Fix column-major parameter ordering in parameterized matmul

CVXPY sends parameter values in Fortran (column-major) order, but the
matmul refresh functions assumed row-major/CSR order via raw memcpy.
This produced incorrect matrix values for non-symmetric matrices.

For sparse matrices, iterate the CSR pattern and index into the
column-major source array. For dense matrices, exploit the fact that
column-major A is row-major A^T to memcpy directly into AT, then
transpose to get A.

Also fixes a latent bug where sparse update_values would blindly copy
the first nnz values from the full d1*d2 parameter array, which is
wrong for matrices with structural zeros.

Adds tests for rectangular (3x2) and sparse (3x3 with zeros) cases.

Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>

Author: Transurgeon

src/atoms/affine/left_matmul.c

@@ -170,24 +170,39 @@ static void refresh_sparse_left(left_matmul_expr *lnode)
 {
     Sparse_Matrix *sm_A = (Sparse_Matrix *) lnode->A;
     Sparse_Matrix *sm_AT = (Sparse_Matrix *) lnode->AT;
-    lnode->A->update_values(lnode->A, lnode->param_source->value);
+    CSR_Matrix *csr = sm_A->csr;
+    const double *vals = lnode->param_source->value;
+    int d1 = lnode->param_source->d1;
+
+    /* Parameter values are column-major; extract into CSR data order */
+    for (int row = 0; row < csr->m; row++)
+    {
+        for (int k = csr->p[row]; k < csr->p[row + 1]; k++)
+        {
+            int col = csr->i[k];
+            csr->x[k] = vals[col * d1 + row];
+        }
+    }
     /* Recompute AT values from A */
     AT_fill_values(sm_A->csr, sm_AT->csr, lnode->base.work->iwork);
 }
 
 static void refresh_dense_left(left_matmul_expr *lnode)
 {
     Dense_Matrix *dm_A = (Dense_Matrix *) lnode->A;
+    Dense_Matrix *dm_AT = (Dense_Matrix *) lnode->AT;
     int m = dm_A->base.m;
     int n = dm_A->base.n;
-    lnode->A->update_values(lnode->A, lnode->param_source->value);
-    /* Recompute AT data (transpose of row-major A) */
-    Dense_Matrix *dm_AT = (Dense_Matrix *) lnode->AT;
+    const double *vals = lnode->param_source->value;
+
+    /* Column-major A is row-major AT; copy directly into AT */
+    memcpy(dm_AT->x, vals, m * n * sizeof(double));
+    /* Transpose AT to get row-major A */
     for (int i = 0; i < m; i++)
     {
         for (int j = 0; j < n; j++)
         {
-            dm_AT->x[j * m + i] = dm_A->x[i * n + j];
+            dm_A->x[i * n + j] = dm_AT->x[j * m + i];
         }
     }
 }

src/atoms/affine/right_matmul.c

@@ -40,8 +40,19 @@ static void refresh_sparse_right(left_matmul_expr *lnode)
 {
     Sparse_Matrix *sm_AT_inner = (Sparse_Matrix *) lnode->A;
     Sparse_Matrix *sm_A_inner = (Sparse_Matrix *) lnode->AT;
-    /* lnode->AT holds the original A; update its values from param */
-    lnode->AT->update_values(lnode->AT, lnode->param_source->value);
+    CSR_Matrix *csr_A = sm_A_inner->csr;
+    const double *vals = lnode->param_source->value;
+    int d1 = lnode->param_source->d1;
+
+    /* Parameter values are column-major; extract into CSR data order */
+    for (int row = 0; row < csr_A->m; row++)
+    {
+        for (int k = csr_A->p[row]; k < csr_A->p[row + 1]; k++)
+        {
+            int col = csr_A->i[k];
+            csr_A->x[k] = vals[col * d1 + row];
+        }
+    }
     /* Recompute A^T (lnode->A) from A (lnode->AT) */
     AT_fill_values(sm_A_inner->csr, sm_AT_inner->csr, lnode->base.work->iwork);
 }
@@ -50,16 +61,19 @@ static void refresh_dense_right(left_matmul_expr *lnode)
 {
     Dense_Matrix *dm_AT_inner = (Dense_Matrix *) lnode->A;
     Dense_Matrix *dm_A_inner = (Dense_Matrix *) lnode->AT;
-    int m_orig = dm_A_inner->base.m; /* original A is m x n */
+    int m_orig = dm_A_inner->base.m;
     int n_orig = dm_A_inner->base.n;
-    /* Update original A (inner's AT) from param values */
-    lnode->AT->update_values(lnode->AT, lnode->param_source->value);
-    /* Recompute A^T (inner's A) from A */
+    const double *vals = lnode->param_source->value;
+
+    /* Column-major A is row-major AT; copy directly into inner A (= A^T) */
+    memcpy(dm_AT_inner->x, vals, m_orig * n_orig * sizeof(double));
+    /* Transpose to get row-major A (inner's AT = original A) */
     for (int i = 0; i < m_orig; i++)
     {
         for (int j = 0; j < n_orig; j++)
         {
-            dm_AT_inner->x[j * m_orig + i] = dm_A_inner->x[i * n_orig + j];
+            dm_A_inner->x[i * n_orig + j] =
+                dm_AT_inner->x[j * m_orig + i];
         }
     }
 }

tests/all_tests.c

@@ -354,6 +354,8 @@ int main(void)
     mu_run_test(test_param_vector_mult_problem, tests_run);
     mu_run_test(test_param_left_matmul_problem, tests_run);
     mu_run_test(test_param_right_matmul_problem, tests_run);
+    mu_run_test(test_param_left_matmul_rectangular, tests_run);
+    mu_run_test(test_param_left_matmul_sparse, tests_run);
     mu_run_test(test_param_fixed_skip_in_update, tests_run);
 #endif /* PROFILE_ONLY */
 

tests/problem/test_param_prob.h

@@ -164,14 +164,14 @@ const char *test_param_vector_mult_problem(void)
  * Test 3: left_param_matmul in constraint
  *
  * Problem: minimize sum(x), subject to A @ x, x size 2, A is 2x2
- *   A is a 2x2 matrix parameter (param_id=0, size=4, CSR data order)
- *   A = [[1,2],[3,4]] → CSR data order theta = [1,2,3,4]
+ *   A is a 2x2 matrix parameter (param_id=0, size=4, column-major order)
+ *   A = [[1,2],[3,4]] → column-major theta = [1,3,2,4]
  *
  * At x=[1,2]:
  *   constraint_values = [1*1+2*2, 3*1+4*2] = [5, 11]
  *   jacobian = [[1,2],[3,4]]
  *
- * After update A = [[5,6],[7,8]] → theta = [5,6,7,8]:
+ * After update A = [[5,6],[7,8]] → column-major theta = [5,7,6,8]:
  *   constraint_values = [5*1+6*2, 7*1+8*2] = [17, 23]
  *   jacobian = [[5,6],[7,8]]
  */
@@ -208,8 +208,8 @@ const char *test_param_left_matmul_problem(void)
     problem_register_params(prob, param_nodes, 1);
     problem_init_derivatives(prob);
 
-    /* Set A = [[1,2],[3,4]], CSR data order: [1,2,3,4] */
-    double theta[4] = {1.0, 2.0, 3.0, 4.0};
+    /* Set A = [[1,2],[3,4]], column-major: [1,3,2,4] */
+    double theta[4] = {1.0, 3.0, 2.0, 4.0};
     problem_update_params(prob, theta);
 
     double u[2] = {1.0, 2.0};
@@ -233,8 +233,8 @@ const char *test_param_left_matmul_problem(void)
     double expected_x[4] = {1.0, 2.0, 3.0, 4.0};
     mu_assert("jac->x wrong (A1)", cmp_double_array(jac->x, expected_x, 4));
 
-    /* Update A = [[5,6],[7,8]], CSR data order: [5,6,7,8] */
-    double theta2[4] = {5.0, 6.0, 7.0, 8.0};
+    /* Update A = [[5,6],[7,8]], column-major: [5,7,6,8] */
+    double theta2[4] = {5.0, 7.0, 6.0, 8.0};
     problem_update_params(prob, theta2);
 
     problem_constraint_forward(prob, u);
@@ -256,14 +256,14 @@ const char *test_param_left_matmul_problem(void)
  * Test 4: right_param_matmul in constraint
  *
  * Problem: minimize sum(x), subject to x @ A, x size 1x2, A is 2x2
- *   A is a 2x2 matrix parameter (param_id=0, size=4, CSR data order)
- *   A = [[1,2],[3,4]] → CSR data order theta = [1,2,3,4]
+ *   A is a 2x2 matrix parameter (param_id=0, size=4, column-major order)
+ *   A = [[1,2],[3,4]] → column-major theta = [1,3,2,4]
  *
  * At x=[1,2]:
  *   constraint_values = [1*1+2*3, 1*2+2*4] = [7, 10]
  *   jacobian = [[1,3],[2,4]] = A^T
  *
- * After update A = [[5,6],[7,8]] → theta = [5,6,7,8]:
+ * After update A = [[5,6],[7,8]] → column-major theta = [5,7,6,8]:
  *   constraint_values = [1*5+2*7, 1*6+2*8] = [19, 22]
  *   jacobian = [[5,7],[6,8]] = A^T
  */
@@ -300,8 +300,8 @@ const char *test_param_right_matmul_problem(void)
     problem_register_params(prob, param_nodes, 1);
     problem_init_derivatives(prob);
 
-    /* Set A = [[1,2],[3,4]], CSR data order: [1,2,3,4] */
-    double theta[4] = {1.0, 2.0, 3.0, 4.0};
+    /* Set A = [[1,2],[3,4]], column-major: [1,3,2,4] */
+    double theta[4] = {1.0, 3.0, 2.0, 4.0};
     problem_update_params(prob, theta);
 
     double u[2] = {1.0, 2.0};
@@ -325,8 +325,8 @@ const char *test_param_right_matmul_problem(void)
     double expected_x[4] = {1.0, 3.0, 2.0, 4.0};
     mu_assert("jac->x wrong (A1)", cmp_double_array(jac->x, expected_x, 4));
 
-    /* Update A = [[5,6],[7,8]], CSR data order: [5,6,7,8] */
-    double theta2[4] = {5.0, 6.0, 7.0, 8.0};
+    /* Update A = [[5,6],[7,8]], column-major: [5,7,6,8] */
+    double theta2[4] = {5.0, 7.0, 6.0, 8.0};
     problem_update_params(prob, theta2);
 
     problem_constraint_forward(prob, u);
@@ -425,4 +425,137 @@ const char *test_param_fixed_skip_in_update(void)
     return 0;
 }
 
+/*
+ * Test 6: left_param_matmul with rectangular 3x2 matrix
+ *
+ * Problem: minimize sum(x), subject to A @ x, x size 2, A is 3x2
+ *   A = [[1,2],[3,4],[5,6]] → column-major theta = [1,3,5,2,4,6]
+ *
+ * At x=[1,2]:
+ *   constraint_values = [1+4, 3+8, 5+12] = [5, 11, 17]
+ *   jacobian = [[1,2],[3,4],[5,6]]
+ */
+const char *test_param_left_matmul_rectangular(void)
+{
+    int n_vars = 2;
+
+    /* Objective: sum(x) */
+    expr *x_obj = new_variable(2, 1, 0, n_vars);
+    expr *objective = new_sum(x_obj, -1);
+
+    /* Constraint: A @ x, A is 3x2 */
+    expr *x_con = new_variable(2, 1, 0, n_vars);
+    expr *A_param = new_parameter(3, 2, 0, n_vars, NULL);
+
+    /* Dense 3x2 CSR */
+    CSR_Matrix *A = new_csr_matrix(3, 2, 6);
+    int Ap[4] = {0, 2, 4, 6};
+    int Ai[6] = {0, 1, 0, 1, 0, 1};
+    double Ax[6] = {0, 0, 0, 0, 0, 0};
+    memcpy(A->p, Ap, 4 * sizeof(int));
+    memcpy(A->i, Ai, 6 * sizeof(int));
+    memcpy(A->x, Ax, 6 * sizeof(double));
+
+    expr *constraint = new_left_matmul(A_param, x_con, A);
+    free_csr_matrix(A);
+
+    expr *constraints[1] = {constraint};
+
+    problem *prob = new_problem(objective, constraints, 1, false);
+
+    expr *param_nodes[1] = {A_param};
+    problem_register_params(prob, param_nodes, 1);
+    problem_init_derivatives(prob);
+
+    /* Set A = [[1,2],[3,4],[5,6]], column-major: [1,3,5,2,4,6] */
+    double theta[6] = {1.0, 3.0, 5.0, 2.0, 4.0, 6.0};
+    problem_update_params(prob, theta);
+
+    double u[2] = {1.0, 2.0};
+    problem_constraint_forward(prob, u);
+    problem_jacobian(prob);
+
+    double expected_cv[3] = {5.0, 11.0, 17.0};
+    mu_assert("rect constraint values wrong",
+              cmp_double_array(prob->constraint_values, expected_cv, 3));
+
+    CSR_Matrix *jac = prob->jacobian;
+    double expected_x[6] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0};
+    mu_assert("rect jac->x wrong", cmp_double_array(jac->x, expected_x, 6));
+
+    free_problem(prob);
+
+    return 0;
+}
+
+/*
+ * Test 7: left_param_matmul with sparse matrix (structural zeros)
+ *
+ * Problem: minimize sum(x), subject to A @ x, x size 3, A is 3x3 sparse
+ *   A = [[1, 0, 2],
+ *        [0, 3, 0],
+ *        [4, 0, 5]]
+ *   CSR: p=[0,2,3,5], i=[0,2, 1, 0,2], nnz=5
+ *   Column-major theta (full 9 values): [1,0,4, 0,3,0, 2,0,5]
+ *
+ * At x=[1,2,3]:
+ *   constraint_values = [1+6, 6, 4+15] = [7, 6, 19]
+ *   jacobian = A
+ */
+const char *test_param_left_matmul_sparse(void)
+{
+    int n_vars = 3;
+
+    /* Objective: sum(x) */
+    expr *x_obj = new_variable(3, 1, 0, n_vars);
+    expr *objective = new_sum(x_obj, -1);
+
+    /* Constraint: A @ x, A is 3x3 sparse */
+    expr *x_con = new_variable(3, 1, 0, n_vars);
+    expr *A_param = new_parameter(3, 3, 0, n_vars, NULL);
+
+    CSR_Matrix *A = new_csr_matrix(3, 3, 5);
+    int Ap[4] = {0, 2, 3, 5};
+    int Ai[5] = {0, 2, 1, 0, 2};
+    double Ax[5] = {0, 0, 0, 0, 0};
+    memcpy(A->p, Ap, 4 * sizeof(int));
+    memcpy(A->i, Ai, 5 * sizeof(int));
+    memcpy(A->x, Ax, 5 * sizeof(double));
+
+    expr *constraint = new_left_matmul(A_param, x_con, A);
+    free_csr_matrix(A);
+
+    expr *constraints[1] = {constraint};
+
+    problem *prob = new_problem(objective, constraints, 1, false);
+
+    expr *param_nodes[1] = {A_param};
+    problem_register_params(prob, param_nodes, 1);
+    problem_init_derivatives(prob);
+
+    /* A = [[1,0,2],[0,3,0],[4,0,5]], column-major: [1,0,4, 0,3,0, 2,0,5] */
+    double theta[9] = {1.0, 0.0, 4.0, 0.0, 3.0, 0.0, 2.0, 0.0, 5.0};
+    problem_update_params(prob, theta);
+
+    double u[3] = {1.0, 2.0, 3.0};
+    problem_constraint_forward(prob, u);
+    problem_jacobian(prob);
+
+    double expected_cv[3] = {7.0, 6.0, 19.0};
+    mu_assert("sparse constraint values wrong",
+              cmp_double_array(prob->constraint_values, expected_cv, 3));
+
+    CSR_Matrix *jac = prob->jacobian;
+    mu_assert("sparse jac nnz wrong", jac->nnz == 5);
+
+    /* CSR data order: row0=[1,2], row1=[3], row2=[4,5] */
+    double expected_x[5] = {1.0, 2.0, 3.0, 4.0, 5.0};
+    mu_assert("sparse jac->x wrong",
+              cmp_double_array(jac->x, expected_x, 5));
+
+    free_problem(prob);
+
+    return 0;
+}
+
 #endif /* TEST_PARAM_PROB_H */