*** empty log message ***

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<T> &x, const int pos) const;
 	void put(const NRVec<T> &x, const int pos);
         void get(NRMat<T> &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<typename T>
+AuxStorage<T>::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<typename T>
 void AuxStorage<T>::get(NRVec<T> &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<typename T>
+void AuxStorage<T>::get(NRMat<T> &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<typename T>
+void AuxStorage<T>::put(const NRMat<T> &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<typename T>
+void AuxStorage<T>::get(NRSMat<T> &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<typename T>
+void AuxStorage<T>::put(const NRSMat<T> &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

davidson.h

@@ -10,6 +10,6 @@
 //does not even have to be explicitly stored - direct CI
 
 export template <typename T, typename Matrix>
-extern void davidson(const Matrix &bigmat, NRVec<T> &eivals, void *eivecs, 
+extern void davidson(const Matrix &bigmat, NRVec<T> &eivals, NRVec<T> *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);

gmres.h

@@ -10,4 +10,4 @@
 //does not even have to be explicitly stored
 
 template<typename T, typename Matrix>
-extern void gmres(const Matrix &bigmat, const NRVec<T> &b, NRVec<T> &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<T> &b, NRVec<T> &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);

mat.h

@@ -496,6 +496,7 @@ void NRMat<T>::resize(const int n, const int m)
 
 
 
+
 // I/O
 template <typename T> extern ostream& operator<<(ostream &s, const NRMat<T> &x);
 template <typename T> extern istream& operator>>(istream  &s, NRMat<T> &x);

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<class T>
 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<class MAT>
-const typename NRMat_traits<MAT>::elementtype determinant(MAT a)//again passed by value
-{
-typename NRMat_traits<MAT>::elementtype det;
-if(a.nrows()!=a.ncols()) laerror("determinant of non-square matrix");
-linear_solve(a,NULL,&det);
-return det;
-}
 
 
 template<class M, class V>
@@ -257,3 +251,5 @@ for(int i=1; i<=(1<<power); ++i) //unfortunatelly, here we have to repeat it man
 
 return result;
 }
+
+#endif

nonclass.cc

@@ -608,3 +608,31 @@ for(int i=n-1; i>=0; --i)//component loop
 }
 
 }
+
+//auxiliary routine to adjust eigenvectors to guarantee real logarithm
+//at the moment not rigorous yet
+void adjustphases(NRMat<double> &v)
+{
+int n=v.nrows();
+double det=determinant(v);
+int nchange=0;
+for(int i=0; i<n;++i) if(v[i][i]<0.)
+        {
+        cblas_dscal(n,-1.,v[i],1); 
+        nchange++;
+        } 
+if(det<0) nchange++;
+if(nchange&1)//still adjust to get determinant=1
+        {
+        int imin=-1; double min=1e200;
+        for(int i=0; i<n;++i)
+                if(abs(v[i][i])<min)
+                        {
+                        imin=i;
+                        min=abs(v[i][i]);
+                        } 
+        cblas_dscal(n,-1.,v[imin],1);
+        }
+}
+
+

nonclass.h

@@ -4,6 +4,7 @@
 #include "vec.h"
 #include "smat.h"
 #include "mat.h"
+#include "la_traits.h"
 
 
 //MISC
@@ -104,4 +105,18 @@ const NRMat<T> inverse(NRMat<T> a, T *det=0)
 	return result;
 }
 
+//general determinant
+template<class MAT>
+const typename NRMat_traits<MAT>::elementtype determinant(MAT a)//again passed by value
+{
+typename NRMat_traits<MAT>::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<double> &v);
+
 #endif

t.cc

@@ -805,11 +805,11 @@ for(int i=0;i<n;++i) for(int j=0;j<=i;++j)
 	if(i==j) a(i,i)+= .5*(i-n);
         }
 
-NRMat<double> aa=a;
+NRMat<double> aa;
-diagonalize(aa,rr);
+aa=a; diagonalize(aa,rr);
 NRVec<double> r(m);
 NRVec<double> *eivecs = new NRVec<double>[m];
-davidson(a,r,eivecs,m,1);
+davidson(a,r,eivecs,NULL,m,1,1e-6,1,200);
 
 cout <<"Davidson energies " <<r;
 cout <<"Exact energies " ;
@@ -824,12 +824,26 @@ for(int i=0; i<m; ++i)
 	cout <<endl;
 	}
 
-
-
 }
 
+//davidson of large very sparse matrix (10n/n^2)
 
-if(1)
+#undef sparsity
+#define sparsity (n/2)
+if(0)
+{
+int n,m;
+cin >>n>>m;
+	SparseMat<double> aa(n,n);
+	aa.setsymmetric();
+	for(int i=0; i<sparsity;i++) aa.add(randind(n),randind(n),random()/(1.+RAND_MAX));
+	for(int i=0; i<n; ++i) aa.add(i,i,500*random()/(1.+RAND_MAX));
+NRVec<double> r(m);
+davidson(aa,r,(NRVec<double> *)NULL,"eivecs",m,1,1e-5,0,300,300);
+cout <<r;
+}
+
+if(0)
 {
 int n,m;
 cin>>n >>m;
@@ -851,8 +865,8 @@ NRVec<double> bb=b;
 NRMat<double>aa=a;
 linear_solve(aa,bb);
 //cout <<bb;
-//gmres   (a,b,x,1,1e-10,100,1,0,1,0);
+gmres   (a,b,x,1,1e-10,100,1,0,1,0);
-  conjgrad(a,b,x,1,1e-10,200,1,0,1);
+//conjgrad(a,b,x,1,1e-10,200,1,0,1);
 
 cout <<"\nsolution compare:\n";
 for(int i=0; i<m; ++i) 
@@ -864,9 +878,25 @@ for(int i=0; i<m; ++i)
 	cout <<endl;
 	}
 
-
-
 }
 
 
+if(1)
+{
+int n,m;
+cin>>n >>m;
+SparseMat<double> aa(n,m);
+NRVec<double> b(n);
+NRVec<double> x(m);
+
+//tridiagonal
+        for(int i=0; i<n; ++i) aa.add(i,i,random()/(1.+RAND_MAX));
+	for(int i=0; i<n-1; i+=1) aa.add(i,i+1,.002*random()/(1.+RAND_MAX));
+	for(int i=0; i<n-1; i+=1) aa.add(i+1,i,.002*random()/(1.+RAND_MAX));
+
+for(int i=0;i<n;++i) b[i] = i+1;
+gmres(aa,b,x,1,1e-20,20,1,1,1,1000,1);
+//conjgrad(aa,b,x,1,1e-10,1000,1,0,1);
+}
+
 }