From 79a0cf9247af88254b5396691ddd270f559f314c Mon Sep 17 00:00:00 2001 From: Jiri Pittner Date: Tue, 27 Oct 2026 23:18:40 +0100 Subject: [PATCH] bugfixes in sparsemat and sparsesmat --- sparsemat.cc | 10 ++++--- sparsemat.h | 8 +++--- sparsesmat.cc | 74 +++++++++++++++++++++++-------------------------- sparsesmat.h | 77 +++++++++++++++++++++++++++++++++++++++------------ t.cc | 6 ++-- 5 files changed, 107 insertions(+), 68 deletions(-) diff --git a/sparsemat.cc b/sparsemat.cc index 678cee9..342eb05 100644 --- a/sparsemat.cc +++ b/sparsemat.cc @@ -44,8 +44,9 @@ static inline void SWAP(T &a, T &b) template -void SparseMat::get(int fd, bool dimen, bool transp) +void SparseMat::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 -void SparseMat::put(int fd,bool dimen, bool transp) const +void SparseMat::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 & SparseMat::oplusequal(const SparseMat &rhs);\ template SparseMat & SparseMat::oplusequal(const NRMat &rhs);\ template SparseMat & SparseMat::oplusequal(const NRSMat &rhs);\ -template void SparseMat::get(int fd, bool dimen, bool transp); \ -template void SparseMat::put(int fd, bool dimen, bool transp) const; \ +template void SparseMat::get(int fd, bool dimen, bool transp, bool orcaformat); \ +template void SparseMat::put(int fd, bool dimen, bool transp, bool orcaformat) const; \ template void SparseMat::copyonwrite(); \ template void SparseMat::unsort(); \ template void SparseMat::resize(const SPMatindex n, const SPMatindex m); \ diff --git a/sparsemat.h b/sparsemat.h index aaab3e0..452e38c 100644 --- a/sparsemat.h +++ b/sparsemat.h @@ -118,8 +118,8 @@ public: void setlist(matel *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 put(int fd, bool dimensions=1, bool transposed=false) const; + 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, 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 diagonal; SparseMat 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 diff --git a/sparsesmat.cc b/sparsesmat.cc index cac53ba..7663c96 100644 --- a/sparsesmat.cc +++ b/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 SparseSMat NRMat::operator*(const SparseSMat &rhs) const { +if(mm!=(int)rhs.nrows()) laerror("incompatible sizes in NRMat*SparseSMat"); SparseSMat r(nn,rhs.ncols()); -if(mm!=rhs.nrows()) laerror("incompatible sizes in NRMat*SparseSMat"); -for(SPMatindex k=0; k * kl = rhs.line(k); - if(kl) +for(int i=0; i *ri=NULL; + for(int k=0; k::iterator p; - int i,j; - NRVec kline(kl->size()); - NRVec klineind(kl->size()); - for(p=kl->begin(), i=0; p!=kl->end(); ++p,++i) - { - klineind[i] = p->first; - kline[i] = p->second; - } - NRVec kcol = column(k); - - //multiply - NRMat prod=kcol.otimes(kline,false,1.); - - //scatter the results - for(i=0; i *kl = rhs.line(k); + if(aik==(T)0 || !kl) continue; + if(!ri) ri = r[i] = new std::map; + typename std::map::const_iterator p; + for(p=kl->begin(); p!=kl->end(); ++p) (*ri)[p->first] += aik * p->second; //new map element is value-initialized to 0 } } 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; kempty() && !b.v[k]->empty()) //nonempty in both (allocated rows can be empty, e.g. after cholesky) { NRVec av(a.v[k]->size()); NRVec bv(b.v[k]->size()); @@ -148,7 +136,7 @@ void SparseSMat::gemv(const T beta, NRVec &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::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::iterator p(rhs); for(; p.notend(); ++p) add(p->row+fromrow, p->col+fromcol, p->elem, false); } template -void SparseSMat::get(int fd, bool dimen, bool transp) { +void SparseSMat::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::IOtype tmp; LA_traits::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 -void SparseSMat::put(int fd, bool dimen, bool transp) const { +void SparseSMat::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,&mm,sizeof(SPMatindex))) laerror("cannot write() 2 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,&(transp ? nn : mm),sizeof(SPMatindex))) laerror("cannot write() 2 in SparseSMat::put"); } typename SparseSMat::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,&(p->col),sizeof(SPMatindex))) laerror("cannot write() 4 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,&(transp ? p->row : p->col),sizeof(SPMatindex))) laerror("cannot write() 4 in SparseSMat::put"); typename LA_traits_io::IOtype tmp = p->elem; LA_traits::put(fd,tmp,dimen,transp); // general way to work when elem is some non-scalar class again } @@ -347,8 +341,8 @@ template SparseSMat & SparseSMat::operator+=(const T &a); \ template SparseSMat & SparseSMat::operator-=(const T &a); \ template LA_traits::normtype SparseSMat::norm(const T scalar) const; \ template const T* SparseSMat::diagonalof(NRVec &r, const bool divide, bool cache) const; \ -template void SparseSMat::get(int fd, bool dimen, bool transp); \ -template void SparseSMat::put(int fd, bool dimen, bool transp) const; \ +template void SparseSMat::get(int fd, bool dimen, bool transp, bool orcaformat); \ +template void SparseSMat::put(int fd, bool dimen, bool transp, bool orcaformat) const; \ INSTANTIZE(double) @@ -358,10 +352,12 @@ INSTANTIZE(complex) //// forced instantization of functions in the header in the corresponding object file template class SparseSMat; template class SparseSMat >; +template SparseSMat NRMat::operator*(const SparseSMat &rhs) const; +template SparseSMat > NRMat >::operator*(const SparseSMat > &rhs) const; /*activate this if needed -template void SparseSMat >::put(int fd, bool dimen, bool transp) const; -template void SparseSMat >::get(int fd, bool dimen, bool transp); +template void SparseSMat >::put(int fd, bool dimen, bool transp, bool orcaformat) const; +template void SparseSMat >::get(int fd, bool dimen, bool transp, bool orcaformat); */ diff --git a/sparsesmat.h b/sparsesmat.h index 7ac6bf5..be1ee7d 100644 --- a/sparsesmat.h +++ b/sparsesmat.h @@ -40,9 +40,15 @@ namespace LA { -//symmetric sparse matrix class with a representation for efficient exponentiatiation -//in particular we need a unitary symmetric complex matrix as exp(iH) with H real symmetric +//sparse matrix stored by rows, each row a std::map from column index to value, for efficient exponentiation and products +//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 @@ -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 *line(SPMatindex n) const {return v[n];}; + inline const std::map *line(const SPMatindex n) const {return v[n];}; //read-only row access typedef std::map *ROWTYPE; inline typename SparseSMat::ROWTYPE & operator[](const SPMatindex i) {return v[i];}; void clear() {resize(nn,mm);} @@ -90,15 +96,15 @@ public: void gemv(const T beta, NRVec &r, const char trans, const T alpha, const NRVec &x) const; inline const NRVec operator*(const NRVec &rhs) const {NRVec result(nn); this->gemv((T)0,result,'n',(T)1,rhs); return result;}; typename LA_traits::normtype norm(const T scalar=(T)0) const; - inline const SparseSMat operator*(const SparseSMat &rhs) const {SparseSMat 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 - 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 + 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*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 put(int fd, bool dimen, bool transp) const; + 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=1, bool transp=0, bool orcaformat=false) const; int nrows() const {return nn;} int ncols() const {return mm;} SparseSMat cholesky(void) const; @@ -197,12 +203,27 @@ int i,j; for(i=0; iSPARSEEPSILON) (*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 +SparseSMat::SparseSMat(const SparseMat &rhs) +:nn(rhs.nrows()), mm(rhs.ncols()), +count(new int(1)) +{ +v= new std::map * [nn]; +memset(v,0,nn*sizeof(std::map *)); +const bool sym=rhs.issymmetric(); +for(const matel *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 SparseSMat::SparseSMat(const NRMat &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 * [nn]; memset(v,0,nn*sizeof(std::map *)); int i,j; @@ -234,7 +255,7 @@ count = new int(1); v=new T[nn2]; memset(v,0,nn2*sizeof(T)); typename SparseSMat::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,14 +328,13 @@ if(*count > 1) //it was shared } else //it was not shared { - mm=m; //delete all trees for(SPMatindex i=0; inn) {delete[] v; v= new std::map * [n];} //shrinking keeps the longer array + nn=n; + mm=m; + for(SPMatindex i=0; i::copyonwrite() } +//general product A B row by row: row i of the result is sum_k A_ik (row k of B) +template +const SparseSMat SparseSMat::operator*(const SparseSMat &rhs) const +{ +if(mm!=rhs.nn) laerror("incompatible sizes in SparseSMat::operator*"); +SparseSMat r(nn,rhs.mm); +for(SPMatindex i=0; i *ri=NULL; + typename std::map::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; + 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 void SparseSMat::add(const SPMatindex n, const SPMatindex m, const T elem, const bool both) { @@ -375,11 +416,11 @@ typename std::map::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; 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::conjugate(elem); else (*v[m])[n] = LA_traits::conjugate(elem); } //addition can lead to zero, but do not treat it now, make a simplify } @@ -585,6 +626,7 @@ return r; template SparseSMat otimes_sparse(const NRVec &lhs, const NRVec &rhs, const bool conjugate=false, const T &scale=1) { +if(lhs.size()!=rhs.size()) laerror("otimes_sparse requires vectors of equal length"); SparseSMat r(lhs.size(),rhs.size()); for(SPMatindex i=0; i::conjugate(rhs[j]):rhs[j])*scale; + if(i==j) x+=LA_traits::conjugate(x); //diagonal of M + M^\dagger, add() mirrors only off-diagonal elements if(std::abs(x)>SPARSEEPSILON) r.add(i,j,x); } } diff --git a/t.cc b/t.cc index fb30c6d..eb93a4b 100644 --- a/t.cc +++ b/t.cc @@ -1696,11 +1696,9 @@ cout <<"SparseSMat exp time "< > h3(h); -NRMat > h2(h3); +NRMat > h2(h); //not via NRSMat: complex NRSMat is hermitian, while i*h1 is complex symmetric NRMat >r2 = exp(h2); -cout <<"error = "<<(r2-NRMat >(NRSMat >(r))).norm()<