diff --git a/la_traits.h b/la_traits.h index 207d6f1..fbbb756 100644 --- a/la_traits.h +++ b/la_traits.h @@ -468,6 +468,8 @@ template<> struct LA_precision static double eigensolver_eps() {return 1e-6;} //default convergence threshold of davidson and lanczos static double newton_threshold() {return 1e-14;} //convergence of Polynomial::newton static double det_tiny() {return 1e-300;} //determinant underflow guard in linear_solve + static double sparse_epsilon() {return 1e-15;} //elements treated as zero in the sparse matrices (SPARSEEPSILON) + static double cholesky_epsilon() {return 1e-16;} //pivot threshold of the sparse Cholesky (CHOLESKYEPSILON) }; template<> struct LA_precision { @@ -478,6 +480,8 @@ template<> struct LA_precision static float eigensolver_eps() {return 1e-4f;} static float newton_threshold() {return 1e-6f;} static float det_tiny() {return 1e-37f;} + static float sparse_epsilon() {return 1e-7f;} + static float cholesky_epsilon() {return 1e-8f;} }; }//namespace diff --git a/nonclass.cc b/nonclass.cc index 7e796a1..d747fc6 100644 --- a/nonclass.cc +++ b/nonclass.cc @@ -2691,6 +2691,18 @@ if(a.issymmetric()) result.setsymmetric(); return result; } +template<> +const SparseMat realpart > >(const SparseMat< std::complex > &a) +{ +SparseMat result(a.nrows(), a.ncols()); +for(SparseMat< std::complex >::iterator p=a.begin(); p!=a.end(); ++p) + { + result.add(p->row,p->col,p->elem.real()); + } +if(a.issymmetric()) result.setsymmetric(); +return result; +} + template<> const SparseMat imagpart > >(const SparseMat< std::complex > &a) { @@ -2703,6 +2715,18 @@ if(a.issymmetric()) result.setsymmetric(); return result; } +template<> +const SparseMat imagpart > >(const SparseMat< std::complex > &a) +{ +SparseMat result(a.nrows(), a.ncols()); +for(SparseMat< std::complex >::iterator p=a.begin(); p!=a.end(); ++p) + { + result.add(p->row,p->col,p->elem.imag()); + } +if(a.issymmetric()) result.setsymmetric(); +return result; +} + template<> const NRMat realpart > >(const NRMat< std::complex > &a) { @@ -3044,6 +3068,27 @@ const SparseSMat< std::complex > complexmatrix >(cons return result; } +template<> +const SparseSMat< std::complex > complexmatrix >(const SparseSMat &re, const SparseSMat &im) { + if(re.nrows()!=im.nrows() || re.ncols() != im.ncols()) laerror("incompatible sizes of real and imaginary parts"); + SparseSMat< std::complex > result(re.nrows(),re.ncols()); + std::complex tmp; + + SparseSMat::iterator pre(re); + for(; pre.notend(); ++pre) { + tmp = pre->elem; + result.add(pre->row,pre->col,tmp,false); + } + + SparseSMat::iterator pim(im); + for(; pim.notend(); ++pim) { + tmp = std::complex(0,1)*(pim->elem); + result.add(pim->row,pim->col,tmp,false); + } + + return result; +} + template<> const SparseSMat< std::complex > realmatrix >(const SparseSMat &re) { SparseSMat< std::complex > result(re.nrows(),re.ncols()); @@ -3058,6 +3103,20 @@ const SparseSMat< std::complex > realmatrix >(const S return result; } +template<> +const SparseSMat< std::complex > realmatrix >(const SparseSMat &re) { + SparseSMat< std::complex > result(re.nrows(),re.ncols()); + std::complex tmp; + + SparseSMat::iterator pre(re); + for(; pre.notend(); ++pre) { + tmp = pre->elem; + result.add(pre->row,pre->col,tmp,false); + } + + return result; +} + template<> const SparseSMat< std::complex > imagmatrix >(const SparseSMat &im) { SparseSMat< std::complex > result(im.nrows(),im.ncols()); @@ -3073,6 +3132,21 @@ const SparseSMat< std::complex > imagmatrix >(const S return result; } +template<> +const SparseSMat< std::complex > imagmatrix >(const SparseSMat &im) { + SparseSMat< std::complex > result(im.nrows(),im.ncols()); + std::complex tmp; + + + SparseSMat::iterator pim(im); + for(; pim.notend(); ++pim) { + tmp = std::complex(0,1)*(pim->elem); + result.add(pim->row,pim->col,tmp,false); + } + + return result; +} + diff --git a/sparsemat.cc b/sparsemat.cc index 90162ed..3dca790 100644 --- a/sparsemat.cc +++ b/sparsemat.cc @@ -334,7 +334,7 @@ return nonzero; //number of (in principle) nonzero elements template -void SparseMat::simplify(const double sparseepsilon) +void SparseMat::simplify(const typename LA_traits::normtype sparseepsilon) { unsigned int n; if(!list) return; @@ -431,7 +431,7 @@ void SparseMat::addsafe(const SPMatindex n, const SPMatindex m, const T elem) #ifdef debug if(n<0||n>=nn||m<0||m>=mm) laerror("SparseMat out of range"); #endif -if(std::abs(elem)<=SPARSEEPSILON) return; +if(std::abs(elem)<=LA_precision::normtype>::sparse_epsilon()) return; if(!count) {count=new int; *count=1; list=NULL;} //blank new matrix else copyonwrite(); //makes also unsort add(n,m,elem); @@ -482,7 +482,7 @@ register matel *l=rhs.list; while(l) { if(rhs.symmetric || lower && l->row <=l->col || !lower && l->row >=l->col) - if(std::abs(l->elem)>SPARSEEPSILON) + if(std::abs(l->elem)>LA_precision::normtype>::sparse_epsilon()) add( l->row,l->col,sign=='+'?l->elem:- l->elem); l=l->next; } @@ -501,14 +501,14 @@ register matel *l=rhs.list; if(symmetrize) while(l) { - if(std::abs(l->elem)>SPARSEEPSILON) + if(std::abs(l->elem)>LA_precision::normtype>::sparse_epsilon()) {add( l->row,l->col,l->elem); if( l->row!=l->col) add( l->col,l->row,l->elem);} l=l->next; } else while(l) { - if(std::abs(l->elem)>SPARSEEPSILON) + if(std::abs(l->elem)>LA_precision::normtype>::sparse_epsilon()) add( l->row,l->col,l->elem); l=l->next; } @@ -528,14 +528,14 @@ register matel *l=rhs.list; if(symmetrize) while(l) { - if(std::abs(l->elem)>SPARSEEPSILON) + if(std::abs(l->elem)>LA_precision::normtype>::sparse_epsilon()) {add( l->row,l->col,- l->elem); if( l->row!=l->col) add( l->col,l->row,- l->elem);} l=l->next; } else while(l) { - if(std::abs(l->elem)>SPARSEEPSILON) + if(std::abs(l->elem)>LA_precision::normtype>::sparse_epsilon()) add( l->row,l->col,- l->elem); l=l->next; } @@ -559,7 +559,7 @@ SPMatindex i,j; for(i=0;iSPARSEEPSILON) + if( std::abs(t)>LA_precision::normtype>::sparse_epsilon()) add(i,j,t); } } @@ -598,12 +598,12 @@ if(nn==mm) //square else //diagonal of A^TA, assuming it has been simplified (only one entry per position), will be used as preconditioner only anyway while(l) { - (*rr)[l->col] += l->elem*l->elem *(l->col!=l->row && symmetric?2.:1.); + (*rr)[l->col] += l->elem*l->elem *(T)(l->col!=l->row && symmetric?2:1); l= l->next; } if(divide) { - for(unsigned int i=0; i & SparseMat::operator=(const T &a) if(!count ||nn<=0||mm<=0) laerror("assignment of scalar to undefined sparse matrix"); if(nn!=mm) laerror("assignment of scalar to non-square sparse matrix"); resize(nn,mm);//clear -if(std::abs(a)<=SPARSEEPSILON) return *this; +if(std::abs(a)<=LA_precision::normtype>::sparse_epsilon()) return *this; SPMatindex i; for(i=0;iSPARSEEPSILON + std::abs(t=rhs(i,j))>LA_precision::normtype>::sparse_epsilon() ) add(i,j,t); } } @@ -774,7 +774,7 @@ matel *l=list; while(l) //include the mirror picture of elements into the list { if( - std::abs(l->elem)>SPARSEEPSILON && + std::abs(l->elem)>LA_precision::normtype>::sparse_epsilon() && l->row!=l->col) add(l->col,l->row,l->elem); //not OK for complex-hermitean l=l->next; } @@ -812,12 +812,32 @@ while(l) return r; } +template<> +const float SparseMat::dot(const NRMat &rhs) const +{ +float r=0; +matel *l=list; +while(l) + { + r+= l->elem*rhs[l->row][l->col]; + if(symmetric&&l->row!=l->col) r+=l->elem*rhs[l->col][l->row]; + l=l->next; + } +return r; +} + template<> void NRMat >::gemm(const std::complex &beta, const SparseMat > &a, const char transa, const NRMat > &b, const char transb, const std::complex &alpha) { laerror("not implemented yet"); } +template<> +void NRMat >::gemm(const std::complex &beta, const SparseMat > &a, const char transa, const NRMat > &b, const char transb, const std::complex &alpha) +{ +laerror("not implemented yet"); +} + template<> void NRMat::gemm(const double &beta, const SparseMat &a, const char transa, const NRMat &b, const char transb, const double &alpha) @@ -867,6 +887,54 @@ else } +template<> +void NRMat::gemm(const float &beta, const SparseMat &a, const char transa, const NRMat &b, const char transb, const float &alpha) +{ +bool transpa = tolower(transa)!='n'; //not OK for complex +bool transpb = tolower(transb)!='n'; //not OK for complex +if(nn!=(int)(transpa?a.ncols():a.nrows()) | mm!= (transpb?b.nrows():b.ncols()) || (int)(transpa?a.nrows():a.ncols()) != (transpb?b.ncols():b.nrows())) laerror("incompatible sizes in gemm"); +copyonwrite(); +if(beta!=(float)0) (*this) *= beta; +else memset(v,0,nn*mm*sizeof(float)); + +matel *l=a.getlist(); +if(alpha==(float)0 || !l) return; + +//raw pointers to the full matrices +const float *bp= b; +float *p= *this; + +int ldb=(transpb?b.ncols():1); +int bstep=(transpb?1:b.ncols()); +int len=(transpb?b.nrows():b.ncols()); + +if(a.issymmetric()) + { + while(l) + { + cblas_saxpy(len, alpha*l->elem , bp+l->row*bstep, ldb, p+l->col*mm, 1); + if(l->row!=l->col) cblas_saxpy(len, alpha*l->elem , bp+l->col*bstep, ldb, p+l->row*mm, 1); + l=l->next; + } + } +else + { + if(transpa) + while(l) + { + cblas_saxpy(len, alpha*l->elem , bp+l->row*bstep, ldb, p+l->col*mm, 1); + l=l->next; + } + else + while(l) + { + cblas_saxpy(len, alpha*l->elem , bp+l->col*bstep, ldb, p+l->row*mm, 1); + l=l->next; + } + } + +} + template @@ -1036,7 +1104,7 @@ if(transp) while(l) { register T t=alpha*l->elem; - if(std::abs(t)>SPARSEEPSILON) + if(std::abs(t)>LA_precision::normtype>::sparse_epsilon()) add( l->col,l->row,t); l=l->next; } @@ -1044,7 +1112,7 @@ else while(l) { register T t=alpha*l->elem; - if(std::abs(t)>SPARSEEPSILON) + if(std::abs(t)>LA_precision::normtype>::sparse_epsilon()) add( l->row,l->col,t); l=l->next; } @@ -1329,7 +1397,7 @@ incsize(rhs.nrows(), rhs.ncols()); T tmp; for(SPMatindex i=0; i<(SPMatindex)rhs.nrows(); ++i) for(SPMatindex j=0; j<(SPMatindex)rhs.ncols(); ++j) - if(std::abs(tmp=rhs(i,j))>SPARSEEPSILON) add(n0+i,m0+j,tmp); + if(std::abs(tmp=rhs(i,j))>LA_precision::normtype>::sparse_epsilon()) add(n0+i,m0+j,tmp); return *this; } @@ -1345,7 +1413,7 @@ T tmp; incsize(rhs.nrows(), rhs.ncols()); for(SPMatindex i=0; i<(SPMatindex)rhs.nrows(); ++i) for(SPMatindex j=0; j<=i; ++j) //symmetric storage: one triangle only - if(std::abs(tmp=rhs(i,j))>SPARSEEPSILON) add(n0+i,m0+j,tmp); + if(std::abs(tmp=rhs(i,j))>LA_precision::normtype>::sparse_epsilon()) add(n0+i,m0+j,tmp); return *this; } @@ -1364,7 +1432,7 @@ incsize(rhs.nrows(), rhs.ncols()); register matel *l=rhs.list; while(l) { - if(std::abs(l->elem)>SPARSEEPSILON) + if(std::abs(l->elem)>LA_precision::normtype>::sparse_epsilon()) add(n0+l->row,m0+l->col,l->elem); l=l->next; } @@ -1429,8 +1497,12 @@ INSTANTIZE(std::complex) //some functions are not OK for hermitean matri //// forced instantization in the corresponding object file template class SparseMat; template class SparseMat >; +template class SparseMat; +template class SparseMat >; template class SemiSparseMat; template class SemiSparseMat >; +template class SemiSparseMat; +template class SemiSparseMat >; #define INSTANTIZE(T) \ @@ -1443,10 +1515,8 @@ template void NRVec::gemv(const T beta, const SparseMat &a, const char tra INSTANTIZE(double) INSTANTIZE(std::complex) - -//needed by the instantiation of NRVec and NRVec > -template void NRVec::gemv(const float beta, const SparseMat &a, const char trans, const float alpha, const NRVec &x, const bool treat_as_symmetric); -template void NRVec >::gemv(const std::complex beta, const SparseMat > &a, const char trans, const std::complex alpha, const NRVec > &x, const bool treat_as_symmetric); +INSTANTIZE(float) +INSTANTIZE(std::complex) }//namespace diff --git a/sparsemat.h b/sparsemat.h index 452e38c..bbf3b0c 100644 --- a/sparsemat.h +++ b/sparsemat.h @@ -141,7 +141,7 @@ public: unsigned int length() const; void copyonwrite(bool detachonly=false, bool deep=true); void clear() {bool sym=symmetric; resize(nn,mm); symmetric=sym;} //zero matrix of the same dimensions and symmetry - void simplify(const double sparseepsilon=SPARSEEPSILON); + void simplify(const typename LA_traits::normtype sparseepsilon=LA_precision::normtype>::sparse_epsilon()); const T trace() const; const typename LA_traits::normtype norm(const T scalar=(T)0) const; //is const only mathematically, not in internal implementation - we have to simplify first unsigned int sort(int type) const; diff --git a/sparsesmat.cc b/sparsesmat.cc index 7663c96..738c4ce 100644 --- a/sparsesmat.cc +++ b/sparsesmat.cc @@ -239,7 +239,7 @@ for(SPMatindex i=0; i) //// forced instantization of functions in the header in the corresponding object file template class SparseSMat; template class SparseSMat >; +template class SparseSMat; +template class SparseSMat >; template SparseSMat NRMat::operator*(const SparseSMat &rhs) const; +template SparseSMat NRMat::operator*(const SparseSMat &rhs) const; 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, bool orcaformat) const; +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); +template void SparseSMat >::get(int fd, bool dimen, bool transp, bool orcaformat); */ diff --git a/sparsesmat.h b/sparsesmat.h index be1ee7d..a0c12c7 100644 --- a/sparsesmat.h +++ b/sparsesmat.h @@ -200,7 +200,7 @@ count(new int(1)) v= new std::map * [nn]; memset(v,0,nn*sizeof(std::map *)); int i,j; -for(i=0; iSPARSEEPSILON) (*this).add(i,j,rhs(i,j),true); +for(i=0; iLA_precision::normtype>::sparse_epsilon()) (*this).add(i,j,rhs(i,j),true); } //from SparseMat: duplicates summed; a symmetric SparseMat (one triangle stored, complex symmetric) is unfolded without conjugation @@ -227,7 +227,7 @@ count(new int(1)) v= new std::map * [nn]; memset(v,0,nn*sizeof(std::map *)); int i,j; -for(i=0; iSPARSEEPSILON) (*this).add(i,j,rhs(i,j),false); +for(i=0; iLA_precision::normtype>::sparse_epsilon()) (*this).add(i,j,rhs(i,j),false); } @@ -437,7 +437,7 @@ for(SPMatindex i=0; i l; typename std::map::iterator p; for(p=v[i]->begin(); p!=v[i]->end(); ++p) - if(std::abs(p->second) < SPARSEEPSILON) l.push_front(p->first); else ++count; + if(std::abs(p->second) < LA_precision::normtype>::sparse_epsilon()) l.push_front(p->first); else ++count; typename std::list::iterator q; for(q=l.begin(); q!=l.end(); ++q) v[i]->erase(*q); if(v[i]->size() == 0) {delete v[i]; v[i]=NULL;} @@ -548,7 +548,7 @@ for(SPMatindex i=0; i r0(r);r.copyonwrite(); -//perform complex, positive semidefinite Cholesky with thresholding by SPARSEEPSILON +//perform complex, positive semidefinite Cholesky with thresholding by LA_precision<>::sparse_epsilon() SPMatindex i,j,k; int rank=0; for(k=0; k::iterator p; p= r.v[k]->find(k); if(p==r.v[k]->end()) continue; //must not break due to the a priori pivoting - if(LA_traits::realpart(p->second) < CHOLESKYEPSILON) continue; //must not break due to the a priori pivoting + if(LA_traits::realpart(p->second) < LA_precision::normtype>::cholesky_epsilon()) continue; //must not break due to the a priori pivoting ++rank; typename LA_traits::normtype tmp = std::sqrt(LA_traits::realpart(p->second)); p->second = tmp; @@ -572,7 +572,7 @@ for(k=0; k::conjugate(linek[j]); NRVec linek_used(0,nn); - if(std::abs(akj)>SPARSEEPSILON) + if(std::abs(akj)>LA_precision::normtype>::sparse_epsilon()) { for(p=r.v[j]->begin(); p!=r.v[j]->end(); ++p) { @@ -581,7 +581,7 @@ for(k=0; k SPARSEEPSILON) + if(!linek_used[i] && std::abs(linek[i]) > LA_precision::normtype>::sparse_epsilon()) { (*r.v[j])[i] = -akj * linek[i]; } @@ -611,7 +611,7 @@ for(k=0; k *vnew = new typename std::map; for(p=r.v[k]->begin(); p!=r.v[k]->end(); ++p) { - if(std::abs(p->second) > SPARSEEPSILON) (*vnew)[pivot[p->first]] = p->second; + if(std::abs(p->second) > LA_precision::normtype>::sparse_epsilon()) (*vnew)[pivot[p->first]] = p->second; } delete r.v[k]; r.v[k]=vnew; @@ -636,7 +636,7 @@ 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); + if(std::abs(x)>LA_precision::normtype>::sparse_epsilon()) r.add(i,j,x); } } return r;