00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019 #ifndef _SPARSESMAT_H_
00020 #define _SPARSESMAT_H_
00021
00022 #include <string>
00023 #include <cmath>
00024 #include <stdlib.h>
00025 #include <sys/types.h>
00026 #include <sys/stat.h>
00027 #include <fcntl.h>
00028 #include <errno.h>
00029 #include "la_traits.h"
00030 #include "sparsemat.h"
00031 #include "vec.h"
00032 #include "mat.h"
00033 #include "smat.h"
00034 #include "qsort.h"
00035
00036 #include <map>
00037 #include <list>
00038
00039 #define CHOLESKYEPSILON 1e-16
00040
00041 namespace LA {
00042
00043
00044
00045
00046
00047
00048 template <typename T>
00049 class SparseSMat
00050 {
00051 protected:
00052 SPMatindex nn;
00053 SPMatindex mm;
00054 std::map<SPMatindex,T> **v;
00055 int *count;
00056 public:
00057 SparseSMat() : nn(0), mm(0), v(NULL), count(NULL) {};
00058 explicit SparseSMat(const SPMatindex n, const SPMatindex m);
00059 explicit SparseSMat(const SPMatindex n);
00060 SparseSMat(const SparseSMat &rhs);
00061 explicit SparseSMat(const SparseMat<T> &rhs);
00062 explicit SparseSMat(const NRSMat<T> &rhs);
00063 explicit SparseSMat(const NRMat<T> &rhs);
00064 SparseSMat & operator=(const SparseSMat &rhs);
00065 void copyonwrite();
00066 void resize(const SPMatindex nn, const SPMatindex mm);
00067 inline void setcoldim(int i) {mm=(SPMatindex)i;};
00068
00069 typedef std::map<SPMatindex,T> *ROWTYPE;
00070 inline typename SparseSMat<T>::ROWTYPE & operator[](const SPMatindex i) {return v[i];};
00071 void clear() {resize(nn,mm);}
00072 unsigned long long simplify();
00073 ~SparseSMat();
00074 inline int getcount() const {return count?*count:0;}
00075 SparseSMat & operator*=(const T &a);
00076 inline const SparseSMat operator*(const T &rhs) const {return SparseSMat(*this) *= rhs;}
00077 inline const SparseSMat operator+(const T &rhs) const {return SparseSMat(*this) += rhs;}
00078 inline const SparseSMat operator-(const T &rhs) const {return SparseSMat(*this) -= rhs;}
00079 inline const SparseSMat operator+(const SparseSMat &rhs) const {return SparseSMat(*this) += rhs;}
00080 inline const SparseSMat operator-(const SparseSMat &rhs) const {return SparseSMat(*this) -= rhs;}
00081 SparseSMat & operator=(const T &a);
00082 SparseSMat & operator+=(const T &a);
00083 SparseSMat & operator-=(const T &a);
00084 void axpy(const T alpha, const SparseSMat &x, const bool transp=0);
00085 inline SparseSMat & operator+=(const SparseSMat &rhs) {axpy((T)1,rhs); return *this;};
00086 inline SparseSMat & operator-=(const SparseSMat &rhs) {axpy((T)-1,rhs); return *this;};
00087 const T* diagonalof(NRVec<T> &, const bool divide=0, bool cache=false) const;
00088 void gemv(const T beta, NRVec<T> &r, const char trans, const T alpha, const NRVec<T> &x) const;
00089 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;};
00090 typename LA_traits<T>::normtype norm(const T scalar=(T)0) const;
00091 inline const SparseSMat operator*(const SparseSMat &rhs) const {SparseSMat<T> r(nn,mm); r.gemm(0,*this,'n',rhs,'n',1); return r;};
00092 void gemm(const T beta, const SparseSMat &a, const char transa, const SparseSMat &b, const char transb, const T alpha);
00093 inline void add(const SPMatindex n, const SPMatindex m, const T elem, const bool both=true);
00094 inline unsigned long long length() {return simplify();};
00095 void transposeme() const {laerror("in-place transposition not necessary/implemented for SparseSMat");};
00096 SparseSMat transpose(bool conj=false) const;
00097 inline bool issymmetric() const {return true;}
00098 void get(int fd, bool dimen, bool transp);
00099 void put(int fd, bool dimen, bool transp) const;
00100 int nrows() const {return nn;}
00101 int ncols() const {return mm;}
00102 SparseSMat<T> cholesky(void) const;
00103
00104 class iterator {
00105 private:
00106 matel<T> *p;
00107 matel<T> my;
00108 SPMatindex row;
00109 typename std::map<SPMatindex,T>::iterator *col;
00110 typename std::map<SPMatindex,T>::iterator mycol;
00111 SPMatindex mynn;
00112 SPMatindex mymm;
00113 std::map<SPMatindex,T> **myv;
00114
00115
00116 public:
00117
00118
00119 iterator(): p(NULL),row(0),col(NULL),mynn(0),mymm(0),myv(NULL) {};
00120 iterator(const SparseSMat &rhs) : mynn(rhs.nn), mymm(rhs.mm), myv(rhs.v), col(NULL) {row=0; p= &my; operator++();}
00121 iterator operator++() {
00122 if(col)
00123 {
00124 ++mycol;
00125 if(mycol != myv[row]->end())
00126 {
00127 p->row = row;
00128 p->col = mycol->first;
00129 p->elem = mycol->second;
00130 return *this;
00131 }
00132 else
00133 {
00134 col=NULL;
00135 ++row; if(row==mynn) {p=NULL; return *this;}
00136 }
00137 }
00138 nextrow:
00139 while(myv[row]==NULL) {++row; if(row==mynn) {p=NULL; return *this;}}
00140
00141
00142 col = &mycol;
00143 mycol = myv[row]->begin();
00144 if(mycol == myv[row]->end()) {col=NULL;
00145 ++row;
00146 if(row==mynn) {p=NULL; return *this;} else goto nextrow;
00147 }
00148
00149 p->row = row;
00150 p->col = mycol->first;
00151 p->elem = mycol->second;
00152 return *this;
00153 };
00154 iterator(matel<T> *q) {p=q; col=NULL;}
00155
00156 bool operator!=(const iterator &rhs) const {return p!=rhs.p;}
00157 bool operator==(const iterator &rhs) const {return p==rhs.p;}
00158 matel<T> & operator*() const {return *p;}
00159 matel<T> * operator->() const {return p;}
00160 bool notend() const {return (bool)p;};
00161 };
00162 iterator begin() const {return iterator(*this);}
00163 iterator end() const {return iterator(NULL);}
00164 };
00165
00166 template <typename T>
00167 SparseSMat<T>::SparseSMat(const SPMatindex n)
00168 :nn(n), mm(n),
00169 count(new int(1))
00170 {
00171 v= new std::map<SPMatindex,T> * [n];
00172 memset(v,0,nn*sizeof(std::map<SPMatindex,T> *));
00173 }
00174
00175
00176 template <typename T>
00177 SparseSMat<T>::SparseSMat(const SPMatindex n, const SPMatindex m)
00178 :nn(n), mm(m),
00179 count(new int(1))
00180 {
00181 v= new std::map<SPMatindex,T> * [n];
00182 memset(v,0,nn*sizeof(std::map<SPMatindex,T> *));
00183 }
00184
00185 template <typename T>
00186 SparseSMat<T>::SparseSMat(const NRSMat<T> &rhs)
00187 :nn(rhs.nrows()), mm(rhs.ncols()),
00188 count(new int(1))
00189 {
00190 v= new std::map<SPMatindex,T> * [nn];
00191 memset(v,0,nn*sizeof(std::map<SPMatindex,T> *));
00192 int i,j;
00193 for(i=0; i<nn; ++i) for(j=0; j<=i; ++j) if(std::abs(rhs(i,j))>SPARSEEPSILON) (*this).add(i,j,rhs(i,j),true);
00194 }
00195
00196 template <typename T>
00197 SparseSMat<T>::SparseSMat(const NRMat<T> &rhs)
00198 :nn(rhs.nrows()), mm(rhs.ncols()),
00199 count(new int(1))
00200 {
00201 if(rhs.nrows()!=rhs.ncols()) laerror("non-square matrix in SparseSMat constructor from NRMat");
00202 v= new std::map<SPMatindex,T> * [nn];
00203 memset(v,0,nn*sizeof(std::map<SPMatindex,T> *));
00204 int i,j;
00205 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);
00206 }
00207
00208
00209 template <typename T>
00210 SparseSMat<T>::SparseSMat(const SparseSMat &rhs)
00211 {
00212 v = rhs.v;
00213 nn = rhs.nn;
00214 mm = rhs.mm;
00215 count = rhs.count;
00216 if(count) (*count)++;
00217 }
00218
00219
00220 #define nn2 (nn*(nn+1)/2)
00221 template <typename T>
00222 NRSMat<T>::NRSMat(const SparseSMat<T> &rhs)
00223 : nn(rhs.nrows())
00224 {
00225 if(rhs.nrows()!=rhs.ncols()) laerror("cannot transform rectangular matrix to NRSMat");
00226 count = new int(1);
00227 v=new T[nn2];
00228 memset(v,0,nn2*sizeof(T));
00229 typename SparseSMat<T>::iterator p(rhs);
00230 for(; p.notend(); ++p) if(p->row <= p->col) (*this)(p->row,p->col)=p->elem;
00231 }
00232 #undef nn2
00233
00234
00235
00236 template <typename T>
00237 NRMat<T>::NRMat(const SparseSMat<T> &rhs) :
00238 nn(rhs.nrows()),
00239 mm(rhs.ncols()),
00240 count(new int(1))
00241 {
00242 #ifdef MATPTR
00243 v = new T*[nn];
00244 v[0] = new T[mm*nn];
00245 for (int i=1; i<nn; i++) v[i] = v[i-1] + mm;
00246 #else
00247 v = new T[mm*nn];
00248 #endif
00249 memset(&(*this)(0,0),0,mm*nn*sizeof(T));
00250 typename SparseSMat<T>::iterator p(rhs);
00251 for(; p.notend(); ++p) (*this)(p->row,p->col)= p->elem;
00252 }
00253
00254
00255
00256 template <typename T>
00257 SparseSMat<T>::~SparseSMat()
00258 {
00259 if(!count) return;
00260 if(--(*count) <= 0) {
00261 if(v)
00262 {
00263 for(SPMatindex i=0; i<nn; ++i) if(v[i]) delete v[i];
00264 delete[] (v);
00265 }
00266 delete count;
00267 }
00268 }
00269
00270
00271 template <typename T>
00272 void SparseSMat<T>::resize(const SPMatindex n, const SPMatindex m)
00273 {
00274 if(!count)
00275 {
00276 if(n==0) return;
00277 count = new int(1);
00278 nn=n;
00279 mm=m;
00280 v= new std::map<SPMatindex,T> * [nn];
00281 for(SPMatindex i=0; i<nn; ++i) v[i]=NULL;
00282 return;
00283 }
00284
00285 if(*count > 1)
00286 {
00287 (*count)--;
00288 if(n)
00289 {
00290 count = new int(1);
00291 nn=n;
00292 mm=m;
00293 v= new std::map<SPMatindex,T> * [nn];
00294 for(SPMatindex i=0; i<nn; ++i) v[i]=NULL;
00295 }
00296 else {v=NULL; nn=0; mm=0; count=NULL;}
00297 }
00298 else
00299 {
00300 mm=m;
00301
00302 for(SPMatindex i=0; i<nn; ++i) if(v[i]) {delete v[i]; v[i]=NULL;}
00303 if(n!=nn)
00304 {
00305 nn=n;
00306 for(SPMatindex i=0; i<nn; ++i) v[i]=NULL;
00307 }
00308 }
00309 }
00310
00311
00312 template <typename T>
00313 SparseSMat<T> & SparseSMat<T>::operator=(const SparseSMat &rhs)
00314 {
00315 if (this != &rhs)
00316 {
00317 if(count)
00318 if(--(*count) == 0)
00319 {
00320 if(v)
00321 {
00322 for(SPMatindex i=0; i<nn; ++i) if(v[i]) delete v[i];
00323 delete[] (v);
00324 }
00325 delete count;
00326 }
00327 v = rhs.v;
00328 nn = rhs.nn;
00329 mm = rhs.mm;
00330 count = rhs.count;
00331 if(count) (*count)++;
00332 }
00333 return *this;
00334 }
00335
00336
00337 template <typename T>
00338 void SparseSMat<T>::copyonwrite()
00339 {
00340 if(!count) laerror("SparseSmat::copyonwrite() of an undefined object");
00341 if(*count > 1)
00342 {
00343 (*count)--;
00344 count = new int;
00345 *count = 1;
00346 typename std::map<SPMatindex,T> **newv= new std::map<SPMatindex,T> * [nn];
00347 for(SPMatindex i=0; i<nn; ++i) if(v[i])
00348 newv[i]= new typename std::map<SPMatindex,T>(*(v[i]));
00349 else
00350 newv[i]= NULL;
00351 v = newv;
00352 }
00353 }
00354
00355
00356 template <typename T>
00357 void SparseSMat<T>::add(const SPMatindex n, const SPMatindex m, const T elem, const bool both)
00358 {
00359 #ifdef DEBUG
00360 if(n>=nn || m>=mm) laerror("illegal index in SparseSMat::add()");
00361 #endif
00362 if(!v[n]) v[n] = new std::map<SPMatindex,T>;
00363
00364 typename std::map<SPMatindex,T>::iterator p;
00365
00366 p= v[n]->find(m);
00367 if(p!=v[n]->end()) p->second+=elem; else (*v[n])[m] = elem;
00368 if(n!=m && both)
00369 {
00370 if(!v[m]) v[m] = new std::map<SPMatindex,T>;
00371 p= v[m]->find(n);
00372 if(p!=v[m]->end()) p->second+=elem; else (*v[m])[n] = elem;
00373 }
00374
00375 }
00376
00377
00378 template <typename T>
00379 unsigned long long SparseSMat<T>::simplify()
00380 {
00381 unsigned long long count=0;
00382 for(SPMatindex i=0; i<nn; ++i) if(v[i])
00383 {
00384
00385
00386 std::list<SPMatindex> l;
00387 typename std::map<SPMatindex,T>::iterator p;
00388 for(p=v[i]->begin(); p!=v[i]->end(); ++p)
00389 if(std::abs(p->second) < SPARSEEPSILON) l.push_front(p->first); else ++count;
00390 typename std::list<SPMatindex>::iterator q;
00391 for(q=l.begin(); q!=l.end(); ++q) v[i]->erase(*q);
00392 if(v[i]->size() == 0) {delete v[i]; v[i]=NULL;}
00393 }
00394 return count;
00395 }
00396
00397 template <typename T>
00398 std::ostream & operator<<(std::ostream &s, const SparseSMat<T> &x)
00399 {
00400 SPMatindex n;
00401
00402 s << x.nrows() << " "<< x.ncols()<< std::endl;
00403
00404 typename SparseSMat<T>::iterator p(x);
00405 for(; p.notend(); ++p) s << (int)p->row << ' ' << (int)p->col << ' ' << (typename LA_traits_io<T>::IOtype) p->elem << '\n';
00406 s << "-1 -1\n";
00407 return s;
00408 }
00409
00410 template <class T>
00411 std::istream& operator>>(std::istream &s, SparseSMat<T> &x)
00412 {
00413 SPMatindex n,m;
00414 long i,j;
00415 s >> n >> m;
00416 if(n!=m) laerror("SparseSMat must be square");
00417 x.resize(n,m);
00418 s >> i >> j;
00419 typename LA_traits_io<T>::IOtype tmp;
00420 while(i>=0 && j>=0)
00421 {
00422 s>>tmp;
00423 if(i>=n||j>=m) laerror("bad index in SparseSMat::operator>>");
00424 x.add(i,j,tmp,false);
00425 s >> i >> j;
00426 }
00427 return s;
00428 }
00429
00430
00431 template <typename T>
00432 SparseSMat<T> SparseSMat<T>::transpose(bool conj) const
00433 {
00434 SparseSMat<T> r(mm,nn);
00435 typename SparseSMat<T>::iterator p(*this);
00436 for(; p.notend(); ++p) r.add(p->col, p->row, (conj?LA_traits<T>::conjugate(p->elem):p->elem), false);
00437 return r;
00438 }
00439
00440
00441
00442
00443
00444
00445
00446
00447 template <typename T>
00448 SparseSMat<T> SparseSMat<T>::cholesky(void) const
00449 {
00450 if(nn!=mm) laerror("Cholesky defined only for square matrices");
00451
00452 NRVec<typename LA_traits<T>::normtype> diagreal(nn);
00453 {
00454 NRVec<T> diag(nn); diagonalof(diag);
00455 for(SPMatindex i=0; i<nn; ++i) diagreal[i]=LA_traits<T>::realpart(diag[i]);
00456 }
00457
00458 NRVec<int> pivot(nn);
00459 for(int i=0; i<nn; ++i) pivot[i]=i;
00460
00461
00463
00464
00465
00466 diagreal.sort(1,0,nn-1,pivot);
00467
00468
00469 NRVec<int> invpivot(nn);
00470 for(int i=0; i<nn; ++i) invpivot[pivot[i]]=i;
00471
00472
00473
00474
00475
00476 SparseSMat<T> r;
00477 r.nn=nn;
00478 r.mm=nn;
00479 r.count = new int(1);
00480 r.v = new std::map<SPMatindex,T> * [nn];
00481 for(SPMatindex i=0; i<nn; ++i)
00482 if(v[pivot[i]])
00483 {
00484 r.v[i]= new typename std::map<SPMatindex,T>;
00485 typename std::map<SPMatindex,T>::iterator p;
00486 for(p=v[pivot[i]]->begin(); p!=v[pivot[i]]->end(); ++p)
00487 {
00488 if(invpivot[p->first] >= i)
00489 {
00490 (*r.v[i])[invpivot[p->first]] = p->second;
00491 }
00492 }
00493 }
00494 else
00495 r.v[i]= NULL;
00496
00497
00498
00499
00500
00501 SPMatindex i,j,k;
00502 int rank=0;
00503 for(k=0; k<nn; ++k)
00504 if(r.v[k])
00505 {
00506 typename std::map<SPMatindex,T>::iterator p;
00507 p= r.v[k]->find(k);
00508 if(p==r.v[k]->end()) continue;
00509 if(LA_traits<T>::realpart(p->second) < CHOLESKYEPSILON) continue;
00510 ++rank;
00511 typename LA_traits<T>::normtype tmp = std::sqrt(LA_traits<T>::realpart(p->second));
00512 p->second = tmp;
00513 NRVec<T> linek(0.,nn);
00514 for(p=r.v[k]->begin(); p!=r.v[k]->end(); ++p)
00515 {
00516 if(p->first > k) p->second /= tmp;
00517 linek[p->first] = p->second;
00518 }
00519 for(j=k+1; j<nn; ++j)
00520 if(r.v[j])
00521 {
00522 T akj = LA_traits<T>::conjugate(linek[j]);
00523 NRVec<int> linek_used(0,nn);
00524 if(std::abs(akj)>SPARSEEPSILON)
00525 {
00526 for(p=r.v[j]->begin(); p!=r.v[j]->end(); ++p)
00527 {
00528 p->second -= akj * linek[p->first];
00529 linek_used[p->first]=1;
00530 }
00531
00532 for(i=j; i<nn; ++i)
00533 if(!linek_used[i] && std::abs(linek[i]) > SPARSEEPSILON)
00534 {
00535 (*r.v[j])[i] = -akj * linek[i];
00536 }
00537 }
00538 }
00539 }
00540
00541
00542
00543
00544
00545
00546
00547
00548
00549
00550
00551
00552
00553
00554
00555
00556 for(k=0; k<nn; ++k)
00557 if(r.v[k])
00558 {
00559 typename std::map<SPMatindex,T>::iterator p;
00560 typename std::map<SPMatindex,T> *vnew = new typename std::map<SPMatindex,T>;
00561 for(p=r.v[k]->begin(); p!=r.v[k]->end(); ++p)
00562 {
00563 if(std::abs(p->second) > SPARSEEPSILON) (*vnew)[pivot[p->first]] = p->second;
00564 }
00565 delete r.v[k];
00566 r.v[k]=vnew;
00567 }
00568
00569 return r;
00570 }
00571
00572
00573
00574
00575 template<typename T>
00576 SparseSMat<T> otimes_sparse(const NRVec<T> &lhs, const NRVec<T> &rhs, const bool conjugate=false, const T &scale=1)
00577 {
00578 SparseSMat<T> r(lhs.size(),rhs.size());
00579 for(SPMatindex i=0; i<lhs.size(); ++i)
00580 if(lhs[i]!=(T)0)
00581 {
00582 for(SPMatindex j=0; j<rhs.size(); ++j)
00583 if(rhs[j]!=(T)0)
00584 {
00585 T x=lhs[i]*(conjugate?LA_traits<T>::conjugate(rhs[j]):rhs[j])*scale;
00586 if(std::abs(x)>SPARSEEPSILON) r.add(i,j,x);
00587 }
00588 }
00589 return r;
00590 }
00591
00592
00593
00594
00595 }
00596 #endif //_SPARSESMAT_H_