sparseMatrix.hh

Go to the documentation of this file.
1 /* @file operator/sparseMatrix.hh
2 
3  A sparse matrix
4  */
5 
6 #ifndef sparseMatrix_hh
7 #define sparseMatrix_hh
8 
9 #include <map>
10 #if __GNUC__ == 2
11 #include <float.h>
12 #define EPS DBL_EPSILON
13 #else
14 #include <limits>
15 #define EPS std::numeric_limits<double>::epsilon()
16 #endif
17 #include <memory>
18 #include "compositions.hh"
19 #include "matrix.hh"
20 #include "CRS.hh"
21 #include "hashedSMatrix.hh"
22 #include "matrixMult.hh"
23 #include "basics/debug.hh"
24 #include "toolbox/sequence.hh"
25 
26 #include "bilinearForm.hh"
27 
28 #include <type_traits>
29 #include <typeinfo>
30 
31 #define SparseAddInto_D 0
32 
33 namespace concepts {
34 
35  // forward declaration
36  template<class F, class G>
37  class BilinearForm;
38 
39  // forward declaration
40  //template<class F, class G>
41  //class BilinearFormContainer;
42 
43  // ********************************************************** SparseMatrix **
44 
64  template<class F>
65  class SparseMatrix : public Matrix<F>, public CRSConvertable<F> {
66  public:
68  typedef typename Realtype<F>::type r_type;
70  typedef typename Cmplxtype<F>::type c_type;
72  typedef typename std::conditional<std::is_same<typename Realtype<F>::type, F>::value ,
73  typename Realtype<F>::type, typename Cmplxtype<F>::type >::type d_type;
74 
77 
80  template<class G>
81  SparseMatrix(const Space<G>& spcX, const Space<G>& spcY)
82  : Matrix<F>(spcX.dim(), spcY.dim())
83  , nX_(spcX.dim())
84  , nY_(spcY.dim()),
85  m_(new HashedSparseMatrix<F>(spcX.dim(), spcY.dim(), 2))
86  { }
87 
90  SparseMatrix(uint nofrows, uint nofcols)
91  : Matrix<F>(nofrows, nofcols)
92  , nX_(nofrows)
93  , nY_(nofcols),
94  m_(new HashedSparseMatrix<F>(nofrows, nofcols, 2))
95  { }
96  SparseMatrix(uint dim = 0)
97  : Matrix<F>(dim, dim)
98  , nX_(dim)
99  , nY_(dim),
100  m_(new HashedSparseMatrix<F>(dim, dim, 2))
101  { }
102 
128  template<class G>
129  SparseMatrix(const Space<G>& spcX, const Space<G>& spcY,
130  const BilinearForm<F,G>& bf, const Real eps = 0.0,
131  const bool single = false);
132 
160  template<class G>
161  SparseMatrix(const Space<G>& spcX, const Space<G>& spcY,
162  const Sequence<bool> & seq,
163  const BilinearForm<F,G>& bf, const Real eps = 0.0,
164  const bool single = false);
165 
192  template<class G>
193  SparseMatrix(const Space<G>& spcX, const Space<G>& spcY,
194  const BilinearForm<F,G>& bf,
195  const Sequence<ElementWithCell<G> * > & seq,
196  const Real eps = 0.0);
197 
208  template<class G>
209  SparseMatrix(const Space<G>& spc, const BilinearForm<F,G>& bf,
210  const Real eps = 0.0);
211 
212 
223  const Real eps = 0.0);
224 
235  template<class G>
236  SparseMatrix(const Space<G>& spc,
237  const Sequence<bool> & seq, const BilinearForm<F,G>& bf,
238  const Real eps = 0.0);
239 
250  template<class G>
251  SparseMatrix(const Space<G>& spc, const BilinearForm<F,G>& bf,
252  const Sequence<ElementWithCell<G> * > & seq,
253  const Real eps = 0.0);
254 
258  SparseMatrix(const SparseMatrix<F>& m, bool t = false);
259 
267  template<class G>
268  SparseMatrix(const Space<G>& spc, const F* const v, const int i = -1);
269 
276  template<class G>
277  SparseMatrix(const Space<G>& spc,
278  const Vector<F>& x, const Vector<F>& y);
279 
289  template<class H>
290  SparseMatrix(Operator<H>& A, bool slow = false);
291 
293 
299  template<class H>
300  SparseMatrix(const SparseMatrix<H>& fncX, F fnc(const H&));
301  template<class H>
302  SparseMatrix(const SparseMatrix<H>& fncX, const F& fnc(const H&));
303 
309  template<class H>
311 
312 
313  virtual ~SparseMatrix(){};
314 
315  virtual void operator()(const Function<r_type>& fncY,
316  Function<F>& fncX);
317  virtual void operator()(const Function<c_type>& fncY,
318  Function<c_type>& fncX);
319 
321  template<class H, class I>
322  void operator()(const Vector<H>& fncY, Vector<I>& fncX) {
323  (*m_)((const H*)fncY, (I*)fncX);
324  }
325 
332  const_iterator begin(uint row = 0) const;
339  iterator begin(uint row = 0);
342 
346  virtual void transpMult(const Vector<r_type>& fncY,
347  Vector<F>& fncX);
348  virtual void transpMult(const Vector<c_type>& fncY,
349  Vector<c_type>& fncX);
350 
351  virtual F operator()(const uint i, const uint j) const {
352  return (const_cast<const HashedSparseMatrix<F>&>(*m_))(i,j); }
353  virtual F& operator()(const uint i, const uint j) { return (*m_)(i,j); }
354 
355  SparseMatrix<F>& operator*=(const F factor) {
356  (*m_) *= factor;
357  return *this;
358  }
359 
361  virtual void resize(uint m, uint n) {
362  this->dimX_ = nX_ = m; this->dimY_ = nY_ = n;
363  m_.reset(new HashedSparseMatrix<F>(m, n, 2));
364  }
365 
369  void write(char* fname) const;
370 
374  bool storeMatlab(const std::string filename, const std::string name = "",
375  bool append = false) const;
376 
378  const HashedSparseMatrix<F>* m() const {
379  return const_cast<const HashedSparseMatrix<F>*>(m_.get()); }
380 
394  template<class H, class I>
395  void addInto(Matrix<H>& dest, const I fact,
396  const uint rowoffset = 0, const uint coloffset = 0) const;
405  template<class H, class I>
406  void addIntoT(Matrix<H>& dest, const I fact,
407  const uint rowoffset = 0, const uint coloffset = 0) const;
408 
412  template<class H, class I>
413  void add(const Vector<H>& v, const I fact,
414  const uint rowoffset = 0, const uint coloffset = 0);
415 
419  template<class H, class I>
420  void addT(const Vector<H>& v, const I fact,
421  const uint rowoffset = 0, const uint coloffset = 0);
422 
424  void multiply(const SparseMatrix<F>& fact, Matrix<F>& dest) const {
425  m_->multiply(fact.m(), dest);
426  }
428  template<class H>
429  void multiply(const H& fact, Matrix<F>& dest) const {
430  matrixMultiplyRowSorting(*this, fact, dest);
431  }
433  virtual uint used() const { return m_->used(); }
435  float memory() const { return sizeof(SparseMatrix<F>) + m_->memory(); }
437  void symmetrize();
438 
444  void compress(Real threshold = EPS) {
445  m_->compress(threshold);
446  }
447 
449  void histogram(std::map<int, uint>& hist) const;
450 
452  void copy(const SparseMatrix<F>& n);
453 
454  virtual void convertCRS(F* a, int* asub, int* xa) const;
455  virtual void convertCCS(F* a, int* asub, int* xa) const;
456  virtual void convertIJK(F* , int* , int* ) const;
457  protected:
458  virtual std::ostream& info(std::ostream& os) const;
459  private:
461  uint nX_;
462 
464  uint nY_;
465 
467  std::unique_ptr<HashedSparseMatrix<F> > m_;
468 
470  static uint storeMatlabCounter_;
471 
472  void copyConstructor_(const SparseMatrix<F>& m, bool t);
473 
474  virtual void apply_(const Vector<F>& fncY, Vector<F>& fncX) {
476  (*m_)((const F*)fncY, (F*)fncX);
477  }
478  };
479 
480  template<class F>
481  template<class H>
482  SparseMatrix<F>::SparseMatrix(const SparseMatrix<H>& m, F fnc(const H&))
483  : Matrix<F>(m.dimX(), m.dimY())
484  , nX_(m.dimX()), nY_(m.dimY())
485  , m_(new HashedSparseMatrix<F>(nX_, nY_, 2))
486  {
487  int hashBits = m.m()->HashBits();
488  int p = 1 << hashBits;
489  typename HashedSparseMatrix<H>::Value** matrix = m.m()->Matrix();
490 
491  for (uint r = 0; r < nX_; ++r)
492  for (int c_mod = 0; c_mod < p; ++c_mod)
493  for (typename HashedSparseMatrix<H>::Value* v =
494  matrix[(r<<hashBits) + c_mod]; v != NULL; v = v->lnk)
495  (*m_)(r, v->idx) = fnc(v->val);
496  }
497 
498  template<class F>
499  template<class H>
501  : Matrix<F>(m.dimX(), m.dimY())
502  , nX_(m.dimX()), nY_(m.dimY())
503  , m_(new HashedSparseMatrix<F>(nX_, nY_, 2))
504  {
505  int hashBits = m.m()->HashBits();
506  int p = 1 << hashBits;
507  typename HashedSparseMatrix<H>::Value** matrix = m.m()->Matrix();
508 
509  for (uint r = 0; r < nX_; ++r)
510  for (int c_mod = 0; c_mod < p; ++c_mod)
511  for (typename HashedSparseMatrix<H>::Value* v =
512  matrix[(r<<hashBits) + c_mod]; v != NULL; v = v->lnk)
513  (*m_)(r, v->idx) = (F)(v->val);
514  }
515 
516  template<class F>
517  template<class H>
519  const F& fnc(const H&))
520  : Matrix<F>(m.dimX(), m.dimY())
521  , nX_(m.dimX()), nY_(m.dimY())
522  , m_(new HashedSparseMatrix<F>(nX_, nY_, 2))
523  {
524  int hashBits = m.m()->HashBits();
525  int p = 1 << hashBits;
526  typename HashedSparseMatrix<H>::Value** matrix = m.m()->Matrix();
527 
528  for (uint r = 0; r < nX_; ++r)
529  for (int c_mod = 0; c_mod < p; ++c_mod)
530  for (typename HashedSparseMatrix<H>::Value* v =
531  matrix[(r<<hashBits) + c_mod]; v != NULL; v = v->lnk)
532  (*m_)(r, v->idx) = fnc(v->val);
533  }
534 
535 
536  template<class F>
537  template<class H, class I>
538  void SparseMatrix<F>::addInto(Matrix<H>& dest, const I fact,
539  const uint rowoffset,
540  const uint coloffset) const
541  {
542  conceptsAssert(dest.dimX() >= this->dimX() + rowoffset, Assertion());
543  conceptsAssert3(dest.dimY() >= this->dimY() + coloffset, Assertion(),
544  "dest.dimY() = " << dest.dimY() << ", dimY() = " <<
545  this->dimY() << ", coloffset = " << coloffset);
546  DEBUGL(SparseAddInto_D, "fact = " << fact <<
547  ", rowoffset = " << rowoffset <<
548  ", coloffset = " << coloffset);
549  if (fact != 0.0) {
550 
551  int hashBits = m_->HashBits();
552  int p = 1 << hashBits;
553  uint n = this->dimX();
554  typename HashedSparseMatrix<F>::Value** matrix = m_->Matrix();
555 
556  for (uint r = 0; r < n; ++r)
557  for (int c_mod = 0; c_mod < p; ++c_mod)
558  for (typename HashedSparseMatrix<F>::Value* v =
559  matrix[(r<<hashBits) + c_mod]; v != NULL; v = v->lnk)
560  if (v->val != 0.0)
561  dest(r + rowoffset, v->idx + coloffset) += (v->val * fact);
562  }
563  }
564 
565  template<class F>
566  template<class H, class I>
567  void SparseMatrix<F>::addIntoT(Matrix<H>& dest, const I fact,
568  const uint rowoffset,
569  const uint coloffset) const
570  {
571  conceptsAssert3(dest.dimY() >= this->dimX() + coloffset, Assertion(),
572  "dest.dimY() = " << dest.dimX() << ", dimX() = " <<
573  this->dimX() << ", coloffset = " << coloffset);
574  conceptsAssert3(dest.dimX() >= this->dimY() + rowoffset, Assertion(),
575  "dest.dimX() = " << dest.dimX() << ", dimY() = " <<
576  this->dimY() << ", rowoffset = " << rowoffset);
577  DEBUGL(SparseAddInto_D, "fact = " << fact <<
578  ", rowoffset = " << rowoffset <<
579  ", coloffset = " << coloffset);
580  if (fact != 0.0) {
581  int hashBits = m_->HashBits();
582  int p = 1 << hashBits;
583  uint n = this->dimX();
584  typename HashedSparseMatrix<F>::Value** matrix = m_->Matrix();
585 
586  for (uint r = 0; r < n; ++r)
587  for (int c_mod = 0; c_mod < p; ++c_mod)
588  for (typename HashedSparseMatrix<F>::Value* v =
589  matrix[(r<<hashBits) + c_mod]; v != NULL; v = v->lnk)
590  if (v->val != 0.0)
591  dest(v->idx + rowoffset, r + coloffset) += (v->val*fact);
592  }
593  }
594 
595  template<class F>
596  template<class H, class I>
597  void SparseMatrix<F>::add(const Vector<H>& v, const I fact,
598  const uint rowoffset, const uint coloffset)
599  {
600  conceptsAssert(this->dimX() >= v.dim() + rowoffset, Assertion());
601  conceptsAssert(this->dimY() >= 1 + coloffset, Assertion());
602 
603  if (fact != 0.0) {
604  const H* h = (const H*)v;
605  uint row = rowoffset;
606  for(uint i = v.dim(); i--; )
607  (*this)(row++, coloffset) += *h++ * fact;
608  }
609  }
610 
611  template<class F>
612  template<class H, class I>
613  void SparseMatrix<F>::addT(const Vector<H>& v, const I fact,
614  const uint rowoffset, const uint coloffset)
615  {
616  conceptsAssert(this->dimX() >= 1 + rowoffset, Assertion());
617  conceptsAssert(this->dimY() >= v.dim() + coloffset, Assertion());
618 
619  if (fact != 0.0) {
620  const H* h = (const H*)v;
621  uint col = coloffset;
622  for(uint i = v.dim(); i--; )
623  (*this)(rowoffset, col++) += *h++ * fact;
624  }
625  }
626 
627 } // namespace concepts
628 
629 #endif // sparseMatrix_hh
SparseMatrix(const Space< G > &spc, const Vector< F > &x, const Vector< F > &y)
Constructor of rank 1 matrix.
std::unique_ptr< HashedSparseMatrix< F > > m_
The matrix.
SparseMatrix(uint nofrows, uint nofcols)
Constructor.
Definition: sparseMatrix.hh:90
SparseMatrix(uint dim=0)
Definition: sparseMatrix.hh:96
void multiply(const H &fact, Matrix< F > &dest) const
Multiplies this matrix with fact and adds the result to dest.
uint nX_
Dimension of image space (spcX_)
bool storeMatlab(const std::string filename, const std::string name="", bool append=false) const
Stores the matrix in a Matlab sparse matrix.
virtual void transpMult(const Vector< c_type > &fncY, Vector< c_type > &fncX)
std::conditional< std::is_same< typename Realtype< F >::type, F >::value, typename Realtype< F >::type, typename Cmplxtype< F >::type >::type d_type
Data type, depending if F is real or complex.
Definition: sparseMatrix.hh:73
SparseMatrix(const Space< G > &spc, const Sequence< bool > &seq, const BilinearForm< F, G > &bf, const Real eps=0.0)
Constructor.
const HashedSparseMatrix< F > * m() const
Returns the sparse matrix itself.
Abstract class for a function.
Definition: basis.hh:21
virtual void convertCCS(F *a, int *asub, int *xa) const
void operator=(const SparseMatrix< F > &)
SparseMatrix(const Space< G > &spcX, const Space< G > &spcY)
Constructor.
Definition: sparseMatrix.hh:81
void add(const Vector< H > &v, const I fact, const uint rowoffset=0, const uint coloffset=0)
Copies the vector v multiplied by fact on position (rowoffset, coloffset)
void compress(Real threshold=EPS)
Compresses the matrix by dropping small entries.
#define conceptsAssert(cond, exc)
Assert that a certain condition is fulfilled.
Definition: exceptions.hh:394
void histogram(std::map< int, uint > &hist) const
Creates a histogram of the matrix entries.
void copyConstructor_(const SparseMatrix< F > &m, bool t)
void addInto(Matrix< H > &dest, const I fact, const uint rowoffset=0, const uint coloffset=0) const
This matrix is added as block to the given matrix dest.
Element with cell.
#define DEBUGL(doit, msg)
Abstract function class to evaluate a bilinear form.
Definition: bilinearForm.hh:33
void symmetrize()
Makes sure a theoretically symmetric matrix is symmetric in memory too.
#define EPS
Definition: sparseMatrix.hh:15
virtual void resize(uint m, uint n)
Sets a new size, previous data might be lost
virtual uint used() const
Returns the number of used entries in the matrix.
void operator()(const Vector< H > &fncY, Vector< I > &fncX)
Multiplies the matrix with fncY. The result is fncX.
_HashedSMatrix_iterator< F, F &, F * > iterator
Definition: sparseMatrix.hh:75
Cmplxtype< F >::type c_type
Complex type of data type.
Definition: sparseMatrix.hh:70
uint nY_
Dimension of source space (spcY_)
virtual void apply_(const Vector< F > &fncY, Vector< F > &fncX)
Exception class for assertions.
Definition: exceptions.hh:258
SparseMatrix(const SparseMatrix< H > &fncX)
Constructor.
Realtype< F >::type r_type
Real type of data type.
Definition: sparseMatrix.hh:68
SparseMatrix(const Space< G > &spcX, const Space< G > &spcY, const Sequence< bool > &seq, const BilinearForm< F, G > &bf, const Real eps=0.0, const bool single=false)
Constructor.
Abstract class for an operator.
Definition: compositions.hh:31
#define SparseAddInto_D
Definition: sparseMatrix.hh:31
void copy(const SparseMatrix< F > &n)
Copies n to this matrix.
SparseMatrix(const SparseMatrix< H > &fncX, F fnc(const H &))
Constructor.
_HashedSMatrix_iterator< F, const F &, const F * > const_iterator
Definition: sparseMatrix.hh:76
SparseMatrix(const SparseMatrix< H > &fncX, const F &fnc(const H &))
virtual void convertIJK(F *, int *, int *) const
virtual void operator()(const Function< r_type > &fncY, Function< F > &fncX)
Abstract class for an operator.
Definition: ARPACK.hh:16
void addT(const Vector< H > &v, const I fact, const uint rowoffset=0, const uint coloffset=0)
Copies the transpose of the vector v multiplied by fact on position (rowoffset, coloffset)
Sequence with operations, output operator, and method of the particular element types.
Definition: sequence.hh:39
SparseMatrix(const Space< G > &spc, const BilinearForm< F, G > &bf, const Sequence< ElementWithCell< G > * > &seq, const Real eps=0.0)
Constructor.
virtual void convertCRS(F *a, int *asub, int *xa) const
static uint storeMatlabCounter_
Counts number of Matlab outputs (used to uniquely name the matrices)
iterator begin(uint row=0)
Iterator over the elements, standing at position (row,c), where row is the given row number and c the...
SparseMatrix(const SparseMatrix< F > &m, bool t=false)
Copy constructor.
SparseMatrix< F > & operator*=(const F factor)
virtual void operator()(const Function< c_type > &fncY, Function< c_type > &fncX)
A matrix in sparse storage using hashes.
Definition: compositions.hh:34
void multiply(const SparseMatrix< F > &fact, Matrix< F > &dest) const
Multiplies this matrix with fact and adds the result to dest.
SparseMatrix(Operator< H > &A, bool slow=false)
Constructor.
float memory() const
Memory usage in byte.
SparseMatrix(const Space< G > &spcX, const Space< G > &spcY, const BilinearForm< F, G > &bf, const Real eps=0.0, const bool single=false)
Constructor.
const_iterator begin(uint row=0) const
Constant iterator over the elements, standing at position (row,c), where row is the given row number ...
void addIntoT(Matrix< H > &dest, const I fact, const uint rowoffset=0, const uint coloffset=0) const
The transposed of this matrix is added as block to the given matrix.
Values of a matrix in sparse notation.
virtual std::ostream & info(std::ostream &os) const
std::complex< F > type
Definition: typedefs.hh:116
virtual void transpMult(const Vector< r_type > &fncY, Vector< F > &fncX)
Multiplies the transpose of the matrix with fncY and adds the results to fncX.
STL like iterator for hashed sparse matrices.
SparseMatrix(const Space< G > &spc, const BilinearForm< F, G > &bf, const Real eps=0.0)
Constructor.
#define conceptsAssert3(cond, exc, msg)
Assert that a certain condition is fulfilled.
Definition: exceptions.hh:442
const_iterator end() const
Last entrance of the particular order.
void matrixMultiplyRowSorting(const F &factL, const G &factR, Matrix< H > &dest)
Multiplies two matrices, which deliver at least a row sorted iterator, and adds (!...
Definition: matrixMult.hh:28
void write(char *fname) const
Writes the matrix to a file.
virtual F operator()(const uint i, const uint j) const
double Real
Type normally used for a floating point number.
Definition: typedefs.hh:36
virtual F & operator()(const uint i, const uint j)
SparseMatrix(const Space< typename Realtype< F >::type > &spc, const BilinearFormContainer< F > bf, const Real eps=0.0)
Constructor.
Basic namespace for Concepts-2.
Definition: pml_formula.h:16
SparseMatrix(const Space< G > &spc, const F *const v, const int i=-1)
Constructor of partial rank 1 matrix.
SparseMatrix(const Space< G > &spcX, const Space< G > &spcY, const BilinearForm< F, G > &bf, const Sequence< ElementWithCell< G > * > &seq, const Real eps=0.0)
Constructor.
Page URL: http://wiki.math.ethz.ch/bin/view/Concepts/WebHome
21 August 2020
© 2020 Eidgenössische Technische Hochschule Zürich