Hex  1.0
Hydrogen-electron collision solver
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
matrix.h
Go to the documentation of this file.
1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *\
2  * *
3  * / / / / __ \ \ / / *
4  * / /__ / / / _ \ \ \/ / *
5  * / ___ / | |/_/ / /\ \ *
6  * / / / / \_\ / / \ \ *
7  * *
8  * Jakub Benda (c) 2014 *
9  * Charles University in Prague *
10  * *
11 \* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
12 
13 #ifndef HEX_SPMATRIX
14 #define HEX_SPMATRIX
15 
16 #include <cassert>
17 #include <complex>
18 #include <cstring>
19 #include <iostream>
20 #include <vector>
21 
22 #ifndef NO_PNG
23 #include <png++/png.hpp>
24 #endif
25 
26 #include <umfpack.h>
27 
28 #include "arrays.h"
29 #include "complex.h"
30 
31 // forward declaration of classes in order to enable them as return types
32 // of other classes before proper definition
33 template <class T> class RowMatrix;
34 template <class T> class ColMatrix;
35 class CooMatrix;
36 class CscMatrix;
37 class CsrMatrix;
38 class SymDiaMatrix;
39 
50 template <class T> class DenseMatrix
51 {
52  public:
53 
54  // constructors
56  : rows_(0), cols_(0), data_() {}
57  DenseMatrix(int rows, int cols)
58  : rows_(rows), cols_(cols), data_(rows * cols) {}
60  : rows_(rows), cols_(cols), data_(data) { assert(data.size() == (size_t)(rows * cols)); }
61 
62  // explicit conversion to cArrayView
63  ArrayView<T> data () { return data_; }
64  const ArrayView<T> data () const { return data_; }
65 
66  // getters
67  size_t size () const { return rows_ * cols_; }
68  int cols () const { return cols_; }
69  int rows () const { return rows_; }
70  T * begin () { return data_.begin(); }
71  T const * begin () const { return data_.begin(); }
72 
73  private:
74 
76  int rows_;
77 
79  int cols_;
80 
82  NumberArray<T> data_;
83 };
84 
92 template <class Type> class ColMatrix : public DenseMatrix<Type>
93 {
94  public:
95 
96  // constructor
98  : DenseMatrix<Type>() {}
100  : DenseMatrix<Type>(size, size) {}
101  ColMatrix (int rows, int cols)
102  : DenseMatrix<Type>(rows, cols) {}
104  : DenseMatrix<Type>(rows, cols, data) {}
105 
106  explicit ColMatrix (DenseMatrix<Type> const & M)
107  : DenseMatrix<Type>(M) {}
108 
115  template <class Function> void populate (Function f)
116  {
117  // pointer to data
118  Type * ptr = this->begin();
119 
120  // fill the matrix
121  for (int icol = 0; icol < this->cols(); icol++)
122  for (int irow = 0; irow < this->rows(); irow++)
123  *(ptr++) = f(irow,icol);
124  }
125 
133  {
134  return RowMatrix<Type>
135  (
136  this->cols(),
137  this->rows(),
138  this->data()
139  );
140  }
141 
149  {
150  return ArrayView<Type>
151  (
152  this->data(),
153  i * this->rows(),
154  this->rows()
155  );
156  }
157  const ArrayView<Type> col (int i) const
158  {
159  return ArrayView<Type>
160  (
161  this->data(),
162  i * this->rows(),
163  this->rows()
164  );
165  }
167 
169 
170  Type operator() (int i, int j) const { return col(j)[i]; }
171  Type & operator() (int i, int j) { return col(j)[i]; }
173 };
174 
181 template <class Type> class RowMatrix : public DenseMatrix<Type>
182 {
183  public:
184 
185  // constructor
187  : DenseMatrix<Type>() {}
189  : DenseMatrix<Type>(size, size) {}
190  RowMatrix (int rows, int cols)
191  : DenseMatrix<Type>(rows, cols) {}
193  : DenseMatrix<Type>(rows, cols, data) {}
195  : DenseMatrix<Type>(m.rows(), m.cols(), m.data()) { reorder_(); }
196 
197  explicit RowMatrix (DenseMatrix<Type> const & M)
198  : DenseMatrix<Type>(M) {}
199 
209  template <class Function> void populate (Function f)
210  {
211  // pointer to data
212  Type * ptr = this->begin();
213 
214  // fill the matrix
215  for (int irow = 0; irow < this->rows(); irow++)
216  for (int icol = 0; icol < this->cols(); icol++)
217  *(ptr++) = f(irow,icol);
218  }
219 
229  template <class Function> void transform (Function f)
230  {
231  // pointer to data
232  Type * ptr = this->begin();
233 
234  // transform the matrix
235  for (int irow = 0; irow < this->rows(); irow++)
236  for (int icol = 0; icol < this->cols(); icol++)
237  f(irow,icol,*(ptr++));
238  }
239 
247  {
248  return ColMatrix<Type>
249  (
250  this->cols(),
251  this->rows(),
252  this->data()
253  );
254  }
255 
263  {
264  return ArrayView<Type>
265  (
266  this->data(),
267  i * this->cols(),
268  this->cols()
269  );
270  }
271  const ArrayView<Type> row (int i) const
272  {
273  return ArrayView<Type>
274  (
275  this->data(),
276  i * this->cols(),
277  this->cols()
278  );
279  }
281 
283 
284  Type operator() (int i, int j) const { return row(i)[j]; }
285  Type & operator() (int i, int j) { return row(i)[j]; }
287 
296  Type operator() (const ArrayView<Type> u, const ArrayView<Type> v) const
297  {
298  // return value
299  Type sum = 0;
300 
301  // pointer to matrix elements
302  Type const * restrict pA = this->data().data();
303 
304  // pointers to vector elements
305  Type const * const restrict pu = u.data();
306  Type const * const restrict pv = v.data();
307 
308  // for all elements
309  for (int i = 0; i < this->rows(); i++)
310  for (int j = 0; j < this->cols(); j++)
311  sum += pu[i] * (*(pA++)) * pv[j];
312 
313  return sum;
314  }
315 
317  static RowMatrix<Type> Eye (int size)
318  {
319  RowMatrix<Type> M (size, size);
320 
321  for (int i = 0; i < 0; i++)
322  M.data()[i * size + i] = Type(1);
323 
324  return M;
325  }
326 
327  // arithmetic operators
329  {
330  assert (this->rows() == A.rows());
331  assert (this->cols() == A.cols());
332 
333  this->data() += A.data();
334 
335  return *this;
336  }
338  {
339  assert (this->rows() == A.rows());
340  assert (this->cols() == A.cols());
341 
342  this->data() -= A.data();
343 
344  return *this;
345  }
347  {
348  this->data() *= x;
349  return *this;
350  }
351 
371  void write (std::ostream & out, std::string const & pre = "", std::string const & pos = "") const
372  {
373  // data pointer
374  Type const * ptr = this->data().begin();
375 
376  for (int irow = 0; irow < this->rows(); irow++)
377  {
378  out << pre;
379  for (int icol = 0; icol < this->cols(); icol++)
380  {
381  if (*ptr >= 0) out << " "; out << *ptr++ << " ";
382  }
383 
384  out << pos << "\n";
385  }
386  }
387 
388 #ifndef NO_PNG
389 
399  void plot_abs (std::ofstream & out) const
400  {
401  // create empty gray-scale 1-bit image
402  png::image<png::gray_pixel_16> image (this->cols(), this->rows());
403 
404  // skip empty matrix
405  if (this->data().size() > 0)
406  {
407  // find minimal and maximal value
408  double min = std::abs(this->data()[0]), max = std::abs(this->data()[0]);
409  for (Type const & x : this->data())
410  {
411  double absx = std::abs(x);
412  if (absx < min) min = absx;
413  if (absx > max) max = absx;
414  }
415 
416  // for all elements
417  for (int y = 0; y < this->rows(); y++)
418  {
419  for (int x = 0; x < this->cols(); x++)
420  {
421  double fraction = 1. - (std::abs(this->data()[y * this->cols() + x]) - min) / (max - min);
422  image.set_pixel(x, y, std::ceil(65535 * fraction));
423  }
424  }
425  }
426 
427  // save image
428  image.write_stream(out);
429  }
430 #endif
431 
432  private:
433 
435  void reorder_ ()
436  {
437  NumberArray<Type> new_data (this->data().size());
438 
439  for (int irow = 0; irow < this->rows(); irow++)
440  for (int icol = 0; icol < this->cols(); icol++)
441  new_data[irow * this->cols() + icol] = this->data()[icol * this->rows() + irow];
442 
443  this->data() = new_data;
444  }
445 };
446 
447 template <class T> RowMatrix<T> operator + (RowMatrix<T> const & A, RowMatrix<T> const & B) { RowMatrix<T> C(A); C.data() += B.data(); return C; }
448 template <class T> RowMatrix<T> operator - (RowMatrix<T> const & A, RowMatrix<T> const & B) { RowMatrix<T> C(A); C.data() -= B.data(); return C; }
449 template <class T> RowMatrix<T> operator * (T x, RowMatrix<T> const & A) { RowMatrix<T> B(A); B.data() *= x; return B; }
450 template <class T> RowMatrix<T> operator * (RowMatrix<T> const & A, T x) { RowMatrix<T> B(A); B.data() *= x; return B; }
451 template <class T> RowMatrix<T> operator / (RowMatrix<T> const & A, T x) { RowMatrix<T> B(A); B.data() /= x; return B; }
452 
453 // matrix multiplication
454 template <class T> RowMatrix<T> operator * (RowMatrix<T> const & A, ColMatrix<T> const & B);
455 
459 template <class Type> RowMatrix<Type> operator * (RowMatrix<Type> const & A, ColMatrix<Type> const & B)
460 {
461  assert(A.cols() == B.rows());
462 
463  // sizes
464  int rows = A.rows();
465  int cols = B.cols();
466  int comm = A.cols();
467 
468  // output matrix, initialized to zero
469  RowMatrix<Type> C(A.rows(), B.cols());
470 
471  // data pointers
472  Type const * restrict pA = A.begin();
473  Type const * restrict pB = B.begin();
474  Type * restrict pC = C.begin();
475 
476  // for all rows of A
477  for (int irow = 0; irow < rows; irow++)
478  {
479  // for all columns of B
480  for (int icol = 0; icol < cols; icol++)
481  {
482  // compute the scalar product "row * column"
483  for (int k = 0; k < comm; k++)
484  (*pC) += (*(pA + k)) * (*pB++);
485 
486  // move to next element of C
487  pC++;
488  }
489 
490  // move to the next row of A
491  pA += A.cols();
492 
493  // reset B's data pointer
494  pB = B.begin();
495  }
496 
497  // return result
498  return C;
499 }
500 
513 typedef enum {
514  none = 0, // b000,
515  strict_lower = 1, // b001, strict_lower
516  strict_upper = 2, // b010, strict_upper
517  strict_both = 3, // b011, strict_lower | strict_upper
518  diagonal = 4, // b100, diagonal
519  lower = 5, // b101, strict_lower | diagonal
520  upper = 6, // b110, diagonal | strict_upper
521  both = 7 // b111, strict_lower | diagonal | strict_upper
523 
536 {
537 
538 public:
539 
540  // Constructors
541 
543  : m_(0), n_(0) {}
544  CscMatrix(size_t m, size_t n)
545  : m_(m), n_(n) {}
546  CscMatrix(CscMatrix const & A)
547  : m_(A.m_), n_(A.n_), p_(A.p_), i_(A.i_), x_(A.x_) {}
548  CscMatrix(size_t m, size_t n, const lArrayView p, const lArrayView i, const cArrayView x)
549  : m_(m), n_(n), p_(p), i_(i), x_(x) {}
550 
551  // Destructor
552 
554 
558  CooMatrix tocoo() const;
559 
566  cArray dotT(const cArrayView b) const;
567 
568  // Getters
569 
570  size_t size() const { return i_.size(); }
571  size_t rows() const { return m_; }
572  size_t cols() const { return n_; }
573  lArray const & p() const { return p_; }
574  lArray const & i() const { return i_; }
575  cArray const & x() const { return x_; }
576 
580  CscMatrix & operator *= (double r);
581 
587  CscMatrix & operator &= (const CscMatrix& B);
588 
594  CscMatrix& operator ^= (const CscMatrix& B);
595 
600  bool hdfsave(const char* name) const;
601 
606  bool hdfload(const char* name);
607 
608 private:
609 
610  // dimensions
611  long m_;
612  long n_;
613 
614  // representation
615  lArray p_;
616  lArray i_;
617  cArray x_;
618 
619 };
620 
636 class CsrMatrix {
637 
638 public:
639 
640  // Constructors
641 
643  : m_(0), n_(0), name_() {}
644  CsrMatrix(size_t m, size_t n)
645  : m_(m), n_(n), name_() {}
646  CsrMatrix(CsrMatrix const & A)
647  : m_(A.m_), n_(A.n_), p_(A.p_), i_(A.i_), x_(A.x_), name_() {}
648  CsrMatrix(size_t m, size_t n, lArrayView const & p, lArrayView const & i, cArrayView const & x)
649  : m_(m), n_(n), p_(p), i_(i), x_(x), name_() {}
650 
651  // Destructor
652 
654 
658  void drop ()
659  {
660  m_ = n_ = 0;
661  p_.drop();
662  i_.drop();
663  x_.drop();
664  }
665 
670  cArray dot(cArrayView const & b) const;
671 
672  // Getters
673 
674  size_t size() const { return i_.size(); }
675  size_t rows() const { return m_; }
676  size_t cols() const { return n_; }
677  lArray const & p() const { return p_; }
678  lArray const & i() const { return i_; }
679  cArray const & x() const { return x_; }
680 
681  // return absolute value of the (in absolute value) largest element
682  double norm() const;
683 
694  void write(const char* filename) const;
695 
696 #ifndef NO_PNG
697 
703  class PngGenerator : public png::generator<png::gray_pixel_1,PngGenerator>
704  {
705  typedef png::generator<png::gray_pixel_1,PngGenerator> base_t;
706  typedef png::packed_pixel_row<png::gray_pixel_1> row;
707  typedef png::row_traits<row> row_traits;
708 
709  public:
710 
711  PngGenerator (CsrMatrix const * mat, double threshold);
712  ~PngGenerator ();
713 
714  png::byte* get_next_row (size_t pos);
715 
716  private:
717 
718  CsrMatrix const * M_;
719  row buff_;
720  double threshold_;
721  };
722 
728  void plot (const char* filename, double threshold = 0.) const;
729 #endif
730 
731 #ifndef NO_UMFPACK
732 
741  class LUft
742  {
743  public:
744 
746  LUft ()
747  : numeric_(nullptr), matrix_(nullptr), filename_(), info_(UMFPACK_INFO) {}
748 
750  LUft (LUft const & lu)
751  : numeric_(lu.numeric_), matrix_(lu.matrix_), filename_(lu.filename_), info_(lu.info_) {}
752 
754  LUft (LUft && lu)
755  : numeric_(lu.numeric_), matrix_(lu.matrix_), filename_(lu.filename_), info_(lu.info_) { lu.numeric_ = nullptr; }
756 
758  LUft (const CsrMatrix* matrix, void* numeric)
759  : numeric_(numeric), matrix_(matrix), filename_(), info_(UMFPACK_INFO) {}
760 
762  ~LUft ()
763  {
764  drop ();
765  }
766 
773  void transfer (LUft && lu)
774  {
775  // clean this object
776  drop();
777 
778  // transfer data
779  numeric_ = lu.numeric_;
780  matrix_ = lu.matrix_;
781  filename_ = lu.filename_;
782  info_ = lu.info_;
783 
784  // clean source
785  lu.numeric_ = nullptr;
786  }
787 
794  size_t size () const
795  {
796  long lnz, unz, m, n, nz_udiag;
797  long status = umfpack_zl_get_lunz (
798  &lnz, &unz, &m, &n, &nz_udiag, numeric_
799  );
800  return status == 0 ? (lnz + unz) * 16 : 0; // Byte count
801  }
802 
812  cArray solve (const cArrayView b, unsigned eqs = 1) const
813  {
814  // reserve space for the solution
815  cArray x(b.size());
816 
817  // solve
818  solve(b, x, eqs);
819 
820  // return the result
821  return x;
822  }
823  void solve (const cArrayView b, cArrayView x, int eqs = 1) const
824  {
825  // check sizes
826  assert (eqs * matrix_->n_ == (int)x.size());
827  assert (eqs * matrix_->n_ == (int)b.size());
828 
829  // solve for all RHSs
830  for (int eq = 0; eq < eqs; eq++)
831  {
832  // solve for current RHS
833  long status = umfpack_zl_solve (
834  UMFPACK_Aat,
835  matrix_->p_.data(), matrix_->i_.data(),
836  reinterpret_cast<const double*>(matrix_->x_.data()), nullptr,
837  reinterpret_cast<double*>(&x[0] + eq * matrix_->n_), nullptr,
838  reinterpret_cast<const double*>(&b[0] + eq * matrix_->n_), nullptr,
839  numeric_, nullptr, &info_[0]
840  );
841 
842  // check output
843  if (status != UMFPACK_OK)
844  {
845  std::cerr << "\n[CsrMatrix::LUft::solve] Exit status " << status << "\n";
846  umfpack_zl_report_status(0, status);
847  }
848  }
849  }
851 
857  rArray const & info() const { return info_; }
858 
866  void link (std::string name) { filename_ = name; }
867 
874  void save (std::string name) const
875  {
876  long err = umfpack_zl_save_numeric (numeric_, const_cast<char*>(filename_.c_str()));
877 
878  if (err == UMFPACK_ERROR_invalid_Numeric_object)
879  throw exception ("[LUft::save] Invalid numeric object.");
880 
881  if (err == UMFPACK_ERROR_file_IO)
882  throw exception ("[LUft::save] I/O error.");
883  }
884  void save () const
885  {
886  save (filename_);
887  }
888 
895  void load (std::string name)
896  {
897  long err = umfpack_zl_load_numeric (&numeric_, const_cast<char*>(filename_.c_str()));
898 
899  if (err == UMFPACK_ERROR_out_of_memory)
900  throw exception ("[LUft::load] Out of memory.");
901 
902  if (err == UMFPACK_ERROR_file_IO)
903  throw exception ("[LUft::save] I/O error.");
904  }
905  void load ()
906  {
907  load (filename_);
908  }
910 
916  void drop ()
917  {
918  if (numeric_ != nullptr)
919  {
920  umfpack_zl_free_numeric (&numeric_);
921  numeric_ = nullptr;
922  }
923  }
924 
925  private:
926 
928  void* numeric_;
929 
931  const CsrMatrix* matrix_;
932 
934  std::string filename_;
935 
937  mutable rArray info_;
938 
939  // Disable bitwise copy
940  LUft const & operator= (LUft const &);
941  };
942 
956  LUft factorize (double droptol = 0) const;
957 #endif
958 
969  cArray solve (const cArrayView b, size_t eqs = 1) const;
970 
977  void link (std::string name);
978 
983  bool hdfsave () const { return hdfsave(name_); }
984  bool hdfsave (std::string name) const;
986 
991  bool hdfload () { return hdfload(name_); }
992  bool hdfload (std::string name);
994 
999  Complex operator() (unsigned i, unsigned j) const;
1000 
1005 
1011  CsrMatrix & operator &= (CsrMatrix const & B);
1012 
1018  CsrMatrix& operator ^= (CsrMatrix const & B);
1019 
1026  CsrMatrix sparse_like (CsrMatrix const & B) const;
1027 
1031  cArray diag() const;
1032 
1036  CooMatrix tocoo() const;
1037 
1041  RowMatrix<Complex> torow() const;
1042 
1049  cArray upperSolve (cArrayView const & b) const;
1050 
1057  cArray lowerSolve (cArrayView const & b) const;
1058 
1066  template <class Functor> CsrMatrix nzTransform (Functor f) const
1067  {
1068  // create output matrix
1069  CsrMatrix A = *this;
1070 
1071  // loop over rows
1072  for (size_t row = 0; row < rows(); row++)
1073  {
1074  // get relevant columns
1075  size_t idx1 = p_[row];
1076  size_t idx2 = p_[row + 1];
1077 
1078  // loop over columns
1079  for (size_t idx = idx1; idx < idx2; idx++)
1080  {
1081  // which column is this?
1082  size_t col = i_[idx];
1083 
1084  // transform the output matrix
1085  A.x_[idx] = f(row, col, x_[idx]);
1086  }
1087  }
1088 
1089  return A;
1090  }
1091 
1092 private:
1093 
1094  // dimensions
1095  long m_;
1096  long n_;
1097 
1098  // representation
1099  lArray p_;
1100  lArray i_;
1101  cArray x_;
1102 
1103  // linked HDF file
1104  std::string name_;
1105 };
1106 
1107 
1120 class CooMatrix {
1121 
1122 public:
1123 
1124  // Empty constructors
1125 
1127  : m_(0), n_(0), sorted_(true) {}
1128  CooMatrix(size_t m, size_t n)
1129  : m_(m), n_(n), sorted_(true) {}
1131  : m_(A.m_), n_(A.n_), i_(A.i_), j_(A.j_), x_(A.x_), sorted_(false) {}
1132  CooMatrix(size_t m, size_t n, NumberArray<long> const & i, NumberArray<long> const & j, NumberArray<Complex> const & x)
1133  : m_(m), n_(n), i_(i), j_(j), x_(x), sorted_(false) {}
1134 
1142  template <class T> CooMatrix (size_t m, size_t n, T a) : m_(m), n_(n), sorted_(false)
1143  {
1144  // initialize from column-major formatted input
1145  size_t i = 0;
1146  for (size_t col = 0; col < n; col++)
1147  {
1148  for (size_t row = 0; row < m; row++)
1149  {
1150  // get element from array
1151  Complex val = *(a + i);
1152 
1153  // if nonzero, store
1154  if (val != 0.)
1155  {
1156  i_.push_back(row);
1157  j_.push_back(col);
1158  x_.push_back(val);
1159  }
1160 
1161  // move to next element
1162  i++;
1163  }
1164  }
1165  }
1166 
1167  // Destructor
1169 
1171  operator Complex () const
1172  {
1173  if (m_ == 1 and n_ == 1)
1174  {
1175  // get number of nonzero matrix elements
1176  size_t elems = this->shake().x_.size();
1177  if (elems > 1)
1178  throw "[CooMatrix::operator Complex] more elements than nominal volume!\n";
1179  else if (elems == 1)
1180  return this->shake().x_[0];
1181  else
1182  return 0.;
1183  }
1184  else
1185  throw "[CooMatrix::operator Complex] matrix is not 1×1!\n";
1186  }
1187 
1188  // Getters
1189 
1190  size_t rows() const { return m_; }
1191  size_t cols() const { return n_; }
1192  size_t size() const { return i_.size(); }
1193  lArray const & i() const { return i_; }
1194  lArray const & j() const { return j_; }
1195  cArray const & v() const { return x_; }
1196 
1198  Complex operator() (size_t ix, size_t iy) const
1199  {
1200  for (size_t n = 0; n < i_.size(); n++)
1201  if ((size_t)i_[n] == ix and (size_t)j_[n] == iy)
1202  return x_[n];
1203  return 0.;
1204  }
1205 
1218  template <class Functor> CooMatrix& symm_populate_band (size_t d, Functor f)
1219  {
1220  Complex val;
1221 
1222  for (size_t row = 0; row < m_; row++)
1223  {
1224  for (size_t col = row; col < n_ and col - row <= d; col++)
1225  {
1226  val = f(row,col);
1227 
1228  if (val != 0.)
1229  {
1230  i_.push_back(row);
1231  j_.push_back(col);
1232  x_.push_back(val);
1233 
1234  if (row != col)
1235  {
1236  i_.push_back(col);
1237  j_.push_back(row);
1238  x_.push_back(val);
1239  }
1240  }
1241  }
1242  }
1243 
1244  sorted_ = false;
1245 
1246  return *this;
1247  }
1248 
1258  template <class Functor> CooMatrix& populate (Functor f)
1259  {
1260  Complex val;
1261 
1262  for (size_t row = 0; row < m_; row++)
1263  {
1264  for (size_t col = 0; col < n_; col++)
1265  {
1266  val = f(row,col);
1267 
1268  if (val != 0.)
1269  {
1270  i_.push_back(row);
1271  j_.push_back(col);
1272  x_.push_back(val);
1273  }
1274  }
1275  }
1276 
1277  sorted_ = false;
1278 
1279  return *this;
1280  }
1281 
1282  // Assignment
1284  {
1285  m_ = A.m_;
1286  n_ = A.n_;
1287  i_ = A.i_;
1288  j_ = A.j_;
1289  x_ = A.x_;
1290 
1291  sorted_ = A.sorted_;
1292 
1293  return *this;
1294  }
1295 
1302  void add (long i, long j, Complex v)
1303  {
1304  i_.push_back(i);
1305  j_.push_back(j);
1306  x_.push_back(v);
1307 
1308  sorted_ = false;
1309  }
1310 
1313  {
1314  CooMatrix tr;
1315 
1316  tr.m_ = n_;
1317  tr.n_ = m_;
1318  tr.i_ = j_;
1319  tr.j_ = i_;
1320  tr.x_ = x_;
1321 
1322  tr.sorted_ = false;
1323 
1324  return tr;
1325  }
1326 
1329  {
1330  assert(m_ == A.m_);
1331  assert(n_ == A.n_);
1332 
1333  i_.append(A.i_.begin(), A.i_.end());
1334  j_.append(A.j_.begin(), A.j_.end());
1335  x_.append(A.x_.begin(), A.x_.end());
1336 
1337  sorted_ = false;
1338 
1339  return *this;
1340  }
1341 
1344  {
1345  assert(m_ == A.m_);
1346  assert(n_ == A.n_);
1347 
1348  size_t prev_size = x_.size();
1349 
1350  i_.append(A.i_.begin(), A.i_.end());
1351  j_.append(A.j_.begin(), A.j_.end());
1352  x_.append(A.x_.begin(), A.x_.end());
1353 
1354  // negate the newly added elements
1355  for (size_t i = prev_size; i < x_.size(); i++)
1356  x_[i] = -x_[i];
1357 
1358  sorted_ = false;
1359 
1360  return *this;
1361  }
1362 
1365  {
1366  long nz = i_.size();
1367  for (long i = 0; i < nz; i++)
1368  x_[i] *= c;
1369 
1370  return *this;
1371  }
1372 
1374  CooMatrix dot (const cArrayView B) const;
1375 
1383  Complex ddot (CooMatrix const & B) const;
1384 
1393  CooMatrix& operator *= (const cArrayView B);
1394 
1397  {
1398  return *this *= 1./c;
1399  }
1400 
1405  void resize (size_t m, size_t n)
1406  {
1407  m_ = m;
1408  n_ = n;
1409  }
1410 
1421  CooMatrix reshape (size_t m, size_t n) const;
1422 
1424  cArray todense () const;
1425 
1427  void sort();
1428  bool sorted() const { return sorted_; }
1429 
1431  CscMatrix tocsc() const;
1432 
1434  CsrMatrix tocsr() const;
1435 
1437  template <typename DenseMatrixType> DenseMatrixType todense() const
1438  {
1439  DenseMatrixType M (rows(), cols());
1440  for (unsigned idx = 0; idx < x_.size(); idx++)
1441  M (i_[idx], j_[idx]) = x_[idx];
1442  return M;
1443  }
1444 
1447  {
1448  return todense<RowMatrix<Complex>>();
1449  }
1450 
1453  {
1454  return todense<ColMatrix<Complex>>();
1455  }
1456 
1464  SymDiaMatrix todia(MatrixTriangle triangle = lower) const;
1465 
1475  cArray solve (const cArrayView b, size_t eqs = 1) const
1476  {
1477  // COO format is not optimal for solving -> covert to CSC
1478  return tocsr().solve(b, eqs);
1479  }
1480 
1491  void write(const char* filename) const;
1492 
1494  CooMatrix shake() const;
1495 
1500  bool hdfsave(const char* name) const;
1501 
1506  bool hdfload(const char* name);
1507 
1508 private:
1509 
1510  // dimensions
1511  size_t m_, n_;
1512 
1513  // ijv-representation
1514  lArray i_, j_;
1515  cArray x_;
1516 
1517  bool sorted_;
1518 };
1519 
1548 
1549 public:
1550  //
1551  // Constructors
1552  //
1553 
1555  SymDiaMatrix();
1556 
1558  SymDiaMatrix(int n);
1559 
1567  SymDiaMatrix(int n, const iArrayView id, const cArrayView v);
1568 
1570  SymDiaMatrix(SymDiaMatrix const & A);
1571 
1573  SymDiaMatrix(SymDiaMatrix && A);
1574 
1589  template <class Functor> SymDiaMatrix & populate(unsigned d, Functor f)
1590  {
1591  // throw away old data
1592  elems_.resize(0);
1593 
1594  // for all diagonals
1595  for (size_t id = 0; id <= d; id++)
1596  {
1597  // add this diagonal
1598  idiag_.push_back(id);
1599 
1600  // for all elements of the diagonal
1601  for (int icol = id; icol < n_; icol++)
1602  {
1603  // get the row index, too
1604  int irow = icol - id;
1605 
1606  // evaluate the element
1607  elems_.push_back(f(irow,icol));
1608  }
1609  }
1610 
1611  return *this;
1612  }
1613 
1614  //
1615  // Destructor
1616  //
1617 
1619 
1621  void drop ()
1622  {
1623  n_ = 0;
1624  elems_.drop();
1625  idiag_.drop();
1626  dptrs_.clear();
1627  }
1628 
1629  //
1630  // Getters
1631  //
1632 
1642  iArray const & diag() const { return idiag_; }
1643  iArray & diag() { return idiag_; }
1645 
1653  int diag(int i) const { return idiag_[i]; }
1654 
1661  cArrayView main_diagonal() const { return cArrayView(elems_, 0, n_); }
1662  cArrayView main_diagonal() { return cArrayView(elems_, 0, n_); }
1664 
1671  cArray const & data() const { return elems_; }
1672  cArray & data() { return elems_; }
1674 
1683  Complex const * dptr(int i) const { return dptrs_[i]; }
1684  Complex * dptr(int i) { return dptrs_[i]; }
1686 
1693  size_t size() const { return n_; }
1694 
1701  int bandwidth() const { return 1 + 2 * idiag_.back(); }
1702 
1710  bool is_compatible (SymDiaMatrix const & B) const;
1711 
1712  //
1713  // Arithmetic and other operators
1714  //
1715 
1716  SymDiaMatrix const & operator = (SymDiaMatrix && A);
1717  SymDiaMatrix const & operator = (SymDiaMatrix const & A);
1718 
1719  SymDiaMatrix const & operator += (SymDiaMatrix const & B);
1720  SymDiaMatrix const & operator -= (SymDiaMatrix const & B);
1721 
1737  cArray dot (const cArrayView B, MatrixTriangle triangle = both, bool parallelize = false) const;
1738 
1747  cArray lowerSolve (const cArrayView b) const;
1748 
1757  cArray upperSolve (const cArrayView b) const;
1758 
1764  SymDiaMatrix kron (SymDiaMatrix const & B) const;
1765 
1766  //
1767  // HDF interface
1768  //
1769 
1771  void link (std::string name);
1772 
1774  std::string linkedto () const { return name_; }
1775 
1784  bool hdfload () { return hdfload (name_); }
1785  bool hdfload (std::string name);
1787 
1800  bool hdfsave () const { return hdfsave (name_); }
1801  bool hdfsave (std::string name, bool docompress = false, int consec = 10) const;
1803 
1804  //
1805  // Conversions to other formats
1806  //
1807 
1867  cArray toPaddedRows () const;
1868 
1875  cArray toPaddedCols () const;
1876 
1878  CooMatrix tocoo (MatrixTriangle triangle = both) const;
1879 
1881  RowMatrix<Complex> torow (MatrixTriangle triangle = both) const;
1882 
1884  friend std::ostream & operator << (std::ostream & out, SymDiaMatrix const & A);
1885 
1886 private:
1887 
1888  // dimension (only square matrices allowed)
1889  int n_;
1890 
1891  // diagonals: concatenated diagonals starting from the longest
1892  // to the shortest (i.e. with rising right index)
1893  cArray elems_;
1894 
1895  // upper diagonal indices starting from zero
1896  iArray idiag_;
1897 
1898  // name of linked HDF file
1899  std::string name_;
1900 
1901  // diagonal data pointers
1902  std::vector<Complex*> dptrs_;
1903 
1917  void setup_dptrs_();
1918 };
1919 
1920 // --------------------------------------------------------------------------//
1921 // ---- Binary arithmetic operators ---------------------------------------- //
1922 // --------------------------------------------------------------------------//
1923 
1927 inline CsrMatrix operator & (const CsrMatrix& A, const CsrMatrix& B)
1928 {
1929  CsrMatrix C = A;
1930  return C &= B;
1931 }
1932 
1936 inline CsrMatrix operator ^ (const CsrMatrix& A, const CsrMatrix& B)
1937 {
1938  CsrMatrix C = A;
1939  return C ^= B;
1940 }
1941 
1946 {
1947  CsrMatrix C = B;
1948  return C *= z;
1949 }
1950 
1954 inline CsrMatrix operator * (const CsrMatrix& A, double r)
1955 {
1956  CsrMatrix C = A;
1957  return C *= r;
1958 }
1959 
1963 inline CscMatrix operator & (const CscMatrix& A, const CscMatrix& B)
1964 {
1965  CscMatrix C = A;
1966  return C &= B;
1967 }
1968 
1972 inline CscMatrix operator ^ (const CscMatrix& A, const CscMatrix& B)
1973 {
1974  CscMatrix C = A;
1975  return C ^= B;
1976 }
1977 
1981 inline CscMatrix operator * (double r, const CscMatrix& B)
1982 {
1983  CscMatrix C = B;
1984  return C *= r;
1985 }
1986 
1990 inline CscMatrix operator * (const CscMatrix& A, double r)
1991 {
1992  CscMatrix C = A;
1993  return C *= r;
1994 }
1995 
2001 inline CooMatrix operator + (const CooMatrix& A, const CooMatrix& B)
2002 {
2003  CooMatrix C = A;
2004  return C += B;
2005 }
2006 
2012 inline CooMatrix operator - (const CooMatrix& A, const CooMatrix& B)
2013 {
2014  CooMatrix C = A;
2015  return C -= B;
2016 }
2017 
2023 inline CooMatrix operator * (const Complex& z, const CooMatrix& B)
2024 {
2025  CooMatrix C = B;
2026  return C *= z;
2027 }
2028 
2034 inline CooMatrix operator * (CooMatrix const & A, const cArrayView B)
2035 {
2036  CooMatrix C = A;
2037  return C *= B;
2038 }
2039 
2045 inline CooMatrix operator * (const CooMatrix& A, const Complex& z)
2046 {
2047  CooMatrix C = A;
2048  return C *= z;
2049 }
2050 
2051 SymDiaMatrix operator + (SymDiaMatrix const & A, SymDiaMatrix const & B);
2052 SymDiaMatrix operator - (SymDiaMatrix const & A, SymDiaMatrix const & B);
2053 SymDiaMatrix operator * (SymDiaMatrix const & A, SymDiaMatrix const & B);
2054 SymDiaMatrix operator * (double z, SymDiaMatrix const & A);
2056 
2074 CooMatrix kron(const CooMatrix& A, const CooMatrix& B);
2075 
2080 CooMatrix eye(size_t N);
2081 
2086 CooMatrix stairs(size_t N);
2087 
2088 #endif
Definition: matrix.h:516
virtual T * data()
Data pointer.
Definition: arrays.h:757
static RowMatrix< Type > Eye(int size)
Identity matrix.
Definition: matrix.h:317
CooMatrix(size_t m, size_t n)
Definition: matrix.h:1128
Array view.
Definition: arrays.h:186
DenseMatrix()
Definition: matrix.h:55
ArrayView< Type > row(int i)
Matrix row.
Definition: matrix.h:262
CooMatrix tocoo() const
Definition: matrix.cpp:756
CooMatrix & symm_populate_band(size_t d, Functor f)
Symmetrical band populator.
Definition: matrix.h:1218
RowMatrix< Type > const & operator-=(DenseMatrix< Type > const &A)
Definition: matrix.h:337
Dense (column-oriented) matrix.
Definition: matrix.h:34
Definition: matrix.h:514
SymDiaMatrix operator+(SymDiaMatrix const &A, SymDiaMatrix const &B)
Definition: matrix.cpp:1478
lArray const & i() const
Definition: matrix.h:574
const ArrayView< Type > row(int i) const
Definition: matrix.h:271
void drop()
Free all fields, set dimensions to zero.
Definition: matrix.h:1621
RowMatrix(ColMatrix< Type > const &m)
Definition: matrix.h:194
bool sorted() const
Definition: matrix.h:1428
const ArrayView< T > data() const
Definition: matrix.h:64
DenseMatrix(int rows, int cols, const ArrayView< T > data)
Definition: matrix.h:59
RowMatrix(int rows, int cols)
Definition: matrix.h:190
size_t rows() const
Definition: matrix.h:1190
CscMatrix(CscMatrix const &A)
Definition: matrix.h:546
RowMatrix< Complex > torow() const
Definition: matrix.cpp:801
rArray threshold(const rArrayView a, double eps)
Drop small elements of array (replace by zero).
Definition: arrays.cpp:136
void transfer(LUft &&lu)
Transfer data from another factorization object.
Definition: matrix.h:773
CooMatrix kron(const CooMatrix &A, const CooMatrix &B)
Kronecker product.
Definition: matrix.cpp:83
cArrayView main_diagonal()
Definition: matrix.h:1662
Dense (row-oriented) matrix.
Definition: matrix.h:33
RowMatrix(int rows, int cols, const ArrayView< Type > data)
Definition: matrix.h:192
Definition: matrix.h:519
CooMatrix reshape(size_t m, size_t n) const
Change matrix shape.
Definition: matrix.cpp:1111
T const * begin() const
Definition: matrix.h:71
int bandwidth() const
Bandwidth.
Definition: matrix.h:1701
virtual size_t resize(size_t n)
Resize array.
Definition: arrays.h:685
cArray dotT(const cArrayView b) const
Definition: matrix.cpp:207
RowMatrix< Complex > torow() const
Convert to dense matrix (row-ordered).
Definition: matrix.h:1446
T min(const ArrayView< T > a)
Minimal element of array.
Definition: arrays.h:1663
CooMatrix transpose() const
Transposition, implemented as an interchange of &quot;i&quot; and &quot;j&quot; data.
Definition: matrix.h:1312
RowMatrix< Type > const & operator*=(Type x)
Definition: matrix.h:346
DenseMatrixType todense() const
Convert to dense matrix of a given underlying type.
Definition: matrix.h:1437
~PngGenerator()
Definition: matrix.cpp:431
CooMatrix tocoo(MatrixTriangle triangle=both) const
Convert matrix part to CooMatrix.
Definition: matrix.cpp:1670
cArray const & x() const
Definition: matrix.h:575
#define restrict
Definition: misc.h:88
LUft(LUft &&lu)
Move constructor.
Definition: matrix.h:754
CsrMatrix operator^(const CsrMatrix &A, const CsrMatrix &B)
Definition: matrix.h:1936
bool hdfload()
Load from file.
Definition: matrix.h:1784
Complex CSC matrix.
Definition: matrix.h:535
Complex CSR matrix.
Definition: matrix.h:636
void resize(size_t m, size_t n)
Change dimension of the matrix.
Definition: matrix.h:1405
CooMatrix dot(const cArrayView B) const
SpMV multiplication.
Definition: matrix.cpp:1162
friend std::ostream & operator<<(std::ostream &out, SymDiaMatrix const &A)
Output to a text stream.
Definition: matrix.cpp:1883
LUft()
Default constructor.
Definition: matrix.h:746
Complex ddot(CooMatrix const &B) const
Double inner matrix-matrix product.
Definition: matrix.cpp:1200
CooMatrix stairs(size_t N)
Definition: matrix.cpp:143
RowMatrix< Complex > torow(MatrixTriangle triangle=both) const
Convert matrix part to RowMatrix.
Definition: matrix.cpp:1899
void write(std::ostream &out, std::string const &pre="", std::string const &pos="") const
Output to file.
Definition: matrix.h:371
lArray const & i() const
Definition: matrix.h:678
CooMatrix eye(size_t N)
Definition: matrix.cpp:133
cArray toPaddedRows() const
Zero-pad rows.
Definition: matrix.cpp:1931
Type operator()(int i, int j) const
Element access.
Definition: matrix.h:170
RowMatrix(int size)
Definition: matrix.h:188
cArray & data()
Definition: matrix.h:1672
void link(std::string name)
Link matrix to a disk file.
Definition: matrix.cpp:2009
SymDiaMatrix const & operator=(SymDiaMatrix &&A)
Definition: matrix.cpp:1460
CooMatrix shake() const
Shake the content, i.e. sum same element entries.
Definition: matrix.cpp:1241
void plot(const char *filename, double threshold=0.) const
Definition: matrix.cpp:455
lArray const & j() const
Definition: matrix.h:1194
void link(std::string name)
Link to a disk file.
Definition: matrix.h:866
int cols() const
Definition: matrix.h:68
T & back(int i=0)
Definition: arrays.h:779
bool hdfload(const char *name)
Definition: matrix.cpp:317
CsrMatrix(size_t m, size_t n, lArrayView const &p, lArrayView const &i, cArrayView const &x)
Definition: matrix.h:648
iterator begin()
Definition: arrays.h:288
rArray const & info() const
Get info array.
Definition: matrix.h:857
ArrayView< T > data()
Definition: matrix.h:63
iArray const & diag() const
Diagonal indices.
Definition: matrix.h:1642
SymDiaMatrix operator-(SymDiaMatrix const &A, SymDiaMatrix const &B)
Definition: matrix.cpp:1484
RowMatrix()
Definition: matrix.h:186
CscMatrix & operator*=(double r)
Definition: matrix.cpp:161
iterator end()
Definition: arrays.h:774
CooMatrix(size_t m, size_t n, T a)
Definition: matrix.h:1142
size_t size() const
Size of the numerical data.
Definition: matrix.h:794
cArray const & x() const
Definition: matrix.h:679
cArray lowerSolve(const cArrayView b) const
Back-substitution (lower).
Definition: matrix.cpp:1801
const ArrayView< Type > col(int i) const
Definition: matrix.h:157
bool hdfsave() const
Definition: matrix.h:983
size_t rows() const
Definition: matrix.h:571
size_t size() const
Definition: matrix.h:1192
void add(long i, long j, Complex v)
Addition of an element to matrix.
Definition: matrix.h:1302
T * begin()
Definition: matrix.h:70
void save() const
Definition: matrix.h:884
~SymDiaMatrix()
Definition: matrix.h:1618
cArray lowerSolve(cArrayView const &b) const
Definition: matrix.cpp:696
CsrMatrix tocsr() const
Convert to CSR matrix.
Definition: matrix.cpp:992
LU factorization object.
Definition: matrix.h:741
Type operator()(int i, int j) const
Element access.
Definition: matrix.h:284
CsrMatrix sparse_like(CsrMatrix const &B) const
Definition: matrix.cpp:821
size_t cols() const
Definition: matrix.h:572
DenseMatrix.
Definition: matrix.h:50
~CsrMatrix()
Definition: matrix.h:653
LUft(LUft const &lu)
Copy constructor.
Definition: matrix.h:750
ColMatrix< Complex > tocol() const
Convert to dense matrix (column-ordered).
Definition: matrix.h:1452
lArray const & i() const
Definition: matrix.h:1193
CooMatrix tocoo() const
Definition: matrix.cpp:233
Complex * dptr(int i)
Definition: matrix.h:1684
cArray solve(const cArrayView b, size_t eqs=1) const
Solve the Ax = b problem, where &quot;b&quot; can be a matrix.
Definition: matrix.cpp:517
png::byte * get_next_row(size_t pos)
Definition: matrix.cpp:435
CooMatrix & operator+=(CooMatrix const &A)
Addition.
Definition: matrix.h:1328
ColMatrix< Type > T() const
Matrix transpose.
Definition: matrix.h:246
Definition: matrix.h:518
SymDiaMatrix kron(SymDiaMatrix const &B) const
Kronecker product.
Definition: matrix.cpp:1795
Definition: matrix.h:517
void transform(Function f)
Transformator.
Definition: matrix.h:229
bool hdfsave(const char *name) const
Save matrix to HDF file.
Definition: matrix.cpp:1248
LUft(const CsrMatrix *matrix, void *numeric)
Initialize the structure using the matrix and its numeric decomposition.
Definition: matrix.h:758
size_t size() const
Length of the array (number of elements).
Definition: arrays.h:276
~LUft()
Destructor.
Definition: matrix.h:762
CscMatrix(size_t m, size_t n)
Definition: matrix.h:544
CscMatrix & operator^=(const CscMatrix &B)
Definition: matrix.cpp:189
CooMatrix & operator=(CooMatrix const &A)
Definition: matrix.h:1283
void drop()
Free memory.
Definition: matrix.h:916
SymDiaMatrix()
Empty constructor.
Definition: matrix.cpp:1326
CooMatrix & operator-=(CooMatrix const &A)
Subtraction.
Definition: matrix.h:1343
MatrixTriangle
Matrix parts.
Definition: matrix.h:513
T sum(const ArrayView< T > v)
Sum elements in array.
Definition: arrays.h:1770
CooMatrix & operator/=(Complex c)
Element-wise divide by a complex number.
Definition: matrix.h:1396
bool hdfsave(const char *name) const
Definition: matrix.cpp:279
CsrMatrix()
Definition: matrix.h:642
iterator begin()
Definition: arrays.h:771
cArray todense() const
Convert matrix to dense column-major ordered 1D-array.
Definition: matrix.cpp:1142
cArray dot(cArrayView const &b) const
Definition: matrix.cpp:401
void plot_abs(std::ofstream &out) const
Plot to PNG file.
Definition: matrix.h:399
ArrayView< Complex > cArrayView
Definition: arrays.h:1615
Definition: matrix.h:703
void solve(const cArrayView b, cArrayView x, int eqs=1) const
Definition: matrix.h:823
cArray upperSolve(const cArrayView b) const
Back-substitution (upper).
Definition: matrix.cpp:1842
auto operator/(std::complex< C1 > a, std::complex< C2 > b) -> std::complex< decltype(C1(0.)/C2(0.))>
Definition: complex.h:29
PngGenerator(CsrMatrix const *mat, double threshold)
Definition: matrix.cpp:426
~CooMatrix()
Definition: matrix.h:1168
ColMatrix(int size)
Definition: matrix.h:99
cArray diag() const
Definition: matrix.cpp:743
cArray solve(const cArrayView b, size_t eqs=1) const
Solve matrix equation.
Definition: matrix.h:1475
lArray const & p() const
Definition: matrix.h:677
cArray solve(const cArrayView b, unsigned eqs=1) const
Solve equations.
Definition: matrix.h:812
size_t size() const
Definition: matrix.h:570
SymDiaMatrix const & operator+=(SymDiaMatrix const &B)
Definition: matrix.cpp:1352
void link(std::string name)
Link to a disk file.
Definition: matrix.cpp:548
size_t cols() const
Definition: matrix.h:1191
LUft factorize(double droptol=0) const
Compute (incomplete) LU factorization.
Definition: matrix.cpp:469
CscMatrix()
Definition: matrix.h:542
RowMatrix< Type > const & operator+=(DenseMatrix< Type > const &A)
Definition: matrix.h:328
void load()
Definition: matrix.h:905
cArray const & v() const
Definition: matrix.h:1195
int rows() const
Definition: matrix.h:69
ColMatrix(int rows, int cols)
Definition: matrix.h:101
CsrMatrix(CsrMatrix const &A)
Definition: matrix.h:646
void populate(Function f)
Populator.
Definition: matrix.h:209
RowMatrix(DenseMatrix< Type > const &M)
Definition: matrix.h:197
CsrMatrix & operator*=(Complex r)
Definition: matrix.cpp:362
void write(const char *filename) const
Write the matrix data to a file.
Definition: matrix.cpp:1100
Complex operator()(unsigned i, unsigned j) const
Definition: matrix.cpp:848
size_t size() const
Matrix dimension.
Definition: matrix.h:1693
CooMatrix()
Definition: matrix.h:1126
void sort()
Sort indices (by i_, then by j_)
void append(InputIterator first, InputIterator last)
Append a range of values at end.
Definition: arrays.h:877
cArray upperSolve(cArrayView const &b) const
Definition: matrix.cpp:648
auto operator*(std::complex< C1 > a, std::complex< C2 > b) -> std::complex< decltype(C1(0.)*C2(0.))>
Definition: complex.h:24
CLArrayView exception
iArray & diag()
Definition: matrix.h:1643
SymDiaMatrix & populate(unsigned d, Functor f)
Plain symmetrical populator.
Definition: matrix.h:1589
size_t rows() const
Definition: matrix.h:675
void save(std::string name) const
Save Numeric object to a disk file.
Definition: matrix.h:874
size_t size() const
Definition: matrix.h:67
CsrMatrix & operator^=(CsrMatrix const &B)
Definition: matrix.cpp:387
virtual T * data()
Pointer to the data.
Definition: arrays.h:280
ArrayView< Type > col(int i)
Matrix column.
Definition: matrix.h:148
bool hdfsave() const
Save data to file.
Definition: matrix.h:1800
double norm() const
Definition: matrix.cpp:629
std::string linkedto() const
Return the name of the linked disk file.
Definition: matrix.h:1774
cArray dot(const cArrayView B, MatrixTriangle triangle=both, bool parallelize=false) const
Dot product.
Definition: matrix.cpp:1722
Definition: matrix.h:520
Complex const * dptr(int i) const
Pointer to diagonal data.
Definition: matrix.h:1683
lArray const & p() const
Definition: matrix.h:573
ColMatrix(int rows, int cols, const ArrayView< Type > data)
Definition: matrix.h:103
CooMatrix(size_t m, size_t n, NumberArray< long > const &i, NumberArray< long > const &j, NumberArray< Complex > const &x)
Definition: matrix.h:1132
CsrMatrix operator&(const CsrMatrix &A, const CsrMatrix &B)
Definition: matrix.h:1927
CsrMatrix(size_t m, size_t n)
Definition: matrix.h:644
cArray const & data() const
Data pointer.
Definition: matrix.h:1671
~CscMatrix()
Definition: matrix.h:553
SymDiaMatrix const & operator-=(SymDiaMatrix const &B)
Definition: matrix.cpp:1405
Complex COO matrix.
Definition: matrix.h:1120
CscMatrix(size_t m, size_t n, const lArrayView p, const lArrayView i, const cArrayView x)
Definition: matrix.h:548
CooMatrix & populate(Functor f)
Full populator.
Definition: matrix.h:1258
Definition: matrix.h:521
void load(std::string name)
Load Numeric object from a disk file.
Definition: matrix.h:895
Definition: matrix.h:515
cArrayView main_diagonal() const
Main diagonal.
Definition: matrix.h:1661
void drop()
Clear data.
Definition: matrix.h:658
virtual void push_back(T const &a)
Add element to end.
Definition: arrays.h:833
rArray abs(const cArrayView u)
Definition: arrays.cpp:19
bool is_compatible(SymDiaMatrix const &B) const
Check compatibility of matrices.
Definition: matrix.cpp:1566
std::complex< double > Complex
Definition: complex.h:20
ColMatrix(DenseMatrix< Type > const &M)
Definition: matrix.h:106
bool hdfload()
Definition: matrix.h:991
Symmetric diagonal matrix.
Definition: matrix.h:1547
void write(const char *filename) const
Definition: matrix.cpp:533
ColMatrix()
Definition: matrix.h:97
T max(const ArrayView< T > a)
Maximal element of array.
Definition: arrays.h:1673
Complex operator()(size_t ix, size_t iy) const
Index operator. Returns the existing value or zero.
Definition: matrix.h:1198
cArray toPaddedCols() const
Zero-pad columns.
Definition: matrix.cpp:1970
CscMatrix tocsc() const
Convert to CSC matrix.
Definition: matrix.cpp:949
CooMatrix & operator*=(Complex c)
Element-wise multiplication by complex number.
Definition: matrix.h:1364
size_t size() const
Definition: matrix.h:674
DenseMatrix(int rows, int cols)
Definition: matrix.h:57
void drop()
Reset array: deallocate everything, resize to zero.
Definition: arrays.h:918
bool hdfload(const char *name)
Load matrix from HDF file.
Definition: matrix.cpp:1286
SymDiaMatrix todia(MatrixTriangle triangle=lower) const
Convert to symmetric DIA format.
Definition: matrix.cpp:1044
size_t cols() const
Definition: matrix.h:676
CooMatrix(CooMatrix const &A)
Definition: matrix.h:1130
int diag(int i) const
Diagonal index.
Definition: matrix.h:1653
size_t size() const
Item count.
Definition: arrays.h:673
CsrMatrix & operator&=(CsrMatrix const &B)
Definition: matrix.cpp:372
void populate(Function f)
Populator.
Definition: matrix.h:115
CscMatrix & operator&=(const CscMatrix &B)
Definition: matrix.cpp:171
CsrMatrix nzTransform(Functor f) const
Definition: matrix.h:1066
RowMatrix< Type > T() const
Matrix transpose.
Definition: matrix.h:132