00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018 #ifndef _SPARSEMAT_H_
00019 #define _SPARSEMAT_H_
00020 #include "la_traits.h"
00021
00022 namespace LA {
00023
00024
00025
00026 const double SPARSEEPSILON=1e-35;
00027
00028 typedef unsigned int SPMatindex;
00029 typedef int SPMatindexdiff;
00030
00031
00032 template<typename T>
00033 struct matel
00034 {
00035 T elem;
00036 SPMatindex row;
00037 SPMatindex col;
00038 matel *next;
00039 };
00040
00041
00042 template <typename T>
00043 class SparseMat {
00044 protected:
00045 SPMatindex nn;
00046 SPMatindex mm;
00047 bool symmetric;
00048 unsigned int nonzero;
00049 int *count;
00050 matel<T> *list;
00051 private:
00052 matel<T> **rowsorted;
00053 matel<T> **colsorted;
00054 void unsort();
00055 void deletelist();
00056 void copylist(const matel<T> *l);
00057 public:
00058
00059 class iterator {
00060 private:
00061 matel<T> *p;
00062 public:
00063 iterator() {};
00064 ~iterator() {};
00065 iterator(matel<T> *list): p(list) {};
00066 bool operator==(const iterator rhs) const {return p==rhs.p;}
00067 bool operator!=(const iterator rhs) const {return p!=rhs.p;}
00068 iterator operator++() {return p=p->next;}
00069 iterator operator++(int) {matel<T> *q=p; p=p->next; return q;}
00070 matel<T> & operator*() const {return *p;}
00071 matel<T> * operator->() const {return p;}
00072 };
00073 iterator begin() const {return list;}
00074 iterator end() const {return NULL;}
00075
00076
00077 inline SparseMat() :nn(0),mm(0),symmetric(0),nonzero(0),count(NULL),list(NULL),rowsorted(NULL),colsorted(NULL) {};
00078 inline SparseMat(const SPMatindex n, const SPMatindex m) :nn(n),mm(m),symmetric(0),nonzero(0),count(new int(1)),list(NULL),rowsorted(NULL),colsorted(NULL) {};
00079 SparseMat(const SparseMat &rhs);
00080 inline int getcount() const {return count?*count:0;}
00081 explicit SparseMat(const NRMat<T> &rhs);
00082 explicit SparseMat(const NRSMat<T> &rhs);
00083 SparseMat & operator=(const SparseMat &rhs);
00084 SparseMat & operator=(const T &a);
00085 SparseMat & operator+=(const T &a);
00086 SparseMat & operator-=(const T &a);
00087 SparseMat & operator*=(const T &a);
00088 SparseMat & operator+=(const SparseMat &rhs);
00089 SparseMat & addtriangle(const SparseMat &rhs, const bool lower, const char sign);
00090 SparseMat & join(SparseMat &rhs);
00091 SparseMat & operator-=(const SparseMat &rhs);
00092 inline const SparseMat operator+(const T &rhs) const {return SparseMat(*this) += rhs;}
00093 inline const SparseMat operator-(const T &rhs) const {return SparseMat(*this) -= rhs;}
00094 inline const SparseMat operator*(const T &rhs) const {return SparseMat(*this) *= rhs;}
00095 inline const SparseMat operator+(const SparseMat &rhs) const {return SparseMat(*this) += rhs;}
00096 inline const SparseMat operator-(const SparseMat &rhs) const {return SparseMat(*this) -= rhs;}
00097 inline const NRVec<T> operator*(const NRVec<T> &rhs) const;
00098 inline const NRMat<T> operator*(const NRMat<T> &rhs) const;
00099 const T* diagonalof(NRVec<T> &, const bool divide=0, bool cache=false) const;
00100 void gemv(const T beta, NRVec<T> &r, const char trans, const T alpha, const NRVec<T> &x, bool treat_as_symmetric=false) const {r.gemv(beta,*this,trans,alpha,x,treat_as_symmetric);};
00101 const SparseMat operator*(const SparseMat &rhs) const;
00102 SparseMat & oplusequal(const SparseMat &rhs);
00103 SparseMat & oplusequal(const NRMat<T> &rhs);
00104 SparseMat & oplusequal(const NRSMat<T> &rhs);
00105 const SparseMat oplus(const SparseMat &rhs) const {return SparseMat(*this).oplusequal(rhs);};
00106 const SparseMat oplus(const NRMat<T> &rhs) const {return SparseMat(*this).oplusequal(rhs);};
00107 const SparseMat oplus(const NRSMat<T> &rhs) const {return SparseMat(*this).oplusequal(rhs);};
00108 const SparseMat otimes(const SparseMat &rhs) const;
00109 const SparseMat otimes(const NRMat<T> &rhs) const;
00110 const SparseMat otimes(const NRSMat<T> &rhs) const;
00111 void gemm(const T beta, const SparseMat &a, const char transa, const SparseMat &b, const char transb, const T alpha);
00112 const T dot(const SparseMat &rhs) const;
00113 const T dot(const NRMat<T> &rhs) const;
00114 inline ~SparseMat();
00115 void axpy(const T alpha, const SparseMat &x, const bool transp=0);
00116 inline matel<T> *getlist() const {return list;}
00117 void setlist(matel<T> *l) {list=l;}
00118 inline SPMatindex nrows() const {return nn;}
00119 inline SPMatindex ncols() const {return mm;}
00120 void get(int fd, bool dimensions=1, bool transposed=false);
00121 void put(int fd, bool dimensions=1, bool transposed=false) const;
00122 void resize(const SPMatindex n, const SPMatindex m);
00123 void incsize(const SPMatindex n, const SPMatindex m);
00124 void transposeme();
00125 const SparseMat transpose() const;
00126 void permuteindices(const NRVec<SPMatindex> &p);
00127 void permuterows(const NRVec<SPMatindex> &p);
00128 void permutecolumns(const NRVec<SPMatindex> &p);
00129 inline void setsymmetric() {if(nn!=mm) laerror("non-square cannot be symmetric"); symmetric=1;}
00130 inline void defineunsymmetric() {symmetric=0;}
00131 void setunsymmetric();
00132 inline bool issymmetric() const {return symmetric;}
00133 unsigned int length() const;
00134 void copyonwrite();
00135 void clear() {copyonwrite();unsort();deletelist();}
00136 void simplify();
00137 const T trace() const;
00138 const typename LA_traits<T>::normtype norm(const T scalar=(T)0) const;
00139 unsigned int sort(int type) const;
00140 inline void add(const SPMatindex n, const SPMatindex m, const T elem) {matel<T> *ltmp= new matel<T>; ltmp->next=list; list=ltmp; list->row=n; list->col=m; list->elem=elem;}
00141 void addsafe(const SPMatindex n, const SPMatindex m, const T elem);
00142 };
00143
00144 }
00145
00146
00147 #include "vec.h"
00148 #include "smat.h"
00149 #include "mat.h"
00150
00151 namespace LA {
00152
00153 template <typename T>
00154 inline const NRVec<T> SparseMat<T>::operator*(const NRVec<T> &rhs) const
00155 {NRVec<T> result(nn); result.gemv((T)0,*this,'n',(T)1,rhs); return result;};
00156
00157 template <typename T>
00158 inline const NRMat<T> SparseMat<T>::operator*(const NRMat<T> &rhs) const
00159 {NRMat<T> result(nn,rhs.ncols()); result.gemm((T)0,*this,'n',rhs,'n',(T)1); return result;};
00160
00161 template <class T>
00162 std::ostream& operator<<(std::ostream &s, const SparseMat<T> &x)
00163 {
00164 SPMatindex n,m;
00165 n=x.nrows();
00166 m=x.ncols();
00167 s << (int)n << ' ' << (int)m << '\n';
00168 matel<T> *list=x.getlist();
00169 while(list)
00170 {
00171 s << (int)list->row << ' ' << (int)list->col << ' ' << (typename LA_traits_io<T>::IOtype)list->elem << '\n';
00172 list=list->next;
00173 }
00174 s << "-1 -1\n";
00175 return s;
00176 }
00177
00178 template <class T>
00179 std::istream& operator>>(std::istream &s, SparseMat<T> &x)
00180 {
00181 int i,j;
00182 int n,m;
00183 matel<T> *l=NULL;
00184 s >> n >> m;
00185 x.resize(n,m);
00186 s >> i >> j;
00187 while(i>=0 && j>=0)
00188 {
00189 matel<T> *ll = l;
00190 l= new matel<T>;
00191 l->next= ll;
00192 l->row=i;
00193 l->col=j;
00194 typename LA_traits_io<T>::IOtype tmp;
00195 s >> tmp;
00196 l->elem=tmp;
00197 s >> i >> j;
00198 }
00199 x.setlist(l);
00200 return s;
00201 }
00202
00203
00204
00205 template <typename T>
00206 SparseMat<T>::~SparseMat()
00207 {
00208 unsort();
00209 if(!count) return;
00210 if(--(*count)<=0)
00211 {
00212 deletelist();
00213 delete count;
00214 }
00215 }
00216
00217
00218 template <typename T>
00219 SparseMat<T>::SparseMat(const SparseMat<T> &rhs)
00220 {
00221 #ifdef debug
00222 if(! &rhs) laerror("SparseMat copy constructor with NULL argument");
00223 #endif
00224 nn=rhs.nn;
00225 mm=rhs.mm;
00226 symmetric=rhs.symmetric;
00227 if(rhs.list&&!rhs.count) laerror("some inconsistency in SparseMat contructors or assignments");
00228 list=rhs.list;
00229 if(list) {count=rhs.count; (*count)++;} else count=new int(1);
00230 colsorted=rowsorted=NULL;
00231 nonzero=0;
00232 }
00233
00234 template <typename T>
00235 const SparseMat<T> SparseMat<T>::transpose() const
00236 {
00237 if(list&&!count) laerror("some inconsistency in SparseMat transpose");
00238 SparseMat<T> result;
00239 result.nn=mm;
00240 result.mm=nn;
00241 result.symmetric=symmetric;
00242 if(result.symmetric)
00243 {
00244 result.list=list;
00245 if(list) {result.count=count; (*result.count)++;} else result.count=new int(1);
00246 }
00247 else
00248 {
00249 result.count=new int(1);
00250 result.list=NULL;
00251 matel<T> *l =list;
00252 while(l)
00253 {
00254 result.add(l->col,l->row,l->elem);
00255 l=l->next;
00256 }
00257 }
00258 result.colsorted=result.rowsorted=NULL;
00259 result.nonzero=0;
00260 return result;
00261 }
00262
00263
00264
00265 template<typename T>
00266 inline const SparseMat<T> commutator ( const SparseMat<T> &x, const SparseMat<T> &y, const bool trx=0, const bool tryy=0)
00267 {
00268 SparseMat<T> r;
00269 r.gemm((T)0,x,trx?'t':'n',y,tryy?'t':'n',(T)1);
00270 r.gemm((T)1,y,tryy?'t':'n',x,trx?'t':'n',(T)-1);
00271 return r;
00272 }
00273
00274 template<typename T>
00275 inline const SparseMat<T> anticommutator ( const SparseMat<T> &x, const SparseMat<T> &y, const bool trx=0, const bool tryy=0)
00276 {
00277 SparseMat<T> r;
00278 r.gemm((T)0,x,trx?'t':'n',y,tryy?'t':'n',(T)1);
00279 r.gemm((T)1,y,tryy?'t':'n',x,trx?'t':'n',(T)1);
00280 return r;
00281 }
00282
00283
00284
00285 template<typename T>
00286 NRMat<T> & NRMat<T>::operator+=(const SparseMat<T> &rhs)
00287 {
00288 if((unsigned int)nn!=rhs.nrows()||(unsigned int)mm!=rhs.ncols()) laerror("incompatible matrices in +=");
00289 matel<T> *l=rhs.getlist();
00290 bool sym=rhs.issymmetric();
00291 while(l)
00292 {
00293 #ifdef MATPTR
00294 v[l->row][l->col] +=l->elem;
00295 if(sym && l->row!=l->col) v[l->col][l->row] +=l->elem;
00296 #else
00297 v[l->row*mm+l->col] +=l->elem;
00298 if(sym && l->row!=l->col) v[l->col*mm+l->row] +=l->elem;
00299 #endif
00300 l=l->next;
00301 }
00302 return *this;
00303 }
00304
00305 }
00306 #endif