Cleanup testing of bml_add

We should be testing all APIs of `bml_add`. This change adds the missing
API tests.

src/C-interface/dense/bml_add_dense_typed.c

@@ -85,7 +85,6 @@ double TYPED_FUNC(
     shared(B_matrix, A_localRowMin, A_localRowMax) \
     shared(N, myRank) \
     reduction(+:trnorm)
-    //for (int i = 0; i < N * N; i++)
     for (int i = A_localRowMin[myRank] * N; i < A_localRowMax[myRank] * N;
          i++)
     {

tests/CMakeLists.txt

@@ -1,7 +1,7 @@
 include_directories(${CMAKE_SOURCE_DIR}/src/C-interface)
 
 set(SOURCES_TYPED
-  add_matrix_typed.c
+  test_bml_add_typed.c
   adjacency_matrix_typed.c
   adjungate_triangle_matrix_typed.c
   allocate_matrix_typed.c
@@ -39,7 +39,7 @@ set_target_properties(bmltests
   POSITION_INDEPENDENT_CODE yes)
 
 add_executable(bml-test
-  add_matrix.c
+  test_bml_add.c
   adjacency_matrix.c
   adjungate_triangle_matrix.c
   allocate_matrix.c

tests/bml_test.c

@@ -66,7 +66,7 @@ const char *test_description[] = {
 };
 
 const test_function_t testers[] = {
-    test_add,
+    test_bml_add,
     test_adjacency,
     test_adjungate_triangle,
     test_allocate,

tests/bml_test.h

@@ -10,7 +10,6 @@ typedef int (
     const bml_matrix_precision_t matrix_precision,
     const int M);
 
-#include "add_matrix.h"
 #include "adjacency_matrix.h"
 #include "adjungate_triangle_matrix.h"
 #include "allocate_matrix.h"
@@ -27,6 +26,7 @@ typedef int (
 #include "scale_matrix.h"
 #include "set_row.h"
 #include "submatrix_matrix.h"
+#include "test_bml_add.h"
 #include "test_bml_gemm.h"
 #include "threshold_matrix.h"
 #include "trace_matrix.h"

tests/add_matrix.c → tests/test_bml_add.c

@@ -4,7 +4,7 @@
 #include <stdio.h>
 
 int
-test_add(
+test_bml_add(
     const int N,
     const bml_matrix_type_t matrix_type,
     const bml_matrix_precision_t matrix_precision,
@@ -13,18 +13,20 @@ test_add(
     switch (matrix_precision)
     {
         case single_real:
-            return test_add_single_real(N, matrix_type, matrix_precision, M);
+            return test_bml_add_single_real(N, matrix_type, matrix_precision,
+                                            M);
             break;
         case double_real:
-            return test_add_double_real(N, matrix_type, matrix_precision, M);
+            return test_bml_add_double_real(N, matrix_type, matrix_precision,
+                                            M);
             break;
         case single_complex:
-            return test_add_single_complex(N, matrix_type, matrix_precision,
-                                           M);
+            return test_bml_add_single_complex(N, matrix_type,
+                                               matrix_precision, M);
             break;
         case double_complex:
-            return test_add_double_complex(N, matrix_type, matrix_precision,
-                                           M);
+            return test_bml_add_double_complex(N, matrix_type,
+                                               matrix_precision, M);
             break;
         default:
             fprintf(stderr, "unknown matrix precision\n");

tests/add_matrix.h → tests/test_bml_add.h

@@ -1,33 +1,33 @@
-#ifndef __ADD_MATRIX_H
-#define __ADD_MATRIX_H
+#ifndef __TEST_BML_ADD_H
+#define __TEST_BML_ADD_H
 
 #include <bml.h>
 
-int test_add(
+int test_bml_add(
     const int N,
     const bml_matrix_type_t matrix_type,
     const bml_matrix_precision_t matrix_precision,
     const int M);
 
-int test_add_single_real(
+int test_bml_add_single_real(
     const int N,
     const bml_matrix_type_t matrix_type,
     const bml_matrix_precision_t matrix_precision,
     const int M);
 
-int test_add_double_real(
+int test_bml_add_double_real(
     const int N,
     const bml_matrix_type_t matrix_type,
     const bml_matrix_precision_t matrix_precision,
     const int M);
 
-int test_add_single_complex(
+int test_bml_add_single_complex(
     const int N,
     const bml_matrix_type_t matrix_type,
     const bml_matrix_precision_t matrix_precision,
     const int M);
 
-int test_add_double_complex(
+int test_bml_add_double_complex(
     const int N,
     const bml_matrix_type_t matrix_type,
     const bml_matrix_precision_t matrix_precision,

tests/test_bml_add_typed.c

@@ -0,0 +1,195 @@
+#include "bml.h"
+#include "../typed.h"
+#include "../macros.h"
+
+#include <complex.h>
+#include <math.h>
+#include <stdlib.h>
+#include <stdio.h>
+#if defined(SINGLE_REAL) || defined(SINGLE_COMPLEX)
+#define REL_TOL 1e-6
+#else
+#define REL_TOL 1e-12
+#endif
+
+int TYPED_FUNC(
+    test_bml_add) (
+    const int N,
+    const bml_matrix_type_t matrix_type,
+    const bml_matrix_precision_t matrix_precision,
+    const int M)
+{
+    bml_matrix_t *A = NULL;
+    bml_matrix_t *B = NULL;
+    bml_matrix_t *C = NULL;
+
+    REAL_T *A_dense = NULL;
+    REAL_T *B_dense = NULL;
+    REAL_T *C_dense = NULL;
+
+    double norm;
+    double expected_norm;
+
+    double alpha = 1.2;
+    double beta = 0.8;
+    double threshold = 0.0;
+
+    LOG_DEBUG("rel. tolerance = %e\n", REL_TOL);
+
+    A = bml_random_matrix(matrix_type, matrix_precision, N, M, sequential);
+    A_dense = bml_convert_to_dense(A, dense_row_major);
+
+    C = bml_random_matrix(matrix_type, matrix_precision, N, M, sequential);
+    C_dense = bml_convert_to_dense(C, dense_row_major);
+
+    LOG_INFO("Testing bml_add()\n");
+
+    B = bml_copy_new(A);
+    bml_add(B, C, alpha, beta, threshold);
+    B_dense = bml_convert_to_dense(B, dense_row_major);
+
+    bml_print_dense_matrix(N, matrix_precision, dense_row_major, A_dense, 0,
+                           N, 0, N);
+    bml_print_dense_matrix(N, matrix_precision, dense_row_major, B_dense, 0,
+                           N, 0, N);
+    bml_print_dense_matrix(N, matrix_precision, dense_row_major, C_dense, 0,
+                           N, 0, N);
+    for (int i = 0; i < N * N; i++)
+    {
+        double expected = alpha * A_dense[i] + beta * C_dense[i];
+        double rel_diff = ABS((expected - B_dense[i]) / expected);
+        if (rel_diff > REL_TOL)
+        {
+            LOG_ERROR
+                ("matrices are not identical; expected[%d] = %e, B[%d] = %e\n",
+                 i, expected, i, B_dense[i]);
+            return -1;
+        }
+    }
+    bml_free_memory(B_dense);
+    bml_deallocate(&B);
+
+    LOG_INFO("Testing bml_add_norm()\n");
+
+    B = bml_copy_new(A);
+    norm = bml_add_norm(B, C, alpha, beta, threshold);
+    B_dense = bml_convert_to_dense(B, dense_row_major);
+
+    bml_print_dense_matrix(N, matrix_precision, dense_row_major, A_dense, 0,
+                           N, 0, N);
+    bml_print_dense_matrix(N, matrix_precision, dense_row_major, B_dense, 0,
+                           N, 0, N);
+    bml_print_dense_matrix(N, matrix_precision, dense_row_major, C_dense, 0,
+                           N, 0, N);
+    expected_norm = 0;
+    for (int i = 0; i < N * N; i++)
+    {
+        expected_norm += C_dense[i] * C_dense[i];
+        double expected = alpha * A_dense[i] + beta * C_dense[i];
+        double rel_diff = ABS((expected - B_dense[i]) / expected);
+        if (rel_diff > REL_TOL)
+        {
+            LOG_ERROR
+                ("matrices are not identical; expected[%d] = %e, B[%d] = %e\n",
+                 i, expected, i, B_dense[i]);
+            return -1;
+        }
+    }
+    if (ABS(expected_norm - norm) / expected_norm > REL_TOL)
+    {
+        LOG_ERROR("norm mismatch\n");
+        return -1;
+    }
+    bml_free_memory(B_dense);
+    bml_deallocate(&B);
+
+    LOG_INFO("Testing bml_add_identity()\n");
+
+    B = bml_copy_new(A);
+    bml_add_identity(B, beta, threshold);
+    B_dense = bml_convert_to_dense(B, dense_row_major);
+
+    bml_print_dense_matrix(N, matrix_precision, dense_row_major, A_dense, 0,
+                           N, 0, N);
+    bml_print_dense_matrix(N, matrix_precision, dense_row_major, B_dense, 0,
+                           N, 0, N);
+    bml_print_dense_matrix(N, matrix_precision, dense_row_major, C_dense, 0,
+                           N, 0, N);
+    for (int i = 0; i < N; i++)
+    {
+        for (int j = 0; j < N; j++)
+        {
+            double expected = 0;
+            if (i == j)
+            {
+                expected = A_dense[ROWMAJOR(i, i, N, N)] + beta;
+            }
+            else
+            {
+                expected = A_dense[ROWMAJOR(i, j, N, N)];
+            }
+
+            double rel_diff =
+                ABS((expected - B_dense[ROWMAJOR(i, j, N, N)]) / expected);
+            if (rel_diff > REL_TOL)
+            {
+                LOG_ERROR
+                    ("matrices are not identical; expected[%d] = %e, B[%d] = %e\n",
+                     i, expected, i, B_dense[ROWMAJOR(i, j, N, N)]);
+                return -1;
+            }
+        }
+    }
+    bml_free_memory(B_dense);
+    bml_deallocate(&B);
+
+    LOG_INFO("Testing bml_scale_add_identity()\n");
+
+    B = bml_copy_new(A);
+    bml_scale_add_identity(B, alpha, beta, threshold);
+    B_dense = bml_convert_to_dense(B, dense_row_major);
+
+    bml_print_dense_matrix(N, matrix_precision, dense_row_major, A_dense, 0,
+                           N, 0, N);
+    bml_print_dense_matrix(N, matrix_precision, dense_row_major, B_dense, 0,
+                           N, 0, N);
+    bml_print_dense_matrix(N, matrix_precision, dense_row_major, C_dense, 0,
+                           N, 0, N);
+    for (int i = 0; i < N; i++)
+    {
+        for (int j = 0; j < N; j++)
+        {
+            double expected = 0;
+            if (i == j)
+            {
+                expected = alpha * A_dense[ROWMAJOR(i, i, N, N)] + beta;
+            }
+            else
+            {
+                expected = alpha * A_dense[ROWMAJOR(i, j, N, N)];
+            }
+
+            double rel_diff =
+                ABS((expected - B_dense[ROWMAJOR(i, j, N, N)]) / expected);
+            if (rel_diff > REL_TOL)
+            {
+                LOG_ERROR
+                    ("matrices are not identical; expected[%d] = %e, B[%d] = %e\n",
+                     i, expected, i, B_dense[ROWMAJOR(i, j, N, N)]);
+                return -1;
+            }
+        }
+    }
+    bml_free_memory(B_dense);
+    bml_deallocate(&B);
+
+    bml_free_memory(A_dense);
+    bml_free_memory(C_dense);
+
+    bml_deallocate(&A);
+    bml_deallocate(&C);
+
+    LOG_INFO("add matrix test passed\n");
+
+    return 0;
+}