*** empty log message ***

davidson.h

@@ -1,9 +1,14 @@
+#include "vec.h"
+#include "smat.h"
+#include "mat.h"
+#include "sparsemat.h"
+#include "nonclass.h"
 
 //Davidson diagonalization of real symmetric matrix
 
-//matrix can be any class which provides nrows(), diagonalof() and operator*(const  NRVec<T>) methods
+//matrix can be any class which has nrows(), ncols(), diagonalof() and NRVec::gemv() available
 //does not even have to be explicitly stored - direct CI
 //n<0 highest eigenvalues, n>0 lowest eigenvalues
 
-template <class T, class matrix>
+export template <typename T, typename Matrix>
-extern void davidson(const matrix &a, NRVec<T> *vecs /*input-output*/, T *vals, const int n=1, const double eps=1e-10);
+extern void davidson(const Matrix &bigmat, NRVec<T> *eivecs /*input-output*/, NRVec<T> &eivals, int nroots=1, const double accur=1e-6, int nits=50, const int ndvdmx = 500, const bool incore=1);

mat.cc

@@ -34,6 +34,25 @@ NRMat<T> & NRMat<T>::operator=(const T &a)
 
 
 
+//get diagonal; for compatibility with large matrices do not return newly created object
+template <typename T>
+void NRMat<T>::diagonalof(NRVec<T> &r) const
+{
+#ifdef DEBUG
+        if (nn != mm) laerror("diagonalof() non-square matrix");
+	if (r.size() != nn) laerror("diagonalof() incompatible vector");
+#endif
+
+#ifdef MATPTR
+         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];}
+#endif
+
+}
+
+
+
 // M += a
 template <typename T>
 NRMat<T> & NRMat<T>::operator+=(const T &a)
@@ -750,7 +769,7 @@ INSTANTIZE(int)
 INSTANTIZE(char)
 
 
-export template <class T>
+export template <typename T>
 ostream& operator<<(ostream &s, const NRMat<T> &x)
                 {
                 int i,j,n,m;
@@ -764,7 +783,7 @@ ostream& operator<<(ostream &s, const NRMat<T> &x)
                 return s;
                 }
 
-export template <class T>
+export template <typename T>
 istream& operator>>(istream  &s, NRMat<T> &x)
                 {
                 int i,j,n,m;

mat.h

@@ -58,6 +58,7 @@ public:
 	const NRVec<T> operator*(const NRVec<T> &rhs) const; // Mat * Vec
 	const NRVec<T> rsum() const; //sum of rows
 	const NRVec<T> csum() const; //sum of columns
+	void diagonalof(NRVec<T> &) 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
@@ -348,21 +349,22 @@ NRMat<T>::~NRMat()
 template <typename T>
 NRMat<T> & NRMat<T>::operator=(const NRMat<T> &rhs)
 {
-        if (this == &rhs) return *this;
+        if (this !=&rhs)
-        if (count) {
+		{
-                if (--(*count) ==0 ) {
+        	if (count) 
+                    if (--(*count) ==0 ) {
 #ifdef MATPTR
                         delete[] (v[0]);
 #endif
                         delete[] v;
                         delete count;
-                }
+                	}
                 v = rhs.v;
                 nn = rhs.nn;
                 mm = rhs.mm;
                 count = rhs.count;
                 if (count) (*count)++;
-        }
+		}
         return *this;
 }
 

nonclass.cc

@@ -186,12 +186,14 @@ extern "C" void FORNAME(dsyev)(const char *JOBZ, const char *UPLO, const int *N,
 
 // a will contain eigenvectors (columns if corder==1), w eigenvalues
 void diagonalize(NRMat<double> &a, NRVec<double> &w, const bool eivec, 
-		const bool corder)
+		const bool corder, int n)
 {
-	int n = a.nrows();
+	int m = a.nrows();
-	if (n != a.ncols()) laerror("diagonalize() call with non-square matrix");
+	if (m != a.ncols()) laerror("diagonalize() call with non-square matrix");
 	if (a.nrows() != w.size()) 
 		laerror("inconsistent dimension of eigenvalue vector in diagonalize()");
+	if(n==0) n=m;
+	if(n<0||n>m) laerror("actual dimension out of range in diagonalize");
 
 	a.copyonwrite();
 	w.copyonwrite();
@@ -204,10 +206,10 @@ void diagonalize(NRMat<double> &a, NRVec<double> &w, const bool eivec,
 	double WORKX;
 
 	// First call is to determine size of workspace
-	FORNAME(dsyev)(&vectors, &U, &n, a, &n, w, (double *)&WORKX, &LWORK, &r );
+	FORNAME(dsyev)(&vectors, &U, &n, a, &m, w, (double *)&WORKX, &LWORK, &r );
 	LWORK = (int)WORKX;
 	double *WORK = new double[LWORK];
-	FORNAME(dsyev)(&vectors, &U, &n, a, &n, w, WORK, &LWORK, &r );
+	FORNAME(dsyev)(&vectors, &U, &n, a, &m, w, WORK, &LWORK, &r );
 	delete[] WORK;
 	if (vectors == 'V' && corder) a.transposeme();
 
@@ -216,6 +218,7 @@ void diagonalize(NRMat<double> &a, NRVec<double> &w, const bool eivec,
 }
 
 
+
 extern "C" void FORNAME(dspev)(const char *JOBZ, const char *UPLO, const int *N,
 		double *AP, double *W, double *Z, const int *LDZ, double *WORK, int *INFO);
 
@@ -526,18 +529,23 @@ double trace2(const NRSMat<double> &a, const NRSMat<double> &b,
 
 
 //generalized diagonalization, eivecs will be in columns of a
-void gendiagonalize(NRMat<double> &a, NRVec<double> &w, NRMat<double> b)
+//counts with actual dimension smaller than allocated dimension
+
+void gendiagonalize(NRMat<double> &a, NRVec<double> &w, NRMat<double> b, int n)
 {
 if(a.nrows()!=a.ncols() || a.nrows()!=w.size() || a.nrows()!=b.nrows()  || b.nrows()!=b.ncols() ) laerror("incompatible Mats in gendiagonalize");
 
 a.copyonwrite();
 w.copyonwrite();
 b.copyonwrite();
-int n=w.size();
+int m=w.size();
-NRVec<double> dl(n);
+NRVec<double> dl(m);
 int i,j;
 double x;
 
+if(n==0) n=m;
+if(n<0 || n>m) laerror("actual dimension in gendiagonalize out of range");
+
 //transform the problem to usual diagonalization
 /*
 c
@@ -620,7 +628,7 @@ for(j=0; j<n ; ++j)
 {
 	for(i=j;i<n;++i)
 	{
-        x = a(i,j) - cblas_ddot(i-j,&a(j,j),n,&b(i,j),1)
+        x = a(i,j) - cblas_ddot(i-j,&a(j,j),m,&b(i,j),1)
 		- cblas_ddot(j,&a(j,0),1,&b(i,0),1);
         a(i,j) = x/dl[i];
 	}
@@ -630,7 +638,7 @@ for(j=0; j<n ; ++j)
 for(i=1;i<n;++i) for(j=0; j<i; ++j) a(j,i)=a(i,j);
 
 //diagonalize by a standard procedure
-diagonalize(a,w,1,1);
+diagonalize(a,w,1,1,n);
 
 //transform the eigenvectors back
 /*
@@ -669,7 +677,7 @@ for(j=0; j<n; ++j)//eigenvector loop
 {
 for(int i=n-1; i>=0; --i)//component loop
 	{
-        if(i<n-1) a(i,j) -= cblas_ddot(n-1-i,&b(i+1,i),n,&a(i+1,j),n);
+        if(i<n-1) a(i,j) -= cblas_ddot(n-1-i,&b(i+1,i),m,&a(i+1,j),m);
         a(i,j) /= dl[i];
 	}
 }

nonclass.h

@@ -1,3 +1,6 @@
+#ifndef _LA_NONCLASS_H_
+#define _LA_NONCLASS_H_
+
 #include "vec.h"
 #include "smat.h"
 #include "mat.h"
@@ -48,7 +51,7 @@ extern T trace2(const NRMat<T> &a, const NRMat<T> &b, bool trb=0); \
 extern T trace2(const NRSMat<T> &a, const NRSMat<T> &b, const bool diagscaled=0);\
 extern void linear_solve(NRMat<T> &a, NRMat<T> *b, double *det=0); \
 extern void linear_solve(NRSMat<T> &a, NRMat<T> *b, double *det=0); \
-extern void diagonalize(NRMat<T> &a, NRVec<T> &w, const bool eivec=1, const bool corder=1); \
+extern void diagonalize(NRMat<T> &a, NRVec<T> &w, const bool eivec=1, const bool corder=1, int n=0); \
 extern void diagonalize(NRSMat<T> &a, NRVec<T> &w, NRMat<T> *v, const bool corder=1);\
 extern void singular_decomposition(NRMat<T> &a, NRMat<T> *u, NRVec<T> &s,\
 		NRMat<T> *v, const bool corder=1);
@@ -72,7 +75,7 @@ extern NRMat<double> matrixfunction(NRMat<double> a, complex<double> (*f)(const
 
 
 //generalized diagonalization of symmetric matrix with symmetric positive definite metric matrix b
-extern void gendiagonalize(NRMat<double> &a, NRVec<double> &w, NRMat<double> b);
+extern void gendiagonalize(NRMat<double> &a, NRVec<double> &w, NRMat<double> b, const int n=0);
 
 //functions on matrices
 inline NRMat<double>  sqrt(const NRSMat<double> &a) { return matrixfunction(a,&sqrt); }
@@ -98,3 +101,4 @@ const NRMat<T> inverse(NRMat<T> a, T *det=0)
 	return result;
 }
 
+#endif

smat.cc

@@ -34,7 +34,7 @@ NRSMat<T>::NRSMat(const NRMat<T> &rhs)
 
 
 
-// assing to diagonal
+// assign to diagonal
 template <typename T>
 NRSMat<T> & NRSMat<T>::operator=(const T &a)
 {
@@ -43,6 +43,17 @@ NRSMat<T> & NRSMat<T>::operator=(const T &a)
 	return *this;
 }
 
+//get diagonal
+template <typename T>
+void NRSMat<T>::diagonalof(NRVec<T> &r) const
+{
+#ifdef DEBUG
+if(r.size()!=nn) laerror("incompatible vector in diagonalof()");
+#endif
+        for (int i=0; i<nn; i++) r[i] = v[i*(i+1)/2+i];
+}
+
+
 // unary minus
 template <typename T>
 const NRSMat<T> NRSMat<T>::operator-() const
@@ -71,7 +82,7 @@ void NRSMat<T>::fprintf(FILE *file, const char *format, const int modulo) const
 }
 
 // read matrix from the file with specific format
-template <class T>
+template <typename T>
 void NRSMat<T>::fscanf(FILE *f, const char *format)
 {
 	int n, m;
@@ -268,7 +279,7 @@ void NRSMat< complex<double> >::axpy(const complex<double> alpha,
 }
 
 
-export template <class T>
+export template <typename T>
 ostream& operator<<(ostream &s, const NRSMat<T> &x)
                 {
                 int i,j,n;
@@ -282,7 +293,7 @@ ostream& operator<<(ostream &s, const NRSMat<T> &x)
                 }
 
 
-export template <class T>
+export template <typename T>
 istream& operator>>(istream  &s, NRSMat<T> &x)
                 {
                 int i,j,n,m;

smat.h

@@ -43,6 +43,7 @@ public:
 	const NRMat<T> operator*(const NRMat<T> &rhs) const; // SMat*Mat 
 	const T dot(const NRSMat &rhs) const; // Smat.Smat
 	const NRVec<T> operator*(const NRVec<T> &rhs) const; 
+	void diagonalof(NRVec<T> &) 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;

sparsemat.cc

@@ -484,6 +484,24 @@ for(i=0;i<nn;++i)
 		}
 }
 
+
+//get diagonal, do not construct a new object, but store in existing one, important for huge CI matrices
+template <class T>
+void SparseMat<T>::diagonalof(NRVec<T> &r) const
+{
+#ifdef DEBUG
+if(nn!=mm) laerror("diagonalof() non-square sparse matrix");
+#endif
+matel<T> *l=list;
+r=(T)0;
+while(l)
+	{
+	if(l->row == l->col) r[l->row]+= l->elem;
+	l= l->next;
+	}
+}
+
+
 //constructor dense matrix from sparse
 export template <class T>
 NRMat<T>::NRMat(const SparseMat<T> &rhs)

sparsemat.h

@@ -1,3 +1,6 @@
+#ifndef _SPARSEMAT_H_
+#define _SPARSEMAT_H_
+
 //for vectors and dense matrices we shall need
 #include "la.h"
 
@@ -85,6 +88,7 @@ public:
         inline const SparseMat operator-(const SparseMat &rhs) const {return SparseMat(*this) -= rhs;} //must not be symmetric+general
 	const NRVec<T> multiplyvector(const NRVec<T> &rhs, const bool transp=0) const; //sparse matrix * dense vector optionally transposed
 	inline const NRVec<T> operator*(const NRVec<T> &rhs) const {return multiplyvector(rhs);} //sparse matrix * dense vector
+	void diagonalof(NRVec<T> &) 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
@@ -219,4 +223,4 @@ while(l)
 return *this;
 }
 
-
+#endif

t.cc

@@ -6,6 +6,7 @@
 #include "sparsemat.h"
 #include "matexp.h"
 #include "fourindex.h"
+#include "davidson.h"
 
 
 extern void test(const NRVec<double> &x);
@@ -489,7 +490,7 @@ for(int n=8; n<=1024*1024;n+=n)
 	}
 }
 
-if(1)
+if(0)
 {
 int n;
 cin>>n;
@@ -536,7 +537,7 @@ if(0)
 {
 const int n=3;
 NRMat<double> a(n,n);
-for(int i=0;i<n;++i) for(int j=0;j<i;++j)
+for(int i=0;i<n;++i) for(int j=0;j<=i;++j)
 	{
 	a(i,j)= random()/(1.+RAND_MAX);
 	a(j,i)= -a(i,j);
@@ -772,7 +773,45 @@ cout <<"test "<<a.dot(a)<<endl;
 
 }
 
+if(0)
+{
+NRMat<double> amat,bmat;
+cin >>amat;
+cin >>bmat;
+NRVec<double> v(amat.nrows());
+gendiagonalize(amat,v,bmat,2);
+cout <<amat;
+cout <<v;
 }
 
 
+if(1)
+{
+int n,m;
+cin>>n >>m;
+NRMat<double> a(n,n);
+NRVec<double> rr(n);
 
+for(int i=0;i<n;++i) for(int j=0;j<=i;++j)
+        {
+        a(i,j)= random()/(1.+RAND_MAX);
+        a(j,i)= a(i,j);
+	if(i==j) a(i,i)+= .5*(i-n);
+        }
+
+NRMat<double> aa=a;
+cout <<aa;
+diagonalize(aa,rr);
+cout <<aa;
+cout <<rr;
+
+NRVec<double> r(m);
+NRVec<double> *eivecs = new NRVec<double>[m];
+davidson(a,eivecs,r,m);
+
+cout <<r;
+
+}
+
+
+}

test.cc

@@ -0,0 +1,8 @@
+#include "vec.h"
+
+int main(void)
+{
+NRVec<double> *p = new NRVec<double>[1000];
+NRVec<double> q(10); q=1.;
+p[500]|=q;
+}

vec.h

@@ -521,7 +521,7 @@ void NRVec<T>::resize(const int n)
 }
 
 
-// assignmet with a physical copy
+// assignmet with a physical (deep) copy
 template <typename T>
 NRVec<T> & NRVec<T>::operator|=(const NRVec<T> &rhs)
 {