[gentoo-commits] proj/auto-numerical-bench:unstable commit in: /, btl/actions/, btl/libs/PBLAS/, btl/generic_bench/init/

[Andrea Arteaga]

commit:     073daf866484570b163359dc3e50a1c9f57acfa7
Author:     spiros <andyspiros <AT> gmail <DOT> com>
AuthorDate: Fri Jul 22 22:51:54 2011 +0000
Commit:     Andrea Arteaga <andyspiros <AT> gmail <DOT> com>
CommitDate: Fri Jul 22 22:51:54 2011 +0000
URL:        http://git.overlays.gentoo.org/gitweb/?p=proj/auto-numerical-bench.git;a=commit;h=073daf86

Added cholesky. Solved some issues.

---
 ...el_lu_decomp.hh => action_parallel_cholesky.hh} |   62 +++++++++----------
 btl/actions/action_parallel_lu_decomp.hh           |   36 ++++++------
 btl/generic_bench/init/init_function.hh            |    2 +-
 btl/libs/PBLAS/main.cpp                            |   42 +++++++++++---
 btl/libs/PBLAS/pblas.h                             |    5 ++
 btl/libs/PBLAS/pblas_interface_impl.hh             |   10 +++
 pblas.py                                           |    4 +-
 7 files changed, 99 insertions(+), 62 deletions(-)

diff --git a/btl/actions/action_parallel_lu_decomp.hh b/btl/actions/action_parallel_cholesky.hh
similarity index 51%
copy from btl/actions/action_parallel_lu_decomp.hh
copy to btl/actions/action_parallel_cholesky.hh
index 4882a21..f89eb98 100644
--- a/btl/actions/action_parallel_lu_decomp.hh
+++ b/btl/actions/action_parallel_cholesky.hh
@@ -1,5 +1,5 @@
-#ifndef ACTION_PARALLEL_LU_DECOMP_HH_
-#define ACTION_PARALLEL_LU_DECOMP_HH_
+#ifndef ACTION_PARALLEL_CHOLESKY_HH_
+#define ACTION_PARALLEL_CHOLESKY_HH_
 
 #include "utilities.h"
 #include "init/init_function.hh"
@@ -11,14 +11,15 @@
 #include <string>
 
 template<class Interface>
-class Action_parallel_lu_decomp {
+class Action_parallel_cholesky {
+  typedef lapack_interface<typename Interface::real_type> LapackInterface;
 
 public :
 
   // Constructor
-  BTL_DONT_INLINE Action_parallel_lu_decomp( int size ) : _size(size)
+  BTL_DONT_INLINE Action_parallel_cholesky( int size ) : _size(size)
   {
-    MESSAGE("Action_parallel_lu_decomp Ctor");
+    MESSAGE("Action_parallel_cholesky Ctor");
 
     int myid, procnum;
     blacs_pinfo_(&myid, &procnum);
@@ -26,7 +27,16 @@ public :
 
     // STL matrix and vector initialization
     if (iamroot) {
-      init_vector<pseudo_random>(Global_A_stl, size*size);
+        typename LapackInterface::stl_matrix temp_stl;
+        init_matrix_symm<pseudo_random>(temp_stl, size);
+        Global_A_stl.reserve(size*size);
+        const double add = 5000./size;
+        for (int r = 0; r < size; ++r)
+          for (int c = 0; c < size; ++c)
+            if (r==c)
+              Global_A_stl.push_back((std::abs(temp_stl[r][c])+add)*size);
+            else
+              Global_A_stl.push_back(temp_stl[r][c]);
     }
 
     Interface::scatter_matrix(Global_A_stl, Local_A_stl, desc, size, size, 64, 64);
@@ -37,21 +47,25 @@ public :
     Interface::matrix_from_stl(Local_A_ref, Local_A_stl);
     Interface::matrix_from_stl(Local_A    , Local_A_stl);
 
-    _cost = 2.0*size*size*size/3.0 + static_cast<double>(size*size);
+    _cost = 0;
+    for (int j=0; j<_size; ++j) {
+      double r = std::max(_size - j -1,0);
+      _cost += 2*(r*j+r+j);
+    }
   }
 
 
   // Invalidate copy constructor
-  Action_parallel_lu_decomp(const Action_parallel_lu_decomp&)
+  Action_parallel_cholesky(const Action_parallel_cholesky&)
   {
-    INFOS("illegal call to Action_parallel_lu_decomp copy constructor");
+    INFOS("illegal call to Action_parallel_cholesky copy constructor");
     exit(1);
   }
 
   // Destructor
-  ~Action_parallel_lu_decomp()
+  ~Action_parallel_cholesky()
   {
-    MESSAGE("Action_parallel_lu_decomp destructor");
+    MESSAGE("Action_parallel_cholesky destructor");
 
     // Deallocation
     Interface::free_matrix(Local_A_ref, Local_A_stl.size());
@@ -61,7 +75,7 @@ public :
   // Action name
   static inline std::string name()
   {
-    return "lu_decomp_" + Interface::name();
+    return "cholesky_" + Interface::name();
   }
 
   double nb_op_base()
@@ -76,31 +90,14 @@ public :
 
   BTL_DONT_INLINE void calculate()
   {
-    Interface::parallel_lu_decomp(Local_A, desc);
+    Interface::parallel_cholesky(Local_A, desc);
   }
 
   BTL_DONT_INLINE void check_result()
   {
-    /* Gather */
-    typename Interface::stl_matrix Global_Y;
-    typename Interface::stl_matrix Local_Y(Local_A_stl.size());
-    Interface::matrix_to_stl(Local_A, Local_Y);
-    Interface::gather_matrix(Global_Y, Local_Y, desc);
-
-    if (iamroot) {
-
-      typename Interface::gene_matrix A;
-      Interface::matrix_from_stl(A, Global_A_stl);
-      lapack_interface<typename Interface::real_type>::lu_decomp(&Global_A_stl[0], A, _size);
-      typename Interface::stl_vector correct(A, A+_size*_size);
-      typename Interface::real_type error = STL_interface<typename Interface::real_type>::norm_diff(Global_Y, correct);
-      if (error > 10e-5)
-        cerr << " { !! error : " << error << " !! } ";
-
-      Interface::free_matrix(A, _size*_size);
-    }
   }
 
+
 private:
   int _size, desc[9], LocalRows, LocalCols;
   double _cost;
@@ -112,5 +109,4 @@ private:
   typename Interface::gene_matrix Local_A;
 };
 
-
-#endif /* ACTION_PARALLEL_LU_DECOMP_HH_ */
+#endif /* ACTION_PARALLEL_CHOLESKY_HH_ */

diff --git a/btl/actions/action_parallel_lu_decomp.hh b/btl/actions/action_parallel_lu_decomp.hh
index 4882a21..18b4ac7 100644
--- a/btl/actions/action_parallel_lu_decomp.hh
+++ b/btl/actions/action_parallel_lu_decomp.hh
@@ -81,24 +81,24 @@ public :
 
   BTL_DONT_INLINE void check_result()
   {
-    /* Gather */
-    typename Interface::stl_matrix Global_Y;
-    typename Interface::stl_matrix Local_Y(Local_A_stl.size());
-    Interface::matrix_to_stl(Local_A, Local_Y);
-    Interface::gather_matrix(Global_Y, Local_Y, desc);
-
-    if (iamroot) {
-
-      typename Interface::gene_matrix A;
-      Interface::matrix_from_stl(A, Global_A_stl);
-      lapack_interface<typename Interface::real_type>::lu_decomp(&Global_A_stl[0], A, _size);
-      typename Interface::stl_vector correct(A, A+_size*_size);
-      typename Interface::real_type error = STL_interface<typename Interface::real_type>::norm_diff(Global_Y, correct);
-      if (error > 10e-5)
-        cerr << " { !! error : " << error << " !! } ";
-
-      Interface::free_matrix(A, _size*_size);
-    }
+//    /* Gather */
+//    typename Interface::stl_matrix Global_Y;
+//    typename Interface::stl_matrix Local_Y(Local_A_stl.size());
+//    Interface::matrix_to_stl(Local_A, Local_Y);
+//    Interface::gather_matrix(Global_Y, Local_Y, desc);
+//
+//    if (iamroot) {
+//
+//      typename Interface::gene_matrix A;
+//      Interface::matrix_from_stl(A, Global_A_stl);
+//      lapack_interface<typename Interface::real_type>::lu_decomp(&Global_A_stl[0], A, _size);
+//      typename Interface::stl_vector correct(A, A+_size*_size);
+//      typename Interface::real_type error = STL_interface<typename Interface::real_type>::norm_diff(Global_Y, correct);
+//      if (error > 10e-5)
+//        cerr << " { !! error : " << error << " !! } ";
+//
+//      Interface::free_matrix(A, _size*_size);
+//    }
   }
 
 private:

diff --git a/btl/generic_bench/init/init_function.hh b/btl/generic_bench/init/init_function.hh
index 5401673..7b3bdba 100644
--- a/btl/generic_bench/init/init_function.hh
+++ b/btl/generic_bench/init/init_function.hh
@@ -32,7 +32,7 @@ double simple_function(int index_i, int index_j)
 
 double pseudo_random(int index)
 {
-  return static_cast<int>((std::rand()/double(RAND_MAX)-.5)*20);
+  return std::rand()/double(RAND_MAX);
 }
 
 double pseudo_random(int index_i, int index_j)

diff --git a/btl/libs/PBLAS/main.cpp b/btl/libs/PBLAS/main.cpp
index a2aae2a..e7b636b 100644
--- a/btl/libs/PBLAS/main.cpp
+++ b/btl/libs/PBLAS/main.cpp
@@ -9,9 +9,11 @@
 
 #include "pblas_interface.hh"
 #include "blacsinit.hh"
+
+#include "action_parallel_axpy.hh"
 #include "action_parallel_matrix_vector_product.hh"
 #include "action_parallel_lu_decomp.hh"
-#include "action_parallel_axpy.hh"
+#include "action_parallel_cholesky.hh"
 
 #include <string>
 
@@ -21,13 +23,37 @@ int main(int argc, char **argv)
 {
   bool iamroot = blacsinit(&argc, &argv);
 
-//  distr_bench<Action_parallel_matrix_vector_product<pblas_interface<double> > >(10,MAX_MV,NB_POINT,!iamroot);
-//  distr_bench<Action_parallel_axpy<pblas_interface<REAL_TYPE> > >(10,MAX_MV,NB_POINT,!iamroot);
-  distr_bench<Action_parallel_lu_decomp<pblas_interface<REAL_TYPE> > >(10,MAX_MV,NB_POINT,!iamroot);
-//  Action_parallel_lu_decomp<pblas_interface<REAL_TYPE> > action(8);
-//  action.initialize();
-//  action.calculate();
-//  action.check_result();
+  bool
+  general_solve=false, least_squares=false, lu_decomp=false, cholesky=false,
+  symm_ev=false
+  ;
+
+
+  for (int i = 1; i < argc; ++i) {
+          std::string arg = argv[i];
+          if (arg == "general_solve") general_solve = true;
+          else if (arg == "least_squares") least_squares = true;
+          else if (arg == "lu_decomp") lu_decomp = true;
+          else if (arg == "cholesky") cholesky = true;
+          else if (arg == "symm_ev") symm_ev = true;
+  }
+
+
+//  if (general_solve)
+//  distr_bench<Action_general_solve<lapack_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT, !iamroot);
+
+//  if (least_squares)
+//  distr_bench<Action_least_squares<lapack_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT, !iamroot);
+
+  if (lu_decomp)
+  distr_bench<Action_parallel_lu_decomp<pblas_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT, !iamroot);
+
+  if (cholesky)
+    distr_bench<Action_parallel_cholesky<pblas_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT, !iamroot);
+
+//  if (symm_ev)
+//  distr_bench<Action_symm_ev<lapack_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT, !iamroot);
+
 
   int iZERO = 0;
   blacs_exit_(&iZERO);

diff --git a/btl/libs/PBLAS/pblas.h b/btl/libs/PBLAS/pblas.h
index 2b5860e..973b91c 100644
--- a/btl/libs/PBLAS/pblas.h
+++ b/btl/libs/PBLAS/pblas.h
@@ -46,9 +46,14 @@ extern "C" {
   /*************
    * Scalapack *
    *************/
+  // lu_decomp
   void psgetrf_(const int*, const int*,  float*, const int*, const int*, const int*, int*, int*);
   void pdgetrf_(const int*, const int*, double*, const int*, const int*, const int*, int*, int*);
 
+  // cholesky
+  void pspotrf_(const char*, const int*,  float*, const int*, const int*, const int*, int*);
+  void pdpotrf_(const char*, const int*, double*, const int*, const int*, const int*, int*);
+
 
 #ifdef __cplusplus
 }

diff --git a/btl/libs/PBLAS/pblas_interface_impl.hh b/btl/libs/PBLAS/pblas_interface_impl.hh
index c36c249..1dbf3b9 100644
--- a/btl/libs/PBLAS/pblas_interface_impl.hh
+++ b/btl/libs/PBLAS/pblas_interface_impl.hh
@@ -54,5 +54,15 @@ public:
     PBLAS_FUNC(getrf)(&GlobalRows, &GlobalCols, X, &iONE, &iONE, desc,
         &ipiv[0], &info);
   }
+
+  static inline void parallel_cholesky(gene_matrix& X, const int* desc)
+  {
+    const int N = desc[2], iONE = 1;
+    const char UPLO = 'U';
+    int info;
+    PBLAS_FUNC(potrf)(&UPLO, &N, X, &iONE, &iONE, desc, &info);
+//    if (info != 0)
+//      cerr << " { cholesky error : " << info << " } ";
+  }
 };
 

diff --git a/pblas.py b/pblas.py
index 6f0f6cd..9cd087e 100644
--- a/pblas.py
+++ b/pblas.py
@@ -5,7 +5,7 @@ numproc = 4
 class Module(btlbase.BTLBase):
     def _initialize(self):
         self.libname = "scalapack"
-        self.avail = ['axpy', 'matrix_vector', 'lu_decomp']
+        self.avail = ['axpy', 'matrix_vector', 'lu_decomp', 'cholesky']
     
     def _parse_args(self, args):     
         # Parse arguments
@@ -61,7 +61,7 @@ class PBLASTest(btlbase.BTLTest):
     
     @staticmethod
     def _btl_includes():
-        return ["libs/BLAS", "libs/BLACS", "libs/PBLAS"]
+        return ["libs/BLAS", "libs/BLACS", "libs/PBLAS", "libs/STL"]
     
     def _btl_defines(self):
         return ["PBLASNAME="+self.libname]
\ No newline at end of file