From 2db7e6f8ab20da54cbe4c42b9d8889e003cafc9c Mon Sep 17 00:00:00 2001 From: Yong-Chull Jang Date: Tue, 24 Mar 2020 01:03:24 -0400 Subject: [PATCH 1/5] merge manually Block Lanczos files from Chulwoo's update (last state = commit 731a05 + untracked files) to develop branch; namespace QCD is removed; FIXME: multiple starting vectors result in nan after initial orthogonalization --- .../ImplicitlyRestartedBlockLanczos.h | 1348 +++++++++++++++++ Grid/util/Init.cc | 7 + Grid/util/Init.h | 4 +- tests/lanczos/Test_dwf_block_lanczos.cc | 398 +++++ 4 files changed, 1756 insertions(+), 1 deletion(-) create mode 100644 Grid/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h create mode 100644 tests/lanczos/Test_dwf_block_lanczos.cc diff --git a/Grid/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h b/Grid/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h new file mode 100644 index 00000000..e3afe43c --- /dev/null +++ b/Grid/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h @@ -0,0 +1,1348 @@ + /************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./lib/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h + + Copyright (C) 2015 + +Author: Peter Boyle +Author: Chulwoo Jung +Author: Yong-Chull Jang +Author: Guido Cossu + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License along + with this program; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + + See the full license in the file "LICENSE" in the top level distribution directory + *************************************************************************************/ + /* END LEGAL */ +#ifndef GRID_IRBL_H +#define GRID_IRBL_H + +#include //memset +#ifdef USE_LAPACK +#include +#endif + +#undef USE_LAPACK +#define Glog std::cout << GridLogMessage + +namespace Grid { + +//////////////////////////////////////////////////////////////////////////////// +// Helper class for sorting the evalues AND evectors by Field +// Use pointer swizzle on vectors SHOULD GET RID OF IT SOON! +//////////////////////////////////////////////////////////////////////////////// +template +class SortEigen { + private: + static bool less_lmd(RealD left,RealD right){ + return left > right; + } + static bool less_pair(std::pair& left, + std::pair& right){ + return left.first > (right.first); + } + + public: + void push(std::vector& lmd,std::vector& evec,int N) { + + //////////////////////////////////////////////////////////////////////// + // PAB: FIXME: VERY VERY VERY wasteful: takes a copy of the entire vector set. + // : The vector reorder should be done by pointer swizzle somehow + //////////////////////////////////////////////////////////////////////// + std::vector cpy(lmd.size(),evec[0].Grid()); + for(int i=0;i > emod(lmd.size()); + + for(int i=0;i(lmd[i],&cpy[i]); + + partial_sort(emod.begin(),emod.begin()+N,emod.end(),less_pair); + + typename std::vector >::iterator it = emod.begin(); + for(int i=0;ifirst; + evec[i]=*(it->second); + ++it; + } + } + void push(std::vector& lmd,int N) { + std::partial_sort(lmd.begin(),lmd.begin()+N,lmd.end(),less_lmd); + } + bool saturated(RealD lmd, RealD thrs) { + return fabs(lmd) > fabs(thrs); + } +}; + +enum class LanczosType { irbl, rbl }; + +///////////////////////////////////////////////////////////// +// Implicitly restarted block lanczos +///////////////////////////////////////////////////////////// +template +class ImplicitlyRestartedBlockLanczos { + +private: + + std::string cname = std::string("ImplicitlyRestartedBlockLanczos"); + int MaxIter; // Max iterations + int Nstop; // Number of evecs checked for convergence + int Nu; // Numbeer of vecs in the unit block + int Nk; // Number of converged sought + int Nm; // total number of vectors + int Nblock_k; // Nk/Nu + int Nblock_m; // Nm/Nu + int Nconv_test_interval; // Number of skipped vectors when checking a convergence + RealD eresid; + IRLdiagonalisation diagonalisation; + int split_test; //test split in the first iteration + //////////////////////////////////// + // Embedded objects + //////////////////////////////////// + SortEigen _sort; + LinearOperatorBase &_Linop; + LinearOperatorBase &_SLinop;//for split + OperatorFunction &_poly; + GridRedBlackCartesian * f_grid; + GridRedBlackCartesian * sf_grid; + int mrhs; + + ///////////////////////// + // Constructor + ///////////////////////// +public: + ImplicitlyRestartedBlockLanczos(LinearOperatorBase &Linop, // op + LinearOperatorBase &SLinop, // op + GridRedBlackCartesian * FrbGrid, + GridRedBlackCartesian * SFrbGrid, + int _mrhs, + OperatorFunction & poly, // polynomial + int _Nstop, // really sought vecs + int _Nconv_test_interval, // conv check interval + int _Nu, // vecs in the unit block + int _Nk, // sought vecs + int _Nm, // total vecs + RealD _eresid, // resid in lmd deficit + int _MaxIter, // Max iterations + IRLdiagonalisation _diagonalisation = IRLdiagonaliseWithEigen) + : _Linop(Linop), _SLinop(SLinop), _poly(poly),sf_grid(SFrbGrid),f_grid(FrbGrid), + Nstop(_Nstop), Nconv_test_interval(_Nconv_test_interval), mrhs(_mrhs), + Nu(_Nu), Nk(_Nk), Nm(_Nm), + Nblock_m(_Nm/_Nu), Nblock_k(_Nk/_Nu), + //eresid(_eresid), MaxIter(10), + eresid(_eresid), MaxIter(_MaxIter), + diagonalisation(_diagonalisation),split_test(0) + { assert( (Nk%Nu==0) && (Nm%Nu==0) ); }; + + //////////////////////////////// + // Helpers + //////////////////////////////// + static RealD normalize(Field& v, int if_print=0) + { + RealD nn = norm2(v); + nn = sqrt(nn); +#if 0 + if(if_print && nn < 1e20) + Glog<<"normalize: "<& evec, int k, int if_print=0) + { + typedef typename Field::scalar_type MyComplex; +// MyComplex ip; + ComplexD ip; + + for(int j=0; j 1e-14) + Glog<<"orthogonalize before: "< 1e-14) + Glog<<"orthogonalize after: "<& evec, int k) + { + orthogonalize(w, evec, k,1); + } + + void orthogonalize(std::vector& w, int _Nu, std::vector& evec, int k, int if_print=0) + { + typedef typename Field::scalar_type MyComplex; + MyComplex ip; +// ComplexD ip; + + for(int j=0; j 1e-14) + Glog<<"orthogonalize before: "< 1e-14) + Glog<<"orthogonalize after: "< &inner, std::vector& lhs, int llhs, std::vector& rhs, int lrhs) +{ + typedef typename Field:vector_object vobj; + typedef typename vobj::scalar_type scalar_type; + typedef typename vobj::vector_typeD vector_type; + GridBase *grid = lhs[0]._grid; + assert(grid == rhs[0]._grid; + const int pad = 8; + int total = llhs*lrhs; + assert(inner.size()==total); + int sum_size=grid->SumArraySize(); + +// std::vector inner(total); + Vector sumarray(sum_size*pad*total); + + parallel_for(int thr=0;throSites(),thr,mywork,myoff); + + std::vector< decltype(innerProductD(lhs[0]._odata[0],rhs[0]._odata[0])) > vinner(total,zero); // private to thread; sub summation + for(int ss=myoff;ssGlobalSum(tmp); + inner[i]=tmp; + } +// return inner; +} +#endif + + + void orthogonalize_blockhead(Field& w, std::vector& evec, int k, int Nu) + { + typedef typename Field::scalar_type MyComplex; + MyComplex ip; + + for(int j=0; j& eval, + std::vector& evec, + const std::vector& src, int& Nconv, LanczosType Impl) + { + switch (Impl) { + case LanczosType::irbl: + calc_irbl(eval,evec,src,Nconv); + break; + + case LanczosType::rbl: + calc_rbl(eval,evec,src,Nconv); + break; + } + } + + void calc_irbl(std::vector& eval, + std::vector& evec, + const std::vector& src, int& Nconv) + { + std::string fname = std::string(cname+"::calc_irbl()"); + GridBase *grid = evec[0].Grid(); + assert(grid == src[0].Grid()); + assert( Nu = src.size() ); + + Glog << std::string(74,'*') << std::endl; + Glog << fname + " starting iteration 0 / "<< MaxIter<< std::endl; + Glog << std::string(74,'*') << std::endl; + Glog <<" -- seek Nk = "<< Nk <<" vectors"<< std::endl; + Glog <<" -- accept Nstop = "<< Nstop <<" vectors"<< std::endl; + Glog <<" -- total Nm = "<< Nm <<" vectors"<< std::endl; + Glog <<" -- size of eval = "<< eval.size() << std::endl; + Glog <<" -- size of evec = "<< evec.size() << std::endl; + if ( diagonalisation == IRLdiagonaliseWithEigen ) { + Glog << "Diagonalisation is Eigen "<< std::endl; +#ifdef USE_LAPACK + } else if ( diagonalisation == IRLdiagonaliseWithLAPACK ) { + Glog << "Diagonalisation is LAPACK "<< std::endl; +#endif + } else { + abort(); + } + Glog << std::string(74,'*') << std::endl; + + assert(Nm == evec.size() && Nm == eval.size()); + + std::vector> lmd(Nu,std::vector(Nm,0.0)); + std::vector> lme(Nu,std::vector(Nm,0.0)); + std::vector> lmd2(Nu,std::vector(Nm,0.0)); + std::vector> lme2(Nu,std::vector(Nm,0.0)); + std::vector eval2(Nm); + std::vector resid(Nk); + + Eigen::MatrixXcd Qt = Eigen::MatrixXcd::Zero(Nm,Nm); + Eigen::MatrixXcd Q = Eigen::MatrixXcd::Zero(Nm,Nm); + + std::vector Iconv(Nm); + std::vector B(Nm,grid); // waste of space replicating + + std::vector f(Nu,grid); + std::vector f_copy(Nu,grid); + Field v(grid); + + Nconv = 0; + + RealD beta_k; + + // set initial vector + for (int i=0; i& eval, + std::vector& evec, + const std::vector& src, int& Nconv) + { + std::string fname = std::string(cname+"::calc_rbl()"); + GridBase *grid = evec[0].Grid(); + assert(grid == src[0].Grid()); + assert( Nu = src.size() ); + + int Np = (Nm-Nk); + if (Np > 0 && MaxIter > 1) Np /= MaxIter; + int Nblock_p = Np/Nu; + + Glog << std::string(74,'*') << std::endl; + Glog << fname + " starting iteration 0 / "<< MaxIter<< std::endl; + Glog << std::string(74,'*') << std::endl; + Glog <<" -- seek (min) Nk = "<< Nk <<" vectors"<< std::endl; + Glog <<" -- seek (inc) Np = "<< Np <<" vectors"<< std::endl; + Glog <<" -- seek (max) Nm = "<< Nm <<" vectors"<< std::endl; + Glog <<" -- accept Nstop = "<< Nstop <<" vectors"<< std::endl; + Glog <<" -- size of eval = "<< eval.size() << std::endl; + Glog <<" -- size of evec = "<< evec.size() << std::endl; + if ( diagonalisation == IRLdiagonaliseWithEigen ) { + Glog << "Diagonalisation is Eigen "<< std::endl; +#ifdef USE_LAPACK + } else if ( diagonalisation == IRLdiagonaliseWithLAPACK ) { + Glog << "Diagonalisation is LAPACK "<< std::endl; +#endif + } else { + abort(); + } + Glog << std::string(74,'*') << std::endl; + + assert(Nm == evec.size() && Nm == eval.size()); + + std::vector> lmd(Nu,std::vector(Nm,0.0)); + std::vector> lme(Nu,std::vector(Nm,0.0)); + std::vector> lmd2(Nu,std::vector(Nm,0.0)); + std::vector> lme2(Nu,std::vector(Nm,0.0)); + std::vector eval2(Nk); + std::vector resid(Nm); + + Eigen::MatrixXcd Qt = Eigen::MatrixXcd::Zero(Nm,Nm); + Eigen::MatrixXcd Q = Eigen::MatrixXcd::Zero(Nm,Nm); + + std::vector Iconv(Nm); + std::vector B(Nm,grid); // waste of space replicating + + std::vector f(Nu,grid); + std::vector f_copy(Nu,grid); + Field v(grid); + + Nconv = 0; + + RealD beta_k; + + // set initial vector + for (int i=0; i>& lmd, + std::vector>& lme, + std::vector& evec, + std::vector& w, + std::vector& w_copy, + int b) + { + const RealD tiny = 1.0e-20; + + int Nu = w.size(); + int Nm = evec.size(); + assert( b < Nm/Nu ); +// GridCartesian *grid = evec[0]._grid; + + // converts block index to full indicies for an interval [L,R) + int L = Nu*b; + int R = Nu*(b+1); + + Real beta; + + Glog << "Using split grid"<< std::endl; +// LatticeGaugeField s_Umu(SGrid); + assert((Nu%mrhs)==0); + std::vector in(mrhs,f_grid); + + Field s_in(sf_grid); + Field s_out(sf_grid); + // unnecessary copy. Can or should it be avoided? +int k_start = 0; +while ( k_start < Nu) { + Glog << "k_start= "<0) { + for (int u=0; u0) { + for (int u=0; u0) { + // orthogonalize_blockhead(w[0],evec,b,Nu); + // for (int u=1; u& eval, + std::vector>& lmd, + std::vector>& lme, + int Nu, int Nk, int Nm, + Eigen::MatrixXcd & Qt, // Nm x Nm + GridBase *grid) + { + assert( Nk%Nu == 0 && Nm%Nu == 0 ); + assert( Nk <= Nm ); + Eigen::MatrixXcd BlockTriDiag = Eigen::MatrixXcd::Zero(Nk,Nk); + + for ( int u=0; u eigensolver(BlockTriDiag); + + for (int i = 0; i < Nk; i++) { + eval[Nk-1-i] = eigensolver.eigenvalues()(i); + } + for (int i = 0; i < Nk; i++) { + for (int j = 0; j < Nk; j++) { + Qt(j,Nk-1-i) = eigensolver.eigenvectors()(j,i); + //Qt(Nk-1-i,j) = eigensolver.eigenvectors()(i,j); + //Qt(i,j) = eigensolver.eigenvectors()(i,j); + } + } + } + +#ifdef USE_LAPACK + void diagonalize_lapack(std::vector& eval, + std::vector>& lmd, + std::vector>& lme, + int Nu, int Nk, int Nm, + Eigen::MatrixXcd & Qt, // Nm x Nm + GridBase *grid) + { + Glog << "diagonalize_lapack: Nu= "<_Nprocessors; + int node = grid->_processor; + int interval = (NN/total)+1; + double vl = 0.0, vu = 0.0; + MKL_INT il = interval*node+1 , iu = interval*(node+1); + if (iu > NN) iu=NN; + Glog << "node "<= il-1; i--){ + evals_tmp[i] = evals_tmp[i - (il-1)]; + if (il>1) evals_tmp[i-(il-1)]=0.; + for (int j = 0; j< NN; j++){ + evec_tmp[i*NN+j] = evec_tmp[(i - (il-1))*NN+j]; + if (il>1) { + evec_tmp[(i-(il-1))*NN+j].imag=0.; + evec_tmp[(i-(il-1))*NN+j].real=0.; + } + } + } + } + { + grid->GlobalSumVector(evals_tmp,NN); + grid->GlobalSumVector((double*)evec_tmp,2*NN*NN); + } + } + // Safer to sort instead of just reversing it, + // but the document of the routine says evals are sorted in increasing order. + // qr gives evals in decreasing order. +// for(int i=0;i + ( evec_tmp[i*Nk+j].real, + evec_tmp[i*Nk+j].imag); +// ( evec_tmp[(Nk-1-j)*Nk+Nk-1-i].real, +// evec_tmp[(Nk-1-j)*Nk+Nk-1-i].imag); + + } + } + +if (1){ + Eigen::SelfAdjointEigenSolver eigensolver(BlockTriDiag); + + for (int i = 0; i < Nk; i++) { + Glog << "eval = "<& eval, + std::vector>& lmd, + std::vector>& lme, + int Nu, int Nk, int Nm, + Eigen::MatrixXcd & Qt, + GridBase *grid) + { + Qt = Eigen::MatrixXcd::Identity(Nm,Nm); + if ( diagonalisation == IRLdiagonaliseWithEigen ) { + diagonalize_Eigen(eval,lmd,lme,Nu,Nk,Nm,Qt,grid); +#ifdef USE_LAPACK + } else if ( diagonalisation == IRLdiagonaliseWithLAPACK ) { + diagonalize_lapack(eval,lmd,lme,Nu,Nk,Nm,Qt,grid); +#endif + } else { + assert(0); + } + } + + + void unpackHermitBlockTriDiagMatToEigen( + std::vector>& lmd, + std::vector>& lme, + int Nu, int Nb, int Nk, int Nm, + Eigen::MatrixXcd& M) + { + //Glog << "unpackHermitBlockTriDiagMatToEigen() begin" << '\n'; + assert( Nk%Nu == 0 && Nm%Nu == 0 ); + assert( Nk <= Nm ); + M = Eigen::MatrixXcd::Zero(Nk,Nk); + + // rearrange + for ( int u=0; u>& lmd, + std::vector>& lme, + int Nu, int Nb, int Nk, int Nm, + Eigen::MatrixXcd& M) + { + //Glog << "packHermitBlockTriDiagMatfromEigen() begin" << '\n'; + assert( Nk%Nu == 0 && Nm%Nu == 0 ); + assert( Nk <= Nm ); + + // rearrange + for ( int u=0; u QRD(Mtmp); + Q = QRD.householderQ(); + R = QRD.matrixQR(); // upper triangular part is the R matrix. + // lower triangular part used to represent series + // of Q sequence. + + // equivalent operation of Qprod *= Q + //M = Eigen::MatrixXcd::Zero(Nm,Nm); + + //for (int i=0; i Nm) kmax = Nm; + for (int k=i; ki) M(i,j) = conj(M(j,i)); + // if (i-j > Nu || j-i > Nu) M(i,j) = 0.; + // } + //} + + //Glog << "shiftedQRDecompEigen() end" << endl; + } + + void exampleQRDecompEigen(void) + { + Eigen::MatrixXd A = Eigen::MatrixXd::Zero(3,3); + Eigen::MatrixXd Q = Eigen::MatrixXd::Zero(3,3); + Eigen::MatrixXd R = Eigen::MatrixXd::Zero(3,3); + Eigen::MatrixXd P = Eigen::MatrixXd::Zero(3,3); + + A(0,0) = 12.0; + A(0,1) = -51.0; + A(0,2) = 4.0; + A(1,0) = 6.0; + A(1,1) = 167.0; + A(1,2) = -68.0; + A(2,0) = -4.0; + A(2,1) = 24.0; + A(2,2) = -41.0; + + Glog << "matrix A before ColPivHouseholder" << std::endl; + for ( int i=0; i<3; i++ ) { + for ( int j=0; j<3; j++ ) { + Glog << "A[" << i << "," << j << "] = " << A(i,j) << '\n'; + } + } + Glog << std::endl; + + Eigen::ColPivHouseholderQR QRD(A); + + Glog << "matrix A after ColPivHouseholder" << std::endl; + for ( int i=0; i<3; i++ ) { + for ( int j=0; j<3; j++ ) { + Glog << "A[" << i << "," << j << "] = " << A(i,j) << '\n'; + } + } + Glog << std::endl; + + Glog << "HouseholderQ with sequence lenth = nonzeroPiviots" << std::endl; + Q = QRD.householderQ().setLength(QRD.nonzeroPivots()); + for ( int i=0; i<3; i++ ) { + for ( int j=0; j<3; j++ ) { + Glog << "Q[" << i << "," << j << "] = " << Q(i,j) << '\n'; + } + } + Glog << std::endl; + + Glog << "HouseholderQ with sequence lenth = 1" << std::endl; + Q = QRD.householderQ().setLength(1); + for ( int i=0; i<3; i++ ) { + for ( int j=0; j<3; j++ ) { + Glog << "Q[" << i << "," << j << "] = " << Q(i,j) << '\n'; + } + } + Glog << std::endl; + + Glog << "HouseholderQ with sequence lenth = 2" << std::endl; + Q = QRD.householderQ().setLength(2); + for ( int i=0; i<3; i++ ) { + for ( int j=0; j<3; j++ ) { + Glog << "Q[" << i << "," << j << "] = " << Q(i,j) << '\n'; + } + } + Glog << std::endl; + + Glog << "matrixR" << std::endl; + R = QRD.matrixR(); + for ( int i=0; i<3; i++ ) { + for ( int j=0; j<3; j++ ) { + Glog << "R[" << i << "," << j << "] = " << R(i,j) << '\n'; + } + } + Glog << std::endl; + + Glog << "rank = " << QRD.rank() << std::endl; + Glog << "threshold = " << QRD.threshold() << std::endl; + + Glog << "matrixP" << std::endl; + P = QRD.colsPermutation(); + for ( int i=0; i<3; i++ ) { + for ( int j=0; j<3; j++ ) { + Glog << "P[" << i << "," << j << "] = " << P(i,j) << '\n'; + } + } + Glog << std::endl; + + + Glog << "QR decomposition without column pivoting" << std::endl; + + A(0,0) = 12.0; + A(0,1) = -51.0; + A(0,2) = 4.0; + A(1,0) = 6.0; + A(1,1) = 167.0; + A(1,2) = -68.0; + A(2,0) = -4.0; + A(2,1) = 24.0; + A(2,2) = -41.0; + + Glog << "matrix A before Householder" << std::endl; + for ( int i=0; i<3; i++ ) { + for ( int j=0; j<3; j++ ) { + Glog << "A[" << i << "," << j << "] = " << A(i,j) << '\n'; + } + } + Glog << std::endl; + + Eigen::HouseholderQR QRDplain(A); + + Glog << "HouseholderQ" << std::endl; + Q = QRDplain.householderQ(); + for ( int i=0; i<3; i++ ) { + for ( int j=0; j<3; j++ ) { + Glog << "Q[" << i << "," << j << "] = " << Q(i,j) << '\n'; + } + } + Glog << std::endl; + + Glog << "matrix A after Householder" << std::endl; + for ( int i=0; i<3; i++ ) { + for ( int j=0; j<3; j++ ) { + Glog << "A[" << i << "," << j << "] = " << A(i,j) << '\n'; + } + } + Glog << std::endl; + } + + }; +} +#undef Glog +#undef USE_LAPACK +#endif diff --git a/Grid/util/Init.cc b/Grid/util/Init.cc index 472013f4..4611b209 100644 --- a/Grid/util/Init.cc +++ b/Grid/util/Init.cc @@ -162,6 +162,13 @@ void GridCmdOptionInt(std::string &str,int & val) return; } +// ypj [add] +void GridCmdOptionFloat(std::string &str,double & val) +{ + std::stringstream ss(str); + ss>>val; + return; +} void GridParseLayout(char **argv,int argc, Coordinate &latt_c, diff --git a/Grid/util/Init.h b/Grid/util/Init.h index f7f032ba..cb8b6c49 100644 --- a/Grid/util/Init.h +++ b/Grid/util/Init.h @@ -56,7 +56,9 @@ std::string GridCmdVectorIntToString(const VectorInt & vec); void GridCmdOptionCSL(std::string str,std::vector & vec); template void GridCmdOptionIntVector(std::string &str,VectorInt & vec); - +// ypj [add] +void GridCmdOptionInt(std::string &str,int & val); +void GridCmdOptionFloat(std::string &str,double & val); void GridParseLayout(char **argv,int argc, std::vector &latt, diff --git a/tests/lanczos/Test_dwf_block_lanczos.cc b/tests/lanczos/Test_dwf_block_lanczos.cc new file mode 100644 index 00000000..2aa83be1 --- /dev/null +++ b/tests/lanczos/Test_dwf_block_lanczos.cc @@ -0,0 +1,398 @@ + /************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./tests/Test_dwf_block_lanczos.cc + + Copyright (C) 2015 + +Author: Peter Boyle + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License along + with this program; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + + See the full license in the file "LICENSE" in the top level distribution directory + *************************************************************************************/ + /* END LEGAL */ +#include +#include +#include + +using namespace std; +using namespace Grid; +//using namespace Grid::QCD; + +//typedef typename GparityDomainWallFermionR::FermionField FermionField; +typedef typename ZMobiusFermionR::FermionField FermionField; + +RealD AllZero(RealD x){ return 0.;} + +class CmdJobParams +{ + public: + std::string gaugefile; + + int Ls; + double mass; + double M5; + double mob_b; + std::vector omega; + std::vector boundary_phase; + std::vector mpi_split; + + LanczosType Impl; + int Nu; + int Nk; + int Np; + int Nm; + int Nstop; + int Ntest; + int MaxIter; + double resid; + + double low; + double high; + int order; + + CmdJobParams() + : gaugefile("Hot"), + Ls(8), mass(0.01), M5(1.8), mob_b(1.5), + Impl(LanczosType::irbl),mpi_split(4,1), + Nu(4), Nk(200), Np(200), Nstop(100), Ntest(1), MaxIter(10), resid(1.0e-8), + low(0.2), high(5.5), order(11) + {Nm=Nk+Np;}; + + void Parse(char **argv, int argc); +}; + + +void CmdJobParams::Parse(char **argv,int argc) +{ + std::string arg; + std::vector vi; + double re,im; + int expect, idx; + std::string vstr; + std::ifstream pfile; + + if( GridCmdOptionExists(argv,argv+argc,"--gconf") ){ + gaugefile = GridCmdOptionPayload(argv,argv+argc,"--gconf"); + } + + if( GridCmdOptionExists(argv,argv+argc,"--phase") ){ + arg = GridCmdOptionPayload(argv,argv+argc,"--phase"); + pfile.open(arg); + assert(pfile); + expect = 0; + while( pfile >> vstr ) { + if ( vstr.compare("boundary_phase") == 0 ) { + pfile >> vstr; + GridCmdOptionInt(vstr,idx); + assert(expect==idx); + pfile >> vstr; + GridCmdOptionFloat(vstr,re); + pfile >> vstr; + GridCmdOptionFloat(vstr,im); + boundary_phase.push_back({re,im}); + expect++; + } + } + pfile.close(); + } else { + for (int i=0; i<4; ++i) boundary_phase.push_back({1.,0.}); + } + + if( GridCmdOptionExists(argv,argv+argc,"--omega") ){ + arg = GridCmdOptionPayload(argv,argv+argc,"--omega"); + pfile.open(arg); + assert(pfile); + Ls = 0; + while( pfile >> vstr ) { + if ( vstr.compare("omega") == 0 ) { + pfile >> vstr; + GridCmdOptionInt(vstr,idx); + assert(Ls==idx); + pfile >> vstr; + GridCmdOptionFloat(vstr,re); + pfile >> vstr; + GridCmdOptionFloat(vstr,im); + omega.push_back({re,im}); + Ls++; + } + } + pfile.close(); + } else { + if( GridCmdOptionExists(argv,argv+argc,"--Ls") ){ + arg = GridCmdOptionPayload(argv,argv+argc,"--Ls"); + GridCmdOptionInt(arg,Ls); + } + } + + if( GridCmdOptionExists(argv,argv+argc,"--mass") ){ + arg = GridCmdOptionPayload(argv,argv+argc,"--mass"); + GridCmdOptionFloat(arg,mass); + } + + if( GridCmdOptionExists(argv,argv+argc,"--M5") ){ + arg = GridCmdOptionPayload(argv,argv+argc,"--M5"); + GridCmdOptionFloat(arg,M5); + } + + if( GridCmdOptionExists(argv,argv+argc,"--mob_b") ){ + arg = GridCmdOptionPayload(argv,argv+argc,"--mob_b"); + GridCmdOptionFloat(arg,mob_b); + } + + if( GridCmdOptionExists(argv,argv+argc,"--irbl") ){ + arg = GridCmdOptionPayload(argv,argv+argc,"--irbl"); + GridCmdOptionIntVector(arg,vi); + Nu = vi[0]; + Nk = vi[1]; + Np = vi[2]; + Nstop = vi[3]; + MaxIter = vi[4]; + // ypj[fixme] mode overriding message is needed. + Impl = LanczosType::irbl; + Nm = Nk+Np; + } + + // block Lanczos with explicit extension of its dimensions + if( GridCmdOptionExists(argv,argv+argc,"--rbl") ){ + arg = GridCmdOptionPayload(argv,argv+argc,"--rbl"); + GridCmdOptionIntVector(arg,vi); + Nu = vi[0]; + Nk = vi[1]; + Np = vi[2]; // vector space is enlarged by adding Np vectors + Nstop = vi[3]; + MaxIter = vi[4]; + // ypj[fixme] mode overriding message is needed. + Impl = LanczosType::rbl; + Nm = Nk+Np*MaxIter; + } + +#if 1 + // block Lanczos with explicit extension of its dimensions + if( GridCmdOptionExists(argv,argv+argc,"--split") ){ + arg = GridCmdOptionPayload(argv,argv+argc,"--split"); + GridCmdOptionIntVector(arg,vi); + for(int i=0;i seeds4({1,2,3,4}); + std::vector seeds5({5,6,7,8}); + GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5); + GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4); + // ypj [note] why seed RNG5 again? bug? In this case, run with a default seed(). + GridParallelRNG RNG5rb(FrbGrid); RNG5rb.SeedFixedIntegers(seeds5); + + LatticeGaugeField Umu(UGrid); + std::vector U(4,UGrid); + + if ( JP.gaugefile.compare("Hot") == 0 ) { + SU3::HotConfiguration(RNG4, Umu); + } else { + FieldMetaData header; + NerscIO::readConfiguration(Umu,header,JP.gaugefile); + // ypj [fixme] additional checks for the loaded configuration? + } + + for(int mu=0;mu(Umu,mu); + } + + RealD mass = JP.mass; + RealD M5 = JP.M5; + +// ypj [fixme] flexible support for a various Fermions +// RealD mob_b = JP.mob_b; // Gparity +// std::vector omega; // ZMobius + +// GparityMobiusFermionD ::ImplParams params; +// std::vector twists({1,1,1,0}); +// params.twists = twists; +// GparityMobiusFermionR Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,mob_b,mob_b-1.,params); +// SchurDiagTwoOperator HermOp(Ddwf); + + +// int mrhs = JP.Nu; + int Ndir=4; + auto mpi_layout = GridDefaultMpi(); + std::vector mpi_split (Ndir,1); +#if 0 + int tmp=mrhs, dir=0; + std::cout << GridLogMessage << "dir= "< evec(JP.Nm,FrbGrid); + for(int i=0;i<1;++i){ + std::cout << GridLogMessage << i <<" / "<< JP.Nm <<" grid pointer "<< evec[i].Grid() << std::endl; + }; + + int Nconv; + IRBL.calc(eval,evec,src,Nconv,JP.Impl); + + + Grid_finalize(); +} From 9266b89ad8c8cb62e9f70c527b8f0f28a7cd94e5 Mon Sep 17 00:00:00 2001 From: Yong-Chull Jang Date: Wed, 25 Mar 2020 15:45:50 -0400 Subject: [PATCH 2/5] fix rngs issue; block Lanczos is working --- Grid/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h | 2 +- tests/lanczos/Test_dwf_block_lanczos.cc | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/Grid/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h b/Grid/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h index e3afe43c..5076a527 100644 --- a/Grid/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h +++ b/Grid/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h @@ -180,7 +180,7 @@ public: Glog<<"orthogonalize after: "<& evec, int k) { diff --git a/tests/lanczos/Test_dwf_block_lanczos.cc b/tests/lanczos/Test_dwf_block_lanczos.cc index 2aa83be1..a55c06a4 100644 --- a/tests/lanczos/Test_dwf_block_lanczos.cc +++ b/tests/lanczos/Test_dwf_block_lanczos.cc @@ -363,7 +363,7 @@ int main (int argc, char ** argv) std::vector eval(JP.Nm); std::vector src(JP.Nu,FrbGrid); -if (1) +if (0) { std::cout << GridLogMessage << "Using RNG5"< Date: Mon, 30 Mar 2020 18:36:21 -0400 Subject: [PATCH 3/5] first working version of Gram Schmidt using cublas gemm; explicit data type and site vector size has to be removed --- .gitignore | 2 + .../ImplicitlyRestartedBlockLanczos.h | 150 ++++++++++++++++-- tests/lanczos/Test_dwf_block_lanczos.cc | 4 +- 3 files changed, 144 insertions(+), 12 deletions(-) diff --git a/.gitignore b/.gitignore index 5338acb9..13efd67c 100644 --- a/.gitignore +++ b/.gitignore @@ -10,6 +10,8 @@ *~ *# *.sublime-* +.ctags +tags # Precompiled Headers # ####################### diff --git a/Grid/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h b/Grid/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h index 5076a527..7cc11653 100644 --- a/Grid/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h +++ b/Grid/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h @@ -39,6 +39,10 @@ Author: Guido Cossu #undef USE_LAPACK #define Glog std::cout << GridLogMessage +#ifdef GRID_NVCC +#include "cublas_v2.h" +#endif + namespace Grid { //////////////////////////////////////////////////////////////////////////////// @@ -89,6 +93,12 @@ class SortEigen { enum class LanczosType { irbl, rbl }; +enum IRBLdiagonalisation { + IRBLdiagonaliseWithDSTEGR, + IRBLdiagonaliseWithQR, + IRBLdiagonaliseWithEigen +}; + ///////////////////////////////////////////////////////////// // Implicitly restarted block lanczos ///////////////////////////////////////////////////////////// @@ -107,7 +117,7 @@ private: int Nblock_m; // Nm/Nu int Nconv_test_interval; // Number of skipped vectors when checking a convergence RealD eresid; - IRLdiagonalisation diagonalisation; + IRBLdiagonalisation diagonalisation; int split_test; //test split in the first iteration //////////////////////////////////// // Embedded objects @@ -137,7 +147,7 @@ public: int _Nm, // total vecs RealD _eresid, // resid in lmd deficit int _MaxIter, // Max iterations - IRLdiagonalisation _diagonalisation = IRLdiagonaliseWithEigen) + IRBLdiagonalisation _diagonalisation = IRBLdiagonaliseWithEigen) : _Linop(Linop), _SLinop(SLinop), _poly(poly),sf_grid(SFrbGrid),f_grid(FrbGrid), Nstop(_Nstop), Nconv_test_interval(_Nconv_test_interval), mrhs(_mrhs), Nu(_Nu), Nk(_Nk), Nm(_Nm), @@ -211,7 +221,126 @@ public: #endif }} for(int i=0; i<_Nu; ++i) - normalize(w[i],if_print); + assert(normalize(w[i],if_print) !=0); + } + + + void orthogonalize_blas(std::vector& w, int _Nu, std::vector& evec, int _R, int _print=0) + { +#ifdef GRID_NVCC + Glog << "cuBLAS orthogonalize" << std::endl; + + typedef typename Field::vector_object vobj; + typedef typename vobj::scalar_type scalar_type; + typedef typename vobj::vector_type vector_type; + + typedef typename Field::scalar_type MyComplex; + + GridBase *grid = w[0].Grid(); + //grid->show_decomposition(); + //const uint64_t nsimd = grid->Nsimd(); + const uint64_t sites = grid->lSites(); + + //auto w_v = w[0].View(); + //cuDoubleComplex *z = reinterpret_cast(&w_v._odata[0]); + //cuDoubleComplex *z = w_v._odata._internal; + //thread_for(ss,w_v.size(),{ + // Glog << w_v[ss] << std::endl; + //}); + //w_v[0] + //exit(0); + //scalar_type *z = (scalar_type *)&w_v[0]; // OK + //cuDoubleComplex *z = reinterpret_cast(&w_v[0]); // OK + + cudaError_t cudaStat; + + cuDoubleComplex *w_acc, *evec_acc, *c_acc; + + cudaStat = cudaMallocManaged((void **)&w_acc, _Nu*sites*12*sizeof(cuDoubleComplex)); + Glog << cudaStat << std::endl; + cudaStat = cudaMallocManaged((void **)&evec_acc, _R*sites*12*sizeof(cuDoubleComplex)); + Glog << cudaStat << std::endl; + cudaStat = cudaMallocManaged((void **)&c_acc, _Nu*_R*12*sizeof(cuDoubleComplex)); + Glog << cudaStat << std::endl; + + Glog << "cuBLAS prepare array"<< std::endl; +#if 0 // a trivial test + for (int col=0; col<_Nu; ++col) { + for (size_t row=0; row(&w_v[0]); + for (size_t row=0; row(&evec_v[0]); + for (size_t row=0; rowGlobalSumVector((double*)c_acc,2*_Nu*_R); + + cublasDestroy(handle); + + Glog << "cuBLAS Zgemm done"<< std::endl; + + for (int i=0; i<_Nu; ++i) { + for (size_t j=0; j<_R; ++j) { + cuDoubleComplex z = c_acc[i*_R+j]; + MyComplex ip(z.x,z.y); + if (_print) { + Glog << "[" << j << "," << i << "] = " + << z.x << " + i " << z.y << std::endl; + } + w[i] = w[i] - ip * evec[j]; + } + assert(normalize(w[i],_print)!=0); + } + + cudaFree(w_acc); + cudaFree(evec_acc); + cudaFree(c_acc); + + Glog << "cuBLAS orthogonalize done" << std::endl; +#else + Glog << "BLAS wrapper is not implemented" << std::endl; + exit(1); +#endif } @@ -310,10 +439,10 @@ for( int i =0;i eval(JP.Nm); From ac7090e6d3af4e9d468323bd2a03e2f2437b813f Mon Sep 17 00:00:00 2001 From: Yong-Chull Jang Date: Mon, 30 Mar 2020 22:25:50 -0400 Subject: [PATCH 4/5] block Lanczos cublas buffer is set at the inital step; buffer width is fixed to the block size then cublas Zgemm is called multiple times --- .../ImplicitlyRestartedBlockLanczos.h | 173 ++++++++++-------- 1 file changed, 92 insertions(+), 81 deletions(-) diff --git a/Grid/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h b/Grid/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h index 7cc11653..cd6bd682 100644 --- a/Grid/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h +++ b/Grid/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h @@ -110,7 +110,7 @@ private: std::string cname = std::string("ImplicitlyRestartedBlockLanczos"); int MaxIter; // Max iterations int Nstop; // Number of evecs checked for convergence - int Nu; // Numbeer of vecs in the unit block + int Nu; // Number of vecs in the unit block int Nk; // Number of converged sought int Nm; // total number of vectors int Nblock_k; // Nk/Nu @@ -129,7 +129,15 @@ private: GridRedBlackCartesian * f_grid; GridRedBlackCartesian * sf_grid; int mrhs; - + ///////////////////////// + // BLAS objects + ///////////////////////// +#ifdef GRID_NVCC + cudaError_t cudaStat; + cuDoubleComplex *w_acc, *evec_acc, *c_acc; +#endif + int Nevec_acc; // Number of eigenvectors stored in the buffer evec_acc + ///////////////////////// // Constructor ///////////////////////// @@ -154,7 +162,8 @@ public: Nblock_m(_Nm/_Nu), Nblock_k(_Nk/_Nu), //eresid(_eresid), MaxIter(10), eresid(_eresid), MaxIter(_MaxIter), - diagonalisation(_diagonalisation),split_test(0) + diagonalisation(_diagonalisation),split_test(0), + Nevec_acc(_Nu) { assert( (Nk%Nu==0) && (Nm%Nu==0) ); }; //////////////////////////////// @@ -225,7 +234,7 @@ public: } - void orthogonalize_blas(std::vector& w, int _Nu, std::vector& evec, int _R, int _print=0) + void orthogonalize_blas(std::vector& w, std::vector& evec, int R, int do_print=0) { #ifdef GRID_NVCC Glog << "cuBLAS orthogonalize" << std::endl; @@ -237,105 +246,86 @@ public: typedef typename Field::scalar_type MyComplex; GridBase *grid = w[0].Grid(); - //grid->show_decomposition(); - //const uint64_t nsimd = grid->Nsimd(); const uint64_t sites = grid->lSites(); - - //auto w_v = w[0].View(); - //cuDoubleComplex *z = reinterpret_cast(&w_v._odata[0]); - //cuDoubleComplex *z = w_v._odata._internal; - //thread_for(ss,w_v.size(),{ - // Glog << w_v[ss] << std::endl; - //}); - //w_v[0] - //exit(0); - //scalar_type *z = (scalar_type *)&w_v[0]; // OK - //cuDoubleComplex *z = reinterpret_cast(&w_v[0]); // OK - - cudaError_t cudaStat; + + int Nbatch = R/Nevec_acc; + assert( R%Nevec_acc == 0 ); + Glog << "nBatch, Nevec_acc, R, Nu = " + << Nbatch << "," << Nevec_acc << "," << R << "," << Nu << std::endl; - cuDoubleComplex *w_acc, *evec_acc, *c_acc; - - cudaStat = cudaMallocManaged((void **)&w_acc, _Nu*sites*12*sizeof(cuDoubleComplex)); - Glog << cudaStat << std::endl; - cudaStat = cudaMallocManaged((void **)&evec_acc, _R*sites*12*sizeof(cuDoubleComplex)); - Glog << cudaStat << std::endl; - cudaStat = cudaMallocManaged((void **)&c_acc, _Nu*_R*12*sizeof(cuDoubleComplex)); - Glog << cudaStat << std::endl; - - Glog << "cuBLAS prepare array"<< std::endl; #if 0 // a trivial test - for (int col=0; col<_Nu; ++col) { + for (int col=0; col(&w_v[0]); for (size_t row=0; row(&evec_v[0]); - for (size_t row=0; rowGlobalSumVector((double*)c_acc,2*_Nu*_R); - cublasDestroy(handle); + Glog << "cuBLAS Zgemm"<< std::endl; + + for (int b=0; b(&evec_v[0]); + for (size_t row=0; rowGlobalSumVector((double*)c_acc,2*Nu*Nevec_acc); + + for (int i=0; i[" << j << "," << i << "] = " + << z.x << " + i " << z.y << std::endl; + } + w[i] = w[i] - ip * evec[b*Nevec_acc+j]; + } + //assert(normalize(w[i],do_print)!=0); + } + } + + for (int i=0; i[" << j << "," << i << "] = " - << z.x << " + i " << z.y << std::endl; - } - w[i] = w[i] - ip * evec[j]; - } - assert(normalize(w[i],_print)!=0); - } + cublasDestroy(handle); - cudaFree(w_acc); - cudaFree(evec_acc); - cudaFree(c_acc); - Glog << "cuBLAS orthogonalize done" << std::endl; #else Glog << "BLAS wrapper is not implemented" << std::endl; @@ -411,6 +401,21 @@ for( int i =0;i& evec, const std::vector& src, int& Nconv, LanczosType Impl) { +#ifdef GRID_NVCC + GridBase *grid = src[0].Grid(); + grid->show_decomposition(); + + // set eigenvector buffers for the cuBLAS calls + //const uint64_t nsimd = grid->Nsimd(); + const uint64_t sites = grid->lSites(); + + cudaStat = cudaMallocManaged((void **)&w_acc, Nu*sites*12*sizeof(cuDoubleComplex)); + //Glog << cudaStat << std::endl; + cudaStat = cudaMallocManaged((void **)&evec_acc, Nevec_acc*sites*12*sizeof(cuDoubleComplex)); + //Glog << cudaStat << std::endl; + cudaStat = cudaMallocManaged((void **)&c_acc, Nu*Nevec_acc*sizeof(cuDoubleComplex)); + //Glog << cudaStat << std::endl; +#endif switch (Impl) { case LanczosType::irbl: calc_irbl(eval,evec,src,Nconv); @@ -420,6 +425,12 @@ for( int i =0;i& eval, @@ -919,8 +930,8 @@ if(!split_test){ orthogonalize(w[u],evec,R); } #else - //orthogonalize(w,Nu,evec,R); - orthogonalize_blas(w,Nu,evec,R); + //orthogonalize(w,Nu,evec,R); + orthogonalize_blas(w,evec,R); #endif // QR part for (int u=1; u Date: Tue, 31 Mar 2020 05:39:31 -0400 Subject: [PATCH 5/5] use gemm twice to complete the Gram Schmidt --- .../ImplicitlyRestartedBlockLanczos.h | 23 ++++++++++++++++--- 1 file changed, 20 insertions(+), 3 deletions(-) diff --git a/Grid/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h b/Grid/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h index cd6bd682..95fedbca 100644 --- a/Grid/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h +++ b/Grid/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h @@ -303,7 +303,7 @@ public: //Glog << stat << std::endl; grid->GlobalSumVector((double*)c_acc,2*Nu*Nevec_acc); - +#if 0 for (int i=0; i(&w_v[0]); + for (size_t row=0; row