jiri/LA_library
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

*** empty log message ***

Browse Source
This commit is contained in:
jiri 2010-01-11 10:12:28 +00:00
parent 12c88e6872
commit 8ec7c11a6e
6 changed files with 92 additions and 42 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

10
la_traits.h
View File

@ -219,8 +219,9 @@ static void copy(complex<C> *dest, complex<C> *src, unsigned int n) {memcpy(dest
static void clear(complex<C> *dest, unsigned int n) {memset(dest,0,n*sizeof(complex<C>));} static void clear(complex<C> *dest, unsigned int n) {memset(dest,0,n*sizeof(complex<C>));}
static void copyonwrite(complex<C> &x) {}; static void copyonwrite(complex<C> &x) {};
static void clearme(complex<C> &x) {x=0;}; static void clearme(complex<C> &x) {x=0;};
static inline complex<C> conjugate(complex<C> &x) {return complex<C>(x.real(),-x.imag());}; static inline complex<C> conjugate(const complex<C> &x) {return complex<C>(x.real(),-x.imag());};
static inline C realpart(complex<C> &x) {return x.real();} static inline C realpart(const complex<C> &x) {return x.real();}
static inline C imagpart(const complex<C> &x) {return x.imag();}
}; };
//non-complex scalars //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 clear(C *dest, unsigned int n) {memset(dest,0,n*sizeof(C));}
static void copyonwrite(C &x) {}; static void copyonwrite(C &x) {};
static void clearme(complex<C> &x) {x=0;}; static void clearme(complex<C> &x) {x=0;};
static inline C conjugate(C&x) {return x;}; static inline C conjugate(const C &x) {return x;};
static inline C realpart(C &x) {return x;} static inline C realpart(const C &x) {return x;}
static inline C imagpart(const C &x) {return 0;}
}; };

2
mat.h
View File

@ -128,7 +128,7 @@ public:
const T trace() const; const T trace() const;
//members concerning sparse matrix //members concerning sparse matrix
SparseSMat<T> & operator*(const SparseSMat<T> &rhs) const; SparseSMat<T> operator*(const SparseSMat<T> &rhs) const;
explicit NRMat(const SparseMat<T> &rhs); // dense from sparse explicit NRMat(const SparseMat<T> &rhs); // dense from sparse
explicit NRMat(const SparseSMat<T> &rhs); // dense from sparse explicit NRMat(const SparseSMat<T> &rhs); // dense from sparse
NRMat & operator+=(const SparseMat<T> &rhs); NRMat & operator+=(const SparseMat<T> &rhs);

21
sparsesmat.cc
View File

@ -30,11 +30,11 @@ namespace LA {
//dense times sparse (not necessarily symmetric) //dense times sparse (not necessarily symmetric)
template <typename T> template <typename T>
SparseSMat<T> & NRMat<T>::operator*(const SparseSMat<T> &rhs) const SparseSMat<T> NRMat<T>::operator*(const SparseSMat<T> &rhs) const
{ {
SparseSMat<T> r(nn,rhs.ncols()); SparseSMat<T> r(nn,rhs.ncols());
if(mm!=rhs.nrows()) laerror("incompatible sizes in NRMat*SparseSMat"); if(mm!=rhs.nrows()) laerror("incompatible sizes in NRMat*SparseSMat");
for(SPMatindex k=0; k<nn; ++k) //summation loop for(SPMatindex k=0; k<mm; ++k) //summation loop
{ {
std::map<SPMatindex,T> * kl = rhs.line(k); std::map<SPMatindex,T> * kl = rhs.line(k);
if(kl) if(kl)
@ -71,6 +71,7 @@ void SparseSMat<T>::gemm(const T beta, const SparseSMat &a, const char transa, c
{ {
(*this) *= beta; (*this) *= beta;
if(alpha==(T)0) return; 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"); if(a.nn!=b.nn || a.nn!=nn) laerror("incompatible sizes in SparseSMat::gemm");
copyonwrite(); copyonwrite();
@ -89,13 +90,10 @@ for(SPMatindex k=0; k<nn; ++k) //summation loop
for(p=a.v[k]->begin(), i=0; p!=a.v[k]->end(); ++p,++i) { ai[i] = p->first; av[i] = LA_traits<T>::conjugate(p->second); } for(p=a.v[k]->begin(), i=0; p!=a.v[k]->end(); ++p,++i) { ai[i] = p->first; av[i] = LA_traits<T>::conjugate(p->second); }
else else
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=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<T>::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 //make multiply via blas
NRMat<T> prod=av.otimes(bv,false,alpha); NRMat<T> prod=av.otimes(bv,tolower(transb)=='c',alpha);
//scatter the results -- probably the computational bottleneck //scatter the results -- probably the computational bottleneck
for(i=0; i<prod.nrows(); ++i) for(j=0; j<prod.ncols(); ++j) for(i=0; i<prod.nrows(); ++i) for(j=0; j<prod.ncols(); ++j)
@ -127,7 +125,7 @@ return *this;
template <class T> template <class T>
void SparseSMat<T>::axpy(const T alpha, const SparseSMat &x, const bool transp) void SparseSMat<T>::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; if(alpha==(T)0) return;
copyonwrite(); copyonwrite();
for(SPMatindex i=0; i<nn; ++i) if(x.v[i]) for(SPMatindex i=0; i<nn; ++i) if(x.v[i])
@ -148,7 +146,8 @@ simplify();
template <class T> template <class T>
void SparseSMat<T>::gemv(const T beta, NRVec<T> &r, const char trans, const T alpha, const NRVec<T> &x) const void SparseSMat<T>::gemv(const T beta, NRVec<T> &r, const char trans, const T alpha, const NRVec<T> &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; r *= beta;
if(alpha == (T)0) return; if(alpha == (T)0) return;
r.copyonwrite(); r.copyonwrite();
@ -236,6 +235,7 @@ return std::sqrt(sum);
template <class T> template <class T>
const T* SparseSMat<T>::diagonalof(NRVec<T> &r, const bool divide, bool cache) const const T* SparseSMat<T>::diagonalof(NRVec<T> &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()"); if(nn!=r.size()) laerror("incompatible vector size in diagonalof()");
NRVec<T> *rr; NRVec<T> *rr;
@ -266,8 +266,7 @@ void SparseSMat<T>::get(int fd, bool dimen, bool transp) {
if(dimen) { if(dimen) {
if(2*sizeof(SPMatindex)!=read(fd,&dim,2*sizeof(SPMatindex))) laerror("read() error in SparseSMat::get "); 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],dim[1]);
resize(dim[0]);
} }
else copyonwrite(); else copyonwrite();
@ -288,7 +287,7 @@ void SparseSMat<T>::put(int fd, bool dimen, bool transp) const {
errno=0; errno=0;
if(dimen) { 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() 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<T>::iterator p(*this); typename SparseSMat<T>::iterator p(*this);

90
sparsesmat.h
View File

@ -50,20 +50,25 @@ class SparseSMat
{ {
protected: protected:
SPMatindex nn; SPMatindex nn;
SPMatindex mm;
std::map<SPMatindex,T> **v; std::map<SPMatindex,T> **v;
int *count; int *count;
public: public:
SparseSMat() : nn(0), v(NULL), count(NULL) {}; SparseSMat() : nn(0), mm(0), v(NULL), count(NULL) {};
explicit SparseSMat(const SPMatindex n); //prevent double -> int -> SparseSMat explicit SparseSMat(const SPMatindex n, const SPMatindex m); //prevent double -> int -> SparseSMat
explicit SparseSMat(const SPMatindex n);
SparseSMat(const SparseSMat &rhs); SparseSMat(const SparseSMat &rhs);
explicit SparseSMat(const SparseMat<T> &rhs); explicit SparseSMat(const SparseMat<T> &rhs);
explicit SparseSMat(const NRSMat<T> &rhs); explicit SparseSMat(const NRSMat<T> &rhs);
explicit SparseSMat(const NRMat<T> &rhs); explicit SparseSMat(const NRMat<T> &rhs);
SparseSMat & operator=(const SparseSMat &rhs); SparseSMat & operator=(const SparseSMat &rhs);
void copyonwrite(); void copyonwrite();
void resize(const SPMatindex n); void resize(const SPMatindex nn, const SPMatindex mm);
std::map<SPMatindex,T> *line(SPMatindex n) const {return v[n];}; inline void setcoldim(int i) {mm=(SPMatindex)i;};
void clear() {resize(nn);} //std::map<SPMatindex,T> *line(SPMatindex n) const {return v[n];};
typedef std::map<SPMatindex,T> *ROWTYPE;
inline typename SparseSMat<T>::ROWTYPE & operator[](const SPMatindex i) {return v[i];};
void clear() {resize(nn,mm);}
unsigned long long simplify(); unsigned long long simplify();
~SparseSMat(); ~SparseSMat();
inline int getcount() const {return count?*count:0;} inline int getcount() const {return count?*count:0;}
@ -83,7 +88,7 @@ public:
void gemv(const T beta, NRVec<T> &r, const char trans, const T alpha, const NRVec<T> &x) const; 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;}; 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; typename LA_traits<T>::normtype norm(const T scalar=(T)0) const;
inline const SparseSMat operator*(const SparseSMat &rhs) const {SparseSMat<T> 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<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
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*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 void add(const SPMatindex n, const SPMatindex m, const T elem, const bool both=true);
inline unsigned long long length() {return simplify();}; inline unsigned long long length() {return simplify();};
@ -92,7 +97,7 @@ public:
void get(int fd, bool dimen, bool transp); void get(int fd, bool dimen, bool transp);
void put(int fd, bool dimen, bool transp) const; void put(int fd, bool dimen, bool transp) const;
int nrows() const {return nn;} int nrows() const {return nn;}
int ncols() const {return nn;} int ncols() const {return mm;}
SparseSMat<T> cholesky(void) const; SparseSMat<T> cholesky(void) const;
class iterator {//not efficient, just for output to ostreams class iterator {//not efficient, just for output to ostreams
@ -103,14 +108,15 @@ public:
typename std::map<SPMatindex,T>::iterator *col; typename std::map<SPMatindex,T>::iterator *col;
typename std::map<SPMatindex,T>::iterator mycol; typename std::map<SPMatindex,T>::iterator mycol;
SPMatindex mynn; SPMatindex mynn;
SPMatindex mymm;
std::map<SPMatindex,T> **myv; std::map<SPMatindex,T> **myv;
public: public:
//compiler-generated iterator & operator=(const iterator &rhs); //compiler-generated iterator & operator=(const iterator &rhs);
//compiler-generated iterator(const iterator &rhs); //compiler-generated iterator(const iterator &rhs);
iterator(): p(NULL),row(0),col(NULL),mynn(0),myv(NULL) {}; iterator(): p(NULL),row(0),col(NULL),mynn(0),mymm(0),myv(NULL) {};
iterator(const SparseSMat &rhs) : mynn(rhs.nn), myv(rhs.v), col(NULL) {row=0; p= &my; operator++();} iterator(const SparseSMat &rhs) : mynn(rhs.nn), mymm(rhs.mm), myv(rhs.v), col(NULL) {row=0; p= &my; operator++();}
iterator operator++() { iterator operator++() {
if(col) //finish column list if(col) //finish column list
{ {
@ -156,10 +162,19 @@ public:
iterator end() const {return iterator(NULL);} iterator end() const {return iterator(NULL);}
}; };
template <typename T> template <typename T>
SparseSMat<T>::SparseSMat(const SPMatindex n) SparseSMat<T>::SparseSMat(const SPMatindex n)
:nn(n), :nn(n), mm(n),
count(new int(1))
{
v= new std::map<SPMatindex,T> * [n];
memset(v,0,nn*sizeof(std::map<SPMatindex,T> *));
}
template <typename T>
SparseSMat<T>::SparseSMat(const SPMatindex n, const SPMatindex m)
:nn(n), mm(m),
count(new int(1)) count(new int(1))
{ {
v= new std::map<SPMatindex,T> * [n]; v= new std::map<SPMatindex,T> * [n];
@ -168,7 +183,7 @@ memset(v,0,nn*sizeof(std::map<SPMatindex,T> *));
template <typename T> template <typename T>
SparseSMat<T>::SparseSMat(const NRSMat<T> &rhs) SparseSMat<T>::SparseSMat(const NRSMat<T> &rhs)
:nn(rhs.nrows()), :nn(rhs.nrows()), mm(rhs.ncols()),
count(new int(1)) count(new int(1))
{ {
v= new std::map<SPMatindex,T> * [nn]; v= new std::map<SPMatindex,T> * [nn];
@ -179,14 +194,14 @@ for(i=0; i<nn; ++i) for(j=0; j<=i; ++j) if(std::abs(rhs(i,j))>SPARSEEPSILON) (*t
template <typename T> template <typename T>
SparseSMat<T>::SparseSMat(const NRMat<T> &rhs) SparseSMat<T>::SparseSMat(const NRMat<T> &rhs)
:nn(rhs.nrows()), :nn(rhs.nrows()), mm(rhs.ncols()),
count(new int(1)) count(new int(1))
{ {
if(rhs.nrows()!=rhs.ncols()) laerror("non-square matrix in SparseSMat constructor from NRMat"); if(rhs.nrows()!=rhs.ncols()) laerror("non-square matrix in SparseSMat constructor from NRMat");
v= new std::map<SPMatindex,T> * [nn]; v= new std::map<SPMatindex,T> * [nn];
memset(v,0,nn*sizeof(std::map<SPMatindex,T> *)); memset(v,0,nn*sizeof(std::map<SPMatindex,T> *));
int i,j; int i,j;
for(i=0; i<nn; ++i) for(j=0; j<nn; ++j) if(std::abs(rhs(i,j))>SPARSEEPSILON) (*this).add(i,j,rhs(i,j),false); for(i=0; i<nn; ++i) for(j=0; j<mm; ++j) if(std::abs(rhs(i,j))>SPARSEEPSILON) (*this).add(i,j,rhs(i,j),false);
} }
@ -195,6 +210,7 @@ SparseSMat<T>::SparseSMat(const SparseSMat &rhs)
{ {
v = rhs.v; v = rhs.v;
nn = rhs.nn; nn = rhs.nn;
mm = rhs.mm;
count = rhs.count; count = rhs.count;
if(count) (*count)++; if(count) (*count)++;
} }
@ -205,6 +221,7 @@ template <typename T>
NRSMat<T>::NRSMat(const SparseSMat<T> &rhs) NRSMat<T>::NRSMat(const SparseSMat<T> &rhs)
: nn(rhs.nrows()) : nn(rhs.nrows())
{ {
if(rhs.nrows()!=rhs.ncols()) laerror("cannot transform rectangular matrix to NRSMat");
count = new int(1); count = new int(1);
v=new T[nn2]; v=new T[nn2];
memset(v,0,nn2*sizeof(T)); memset(v,0,nn2*sizeof(T));
@ -251,13 +268,14 @@ SparseSMat<T>::~SparseSMat()
template <typename T> template <typename T>
void SparseSMat<T>::resize(const SPMatindex n) void SparseSMat<T>::resize(const SPMatindex n, const SPMatindex m)
{ {
if(!count) if(!count)
{ {
if(n==0) return; if(n==0) return;
count = new int(1); count = new int(1);
nn=n; nn=n;
mm=m;
v= new std::map<SPMatindex,T> * [nn]; v= new std::map<SPMatindex,T> * [nn];
for(SPMatindex i=0; i<nn; ++i) v[i]=NULL; for(SPMatindex i=0; i<nn; ++i) v[i]=NULL;
return; return;
@ -270,13 +288,15 @@ if(*count > 1) //it was shared
{ {
count = new int(1); count = new int(1);
nn=n; nn=n;
mm=m;
v= new std::map<SPMatindex,T> * [nn]; v= new std::map<SPMatindex,T> * [nn];
for(SPMatindex i=0; i<nn; ++i) v[i]=NULL; for(SPMatindex i=0; i<nn; ++i) v[i]=NULL;
} }
else {v=NULL; nn=0; count=NULL;} else {v=NULL; nn=0; mm=0; count=NULL;}
} }
else //it was not shared else //it was not shared
{ {
mm=m;
//delete all trees //delete all trees
for(SPMatindex i=0; i<nn; ++i) if(v[i]) {delete v[i]; v[i]=NULL;} for(SPMatindex i=0; i<nn; ++i) if(v[i]) {delete v[i]; v[i]=NULL;}
if(n!=nn) if(n!=nn)
@ -305,6 +325,7 @@ SparseSMat<T> & SparseSMat<T>::operator=(const SparseSMat &rhs)
} }
v = rhs.v; v = rhs.v;
nn = rhs.nn; nn = rhs.nn;
mm = rhs.mm;
count = rhs.count; count = rhs.count;
if(count) (*count)++; if(count) (*count)++;
} }
@ -335,7 +356,7 @@ template <typename T>
void SparseSMat<T>::add(const SPMatindex n, const SPMatindex m, const T elem, const bool both) void SparseSMat<T>::add(const SPMatindex n, const SPMatindex m, const T elem, const bool both)
{ {
#ifdef DEBUG #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 #endif
if(!v[n]) v[n] = new std::map<SPMatindex,T>; if(!v[n]) v[n] = new std::map<SPMatindex,T>;
@ -377,8 +398,7 @@ std::ostream & operator<<(std::ostream &s, const SparseSMat<T> &x)
{ {
SPMatindex n; SPMatindex n;
n = x.nrows(); s << x.nrows() << " "<< x.ncols()<< std::endl;
s << n << " "<< n<< std::endl;
typename SparseSMat<T>::iterator p(x); typename SparseSMat<T>::iterator p(x);
for(; p.notend(); ++p) s << (int)p->row << ' ' << (int)p->col << ' ' << (typename LA_traits_io<T>::IOtype) p->elem << '\n'; for(; p.notend(); ++p) s << (int)p->row << ' ' << (int)p->col << ' ' << (typename LA_traits_io<T>::IOtype) p->elem << '\n';
@ -393,28 +413,31 @@ std::istream& operator>>(std::istream &s, SparseSMat<T> &x)
long i,j; long i,j;
s >> n >> m; s >> n >> m;
if(n!=m) laerror("SparseSMat must be square"); if(n!=m) laerror("SparseSMat must be square");
x.resize(n); x.resize(n,m);
s >> i >> j; s >> i >> j;
typename LA_traits_io<T>::IOtype tmp; typename LA_traits_io<T>::IOtype tmp;
while(i>=0 && j>=0) while(i>=0 && j>=0)
{ {
s>>tmp; s>>tmp;
if(i>=n||j>=m) laerror("bad index in SparseSMat::operator>>");
x.add(i,j,tmp,false); x.add(i,j,tmp,false);
s >> i >> j; s >> i >> j;
} }
return s; return s;
} }
template <typename T> template <typename T>
SparseSMat<T> SparseSMat<T>::transpose(bool conj) const SparseSMat<T> SparseSMat<T>::transpose(bool conj) const
{ {
SparseSMat<T> r(nn); SparseSMat<T> r(mm,nn);
typename SparseSMat<T>::iterator p(*this); typename SparseSMat<T>::iterator p(*this);
for(; p.notend(); ++p) r.add(p->col, p->row, (conj?LA_traits<T>::conjugate(p->elem):p->elem), false); for(; p.notend(); ++p) r.add(p->col, p->row, (conj?LA_traits<T>::conjugate(p->elem):p->elem), false);
return r; return r;
} }
//Cholesky decomposition, pivoted, positive semidefinite, not in place //Cholesky decomposition, pivoted, positive semidefinite, not in place
//it is NOT checked that the input matrix is symmetric/hermitean //it is NOT checked that the input matrix is symmetric/hermitean
//result.transpose(true)*result reproduces the original matrix //result.transpose(true)*result reproduces the original matrix
@ -423,7 +446,7 @@ return r;
template <typename T> template <typename T>
SparseSMat<T> SparseSMat<T>::cholesky(void) const SparseSMat<T> SparseSMat<T>::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 //we need real values for sorting, if T is already real it makes just an unnecessary copy of one vector
NRVec<typename LA_traits<T>::normtype> diagreal(nn); NRVec<typename LA_traits<T>::normtype> diagreal(nn);
{ {
@ -451,6 +474,7 @@ for(int i=0; i<nn; ++i) invpivot[pivot[i]]=i;
//copy-permute upper triangle //copy-permute upper triangle
SparseSMat<T> r; SparseSMat<T> r;
r.nn=nn; r.nn=nn;
r.mm=nn;
r.count = new int(1); r.count = new int(1);
r.v = new std::map<SPMatindex,T> * [nn]; r.v = new std::map<SPMatindex,T> * [nn];
for(SPMatindex i=0; i<nn; ++i) for(SPMatindex i=0; i<nn; ++i)
@ -545,5 +569,27 @@ return r;
} }
//outer product expected to be sparse
template<typename T>
SparseSMat<T> otimes_sparse(const NRVec<T> &lhs, const NRVec<T> &rhs, const bool conjugate=false, const T &scale=1)
{
SparseSMat<T> r(lhs.size(),rhs.size());
for(SPMatindex i=0; i<lhs.size(); ++i)
if(lhs[i])
{
for(SPMatindex j=0; j<rhs.size(); ++j)
if(rhs[j])
{
T x=lhs[i]*(conjugate?LA_traits<T>::conjugate(rhs[j]):rhs[j])*scale;
if(std::abs(x)>SPARSEEPSILON) r.add(i,j,x);
}
}
return r;
}
}//namespace }//namespace
#endif //_SPARSESMAT_H_ #endif //_SPARSESMAT_H_

6
t.cc
View File

@ -1555,17 +1555,17 @@ if(1)
{ {
int n; int n;
cin >>n; cin >>n;
SparseSMat<double> bh(n); SparseSMat<double> bh(n,n);
for(int i=0; i<=n/400; ++i) for(int j=i; j<n; ++j) {if((double)random()/RAND_MAX>0.995 || i==j) for(int i=0; i<=n/400; ++i) for(int j=i; j<n; ++j) {if((double)random()/RAND_MAX>0.995 || i==j)
{bh.add(i,j,(double)random()/RAND_MAX*(i==j?10.:(random()>RAND_MAX/2?1:-1)),false);}} {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"<<bh; if(n<1000) cout <<"Random matrix\n"<<bh;
SparseSMat<double> bb(n); SparseSMat<double> bb(n,n);
bb.gemm(0.,bh,'c',bh,'n',1.); bb.gemm(0.,bh,'c',bh,'n',1.);
if(n<1000) cout <<"Input matrix\n"<<bb; if(n<1000) cout <<"Input matrix\n"<<bb;
cout <<"Original filling = "<<bb.simplify()<<endl; cout <<"Original filling = "<<bb.simplify()<<endl;
SparseSMat<double> b = bb.cholesky(); SparseSMat<double> b = bb.cholesky();
if(n<1000) cout <<"Cholesky result\n"<<b; if(n<1000) cout <<"Cholesky result\n"<<b;
SparseSMat<double> br(n); SparseSMat<double> br(n,n);
br.gemm(0.,b,'c',b,'n',1.); br.gemm(0.,b,'c',b,'n',1.);
if(n<1000) cout <<"Result of reconstruction\n"<<br; if(n<1000) cout <<"Result of reconstruction\n"<<br;
if(n<1000) cout <<"Difference\n"<<br-bb; if(n<1000) cout <<"Difference\n"<<br-bb;

3
vec.h
View File

@ -113,6 +113,8 @@ public:
inline const T dot(const T *a, const int stride=1) const; // ddot with a stride-vector inline const T dot(const T *a, const int stride=1) const; // ddot with a stride-vector
inline T & operator[](const int i); inline T & operator[](const int i);
inline const T & operator[](const int i) const; inline const T & operator[](const int i) const;
typedef T ROWTYPE;
inline void setcoldim(int i) {}; //dummy
inline int size() const; inline int size() const;
inline operator T*(); //get a pointer to the data inline operator T*(); //get a pointer to the data
inline operator const T*() const; //get a pointer to the data inline operator const T*() const; //get a pointer to the data
@ -145,6 +147,7 @@ public:
#include "mat.h" #include "mat.h"
#include "smat.h" #include "smat.h"
#include "sparsemat.h" #include "sparsemat.h"
#include "sparsesmat.h"
namespace LA { namespace LA {
// formatted I/O // formatted I/O