* [gentoo-commits] proj/auto-numerical-bench:unstable commit in: /, btl/actions/, btl/libs/PBLAS/
2011-08-02 18:45 Andrea Arteaga
@ 2011-07-23 11:46 ` Andrea Arteaga
0 siblings, 0 replies; 4+ messages in thread
From: Andrea Arteaga @ 2011-07-23 11:46 UTC (permalink / raw
To: gentoo-commits
commit: c51172bc46f4b95af6282d8782e4b145911c7afe
Author: spiros <andyspiros <AT> gmail <DOT> com>
AuthorDate: Sat Jul 23 11:45:29 2011 +0000
Commit: Andrea Arteaga <andyspiros <AT> gmail <DOT> com>
CommitDate: Sat Jul 23 11:45:29 2011 +0000
URL: http://git.overlays.gentoo.org/gitweb/?p=proj/auto-numerical-bench.git;a=commit;h=c51172bc
Added working QR decomposition; added working symm_ev (but some negative
MFlops).
---
btl/actions/action_parallel_cholesky.hh | 3 +-
btl/actions/action_parallel_lu_decomp.hh | 3 +-
...el_cholesky.hh => action_parallel_qr_decomp.hh} | 48 +++-----
btl/actions/action_parallel_symm_ev.hh | 121 ++++++++++++++++++++
btl/libs/PBLAS/main.cpp | 14 ++-
btl/libs/PBLAS/pblas.h | 8 ++
btl/libs/PBLAS/pblas_interface_impl.hh | 54 +++++++++-
pblas.py | 3 +-
8 files changed, 213 insertions(+), 41 deletions(-)
diff --git a/btl/actions/action_parallel_cholesky.hh b/btl/actions/action_parallel_cholesky.hh
index f89eb98..05ef3ef 100644
--- a/btl/actions/action_parallel_cholesky.hh
+++ b/btl/actions/action_parallel_cholesky.hh
@@ -39,7 +39,8 @@ public :
Global_A_stl.push_back(temp_stl[r][c]);
}
- Interface::scatter_matrix(Global_A_stl, Local_A_stl, desc, size, size, 64, 64);
+ const int blocksize = std::max(std::min(size/4, 64), 2);
+ Interface::scatter_matrix(Global_A_stl, Local_A_stl, desc, size, size, blocksize, blocksize);
LocalRows = desc[8];
LocalCols = Local_A_stl.size()/desc[8];
diff --git a/btl/actions/action_parallel_lu_decomp.hh b/btl/actions/action_parallel_lu_decomp.hh
index 18b4ac7..d3dc620 100644
--- a/btl/actions/action_parallel_lu_decomp.hh
+++ b/btl/actions/action_parallel_lu_decomp.hh
@@ -29,7 +29,8 @@ public :
init_vector<pseudo_random>(Global_A_stl, size*size);
}
- Interface::scatter_matrix(Global_A_stl, Local_A_stl, desc, size, size, 64, 64);
+ const int blocksize = std::max(std::min(size/4, 64), 2);
+ Interface::scatter_matrix(Global_A_stl, Local_A_stl, desc, size, size, blocksize, blocksize);
LocalRows = desc[8];
LocalCols = Local_A_stl.size()/desc[8];
diff --git a/btl/actions/action_parallel_cholesky.hh b/btl/actions/action_parallel_qr_decomp.hh
similarity index 55%
copy from btl/actions/action_parallel_cholesky.hh
copy to btl/actions/action_parallel_qr_decomp.hh
index f89eb98..a41414c 100644
--- a/btl/actions/action_parallel_cholesky.hh
+++ b/btl/actions/action_parallel_qr_decomp.hh
@@ -1,5 +1,5 @@
-#ifndef ACTION_PARALLEL_CHOLESKY_HH_
-#define ACTION_PARALLEL_CHOLESKY_HH_
+#ifndef ACTION_PARALLEL_QR_DECOMP_HH_
+#define ACTION_PARALLEL_QR_DECOMP_HH_
#include "utilities.h"
#include "init/init_function.hh"
@@ -9,17 +9,17 @@
#include "STL_interface.hh"
#include <string>
+#include <algorithm>
template<class Interface>
-class Action_parallel_cholesky {
- typedef lapack_interface<typename Interface::real_type> LapackInterface;
+class Action_parallel_qr_decomp {
public :
// Constructor
- BTL_DONT_INLINE Action_parallel_cholesky( int size ) : _size(size)
+ BTL_DONT_INLINE Action_parallel_qr_decomp( int size ) : _size(size)
{
- MESSAGE("Action_parallel_cholesky Ctor");
+ MESSAGE("Action_parallel_qr_decomp Ctor");
int myid, procnum;
blacs_pinfo_(&myid, &procnum);
@@ -27,19 +27,11 @@ public :
// STL matrix and vector initialization
if (iamroot) {
- 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]);
+ init_vector<pseudo_random>(Global_A_stl, size*size);
}
- Interface::scatter_matrix(Global_A_stl, Local_A_stl, desc, size, size, 64, 64);
+ const int blocksize = std::max(std::min(size/4, 64), 2);
+ Interface::scatter_matrix(Global_A_stl, Local_A_stl, desc, size, size, blocksize, blocksize);
LocalRows = desc[8];
LocalCols = Local_A_stl.size()/desc[8];
@@ -47,25 +39,21 @@ public :
Interface::matrix_from_stl(Local_A_ref, Local_A_stl);
Interface::matrix_from_stl(Local_A , Local_A_stl);
- _cost = 0;
- for (int j=0; j<_size; ++j) {
- double r = std::max(_size - j -1,0);
- _cost += 2*(r*j+r+j);
- }
+ _cost = 2.0*size*size*size;
}
// Invalidate copy constructor
- Action_parallel_cholesky(const Action_parallel_cholesky&)
+ Action_parallel_qr_decomp(const Action_parallel_qr_decomp&)
{
- INFOS("illegal call to Action_parallel_cholesky copy constructor");
+ INFOS("illegal call to Action_parallel_qr_decomp copy constructor");
exit(1);
}
// Destructor
- ~Action_parallel_cholesky()
+ ~Action_parallel_qr_decomp()
{
- MESSAGE("Action_parallel_cholesky destructor");
+ MESSAGE("Action_parallel_qr_decomp destructor");
// Deallocation
Interface::free_matrix(Local_A_ref, Local_A_stl.size());
@@ -75,7 +63,7 @@ public :
// Action name
static inline std::string name()
{
- return "cholesky_" + Interface::name();
+ return "qr_decomp_" + Interface::name();
}
double nb_op_base()
@@ -90,14 +78,13 @@ public :
BTL_DONT_INLINE void calculate()
{
- Interface::parallel_cholesky(Local_A, desc);
+ Interface::parallel_qr_decomp(Local_A, desc);
}
BTL_DONT_INLINE void check_result()
{
}
-
private:
int _size, desc[9], LocalRows, LocalCols;
double _cost;
@@ -109,4 +96,5 @@ private:
typename Interface::gene_matrix Local_A;
};
-#endif /* ACTION_PARALLEL_CHOLESKY_HH_ */
+
+#endif /* ACTION_PARALLEL_QR_DECOMP_HH_ */
diff --git a/btl/actions/action_parallel_symm_ev.hh b/btl/actions/action_parallel_symm_ev.hh
new file mode 100644
index 0000000..f0af0e3
--- /dev/null
+++ b/btl/actions/action_parallel_symm_ev.hh
@@ -0,0 +1,121 @@
+#ifndef ACTION_PARALLEL_SYMM_EV_HH_
+#define ACTION_PARALLEL_SYMM_EV_HH_
+
+#include "utilities.h"
+#include "init/init_function.hh"
+#include "init/init_vector.hh"
+
+#include "lapack_interface.hh"
+#include "STL_interface.hh"
+
+#include <string>
+
+template<class Interface>
+class Action_parallel_symm_ev {
+
+public :
+
+ // Constructor
+ BTL_DONT_INLINE Action_parallel_symm_ev( int size ) : _size(size)
+ {
+ MESSAGE("Action_parallel_symm_ev constructor");
+
+ int myid, procnum;
+ blacs_pinfo_(&myid, &procnum);
+ iamroot = (myid == 0);
+
+ // STL matrix and vector initialization
+ if (iamroot) {
+ init_vector<pseudo_random>(Global_A_stl, size*size);
+ init_vector<null_function>(Global_Z_stl, size*size);
+ }
+ init_vector<null_function>(Local_w_stl, size);
+
+ const int blocksize = std::max(std::min(size/4, 64), 2);
+ Interface::scatter_matrix(Global_A_stl, Local_A_stl, descA, size, size, blocksize, blocksize);
+ Interface::scatter_matrix(Global_Z_stl, Local_Z_stl, descZ, size, size, blocksize, blocksize);
+ LocalRows = descA[8];
+ LocalCols = Local_A_stl.size()/descA[8];
+
+ // Generic local matrix and vectors initialization
+ Interface::matrix_from_stl(Local_A_ref, Local_A_stl);
+ Interface::matrix_from_stl(Local_A , Local_A_stl);
+ Interface::matrix_from_stl(Local_Z_ref, Local_Z_stl);
+ Interface::matrix_from_stl(Local_Z , Local_Z_stl);
+ Interface::vector_from_stl(Local_w , Local_w_stl);
+ Interface::vector_from_stl(Local_w_ref, Local_w_stl);
+
+ _cost = size*size*size;
+ }
+
+
+ // Invalidate copy constructor
+ Action_parallel_symm_ev(const Action_parallel_symm_ev&)
+ {
+ INFOS("illegal call to Action_parallel_symm_ev copy constructor");
+ exit(1);
+ }
+
+ // Destructor
+ ~Action_parallel_symm_ev()
+ {
+ MESSAGE("Action_parallel_lu_decomp destructor");
+
+ // Deallocation
+ Interface::free_matrix(Local_A_ref, Local_A_stl.size());
+ Interface::free_matrix(Local_A , Local_A_stl.size());
+ Interface::free_matrix(Local_Z_ref, Local_Z_stl.size());
+ Interface::free_matrix(Local_Z , Local_Z_stl.size());
+ Interface::free_vector(Local_w_ref);
+ Interface::free_vector(Local_w );
+ }
+
+ // Action name
+ static inline std::string name()
+ {
+ return "symm_ev_" + Interface::name();
+ }
+
+ double nb_op_base()
+ {
+ return _cost;
+ }
+
+ BTL_DONT_INLINE void initialize()
+ {
+ Interface::copy_matrix(Local_A_ref, Local_A, Local_A_stl.size());
+ Interface::copy_matrix(Local_Z_ref, Local_Z, Local_Z_stl.size());
+ Interface::copy_vector(Local_w_ref, Local_w, Local_w_stl.size());
+ }
+
+ BTL_DONT_INLINE void calculate()
+ {
+ Interface::parallel_symm_ev(Local_A, descA, Local_w, Local_Z, descZ);
+ }
+
+ BTL_DONT_INLINE void check_result()
+ {
+ }
+
+private:
+ int _size, descA[9], descZ[9], LocalRows, LocalCols;
+ double _cost;
+ bool iamroot;
+
+ typename Interface::stl_matrix Global_A_stl;
+ typename Interface::stl_matrix Local_A_stl;
+ typename Interface::gene_matrix Local_A_ref;
+ typename Interface::gene_matrix Local_A;
+
+ typename Interface::stl_matrix Global_Z_stl;
+ typename Interface::stl_matrix Local_Z_stl;
+ typename Interface::gene_matrix Local_Z_ref;
+ typename Interface::gene_matrix Local_Z;
+
+ typename Interface::stl_vector Local_w_stl;
+ typename Interface::gene_vector Local_w_ref;
+ typename Interface::gene_vector Local_w;
+};
+
+
+#endif /* ACTION_PARALLEL_LU_DECOMP_HH_ */
diff --git a/btl/libs/PBLAS/main.cpp b/btl/libs/PBLAS/main.cpp
index e7b636b..c209afe 100644
--- a/btl/libs/PBLAS/main.cpp
+++ b/btl/libs/PBLAS/main.cpp
@@ -14,6 +14,8 @@
#include "action_parallel_matrix_vector_product.hh"
#include "action_parallel_lu_decomp.hh"
#include "action_parallel_cholesky.hh"
+#include "action_parallel_qr_decomp.hh"
+#include "action_parallel_symm_ev.hh"
#include <string>
@@ -24,7 +26,7 @@ int main(int argc, char **argv)
bool iamroot = blacsinit(&argc, &argv);
bool
- general_solve=false, least_squares=false, lu_decomp=false, cholesky=false,
+ general_solve=false, qr_decomp=false, lu_decomp=false, cholesky=false,
symm_ev=false
;
@@ -32,7 +34,7 @@ int main(int argc, char **argv)
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 == "qr_decomp") qr_decomp = true;
else if (arg == "lu_decomp") lu_decomp = true;
else if (arg == "cholesky") cholesky = true;
else if (arg == "symm_ev") symm_ev = true;
@@ -42,8 +44,8 @@ int main(int argc, char **argv)
// 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 (qr_decomp)
+ distr_bench<Action_parallel_qr_decomp<pblas_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);
@@ -51,8 +53,8 @@ int main(int argc, char **argv)
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);
+ if (symm_ev)
+ distr_bench<Action_parallel_symm_ev<pblas_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT, !iamroot);
int iZERO = 0;
diff --git a/btl/libs/PBLAS/pblas.h b/btl/libs/PBLAS/pblas.h
index 973b91c..a6cbeb2 100644
--- a/btl/libs/PBLAS/pblas.h
+++ b/btl/libs/PBLAS/pblas.h
@@ -54,6 +54,14 @@ extern "C" {
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*);
+ // qr_decomp
+ void psgeqpf_(const int*, const int*, float*, const int*, const int*, const int*, int*, float*, float*, const int*, int*);
+ void pdgeqpf_(const int*, const int*, double*, const int*, const int*, const int*, int*, double*, double*, const int*, int*);
+
+ // symm_ev
+ void pssyevd_(const char*, const char*, const int*, float*, const int*, const int*, const int*, float*, float*, const int*, const int*, const int*, float*, const int*, int*, const int*, int*);
+ void pdsyevd_(const char*, const char*, const int*, double*, const int*, const int*, const int*, double*, double*, const int*, const int*, const int*, double*, const int*, 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 1dbf3b9..4522946 100644
--- a/btl/libs/PBLAS/pblas_interface_impl.hh
+++ b/btl/libs/PBLAS/pblas_interface_impl.hh
@@ -61,8 +61,58 @@ public:
const char UPLO = 'U';
int info;
PBLAS_FUNC(potrf)(&UPLO, &N, X, &iONE, &iONE, desc, &info);
+ if (info != 0)
+ cerr << " { cholesky error : " << info << " } ";
+ }
+
+ static inline void parallel_qr_decomp(gene_matrix& X, const int* desc)
+ {
+ const int GlobalRows = desc[2], GlobalCols = desc[3],
+ BlockRows = desc[4], BlockCols = desc[5],
+ ctxt = desc[1];
+
+ int myrow, mycol, nprow, npcol, lwork;
+ SCALAR lworkd;
+ blacs_gridinfo_(&ctxt, &nprow, &npcol, &myrow, &mycol);
+
+ const int iONE = 1, iZERO = 0, imONE = -1,
+ ipivdim = numroc_(&GlobalCols, &BlockCols, &mycol, &iZERO, &npcol);
+ int info;
+ std::vector<int> ipiv(ipivdim);
+ std::vector<SCALAR> tau(ipivdim);
+
+ // Retrieve LWORK
+ PBLAS_FUNC(geqpf)(&GlobalRows, &GlobalCols, X, &iONE, &iONE, desc, &ipiv[0], &tau[0], &lworkd, &imONE, &info);
+ lwork = static_cast<int>(lworkd);
+// if (info != 0)
+// cerr << " { qr_decomp lwork error } ";
+
+ std::vector<SCALAR> work(lwork);
+ PBLAS_FUNC(geqpf)(&GlobalRows, &GlobalCols, X, &iONE, &iONE, desc, &ipiv[0], &tau[0], &work[0], &lwork, &info);
// if (info != 0)
-// cerr << " { cholesky error : " << info << " } ";
+// cerr << " { qr_decomp computation error } ";
}
-};
+ static inline void parallel_symm_ev(gene_matrix& A, const int* descA, gene_vector& w, gene_matrix& Z, const int* descZ)
+ {
+ const char jobz = 'V', uplo = 'u';
+ const int N = descA[2], iONE = 1, iZERO = 0, imONE = -1;
+ std::vector<SCALAR> work;
+ std::vector<int> iwork;
+ int lwork, liwork, info;
+ SCALAR lworkd;
+
+ // Retrieve l(i)work
+ PBLAS_FUNC(syevd)(&jobz, &uplo, &N, A, &iONE, &iONE, descA, w,
+ Z, &iONE, &iONE, descZ, &lworkd, &imONE, &liwork, &imONE, &info);
+ lwork = static_cast<int>(lworkd);
+ work.resize(lwork); iwork.resize(liwork);
+// if (info != 0)
+// cerr << " { symm_ev l(i)work error } ";
+
+ PBLAS_FUNC(syevd)(&jobz, &uplo, &N, A, &iONE, &iONE, descA, w,
+ Z, &iONE, &iONE, descZ, &work[0], &lwork, &iwork[0], &liwork, &info);
+// if (info != 0)
+// cerr << " { symm_ev computation error } ";
+ }
+};
diff --git a/pblas.py b/pblas.py
index 9cd087e..792f343 100644
--- a/pblas.py
+++ b/pblas.py
@@ -5,7 +5,8 @@ numproc = 4
class Module(btlbase.BTLBase):
def _initialize(self):
self.libname = "scalapack"
- self.avail = ['axpy', 'matrix_vector', 'lu_decomp', 'cholesky']
+ self.avail = ['axpy', 'matrix_vector', 'lu_decomp', 'cholesky',
+ 'qr_decomp', 'symm_ev']
def _parse_args(self, args):
# Parse arguments
^ permalink raw reply related [flat|nested] 4+ messages in thread
* [gentoo-commits] proj/auto-numerical-bench:unstable commit in: /, btl/actions/, btl/libs/PBLAS/
@ 2011-07-23 22:59 Andrea Arteaga
2011-08-02 18:45 ` Andrea Arteaga
0 siblings, 1 reply; 4+ messages in thread
From: Andrea Arteaga @ 2011-07-23 22:59 UTC (permalink / raw
To: gentoo-commits
commit: 540fa4b69dcc294f9462bf7ab5dc852b79944ee8
Author: spiros <andyspiros <AT> gmail <DOT> com>
AuthorDate: Sat Jul 23 22:59:15 2011 +0000
Commit: Andrea Arteaga <andyspiros <AT> gmail <DOT> com>
CommitDate: Sat Jul 23 22:59:15 2011 +0000
URL: http://git.overlays.gentoo.org/gitweb/?p=proj/auto-numerical-bench.git;a=commit;h=540fa4b6
Added SVD decomposition. Work on eigenvalues action and mat-vec
multiply.
---
.../action_parallel_matrix_vector_product.hh | 52 +++-----
btl/actions/action_parallel_svd_decomp.hh | 134 ++++++++++++++++++++
btl/actions/action_parallel_symm_ev.hh | 4 +-
btl/libs/PBLAS/main.cpp | 39 ++++--
btl/libs/PBLAS/pblas.h | 4 +
btl/libs/PBLAS/pblas_interface_impl.hh | 28 ++++-
pblas.py | 2 +-
7 files changed, 215 insertions(+), 48 deletions(-)
diff --git a/btl/actions/action_parallel_matrix_vector_product.hh b/btl/actions/action_parallel_matrix_vector_product.hh
index 67e64bf..5920115 100644
--- a/btl/actions/action_parallel_matrix_vector_product.hh
+++ b/btl/actions/action_parallel_matrix_vector_product.hh
@@ -22,27 +22,21 @@ public :
MESSAGE("Action_parallel_matrix_vector_product Ctor");
int iZERO = 0, iONE = 1;
- GlobalRows = _size;
- GlobalCols = _size;
- BlockRows = 2;
- BlockCols= 2;
- LocalXCols = 1;
- LocalYCols = 1;
-
int myid, procnum;
blacs_pinfo_(&myid, &procnum);
iamroot = (myid == 0);
// STL matrix and vector initialization
if (iamroot) {
- init_vector<pseudo_random>(Global_A_stl, GlobalRows*GlobalCols);
- init_vector<pseudo_random>(Global_x_stl, GlobalCols);
- init_vector<null_function>(Global_y_stl, GlobalRows);
+ init_vector<pseudo_random>(Global_A_stl, size*size);
+ init_vector<pseudo_random>(Global_x_stl, size);
+ init_vector<null_function>(Global_y_stl, size);
}
- Interface::scatter_matrix(Global_A_stl, Local_A_stl, descA, GlobalRows, GlobalCols, BlockRows, BlockCols);
- Interface::scatter_matrix(Global_x_stl, Local_x_stl, descX, GlobalCols, 1, BlockRows, BlockCols);
- Interface::scatter_matrix(Global_y_stl, Local_y_stl, descY, GlobalRows, 1, BlockRows, BlockCols);
+ const int blocksize = std::max(std::min(size/4, 64), 2);
+ Interface::scatter_matrix(Global_A_stl, Local_A_stl, descA, size, size, blocksize, blocksize);
+ Interface::scatter_matrix(Global_x_stl, Local_x_stl, descX, size, 1, blocksize, blocksize);
+ Interface::scatter_matrix(Global_y_stl, Local_y_stl, descY, size, 1, blocksize, blocksize);
// generic local matrix and vectors initialization
@@ -57,22 +51,22 @@ 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");
+ INFOS("illegal call to Action_parallel_matrix_vector_product copy constructor");
exit(1);
}
// Dtor
- BTL_DONT_INLINE ~Action_parallel_matrix_vector_product( void ){
+ BTL_DONT_INLINE ~Action_parallel_matrix_vector_product(){
- MESSAGE("Action_parallel_matrix_vector_product Dtor");
+ MESSAGE("Action_parallel_matrix_vector_product destructor");
// deallocation
- Interface::free_matrix(Local_A_ref, GlobalRows*GlobalCols);
+ Interface::free_matrix(Local_A_ref, _size*_size);
Interface::free_vector(Local_x_ref);
Interface::free_vector(Local_y_ref);
- Interface::free_matrix(Local_A, GlobalRows*GlobalCols);
+ Interface::free_matrix(Local_A, _size*_size);
Interface::free_vector(Local_x);
Interface::free_vector(Local_y);
@@ -89,35 +83,32 @@ public :
}
BTL_DONT_INLINE void initialize( void ){
- Interface::copy_matrix(Local_A_ref,Local_A,LocalRows*LocalCols);
- Interface::copy_vector(Local_x_ref,Local_x,LocalXRows*LocalXCols);
- Interface::copy_vector(Local_y_ref,Local_y,LocalYRows*LocalYCols);
+ Interface::copy_matrix(Local_A_ref, Local_A, Local_A_stl.size());
+ Interface::copy_vector(Local_x_ref, Local_x, Local_x_stl.size());
+ Interface::copy_vector(Local_y_ref, Local_y, Local_y_stl.size());
}
BTL_DONT_INLINE void calculate( void ) {
- BTL_ASM_COMMENT("#begin matrix_vector_product");
- Interface::parallel_matrix_vector_product(GlobalRows, GlobalCols, Local_A, descA, Local_x, descX, Local_y, descY);
- BTL_ASM_COMMENT("end matrix_vector_product");
+ Interface::parallel_matrix_vector_product(_size, _size, Local_A, descA, Local_x, descX, Local_y, descY);
}
BTL_DONT_INLINE void check_result( void ){
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);
+ Interface::gather_matrix(Global_y_stl, Local_y_stl, descY);
// calculation check
if (iamroot) {
// Compute YTest
- Test_y_stl.resize(GlobalRows);
+ Test_y_stl.resize(_size);
STL_interface<typename Interface::real_type>::matrix_vector_product(Global_A_stl, Global_x_stl, Test_y_stl, _size);
- typename Interface::real_type error =
- STL_interface<typename Interface::real_type>::norm_diff(Global_y_stl, Test_y_stl);
+ typename Interface::real_type error = STL_interface<typename Interface::real_type>::norm_diff(Global_y_stl, Test_y_stl);
if (error > 1e-5)
- std::cerr << "Error: " << error << std::endl;
+ std::cerr << "Error: " << error << " ";
}
}
@@ -142,8 +133,7 @@ private :
typename Interface::gene_vector Local_y;
bool iamroot;
- int _size, GlobalRows, GlobalCols, LocalRows, LocalCols, BlockRows, BlockCols;
- int LocalXRows, LocalXCols, LocalYRows, LocalYCols;
+ int _size;
int descA[9], descX[9], descY[9];
};
diff --git a/btl/actions/action_parallel_svd_decomp.hh b/btl/actions/action_parallel_svd_decomp.hh
new file mode 100644
index 0000000..790ff6d
--- /dev/null
+++ b/btl/actions/action_parallel_svd_decomp.hh
@@ -0,0 +1,134 @@
+#ifndef ACTION_PARALLEL_SVD_DECOMP_HH_
+#define ACTION_PARALLEL_SVD_DECOMP_HH_
+
+#include "utilities.h"
+#include "init/init_function.hh"
+#include "init/init_vector.hh"
+
+#include "lapack_interface.hh"
+#include "STL_interface.hh"
+
+#include <string>
+#include <algorithm>
+
+template<class Interface>
+class Action_parallel_svd_decomp {
+
+public :
+
+ // Constructor
+ BTL_DONT_INLINE Action_parallel_svd_decomp( int size ) : _size(size)
+ {
+ MESSAGE("Action_parallel_svd_decomp Ctor");
+
+ int myid, procnum;
+ blacs_pinfo_(&myid, &procnum);
+ iamroot = (myid == 0);
+
+ // STL matrix and vector initialization
+ if (iamroot) {
+ init_vector<pseudo_random>(Global_A_stl, size*size);
+ init_vector<pseudo_random>(Global_U_stl, size*size);
+ init_vector<pseudo_random>(Global_V_stl, size*size);
+ }
+ init_vector<null_function>(Local_s_stl, size);
+
+ const int blocksize = std::max(std::min(size/4, 64), 2);
+ Interface::scatter_matrix(Global_A_stl, Local_A_stl, descA, size, size, blocksize, blocksize);
+ Interface::scatter_matrix(Global_U_stl, Local_U_stl, descU, size, size, blocksize, blocksize);
+ Interface::scatter_matrix(Global_V_stl, Local_V_stl, descV, size, size, blocksize, blocksize);
+ LocalRows = descA[8];
+ LocalCols = Local_A_stl.size()/descA[8];
+
+ // Generic local matrix and vectors initialization
+ Interface::matrix_from_stl(Local_A_ref, Local_A_stl);
+ Interface::matrix_from_stl(Local_A , Local_A_stl);
+ Interface::matrix_from_stl(Local_U_ref, Local_U_stl);
+ Interface::matrix_from_stl(Local_U , Local_U_stl);
+ Interface::matrix_from_stl(Local_V_ref, Local_V_stl);
+ Interface::matrix_from_stl(Local_V , Local_V_stl);
+ Interface::vector_from_stl(Local_s_ref, Local_s_stl);
+ Interface::vector_from_stl(Local_s , Local_s_stl);
+
+ _cost = 2.0*size*size*size;
+ }
+
+
+ // Invalidate copy constructor
+ Action_parallel_svd_decomp(const Action_parallel_svd_decomp&)
+ {
+ INFOS("illegal call to Action_parallel_svd_decomp copy constructor");
+ exit(1);
+ }
+
+ // Destructor
+ ~Action_parallel_svd_decomp()
+ {
+ MESSAGE("Action_parallel_svd_decomp destructor");
+
+ // Deallocation
+ Interface::free_matrix(Local_A_ref, Local_A_stl.size());
+ Interface::free_matrix(Local_A , Local_A_stl.size());
+ Interface::free_matrix(Local_U_ref, Local_U_stl.size());
+ Interface::free_matrix(Local_U , Local_U_stl.size());
+ Interface::free_matrix(Local_V_ref, Local_V_stl.size());
+ Interface::free_matrix(Local_V , Local_V_stl.size());
+ Interface::free_vector(Local_s_ref);
+ Interface::free_vector(Local_s );
+ }
+
+ // Action name
+ static inline std::string name()
+ {
+ return "svd_decomp_" + Interface::name();
+ }
+
+ double nb_op_base()
+ {
+ return _cost;
+ }
+
+ BTL_DONT_INLINE void initialize()
+ {
+ Interface::copy_matrix(Local_A_ref, Local_A, Local_A_stl.size());
+ Interface::copy_matrix(Local_U_ref, Local_U, Local_U_stl.size());
+ Interface::copy_matrix(Local_V_ref, Local_V, Local_V_stl.size());
+ Interface::copy_vector(Local_s_ref, Local_s, Local_s_stl.size());
+ }
+
+ BTL_DONT_INLINE void calculate()
+ {
+ Interface::parallel_svd_decomp(Local_A, descA, Local_U, descU, Local_V, descV, Local_s);
+ }
+
+ BTL_DONT_INLINE void check_result()
+ {
+ }
+
+private:
+ int _size, descA[9], descU[9], descV[9], LocalRows, LocalCols;
+ double _cost;
+ bool iamroot;
+
+ typename Interface::stl_matrix Global_A_stl;
+ typename Interface::stl_matrix Local_A_stl;
+ typename Interface::gene_matrix Local_A_ref;
+ typename Interface::gene_matrix Local_A;
+
+ typename Interface::stl_matrix Global_U_stl;
+ typename Interface::stl_matrix Local_U_stl;
+ typename Interface::gene_matrix Local_U_ref;
+ typename Interface::gene_matrix Local_U;
+
+ typename Interface::stl_matrix Global_V_stl;
+ typename Interface::stl_matrix Local_V_stl;
+ typename Interface::gene_matrix Local_V_ref;
+ typename Interface::gene_matrix Local_V;
+
+ typename Interface::stl_vector Local_s_stl;
+ typename Interface::gene_vector Local_s_ref;
+ typename Interface::gene_vector Local_s;
+};
+
+
+#endif /* ACTION_PARALLEL_SVD_DECOMP_HH */
diff --git a/btl/actions/action_parallel_symm_ev.hh b/btl/actions/action_parallel_symm_ev.hh
index f0af0e3..a4f8237 100644
--- a/btl/actions/action_parallel_symm_ev.hh
+++ b/btl/actions/action_parallel_symm_ev.hh
@@ -59,7 +59,7 @@ public :
// Destructor
~Action_parallel_symm_ev()
{
- MESSAGE("Action_parallel_lu_decomp destructor");
+ MESSAGE("Action_parallel_symm_ev destructor");
// Deallocation
Interface::free_matrix(Local_A_ref, Local_A_stl.size());
@@ -118,4 +118,4 @@ private:
};
-#endif /* ACTION_PARALLEL_LU_DECOMP_HH_ */
+#endif /* ACTION_PARALLEL_SYMM_EV_HH_ */
diff --git a/btl/libs/PBLAS/main.cpp b/btl/libs/PBLAS/main.cpp
index c209afe..f1f7d69 100644
--- a/btl/libs/PBLAS/main.cpp
+++ b/btl/libs/PBLAS/main.cpp
@@ -15,6 +15,7 @@
#include "action_parallel_lu_decomp.hh"
#include "action_parallel_cholesky.hh"
#include "action_parallel_qr_decomp.hh"
+#include "action_parallel_svd_decomp.hh"
#include "action_parallel_symm_ev.hh"
#include <string>
@@ -26,32 +27,46 @@ int main(int argc, char **argv)
bool iamroot = blacsinit(&argc, &argv);
bool
- general_solve=false, qr_decomp=false, lu_decomp=false, cholesky=false,
+ axpy=false, matrix_vector=false,
+ lu_decomp=false, cholesky=false, qr_decomp=false, svd_decomp=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 == "qr_decomp") qr_decomp = true;
- else if (arg == "lu_decomp") lu_decomp = true;
- else if (arg == "cholesky") cholesky = true;
- else if (arg == "symm_ev") symm_ev = true;
+ std::string arg = argv[i];
+ if (arg == "axpy") axpy = true;
+ else if (arg == "matrix_vector") matrix_vector=true;
+ else if (arg == "lu_decomp") lu_decomp = true;
+ else if (arg == "cholesky") cholesky = true;
+ else if (arg == "qr_decomp") qr_decomp = true;
+ else if (arg == "svd_decomp") svd_decomp = true;
+ else if (arg == "symm_ev") symm_ev = true;
+ else if(iamroot) {
+ cerr << "Argument not recognized: " << arg << endl << "Exit" << endl;
+ return 1;
+ }
}
-// if (general_solve)
-// distr_bench<Action_general_solve<lapack_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT, !iamroot);
- if (qr_decomp)
- distr_bench<Action_parallel_qr_decomp<pblas_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT, !iamroot);
+ if (axpy)
+ distr_bench<Action_parallel_axpy<pblas_interface<REAL_TYPE> > >(MIN_AXPY,MAX_AXPY,NB_POINT, !iamroot);
+
+ if (matrix_vector)
+ distr_bench<Action_parallel_matrix_vector_product<pblas_interface<REAL_TYPE> > >(MIN_MV,MAX_MV,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);
+ distr_bench<Action_parallel_cholesky<pblas_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT, !iamroot);
+
+ if (qr_decomp)
+ distr_bench<Action_parallel_qr_decomp<pblas_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT, !iamroot);
+
+ if (svd_decomp)
+ distr_bench<Action_parallel_svd_decomp<pblas_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT, !iamroot);
if (symm_ev)
distr_bench<Action_parallel_symm_ev<pblas_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT, !iamroot);
diff --git a/btl/libs/PBLAS/pblas.h b/btl/libs/PBLAS/pblas.h
index a6cbeb2..e801b8b 100644
--- a/btl/libs/PBLAS/pblas.h
+++ b/btl/libs/PBLAS/pblas.h
@@ -58,6 +58,10 @@ extern "C" {
void psgeqpf_(const int*, const int*, float*, const int*, const int*, const int*, int*, float*, float*, const int*, int*);
void pdgeqpf_(const int*, const int*, double*, const int*, const int*, const int*, int*, double*, double*, const int*, int*);
+ // svd_decomp
+ void psgesvd_(const char*, const char*, const int*, const int*, float*, const int*, const int*, const int*, float*, float*, const int*, const int*, const int*, float*, const int*, const int*, const int*, float*, const int*, int*);
+ void pdgesvd_(const char*, const char*, const int*, const int*, double*, const int*, const int*, const int*, double*, double*, const int*, const int*, const int*, double*, const int*, const int*, const int*, double*, const int*, int*);
+
// symm_ev
void pssyevd_(const char*, const char*, const int*, float*, const int*, const int*, const int*, float*, float*, const int*, const int*, const int*, float*, const int*, int*, const int*, int*);
void pdsyevd_(const char*, const char*, const int*, double*, const int*, const int*, const int*, double*, double*, const int*, const int*, const int*, double*, const int*, int*, const int*, int*);
diff --git a/btl/libs/PBLAS/pblas_interface_impl.hh b/btl/libs/PBLAS/pblas_interface_impl.hh
index 4522946..d71d61e 100644
--- a/btl/libs/PBLAS/pblas_interface_impl.hh
+++ b/btl/libs/PBLAS/pblas_interface_impl.hh
@@ -53,6 +53,8 @@ public:
std::vector<int> ipiv(desc[8] + desc[4]);
PBLAS_FUNC(getrf)(&GlobalRows, &GlobalCols, X, &iONE, &iONE, desc,
&ipiv[0], &info);
+// if (info != 0)
+// cerr << " { LU error : " << info << " } ";
}
static inline void parallel_cholesky(gene_matrix& X, const int* desc)
@@ -61,8 +63,8 @@ public:
const char UPLO = 'U';
int info;
PBLAS_FUNC(potrf)(&UPLO, &N, X, &iONE, &iONE, desc, &info);
- if (info != 0)
- cerr << " { cholesky error : " << info << " } ";
+// if (info != 0)
+// cerr << " { cholesky error : " << info << " } ";
}
static inline void parallel_qr_decomp(gene_matrix& X, const int* desc)
@@ -115,4 +117,26 @@ public:
// if (info != 0)
// cerr << " { symm_ev computation error } ";
}
+
+ static inline void parallel_svd_decomp(gene_matrix& A, int* descA, gene_matrix& U, int *descU, gene_matrix& V, int *descV, gene_vector& s)
+ {
+ const char job = 'V';
+ const int size = descA[2], iONE = 1, iZERO = 0, imONE = -1;
+ std::vector<SCALAR> work;
+ int info, lwork;
+ SCALAR lworkd;
+
+ // Retrieve lwork
+ PBLAS_FUNC(gesvd)(&job, &job, &size, &size, A, &iONE, &iONE, descA, s,
+ U, &iONE, &iONE, descU, V, &iONE, &iONE, descV, &lworkd, &imONE, &info);
+// if (info != 0)
+// cerr << " { svd_decomp lwork error } ";
+ lwork = static_cast<int>(lworkd);
+ work.resize(lwork);
+
+ PBLAS_FUNC(gesvd)(&job, &job, &size, &size, A, &iONE, &iONE, descA, s,
+ U, &iONE, &iONE, descU, V, &iONE, &iONE, descV, &work[0], &lwork, &info);
+// if (info != 0)
+// cerr << " { svd_decomp computation error } ";
+ }
};
diff --git a/pblas.py b/pblas.py
index 792f343..64d1eb7 100644
--- a/pblas.py
+++ b/pblas.py
@@ -6,7 +6,7 @@ class Module(btlbase.BTLBase):
def _initialize(self):
self.libname = "scalapack"
self.avail = ['axpy', 'matrix_vector', 'lu_decomp', 'cholesky',
- 'qr_decomp', 'symm_ev']
+ 'qr_decomp', 'svd_decomp', 'symm_ev']
def _parse_args(self, args):
# Parse arguments
^ permalink raw reply related [flat|nested] 4+ messages in thread
* [gentoo-commits] proj/auto-numerical-bench:unstable commit in: /, btl/actions/, btl/libs/PBLAS/
2011-07-23 22:59 [gentoo-commits] proj/auto-numerical-bench:unstable commit in: /, btl/actions/, btl/libs/PBLAS/ Andrea Arteaga
@ 2011-08-02 18:45 ` Andrea Arteaga
0 siblings, 0 replies; 4+ messages in thread
From: Andrea Arteaga @ 2011-08-02 18:45 UTC (permalink / raw
To: gentoo-commits
commit: 540fa4b69dcc294f9462bf7ab5dc852b79944ee8
Author: spiros <andyspiros <AT> gmail <DOT> com>
AuthorDate: Sat Jul 23 22:59:15 2011 +0000
Commit: Andrea Arteaga <andyspiros <AT> gmail <DOT> com>
CommitDate: Sat Jul 23 22:59:15 2011 +0000
URL: http://git.overlays.gentoo.org/gitweb/?p=proj/auto-numerical-bench.git;a=commit;h=540fa4b6
Added SVD decomposition. Work on eigenvalues action and mat-vec
multiply.
---
.../action_parallel_matrix_vector_product.hh | 52 +++-----
btl/actions/action_parallel_svd_decomp.hh | 134 ++++++++++++++++++++
btl/actions/action_parallel_symm_ev.hh | 4 +-
btl/libs/PBLAS/main.cpp | 39 ++++--
btl/libs/PBLAS/pblas.h | 4 +
btl/libs/PBLAS/pblas_interface_impl.hh | 28 ++++-
pblas.py | 2 +-
7 files changed, 215 insertions(+), 48 deletions(-)
diff --git a/btl/actions/action_parallel_matrix_vector_product.hh b/btl/actions/action_parallel_matrix_vector_product.hh
index 67e64bf..5920115 100644
--- a/btl/actions/action_parallel_matrix_vector_product.hh
+++ b/btl/actions/action_parallel_matrix_vector_product.hh
@@ -22,27 +22,21 @@ public :
MESSAGE("Action_parallel_matrix_vector_product Ctor");
int iZERO = 0, iONE = 1;
- GlobalRows = _size;
- GlobalCols = _size;
- BlockRows = 2;
- BlockCols= 2;
- LocalXCols = 1;
- LocalYCols = 1;
-
int myid, procnum;
blacs_pinfo_(&myid, &procnum);
iamroot = (myid == 0);
// STL matrix and vector initialization
if (iamroot) {
- init_vector<pseudo_random>(Global_A_stl, GlobalRows*GlobalCols);
- init_vector<pseudo_random>(Global_x_stl, GlobalCols);
- init_vector<null_function>(Global_y_stl, GlobalRows);
+ init_vector<pseudo_random>(Global_A_stl, size*size);
+ init_vector<pseudo_random>(Global_x_stl, size);
+ init_vector<null_function>(Global_y_stl, size);
}
- Interface::scatter_matrix(Global_A_stl, Local_A_stl, descA, GlobalRows, GlobalCols, BlockRows, BlockCols);
- Interface::scatter_matrix(Global_x_stl, Local_x_stl, descX, GlobalCols, 1, BlockRows, BlockCols);
- Interface::scatter_matrix(Global_y_stl, Local_y_stl, descY, GlobalRows, 1, BlockRows, BlockCols);
+ const int blocksize = std::max(std::min(size/4, 64), 2);
+ Interface::scatter_matrix(Global_A_stl, Local_A_stl, descA, size, size, blocksize, blocksize);
+ Interface::scatter_matrix(Global_x_stl, Local_x_stl, descX, size, 1, blocksize, blocksize);
+ Interface::scatter_matrix(Global_y_stl, Local_y_stl, descY, size, 1, blocksize, blocksize);
// generic local matrix and vectors initialization
@@ -57,22 +51,22 @@ 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");
+ INFOS("illegal call to Action_parallel_matrix_vector_product copy constructor");
exit(1);
}
// Dtor
- BTL_DONT_INLINE ~Action_parallel_matrix_vector_product( void ){
+ BTL_DONT_INLINE ~Action_parallel_matrix_vector_product(){
- MESSAGE("Action_parallel_matrix_vector_product Dtor");
+ MESSAGE("Action_parallel_matrix_vector_product destructor");
// deallocation
- Interface::free_matrix(Local_A_ref, GlobalRows*GlobalCols);
+ Interface::free_matrix(Local_A_ref, _size*_size);
Interface::free_vector(Local_x_ref);
Interface::free_vector(Local_y_ref);
- Interface::free_matrix(Local_A, GlobalRows*GlobalCols);
+ Interface::free_matrix(Local_A, _size*_size);
Interface::free_vector(Local_x);
Interface::free_vector(Local_y);
@@ -89,35 +83,32 @@ public :
}
BTL_DONT_INLINE void initialize( void ){
- Interface::copy_matrix(Local_A_ref,Local_A,LocalRows*LocalCols);
- Interface::copy_vector(Local_x_ref,Local_x,LocalXRows*LocalXCols);
- Interface::copy_vector(Local_y_ref,Local_y,LocalYRows*LocalYCols);
+ Interface::copy_matrix(Local_A_ref, Local_A, Local_A_stl.size());
+ Interface::copy_vector(Local_x_ref, Local_x, Local_x_stl.size());
+ Interface::copy_vector(Local_y_ref, Local_y, Local_y_stl.size());
}
BTL_DONT_INLINE void calculate( void ) {
- BTL_ASM_COMMENT("#begin matrix_vector_product");
- Interface::parallel_matrix_vector_product(GlobalRows, GlobalCols, Local_A, descA, Local_x, descX, Local_y, descY);
- BTL_ASM_COMMENT("end matrix_vector_product");
+ Interface::parallel_matrix_vector_product(_size, _size, Local_A, descA, Local_x, descX, Local_y, descY);
}
BTL_DONT_INLINE void check_result( void ){
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);
+ Interface::gather_matrix(Global_y_stl, Local_y_stl, descY);
// calculation check
if (iamroot) {
// Compute YTest
- Test_y_stl.resize(GlobalRows);
+ Test_y_stl.resize(_size);
STL_interface<typename Interface::real_type>::matrix_vector_product(Global_A_stl, Global_x_stl, Test_y_stl, _size);
- typename Interface::real_type error =
- STL_interface<typename Interface::real_type>::norm_diff(Global_y_stl, Test_y_stl);
+ typename Interface::real_type error = STL_interface<typename Interface::real_type>::norm_diff(Global_y_stl, Test_y_stl);
if (error > 1e-5)
- std::cerr << "Error: " << error << std::endl;
+ std::cerr << "Error: " << error << " ";
}
}
@@ -142,8 +133,7 @@ private :
typename Interface::gene_vector Local_y;
bool iamroot;
- int _size, GlobalRows, GlobalCols, LocalRows, LocalCols, BlockRows, BlockCols;
- int LocalXRows, LocalXCols, LocalYRows, LocalYCols;
+ int _size;
int descA[9], descX[9], descY[9];
};
diff --git a/btl/actions/action_parallel_svd_decomp.hh b/btl/actions/action_parallel_svd_decomp.hh
new file mode 100644
index 0000000..790ff6d
--- /dev/null
+++ b/btl/actions/action_parallel_svd_decomp.hh
@@ -0,0 +1,134 @@
+#ifndef ACTION_PARALLEL_SVD_DECOMP_HH_
+#define ACTION_PARALLEL_SVD_DECOMP_HH_
+
+#include "utilities.h"
+#include "init/init_function.hh"
+#include "init/init_vector.hh"
+
+#include "lapack_interface.hh"
+#include "STL_interface.hh"
+
+#include <string>
+#include <algorithm>
+
+template<class Interface>
+class Action_parallel_svd_decomp {
+
+public :
+
+ // Constructor
+ BTL_DONT_INLINE Action_parallel_svd_decomp( int size ) : _size(size)
+ {
+ MESSAGE("Action_parallel_svd_decomp Ctor");
+
+ int myid, procnum;
+ blacs_pinfo_(&myid, &procnum);
+ iamroot = (myid == 0);
+
+ // STL matrix and vector initialization
+ if (iamroot) {
+ init_vector<pseudo_random>(Global_A_stl, size*size);
+ init_vector<pseudo_random>(Global_U_stl, size*size);
+ init_vector<pseudo_random>(Global_V_stl, size*size);
+ }
+ init_vector<null_function>(Local_s_stl, size);
+
+ const int blocksize = std::max(std::min(size/4, 64), 2);
+ Interface::scatter_matrix(Global_A_stl, Local_A_stl, descA, size, size, blocksize, blocksize);
+ Interface::scatter_matrix(Global_U_stl, Local_U_stl, descU, size, size, blocksize, blocksize);
+ Interface::scatter_matrix(Global_V_stl, Local_V_stl, descV, size, size, blocksize, blocksize);
+ LocalRows = descA[8];
+ LocalCols = Local_A_stl.size()/descA[8];
+
+ // Generic local matrix and vectors initialization
+ Interface::matrix_from_stl(Local_A_ref, Local_A_stl);
+ Interface::matrix_from_stl(Local_A , Local_A_stl);
+ Interface::matrix_from_stl(Local_U_ref, Local_U_stl);
+ Interface::matrix_from_stl(Local_U , Local_U_stl);
+ Interface::matrix_from_stl(Local_V_ref, Local_V_stl);
+ Interface::matrix_from_stl(Local_V , Local_V_stl);
+ Interface::vector_from_stl(Local_s_ref, Local_s_stl);
+ Interface::vector_from_stl(Local_s , Local_s_stl);
+
+ _cost = 2.0*size*size*size;
+ }
+
+
+ // Invalidate copy constructor
+ Action_parallel_svd_decomp(const Action_parallel_svd_decomp&)
+ {
+ INFOS("illegal call to Action_parallel_svd_decomp copy constructor");
+ exit(1);
+ }
+
+ // Destructor
+ ~Action_parallel_svd_decomp()
+ {
+ MESSAGE("Action_parallel_svd_decomp destructor");
+
+ // Deallocation
+ Interface::free_matrix(Local_A_ref, Local_A_stl.size());
+ Interface::free_matrix(Local_A , Local_A_stl.size());
+ Interface::free_matrix(Local_U_ref, Local_U_stl.size());
+ Interface::free_matrix(Local_U , Local_U_stl.size());
+ Interface::free_matrix(Local_V_ref, Local_V_stl.size());
+ Interface::free_matrix(Local_V , Local_V_stl.size());
+ Interface::free_vector(Local_s_ref);
+ Interface::free_vector(Local_s );
+ }
+
+ // Action name
+ static inline std::string name()
+ {
+ return "svd_decomp_" + Interface::name();
+ }
+
+ double nb_op_base()
+ {
+ return _cost;
+ }
+
+ BTL_DONT_INLINE void initialize()
+ {
+ Interface::copy_matrix(Local_A_ref, Local_A, Local_A_stl.size());
+ Interface::copy_matrix(Local_U_ref, Local_U, Local_U_stl.size());
+ Interface::copy_matrix(Local_V_ref, Local_V, Local_V_stl.size());
+ Interface::copy_vector(Local_s_ref, Local_s, Local_s_stl.size());
+ }
+
+ BTL_DONT_INLINE void calculate()
+ {
+ Interface::parallel_svd_decomp(Local_A, descA, Local_U, descU, Local_V, descV, Local_s);
+ }
+
+ BTL_DONT_INLINE void check_result()
+ {
+ }
+
+private:
+ int _size, descA[9], descU[9], descV[9], LocalRows, LocalCols;
+ double _cost;
+ bool iamroot;
+
+ typename Interface::stl_matrix Global_A_stl;
+ typename Interface::stl_matrix Local_A_stl;
+ typename Interface::gene_matrix Local_A_ref;
+ typename Interface::gene_matrix Local_A;
+
+ typename Interface::stl_matrix Global_U_stl;
+ typename Interface::stl_matrix Local_U_stl;
+ typename Interface::gene_matrix Local_U_ref;
+ typename Interface::gene_matrix Local_U;
+
+ typename Interface::stl_matrix Global_V_stl;
+ typename Interface::stl_matrix Local_V_stl;
+ typename Interface::gene_matrix Local_V_ref;
+ typename Interface::gene_matrix Local_V;
+
+ typename Interface::stl_vector Local_s_stl;
+ typename Interface::gene_vector Local_s_ref;
+ typename Interface::gene_vector Local_s;
+};
+
+
+#endif /* ACTION_PARALLEL_SVD_DECOMP_HH */
diff --git a/btl/actions/action_parallel_symm_ev.hh b/btl/actions/action_parallel_symm_ev.hh
index f0af0e3..a4f8237 100644
--- a/btl/actions/action_parallel_symm_ev.hh
+++ b/btl/actions/action_parallel_symm_ev.hh
@@ -59,7 +59,7 @@ public :
// Destructor
~Action_parallel_symm_ev()
{
- MESSAGE("Action_parallel_lu_decomp destructor");
+ MESSAGE("Action_parallel_symm_ev destructor");
// Deallocation
Interface::free_matrix(Local_A_ref, Local_A_stl.size());
@@ -118,4 +118,4 @@ private:
};
-#endif /* ACTION_PARALLEL_LU_DECOMP_HH_ */
+#endif /* ACTION_PARALLEL_SYMM_EV_HH_ */
diff --git a/btl/libs/PBLAS/main.cpp b/btl/libs/PBLAS/main.cpp
index c209afe..f1f7d69 100644
--- a/btl/libs/PBLAS/main.cpp
+++ b/btl/libs/PBLAS/main.cpp
@@ -15,6 +15,7 @@
#include "action_parallel_lu_decomp.hh"
#include "action_parallel_cholesky.hh"
#include "action_parallel_qr_decomp.hh"
+#include "action_parallel_svd_decomp.hh"
#include "action_parallel_symm_ev.hh"
#include <string>
@@ -26,32 +27,46 @@ int main(int argc, char **argv)
bool iamroot = blacsinit(&argc, &argv);
bool
- general_solve=false, qr_decomp=false, lu_decomp=false, cholesky=false,
+ axpy=false, matrix_vector=false,
+ lu_decomp=false, cholesky=false, qr_decomp=false, svd_decomp=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 == "qr_decomp") qr_decomp = true;
- else if (arg == "lu_decomp") lu_decomp = true;
- else if (arg == "cholesky") cholesky = true;
- else if (arg == "symm_ev") symm_ev = true;
+ std::string arg = argv[i];
+ if (arg == "axpy") axpy = true;
+ else if (arg == "matrix_vector") matrix_vector=true;
+ else if (arg == "lu_decomp") lu_decomp = true;
+ else if (arg == "cholesky") cholesky = true;
+ else if (arg == "qr_decomp") qr_decomp = true;
+ else if (arg == "svd_decomp") svd_decomp = true;
+ else if (arg == "symm_ev") symm_ev = true;
+ else if(iamroot) {
+ cerr << "Argument not recognized: " << arg << endl << "Exit" << endl;
+ return 1;
+ }
}
-// if (general_solve)
-// distr_bench<Action_general_solve<lapack_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT, !iamroot);
- if (qr_decomp)
- distr_bench<Action_parallel_qr_decomp<pblas_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT, !iamroot);
+ if (axpy)
+ distr_bench<Action_parallel_axpy<pblas_interface<REAL_TYPE> > >(MIN_AXPY,MAX_AXPY,NB_POINT, !iamroot);
+
+ if (matrix_vector)
+ distr_bench<Action_parallel_matrix_vector_product<pblas_interface<REAL_TYPE> > >(MIN_MV,MAX_MV,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);
+ distr_bench<Action_parallel_cholesky<pblas_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT, !iamroot);
+
+ if (qr_decomp)
+ distr_bench<Action_parallel_qr_decomp<pblas_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT, !iamroot);
+
+ if (svd_decomp)
+ distr_bench<Action_parallel_svd_decomp<pblas_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT, !iamroot);
if (symm_ev)
distr_bench<Action_parallel_symm_ev<pblas_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT, !iamroot);
diff --git a/btl/libs/PBLAS/pblas.h b/btl/libs/PBLAS/pblas.h
index a6cbeb2..e801b8b 100644
--- a/btl/libs/PBLAS/pblas.h
+++ b/btl/libs/PBLAS/pblas.h
@@ -58,6 +58,10 @@ extern "C" {
void psgeqpf_(const int*, const int*, float*, const int*, const int*, const int*, int*, float*, float*, const int*, int*);
void pdgeqpf_(const int*, const int*, double*, const int*, const int*, const int*, int*, double*, double*, const int*, int*);
+ // svd_decomp
+ void psgesvd_(const char*, const char*, const int*, const int*, float*, const int*, const int*, const int*, float*, float*, const int*, const int*, const int*, float*, const int*, const int*, const int*, float*, const int*, int*);
+ void pdgesvd_(const char*, const char*, const int*, const int*, double*, const int*, const int*, const int*, double*, double*, const int*, const int*, const int*, double*, const int*, const int*, const int*, double*, const int*, int*);
+
// symm_ev
void pssyevd_(const char*, const char*, const int*, float*, const int*, const int*, const int*, float*, float*, const int*, const int*, const int*, float*, const int*, int*, const int*, int*);
void pdsyevd_(const char*, const char*, const int*, double*, const int*, const int*, const int*, double*, double*, const int*, const int*, const int*, double*, const int*, int*, const int*, int*);
diff --git a/btl/libs/PBLAS/pblas_interface_impl.hh b/btl/libs/PBLAS/pblas_interface_impl.hh
index 4522946..d71d61e 100644
--- a/btl/libs/PBLAS/pblas_interface_impl.hh
+++ b/btl/libs/PBLAS/pblas_interface_impl.hh
@@ -53,6 +53,8 @@ public:
std::vector<int> ipiv(desc[8] + desc[4]);
PBLAS_FUNC(getrf)(&GlobalRows, &GlobalCols, X, &iONE, &iONE, desc,
&ipiv[0], &info);
+// if (info != 0)
+// cerr << " { LU error : " << info << " } ";
}
static inline void parallel_cholesky(gene_matrix& X, const int* desc)
@@ -61,8 +63,8 @@ public:
const char UPLO = 'U';
int info;
PBLAS_FUNC(potrf)(&UPLO, &N, X, &iONE, &iONE, desc, &info);
- if (info != 0)
- cerr << " { cholesky error : " << info << " } ";
+// if (info != 0)
+// cerr << " { cholesky error : " << info << " } ";
}
static inline void parallel_qr_decomp(gene_matrix& X, const int* desc)
@@ -115,4 +117,26 @@ public:
// if (info != 0)
// cerr << " { symm_ev computation error } ";
}
+
+ static inline void parallel_svd_decomp(gene_matrix& A, int* descA, gene_matrix& U, int *descU, gene_matrix& V, int *descV, gene_vector& s)
+ {
+ const char job = 'V';
+ const int size = descA[2], iONE = 1, iZERO = 0, imONE = -1;
+ std::vector<SCALAR> work;
+ int info, lwork;
+ SCALAR lworkd;
+
+ // Retrieve lwork
+ PBLAS_FUNC(gesvd)(&job, &job, &size, &size, A, &iONE, &iONE, descA, s,
+ U, &iONE, &iONE, descU, V, &iONE, &iONE, descV, &lworkd, &imONE, &info);
+// if (info != 0)
+// cerr << " { svd_decomp lwork error } ";
+ lwork = static_cast<int>(lworkd);
+ work.resize(lwork);
+
+ PBLAS_FUNC(gesvd)(&job, &job, &size, &size, A, &iONE, &iONE, descA, s,
+ U, &iONE, &iONE, descU, V, &iONE, &iONE, descV, &work[0], &lwork, &info);
+// if (info != 0)
+// cerr << " { svd_decomp computation error } ";
+ }
};
diff --git a/pblas.py b/pblas.py
index 792f343..64d1eb7 100644
--- a/pblas.py
+++ b/pblas.py
@@ -6,7 +6,7 @@ class Module(btlbase.BTLBase):
def _initialize(self):
self.libname = "scalapack"
self.avail = ['axpy', 'matrix_vector', 'lu_decomp', 'cholesky',
- 'qr_decomp', 'symm_ev']
+ 'qr_decomp', 'svd_decomp', 'symm_ev']
def _parse_args(self, args):
# Parse arguments
^ permalink raw reply related [flat|nested] 4+ messages in thread
* [gentoo-commits] proj/auto-numerical-bench:unstable commit in: /, btl/actions/, btl/libs/PBLAS/
@ 2011-08-02 18:45 Andrea Arteaga
2011-07-23 11:46 ` Andrea Arteaga
0 siblings, 1 reply; 4+ messages in thread
From: Andrea Arteaga @ 2011-08-02 18:45 UTC (permalink / raw
To: gentoo-commits
commit: c51172bc46f4b95af6282d8782e4b145911c7afe
Author: spiros <andyspiros <AT> gmail <DOT> com>
AuthorDate: Sat Jul 23 11:45:29 2011 +0000
Commit: Andrea Arteaga <andyspiros <AT> gmail <DOT> com>
CommitDate: Sat Jul 23 11:45:29 2011 +0000
URL: http://git.overlays.gentoo.org/gitweb/?p=proj/auto-numerical-bench.git;a=commit;h=c51172bc
Added working QR decomposition; added working symm_ev (but some negative
MFlops).
---
btl/actions/action_parallel_cholesky.hh | 3 +-
btl/actions/action_parallel_lu_decomp.hh | 3 +-
...el_cholesky.hh => action_parallel_qr_decomp.hh} | 48 +++-----
btl/actions/action_parallel_symm_ev.hh | 121 ++++++++++++++++++++
btl/libs/PBLAS/main.cpp | 14 ++-
btl/libs/PBLAS/pblas.h | 8 ++
btl/libs/PBLAS/pblas_interface_impl.hh | 54 +++++++++-
pblas.py | 3 +-
8 files changed, 213 insertions(+), 41 deletions(-)
diff --git a/btl/actions/action_parallel_cholesky.hh b/btl/actions/action_parallel_cholesky.hh
index f89eb98..05ef3ef 100644
--- a/btl/actions/action_parallel_cholesky.hh
+++ b/btl/actions/action_parallel_cholesky.hh
@@ -39,7 +39,8 @@ public :
Global_A_stl.push_back(temp_stl[r][c]);
}
- Interface::scatter_matrix(Global_A_stl, Local_A_stl, desc, size, size, 64, 64);
+ const int blocksize = std::max(std::min(size/4, 64), 2);
+ Interface::scatter_matrix(Global_A_stl, Local_A_stl, desc, size, size, blocksize, blocksize);
LocalRows = desc[8];
LocalCols = Local_A_stl.size()/desc[8];
diff --git a/btl/actions/action_parallel_lu_decomp.hh b/btl/actions/action_parallel_lu_decomp.hh
index 18b4ac7..d3dc620 100644
--- a/btl/actions/action_parallel_lu_decomp.hh
+++ b/btl/actions/action_parallel_lu_decomp.hh
@@ -29,7 +29,8 @@ public :
init_vector<pseudo_random>(Global_A_stl, size*size);
}
- Interface::scatter_matrix(Global_A_stl, Local_A_stl, desc, size, size, 64, 64);
+ const int blocksize = std::max(std::min(size/4, 64), 2);
+ Interface::scatter_matrix(Global_A_stl, Local_A_stl, desc, size, size, blocksize, blocksize);
LocalRows = desc[8];
LocalCols = Local_A_stl.size()/desc[8];
diff --git a/btl/actions/action_parallel_cholesky.hh b/btl/actions/action_parallel_qr_decomp.hh
similarity index 55%
copy from btl/actions/action_parallel_cholesky.hh
copy to btl/actions/action_parallel_qr_decomp.hh
index f89eb98..a41414c 100644
--- a/btl/actions/action_parallel_cholesky.hh
+++ b/btl/actions/action_parallel_qr_decomp.hh
@@ -1,5 +1,5 @@
-#ifndef ACTION_PARALLEL_CHOLESKY_HH_
-#define ACTION_PARALLEL_CHOLESKY_HH_
+#ifndef ACTION_PARALLEL_QR_DECOMP_HH_
+#define ACTION_PARALLEL_QR_DECOMP_HH_
#include "utilities.h"
#include "init/init_function.hh"
@@ -9,17 +9,17 @@
#include "STL_interface.hh"
#include <string>
+#include <algorithm>
template<class Interface>
-class Action_parallel_cholesky {
- typedef lapack_interface<typename Interface::real_type> LapackInterface;
+class Action_parallel_qr_decomp {
public :
// Constructor
- BTL_DONT_INLINE Action_parallel_cholesky( int size ) : _size(size)
+ BTL_DONT_INLINE Action_parallel_qr_decomp( int size ) : _size(size)
{
- MESSAGE("Action_parallel_cholesky Ctor");
+ MESSAGE("Action_parallel_qr_decomp Ctor");
int myid, procnum;
blacs_pinfo_(&myid, &procnum);
@@ -27,19 +27,11 @@ public :
// STL matrix and vector initialization
if (iamroot) {
- 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]);
+ init_vector<pseudo_random>(Global_A_stl, size*size);
}
- Interface::scatter_matrix(Global_A_stl, Local_A_stl, desc, size, size, 64, 64);
+ const int blocksize = std::max(std::min(size/4, 64), 2);
+ Interface::scatter_matrix(Global_A_stl, Local_A_stl, desc, size, size, blocksize, blocksize);
LocalRows = desc[8];
LocalCols = Local_A_stl.size()/desc[8];
@@ -47,25 +39,21 @@ public :
Interface::matrix_from_stl(Local_A_ref, Local_A_stl);
Interface::matrix_from_stl(Local_A , Local_A_stl);
- _cost = 0;
- for (int j=0; j<_size; ++j) {
- double r = std::max(_size - j -1,0);
- _cost += 2*(r*j+r+j);
- }
+ _cost = 2.0*size*size*size;
}
// Invalidate copy constructor
- Action_parallel_cholesky(const Action_parallel_cholesky&)
+ Action_parallel_qr_decomp(const Action_parallel_qr_decomp&)
{
- INFOS("illegal call to Action_parallel_cholesky copy constructor");
+ INFOS("illegal call to Action_parallel_qr_decomp copy constructor");
exit(1);
}
// Destructor
- ~Action_parallel_cholesky()
+ ~Action_parallel_qr_decomp()
{
- MESSAGE("Action_parallel_cholesky destructor");
+ MESSAGE("Action_parallel_qr_decomp destructor");
// Deallocation
Interface::free_matrix(Local_A_ref, Local_A_stl.size());
@@ -75,7 +63,7 @@ public :
// Action name
static inline std::string name()
{
- return "cholesky_" + Interface::name();
+ return "qr_decomp_" + Interface::name();
}
double nb_op_base()
@@ -90,14 +78,13 @@ public :
BTL_DONT_INLINE void calculate()
{
- Interface::parallel_cholesky(Local_A, desc);
+ Interface::parallel_qr_decomp(Local_A, desc);
}
BTL_DONT_INLINE void check_result()
{
}
-
private:
int _size, desc[9], LocalRows, LocalCols;
double _cost;
@@ -109,4 +96,5 @@ private:
typename Interface::gene_matrix Local_A;
};
-#endif /* ACTION_PARALLEL_CHOLESKY_HH_ */
+
+#endif /* ACTION_PARALLEL_QR_DECOMP_HH_ */
diff --git a/btl/actions/action_parallel_symm_ev.hh b/btl/actions/action_parallel_symm_ev.hh
new file mode 100644
index 0000000..f0af0e3
--- /dev/null
+++ b/btl/actions/action_parallel_symm_ev.hh
@@ -0,0 +1,121 @@
+#ifndef ACTION_PARALLEL_SYMM_EV_HH_
+#define ACTION_PARALLEL_SYMM_EV_HH_
+
+#include "utilities.h"
+#include "init/init_function.hh"
+#include "init/init_vector.hh"
+
+#include "lapack_interface.hh"
+#include "STL_interface.hh"
+
+#include <string>
+
+template<class Interface>
+class Action_parallel_symm_ev {
+
+public :
+
+ // Constructor
+ BTL_DONT_INLINE Action_parallel_symm_ev( int size ) : _size(size)
+ {
+ MESSAGE("Action_parallel_symm_ev constructor");
+
+ int myid, procnum;
+ blacs_pinfo_(&myid, &procnum);
+ iamroot = (myid == 0);
+
+ // STL matrix and vector initialization
+ if (iamroot) {
+ init_vector<pseudo_random>(Global_A_stl, size*size);
+ init_vector<null_function>(Global_Z_stl, size*size);
+ }
+ init_vector<null_function>(Local_w_stl, size);
+
+ const int blocksize = std::max(std::min(size/4, 64), 2);
+ Interface::scatter_matrix(Global_A_stl, Local_A_stl, descA, size, size, blocksize, blocksize);
+ Interface::scatter_matrix(Global_Z_stl, Local_Z_stl, descZ, size, size, blocksize, blocksize);
+ LocalRows = descA[8];
+ LocalCols = Local_A_stl.size()/descA[8];
+
+ // Generic local matrix and vectors initialization
+ Interface::matrix_from_stl(Local_A_ref, Local_A_stl);
+ Interface::matrix_from_stl(Local_A , Local_A_stl);
+ Interface::matrix_from_stl(Local_Z_ref, Local_Z_stl);
+ Interface::matrix_from_stl(Local_Z , Local_Z_stl);
+ Interface::vector_from_stl(Local_w , Local_w_stl);
+ Interface::vector_from_stl(Local_w_ref, Local_w_stl);
+
+ _cost = size*size*size;
+ }
+
+
+ // Invalidate copy constructor
+ Action_parallel_symm_ev(const Action_parallel_symm_ev&)
+ {
+ INFOS("illegal call to Action_parallel_symm_ev copy constructor");
+ exit(1);
+ }
+
+ // Destructor
+ ~Action_parallel_symm_ev()
+ {
+ MESSAGE("Action_parallel_lu_decomp destructor");
+
+ // Deallocation
+ Interface::free_matrix(Local_A_ref, Local_A_stl.size());
+ Interface::free_matrix(Local_A , Local_A_stl.size());
+ Interface::free_matrix(Local_Z_ref, Local_Z_stl.size());
+ Interface::free_matrix(Local_Z , Local_Z_stl.size());
+ Interface::free_vector(Local_w_ref);
+ Interface::free_vector(Local_w );
+ }
+
+ // Action name
+ static inline std::string name()
+ {
+ return "symm_ev_" + Interface::name();
+ }
+
+ double nb_op_base()
+ {
+ return _cost;
+ }
+
+ BTL_DONT_INLINE void initialize()
+ {
+ Interface::copy_matrix(Local_A_ref, Local_A, Local_A_stl.size());
+ Interface::copy_matrix(Local_Z_ref, Local_Z, Local_Z_stl.size());
+ Interface::copy_vector(Local_w_ref, Local_w, Local_w_stl.size());
+ }
+
+ BTL_DONT_INLINE void calculate()
+ {
+ Interface::parallel_symm_ev(Local_A, descA, Local_w, Local_Z, descZ);
+ }
+
+ BTL_DONT_INLINE void check_result()
+ {
+ }
+
+private:
+ int _size, descA[9], descZ[9], LocalRows, LocalCols;
+ double _cost;
+ bool iamroot;
+
+ typename Interface::stl_matrix Global_A_stl;
+ typename Interface::stl_matrix Local_A_stl;
+ typename Interface::gene_matrix Local_A_ref;
+ typename Interface::gene_matrix Local_A;
+
+ typename Interface::stl_matrix Global_Z_stl;
+ typename Interface::stl_matrix Local_Z_stl;
+ typename Interface::gene_matrix Local_Z_ref;
+ typename Interface::gene_matrix Local_Z;
+
+ typename Interface::stl_vector Local_w_stl;
+ typename Interface::gene_vector Local_w_ref;
+ typename Interface::gene_vector Local_w;
+};
+
+
+#endif /* ACTION_PARALLEL_LU_DECOMP_HH_ */
diff --git a/btl/libs/PBLAS/main.cpp b/btl/libs/PBLAS/main.cpp
index e7b636b..c209afe 100644
--- a/btl/libs/PBLAS/main.cpp
+++ b/btl/libs/PBLAS/main.cpp
@@ -14,6 +14,8 @@
#include "action_parallel_matrix_vector_product.hh"
#include "action_parallel_lu_decomp.hh"
#include "action_parallel_cholesky.hh"
+#include "action_parallel_qr_decomp.hh"
+#include "action_parallel_symm_ev.hh"
#include <string>
@@ -24,7 +26,7 @@ int main(int argc, char **argv)
bool iamroot = blacsinit(&argc, &argv);
bool
- general_solve=false, least_squares=false, lu_decomp=false, cholesky=false,
+ general_solve=false, qr_decomp=false, lu_decomp=false, cholesky=false,
symm_ev=false
;
@@ -32,7 +34,7 @@ int main(int argc, char **argv)
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 == "qr_decomp") qr_decomp = true;
else if (arg == "lu_decomp") lu_decomp = true;
else if (arg == "cholesky") cholesky = true;
else if (arg == "symm_ev") symm_ev = true;
@@ -42,8 +44,8 @@ int main(int argc, char **argv)
// 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 (qr_decomp)
+ distr_bench<Action_parallel_qr_decomp<pblas_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);
@@ -51,8 +53,8 @@ int main(int argc, char **argv)
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);
+ if (symm_ev)
+ distr_bench<Action_parallel_symm_ev<pblas_interface<REAL_TYPE> > >(MIN_MM,MAX_MM,NB_POINT, !iamroot);
int iZERO = 0;
diff --git a/btl/libs/PBLAS/pblas.h b/btl/libs/PBLAS/pblas.h
index 973b91c..a6cbeb2 100644
--- a/btl/libs/PBLAS/pblas.h
+++ b/btl/libs/PBLAS/pblas.h
@@ -54,6 +54,14 @@ extern "C" {
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*);
+ // qr_decomp
+ void psgeqpf_(const int*, const int*, float*, const int*, const int*, const int*, int*, float*, float*, const int*, int*);
+ void pdgeqpf_(const int*, const int*, double*, const int*, const int*, const int*, int*, double*, double*, const int*, int*);
+
+ // symm_ev
+ void pssyevd_(const char*, const char*, const int*, float*, const int*, const int*, const int*, float*, float*, const int*, const int*, const int*, float*, const int*, int*, const int*, int*);
+ void pdsyevd_(const char*, const char*, const int*, double*, const int*, const int*, const int*, double*, double*, const int*, const int*, const int*, double*, const int*, 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 1dbf3b9..4522946 100644
--- a/btl/libs/PBLAS/pblas_interface_impl.hh
+++ b/btl/libs/PBLAS/pblas_interface_impl.hh
@@ -61,8 +61,58 @@ public:
const char UPLO = 'U';
int info;
PBLAS_FUNC(potrf)(&UPLO, &N, X, &iONE, &iONE, desc, &info);
+ if (info != 0)
+ cerr << " { cholesky error : " << info << " } ";
+ }
+
+ static inline void parallel_qr_decomp(gene_matrix& X, const int* desc)
+ {
+ const int GlobalRows = desc[2], GlobalCols = desc[3],
+ BlockRows = desc[4], BlockCols = desc[5],
+ ctxt = desc[1];
+
+ int myrow, mycol, nprow, npcol, lwork;
+ SCALAR lworkd;
+ blacs_gridinfo_(&ctxt, &nprow, &npcol, &myrow, &mycol);
+
+ const int iONE = 1, iZERO = 0, imONE = -1,
+ ipivdim = numroc_(&GlobalCols, &BlockCols, &mycol, &iZERO, &npcol);
+ int info;
+ std::vector<int> ipiv(ipivdim);
+ std::vector<SCALAR> tau(ipivdim);
+
+ // Retrieve LWORK
+ PBLAS_FUNC(geqpf)(&GlobalRows, &GlobalCols, X, &iONE, &iONE, desc, &ipiv[0], &tau[0], &lworkd, &imONE, &info);
+ lwork = static_cast<int>(lworkd);
+// if (info != 0)
+// cerr << " { qr_decomp lwork error } ";
+
+ std::vector<SCALAR> work(lwork);
+ PBLAS_FUNC(geqpf)(&GlobalRows, &GlobalCols, X, &iONE, &iONE, desc, &ipiv[0], &tau[0], &work[0], &lwork, &info);
// if (info != 0)
-// cerr << " { cholesky error : " << info << " } ";
+// cerr << " { qr_decomp computation error } ";
}
-};
+ static inline void parallel_symm_ev(gene_matrix& A, const int* descA, gene_vector& w, gene_matrix& Z, const int* descZ)
+ {
+ const char jobz = 'V', uplo = 'u';
+ const int N = descA[2], iONE = 1, iZERO = 0, imONE = -1;
+ std::vector<SCALAR> work;
+ std::vector<int> iwork;
+ int lwork, liwork, info;
+ SCALAR lworkd;
+
+ // Retrieve l(i)work
+ PBLAS_FUNC(syevd)(&jobz, &uplo, &N, A, &iONE, &iONE, descA, w,
+ Z, &iONE, &iONE, descZ, &lworkd, &imONE, &liwork, &imONE, &info);
+ lwork = static_cast<int>(lworkd);
+ work.resize(lwork); iwork.resize(liwork);
+// if (info != 0)
+// cerr << " { symm_ev l(i)work error } ";
+
+ PBLAS_FUNC(syevd)(&jobz, &uplo, &N, A, &iONE, &iONE, descA, w,
+ Z, &iONE, &iONE, descZ, &work[0], &lwork, &iwork[0], &liwork, &info);
+// if (info != 0)
+// cerr << " { symm_ev computation error } ";
+ }
+};
diff --git a/pblas.py b/pblas.py
index 9cd087e..792f343 100644
--- a/pblas.py
+++ b/pblas.py
@@ -5,7 +5,8 @@ numproc = 4
class Module(btlbase.BTLBase):
def _initialize(self):
self.libname = "scalapack"
- self.avail = ['axpy', 'matrix_vector', 'lu_decomp', 'cholesky']
+ self.avail = ['axpy', 'matrix_vector', 'lu_decomp', 'cholesky',
+ 'qr_decomp', 'symm_ev']
def _parse_args(self, args):
# Parse arguments
^ permalink raw reply related [flat|nested] 4+ messages in thread
end of thread, other threads:[~2011-08-02 18:45 UTC | newest]
Thread overview: 4+ messages (download: mbox.gz follow: Atom feed
-- links below jump to the message on this page --
2011-07-23 22:59 [gentoo-commits] proj/auto-numerical-bench:unstable commit in: /, btl/actions/, btl/libs/PBLAS/ Andrea Arteaga
2011-08-02 18:45 ` Andrea Arteaga
-- strict thread matches above, loose matches on Subject: below --
2011-08-02 18:45 Andrea Arteaga
2011-07-23 11:46 ` Andrea Arteaga
This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox