From 5f5f0343a688d08f738e7121ed270273883a49d4 Mon Sep 17 00:00:00 2001 From: jiri Date: Thu, 7 Jan 2010 16:10:12 +0000 Subject: [PATCH] *** empty log message *** --- la_traits.h | 22 ++++--- mat.h | 3 + nonclass.cc | 44 +++++++++++++ nonclass.h | 4 ++ sparsesmat.cc | 164 ++++++++++++++++++++++++++++++++++++++++++----- sparsesmat.h | 173 ++++++++++++++++++++++++++++++++++++++++++++++++++ t.cc | 134 +++++++++++++++++++++++++++++++++++++- 7 files changed, 517 insertions(+), 27 deletions(-) diff --git a/la_traits.h b/la_traits.h index 7911498..a45817e 100644 --- a/la_traits.h +++ b/la_traits.h @@ -18,7 +18,7 @@ // //for autotools // -#include "config.h" +//#include "config.h" //this would force the user of the library to have config.h //////////////////////////////////////////////////////////////////////////// //LA traits classes and generally needed includes @@ -219,6 +219,8 @@ static void copy(complex *dest, complex *src, unsigned int n) {memcpy(dest static void clear(complex *dest, unsigned int n) {memset(dest,0,n*sizeof(complex));} static void copyonwrite(complex &x) {}; static void clearme(complex &x) {x=0;}; +static inline complex conjugate(complex &x) {return complex(x.real(),-x.imag());}; +static inline C realpart(complex &x) {return x.real();} }; //non-complex scalars @@ -241,6 +243,8 @@ static void copy(C *dest, C *src, unsigned int n) {memcpy(dest,src,n*sizeof(C)); static void clear(C *dest, unsigned int n) {memset(dest,0,n*sizeof(C));} static void copyonwrite(C &x) {}; static void clearme(complex &x) {x=0;}; +static inline C conjugate(C&x) {return x;}; +static inline C realpart(C &x) {return x;} }; @@ -260,10 +264,10 @@ static inline bool bigger(const C &x, const C &y) {return x>y;} \ static inline bool smaller(const C &x, const C &y) {return x &x) {return x.norm();} \ static inline void axpy (X&s, const X &x, const C c) {s.axpy(c,x);} \ -static void put(int fd, const C &x, bool dimensions=1, bool transp=0) {x.put(fd,dimensions,transp);} \ -static void get(int fd, C &x, bool dimensions=1, bool transp=0) {x.get(fd,dimensions,transp);} \ -static void multiput(unsigned int n,int fd, const C *x, bool dimensions=1) {for(unsigned int i=0; i &x, bool dimensions=1, bool transp=0) {x.put(fd,dimensions,transp);} \ +static void get(int fd, X &x, bool dimensions=1, bool transp=0) {x.get(fd,dimensions,transp);} \ +static void multiput(unsigned int n,int fd, const X *x, bool dimensions=1) {for(unsigned int i=0; i *x, bool dimensions=1) {for(unsigned int i=0; i &x) {x.copyonwrite();}\ @@ -292,10 +296,10 @@ static inline bool bigger(const C &x, const C &y) {return x>y;} \ static inline bool smaller(const C &x, const C &y) {return x &x) {return x.norm();} \ static inline void axpy (X&s, const X &x, const C c) {s.axpy(c,x);} \ -static void put(int fd, const C &x, bool dimensions=1, bool transp=0) {x.put(fd,dimensions);} \ -static void get(int fd, C &x, bool dimensions=1, bool transp=0) {x.get(fd,dimensions);} \ -static void multiput(unsigned int n,int fd, const C *x, bool dimensions=1) {for(unsigned int i=0; i &x, bool dimensions=1, bool transp=0) {x.put(fd,dimensions);} \ +static void get(int fd, X &x, bool dimensions=1, bool transp=0) {x.get(fd,dimensions);} \ +static void multiput(unsigned int n,int fd, const X *x, bool dimensions=1) {for(unsigned int i=0; i *x, bool dimensions=1) {for(unsigned int i=0; i &x) {x.copyonwrite();} \ diff --git a/mat.h b/mat.h index 9ed2abd..c6f6822 100644 --- a/mat.h +++ b/mat.h @@ -128,7 +128,9 @@ public: const T trace() const; //members concerning sparse matrix + SparseSMat & operator*(const SparseSMat &rhs) const; explicit NRMat(const SparseMat &rhs); // dense from sparse + explicit NRMat(const SparseSMat &rhs); // dense from sparse NRMat & operator+=(const SparseMat &rhs); NRMat & operator-=(const SparseMat &rhs); void gemm(const T &beta, const SparseMat &a, const char transa, const NRMat &b, const char transb, const T &alpha);//this = alpha*op( A )*op( B ) + beta*this @@ -146,6 +148,7 @@ public: #include "vec.h" #include "smat.h" #include "sparsemat.h" +#include "sparsesmat.h" namespace LA { // ctors diff --git a/nonclass.cc b/nonclass.cc index 0dd38e4..6c3966d 100644 --- a/nonclass.cc +++ b/nonclass.cc @@ -949,6 +949,50 @@ return r; } +//Cholesky interface +extern "C" void FORNAME(dpotrf)(const char *UPLO, const int *N, double *A, const int *LDA, int *INFO); +extern "C" void FORNAME(zpotrf)(const char *UPLO, const int *N, complex *A, const int *LDA, int *INFO); + +void cholesky(NRMat &a, bool upper) +{ +if(a.nrows()!=a.ncols()) laerror("matrix must be square in Cholesky"); +int lda=a.ncols(); +int n=a.nrows(); +char uplo=upper?'U':'L'; +int info; +a.copyonwrite(); +FORNAME(dpotrf)(&uplo, &n, a, &lda, &info); +if(info) {std::cerr << "Lapack error "< > &a, bool upper) +{ +if(a.nrows()!=a.ncols()) laerror("matrix must be square in Cholesky"); +int lda=a.ncols(); +int n=a.nrows(); +char uplo=upper?'U':'L'; +int info; +a.copyonwrite(); +a.transposeme();//switch to Fortran order +FORNAME(zpotrf)(&uplo, &n, a, &lda, &info); +if(info) {std::cerr << "Lapack error "< &a, NRVec &w, NRMat b, int n) diff --git a/nonclass.h b/nonclass.h index 9f7a3ad..740c933 100644 --- a/nonclass.h +++ b/nonclass.h @@ -135,6 +135,10 @@ extern const NRMat< complex > realmatrix (const NRMat&); extern const NRMat< complex > imagmatrix (const NRMat&); extern const NRMat< complex > complexmatrix (const NRMat&, const NRMat&); +//Cholesky decomposition +extern void cholesky(NRMat &a, bool upper=1); +extern void cholesky(NRMat > &a, bool upper=1); + //inverse by means of linear solve, preserving rhs intact template const NRMat inverse(NRMat a, T *det=0) diff --git a/sparsesmat.cc b/sparsesmat.cc index 0ecaa39..bd9f95c 100644 --- a/sparsesmat.cc +++ b/sparsesmat.cc @@ -27,6 +27,45 @@ namespace LA { + +//dense times sparse (not necessarily symmetric) +template +SparseSMat & NRMat::operator*(const SparseSMat &rhs) const +{ +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) + { + //gather the data + typename std::map::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 void SparseSMat::gemm(const T beta, const SparseSMat &a, const char transa, const SparseSMat &b, const char transb, const T alpha) { @@ -46,16 +85,14 @@ for(SPMatindex k=0; k::iterator p; int i,j; - for(p=a.v[k]->begin(), i=0; p!=a.v[k]->end(); ++p,++i) - { - ai[i] = p->first; - av[i] = p->second; - } - for(p=b.v[k]->begin(), i=0; p!=b.v[k]->end(); ++p,++i) - { - bi[i] = p->first; - bv[i] = p->second; - } + if(tolower(transa)=='c') + for(p=a.v[k]->begin(), i=0; p!=a.v[k]->end(); ++p,++i) { ai[i] = p->first; av[i] = LA_traits::conjugate(p->second); } + else + for(p=a.v[k]->begin(), i=0; p!=a.v[k]->end(); ++p,++i) { ai[i] = p->first; av[i] = p->second; } + if(tolower(transb)=='c') + for(p=b.v[k]->begin(), i=0; p!=b.v[k]->end(); ++p,++i) { bi[i] = p->first; bv[i] = LA_traits::conjugate(p->second); } + else + for(p=b.v[k]->begin(), i=0; p!=b.v[k]->end(); ++p,++i) { bi[i] = p->first; bv[i] = p->second; } //make multiply via blas NRMat prod=av.otimes(bv,false,alpha); @@ -96,8 +133,13 @@ copyonwrite(); for(SPMatindex i=0; i; - typename std::map::iterator p; - for(p=x.v[i]->begin(); p!=x.v[i]->end(); ++p) (*v[i])[p->first] = p->second * alpha; + typename std::map::iterator p,q; + for(p=x.v[i]->begin(); p!=x.v[i]->end(); ++p) + { + q=v[i]->find(p->first); + if(q!=v[i]->end()) q->second += p->second * alpha; + else (*v[i])[p->first] = p->second * alpha; + } } simplify(); } @@ -165,26 +207,105 @@ for(SPMatindex i=0; i typename LA_traits::normtype SparseSMat::norm(const T scalar) const { typename LA_traits::normtype sum=0; for(SPMatindex i=0; i::iterator p; - p= v[i]->find(i); - if(p!=v[i]->end()) sum += LA_traits::sqrabs(p->second - scalar); - else sum += LA_traits::sqrabs(scalar); + bool diagonal_present=false; + for(p=v[i]->begin(); p!=v[i]->end(); ++p) //loop over all existing elements + { + if(i==p->first) {diagonal_present=true; sum += LA_traits::sqrabs(p->second - scalar);} + else sum += LA_traits::sqrabs(p->second); + } + if(!diagonal_present) sum += LA_traits::sqrabs(scalar); //there was zero on the diagonal } - else sum += LA_traits::sqrabs(scalar); //missing diagonal element + else sum += LA_traits::sqrabs(scalar); //missing whole line, subtracted diagonal element contributes return std::sqrt(sum); } +//get diagonal, do not construct a new object, but store in existing one +template +const T* SparseSMat::diagonalof(NRVec &r, const bool divide, bool cache) const +{ +if(nn!=r.size()) laerror("incompatible vector size in diagonalof()"); +NRVec *rr; + +r.copyonwrite(); +if(divide) {rr=new NRVec(nn); *rr=(T)0;} +else {r=(T)0; rr=&r;} +for(SPMatindex i=0; i::iterator p; + p= v[i]->find(i); + if(p!=v[i]->end()) (*rr)[i] += p->second; + } +if(divide) + { + for(unsigned int i=0; i +void SparseSMat::get(int fd, bool dimen, bool transp) { + errno=0; + SPMatindex dim[2]; + +if(dimen) { + if(2*sizeof(SPMatindex)!=read(fd,&dim,2*sizeof(SPMatindex))) laerror("read() error in SparseSMat::get "); + if(dim[0]!=dim[1]) laerror("SparseSMat must be square (nonsquare read in ::get)"); + resize(dim[0]); + } +else copyonwrite(); + +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); + } +while(1); +} + + + +template +void SparseSMat::put(int fd, bool dimen, bool transp) const { + errno=0; + if(dimen) { + if(sizeof(SPMatindex)!=write(fd,&nn,sizeof(SPMatindex))) laerror("cannot write() 1 in SparseSMat::put"); + if(sizeof(SPMatindex)!=write(fd,&nn,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"); + 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 + } + + SPMatindex sentinel[2]; + sentinel[0] = sentinel[1] = (SPMatindex) -1; + if(2*sizeof(SPMatindex) != write(fd,&sentinel,2*sizeof(SPMatindex))) laerror("cannot write() 5 in SparseSMat::put"); +} + + + #define INSTANTIZE(T) \ template void SparseSMat::gemm(const T beta, const SparseSMat &a, const char transa, const SparseSMat &b, const char transb, const T alpha); \ @@ -195,6 +316,9 @@ template SparseSMat & SparseSMat::operator=(const T &a); \ 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; \ INSTANTIZE(double) @@ -205,4 +329,10 @@ INSTANTIZE(complex) template class SparseSMat; template class SparseSMat >; +/*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); +*/ + + }//namespace diff --git a/sparsesmat.h b/sparsesmat.h index 58bd828..1cd19e5 100644 --- a/sparsesmat.h +++ b/sparsesmat.h @@ -31,10 +31,13 @@ #include "vec.h" #include "mat.h" #include "smat.h" +#include "qsort.h" #include #include +#define CHOLESKYEPSILON 1e-16 + namespace LA { //symmetric sparse matrix class with a representation for efficient exponentiatiation @@ -55,9 +58,11 @@ public: SparseSMat(const SparseSMat &rhs); explicit SparseSMat(const SparseMat &rhs); explicit SparseSMat(const NRSMat &rhs); + explicit SparseSMat(const NRMat &rhs); SparseSMat & operator=(const SparseSMat &rhs); void copyonwrite(); void resize(const SPMatindex n); + std::map *line(SPMatindex n) const {return v[n];}; void clear() {resize(nn);} unsigned long long simplify(); ~SparseSMat(); @@ -74,15 +79,20 @@ public: void axpy(const T alpha, const SparseSMat &x, const bool transp=0); // this+= a*x inline SparseSMat & operator+=(const SparseSMat &rhs) {axpy((T)1,rhs); return *this;}; inline SparseSMat & operator-=(const SparseSMat &rhs) {axpy((T)-1,rhs); return *this;}; + const T* diagonalof(NRVec &, const bool divide=0, bool cache=false) const; //get diagonal 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); 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 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 {}; + void get(int fd, bool dimen, bool transp); + void put(int fd, bool dimen, bool transp) const; int nrows() const {return nn;} int ncols() const {return nn;} + SparseSMat cholesky(void) const; class iterator {//not efficient, just for output to ostreams private: @@ -145,6 +155,7 @@ public: iterator end() const {return iterator(NULL);} }; + template SparseSMat::SparseSMat(const SPMatindex n) :nn(n), @@ -165,6 +176,19 @@ int i,j; for(i=0; iSPARSEEPSILON) (*this).add(i,j,rhs(i,j),true); } +template +SparseSMat::SparseSMat(const NRMat &rhs) +:nn(rhs.nrows()), +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; +for(i=0; iSPARSEEPSILON) (*this).add(i,j,rhs(i,j),false); +} + + template SparseSMat::SparseSMat(const SparseSMat &rhs) { @@ -189,6 +213,27 @@ for(; p.notend(); ++p) if(p->row <= p->col) (*this)(p->row,p->col)=p->elem; #undef nn2 +//construct dense from sparse +template +NRMat::NRMat(const SparseSMat &rhs) : +nn(rhs.nrows()), +mm(rhs.ncols()), +count(new int(1)) +{ +#ifdef MATPTR + v = new T*[nn]; + v[0] = new T[mm*nn]; + for (int i=1; i::iterator p(rhs); +for(; p.notend(); ++p) (*this)(p->row,p->col)= p->elem; +} + + + template SparseSMat::~SparseSMat() { @@ -360,5 +405,133 @@ std::istream& operator>>(std::istream &s, SparseSMat &x) } +//Cholesky decomposition, pivoted, positive semidefinite, not in place +//it is NOT checked that the input matrix is symmetric/hermitean +//result.transpose(true)*result reproduces the original matrix +//Due to pivoting the result is upper triangular only before permutation +// +template +SparseSMat SparseSMat::cholesky(void) const +{ + +//we need real values for sorting, if T is already real it makes just an unnecessary copy of one vector +NRVec::normtype> diagreal(nn); +{ +NRVec diag(nn); diagonalof(diag); +for(SPMatindex i=0; i::realpart(diag[i]); +} + +NRVec pivot(nn); +for(int i=0; i invpivot(nn); +for(int i=0; i r; +r.nn=nn; +r.count = new int(1); +r.v = new std::map * [nn]; +for(SPMatindex i=0; i; + typename std::map::iterator p; + for(p=v[pivot[i]]->begin(); p!=v[pivot[i]]->end(); ++p) + { + if(invpivot[p->first] >= i) + { + (*r.v[i])[invpivot[p->first]] = p->second; + } + } + } + else + r.v[i]= NULL; + +//std::cout <<"Permuted upper triangle matrix\n"< r0(r);r.copyonwrite(); + +//perform complex, positive semidefinite Cholesky with thresholding by SPARSEEPSILON +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 + ++rank; + typename LA_traits::normtype tmp = std::sqrt(LA_traits::realpart(p->second)); + p->second = tmp; + NRVec linek(0.,nn); + for(p=r.v[k]->begin(); p!=r.v[k]->end(); ++p) + { + if(p->first > k) p->second /= tmp; + linek[p->first] = p->second; + } + for(j=k+1; j::conjugate(linek[j]); + NRVec linek_used(0,nn); + if(std::abs(akj)>SPARSEEPSILON) + { + for(p=r.v[j]->begin(); p!=r.v[j]->end(); ++p) + { + p->second -= akj * linek[p->first]; + linek_used[p->first]=1; + } + //subtract also elements nonzero in line k but non-existent in line j + for(i=j; i SPARSEEPSILON) + { + (*r.v[j])[i] = -akj * linek[i]; + } + } + } + } + +/* +SparseSMat br(nn); +br.gemm(0,r,'c',r,'n',1); +//cancel low triangle from br +for(k=0; k::iterator p; + for(p=br.v[k]->begin(); p!=br.v[k]->end(); ++p) + if(p->first second=0.; + } +std::cout << "Error before permute = " <<(br-r0).norm()<::iterator p; + typename std::map *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; + } + delete r.v[k]; + r.v[k]=vnew; + } + +return r; +} + + }//namespace #endif //_SPARSESMAT_H_ diff --git a/t.cc b/t.cc index 183f900..4ce871d 100644 --- a/t.cc +++ b/t.cc @@ -1424,7 +1424,7 @@ NRMat r2(rx); cout <<"Error = "<<(r2-rd).norm()< > h; cin>>h; @@ -1442,5 +1442,137 @@ cout <<"errorx = "<<(r2-NRSMat >(r)).norm()<>n; +NRMat a(n,n); +a.randomize(1); +for(int i=0; i bb=a.transpose()*a; +NRMat cc(bb); +NRMat b(bb); +NRMat c(cc); +cholesky(b,0); +cout < > b(bb); +//cholesky(b,0); +//if(n<100)cout <<"Dense Cholesky result\n"< > c; +c= cc.cholesky(); +NRMat > cx(c); +if(n<100)cout <<"Sparse pivoted Cholesky result \n"< bb(n); +bb.gemm(0.,bh,'c',bh,'n',1.); +if(n<1000) cout <<"Input matrix\n"< br(n); +br.gemm(0.,b,'c',b,'n',1.); +if(n<1000) cout <<"Result of reconstruction\n"<