From 8ec7c11a6ef622735cf97deedfbe9e15851b6335 Mon Sep 17 00:00:00 2001 From: jiri Date: Mon, 11 Jan 2010 10:12:28 +0000 Subject: [PATCH] *** empty log message *** --- la_traits.h | 10 +++--- mat.h | 2 +- sparsesmat.cc | 23 +++++++------ sparsesmat.h | 90 ++++++++++++++++++++++++++++++++++++++------------- t.cc | 6 ++-- vec.h | 3 ++ 6 files changed, 92 insertions(+), 42 deletions(-) diff --git a/la_traits.h b/la_traits.h index a45817e..60ec83a 100644 --- a/la_traits.h +++ b/la_traits.h @@ -219,8 +219,9 @@ 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();} +static inline complex conjugate(const complex &x) {return complex(x.real(),-x.imag());}; +static inline C realpart(const complex &x) {return x.real();} +static inline C imagpart(const complex &x) {return x.imag();} }; //non-complex scalars @@ -243,8 +244,9 @@ 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;} +static inline C conjugate(const C &x) {return x;}; +static inline C realpart(const C &x) {return x;} +static inline C imagpart(const C &x) {return 0;} }; diff --git a/mat.h b/mat.h index c6f6822..36eabf1 100644 --- a/mat.h +++ b/mat.h @@ -128,7 +128,7 @@ public: const T trace() const; //members concerning sparse matrix - SparseSMat & operator*(const SparseSMat &rhs) const; + 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); diff --git a/sparsesmat.cc b/sparsesmat.cc index bd9f95c..5720d9e 100644 --- a/sparsesmat.cc +++ b/sparsesmat.cc @@ -30,11 +30,11 @@ namespace LA { //dense times sparse (not necessarily symmetric) template -SparseSMat & NRMat::operator*(const SparseSMat &rhs) const +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) @@ -71,6 +71,7 @@ void SparseSMat::gemm(const T beta, const SparseSMat &a, const char transa, c { (*this) *= beta; if(alpha==(T)0) return; +if(a.nn!=a.mm || b.nn!=b.mm || nn!=mm) laerror("SparseSMat::gemm implemented only for square symmetric matrices"); if(a.nn!=b.nn || a.nn!=nn) laerror("incompatible sizes in SparseSMat::gemm"); copyonwrite(); @@ -89,13 +90,10 @@ for(SPMatindex k=0; kbegin(), 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; } + 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); + NRMat prod=av.otimes(bv,tolower(transb)=='c',alpha); //scatter the results -- probably the computational bottleneck for(i=0; i void SparseSMat::axpy(const T alpha, const SparseSMat &x, const bool transp) { -if(nn!=x.nn) laerror("incompatible matrix dimensions in SparseSMat::axpy"); +if(nn!=x.nn || mm!=x.mm) laerror("incompatible matrix dimensions in SparseSMat::axpy"); if(alpha==(T)0) return; copyonwrite(); for(SPMatindex i=0; i void SparseSMat::gemv(const T beta, NRVec &r, const char trans, const T alpha, const NRVec &x) const { -if(nn!=r.size() || nn!= x.size()) laerror("incompatible matrix vector dimensions in SparseSMat::gemv"); +if(nn!=r.size() || mm!= x.size()) laerror("incompatible matrix vector dimensions in SparseSMat::gemv"); +if(trans) laerror("transposition not implemented yet in SparseSMat::gemv"); r *= beta; if(alpha == (T)0) return; r.copyonwrite(); @@ -236,6 +235,7 @@ return std::sqrt(sum); template const T* SparseSMat::diagonalof(NRVec &r, const bool divide, bool cache) const { +if(nn!=mm) laerror("non-square matrix in SparseSMat::diagonalof"); if(nn!=r.size()) laerror("incompatible vector size in diagonalof()"); NRVec *rr; @@ -266,8 +266,7 @@ void SparseSMat::get(int fd, bool dimen, bool transp) { 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]); + resize(dim[0],dim[1]); } else copyonwrite(); @@ -288,7 +287,7 @@ 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"); + if(sizeof(SPMatindex)!=write(fd,&mm,sizeof(SPMatindex))) laerror("cannot write() 2 in SparseSMat::put"); } typename SparseSMat::iterator p(*this); diff --git a/sparsesmat.h b/sparsesmat.h index 25ca17d..66ef6b9 100644 --- a/sparsesmat.h +++ b/sparsesmat.h @@ -50,20 +50,25 @@ class SparseSMat { protected: SPMatindex nn; + SPMatindex mm; std::map **v; int *count; public: - SparseSMat() : nn(0), v(NULL), count(NULL) {}; - explicit SparseSMat(const SPMatindex n); //prevent double -> int -> SparseSMat + SparseSMat() : nn(0), mm(0), v(NULL), count(NULL) {}; + explicit SparseSMat(const SPMatindex n, const SPMatindex m); //prevent double -> int -> SparseSMat + explicit SparseSMat(const SPMatindex n); 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);} + void resize(const SPMatindex nn, const SPMatindex mm); + inline void setcoldim(int i) {mm=(SPMatindex)i;}; + //std::map *line(SPMatindex n) const {return v[n];}; + typedef std::map *ROWTYPE; + inline typename SparseSMat::ROWTYPE & operator[](const SPMatindex i) {return v[i];}; + void clear() {resize(nn,mm);} unsigned long long simplify(); ~SparseSMat(); inline int getcount() const {return count?*count:0;} @@ -83,7 +88,7 @@ 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); r.gemm(0,*this,'n',rhs,'n',1); return r;}; //!!!NOT A GENERAL ROUTINE, JUST FOR THE CASES WHEN THE RESULT STAYS SYMMETRIC + 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 inline void add(const SPMatindex n, const SPMatindex m, const T elem, const bool both=true); inline unsigned long long length() {return simplify();}; @@ -92,7 +97,7 @@ public: 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;} + int ncols() const {return mm;} SparseSMat cholesky(void) const; class iterator {//not efficient, just for output to ostreams @@ -103,14 +108,15 @@ public: typename std::map::iterator *col; typename std::map::iterator mycol; SPMatindex mynn; + SPMatindex mymm; std::map **myv; public: //compiler-generated iterator & operator=(const iterator &rhs); //compiler-generated iterator(const iterator &rhs); - iterator(): p(NULL),row(0),col(NULL),mynn(0),myv(NULL) {}; - iterator(const SparseSMat &rhs) : mynn(rhs.nn), myv(rhs.v), col(NULL) {row=0; p= &my; operator++();} + iterator(): p(NULL),row(0),col(NULL),mynn(0),mymm(0),myv(NULL) {}; + iterator(const SparseSMat &rhs) : mynn(rhs.nn), mymm(rhs.mm), myv(rhs.v), col(NULL) {row=0; p= &my; operator++();} iterator operator++() { if(col) //finish column list { @@ -156,10 +162,19 @@ public: iterator end() const {return iterator(NULL);} }; - template SparseSMat::SparseSMat(const SPMatindex n) -:nn(n), +:nn(n), mm(n), +count(new int(1)) +{ +v= new std::map * [n]; +memset(v,0,nn*sizeof(std::map *)); +} + + +template +SparseSMat::SparseSMat(const SPMatindex n, const SPMatindex m) +:nn(n), mm(m), count(new int(1)) { v= new std::map * [n]; @@ -168,7 +183,7 @@ memset(v,0,nn*sizeof(std::map *)); template SparseSMat::SparseSMat(const NRSMat &rhs) -:nn(rhs.nrows()), +:nn(rhs.nrows()), mm(rhs.ncols()), count(new int(1)) { v= new std::map * [nn]; @@ -179,14 +194,14 @@ for(i=0; iSPARSEEPSILON) (*t template SparseSMat::SparseSMat(const NRMat &rhs) -:nn(rhs.nrows()), +: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; -for(i=0; iSPARSEEPSILON) (*this).add(i,j,rhs(i,j),false); +for(i=0; iSPARSEEPSILON) (*this).add(i,j,rhs(i,j),false); } @@ -195,6 +210,7 @@ SparseSMat::SparseSMat(const SparseSMat &rhs) { v = rhs.v; nn = rhs.nn; +mm = rhs.mm; count = rhs.count; if(count) (*count)++; } @@ -205,6 +221,7 @@ template NRSMat::NRSMat(const SparseSMat &rhs) : nn(rhs.nrows()) { +if(rhs.nrows()!=rhs.ncols()) laerror("cannot transform rectangular matrix to NRSMat"); count = new int(1); v=new T[nn2]; memset(v,0,nn2*sizeof(T)); @@ -251,13 +268,14 @@ SparseSMat::~SparseSMat() template -void SparseSMat::resize(const SPMatindex n) +void SparseSMat::resize(const SPMatindex n, const SPMatindex m) { if(!count) { if(n==0) return; count = new int(1); nn=n; + mm=m; v= new std::map * [nn]; for(SPMatindex i=0; i 1) //it was shared { count = new int(1); nn=n; + mm=m; v= new std::map * [nn]; for(SPMatindex i=0; i & SparseSMat::operator=(const SparseSMat &rhs) } v = rhs.v; nn = rhs.nn; + mm = rhs.mm; count = rhs.count; if(count) (*count)++; } @@ -335,7 +356,7 @@ template void SparseSMat::add(const SPMatindex n, const SPMatindex m, const T elem, const bool both) { #ifdef DEBUG -if(n>=nn || m>=nn) laerror("illegal index in SparseSMat::add()"); +if(n>=nn || m>=mm) laerror("illegal index in SparseSMat::add()"); #endif if(!v[n]) v[n] = new std::map; @@ -377,8 +398,7 @@ std::ostream & operator<<(std::ostream &s, const SparseSMat &x) { SPMatindex n; -n = x.nrows(); -s << n << " "<< n<< std::endl; +s << x.nrows() << " "<< x.ncols()<< std::endl; typename SparseSMat::iterator p(x); for(; p.notend(); ++p) s << (int)p->row << ' ' << (int)p->col << ' ' << (typename LA_traits_io::IOtype) p->elem << '\n'; @@ -393,28 +413,31 @@ std::istream& operator>>(std::istream &s, SparseSMat &x) long i,j; s >> n >> m; if(n!=m) laerror("SparseSMat must be square"); - x.resize(n); + x.resize(n,m); s >> i >> j; typename LA_traits_io::IOtype tmp; while(i>=0 && j>=0) { s>>tmp; + if(i>=n||j>=m) laerror("bad index in SparseSMat::operator>>"); x.add(i,j,tmp,false); s >> i >> j; } return s; } + template SparseSMat SparseSMat::transpose(bool conj) const { -SparseSMat r(nn); +SparseSMat r(mm,nn); typename SparseSMat::iterator p(*this); for(; p.notend(); ++p) r.add(p->col, p->row, (conj?LA_traits::conjugate(p->elem):p->elem), false); return r; } + //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 @@ -423,7 +446,7 @@ return r; template SparseSMat SparseSMat::cholesky(void) const { - +if(nn!=mm) laerror("Cholesky defined only for square matrices"); //we need real values for sorting, if T is already real it makes just an unnecessary copy of one vector NRVec::normtype> diagreal(nn); { @@ -451,6 +474,7 @@ for(int i=0; i r; r.nn=nn; +r.mm=nn; r.count = new int(1); r.v = new std::map * [nn]; for(SPMatindex i=0; i +SparseSMat otimes_sparse(const NRVec &lhs, const NRVec &rhs, const bool conjugate=false, const T &scale=1) +{ +SparseSMat r(lhs.size(),rhs.size()); +for(SPMatindex i=0; i::conjugate(rhs[j]):rhs[j])*scale; + if(std::abs(x)>SPARSEEPSILON) r.add(i,j,x); + } + } +return r; +} + + + + }//namespace #endif //_SPARSESMAT_H_ diff --git a/t.cc b/t.cc index 4ce871d..f683a95 100644 --- a/t.cc +++ b/t.cc @@ -1555,17 +1555,17 @@ if(1) { int n; cin >>n; -SparseSMat bh(n); +SparseSMat bh(n,n); for(int i=0; i<=n/400; ++i) for(int j=i; j0.995 || i==j) {bh.add(i,j,(double)random()/RAND_MAX*(i==j?10.:(random()>RAND_MAX/2?1:-1)),false);}} if(n<1000) cout <<"Random matrix\n"< bb(n); +SparseSMat bb(n,n); bb.gemm(0.,bh,'c',bh,'n',1.); if(n<1000) cout <<"Input matrix\n"< br(n); +SparseSMat br(n,n); br.gemm(0.,b,'c',b,'n',1.); if(n<1000) cout <<"Result of reconstruction\n"<