bugfixes in sparsemat and sparsesmat

sparsemat.cc

@@ -44,8 +44,9 @@ static inline void SWAP(T &a, T &b)
 
 
 template <class T>
-void SparseMat<T>::get(int fd, bool dimen, bool transp)
+void SparseMat<T>::get(int fd, bool dimen, bool transp, bool orcaformat)
 {
+if(orcaformat) laerror("ORCA format not supported for SparseMat");
 errno=0;
 SPMatindex dim[2];
 if(dimen)
@@ -81,8 +82,9 @@ list=l;
 
 
 template <class T>
-void SparseMat<T>::put(int fd,bool dimen, bool transp) const
+void SparseMat<T>::put(int fd,bool dimen, bool transp, bool orcaformat) const
 {
+if(orcaformat) laerror("ORCA format not supported for SparseMat");
 errno=0;
 if(dimen)
 {
@@ -1378,8 +1380,8 @@ return *this;
 template SparseMat<T> & SparseMat<T>::oplusequal(const SparseMat<T> &rhs);\
 template SparseMat<T> & SparseMat<T>::oplusequal(const NRMat<T> &rhs);\
 template SparseMat<T> & SparseMat<T>::oplusequal(const NRSMat<T> &rhs);\
-template void SparseMat<T>::get(int fd, bool dimen, bool transp); \
+template void SparseMat<T>::get(int fd, bool dimen, bool transp, bool orcaformat); \
-template void SparseMat<T>::put(int fd, bool dimen, bool transp) const; \
+template void SparseMat<T>::put(int fd, bool dimen, bool transp, bool orcaformat) const; \
 template void SparseMat<T>::copyonwrite(); \
 template void SparseMat<T>::unsort(); \
 template void SparseMat<T>::resize(const SPMatindex n, const SPMatindex m); \

sparsemat.h

@@ -118,8 +118,8 @@ public:
 	void setlist(matel<T> *l) {list=l;}
 	inline SPMatindex nrows() const {return nn;}
         inline SPMatindex ncols() const {return mm;}
-	void get(int fd, bool dimensions=1, bool transposed=false);
+	void get(int fd, bool dimensions=1, bool transposed=false, bool orcaformat=false); //orcaformat not supported (dense only), for LA_traits
-	void put(int fd, bool dimensions=1, bool transposed=false) const;
+	void put(int fd, bool dimensions=1, bool transposed=false, bool orcaformat=false) const;
 	void resize(const SPMatindex n, const SPMatindex m); //destructive
         void dealloc(void) {resize(0,0);}
 	void incsize(const SPMatindex n, const SPMatindex m); //increase size without destroying the data
@@ -353,7 +353,7 @@ public:
 	NRVec<T> diagonal;
 	SparseMat<T> offdiagonal;
         SemiSparseMat() {};
-        SemiSparseMat(const SPMatindex n, const SPMatindex m) : offdiagonal(n,m) {if(n==m) diagonal.resize(n);};
+        SemiSparseMat(const SPMatindex n, const SPMatindex m) : offdiagonal(n,m) {if(n==m) {diagonal.resize(n); diagonal.clear();}}; //zero matrix, add() accumulates into the diagonal
         SPMatindex nrows() const {return offdiagonal.nrows();}
         SPMatindex ncols() const {return offdiagonal.ncols();}
         SPMatindex issymmetric() const {return offdiagonal.issymmetric();}
@@ -381,7 +381,7 @@ public:
 	void add(const SPMatindex n, const SPMatindex m, const T elem) {if(diagonal.size()>0 && n==m) diagonal[n]+=elem; else offdiagonal.add(n,m,elem);}
 	void clear() {diagonal.clear(); offdiagonal.clear();}
 	void copyonwrite() {diagonal.copyonwrite(); offdiagonal.copyonwrite();}
-	void resize(const SPMatindex n, const SPMatindex m) {offdiagonal.resize(n,m); diagonal.resize(n==m?n:0);}
+	void resize(const SPMatindex n, const SPMatindex m) {offdiagonal.resize(n,m); diagonal.resize(n==m?n:0); if(n==m) diagonal.clear();} //zero matrix as SparseMat::resize
 };
 
 template <class T>

sparsesmat.cc

@@ -29,35 +29,23 @@ namespace LA {
 
 
 //dense times sparse (not necessarily symmetric)
+//result row i accumulated as sum_k A_ik (row k of rhs), skipping zero A_ik
 template <typename T>
 SparseSMat<T>  NRMat<T>::operator*(const SparseSMat<T> &rhs) const
 {
+if(mm!=(int)rhs.nrows())  laerror("incompatible sizes in NRMat*SparseSMat");
 SparseSMat<T> r(nn,rhs.ncols());
-if(mm!=rhs.nrows())  laerror("incompatible sizes in NRMat*SparseSMat");
+for(int i=0; i<nn; ++i)
-for(SPMatindex k=0; k<mm; ++k) //summation loop
 	{
-	std::map<SPMatindex,T> * kl = rhs.line(k);
+	std::map<SPMatindex,T> *ri=NULL;
-	if(kl)
+	for(int k=0; k<mm; ++k)
 		{
-		//gather the data
+		const T aik=(*this)(i,k);
-		typename std::map<SPMatindex,T>::iterator p;
+		const std::map<SPMatindex,T> *kl = rhs.line(k);
-		int i,j;
+		if(aik==(T)0 || !kl) continue;
-		NRVec<T> kline(kl->size());
+		if(!ri) ri = r[i] = new std::map<SPMatindex,T>;
-		NRVec<SPMatindex> klineind(kl->size());
+		typename std::map<SPMatindex,T>::const_iterator p;
-		for(p=kl->begin(), i=0; p!=kl->end(); ++p,++i)
+		for(p=kl->begin(); p!=kl->end(); ++p) (*ri)[p->first] += aik * p->second; //new map element is value-initialized to 0
-       		         {
-       		         klineind[i] = p->first;
-       		         kline[i] = p->second;
-       		         }
-		NRVec<T> kcol = column(k);
-
-		//multiply
-		NRMat<T> prod=kcol.otimes(kline,false,1.);
-
-		//scatter the results
-		for(i=0; i<prod.nrows(); ++i) for(j=0; j<prod.ncols(); ++j)
-                	add(i,klineind[j],prod(i,j),false);
-
 		}
 	}
 r.simplify();
@@ -76,7 +64,7 @@ if(a.nn!=b.nn || a.nn!=nn) laerror("incompatible sizes in SparseSMat::gemm");
 copyonwrite();
 
 for(SPMatindex k=0; k<nn; ++k) //summation loop
-    if(a.v[k] && b.v[k]) //nonempty in both
+    if(a.v[k] && b.v[k] && !a.v[k]->empty() && !b.v[k]->empty()) //nonempty in both (allocated rows can be empty, e.g. after cholesky)
 	{
 	NRVec<T> av(a.v[k]->size());
 	NRVec<T> bv(b.v[k]->size());
@@ -148,7 +136,7 @@ void SparseSMat<T>::gemv(const T beta, NRVec<T> &r, const char trans, const T al
 {
 if(nn!=r.size() || mm!= x.size()) laerror("incompatible matrix vector dimensions in SparseSMat::gemv");
 if(tolower(trans)!='n') laerror("transposition not implemented yet in SparseSMat::gemv");
-r *= beta;
+if(beta==(T)0) r.clear(); else if(beta!=(T)1) r *= beta; //no multiplication by 0, r may be uninitialized
 if(alpha == (T)0) return;
 r.copyonwrite();
 for(SPMatindex i=0; i<nn; ++i) if(v[i])
@@ -283,29 +271,34 @@ void SparseSMat<T>::storesubmatrix(const int fromrow, const int fromcol, const S
         const int torow = fromrow + rhs.nrows() - 1;
 #ifdef DEBUG
         if(fromrow<0 || fromrow>=nn || torow>=nn || fromcol<0 || fromcol>=mm || tocol>=mm) laerror("bad indices in storesubmatrix");
+	if(&rhs==this) laerror("storesubmatrix from itself");
 #endif
+	copyonwrite();
+	//store, not add: remove the old elements of the block first
+	for(int i=fromrow; i<=torow; ++i) if(v[i]) v[i]->erase(v[i]->lower_bound((SPMatindex)fromcol), v[i]->upper_bound((SPMatindex)tocol));
 	typename SparseSMat<T>::iterator p(rhs);
 	for(; p.notend(); ++p) add(p->row+fromrow, p->col+fromcol, p->elem, false);
 }
 
 
 template <class T>
-void SparseSMat<T>::get(int fd, bool dimen, bool transp) {
+void SparseSMat<T>::get(int fd, bool dimen, bool transp, bool orcaformat) {
+  if(orcaformat) laerror("ORCA format not supported for SparseSMat");
   errno=0;
   SPMatindex dim[2];
 
 if(dimen) {
     if(2*sizeof(SPMatindex)!=read(fd,&dim,2*sizeof(SPMatindex))) laerror("read() error in SparseSMat::get ");
-    resize(dim[0],dim[1]);
+    if(transp) resize(dim[1],dim[0]); else resize(dim[0],dim[1]);
   }
-else  copyonwrite(); 
+else  clear(); //the read matrix replaces the old contents
 
 do {
     if(2*sizeof(SPMatindex)!=read(fd,&dim,2*sizeof(SPMatindex))) laerror("read() error 2 in SparseSMat::get");
     if(dim[0]==(SPMatindex) -1 || dim[1]==(SPMatindex) -1) break;
     typename LA_traits_io<T>::IOtype tmp;
     LA_traits<T>::get(fd,tmp,dimen,transp); // general way to work when elem is some complex class again
-    if(transp) add(dim[0],dim[1],tmp,false);  else add(dim[1],dim[0],tmp,false); 
+    if(transp) add(dim[1],dim[0],tmp,false);  else add(dim[0],dim[1],tmp,false); 
   } 
 while(1);
 }
@@ -313,17 +306,18 @@ while(1);
 
 
 template <class T>
-void SparseSMat<T>::put(int fd, bool dimen, bool transp) const {
+void SparseSMat<T>::put(int fd, bool dimen, bool transp, bool orcaformat) const {
+  if(orcaformat) laerror("ORCA format not supported for SparseSMat");
   errno=0;  
   if(dimen) {
-    if(sizeof(SPMatindex)!=write(fd,&nn,sizeof(SPMatindex))) laerror("cannot write() 1 in SparseSMat::put");
+    if(sizeof(SPMatindex)!=write(fd,&(transp ? mm : nn),sizeof(SPMatindex))) laerror("cannot write() 1 in SparseSMat::put");
-    if(sizeof(SPMatindex)!=write(fd,&mm,sizeof(SPMatindex))) laerror("cannot write() 2 in SparseSMat::put");
+    if(sizeof(SPMatindex)!=write(fd,&(transp ? nn : mm),sizeof(SPMatindex))) laerror("cannot write() 2 in SparseSMat::put");
   }
   
   typename SparseSMat<T>::iterator p(*this);
   for(; p.notend(); ++p) {
-    if(sizeof(SPMatindex)!=write(fd,&(p->row),sizeof(SPMatindex))) laerror("cannot write() 3 in SparseSMat::put");
+    if(sizeof(SPMatindex)!=write(fd,&(transp ? p->col : p->row),sizeof(SPMatindex))) laerror("cannot write() 3 in SparseSMat::put");
-    if(sizeof(SPMatindex)!=write(fd,&(p->col),sizeof(SPMatindex))) laerror("cannot write() 4 in SparseSMat::put");
+    if(sizeof(SPMatindex)!=write(fd,&(transp ? p->row : p->col),sizeof(SPMatindex))) laerror("cannot write() 4 in SparseSMat::put");
     typename LA_traits_io<T>::IOtype tmp = p->elem;
     LA_traits<T>::put(fd,tmp,dimen,transp); // general way to work when elem is some non-scalar class again
   }
@@ -347,8 +341,8 @@ template SparseSMat<T> & SparseSMat<T>::operator+=(const T &a); \
 template SparseSMat<T> & SparseSMat<T>::operator-=(const T &a); \
 template LA_traits<T>::normtype SparseSMat<T>::norm(const T scalar) const; \
 template const T* SparseSMat<T>::diagonalof(NRVec<T> &r, const bool divide, bool cache) const; \
-template void SparseSMat<T>::get(int fd, bool dimen, bool transp); \
+template void SparseSMat<T>::get(int fd, bool dimen, bool transp, bool orcaformat); \
-template void SparseSMat<T>::put(int fd, bool dimen, bool transp) const; \
+template void SparseSMat<T>::put(int fd, bool dimen, bool transp, bool orcaformat) const; \
 
 
 INSTANTIZE(double)
@@ -358,10 +352,12 @@ INSTANTIZE(complex<double>)
 //// forced instantization of functions in the header in the corresponding object file
 template class SparseSMat<double>;
 template class SparseSMat<std::complex<double> >;
+template SparseSMat<double> NRMat<double>::operator*(const SparseSMat<double> &rhs) const;
+template SparseSMat<std::complex<double> > NRMat<std::complex<double> >::operator*(const SparseSMat<std::complex<double> > &rhs) const;
 
 /*activate this if needed
-template void SparseSMat<NRMat<double> >::put(int fd, bool dimen, bool transp) const;
+template void SparseSMat<NRMat<double> >::put(int fd, bool dimen, bool transp, bool orcaformat) const;
-template void SparseSMat<NRMat<double> >::get(int fd, bool dimen, bool transp);
+template void SparseSMat<NRMat<double> >::get(int fd, bool dimen, bool transp, bool orcaformat);
 */
 
 

sparsesmat.h

@@ -40,9 +40,15 @@
 
 namespace LA {
 
-//symmetric sparse matrix class with a representation for efficient exponentiatiation
+//sparse matrix stored by rows, each row a std::map from column index to value, for efficient exponentiation and products
-//in particular we need a unitary symmetric complex matrix as exp(iH) with H real symmetric
+//originally written for symmetric matrices, in particular the unitary complex symmetric exp(iH) with H real symmetric,
+//but it need not be symmetric nor square: elements (i,j) and (j,i) are stored independently and operator*, gemv,
+//transpose and the conversions work for general matrices
+//the "S" remains in some conventions: add() with both=true (default) inserts (i,j) and conj at (j,i) (hermitian for complex),
+//as do the constructor from NRSMat and otimes_sparse; gemm assumes A=A^T (A^dagger); issymmetric() returns always true (for davidson)
 //indices are counted from zero
+//TODO: handle it as a general matrix class with a member symmetry_type (0 none, 1 symmetric, 2 hermitian);
+//the conversion from NRSMat, otimes_sparse, gemm (and issymmetric) would then be implemented according to it, in a backward compatible way
 
 
 template <typename T>
@@ -67,7 +73,7 @@ public:
         void resize(const SPMatindex nn, const SPMatindex mm);
    	void dealloc(void) {resize(0,0);}
 	inline void setcoldim(int i) {mm=(SPMatindex)i;};
-	//std::map<SPMatindex,T> *line(SPMatindex n) const {return v[n];};
+	inline const std::map<SPMatindex,T> *line(const SPMatindex n) const {return v[n];}; //read-only row access
 	typedef std::map<SPMatindex,T> *ROWTYPE;
 	inline typename SparseSMat<T>::ROWTYPE & operator[](const SPMatindex i) {return v[i];}; 
         void clear() {resize(nn,mm);}
@@ -90,15 +96,15 @@ public:
 	void gemv(const T beta, NRVec<T> &r, const char trans, const T alpha, const NRVec<T> &x) const;
 	inline const NRVec<T> operator*(const NRVec<T> &rhs) const {NRVec<T> result(nn); this->gemv((T)0,result,'n',(T)1,rhs); return result;};
 	typename LA_traits<T>::normtype norm(const T scalar=(T)0) const;
-        inline const SparseSMat operator*(const SparseSMat &rhs) const {SparseSMat<T> r(nn,mm); r.gemm(0,*this,'n',rhs,'n',1); return r;}; //!!!NOT A GENERAL ROUTINE, JUST FOR THE CASES WHEN THE RESULT STAYS SYMMETRIC
+        const SparseSMat operator*(const SparseSMat &rhs) const; //general product A B of the stored matrices (hermitian, complex symmetric or general)
-	void gemm(const T beta, const SparseSMat &a, const char transa, const SparseSMat &b, const char transb, const T alpha); //this := alpha*op( A )*op( B ) + beta*this !!!NOT A GENERAL ROUTINE, JUST FOR THE CASES WHEN THE RESULT STAYS SYMMETRIC
+	void gemm(const T beta, const SparseSMat &a, const char transa, const SparseSMat &b, const char transb, const T alpha); //this := alpha*A^T*op( B ) + beta*this (A^dagger if transa=='c'), i.e. A*op(B) only for A=A^T (A=A^dagger) !!!NOT A GENERAL ROUTINE
 	inline void add(const SPMatindex n, const SPMatindex m, const T elem, const bool both=true);
 	inline unsigned long long length() {return simplify();};
 	void transposeme() const {laerror("in-place transposition not necessary/implemented for SparseSMat");};
 	SparseSMat transpose(bool conj=false) const; //if we store a non-symmetric matrix there
         inline bool issymmetric() const {return true;} // LV: for davidson
-	void get(int fd, bool dimen, bool transp);
+	void get(int fd, bool dimen=1, bool transp=0, bool orcaformat=false); //orcaformat not supported (dense only), for LA_traits
-        void put(int fd, bool dimen, bool transp) const;
+        void put(int fd, bool dimen=1, bool transp=0, bool orcaformat=false) const;
 	int nrows() const {return nn;}
 	int ncols() const {return mm;}
 	SparseSMat<T>  cholesky(void) const;
@@ -197,12 +203,27 @@ int i,j;
 for(i=0; i<nn; ++i) for(j=0; j<=i; ++j) if(std::abs(rhs(i,j))>SPARSEEPSILON) (*this).add(i,j,rhs(i,j),true);
 }
 
+//from SparseMat: duplicates summed; a symmetric SparseMat (one triangle stored, complex symmetric) is unfolded without conjugation
+template <typename T>
+SparseSMat<T>::SparseSMat(const SparseMat<T> &rhs)
+:nn(rhs.nrows()), mm(rhs.ncols()),
+count(new int(1))
+{
+v= new std::map<SPMatindex,T> * [nn];
+memset(v,0,nn*sizeof(std::map<SPMatindex,T> *));
+const bool sym=rhs.issymmetric();
+for(const matel<T> *l=rhs.getlist(); l; l=l->next)
+	{
+	add(l->row,l->col,l->elem,false);
+	if(sym && l->row!=l->col) add(l->col,l->row,l->elem,false);
+	}
+}
+
 template <typename T>
 SparseSMat<T>::SparseSMat(const NRMat<T> &rhs)
 :nn(rhs.nrows()), mm(rhs.ncols()),
 count(new int(1))
 {
-if(rhs.nrows()!=rhs.ncols()) laerror("non-square matrix in SparseSMat constructor from NRMat");
 v= new std::map<SPMatindex,T> * [nn];
 memset(v,0,nn*sizeof(std::map<SPMatindex,T> *));
 int i,j;
@@ -234,7 +255,7 @@ count = new int(1);
 v=new T[nn2];
 memset(v,0,nn2*sizeof(T));
 typename SparseSMat<T>::iterator p(rhs);
-for(; p.notend(); ++p) if(p->row <= p->col) (*this)(p->row,p->col)=p->elem;
+for(; p.notend(); ++p) if(p->row >= p->col) (*this)(p->row,p->col)=p->elem; //lower triangle as stored in NRSMat
 }
 #undef nn2
 
@@ -307,15 +328,14 @@ if(*count > 1) //it was shared
   }
 else  //it was not shared
 	{
-	mm=m;
 	//delete all trees
 	for(SPMatindex i=0; i<nn; ++i) if(v[i]) {delete v[i]; v[i]=NULL;}
-	if(n!=nn)
+	if(n==0) {delete[] v; v=NULL; delete count; count=NULL; nn=mm=0; return;} //as in the shared case
-		{
+	if(n>nn) {delete[] v; v= new std::map<SPMatindex,T> * [n];} //shrinking keeps the longer array
 	nn=n;
+	mm=m;
 	for(SPMatindex i=0; i<nn; ++i) v[i]=NULL;
 	}
-	}
 }
 
 
@@ -363,6 +383,27 @@ void SparseSMat<T>::copyonwrite()
 }
 
 
+//general product A B row by row: row i of the result is sum_k A_ik (row k of B)
+template <typename T>
+const SparseSMat<T> SparseSMat<T>::operator*(const SparseSMat<T> &rhs) const
+{
+if(mm!=rhs.nn) laerror("incompatible sizes in SparseSMat::operator*");
+SparseSMat<T> r(nn,rhs.mm);
+for(SPMatindex i=0; i<nn; ++i) if(v[i])
+	{
+	std::map<SPMatindex,T> *ri=NULL;
+	typename std::map<SPMatindex,T>::const_iterator p,q;
+	for(p=v[i]->begin(); p!=v[i]->end(); ++p) if(rhs.v[p->first])
+		{
+		if(!ri) ri = r.v[i] = new std::map<SPMatindex,T>;
+		for(q=rhs.v[p->first]->begin(); q!=rhs.v[p->first]->end(); ++q)
+			(*ri)[q->first] += p->second * q->second; //new map element is value-initialized to 0
+		}
+	}
+r.simplify();
+return r;
+}
+
 template <typename T>
 void SparseSMat<T>::add(const SPMatindex n, const SPMatindex m, const T elem, const bool both)
 {
@@ -375,11 +416,11 @@ typename std::map<SPMatindex,T>::iterator p;
 
 p= v[n]->find(m);
 if(p!=v[n]->end()) p->second+=elem; else (*v[n])[m] = elem;
-if(n!=m && both) //add also transposed
+if(n!=m && both) //add also transposed, conjugated (hermitian for complex)
 	{
 	if(!v[m]) v[m] = new std::map<SPMatindex,T>;
 	p= v[m]->find(n);
-	if(p!=v[m]->end()) p->second+=elem; else (*v[m])[n] = elem;
+	if(p!=v[m]->end()) p->second+=LA_traits<T>::conjugate(elem); else (*v[m])[n] = LA_traits<T>::conjugate(elem);
 	}
 //addition can lead to zero, but do not treat it now, make a simplify
 }
@@ -585,6 +626,7 @@ return r;
 template<typename T>
 SparseSMat<T> otimes_sparse(const NRVec<T> &lhs, const NRVec<T> &rhs, const bool conjugate=false, const T &scale=1) 
 {
+if(lhs.size()!=rhs.size()) laerror("otimes_sparse requires vectors of equal length");
 SparseSMat<T> r(lhs.size(),rhs.size());
 for(SPMatindex i=0; i<lhs.size(); ++i)
     if(lhs[i]!=(T)0)
@@ -593,6 +635,7 @@ for(SPMatindex i=0; i<lhs.size(); ++i)
 	    if(rhs[j]!=(T)0)
 		{
 		T x=lhs[i]*(conjugate?LA_traits<T>::conjugate(rhs[j]):rhs[j])*scale;
+		if(i==j) x+=LA_traits<T>::conjugate(x); //diagonal of M + M^\dagger, add() mirrors only off-diagonal elements
 		if(std::abs(x)>SPARSEEPSILON) r.add(i,j,x);
 		}
 	}

t.cc

@@ -1696,11 +1696,9 @@ cout <<"SparseSMat exp time "<<clock()/((double) (CLOCKS_PER_SEC))-t <<"\n";
 cout.flush();
 if(h.nrows()<=1024)
 {
-NRSMat<complex<double> > h3(h);
+NRMat<complex<double> > h2(h); //not via NRSMat: complex NRSMat is hermitian, while i*h1 is complex symmetric
-NRMat<complex<double> > h2(h3);
 NRMat<complex<double> >r2 = exp(h2);
-cout <<"error = "<<(r2-NRMat<complex<double> >(NRSMat<complex<double> >(r))).norm()<<endl;
+cout <<"error = "<<(r2-NRMat<complex<double> >(r)).norm()<<endl;
-cout <<"errorx = "<<(r2-NRSMat<complex<double> >(r)).norm()<<endl;
 }
 }