Add functionalities for distributed2d gershgorin

src/C-interface/bml_import.c

@@ -40,8 +40,8 @@ bml_import_from_dense(
 #ifdef DO_MPI
     if (distrib_mode == distributed)
         return bml_import_from_dense_distributed2d(matrix_type,
-                                                   matrix_precision,
-                                                   order, N, A, threshold, M);
+                                                   matrix_precision, order, N,
+                                                   A, threshold, M);
     else
 #endif
         switch (matrix_type)

src/C-interface/bml_normalize.c

@@ -6,6 +6,9 @@
 #include "ellsort/bml_normalize_ellsort.h"
 #include "ellblock/bml_normalize_ellblock.h"
 #include "csr/bml_normalize_csr.h"
+#ifdef DO_MPI
+#include "distributed2d/bml_normalize_distributed2d.h"
+#endif
 
 #include <stdlib.h>
 
@@ -40,10 +43,44 @@ bml_normalize(
         case csr:
             bml_normalize_csr(A, mineval, maxeval);
             break;
+#ifdef DO_MPI
+        case distributed2d:
+            return bml_normalize_distributed2d(A, mineval, maxeval);
+            break;
+#endif
+        default:
+            LOG_ERROR("unknown matrix type\n");
+            break;
+    }
+}
+
+void *
+bml_accumulate_offdiag(
+    bml_matrix_t * A,
+    int flag)
+{
+    switch (bml_get_type(A))
+    {
+        case dense:
+            return bml_accumulate_offdiag_dense(A, flag);
+            break;
+        case ellpack:
+            return bml_accumulate_offdiag_ellpack(A, flag);
+            break;
+        case ellsort:
+            return bml_accumulate_offdiag_ellsort(A, flag);
+            break;
+        case ellblock:
+            return bml_accumulate_offdiag_ellblock(A, flag);
+            break;
+        case csr:
+            return bml_accumulate_offdiag_csr(A, flag);
+            break;
         default:
             LOG_ERROR("unknown matrix type\n");
             break;
     }
+    return NULL;
 }
 
 /** Calculate Gershgorin bounds.
@@ -75,6 +112,11 @@ bml_gershgorin(
         case csr:
             return bml_gershgorin_csr(A);
             break;
+#ifdef DO_MPI
+        case distributed2d:
+            return bml_gershgorin_distributed2d(A);
+            break;
+#endif
         default:
             LOG_ERROR("unknown matrix type\n");
             break;

src/C-interface/bml_normalize.h

@@ -20,4 +20,7 @@ void *bml_gershgorin_partial(
     bml_matrix_t * A,
     int nrows);
 
+void *bml_accumulate_offdiag(
+    bml_matrix_t * A,
+    int flag);
 #endif

src/C-interface/csr/bml_normalize_csr.c

@@ -7,6 +7,31 @@
 #include <stdlib.h>
 #include <string.h>
 
+void *
+bml_accumulate_offdiag_csr(
+    bml_matrix_csr_t * A,
+    int flag)
+{
+    switch (A->matrix_precision)
+    {
+        case single_real:
+            return bml_accumulate_offdiag_csr_single_real(A, flag);
+            break;
+        case double_real:
+            return bml_accumulate_offdiag_csr_double_real(A, flag);
+            break;
+        case single_complex:
+            return bml_accumulate_offdiag_csr_single_complex(A, flag);
+            break;
+        case double_complex:
+            return bml_accumulate_offdiag_csr_double_complex(A, flag);
+            break;
+        default:
+            LOG_ERROR("unknown precision\n");
+            break;
+    }
+}
+
 /** Normalize csr matrix given gershgorin bounds.
  *
  *  \ingroup normalize_group

src/C-interface/csr/bml_normalize_csr.h

@@ -28,6 +28,26 @@ void bml_normalize_csr_double_complex(
     double mineval,
     double maxeval);
 
+void *bml_accumulate_offdiag_csr(
+    bml_matrix_csr_t * A,
+    int);
+
+void *bml_accumulate_offdiag_csr_single_real(
+    bml_matrix_csr_t * A,
+    int);
+
+void *bml_accumulate_offdiag_csr_double_real(
+    bml_matrix_csr_t * A,
+    int);
+
+void *bml_accumulate_offdiag_csr_single_complex(
+    bml_matrix_csr_t * A,
+    int);
+
+void *bml_accumulate_offdiag_csr_double_complex(
+    bml_matrix_csr_t * A,
+    int);
+
 void *bml_gershgorin_csr(
     bml_matrix_csr_t * A);
 

src/C-interface/csr/bml_normalize_csr_typed.c

@@ -44,6 +44,32 @@ void TYPED_FUNC(
     bml_add_identity_csr(A, gershfact, threshold);
 }
 
+void *TYPED_FUNC(
+    bml_accumulate_offdiag_csr) (
+    bml_matrix_csr_t * A,
+    int include_diag)
+{
+    int N = A->N_;
+    REAL_T *offdiag_sum = calloc(N, sizeof(REAL_T));
+
+    for (int i = 0; i < N; i++)
+    {
+        double radius = 0.0;
+        int *cols = A->data_[i]->cols_;
+        REAL_T *vals = (REAL_T *) A->data_[i]->vals_;
+        const int annz = A->data_[i]->NNZ_;
+        for (int pos = 0; pos < annz; pos++)
+        {
+            const int j = cols[pos];
+            if ((i != j) || include_diag)
+                radius += (double) ABS(vals[pos]);
+        }
+        offdiag_sum[i] = radius;
+    }
+
+    return offdiag_sum;
+}
+
 /** Calculate Gershgorin bounds for an csr matrix.
  *
  *  \ingroup normalize_group

src/C-interface/dense/bml_normalize_dense.c

@@ -10,6 +10,31 @@
 #include <stdlib.h>
 #include <string.h>
 
+void *
+bml_accumulate_offdiag_dense(
+    bml_matrix_dense_t * A,
+    int flag)
+{
+    switch (A->matrix_precision)
+    {
+        case single_real:
+            return bml_accumulate_offdiag_dense_single_real(A, flag);
+            break;
+        case double_real:
+            return bml_accumulate_offdiag_dense_double_real(A, flag);
+            break;
+        case single_complex:
+            return bml_accumulate_offdiag_dense_single_complex(A, flag);
+            break;
+        case double_complex:
+            return bml_accumulate_offdiag_dense_double_complex(A, flag);
+            break;
+        default:
+            LOG_ERROR("unknown precision\n");
+            break;
+    }
+}
+
 /** Normalize dense matrix given Gershgorin bounds.
  *
  *  \ingroup normalize_group

src/C-interface/dense/bml_normalize_dense.h

@@ -28,6 +28,26 @@ void bml_normalize_dense_double_complex(
     double mineval,
     double maxeval);
 
+void *bml_accumulate_offdiag_dense(
+    bml_matrix_dense_t * A,
+    int);
+
+void *bml_accumulate_offdiag_dense_single_real(
+    bml_matrix_dense_t * A,
+    int);
+
+void *bml_accumulate_offdiag_dense_double_real(
+    bml_matrix_dense_t * A,
+    int);
+
+void *bml_accumulate_offdiag_dense_single_complex(
+    bml_matrix_dense_t * A,
+    int);
+
+void *bml_accumulate_offdiag_dense_double_complex(
+    bml_matrix_dense_t * A,
+    int);
+
 void *bml_gershgorin_dense(
     bml_matrix_dense_t * A);
 

src/C-interface/dense/bml_normalize_dense_typed.c

@@ -47,6 +47,30 @@ void TYPED_FUNC(
     TYPED_FUNC(bml_scale_add_identity_dense) (A, alpha, beta);
 }
 
+void *TYPED_FUNC(
+    bml_accumulate_offdiag_dense) (
+    bml_matrix_dense_t * A,
+    int include_diag)
+{
+    int N = A->N;
+    REAL_T *offdiag_sum = calloc(N, sizeof(REAL_T));
+    REAL_T *A_value = (REAL_T *) A->matrix;
+
+    for (int i = 0; i < N; i++)
+    {
+        double radius = 0.0;
+        for (int j = 0; j < N; j++)
+        {
+            int ind = ROWMAJOR(i, j, N, N);
+            if ((i != j) || include_diag)
+                radius += (double) ABS(A_value[ind]);
+        }
+        offdiag_sum[i] = radius;
+    }
+
+    return offdiag_sum;
+}
+
 /** Calculate Gershgorin bounds for a dense matrix.
  *
  *  \ingroup gershgorin_group

src/C-interface/distributed2d/CMakeLists.txt

@@ -9,6 +9,7 @@ set(HEADERS-DISTRIBUTED
   bml_add_distributed2d.h
   bml_introspection_distributed2d.h
   bml_multiply_distributed2d.h
+  bml_normalize_distributed2d.h
   bml_norm_distributed2d.h
   bml_getters_distributed2d.h
   bml_setters_distributed2d.h
@@ -29,6 +30,7 @@ set(SOURCES-DISTRIBUTED
   bml_introspection_distributed2d.c
   bml_import_distributed2d.c
   bml_multiply_distributed2d.c
+  bml_normalize_distributed2d.c
   bml_norm_distributed2d.c
   bml_getters_distributed2d.c
   bml_setters_distributed2d.c
@@ -51,6 +53,7 @@ set(SOURCES-DISTRIBUTED-TYPED
   bml_export_distributed2d_typed.c
   bml_import_distributed2d_typed.c
   bml_multiply_distributed2d_typed.c
+  bml_normalize_distributed2d_typed.c
   bml_diagonalize_distributed2d_typed.c
   bml_getters_distributed2d_typed.c
   bml_setters_distributed2d_typed.c

src/C-interface/distributed2d/bml_allocate_distributed2d.c

@@ -12,6 +12,8 @@
 
 /* MPI communicator for all the distributed2d matrices */
 static MPI_Comm s_comm = MPI_COMM_NULL;
+static MPI_Comm r_comm = MPI_COMM_NULL;
+static MPI_Comm c_comm = MPI_COMM_NULL;
 
 void
 bml_setcomm_distributed2d(
@@ -66,6 +68,9 @@ bml_setup_distributed2d(
     A->n = N / p;
     A->matrix_type = distributed2d;
 
+    MPI_Comm_split(A->comm, A->myprow, A->mypcol, &A->row_comm);
+    MPI_Comm_split(A->comm, A->mypcol, A->myprow, &A->col_comm);
+
     assert(A->n * p == N);
 }
 

src/C-interface/distributed2d/bml_normalize_distributed2d.c

@@ -0,0 +1,76 @@
+#include "../bml_logger.h"
+#include "../bml_normalize.h"
+#include "../bml_types.h"
+#include "../bml_introspection.h"
+#include "bml_normalize_distributed2d.h"
+#include "bml_types_distributed2d.h"
+
+#include <stdlib.h>
+#include <string.h>
+
+/** Normalize distributed2d matrix given gershgorin bounds.
+ *
+ *  \ingroup normalize_group
+ *
+ *  \param A The matrix
+ *  \param mineval Calculated min value
+ *  \param maxeval Calculated max value
+ */
+void
+bml_normalize_distributed2d(
+    bml_matrix_distributed2d_t * A,
+    double mineval,
+    double maxeval)
+{
+    switch (bml_get_precision(A))
+    {
+        case single_real:
+            bml_normalize_distributed2d_single_real(A, mineval, maxeval);
+            break;
+        case double_real:
+            bml_normalize_distributed2d_double_real(A, mineval, maxeval);
+            break;
+        case single_complex:
+            bml_normalize_distributed2d_single_complex(A, mineval, maxeval);
+            break;
+        case double_complex:
+            bml_normalize_distributed2d_double_complex(A, mineval, maxeval);
+            break;
+        default:
+            LOG_ERROR("unknown precision\n");
+            break;
+    }
+}
+
+/** Calculate gershgorin bounds for an distributed2d matrix.
+ *
+ *  \ingroup normalize_group
+ *
+ *  \param A The matrix
+ *  returns mineval Calculated min value
+ *  returns maxeval Calculated max value
+ */
+void *
+bml_gershgorin_distributed2d(
+    bml_matrix_distributed2d_t * A)
+{
+    switch (bml_get_precision(A))
+    {
+        case single_real:
+            return bml_gershgorin_distributed2d_single_real(A);
+            break;
+        case double_real:
+            return bml_gershgorin_distributed2d_double_real(A);
+            break;
+        case single_complex:
+            return bml_gershgorin_distributed2d_single_complex(A);
+            break;
+        case double_complex:
+            return bml_gershgorin_distributed2d_double_complex(A);
+            break;
+        default:
+            LOG_ERROR("unknown precision\n");
+            break;
+    }
+    return NULL;
+}