*** empty log message ***

davidson.h

@@ -11,6 +11,8 @@
 
 //matrix can be any class which has nrows(), ncols(), diagonalof(), issymmetric(), and gemv() available
 //does not even have to be explicitly stored - direct CI
+//therefore the whole implementation must be a template in a header
+//Note that for efficiency in a direct CI case the diagonalof() should cache its result
 
 export template <typename T, typename Matrix>
 extern void davidson(const Matrix &bigmat, NRVec<T> &eivals, NRVec<T> *eivecs, const char *eivecsfile, 
@@ -62,10 +64,10 @@ int oldnroot;
 smallS=0;
 smallH=0;     
 //guess based on lowest diagonal element of the matrix
-bigmat.diagonalof(vec2);
+const T *diagonal = bigmat.diagonalof(vec2,false,true);
 vec1=0;
 {T t=1e100; int i,j;  
-for(i=0, j= -1; i<n; ++i) if(vec2[i]<t) {t=vec2[i]; j=i;}
+for(i=0, j= -1; i<n; ++i) if(diagonal[i]<t) {t=diagonal[i]; j=i;}
 vec1[j]=1;}
 
 //init Krylov matrices
@@ -170,11 +172,11 @@ for(j=0; j<=krylovsize; ++j)
 	if(!incore) s1->get(vec2,j);
         vec1.axpy(smallV(j,nroot),incore?v1[j]:vec2);
 	}
-bigmat.diagonalof(vec2);
+diagonal = bigmat.diagonalof(vec2,false,true);
 eival_n = r[nroot];
 for(i=0; i<n; ++i)
 	{
-	T denom = vec2[i] - eival_n;
+	T denom = diagonal[i] - eival_n;
         denom = denom<0?-max(0.1,abs(denom)):max(0.1,abs(denom));
         vec1[i] /= denom;
 	}

fourindex.h

@@ -12,6 +12,8 @@
 //it is actually not needed for the algorithms here, but may be useful for the
 //user of this class to keep this piece of information along with the data
 
+//when patient enough, make const_casts for piterators to have pbegin() const
+
 template<class I>
 union packed_index {
                 I packed[4];
@@ -565,28 +567,28 @@ for(p=rhs.begin(); p!= rhs.end(); ++p) (*this)(p->index.indiv.i,p->index.indiv.j
 template<class T, class DUMMY>
 T& fourindex_dense<twoelectronrealmullikan,T,DUMMY>::operator() (unsigned int i, unsigned int j, unsigned int k, unsigned int l)
 {
-unsigned long I = i>j? i*(i-1)/2+j-1 : j*(j-1)/2+i-1;
+unsigned long I = SMat_index_1(i,j);
-unsigned long J = k>l? k*(k-1)/2+l-1 : l*(l-1)/2+k-1;
+unsigned long J = SMat_index_1(k,l);
 //I,J act as indices of a NRSmat
 #ifdef DEBUG
      	if (*count != 1) laerror("lval (i,j,k,l) with count > 1 in fourindex_dense");
        	if (I<0 || I>=(unsigned long)nn || J<0 || J>=(unsigned long)nn) laerror("fourindex_dense index out of range");
        	if (!v) laerror("access to unallocated fourindex_dense");
 #endif
-return I>=J ? v[I*(I+1)/2+J] : v[J*(J+1)/2+I];
+return v[SMat_index(I,J)];
 }
 
 template<class T, class DUMMY>
 const T& fourindex_dense<twoelectronrealmullikan,T,DUMMY>::operator() (unsigned int i, unsigned int j, unsigned int k, unsigned int l) const
 {
-unsigned long I = i>j? i*(i-1)/2+j-1 : j*(j-1)/2+i-1;
+unsigned long I = SMat_index_1(i,j);
-unsigned long J = k>l? k*(k-1)/2+l-1 : l*(l-1)/2+k-1;
+unsigned long J = SMat_index_1(k,l);
 //I,J act as indices of a NRSmat
 #ifdef DEBUG
-       	if (I<0 || I>=nn || J<0 || J>=nn) laerror("fourindex_dense index out of range");
+       	if (I<0 || I>=(unsigned long)nn || J<0 || J>=(unsigned long)nn) laerror("fourindex_dense index out of range");
        	if (!v) laerror("access to unallocated fourindex_dense");
 #endif
-return I>=J ? v[I*(I+1)/2+J] : v[J*(J+1)/2+I];
+return v[SMat_index(I,J)];
 }
 
 

mat.cc

@@ -866,12 +866,12 @@ void NRMat<double>::gemm(const double &beta, const NRMat<double> &a,
 		const char transa, const NRMat<double> &b, const char transb, 
 		const double &alpha)
 {
-	int l(transa=='n'?a.nn:a.mm);
 	int k(transa=='n'?a.mm:a.nn);
-	int kk(transb=='n'?b.nn:b.mm);
-	int ll(transb=='n'?b.mm:b.nn);
 
 #ifdef DEBUG
+	int l(transa=='n'?a.nn:a.mm);
+	int kk(transb=='n'?b.nn:b.mm);
+	int ll(transb=='n'?b.mm:b.nn);
 	if (l!=nn || ll!=mm || k!=kk) laerror("incompatible matrices in Mat:gemm()");
 #endif
 	if (alpha==0.0 && beta==1.0) return;
@@ -890,12 +890,12 @@ void NRMat< complex<double> >::gemm(const complex<double> & beta,
 		const NRMat< complex<double> > & b, const char transb, 
 		const complex<double> & alpha)
 {
-	int l(transa=='n'?a.nn:a.mm);
 	int k(transa=='n'?a.mm:a.nn);
-	int kk(transb=='n'?b.nn:b.mm);
-	int ll(transb=='n'?b.mm:b.nn);
 
 #ifdef DEBUG
+	int l(transa=='n'?a.nn:a.mm);
+	int kk(transb=='n'?b.nn:b.mm);
+	int ll(transb=='n'?b.mm:b.nn);
 	if (l!=nn || ll!=mm || k!=kk) laerror("incompatible matrices in Mat:gemm()");
 #endif
 	if (alpha==CZERO && beta==CONE) return;
@@ -1013,7 +1013,7 @@ const complex<double> NRMat< complex<double> >::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 
 template<>
-void NRMat<double>::diagonalof(NRVec<double> &r, const bool divide) const
+const double * NRMat<double>::diagonalof(NRVec<double> &r, const bool divide, bool cache) const
 {
 #ifdef DEBUG
 	if (r.size() != nn) laerror("diagonalof() incompatible vector");
@@ -1047,6 +1047,7 @@ for (int i=0; i< mm; i++)
 	}
 }
 
+return divide?NULL:&r[0];
 }
 
 

mat.h

@@ -71,7 +71,8 @@ public:
 	const NRVec<T> csum() const; //sum of columns
 	const NRVec<T> row(const int i) const; //row of, efficient
 	const NRVec<T> column(const int j) const {NRVec<T> r(nn); for(int i=0; i<nn; ++i) r[i]= (*this)(i,j); return r;}; //column of, general but not very efficient
-	void diagonalof(NRVec<T> &, const bool divide=0) const; //get diagonal 
+	const T* diagonalof(NRVec<T> &, const bool divide=0, bool cache=false) const; //get diagonal 
+	void gemv(const T beta, NRVec<T> &r, const char trans, const T alpha, const NRVec<T> &x) const {r.gemv(beta,*this,trans,alpha,x);};
 	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

smat.cc

@@ -85,7 +85,7 @@ NRSMat<T> & NRSMat<T>::operator=(const T &a)
 
 //get diagonal
 template <typename T>
-void NRSMat<T>::diagonalof(NRVec<T> &r, const bool divide) const
+const T* NRSMat<T>::diagonalof(NRVec<T> &r, const bool divide, bool cache) const
 {
 #ifdef DEBUG
 if(r.size()!=nn) laerror("incompatible vector in diagonalof()");
@@ -97,6 +97,7 @@ if (divide)
 	for (int i=0; i<nn; i++) {T a =v[i*(i+1)/2+i]; if(a!=0.) r[i] /= a;}
 else
         for (int i=0; i<nn; i++) r[i] = v[i*(i+1)/2+i];
+return divide?NULL:&r[0];
 }
 
 
@@ -266,7 +267,31 @@ NRSMat< complex<double> >::dot(const NRSMat< complex<double> > &rhs) const
 	if (nn != rhs.nn) laerror("dot of incompatible SMat's");
 #endif
 	complex<double> dot;
-	cblas_zdotc_sub(nn, (void *)v, 1, (void *)rhs.v, 1, (void *)(&dot));
+	cblas_zdotc_sub(NN2, (void *)v, 1, (void *)rhs.v, 1, (void *)(&dot));
+	return dot;
+}
+
+
+template<>
+const double NRSMat<double>::dot(const NRVec<double> &rhs) const
+{
+#ifdef DEBUG
+	if (NN2 != rhs.nn) laerror("dot of incompatible SMat's");
+#endif
+	return cblas_ddot(NN2, v, 1, rhs.v, 1);
+}
+
+
+
+template<>
+const complex<double> 
+NRSMat< complex<double> >::dot(const NRVec< complex<double> > &rhs) const
+{
+#ifdef DEBUG
+	if (NN2 != rhs.nn) laerror("dot of incompatible SMat's");
+#endif
+	complex<double> dot;
+	cblas_zdotc_sub(NN2, (void *)v, 1, (void *)rhs.v, 1, (void *)(&dot));
 	return dot;
 }
 

smat.h

@@ -43,8 +43,10 @@ public:
 	const NRMat<T> operator*(const NRSMat &rhs) const; // SMat*SMat
 	const NRMat<T> operator*(const NRMat<T> &rhs) const; // SMat*Mat 
 	const T dot(const NRSMat &rhs) const; // Smat.Smat//@@@for complex do conjugate
+	const T dot(const NRVec<T> &rhs) const; //Smat(as vec).vec //@@@for complex do conjugate
 	const NRVec<T> operator*(const NRVec<T> &rhs) const {NRVec<T> result(nn); result.gemv((T)0,*this,'n',(T)1,rhs); return result;}; // Mat * Vec
-	void diagonalof(NRVec<T> &, const bool divide=0) const; //get diagonal
+	const T* diagonalof(NRVec<T> &, const bool divide=0, bool cache=false) const; //get diagonal
+	void gemv(const T beta, NRVec<T> &r, const char trans, const T alpha, const NRVec<T> &x) const {r.gemv(beta,*this,trans,alpha,x);};
 	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;
@@ -52,6 +54,7 @@ public:
 	inline int nrows() const;
 	inline int ncols() const;
 	inline int size() const;
+	inline bool transp(const int i, const int j) const {return i>j;} //this can be used for compact storage of matrices, which are actually not symmetric, but one triangle of them is redundant
 	const double norm(const T scalar=(T)0) const;
 	void axpy(const T alpha, const NRSMat &x); // this+= a*x
 	inline const T amax() const;
@@ -274,6 +277,18 @@ inline const T & NRSMat<T>::operator[](const int ij) const
 	return v[ij];
 }
 
+template<typename T>
+inline T SMat_index(T i, T j)
+{
+return i>=j ? i*(i+1)/2+j : j*(j+1)/2+i;
+}
+
+template<typename T>
+inline T SMat_index_1(T i, T j)
+{
+return i>j? i*(i-1)/2+j-1 : j*(j-1)/2+i-1;
+}
+
 // access the element, 2-dim array case
 template <typename T>
 inline T & NRSMat<T>::operator()(const int i, const int j)
@@ -283,7 +298,7 @@ inline T & NRSMat<T>::operator()(const int i, const int j)
 	if (i<0 || i>=nn || j<0 || j>=nn) laerror("SMat (i,j) out of range");
 	if (!v) laerror("(i,j) for unallocated Smat");
 #endif
-	return i>=j ? v[i*(i+1)/2+j] : v[j*(j+1)/2+i];
+	return v[SMat_index(i,j)];
 }
 template <typename T>
 inline const T & NRSMat<T>::operator()(const int i, const int j) const
@@ -292,7 +307,7 @@ inline const T & NRSMat<T>::operator()(const int i, const int j) const
 	if (i<0 || i>=nn || j<0 || j>=nn) laerror("SMat (i,j) out of range");
 	if (!v) laerror("(i,j) for unallocated Smat");
 #endif
-	return i>=j ? v[i*(i+1)/2+j] : v[j*(j+1)/2+i];
+	return v[SMat_index(i,j)];
 }
 
 // return the number of rows and columns
@@ -471,18 +486,10 @@ for(int i=0; i<rhs.nrows(); ++i)
 return r;
 }
 
-
-
-
-
-
 // I/O
 template <typename T> extern ostream& operator<<(ostream &s, const NRSMat<T> &x);
 template <typename T> extern istream& operator>>(istream  &s, NRSMat<T> &x);
 
-
-
-
 // generate operators: SMat + a, a + SMat, SMat * a
 NRVECMAT_OPER(SMat,+)
 NRVECMAT_OPER(SMat,-)
@@ -491,4 +498,35 @@ NRVECMAT_OPER(SMat,*)
 NRVECMAT_OPER2(SMat,+)
 NRVECMAT_OPER2(SMat,-)
 
+//optional indexing from 1
+//all possible constructors have to be given explicitly, other stuff is inherited
+//with exception of the operator() which differs
+template<typename T>
+class NRSMat_from1 : public NRSMat<T> {
+public:
+        NRSMat_from1(): NRSMat<T>() {};
+        explicit NRSMat_from1(const int n): NRSMat<T>(n) {};
+	NRSMat_from1(const NRSMat<T> &rhs): NRSMat<T>(rhs) {}; //be able to convert the parent class  transparently to this
+        NRSMat_from1(const T &a, const int n): NRSMat<T>(a,n) {}; 
+        NRSMat_from1(const T *a, const int n): NRSMat<T>(a,n) {}; 
+        explicit NRSMat_from1(const NRMat<T> &rhs): NRSMat<T>(rhs) {};
+        explicit NRSMat_from1(const NRVec<T> &rhs, const int n): NRSMat<T>(rhs,n) {};
+
+        inline const T& operator() (const int i, const int j) const 
+		{
+#ifdef DEBUG
+		if(i<=0||j<=0||i>nn||j>nn) laerror("index out of range in NRSMat_from1");
+#endif
+		return v[SMat_index_1(i,j)];
+		}
+        inline  T& operator() (const int i, const int j) 
+		{
+#ifdef DEBUG
+                if(i<=0||j<=0||i>nn||j>nn) laerror("index out of range in NRSMat_from1");
+#endif
+		return v[SMat_index_1(i,j)];
+		}
+};
+
+
 #endif /* _LA_SMAT_H_ */

sparsemat.cc

@@ -579,7 +579,7 @@ for(i=0;i<nn;++i)
 //get diagonal, do not construct a new object, but store in existing one, important for huge CI matrices
 // with the divide option is used as a preconditioner, another choice of preconditioner is possible
 template <class T>
-void SparseMat<T>::diagonalof(NRVec<T> &r, const bool divide) const
+const T* SparseMat<T>::diagonalof(NRVec<T> &r, const bool divide, bool cache) const
 {
 #ifdef DEBUG
 if((int)mm!=r.size()) laerror("incompatible vector size in diagonalof()");
@@ -607,6 +607,7 @@ if(divide)
 	for(unsigned int i=0; i<mm; ++i) if((*rr)[i]!=0.) r[i]/=(*rr)[i];
 	delete rr;
 	}
+return divide?NULL:&r[0];
 }
 
 
@@ -877,7 +878,7 @@ else
 
 
 template<class T>
-void NRVec<T>::gemv(const T beta, const SparseMat<T> &a, const char trans, const T alpha, const NRVec<T> &x)
+void NRVec<T>::gemv(const T beta, const SparseMat<T> &a, const char trans, const T alpha, const NRVec<T> &x, const bool treat_as_symmetric)
 {
 if((trans=='n'?a.ncols():a.nrows())!= (SPMatindex)x.size()) laerror("incompatible sizes in gemv");
 copyonwrite();
@@ -893,7 +894,7 @@ T *vec=x.v;
 
 if(alpha==(T)1)
 {
-        if(a.issymmetric())
+        if(a.issymmetric()||treat_as_symmetric)
         {
         while(l)
                 {
@@ -920,7 +921,7 @@ if(alpha==(T)1)
 }
 else
 {
-	if(a.issymmetric())
+	if(a.issymmetric()||treat_as_symmetric)
 	{
 	while(l)
 	        {
@@ -1350,7 +1351,7 @@ template void SparseMat<T>::axpy(const T alpha, const SparseMat<T> &x, const boo
 template const SparseMat<T> SparseMat<T>::operator*(const SparseMat<T> &rhs) const; \
 template const T SparseMat<T>::dot(const SparseMat<T> &rhs) const; \
 template void SparseMat<T>::gemm(const T beta, const SparseMat<T> &a, const char transa, const SparseMat<T> &b, const char transb, const T alpha); \
-template void NRVec<T>::gemv(const T beta, const SparseMat<T> &a, const char trans, const T alpha, const NRVec<T> &x);\
+template void NRVec<T>::gemv(const T beta, const SparseMat<T> &a, const char trans, const T alpha, const NRVec<T> &x, const bool treat_as_symmetric);\
 template void SparseMat<T>::permuterows(const NRVec<SPMatindex> &p);\
 template void SparseMat<T>::permutecolumns(const NRVec<SPMatindex> &p);\
 template void SparseMat<T>::permuteindices(const NRVec<SPMatindex> &p);\

sparsemat.h

@@ -79,7 +79,8 @@ public:
         inline const SparseMat operator-(const SparseMat &rhs) const {return SparseMat(*this) -= rhs;} //must not be symmetric+general
 	inline const NRVec<T> operator*(const NRVec<T> &rhs) const; // SparseMat * Vec
 	inline const NRMat<T> operator*(const NRMat<T> &rhs) const; // SparseMat * Mat
-	void diagonalof(NRVec<T> &, const bool divide=0) const; //get diagonal
+	const T* diagonalof(NRVec<T> &, const bool divide=0, bool cache=false) const; //get diagonal
+	void gemv(const T beta, NRVec<T> &r, const char trans, const T alpha, const NRVec<T> &x, bool treat_as_symmetric=false) const {r.gemv(beta,*this,trans,alpha,x,treat_as_symmetric);};
 	const SparseMat operator*(const SparseMat &rhs) const; 
 	SparseMat & oplusequal(const SparseMat &rhs); //direct sum
 	SparseMat & oplusequal(const NRMat<T> &rhs);

test.cc

@@ -1,5 +1,7 @@
 #include "bitvector.h"
 #include "qsort.h"
+#include "la.h"
+#include "fourindex.h"
 
 
 int main(void)
@@ -31,4 +33,10 @@ ptrqsortup(&t[0],&t[9],&u[0]);
 cout<<t <<"U= "<<u;
 ptrqsortdown<int,int>(&t[0],&t[9]);
 cout<<t;
+
+NRSMat_from1<double> a(5),b(5),c;
+c=a+b;
+
+fourindex<int,double> f;
+fourindex_dense<twoelectronrealmullikan,double,int> ff(f);
 }

vec.cc

@@ -377,7 +377,7 @@ laerror("not yet implemented");
 template<>
 void NRVec<int>::gemv(const int beta,
                 const SparseMat<int> &A, const char trans,
-                const int alpha, const NRVec &x)
+                const int alpha, const NRVec &x, bool s)
 {
 laerror("not yet implemented");
 }
@@ -385,7 +385,7 @@ laerror("not yet implemented");
 template<>
 void NRVec<short>::gemv(const short beta,
                 const SparseMat<short> &A, const char trans,
-                const short alpha, const NRVec &x)
+                const short alpha, const NRVec &x, bool s)
 {
 laerror("not yet implemented");
 }
@@ -394,7 +394,7 @@ laerror("not yet implemented");
 template<>
 void NRVec<char>::gemv(const char beta,
                 const SparseMat<char> &A, const char trans,
-                const char alpha, const NRVec &x)
+                const char alpha, const NRVec &x, bool s)
 {
 laerror("not yet implemented");
 }
@@ -402,7 +402,7 @@ laerror("not yet implemented");
 template<>
 void NRVec<unsigned long>::gemv(const unsigned long beta,
                 const SparseMat<unsigned long> &A, const char trans,
-                const unsigned long alpha, const NRVec &x)
+                const unsigned long alpha, const NRVec &x, bool s)
 {
 laerror("not yet implemented");
 }
@@ -410,7 +410,7 @@ laerror("not yet implemented");
 template<>
 void NRVec<unsigned char>::gemv(const unsigned char beta,
                 const SparseMat<unsigned char> &A, const char trans,
-                const unsigned char alpha, const NRVec &x)
+                const unsigned char alpha, const NRVec &x, bool s)
 {
 laerror("not yet implemented");
 }

vec.h

@@ -62,6 +62,8 @@ public:
 	const NRVec operator-() const;
 	inline NRVec & operator+=(const NRVec &rhs);
 	inline NRVec & operator-=(const NRVec &rhs);
+	inline NRVec & operator*=(const NRVec &rhs); //elementwise
+	inline NRVec & operator/=(const NRVec &rhs); //elementwise
 	inline NRVec & operator+=(const T &a);
 	inline NRVec & operator-=(const T &a);
 	inline NRVec & operator*=(const T &a);
@@ -75,7 +77,7 @@ public:
 	inline const T dot(const NRVec &rhs) const {return *this * rhs;}; //@@@for complex do conjugate
 	void gemv(const T beta, const NRMat<T> &a, const char trans, const T alpha, const NRVec &x);
 	void gemv(const T beta, const NRSMat<T> &a, const char trans /*just for compatibility*/,  const T alpha, const NRVec &x);
-	void gemv(const T beta, const SparseMat<T> &a, const char trans, const T alpha, const NRVec &x);
+	void gemv(const T beta, const SparseMat<T> &a, const char trans, const T alpha, const NRVec &x,const bool treat_as_symmetric=false);
 	const NRVec operator*(const NRMat<T> &mat) const {NRVec<T> result(mat.ncols()); result.gemv((T)0,mat,'t',(T)1,*this); return result;};
 	const NRVec operator*(const NRSMat<T> &mat) const {NRVec<T> result(mat.ncols()); result.gemv((T)0,mat,'t',(T)1,*this); return result;};
 	const NRVec operator*(const SparseMat<T> &mat) const {NRVec<T> result(mat.ncols()); result.gemv((T)0,mat,'t',(T)1,*this); return result;};
@@ -253,6 +255,32 @@ inline NRVec<T> & NRVec<T>::operator+=(const NRVec<T> &rhs)
         return *this;
 }
 
+//for general type only
+template <typename T>
+inline NRVec<T> & NRVec<T>::operator*=(const NRVec<T> &rhs)
+{
+#ifdef DEBUG
+        if (nn != rhs.nn) laerror("*= of incompatible vectors");
+#endif
+        copyonwrite();
+        int i;
+        for(i=0; i<nn; ++i) v[i]*=rhs.v[i];
+        return *this;
+}
+
+template <typename T>
+inline NRVec<T> & NRVec<T>::operator/=(const NRVec<T> &rhs)
+{
+#ifdef DEBUG
+        if (nn != rhs.nn) laerror("/= of incompatible vectors");
+#endif
+        copyonwrite();
+        int i;
+        for(i=0; i<nn; ++i) v[i]/=rhs.v[i];
+        return *this;
+}
+
+
 
 // x -= x
 template<>
@@ -568,7 +596,7 @@ void NRVec<T>::resize(const int n)
 }
 
 
-// assignmet with a physical (deep) copy
+// assignment with a physical (deep) copy
 template <typename T>
 NRVec<T> & NRVec<T>::operator|=(const NRVec<T> &rhs)
 {