Not using cannons algorithm without MPI

also adding wrapper for the backwords step

Makefile.am

@@ -62,7 +62,7 @@ libelpa@SUFFIX@_private_la_SOURCES = \
   src/elpa2/elpa2.F90 \
   src/elpa_generalized/cannon_forw_real_double.c \
   src/elpa_generalized/cannon_back_real_double.c \
-  src/elpa_generalized/test_c_bindings.c \
+  #src/elpa_generalized/test_c_bindings.c \
   src/helpers/matrix_plot.F90 \
   src/elpa_index.c
 
@@ -601,11 +601,11 @@ single_complex_2stage_banded@SUFFIX@_LDADD = $(test_program_ldadd)
 single_complex_2stage_banded@SUFFIX@_FCFLAGS = $(AM_FCFLAGS) $(FC_MODINC)test_modules $(FC_MODINC)modules
 endif
 
-noinst_PROGRAMS += test_c_bindings@SUFFIX@
-check_SCRIPTS += test_c_bindings@SUFFIX@_default.sh
-test_c_bindings@SUFFIX@_SOURCES = test/Fortran/elpa_generalized/test_bindings.F90
-test_c_bindings@SUFFIX@_LDADD = $(test_program_ldadd) $(FCLIBS)
-test_c_bindings@SUFFIX@_FCFLAGS = $(AM_FCFLAGS) $(FC_MODINC)test_modules $(FC_MODINC)modules
+#noinst_PROGRAMS += test_c_bindings@SUFFIX@
+#check_SCRIPTS += test_c_bindings@SUFFIX@_default.sh
+#test_c_bindings@SUFFIX@_SOURCES = test/Fortran/elpa_generalized/test_bindings.F90
+#test_c_bindings@SUFFIX@_LDADD = $(test_program_ldadd) $(FCLIBS)
+#test_c_bindings@SUFFIX@_FCFLAGS = $(AM_FCFLAGS) $(FC_MODINC)test_modules $(FC_MODINC)modules
 
 #noinst_PROGRAMS += test_c_cannon@SUFFIX@
 #check_SCRIPTS += test_c_cannon@SUFFIX@_default.sh

src/elpa_generalized/cannon_back_real_double.c

@@ -2,10 +2,11 @@
 
 #include <stdio.h>
 #include <stdlib.h>
+#include <math.h>
+
+// most of the file is not compiled if not using MPI
 #ifdef WITH_MPI
 #include <mpi.h>
-#endif
-#include <math.h>
 
 //#include <elpa/elpa.h>
 //#include <elpa/elpa_generated.h>
@@ -435,4 +436,43 @@ void d_cannons_triang_rectangular(double* U, double* B, int np_rows, int np_cols
       free(Buf_U);
 }
 
+#endif
+
+//***********************************************************************************************************
+/*
+!f> interface
+!f>   subroutine cannons_triang_rectangular(U, B, local_rows, local_cols, np_rows, np_cols, my_prow, my_pcol, u_desc, b_desc, &
+!f>                                Res, row_comm, col_comm) &
+!f>                             bind(C, name="d_cannons_triang_rectangular_c")
+!f>     use, intrinsic :: iso_c_binding
+!f>     real(c_double)                        :: U(local_rows, local_cols), B(local_rows, local_cols), Res(local_rows, local_cols)
+!f>     integer(kind=c_int)                   :: u_desc(9), b_desc(9)
+!f>     integer(kind=c_int),value             :: local_rows, local_cols
+!f>     integer(kind=c_int),value     :: np_rows, np_cols, my_prow, my_pcol, row_comm, col_comm
+!f>   end subroutine
+!f> end interface
+*/
+void d_cannons_triang_rectantular_c(double* U, double* B, int local_rows, int local_cols, int np_rows, int np_cols,
+                                     int my_prow, int my_pcol, int* u_desc, int* b_desc, double *Res, int row_comm, int col_comm)
+{
+#ifdef WITH_MPI
+  MPI_Comm c_row_comm = MPI_Comm_f2c(row_comm);
+  MPI_Comm c_col_comm = MPI_Comm_f2c(col_comm);
+  
+  //int c_my_prow, c_my_pcol;
+  //MPI_Comm_rank(c_row_comm, &c_my_prow);
+  //MPI_Comm_rank(c_col_comm, &c_my_pcol);
+  //printf("FORT<->C row: %d<->%d, col: %d<->%d\n", my_prow, c_my_prow, my_pcol, c_my_pcol);
+
+  // BEWARE
+  // in the cannons algorithm, column and row communicators are exchanged
+  // What we usually call row_comm in elpa, is thus passed to col_comm parameter of the function and vice versa
+  // (order is swapped in the following call)
+  // It is a bit unfortunate, maybe it should be changed in the Cannon algorithm to comply with ELPA standard notation?
+  d_cannons_triang_rectangular(U, B, np_rows, np_cols, my_prow, my_pcol, u_desc, b_desc, Res, c_col_comm, c_row_comm);
+#else
+  printf("Internal error: Cannons algorithm should not be called without MPI, stopping...\n");
+  exit(1);
+#endif
+}
 

src/elpa_generalized/cannon_forw_real_double.c

 #include "config-f90.h"
 #include <stdio.h>
 #include <stdlib.h>
+#include <math.h>
+
+// most of the file is not compiled if not using MPI
 #ifdef WITH_MPI
 #include <mpi.h>
-#endif
-#include <math.h>
 
 //#include <elpa/elpa.h>
 //#include <elpa/elpa_generated.h>
@@ -896,13 +897,8 @@ void d_cannons_reduction(double* A, double* U, int np_rows, int np_cols, int my_
    if(ratio != 1)
       free(Buf_A);
    free(U_stored);
-
-   double num = 7.1354857767573481E-003;
-   double eps = 0.0000001;
-   for (int i = 0; i < na_rows * na_cols; i++)
-     if(Res[i] < num + eps && Res[i] > num - eps)
-       printf("end original C code, Res(%d)  %g, size mat %d, %d\n",i,  Res[i], na_rows, na_cols);
 }
+#endif
 
 //***********************************************************************************************************
 /*
@@ -922,8 +918,7 @@ void d_cannons_reduction(double* A, double* U, int np_rows, int np_cols, int my_
 void d_cannons_reduction_c(double* A, double* U, int local_rows, int local_cols, int np_rows, int np_cols, int my_prow, int my_pcol, int* a_desc,
                          double *Res, int ToStore, int row_comm, int col_comm)
 {
-  //printf("%d, %d, %d, %d, %lf, %lf, %lf, %lf, com: %d, %d\n", np_rows, np_cols, my_prow, my_pcol, A[0], A[1], U[0], U[1], row_comm, col_comm);
-
+#ifdef WITH_MPI
   MPI_Comm c_row_comm = MPI_Comm_f2c(row_comm);
   MPI_Comm c_col_comm = MPI_Comm_f2c(col_comm);
   
@@ -938,5 +933,9 @@ void d_cannons_reduction_c(double* A, double* U, int local_rows, int local_cols,
   // (order is swapped in the following call)
   // It is a bit unfortunate, maybe it should be changed in the Cannon algorithm to comply with ELPA standard notation?
   d_cannons_reduction(A, U, np_rows, np_cols, my_prow, my_pcol, a_desc, Res, ToStore, c_col_comm, c_row_comm);
+#else
+  printf("Internal error: Cannons algorithm should not be called without MPI, stopping...\n");
+  exit(1);
+#endif
 }
 

src/elpa_impl.F90

@@ -53,9 +53,7 @@ module elpa_impl
   use elpa2_impl
   use elpa1_impl
   use elpa1_auxiliary_impl
-#ifdef WITH_MPI
   use elpa_mpi
-#endif
   use elpa_generated_fortran_interfaces
   use elpa_utilities, only : error_unit
 

src/elpa_impl_generalized_transform_template.F90

@@ -42,6 +42,14 @@
      use_cannon = 0
 #endif
 
+#if !defined(WITH_MPI)
+     if(my_p == 0) then
+       write(*,*) "Cannons algorithm can be used with MPI"
+       write(*,*) "Switching to elpa Hermitian and scalapack"
+     end if
+     use_cannon = 0
+#endif
+
      if (mod(np_cols, np_rows) /= 0) then
        if(my_p == 0) then
          write(*,*) "To use Cannons algorithm, np_cols must be a multiple of np_rows."