From 6e2727f59577dd25fecfcfa3b8e0c480fd3a562c Mon Sep 17 00:00:00 2001 From: jiri Date: Fri, 4 Feb 2005 14:31:42 +0000 Subject: [PATCH] *** empty log message *** --- conjgrad.h | 14 +++++++++++ gmres.h | 2 +- laerror.cc | 16 ++++++++++++ laerror.h | 14 +++++++++++ mat.cc | 20 +++++++++++---- mat.h | 2 +- nonclass.cc | 17 ------------- nonclass.h | 3 +++ smat.cc | 8 +++++- smat.h | 2 +- sparsemat.cc | 18 +++++++++++--- sparsemat.h | 2 +- t.cc | 69 ++++++++++++++++++++++++++++++++++++++++++++++------ vec.h | 3 +-- 14 files changed, 150 insertions(+), 40 deletions(-) create mode 100644 conjgrad.h create mode 100644 laerror.cc create mode 100644 laerror.h diff --git a/conjgrad.h b/conjgrad.h new file mode 100644 index 0000000..d16ca62 --- /dev/null +++ b/conjgrad.h @@ -0,0 +1,14 @@ +#include "vec.h" +#include "smat.h" +#include "mat.h" +#include "sparsemat.h" +#include "nonclass.h" + +//conjugate gradient solution of a linear system + +//matrix can be any class which has nrows(), ncols(), diagonalof() and NRVec::gemv() available +//does not even have to be explicitly stored + + +template +extern void conjgrad(const Matrix &bigmat, const NRVec &b, NRVec &x, const bool doguess, const double tol, const int itmax, const bool verbose, bool issquare,const bool precondition); diff --git a/gmres.h b/gmres.h index e4dc6e6..04356a1 100644 --- a/gmres.h +++ b/gmres.h @@ -10,4 +10,4 @@ //does not even have to be explicitly stored template -extern void gmres(const Matrix &bigmat, const NRVec &b, NRVec &x, const bool doguess=1, const double eps=1e-7, const int MAXIT=50, int neustart=0, const bool verbose=1, bool square=1,const bool precondition=1); +extern void gmres(const Matrix &bigmat, const NRVec &b, NRVec &x, const bool doguess=1, const double eps=1e-7, const int MAXIT=50, const bool verbose=1, bool square=1,const bool precondition=1, int neustart=0); diff --git a/laerror.cc b/laerror.cc new file mode 100644 index 0000000..168c7ef --- /dev/null +++ b/laerror.cc @@ -0,0 +1,16 @@ +// LA errorr handler +#include +#include "laerror.h" +void laerror(const char *s1) +{ + std::cerr << "LA:ERROR - "; + if(!s1) + std::cerr << "udefined.\n"; + else { + if(s1) std::cerr << s1; + std::cerr << "\n"; + } + throw LAerror(s1); + exit(1); +} + diff --git a/laerror.h b/laerror.h new file mode 100644 index 0000000..b78f8cb --- /dev/null +++ b/laerror.h @@ -0,0 +1,14 @@ +#ifndef _LAERROR_H_ +#define _LAERROR_H_ + +//exception class for laerror +class LAerror + { + const char *msg; + public: + LAerror(const char *s) {msg=s;}; + }; + +extern void laerror(const char *); + +#endif diff --git a/mat.cc b/mat.cc index 1007f0d..2958b30 100644 --- a/mat.cc +++ b/mat.cc @@ -743,28 +743,38 @@ const complex NRMat< complex >::trace() const //get diagonal; for compatibility with large matrices do not return newly created object //for non-square get diagonal of A^TA, will be used as preconditioner -void NRMat::diagonalof(NRVec &r) const +void NRMat::diagonalof(NRVec &r, const bool divide) const { #ifdef DEBUG if (r.size() != nn) laerror("diagonalof() incompatible vector"); #endif +double a; + +r.copyonwrite(); + if(nn==mm) { #ifdef MATPTR - for (int i=0; i< nn; i++) r[i] = v[i][i]; +if(divide) for (int i=0; i< nn; i++) if((a=v[i][i])) r[i]/=a; +else for (int i=0; i< nn; i++) r[i] = v[i][i]; #else - {int i,j; for (i=j=0; j< nn; ++j, i+=nn+1) r[j] = v[i];} +if(divide) {int i,j; for (i=j=0; j< nn; ++j, i+=nn+1) if((a=v[i])) r[j] /=a;} +else {int i,j; for (i=j=0; j< nn; ++j, i+=nn+1) r[j] = v[i];} #endif } else //non-square { for (int i=0; i< mm; i++) + { #ifdef MATPTR - r[i] = cblas_ddot(nn,v[0]+i,mm,v[0]+i,mm); + a= cblas_ddot(nn,v[0]+i,mm,v[0]+i,mm); #else - r[i] = cblas_ddot(nn,v+i,mm,v+i,mm); + a=cblas_ddot(nn,v+i,mm,v+i,mm); #endif + if(divide) {if(a) r[i]/=a;} + else r[i] = a; + } } } diff --git a/mat.h b/mat.h index ecd5310..d463a02 100644 --- a/mat.h +++ b/mat.h @@ -58,7 +58,7 @@ public: const NRVec operator*(const NRVec &rhs) const; // Mat * Vec const NRVec rsum() const; //sum of rows const NRVec csum() const; //sum of columns - void diagonalof(NRVec &) const; //get diagonal + void diagonalof(NRVec &, const bool divide=0) const; //get diagonal inline T* operator[](const int i); //subscripting: pointer to row i inline const T* operator[](const int i) const; inline T& operator()(const int i, const int j); // (i,j) subscripts diff --git a/nonclass.cc b/nonclass.cc index a5e5e60..7a7959e 100644 --- a/nonclass.cc +++ b/nonclass.cc @@ -98,23 +98,6 @@ void lawritemat(FILE *file,const T *a,int r,int c,const char *form0, } } -// LA errorr handler -void laerror(const char *s1, const char *s2, const char *s3, const char *s4) -{ - std::cerr << "LA:ERROR - "; - if(!s1) - std::cerr << "udefined."; - else { - if(s1) std::cerr << s1; - if(s2) std::cerr << s2; - if(s3) std::cerr << s3; - if(s4) std::cerr << s4; - } - std::cerr << endl; -//@@@throw - exit(1); -} - ////////////////////// // LAPACK interface // ////////////////////// diff --git a/nonclass.h b/nonclass.h index fba68b6..959c467 100644 --- a/nonclass.h +++ b/nonclass.h @@ -5,6 +5,7 @@ #include "smat.h" #include "mat.h" + //MISC export template const NRMat diagonalmatrix(const NRVec &x) @@ -51,6 +52,8 @@ extern T trace2(const NRMat &a, const NRMat &b, bool trb=0); \ extern T trace2(const NRSMat &a, const NRSMat &b, const bool diagscaled=0);\ extern void linear_solve(NRMat &a, NRMat *b, double *det=0); \ extern void linear_solve(NRSMat &a, NRMat *b, double *det=0); \ +extern void linear_solve(NRMat &a, NRVec &b, double *det=0); \ +extern void linear_solve(NRSMat &a, NRVec &b, double *det=0); \ extern void diagonalize(NRMat &a, NRVec &w, const bool eivec=1, const bool corder=1, int n=0); \ extern void diagonalize(NRSMat &a, NRVec &w, NRMat *v, const bool corder=1);\ extern void singular_decomposition(NRMat &a, NRMat *u, NRVec &s,\ diff --git a/smat.cc b/smat.cc index 1e04085..404c0a9 100644 --- a/smat.cc +++ b/smat.cc @@ -45,11 +45,17 @@ NRSMat & NRSMat::operator=(const T &a) //get diagonal template -void NRSMat::diagonalof(NRVec &r) const +void NRSMat::diagonalof(NRVec &r, const bool divide) const { #ifdef DEBUG if(r.size()!=nn) laerror("incompatible vector in diagonalof()"); #endif + +r.copyonwrite(); + +if (divide) + for (int i=0; i operator*(const NRMat &rhs) const; // SMat*Mat const T dot(const NRSMat &rhs) const; // Smat.Smat const NRVec operator*(const NRVec &rhs) const; - void diagonalof(NRVec &) const; //get diagonal + void diagonalof(NRVec &, const bool divide=0) const; //get diagonal inline const T& operator[](const int ij) const; inline T& operator[](const int ij); inline const T& operator()(const int i, const int j) const; diff --git a/sparsemat.cc b/sparsemat.cc index 5c2c789..cd631cb 100644 --- a/sparsemat.cc +++ b/sparsemat.cc @@ -486,26 +486,36 @@ for(i=0;i -void SparseMat::diagonalof(NRVec &r) const +void SparseMat::diagonalof(NRVec &r, const bool divide) const { #ifdef DEBUG if((int)mm!=r.size()) laerror("incompatible vector size in diagonalof()"); #endif matel *l=list; -r=(T)0; +NRVec *rr; + +r.copyonwrite(); +if(divide) {rr=new NRVec(mm); *rr=(T)0;} +else {r=(T)0; rr=&r;} if(nn==mm) //square while(l) { - if(l->row == l->col) r[l->row]+= l->elem; + if(l->row == l->col) (*rr)[l->row]+= l->elem; l= l->next; } else //diagonal of A^TA, assuming it has been simplified (only one entry per position), will be used as preconditioner only anyway while(l) { - r[l->col] += l->elem*l->elem *(l->col!=l->row && symmetric?2.:1.); + (*rr)[l->col] += l->elem*l->elem *(l->col!=l->row && symmetric?2.:1.); l= l->next; } +if(divide) + { + for(unsigned int i=0; i multiplyvector(const NRVec &rhs, const bool transp=0) const; //sparse matrix * dense vector optionally transposed inline const NRVec operator*(const NRVec &rhs) const {return multiplyvector(rhs);} //sparse matrix * dense vector - void diagonalof(NRVec &) const; //get diagonal + void diagonalof(NRVec &, const bool divide=0) const; //get diagonal const SparseMat operator*(const SparseMat &rhs) const; void gemm(const T beta, const SparseMat &a, const char transa, const SparseMat &b, const char transb, const T alpha);//this := alpha*op( A )*op( B ) + beta*this, if this is symemtric, only half will be added onto it const T dot(const SparseMat &rhs) const; //supervector dot product diff --git a/t.cc b/t.cc index 262ed7b..6824ba5 100644 --- a/t.cc +++ b/t.cc @@ -7,6 +7,8 @@ #include "matexp.h" #include "fourindex.h" #include "davidson.h" +#include "gmres.h" +#include "conjgrad.h" extern void test(const NRVec &x); @@ -42,6 +44,10 @@ NRVec x(1.,10); NRVec y(2.,10); NRVec z(-2.,10); +cout.setf(ios::fixed); +cout.precision(12); + + if(0) test(x); y.axpy(3,x); @@ -785,7 +791,7 @@ cout <>n >>m; @@ -800,16 +806,65 @@ for(int i=0;i aa=a; -cout < r(m); NRVec *eivecs = new NRVec[m]; -davidson(a,eivecs,r,m); +davidson(a,r,eivecs,m,1); + +cout <<"Davidson energies " <>n >>m; +NRMat a(n,m); +NRVec b(n); +NRVec x(m); + +for(int i=0;i bb=b; +//cout <aa=a; +linear_solve(aa,bb); +//cout < class SparseMat; ////////////////////////////////////////////////////////////////////////////// // Forward declarations -void laerror(const char *s1=0, const char *s2=0, const char *s3=0, const char *s4=0); template void lawritemat(FILE *file,const T *a,int r,int c, const char *form0,int nodim,int modulo, int issym);