add test

Author: dance858

src/bivariate/right_matmul.c

@@ -54,7 +54,7 @@ expr *new_right_matmul_dense(expr *u, int m, int n, const double *data)
     }
 
     expr *u_transpose = new_transpose(u);
-    expr *left_matmul_node = new_left_matmul_dense(u_transpose, m, n, AT);
+    expr *left_matmul_node = new_left_matmul_dense(u_transpose, n, m, AT);
     expr *node = new_transpose(left_matmul_node);
 
     free(AT);

src/utils/dense_matrix.c

@@ -30,14 +30,22 @@ static void dense_block_left_mult_vec(const Matrix *A, const double *x, double *
 
     /* y = kron(I_p, A) @ x via a single dgemm call:
        Treat x as n x p (column-major blocks) and y as m x p.
-       But x and y are stored as p blocks of length n and m respectively
-       (i.e. block-interleaved). This is the same as treating them as
-       row-major matrices of shape p x n and p x m, so:
+       But x and y are stored as p blocks of length n and m
+       respectively (i.e. block-interleaved). This is the same as
+       treating them as row-major matrices of shape p x n and
+       p x m, so:
        y (p x m) = x (p x n) * A^T (n x m), all row-major.
        cblas with RowMajor: C = alpha * A * B + beta * C
        where A = x (p x n), B = A^T (n x m), C = y (p x m). */
-    cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasTrans, p, m, n, 1.0, x, n, dm->x,
-                n, 0.0, y, m);
+    /* cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasTrans,
+                   p, m, n, 1.0, x, n, dm->x,
+                   n, 0.0, y, m); */
+    for (int b = 0; b < p; b++)
+    {
+        cblas_dgemv(CblasRowMajor, CblasNoTrans, m, n, 1.0,
+                    dm->x, n, x + b * n, 1,
+                    0.0, y + b * m, 1);
+    }
 }
 
 static CSC_Matrix *dense_block_left_mult_sparsity(const Matrix *A,

tests/all_tests.c

@@ -15,6 +15,7 @@
 #include "forward_pass/composite/test_composite.h"
 #include "forward_pass/elementwise/test_exp.h"
 #include "forward_pass/elementwise/test_log.h"
+#include "forward_pass/test_left_matmul_dense.h"
 #include "forward_pass/test_matmul.h"
 #include "forward_pass/test_prod_axis_one.h"
 #include "forward_pass/test_prod_axis_zero.h"
@@ -112,6 +113,7 @@ int main(void)
     mu_run_test(test_forward_prod_axis_zero, tests_run);
     mu_run_test(test_forward_prod_axis_one, tests_run);
     mu_run_test(test_matmul, tests_run);
+    mu_run_test(test_left_matmul_dense, tests_run);
 
     printf("\n--- Jacobian Tests ---\n");
     mu_run_test(test_neg_jacobian, tests_run);

tests/forward_pass/test_left_matmul_dense.h

@@ -0,0 +1,58 @@
+#include <stdio.h>
+#include <stdlib.h>
+
+#include "bivariate.h"
+#include "expr.h"
+#include "minunit.h"
+#include "test_helpers.h"
+
+const char *test_left_matmul_dense(void)
+{
+    /* Test: Z = A @ X where
+     * A is 3x3 (row-major): [1 2 3; 4 5 6; 7 8 9]
+     * X is 3x3 variable (col-major): [1 4 7; 2 5 8; 3 6 9]
+     *
+     * Z = A @ X = [14   32   50]
+     *             [32   77  122]
+     *             [50  122  194]
+     */
+
+    /* Create X variable (3 x 3) */
+    expr *X = new_variable(3, 3, 0, 9);
+
+    /* Constant matrix A in row-major order */
+    double A_data[9] = {1.0, 2.0, 3.0,
+                        4.0, 5.0, 6.0,
+                        7.0, 8.0, 9.0};
+
+    /* Build expression Z = A @ X */
+    expr *Z = new_left_matmul_dense(X, 3, 3, A_data);
+
+    /* Variable values in column-major order */
+    double u[9] = {1.0, 2.0, 3.0,   /* first column */
+                   4.0, 5.0, 6.0,   /* second column */
+                   7.0, 8.0, 9.0};  /* third column */
+
+    /* Evaluate forward pass */
+    Z->forward(Z, u);
+
+    /* Expected result (3 x 3) in column-major order */
+    double expected[9] = {14.0,  32.0,  50.0,   /* first column */
+                          32.0,  77.0,  122.0,  /* second column */
+                          50.0,  122.0, 194.0}; /* third column */
+
+    /* Verify dimensions */
+    mu_assert("left_matmul_dense result should have d1=3",
+              Z->d1 == 3);
+    mu_assert("left_matmul_dense result should have d2=3",
+              Z->d2 == 3);
+    mu_assert("left_matmul_dense result should have size=9",
+              Z->size == 9);
+
+    /* Verify values */
+    mu_assert("Left matmul dense forward pass test failed",
+              cmp_double_array(Z->value, expected, 9));
+
+    free_expr(Z);
+    return 0;
+}