Make random matrices more complex

src/C-interface/bml_transpose.c

@@ -70,3 +70,25 @@ bml_transpose(
             break;
     }
 }
+
+/** Complex conjugate of a matrix.
+ *
+ * \ingroup transpose_group_C
+ *
+ * \param A Matrix to be complex conjugated
+ * \return Complex conjugate of A
+ */
+void
+bml_complex_conjugate(
+    bml_matrix_t * A)
+{
+    switch (bml_get_type(A))
+    {
+        case dense:
+            bml_complex_conjugate_dense(A);
+            break;
+        default:
+            LOG_ERROR("unknown matrix type\n");
+            break;
+    }
+}

src/C-interface/bml_transpose.h

@@ -5,11 +5,12 @@
 
 #include "bml_types.h"
 
-// Transpose A - B = A^T
 bml_matrix_t *bml_transpose_new(
     bml_matrix_t * A);
 
 void bml_transpose(
     bml_matrix_t * A);
 
+void bml_complex_conjugate(bml_matrix_t * A);
+
 #endif

src/C-interface/dense/bml_allocate_dense_typed.c

@@ -11,6 +11,9 @@
 #include "magma_v2.h"
 #endif
 
+#define __USE_MISC
+#include <math.h>
+
 #include <complex.h>
 #include <string.h>
 #include <stdio.h>
@@ -158,11 +161,16 @@ bml_matrix_dense_t *TYPED_FUNC(
     {
         for (int j = 0; j < N; j++)
         {
+            double angle = rand() / (double) RAND_MAX * 2 * M_PI;
 #ifdef BML_USE_MAGMA
             A_dense[ROWMAJOR(i, j, N, N)] =
-                MAGMACOMPLEX(MAKE) (rand() / (double) RAND_MAX, 0.);
+                MAGMACOMPLEX(MAKE) (cos(angle), sin(angle));
 #else
-            A_dense[ROWMAJOR(i, j, N, N)] = rand() / (double) RAND_MAX;
+#if defined(BML_COMPLEX) && (defined(SINGLE_COMPLEX) || defined(DOUBLE_COMPLEX))
+            A_dense[ROWMAJOR(i, j, N, N)] = cos(angle) + sin(angle) * I;
+#else
+            A_dense[ROWMAJOR(i, j, N, N)] = cos(angle);
+#endif
 #endif
         }
     }

src/C-interface/dense/bml_transpose_dense.c

@@ -73,3 +73,34 @@ bml_transpose_dense(
             break;
     }
 }
+
+/** Complex conjugate a matrix in place.
+ *
+ * \ingroup transpose_group
+ *
+ * \param A The matrix to be complex conjugated
+ * \return Complex conjugate of A
+ */
+void
+bml_complex_conjugate_dense(
+    bml_matrix_dense_t * A)
+{
+    switch (A->matrix_precision)
+    {
+        case single_real:
+            bml_complex_conjugate_dense_single_real(A);
+            break;
+        case double_real:
+            bml_complex_conjugate_dense_double_real(A);
+            break;
+        case single_complex:
+            bml_complex_conjugate_dense_single_complex(A);
+            break;
+        case double_complex:
+            bml_complex_conjugate_dense_double_complex(A);
+            break;
+        default:
+            LOG_ERROR("unknown precision\n");
+            break;
+    }
+}

src/C-interface/dense/bml_transpose_dense.h

@@ -33,4 +33,19 @@ void bml_transpose_dense_single_complex(
 void bml_transpose_dense_double_complex(
     bml_matrix_dense_t * A);
 
+void bml_complex_conjugate_dense(
+    bml_matrix_dense_t * A);
+
+void bml_complex_conjugate_dense_single_real(
+    bml_matrix_dense_t * A);
+
+void bml_complex_conjugate_dense_double_real(
+    bml_matrix_dense_t * A);
+
+void bml_complex_conjugate_dense_single_complex(
+    bml_matrix_dense_t * A);
+
+void bml_complex_conjugate_dense_double_complex(
+    bml_matrix_dense_t * A);
+
 #endif

src/C-interface/dense/bml_transpose_dense_typed.c

@@ -105,3 +105,32 @@ void TYPED_FUNC(
     }
 #endif
 }
+
+/** Complex conjugate of a matrix.
+ *
+ * \ingroup transpose_group_C
+ *
+ * \param A Matrix to be complex conjugated
+ * \return Complex conjugate of A
+ */
+void TYPED_FUNC(
+    bml_complex_conjugate_dense) (
+    bml_matrix_dense_t * A)
+{
+    int N = A->N;
+
+    REAL_T *A_matrix = A->matrix;
+    REAL_T tmp;
+
+#if defined(SINGLE_REAL) || defined(DOUBLE_REAL)
+    return;
+#else
+#pragma omp parallel for                        \
+  private(tmp)                                  \
+  shared(N, A_matrix)
+    for (int i = 0; i < N * N; i++)
+    {
+        A_matrix[i] = COMPLEX_CONJUGATE(A_matrix[i]);
+    }
+#endif
+}

src/C-interface/ellblock/bml_allocate_ellblock_typed.c

@@ -6,8 +6,10 @@
 #include "bml_setters_ellblock.h"
 #include "bml_types_ellblock.h"
 
-#include <complex.h>
+#define __USE_MISC
 #include <math.h>
+
+#include <complex.h>
 #include <stdlib.h>
 #include <string.h>
 #include <assert.h>
@@ -212,8 +214,14 @@ bml_matrix_ellblock_t *TYPED_FUNC(
             for (int ii = 0; ii < bsize[ib]; ii++)
                 for (int jj = 0; jj < bsize[jb]; jj++)
                 {
+                    double angle = rand() / (double) RAND_MAX * 2 * M_PI;
+#if defined(BML_COMPLEX) && (defined(SINGLE_COMPLEX) || defined(DOUBLE_COMPLEX))
                     A_value[ROWMAJOR(ii, jj, bsize[ib], bsize[jb])] =
-                        rand() / (REAL_T) RAND_MAX;
+                        cos(angle) + sin(angle) * I;
+#else
+                    A_value[ROWMAJOR(ii, jj, bsize[ib], bsize[jb])] =
+                        cos(angle);
+#endif
                 }
         }
         A_nnzb[ib] = MB;

src/C-interface/ellpack/bml_allocate_ellpack_typed.c

@@ -5,8 +5,10 @@
 #include "bml_allocate_ellpack.h"
 #include "bml_types_ellpack.h"
 
-#include <complex.h>
+#define __USE_MISC
 #include <math.h>
+
+#include <complex.h>
 #include <stdlib.h>
 #include <string.h>
 
@@ -203,13 +205,18 @@ bml_matrix_ellpack_t *TYPED_FUNC(
     REAL_T *A_value = A->value;
     int *A_index = A->index;
     int *A_nnz = A->nnz;
-    const REAL_T INV_RAND_MAX = 1.0 / (REAL_T) RAND_MAX;
+
     for (int i = 0; i < N; i++)
     {
         int jind = 0;
         for (int j = 0; j < M; j++)
         {
-            A_value[ROWMAJOR(i, jind, N, M)] = rand() * INV_RAND_MAX;
+            double angle = rand() / (double) RAND_MAX * 2 * M_PI;
+#if defined(BML_COMPLEX) && (defined(SINGLE_COMPLEX) || defined(DOUBLE_COMPLEX))
+            A_value[ROWMAJOR(i, jind, N, M)] = cos(angle) + sin(angle) * I;
+#else
+            A_value[ROWMAJOR(i, jind, N, M)] = cos(angle);
+#endif
             A_index[ROWMAJOR(i, jind, N, M)] = j;
             jind++;
         }

src/C-interface/ellsort/bml_allocate_ellsort_typed.c

@@ -5,8 +5,10 @@
 #include "bml_allocate_ellsort.h"
 #include "bml_types_ellsort.h"
 
-#include <complex.h>
+#define __USE_MISC
 #include <math.h>
+
+#include <complex.h>
 #include <stdlib.h>
 #include <string.h>
 
@@ -134,15 +136,19 @@ bml_matrix_ellsort_t *TYPED_FUNC(
     int *A_index = A->index;
     int *A_nnz = A->nnz;
 
-    const REAL_T INV_RAND_MAX = 1.0 / (REAL_T) RAND_MAX;
 #pragma omp parallel for shared(A_value, A_index, A_nnz)
     for (int i = 0; i < N; i++)
     {
         int jind = 0;
         for (int j = (i - M / 2 >= 0 ? i - M / 2 : 0);
              j < (i - M / 2 + M <= N ? i - M / 2 + M : N); j++)
         {
-            A_value[ROWMAJOR(i, jind, N, M)] = rand() * INV_RAND_MAX;
+            double angle = rand() / (double) RAND_MAX * 2 * M_PI;
+#if defined(BML_COMPLEX) && (defined(SINGLE_COMPLEX) || defined(DOUBLE_COMPLEX))
+            A_value[ROWMAJOR(i, jind, N, M)] = cos(angle) + sin(angle) * I;
+#else
+            A_value[ROWMAJOR(i, jind, N, M)] = cos(angle);
+#endif
             A_index[ROWMAJOR(i, jind, N, M)] = j;
             jind++;
         }

tests/C-tests/add_matrix_typed.c

@@ -6,7 +6,7 @@
 #include <stdlib.h>
 #include <stdio.h>
 #if defined(SINGLE_REAL) || defined(SINGLE_COMPLEX)
-#define REL_TOL 1e-6
+#define REL_TOL 2e-6
 #else
 #define REL_TOL 1e-12
 #endif
@@ -30,33 +30,44 @@ int TYPED_FUNC(
     double beta = 0.8;
     double threshold = 0.0;
 
-    LOG_DEBUG("rel. tolerance = %e\n", REL_TOL);
+    LOG_INFO("rel. tolerance = %e\n", REL_TOL);
 
     A = bml_random_matrix(matrix_type, matrix_precision, N, M, sequential);
     B = bml_copy_new(A);
     C = bml_random_matrix(matrix_type, matrix_precision, N, M, sequential);
 
+    LOG_INFO("A\n");
+    bml_print_bml_matrix(A, 0, N, 0, N);
+    LOG_INFO("B\n");
+    bml_print_bml_matrix(B, 0, N, 0, N);
+    LOG_INFO("C\n");
+    bml_print_bml_matrix(C, 0, N, 0, N);
+
     bml_add(B, C, alpha, beta, threshold);
 
     A_dense = bml_export_to_dense(A, dense_row_major);
     B_dense = bml_export_to_dense(B, dense_row_major);
     C_dense = bml_export_to_dense(C, dense_row_major);
+
+    LOG_INFO("A_dense\n");
     bml_print_dense_matrix(N, matrix_precision, dense_row_major, A_dense, 0,
                            N, 0, N);
+    LOG_INFO("B_dense\n");
     bml_print_dense_matrix(N, matrix_precision, dense_row_major, B_dense, 0,
                            N, 0, N);
+    LOG_INFO("C_dense\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_val = (expected - B_dense[i]) / expected;
-        double rel_diff = fabs(rel_diff_val);
+        double abs_diff = fabs(expected - B_dense[i]);
+        double rel_diff = abs_diff / 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]);
+                ("matrices are not identical; expected[%d] = %e, B[%d] = %e, abs. diff = %e, |rel. diff| = %e\n",
+                 i, expected, i, B_dense[i], abs_diff, rel_diff);
             return -1;
         }
     }

tests/C-tests/diagonalize_matrix_typed.c

@@ -36,6 +36,7 @@ int TYPED_FUNC(
     A = bml_random_matrix(matrix_type, matrix_precision, N, M, sequential);
 
     A_t = bml_transpose_new(A);
+    bml_complex_conjugate(A_t);
     bml_add(A, A_t, 0.5, 0.5, 0.0);
     bml_print_bml_matrix(A, 0, N, 0, N);
 
@@ -118,7 +119,7 @@ int TYPED_FUNC(
 
     printf("%s\n", "eigenvalues");
     for (int i = 0; i < N; i++)
-        printf("val = %e  i%e\n", REAL_PART(eigenvalues[i]),
+        printf("val = % e%+ei\n", REAL_PART(eigenvalues[i]),
                IMAGINARY_PART(eigenvalues[i]));
 
     aux = bml_transpose_new(eigenvectors);
@@ -129,7 +130,7 @@ int TYPED_FUNC(
         REAL_T *val = bml_get(aux2, i, i);
         if (ABS(*val - (REAL_T) 1.0) > REL_TOL)
         {
-            printf("i = %d, val = %e  i%e\n", i, REAL_PART(*val),
+            printf("i = %d, val = % e%+ei\n", i, REAL_PART(*val),
                    IMAGINARY_PART(*val));
             LOG_ERROR
                 ("Error in matrix diagonalization; eigenvector not normalized\n");

update_tags.sh

@@ -1,8 +1,10 @@
 #!/bin/bash
 
-basedir=$(dirname $0)
-readarray -t files < <(git ls-files '*.c' '*.h' '*.F90')
+set -x
+
+basedir=$(readlink --canonicalize $(dirname $0))
 pushd "${basedir}" || exit
+readarray -t files < <(git ls-files '*.c' '*.h' '*.F90')
 ctags "${files[@]}"
 etags "${files[@]}"
 popd || exit