portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2025-04-10 06:00:45 +01:00
Code Releases Activity

Commenting and clean up

Browse Source
This commit is contained in:
paboyle 2017-04-10 20:38:20 +09:00
parent d80d802f9d
commit b12dc89d26
2 changed files with 21 additions and 2 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

20
lib/algorithms/iterative/BlockConjugateGradient.h
View File

@ -69,8 +69,14 @@ static void sliceMaddMatrix (Lattice<vobj> &R,Eigen::MatrixXcd &aa,const Lattice
Lattice<vobj> Xslice(SliceGrid); Lattice<vobj> Xslice(SliceGrid);
Lattice<vobj> Rslice(SliceGrid); Lattice<vobj> Rslice(SliceGrid);
// FIXME: Implementation is slow
// If we based this on Cshift it would work for spread out // If we based this on Cshift it would work for spread out
// but it would be even slower // but it would be even slower
//
// Repeated extract slice is inefficient
//
// Best base the linear combination by constructing a
// set of vectors of size grid->_rdimensions[Orthog].
for(int i=0;i<Nblock;i++){ for(int i=0;i<Nblock;i++){
ExtractSlice(Rslice,Y,i,Orthog); ExtractSlice(Rslice,Y,i,Orthog);
for(int j=0;j<Nblock;j++){ for(int j=0;j<Nblock;j++){
@ -84,6 +90,9 @@ template<class vobj>
static void sliceMaddVector (Lattice<vobj> &R,std::vector<RealD> &a,const Lattice<vobj> &X,const Lattice<vobj> &Y, static void sliceMaddVector (Lattice<vobj> &R,std::vector<RealD> &a,const Lattice<vobj> &X,const Lattice<vobj> &Y,
int Orthog,RealD scale=1.0) int Orthog,RealD scale=1.0)
{ {
// FIXME: Implementation is slow
// Best base the linear combination by constructing a
// set of vectors of size grid->_rdimensions[Orthog].
typedef typename vobj::scalar_object sobj; typedef typename vobj::scalar_object sobj;
typedef typename vobj::scalar_type scalar_type; typedef typename vobj::scalar_type scalar_type;
typedef typename vobj::vector_type vector_type; typedef typename vobj::vector_type vector_type;
@ -107,6 +116,8 @@ static void sliceMaddVector (Lattice<vobj> &R,std::vector<RealD> &a,const Lattic
template<class vobj> template<class vobj>
static void sliceInnerProductMatrix( Eigen::MatrixXcd &mat, const Lattice<vobj> &lhs,const Lattice<vobj> &rhs,int Orthog) static void sliceInnerProductMatrix( Eigen::MatrixXcd &mat, const Lattice<vobj> &lhs,const Lattice<vobj> &rhs,int Orthog)
{ {
// FIXME: Implementation is slow
// Not sure of best solution.. think about it
typedef typename vobj::scalar_object sobj; typedef typename vobj::scalar_object sobj;
typedef typename vobj::scalar_type scalar_type; typedef typename vobj::scalar_type scalar_type;
typedef typename vobj::vector_type vector_type; typedef typename vobj::vector_type vector_type;
@ -141,6 +152,9 @@ static void sliceInnerProductMatrix( Eigen::MatrixXcd &mat, const Lattice<vobj>
template<class vobj> template<class vobj>
static void sliceInnerProductVector( std::vector<ComplexD> & vec, const Lattice<vobj> &lhs,const Lattice<vobj> &rhs,int Orthog) static void sliceInnerProductVector( std::vector<ComplexD> & vec, const Lattice<vobj> &lhs,const Lattice<vobj> &rhs,int Orthog)
{ {
// FIXME: Implementation is slow
// Look at localInnerProduct implementation,
// and do inside a site loop with block strided iterators
typedef typename vobj::scalar_object sobj; typedef typename vobj::scalar_object sobj;
typedef typename vobj::scalar_type scalar_type; typedef typename vobj::scalar_type scalar_type;
typedef typename vobj::vector_type vector_type; typedef typename vobj::vector_type vector_type;
@ -257,6 +271,10 @@ void operator()(LinearOperatorBase<Field> &Linop, const Field &Src, Field &Psi)
Field AP(Src); Field AP(Src);
Field R(Src); Field R(Src);
GridStopWatch LinalgTimer;
GridStopWatch MatrixTimer;
GridStopWatch SolverTimer;
Eigen::MatrixXcd m_pAp = Eigen::MatrixXcd::Identity(Nblock,Nblock); Eigen::MatrixXcd m_pAp = Eigen::MatrixXcd::Identity(Nblock,Nblock);
Eigen::MatrixXcd m_pAp_inv= Eigen::MatrixXcd::Identity(Nblock,Nblock); Eigen::MatrixXcd m_pAp_inv= Eigen::MatrixXcd::Identity(Nblock,Nblock);
Eigen::MatrixXcd m_rr = Eigen::MatrixXcd::Zero(Nblock,Nblock); Eigen::MatrixXcd m_rr = Eigen::MatrixXcd::Zero(Nblock,Nblock);
@ -339,8 +357,10 @@ void operator()(LinearOperatorBase<Field> &Linop, const Field &Src, Field &Psi)
} }
if ( max_resid < Tolerance*Tolerance ) { if ( max_resid < Tolerance*Tolerance ) {
std::cout << GridLogMessage<<" Block solver has converged in " std::cout << GridLogMessage<<" Block solver has converged in "
<<k<<" iterations; max residual is "<<std::sqrt(max_resid)<<std::endl; <<k<<" iterations; max residual is "<<std::sqrt(max_resid)<<std::endl;
for(int b=0;b<Nblock;b++){ for(int b=0;b<Nblock;b++){
std::cout << GridLogMessage<< " block "<<b<<" resid "<< std::sqrt(real(m_rr(b,b))/ssq[b])<<std::endl; std::cout << GridLogMessage<< " block "<<b<<" resid "<< std::sqrt(real(m_rr(b,b))/ssq[b])<<std::endl;
} }

3
lib/algorithms/iterative/ConjugateGradient.h
View File

@ -99,8 +99,7 @@ class ConjugateGradient : public OperatorFunction<Field> {
} }
std::cout << GridLogIterative << std::setprecision(4) std::cout << GridLogIterative << std::setprecision(4)
<< "ConjugateGradient: k=0 residual " << cp << " target " << rsq << "ConjugateGradient: k=0 residual " << cp << " target " << rsq << std::endl;
<< std::endl;
GridStopWatch LinalgTimer; GridStopWatch LinalgTimer;
GridStopWatch MatrixTimer; GridStopWatch MatrixTimer;