00001
00002
00003
00004 #ifndef sparseMatrix_hh
00005 #define sparseMatrix_hh
00006
00007 #include <map>
00008 #if __GNUC__ == 2
00009 #include <float.h>
00010 #define EPS DBL_EPSILON
00011 #else
00012 #include <limits>
00013 #define EPS std::numeric_limits<double>::epsilon()
00014 #endif
00015 #include <memory>
00016 #include "compositions.hh"
00017 #include "matrix.hh"
00018 #include "CRS.hh"
00019 #include "hashedSMatrix.hh"
00020 #include "matrixMult.hh"
00021 #include "basics/debug.hh"
00022
00023 #define SparseAddInto_D 0
00024
00025 namespace concepts {
00026
00027
00028 template<class F, class G>
00029 class BilinearForm;
00030
00031
00032
00052 template<class F>
00053 class SparseMatrix : public Matrix<F>, public CRSConvertable<F> {
00054 public:
00056 typedef typename Realtype<F>::type r_type;
00058 typedef typename Cmplxtype<F>::type c_type;
00059
00060 typedef _HashedSMatrix_iterator<F, F&, F*> iterator;
00061 typedef _HashedSMatrix_iterator<F, const F&, const F*> const_iterator;
00062
00065 template<class G>
00066 SparseMatrix(const Space<G>& spcX, const Space<G>& spcY)
00067 : Matrix<F>(spcX.dim(), spcY.dim())
00068 , nX_(spcX.dim())
00069 , nY_(spcY.dim()),
00070 m_(new HashedSparseMatrix<F>(spcX.dim(), spcY.dim(), 2))
00071 { }
00072
00075 SparseMatrix(uint nofrows, uint nofcols)
00076 : Matrix<F>(nofrows, nofcols)
00077 , nX_(nofrows)
00078 , nY_(nofcols),
00079 m_(new HashedSparseMatrix<F>(nofrows, nofcols, 2))
00080 { }
00081 SparseMatrix(uint dim)
00082 : Matrix<F>(dim, dim)
00083 , nX_(dim)
00084 , nY_(dim),
00085 m_(new HashedSparseMatrix<F>(dim, dim, 2))
00086 { }
00087
00088
00114 template<class G>
00115 SparseMatrix(const Space<G>& spcX, const Space<G>& spcY,
00116 BilinearForm<F,G>& bf, const Real threshold = 0.0);
00117
00127 template<class G>
00128 SparseMatrix(const Space<G>& spc, BilinearForm<F,G>& bf,
00129 const Real threshold = 0.0);
00130
00134 SparseMatrix(const SparseMatrix<F>& m, bool t = false);
00135
00143 template<class G>
00144 SparseMatrix(const Space<G>& spc, const F* const v, const int i = -1);
00145
00152 template<class G>
00153 SparseMatrix(const Space<G>& spc,
00154 const Vector<F>& x, const Vector<F>& y);
00155
00165 template<class H>
00166 SparseMatrix(Operator<H>& A, bool slow = false);
00167
00168 void operator=(const SparseMatrix<F>&);
00169
00175 template<class H>
00176 SparseMatrix(const SparseMatrix<H>& fncX, F fnc(const H&));
00177 template<class H>
00178 SparseMatrix(const SparseMatrix<H>& fncX, const F& fnc(const H&));
00179
00185 template<class H>
00186 SparseMatrix(const SparseMatrix<H>& fncX);
00187
00188 virtual void operator()(const Function<r_type>& fncY,
00189 Function<F>& fncX);
00190 virtual void operator()(const Function<c_type>& fncY,
00191 Function<c_type>& fncX);
00192
00194 template<class H, class I>
00195 void operator()(const Vector<H>& fncY, Vector<I>& fncX) {
00196 (*m_)((const H*)fncY, (I*)fncX);
00197 }
00198
00205 const_iterator begin(uint row = 0) const;
00212 iterator begin(uint row = 0);
00214 const_iterator end() const;
00215
00219 virtual void transpMult(const Vector<r_type>& fncY,
00220 Vector<F>& fncX);
00221 virtual void transpMult(const Vector<c_type>& fncY,
00222 Vector<c_type>& fncX);
00223
00224 virtual F operator()(const uint i, const uint j) const {
00225 return (const_cast<const HashedSparseMatrix<F>&>(*m_))(i,j); }
00226 virtual F& operator()(const uint i, const uint j) { return (*m_)(i,j); }
00227
00228 SparseMatrix<F>& operator*=(const F factor) {
00229 (*m_) *= factor;
00230 return *this;
00231 }
00232
00236 void write(char* fname) const;
00237
00241 bool storeMatlab(const std::string filename, const std::string name = "",
00242 bool append = false) const;
00243
00245 const HashedSparseMatrix<F>* m() const {
00246 return const_cast<const HashedSparseMatrix<F>*>(m_.get()); }
00247
00261 template<class H, class I>
00262 void addInto(Matrix<H>& dest, const I fact,
00263 const uint rowoffset = 0, const uint coloffset = 0) const;
00272 template<class H, class I>
00273 void addIntoT(Matrix<H>& dest, const I fact,
00274 const uint rowoffset = 0, const uint coloffset = 0) const;
00275
00277 void multiply(const SparseMatrix<F>& fact, Matrix<F>& dest) const {
00278 m_->multiply(fact.m(), dest);
00279 }
00281 template<class H>
00282 void multiply(const H& fact, Matrix<F>& dest) const {
00283 matrixMultiplyRowSorting(*this, fact, dest);
00284 }
00286 virtual uint used() const { return m_->used(); }
00288 float memory() const { return sizeof(SparseMatrix<F>) + m_->memory(); }
00290 void symmetrize();
00291
00297 void compress(Real threshold = EPS) {
00298 m_->compress(threshold);
00299 }
00300
00302 void histogram(std::map<int, uint>& hist) const;
00303
00305 void copy(const SparseMatrix<F>& n);
00306
00307 virtual void convertCRS(F* a, int* asub, int* xa) const;
00308 virtual void convertCCS(F* a, int* asub, int* xa) const;
00309 protected:
00310 virtual std::ostream& info(std::ostream& os) const;
00311 private:
00313 uint nX_;
00314
00316 uint nY_;
00317
00319 std::auto_ptr<HashedSparseMatrix<F> > m_;
00320
00322 static uint storeMatlabCounter_;
00323
00324 void copyConstructor_(const SparseMatrix<F>& m, bool t);
00325
00326 virtual void apply_(const Vector<F>& fncY, Vector<F>& fncX) {
00327 conceptsAssert(0, Assertion());
00328 (*m_)((const F*)fncY, (F*)fncX);
00329 }
00330 };
00331
00332 template<class F>
00333 template<class H>
00334 SparseMatrix<F>::SparseMatrix(const SparseMatrix<H>& m, F fnc(const H&))
00335 : Matrix<F>(m.dimX(), m.dimY())
00336 , nX_(m.dimX()), nY_(m.dimY())
00337 , m_(new HashedSparseMatrix<F>(nX_, nY_, 2))
00338 {
00339 int hashBits = m.m()->HashBits();
00340 int p = 1 << hashBits;
00341 typename HashedSparseMatrix<H>::Value** matrix = m.m()->Matrix();
00342
00343 for (uint r = 0; r < nX_; ++r)
00344 for (int c_mod = 0; c_mod < p; ++c_mod)
00345 for (typename HashedSparseMatrix<H>::Value* v =
00346 matrix[(r<<hashBits) + c_mod]; v != NULL; v = v->lnk)
00347 (*m_)(v->idx, r) = fnc(v->val);
00348 }
00349
00350 template<class F>
00351 template<class H>
00352 SparseMatrix<F>::SparseMatrix(const SparseMatrix<H>& m)
00353 : Matrix<F>(m.dimX(), m.dimY())
00354 , nX_(m.dimX()), nY_(m.dimY())
00355 , m_(new HashedSparseMatrix<F>(nX_, nY_, 2))
00356 {
00357 int hashBits = m.m()->HashBits();
00358 int p = 1 << hashBits;
00359 typename HashedSparseMatrix<H>::Value** matrix = m.m()->Matrix();
00360
00361 for (uint r = 0; r < nX_; ++r)
00362 for (int c_mod = 0; c_mod < p; ++c_mod)
00363 for (typename HashedSparseMatrix<H>::Value* v =
00364 matrix[(r<<hashBits) + c_mod]; v != NULL; v = v->lnk)
00365 (*m_)(r, v->idx) = (F)(v->val);
00366 }
00367
00368 template<class F>
00369 template<class H>
00370 SparseMatrix<F>::SparseMatrix(const SparseMatrix<H>& m,
00371 const F& fnc(const H&))
00372 : Matrix<F>(m.dimX(), m.dimY())
00373 , nX_(m.dimX()), nY_(m.dimY())
00374 , m_(new HashedSparseMatrix<F>(nX_, nY_, 2))
00375 {
00376 int hashBits = m.m()->HashBits();
00377 int p = 1 << hashBits;
00378 typename HashedSparseMatrix<H>::Value** matrix = m.m()->Matrix();
00379
00380 for (uint r = 0; r < nX_; ++r)
00381 for (int c_mod = 0; c_mod < p; ++c_mod)
00382 for (typename HashedSparseMatrix<H>::Value* v =
00383 matrix[(r<<hashBits) + c_mod]; v != NULL; v = v->lnk)
00384 (*m_)(v->idx, r) = fnc(v->val);
00385 }
00386
00387
00388 template<class F>
00389 template<class H, class I>
00390 void SparseMatrix<F>::addInto(Matrix<H>& dest, const I fact,
00391 const uint rowoffset,
00392 const uint coloffset) const
00393 {
00394 conceptsAssert(dest.dimX() >= this->dimX() + rowoffset, Assertion());
00395 conceptsAssert3(dest.dimY() >= this->dimY() + coloffset, Assertion(),
00396 "dest.dimY() = " << dest.dimY() << ", dimY() = " <<
00397 this->dimY() << ", coloffset = " << coloffset);
00398 DEBUGL(SparseAddInto_D, "fact = " << fact <<
00399 ", rowoffset = " << rowoffset <<
00400 ", coloffset = " << coloffset);
00401 if (fact != 0.0) {
00402
00403 int hashBits = m_->HashBits();
00404 int p = 1 << hashBits;
00405 uint n = this->dimX();
00406 typename HashedSparseMatrix<F>::Value** matrix = m_->Matrix();
00407
00408 for (uint r = 0; r < n; ++r)
00409 for (int c_mod = 0; c_mod < p; ++c_mod)
00410 for (typename HashedSparseMatrix<F>::Value* v =
00411 matrix[(r<<hashBits) + c_mod]; v != NULL; v = v->lnk)
00412 if (v->val != 0.0)
00413 dest(r + rowoffset, v->idx + coloffset) += (v->val * fact);
00414 }
00415 }
00416
00417 template<class F>
00418 template<class H, class I>
00419 void SparseMatrix<F>::addIntoT(Matrix<H>& dest, const I fact,
00420 const uint rowoffset,
00421 const uint coloffset) const
00422 {
00423 conceptsAssert3(dest.dimY() >= this->dimX() + coloffset, Assertion(),
00424 "dest.dimY() = " << dest.dimX() << ", dimX() = " <<
00425 this->dimX() << ", coloffset = " << coloffset);
00426 conceptsAssert3(dest.dimX() >= this->dimY() + rowoffset, Assertion(),
00427 "dest.dimX() = " << dest.dimX() << ", dimY() = " <<
00428 this->dimY() << ", rowoffset = " << rowoffset);
00429 DEBUGL(SparseAddInto_D, "fact = " << fact <<
00430 ", rowoffset = " << rowoffset <<
00431 ", coloffset = " << coloffset);
00432 if (fact != 0.0) {
00433 int hashBits = m_->HashBits();
00434 int p = 1 << hashBits;
00435 uint n = this->dimX();
00436 typename HashedSparseMatrix<F>::Value** matrix = m_->Matrix();
00437
00438 for (uint r = 0; r < n; ++r)
00439 for (int c_mod = 0; c_mod < p; ++c_mod)
00440 for (typename HashedSparseMatrix<F>::Value* v =
00441 matrix[(r<<hashBits) + c_mod]; v != NULL; v = v->lnk)
00442 if (v->val != 0.0)
00443 dest(v->idx + rowoffset, r + coloffset) += (v->val*fact);
00444 }
00445 }
00446
00447 }
00448
00449 #endif // sparseMatrix_hh