Support for result checking in parallel matrix-vector. It works!

4 files changed, 58 insertions, 18 deletions

diff --git a/btl/actions/action_parallel_matrix_vector_product.hh b/btl/actions/action_parallel_matrix_vector_product.hh
index 07886a2..bac8bc3 100644
--- a/btl/actions/action_parallel_matrix_vector_product.hh
+++ b/btl/actions/action_parallel_matrix_vector_product.hh
@@ -27,6 +27,11 @@
 #include "init/init_vector.hh"
 #include "init/init_matrix.hh"
 
+#include <complex>
+extern "C" {
+#include "blas.h"
+}
+
 using namespace std;
 
 template<class Interface>
@@ -39,7 +44,7 @@ public :
   BTL_DONT_INLINE Action_parallel_matrix_vector_product( int size ):_size(size)
   {
     MESSAGE("Action_parallel_matrix_vector_product Ctor");
-    const int iZERO = 0, iONE = 1;
+    int iZERO = 0, iONE = 1;
 
     GlobalRows = _size;
     GlobalCols = _size;
@@ -59,6 +64,14 @@ public :
       init_vector<null_function>(Global_y_stl, GlobalRows);
 
       // Compute YTest (?)
+      Test_y_stl.resize(GlobalRows);
+      double alpha = 1., beta = 0.;
+      char notrans = 'N';
+      dgemv_(&notrans, &GlobalRows, &GlobalCols,
+          &alpha, &Global_A_stl[0], &GlobalRows,
+                  &Global_x_stl[0], &iONE,
+           &beta,   &Test_y_stl[0], &iONE
+      );
     }
 
     Interface::scatter_matrix(Global_A_stl, Local_A_stl, GlobalRows, GlobalCols, BlockRows, BlockCols,  LocalRows, LocalCols);
@@ -87,7 +100,6 @@ public :
   }
 
   // invalidate copy ctor
-
   Action_parallel_matrix_vector_product( const  Action_parallel_matrix_vector_product & )
   {
     INFOS("illegal call to Action_parallel_matrix_vector_product Copy Ctor");
@@ -95,7 +107,6 @@ public :
   }
 
   // Dtor
-
   BTL_DONT_INLINE ~Action_parallel_matrix_vector_product( void ){
 
     MESSAGE("Action_parallel_matrix_vector_product Dtor");
@@ -113,7 +124,6 @@ public :
   }
 
   // action name
-
   static inline std::string name( void )
   {
     return "parallel_matrix_vector_" + Interface::name();
@@ -138,10 +148,21 @@ public :
   }
 
   BTL_DONT_INLINE void check_result( void ){
+    int iONE = 1;
+    double dmONE = -1.;
+    int GlobalYCols;
+    Interface::vector_to_stl(Local_y, Local_y_stl);
+
+    Interface::gather_matrix(Global_y_stl, Local_y_stl, GlobalRows, GlobalYCols, BlockRows, BlockCols, LocalYRows, LocalYCols);
 
     // calculation check
+    if (iamroot) {
+      daxpy_(&GlobalRows, &dmONE, &Global_y_stl[0], &iONE, &Test_y_stl[0], &iONE);
+      double nrm = dnrm2_(&GlobalRows, &Test_y_stl[0], &iONE);
 
-    // TODO
+      if (nrm > 1e-5)
+        std::cerr << "Error: " << nrm << std::endl;
+    }
 
   }
 
@@ -150,6 +171,7 @@ private :
   typename Interface::stl_matrix Global_A_stl;
   typename Interface::stl_vector Global_x_stl;
   typename Interface::stl_vector Global_y_stl;
+  typename Interface::stl_vector Test_y_stl;
 
   typename Interface::stl_matrix Local_A_stl;
   typename Interface::stl_vector Local_x_stl;
diff --git a/btl/libs/BLACS/blacsinit.hh b/btl/libs/BLACS/blacsinit.hh
new file mode 100644
index 0000000..6cfeee3
--- /dev/null
+++ b/btl/libs/BLACS/blacsinit.hh
@@ -0,0 +1,21 @@
+#ifndef BLACSINIT_HH
+#define BLACSINIT_HH
+
+#include <cmath>
+
+bool blacsinit(int *argc, char ***argv)
+{
+  int context, myid, numproc, prows, pcols, iZERO = 0;
+  MPI_Init(argc, argv);
+  blacs_pinfo_(&myid, &numproc);
+  blacs_get_(&iZERO, &iZERO, &context);
+
+  int p = static_cast<double>(std::sqrt(static_cast<double>(numproc)) + 1.);
+  while (numproc % p) --p;
+  prows = p; pcols = numproc/p;
+
+  blacs_gridinit_(&context, "Row-major", &prows, &pcols);
+  return (myid == 0);
+}
+
+#endif /* BLACSINIT_HH */
diff --git a/btl/libs/BLACS/gather_impl.h b/btl/libs/BLACS/gather_impl.h
index 7683454..92a52f3 100644
--- a/btl/libs/BLACS/gather_impl.h
+++ b/btl/libs/BLACS/gather_impl.h
@@ -55,9 +55,7 @@ inline void gather(
     minfo[0] = LocalRows; minfo[1] = LocalCols;
     igsum2d_(&context, "Col", " ", &iONE, &iONE, minfo, &iONE, &imONE, &imONE);
     igsum2d_(&context, "Row", " ", &iONE, &iONE, minfo+1, &iONE, &imONE, &imONE);
-    if(iamroot) {
-    	GlobalRows = minfo[0]; GlobalCols = minfo[1];
-    }
+    GlobalRows = minfo[0]; GlobalCols = minfo[1];
 
     
     /* Reserve space on root */
diff --git a/btl/libs/PBLAS/main.cpp b/btl/libs/PBLAS/main.cpp
index bca245c..33a4f96 100644
--- a/btl/libs/PBLAS/main.cpp
+++ b/btl/libs/PBLAS/main.cpp
@@ -5,6 +5,7 @@
 #include <iostream>
 //using namespace std;
 
+#include "blacsinit.hh"
 #include "pblas_interface.hh"
 #include "action_parallel_matrix_vector_product.hh"
 
@@ -14,18 +15,16 @@ BTL_MAIN;
 
 int main(int argc, char **argv)
 {
-  MPI_Init(&argc, &argv);
+  bool iamroot = blacsinit(&argc, &argv);
 
-  int myid, numproc, context, procrows = 2, proccols = 2, iZERO = 0;
-  blacs_pinfo_(&myid, &numproc);
-  blacs_get_(&iZERO, &iZERO, &context);
-  blacs_gridinit_(&context, "Row-major", &procrows, &proccols);
-  bool iamroot = (myid == 0);
+  distr_bench<Action_parallel_matrix_vector_product<pblas_interface<double> > >(10,MAX_MV,NB_POINT,!iamroot);
+//  Action_parallel_matrix_vector_product<pblas_interface<double> > action(3000);
+//  action.initialize();
+//  action.calculate();
+//  action.check_result();
 
-  distr_bench<Action_parallel_matrix_vector_product<pblas_interface<REAL_TYPE> > >(10,MAX_MV,NB_POINT,!iamroot);
-
-//  blacs_exit_(&iZERO);
-  MPI_Finalize();
+  int iZERO = 0;
+  blacs_exit_(&iZERO);
   return 0;
 }