From ac8afe89ad73734e4d0ff4b15ba8b074505f8045 Mon Sep 17 00:00:00 2001 From: jiri Date: Sun, 6 Feb 2005 14:01:27 +0000 Subject: [PATCH] *** empty log message *** --- auxstorage.h | 56 +++++++++++++++++++++++++++++++++++++++++++++++++++- davidson.h | 2 +- gmres.h | 2 +- mat.h | 1 + matexp.h | 14 +++++-------- nonclass.cc | 28 ++++++++++++++++++++++++++ nonclass.h | 15 ++++++++++++++ t.cc | 50 ++++++++++++++++++++++++++++++++++++---------- 8 files changed, 146 insertions(+), 22 deletions(-) diff --git a/auxstorage.h b/auxstorage.h index c8d668c..422747d 100644 --- a/auxstorage.h +++ b/auxstorage.h @@ -21,10 +21,12 @@ class AuxStorage { char filename[32]; int fd; + bool deleteme; size_t recl; public: AuxStorage(void); - ~AuxStorage(void) {close(fd); unlink(filename);}; + AuxStorage(const char *name); + ~AuxStorage(void) {close(fd); if(deleteme) unlink(filename);}; void get(NRVec &x, const int pos) const; void put(const NRVec &x, const int pos); void get(NRMat &x, const int pos) const; @@ -42,8 +44,22 @@ unlink(filename); fd=open(filename,O_CREAT|O_LARGEFILE|O_RDWR,0777); if(fd<0) {perror(""); laerror("open failed in AuxStorage");} recl=0; +deleteme=1; } +template +AuxStorage::AuxStorage(const char *name) +{ +strcpy(filename,name); +unlink(filename); +fd=open(filename,O_CREAT|O_LARGEFILE|O_RDWR,0777); +if(fd<0) {perror(""); laerror("open failed in AuxStorage");} +recl=0; +deleteme=0; +} + + +//vectors template void AuxStorage::get(NRVec &x, const int pos) const { @@ -61,5 +77,43 @@ if((off64_t)-1 == lseek64(fd,pos*((off64_t)recl),SEEK_SET)) {perror(""); laerror if(0>write(fd,&x[0],recl)) {perror(""); laerror("write failed in AuxStorage");} } +//matrices +template +void AuxStorage::get(NRMat &x, const int pos) const +{ +if(recl==0) laerror("get from an empty file in AuxStorage"); +if((off64_t)-1 == lseek64(fd,pos*((off64_t)recl),SEEK_SET)) {perror(""); laerror("seek failed in AuxStorage");} +if((ssize_t)recl!=read(fd,&x(0,0),recl)) {perror(""); laerror("read failed in AuxStorage");} +} + +template +void AuxStorage::put(const NRMat &x, const int pos) +{ +if(recl) {if(recl!=x.nrows()*x.ncols()*sizeof(T)) laerror("attempt to write objects of different size to one AuxStorage");} +else recl=x.nrows()*x.ncols()*sizeof(T); +if((off64_t)-1 == lseek64(fd,pos*((off64_t)recl),SEEK_SET)) {perror(""); laerror("seek failed in AuxStorage");} +if(0>write(fd,&x(0,0),recl)) {perror(""); laerror("write failed in AuxStorage");} +} + +//packed symmetric matrices +template +void AuxStorage::get(NRSMat &x, const int pos) const +{ +if(recl==0) laerror("get from an empty file in AuxStorage"); +if((off64_t)-1 == lseek64(fd,pos*((off64_t)recl),SEEK_SET)) {perror(""); laerror("seek failed in AuxStorage");} +if((ssize_t)recl!=read(fd,&x(0,0),recl)) {perror(""); laerror("read failed in AuxStorage");} +} + +template +void AuxStorage::put(const NRSMat &x, const int pos) +{ +if(recl) {if(recl!=x.nrows()*(x.nrows()+1)/2*sizeof(T)) laerror("attempt to write objects of different size to one AuxStorage");} +else recl=x.nrows()*(x.nrows()+1)/2*sizeof(T); +if((off64_t)-1 == lseek64(fd,pos*((off64_t)recl),SEEK_SET)) {perror(""); laerror("seek failed in AuxStorage");} +if(0>write(fd,&x(0,0),recl)) {perror(""); laerror("write failed in AuxStorage");} +} + + + #endif diff --git a/davidson.h b/davidson.h index 5840b3f..0103ad8 100644 --- a/davidson.h +++ b/davidson.h @@ -10,6 +10,6 @@ //does not even have to be explicitly stored - direct CI export template -extern void davidson(const Matrix &bigmat, NRVec &eivals, void *eivecs, +extern void davidson(const Matrix &bigmat, NRVec &eivals, NRVec *eivecs, const char *eivecsfile, int nroots=1, const bool verbose=0, const double eps=1e-6, const bool incore=1, int maxit=100, const int maxkrylov = 500); diff --git a/gmres.h b/gmres.h index 04356a1..bc90540 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, const bool verbose=1, bool square=1,const bool precondition=1, int neustart=0); +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, const int incore=1); diff --git a/mat.h b/mat.h index d463a02..3225205 100644 --- a/mat.h +++ b/mat.h @@ -496,6 +496,7 @@ void NRMat::resize(const int n, const int m) + // I/O template extern ostream& operator<<(ostream &s, const NRMat &x); template extern istream& operator>>(istream &s, NRMat &x); diff --git a/matexp.h b/matexp.h index ffb9bfc..9b0cc71 100644 --- a/matexp.h +++ b/matexp.h @@ -1,3 +1,5 @@ +#ifndef _MATEXP_H_ +#define _MATEXP_H_ //general routine for polynomial of a matrix, tuned to minimize the number //of matrix-matrix multiplications on cost of additions and memory // the polynom and exp routines will work on any type, for which traits class @@ -123,7 +125,7 @@ template const T ipow( const T &x, int i) { if(i<0) laerror("negative exponent in ipow"); -if(i==0) {T r=x; r=1.; return r;}//trick for matrix dimension +if(i==0) {T r=x; r=(T)1; return r;}//trick for matrix dimension if(i==1) return x; T y,z; z=x; @@ -223,14 +225,6 @@ return r; } -template -const typename NRMat_traits::elementtype determinant(MAT a)//again passed by value -{ -typename NRMat_traits::elementtype det; -if(a.nrows()!=a.ncols()) laerror("determinant of non-square matrix"); -linear_solve(a,NULL,&det); -return det; -} template @@ -257,3 +251,5 @@ for(int i=1; i<=(1<=0; --i)//component loop } } + +//auxiliary routine to adjust eigenvectors to guarantee real logarithm +//at the moment not rigorous yet +void adjustphases(NRMat &v) +{ +int n=v.nrows(); +double det=determinant(v); +int nchange=0; +for(int i=0; i inverse(NRMat a, T *det=0) return result; } +//general determinant +template +const typename NRMat_traits::elementtype determinant(MAT a)//again passed by value +{ +typename NRMat_traits::elementtype det; +if(a.nrows()!=a.ncols()) laerror("determinant of non-square matrix"); +linear_solve(a,NULL,&det); +return det; +} + + +//auxiliary routine to adjust eigenvectors to guarantee real logarithm +extern void adjustphases(NRMat &v); + #endif diff --git a/t.cc b/t.cc index 6824ba5..7eaf0fe 100644 --- a/t.cc +++ b/t.cc @@ -805,11 +805,11 @@ for(int i=0;i aa=a; -diagonalize(aa,rr); +NRMat aa; +aa=a; diagonalize(aa,rr); NRVec r(m); NRVec *eivecs = new NRVec[m]; -davidson(a,r,eivecs,m,1); +davidson(a,r,eivecs,NULL,m,1,1e-6,1,200); cout <<"Davidson energies " <>n>>m; + SparseMat aa(n,n); + aa.setsymmetric(); + for(int i=0; i r(m); +davidson(aa,r,(NRVec *)NULL,"eivecs",m,1,1e-5,0,300,300); +cout <>n >>m; @@ -851,8 +865,8 @@ NRVec bb=b; NRMataa=a; linear_solve(aa,bb); //cout <>n >>m; +SparseMat aa(n,m); +NRVec b(n); +NRVec x(m); + +//tridiagonal + for(int i=0; i