ergo
MatrixSymmetric.h
Go to the documentation of this file.
1 /* Ergo, version 3.8, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2019 Elias Rudberg, Emanuel H. Rubensson, Pawel Salek,
4  * and Anastasia Kruchinina.
5  *
6  * This program is free software: you can redistribute it and/or modify
7  * it under the terms of the GNU General Public License as published by
8  * the Free Software Foundation, either version 3 of the License, or
9  * (at your option) any later version.
10  *
11  * This program is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  * GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License
17  * along with this program. If not, see <http://www.gnu.org/licenses/>.
18  *
19  * Primary academic reference:
20  * Ergo: An open-source program for linear-scaling electronic structure
21  * calculations,
22  * Elias Rudberg, Emanuel H. Rubensson, Pawel Salek, and Anastasia
23  * Kruchinina,
24  * SoftwareX 7, 107 (2018),
25  * <http://dx.doi.org/10.1016/j.softx.2018.03.005>
26  *
27  * For further information about Ergo, see <http://www.ergoscf.org>.
28  */
29 
38 #ifndef MAT_MatrixSymmetric
39 #define MAT_MatrixSymmetric
40 
41 #include <algorithm>
42 
43 #include "MatrixBase.h"
44 #include "Interval.h"
46 #include "mat_utils.h"
47 #include "truncation.h"
48 
49 namespace mat {
67  template<typename Treal, typename Tmatrix>
68  class MatrixSymmetric : public MatrixBase<Treal, Tmatrix> {
69  public:
71  typedef Treal real;
72 
74  :MatrixBase<Treal, Tmatrix>() {}
75  /* In the code we are using std::vector<MatrixSymmetric> which in the c++ standard before c++11 requires move operation like T x_new = x which calls implicitly the copy constructor. To make it work with g++ versions without c++11 support we remove the keyword explicit. */
76 #if __cplusplus >= 201103L
78  :MatrixBase<Treal, Tmatrix>(symm) {}
79 #else
81  :MatrixBase<Treal, Tmatrix>(symm) {}
82 #endif
83  explicit MatrixSymmetric(const XY<Treal, MatrixSymmetric<Treal, Tmatrix> >& sm)
84  :MatrixBase<Treal, Tmatrix>() { *this = sm.A * sm.B; }
86  :MatrixBase<Treal, Tmatrix>(matr) {
87  this->matrixPtr->nosymToSym();
88  }
91 #if 0
92  template<typename Tfull>
93  inline void assign_from_full
94  (Tfull const* const fullmatrix,
95  int const totnrows, int const totncols) {
96  assert(totnrows == totncols);
97  this->matrixPtr->assign_from_full(fullmatrix, totnrows, totncols);
98  this->matrixPtr->nosym_to_sym();
99  }
100  inline void assign_from_full
101  (Treal const* const fullmatrix,
102  int const totnrows, int const totncols) {
103  assert(totnrows == totncols);
104  this->matrixPtr->assign_from_full(fullmatrix, totnrows, totncols);
105  this->matrixPtr->nosym_to_sym();
106  }
107 #endif
108 
109  inline void assignFromFull
110  (std::vector<Treal> const & fullMat) {
111  assert((int)fullMat.size() == this->get_nrows() * this->get_ncols());
112  this->matrixPtr->assignFromFull(fullMat);
113  this->matrixPtr->nosymToSym();
114  }
115 
116  inline void assignFromFull
117  (std::vector<Treal> const & fullMat,
118  std::vector<int> const & rowPermutation,
119  std::vector<int> const & colPermutation) {
120  assert((int)fullMat.size() == this->get_nrows() * this->get_ncols());
121  std::vector<int> rowind(fullMat.size());
122  std::vector<int> colind(fullMat.size());
123  int ind = 0;
124  for (int col = 0; col < this->get_ncols(); ++col)
125  for (int row = 0; row < this->get_nrows(); ++row) {
126  rowind[ind] = row;
127  colind[ind] = col;
128  ++ind;
129  }
130  this->assign_from_sparse(rowind,
131  colind,
132  fullMat,
133  rowPermutation,
134  colPermutation);
135  this->matrixPtr->nosymToSym();
136  }
137 
138  inline void fullMatrix(std::vector<Treal> & fullMat) const {
139  this->matrixPtr->syFullMatrix(fullMat);
140  }
147  inline void fullMatrix
148  (std::vector<Treal> & fullMat,
149  std::vector<int> const & rowInversePermutation,
150  std::vector<int> const & colInversePermutation) const {
151  std::vector<int> rowind;
152  std::vector<int> colind;
153  std::vector<Treal> values;
154  get_all_values(rowind, colind, values,
155  rowInversePermutation,
156  colInversePermutation);
157  fullMat.assign(this->get_nrows()*this->get_ncols(),0);
158  assert(rowind.size() == values.size());
159  assert(rowind.size() == colind.size());
160  for (unsigned int ind = 0; ind < values.size(); ++ind) {
161  assert(rowind[ind] + this->get_nrows() * colind[ind] <
162  this->get_nrows()*this->get_ncols());
163  fullMat[rowind[ind] + this->get_nrows() * colind[ind]] =
164  values[ind];
165  fullMat[colind[ind] + this->get_nrows() * rowind[ind]] =
166  values[ind];
167  }
168  }
175  inline void assign_from_sparse
176  (std::vector<int> const & rowind,
177  std::vector<int> const & colind,
178  std::vector<Treal> const & values) {
179  this->matrixPtr->syAssignFromSparse(rowind, colind, values);
180  }
192  inline void assign_from_sparse
193  (std::vector<int> const & rowind,
194  std::vector<int> const & colind,
195  std::vector<Treal> const & values,
196  SizesAndBlocks const & newRows,
197  SizesAndBlocks const & newCols) {
198  this->resetSizesAndBlocks(newRows, newCols);
199  this->matrixPtr->syAssignFromSparse(rowind, colind, values);
200  }
210  inline void assign_from_sparse
211  (std::vector<int> const & rowind,
212  std::vector<int> const & colind,
213  std::vector<Treal> const & values,
214  std::vector<int> const & rowPermutation,
215  std::vector<int> const & colPermutation) {
216  std::vector<int> newRowind;
217  std::vector<int> newColind;
218  this->getPermutedAndSymmetrized(rowind, rowPermutation, newRowind,
219  colind, colPermutation, newColind);
220 
221  this->matrixPtr->syAssignFromSparse(newRowind, newColind, values);
222  }
228  inline void assign_from_sparse
229  (std::vector<int> const & rowind,
230  std::vector<int> const & colind,
231  std::vector<Treal> const & values,
232  SizesAndBlocks const & newRows,
233  SizesAndBlocks const & newCols,
234  std::vector<int> const & rowPermutation,
235  std::vector<int> const & colPermutation) {
236  this->resetSizesAndBlocks(newRows, newCols);
237  assign_from_sparse(rowind, colind, values,
238  rowPermutation, colPermutation);
239  }
247  inline void add_values
248  (std::vector<int> const & rowind,
249  std::vector<int> const & colind,
250  std::vector<Treal> const & values) {
251  this->matrixPtr->syAddValues(rowind, colind, values);
252  }
253 
257  inline void add_values
258  (std::vector<int> const & rowind,
259  std::vector<int> const & colind,
260  std::vector<Treal> const & values,
261  std::vector<int> const & rowPermutation,
262  std::vector<int> const & colPermutation) {
263  std::vector<int> newRowind;
264  std::vector<int> newColind;
265  this->getPermutedAndSymmetrized(rowind, rowPermutation, newRowind,
266  colind, colPermutation, newColind);
267  this->matrixPtr->syAddValues(newRowind, newColind, values);
268  }
269 
270 
271 
272  inline void get_values
273  (std::vector<int> const & rowind,
274  std::vector<int> const & colind,
275  std::vector<Treal> & values) const {
276  this->matrixPtr->syGetValues(rowind, colind, values);
277  }
285  inline void get_values
286  (std::vector<int> const & rowind,
287  std::vector<int> const & colind,
288  std::vector<Treal> & values,
289  std::vector<int> const & rowPermutation,
290  std::vector<int> const & colPermutation) const {
291  std::vector<int> newRowind;
292  std::vector<int> newColind;
293  this->getPermutedAndSymmetrized(rowind, rowPermutation, newRowind,
294  colind, colPermutation, newColind);
295  this->matrixPtr->syGetValues(newRowind, newColind, values);
296  }
302  inline void get_all_values
303  (std::vector<int> & rowind,
304  std::vector<int> & colind,
305  std::vector<Treal> & values) const {
306  rowind.resize(0);
307  colind.resize(0);
308  values.resize(0);
309  rowind.reserve(nnz());
310  colind.reserve(nnz());
311  values.reserve(nnz());
312  this->matrixPtr->syGetAllValues(rowind, colind, values);
313  }
319  inline void get_all_values
320  (std::vector<int> & rowind,
321  std::vector<int> & colind,
322  std::vector<Treal> & values,
323  std::vector<int> const & rowInversePermutation,
324  std::vector<int> const & colInversePermutation) const {
325  std::vector<int> tmpRowind;
326  std::vector<int> tmpColind;
327  tmpRowind.reserve(rowind.capacity());
328  tmpColind.reserve(colind.capacity());
329  values.resize(0);
330  this->matrixPtr->syGetAllValues(tmpRowind, tmpColind, values);
331  this->getPermutedAndSymmetrized(tmpRowind, rowInversePermutation, rowind,
332  tmpColind, colInversePermutation, colind);
333  }
347  return *this;
348  }
352  this->matrixPtr->nosymToSym();
353  return *this;
354  }
356  *this->matrixPtr = k;
357  return *this;
358  }
359 
360  inline Treal frob() const {
361  return this->matrixPtr->syFrob();
362  }
363  Treal mixed_norm(Treal const requestedAccuracy,
364  int maxIter = -1) const;
365  Treal eucl(Treal const requestedAccuracy,
366  int maxIter = -1) const;
367 
368  void quickEuclBounds(Treal & euclLowerBound,
369  Treal & euclUpperBound) const {
370  Treal frobTmp = frob();
371  euclLowerBound = frobTmp / template_blas_sqrt( (Treal)this->get_nrows() );
372  euclUpperBound = frobTmp;
373  }
374 
375 
381  static Interval<Treal>
384  normType const norm,
385  Treal const requestedAccuracy);
393  static Interval<Treal>
396  normType const norm,
397  Treal const requestedAccuracy,
398  Treal const maxAbsVal);
402  static inline Treal frob_diff
405  return Tmatrix::syFrobDiff(*A.matrixPtr, *B.matrixPtr);
406  }
407 
411  static Treal eucl_diff
414  Treal const requestedAccuracy,
415  int maxIter = -1);
416 
420  static Treal mixed_diff
423  Treal const requestedAccuracy);
424 
434  Treal const requestedAccuracy,
435  Treal const maxAbsVal,
436  VectorType * const eVecPtr = 0);
437 
438 
449  Treal thresh(Treal const threshold,
450  normType const norm);
451 
458  inline Treal frob_thresh(Treal const threshold) {
459  return 2.0 * this->matrixPtr->frob_thresh(threshold / 2);
460  }
461 
462  Treal eucl_thresh(Treal const threshold,
463  MatrixTriangular<Treal, Tmatrix> const * const Zptr = NULL);
464 
465  Treal eucl_element_level_thresh(Treal const threshold);
466 
468  ( SizesAndBlocks & rows_new, SizesAndBlocks & cols_new ) const;
469 
470  Treal mixed_norm_thresh(Treal const threshold);
471 
472  void simple_blockwise_frob_thresh(Treal const threshold) {
473  this->matrixPtr->frobThreshLowestLevel(threshold*threshold, 0);
474  }
475 
476  inline void gershgorin(Treal& lmin, Treal& lmax) const {
477  this->matrixPtr->sy_gershgorin(lmin, lmax);
478  }
479  static inline Treal trace_ab
482  return Tmatrix::sy_trace_ab(*A.matrixPtr, *B.matrixPtr);
483  }
484  inline size_t nnz() const { /* Note: size_t instead of int here to avoid integer overflow. */
485  return this->matrixPtr->sy_nnz();
486  }
487  inline size_t nvalues() const { /* Note: size_t instead of int here to avoid integer overflow. */
488  return this->matrixPtr->sy_nvalues();
489  }
490  inline void write_to_buffer(void* buffer, const int n_bytes) const {
491  this->write_to_buffer_base(buffer, n_bytes, matrix_symm);
492  }
493  inline void read_from_buffer(void* buffer, const int n_bytes) {
494  this->read_from_buffer_base(buffer, n_bytes, matrix_symm);
495  }
496 
497 
503  (const XYZpUV<Treal,
506  Treal,
510  (const XYZ<Treal,
515  (const XYZ<Treal,
520  (const XYZpUV<Treal,
523  Treal,
527  (const XYZ<Treal,
532  (const XYZ<Treal,
535 
543 
548  static void ssmmUpperTriangleOnly(const Treal alpha,
551  const Treal beta,
553 
554 
555  /* Addition */
559  MatrixSymmetric<Treal, Tmatrix> > const & mpm);
563  MatrixSymmetric<Treal, Tmatrix> > const & mm);
570 
574 
578 
579  template<typename Top>
580  Treal accumulateWith(Top & op) {
581  return this->matrixPtr->syAccumulateWith(op);
582  }
583 
584  void random() {
585  this->matrixPtr->syRandom();
586  }
587 
588  void randomZeroStructure(Treal probabilityBeingZero) {
589  this->matrixPtr->syRandomZeroStructure(probabilityBeingZero);
590  }
591 
596  template<typename TRule>
597  void setElementsByRule(TRule & rule) {
598  this->matrixPtr->sySetElementsByRule(rule);
599  return;
600  }
601 
605  ValidPtr<Tmatrix>::swap( this->matrixPtr, dest.matrixPtr );
606  // *this now contains previous content of dest
607  this->clear();
608  }
609 
610  template<typename Tvector>
611  void matVecProd(Tvector & y, Tvector const & x) const {
612  Treal const ONE = 1.0;
613  y = (ONE * (*this) * x);
614  }
615 
616 
617  std::string obj_type_id() const {return "MatrixSymmetric";}
618  protected:
619  inline void writeToFileProt(std::ofstream & file) const {
620  this->writeToFileBase(file, matrix_symm);
621  }
622  inline void readFromFileProt(std::ifstream & file) {
623  this->readFromFileBase(file, matrix_symm);
624  }
625 
632  (std::vector<int> const & rowind,
633  std::vector<int> const & rowPermutation,
634  std::vector<int> & newRowind,
635  std::vector<int> const & colind,
636  std::vector<int> const & colPermutation,
637  std::vector<int> & newColind) {
639  getPermutedIndexes(rowind, rowPermutation, newRowind);
641  getPermutedIndexes(colind, colPermutation, newColind);
642  int tmp;
643  for (unsigned int i = 0; i < newRowind.size(); ++i) {
644  if (newRowind[i] > newColind[i]) {
645  tmp = newRowind[i];
646  newRowind[i] = newColind[i];
647  newColind[i] = tmp;
648  }
649  }
650  }
651  private:
652  };
653 
654  template<typename Treal, typename Tmatrix>
656  mixed_norm(Treal const requestedAccuracy,
657  int maxIter) const {
658  // Construct SizesAndBlocks for frobNormMat
659  SizesAndBlocks rows_new;
660  SizesAndBlocks cols_new;
661  this->getSizesAndBlocksForFrobNormMat( rows_new, cols_new );
662  // Now we can construct an empty matrix where the Frobenius norms
663  // of lowest level nonzero submatrices will be stored
665  frobNormMat.resetSizesAndBlocks(rows_new, cols_new);
666  frobNormMat.getMatrix().syAssignFrobNormsLowestLevel(this->getMatrix());
667  return frobNormMat.eucl(requestedAccuracy, maxIter);
668  }
669 
670 
671  template<typename Treal, typename Tmatrix>
673  eucl(Treal const requestedAccuracy,
674  int maxIter) const {
675  assert(requestedAccuracy >= 0);
676  /* Check if norm is really small, in that case quick return */
677  Treal frobTmp = this->frob();
678  if (frobTmp < requestedAccuracy)
679  return (Treal)0.0;
680  if (maxIter < 0)
681  maxIter = this->get_nrows() * 100;
682  VectorType guess;
684  this->getCols(cols);
685  guess.resetSizesAndBlocks(cols);
686  guess.rand();
687  // Elias note 2010-03-26: changed this back from "new code" to "old code" to reduce memory usage.
688 #if 0 // "new code"
690  Copy.frob_thresh(requestedAccuracy / 2.0);
693  lan(Copy, guess, maxIter);
694  lan.setAbsTol( requestedAccuracy / 2.0 );
695 #else // "old code"
698  lan(*this, guess, maxIter);
699  lan.setAbsTol( requestedAccuracy );
700 #endif
701  lan.run();
702  Treal eVal = 0;
703  Treal acc = 0;
704  lan.getLargestMagnitudeEig(eVal, acc);
705  return template_blas_fabs(eVal);
706  }
707 
708  template<typename Treal, typename Tmatrix>
712  normType const norm, Treal const requestedAccuracy) {
713  Treal diff;
714  Treal eNMin;
715  switch (norm) {
716  case frobNorm:
717  diff = frob_diff(A, B);
718  return Interval<Treal>(diff / template_blas_sqrt((Treal)A.get_nrows()), diff);
719  break;
720  case euclNorm:
721  diff = eucl_diff(A, B, requestedAccuracy);
722  eNMin = diff - requestedAccuracy;
723  eNMin = eNMin >= 0 ? eNMin : 0;
724  return Interval<Treal>(eNMin, diff + requestedAccuracy);
725  break;
726  default:
727  throw Failure("MatrixSymmetric<Treal, Tmatrix>::"
728  "diff(const MatrixSymmetric<Treal, Tmatrix>&, "
729  "const MatrixSymmetric<Treal, Tmatrix>&, "
730  "normType const, Treal): "
731  "Diff not implemented for selected norm");
732  }
733  }
734 
735 #if 1
736  template<typename Treal, typename Tmatrix>
740  normType const norm,
741  Treal const requestedAccuracy,
742  Treal const maxAbsVal) {
743  Treal diff;
744  switch (norm) {
745  case frobNorm:
746  {
747  diff = frob_diff(A, B);
748  return Interval<Treal>(diff / template_blas_sqrt((Treal)A.get_nrows()), diff);
749  }
750  break;
751  case euclNorm:
752  return euclDiffIfSmall(A, B, requestedAccuracy, maxAbsVal);
753  break;
754  default:
755  throw Failure("MatrixSymmetric<Treal, Tmatrix>::"
756  "diffIfSmall"
757  "(const MatrixSymmetric<Treal, Tmatrix>&, "
758  "const MatrixSymmetric<Treal, Tmatrix>&, "
759  "normType const, Treal const, Treal const): "
760  "Diff not implemented for selected norm");
761  }
762  }
763 
764 #endif
765 
766 
767  template<typename Treal, typename Tmatrix>
771  Treal const requestedAccuracy,
772  int maxIter) {
773  // DiffMatrix is a lightweight proxy object:
775  Treal maxAbsVal = 2 * frob_diff(A,B);
776  // Now, maxAbsVal should be larger than the Eucl norm
777  // Note that mat::euclIfSmall lies outside this class
778  Treal relTol = template_blas_sqrt(template_blas_sqrt(mat::getMachineEpsilon<Treal>()));
779  VectorType * const eVecPtrNotUsed = 0;
780  Interval<Treal> euclInt =
781  mat::euclIfSmall(Diff, requestedAccuracy, relTol, maxAbsVal, eVecPtrNotUsed, maxIter);
782  return euclInt.midPoint();
783  }
784 
785  template<typename Treal, typename Tmatrix>
789  Treal const requestedAccuracy) {
791  {
792  SizesAndBlocks rows_new;
793  SizesAndBlocks cols_new;
794  A.getSizesAndBlocksForFrobNormMat( rows_new, cols_new );
795  frobNormMat.resetSizesAndBlocks(rows_new, cols_new);
796  frobNormMat.getMatrix().syAssignDiffFrobNormsLowestLevel(A.getMatrix(),B.getMatrix());
797  }
798  return frobNormMat.eucl(requestedAccuracy);
799  }
800 
801 #if 1
802  template<typename Treal, typename Tmatrix>
806  Treal const requestedAccuracy,
807  Treal const maxAbsVal,
808  VectorType * const eVecPtr) {
809  // DiffMatrix is a lightweight proxy object:
811  // Note that this function lies outside this class
812  Treal relTol = template_blas_sqrt(template_blas_sqrt(mat::getMachineEpsilon<Treal>()));
813  Interval<Treal> tmpInterval = mat::euclIfSmall(Diff, requestedAccuracy, relTol, maxAbsVal, eVecPtr);
814  // Emanuel note: Ugly fix to make certain tests pass, we expand
815  // the interval up to the requested accuracy. Note that larger
816  // intervals may occur if the norm is not 'small'. It happens that
817  // Lanczos misconverges to for example the second largest
818  // eigenvalue. This happens in particular when the first and second
819  // eigenvalues are very close (of the order of the requested
820  // accuracy). Expanding the interval makes the largest eigenvalue
821  // (at least for certain cases) end up inside the interval even
822  // though Lanczos has misconverged.
823  if ( tmpInterval.length() < 2*requestedAccuracy )
824  return Interval<Treal>( tmpInterval.midPoint()-requestedAccuracy, tmpInterval.midPoint()+requestedAccuracy );
825  return tmpInterval;
826  }
827 
828 #endif
829 
830 
831  template<typename Treal, typename Tmatrix>
833  thresh(Treal const threshold,
834  normType const norm) {
835  switch (norm) {
836  case frobNorm:
837  return this->frob_thresh(threshold);
838  break;
839  case euclNorm:
840  return this->eucl_thresh(threshold);
841  break;
842  case mixedNorm:
843  return this->mixed_norm_thresh(threshold);
844  break;
845  default:
846  throw Failure("MatrixSymmetric<Treal, Tmatrix>::"
847  "thresh(Treal const, "
848  "normType const): "
849  "Thresholding not implemented for selected norm");
850  }
851  }
852 
853 #if 1
854 
855  template<typename Treal, typename Tmatrix>
857  eucl_thresh(Treal const threshold,
858  MatrixTriangular<Treal, Tmatrix> const * const Zptr) {
859  if ( Zptr == NULL ) {
861  return TruncObj.run( threshold );
862  }
864  return TruncObj.run( threshold );
865  }
866 
867 #endif
868 
869 
870  template<typename Treal, typename Tmatrix>
872  ( SizesAndBlocks & rows_new, SizesAndBlocks & cols_new ) const {
873  std::vector<int> rows_block_sizes;
874  std::vector<int> cols_block_sizes;
875  int n_rows;
876  int n_cols;
877  {
880  this->getRows(rows);
881  this->getCols(cols);
882  rows.getBlockSizeVector( rows_block_sizes );
883  cols.getBlockSizeVector( cols_block_sizes );
884  rows_block_sizes.pop_back(); // Remove the '1' at the end
885  cols_block_sizes.pop_back(); // Remove the '1' at the end
886  n_rows = rows.getNTotalScalars();
887  n_cols = cols.getNTotalScalars();
888  int factor_rows = rows_block_sizes[rows_block_sizes.size()-1];
889  int factor_cols = cols_block_sizes[cols_block_sizes.size()-1];
890  for (unsigned int ind = 0; ind < rows_block_sizes.size(); ++ind)
891  rows_block_sizes[ind] = rows_block_sizes[ind] / factor_rows;
892  for (unsigned int ind = 0; ind < cols_block_sizes.size(); ++ind)
893  cols_block_sizes[ind] = cols_block_sizes[ind] / factor_cols;
894  // Now set the number of (scalar) rows and cols, should be equal
895  // to the number of blocks at the lowest level of the original
896  // matrix
897  if (n_rows % factor_rows)
898  n_rows = n_rows / factor_rows + 1;
899  else
900  n_rows = n_rows / factor_rows;
901  if (n_cols % factor_cols)
902  n_cols = n_cols / factor_cols + 1;
903  else
904  n_cols = n_cols / factor_cols;
905  }
906  rows_new = SizesAndBlocks( rows_block_sizes, n_rows );
907  cols_new = SizesAndBlocks( cols_block_sizes, n_cols );
908  }
909 
910 
911  template<typename Treal, typename Tmatrix>
913  mixed_norm_thresh(Treal const threshold) {
914  assert(threshold >= (Treal)0.0);
915  if (threshold == (Treal)0.0)
916  return (Treal)0;
917  // Construct SizesAndBlocks for frobNormMat
918  SizesAndBlocks rows_new;
919  SizesAndBlocks cols_new;
920  this->getSizesAndBlocksForFrobNormMat( rows_new, cols_new );
921  // Now we can construct an empty matrix where the Frobenius norms
922  // of lowest level nonzero submatrices will be stored
924  frobNormMat.resetSizesAndBlocks(rows_new, cols_new);
925  // We want the following step to dominate the mixed_norm truncation (this
926  // is where Frobenius norms of submatrices are computed, i.e. it
927  // is here we loop over all matrix elements.)
928  frobNormMat.getMatrix().syAssignFrobNormsLowestLevel(this->getMatrix());
930  Treal mixed_norm_result = TruncObj.run( threshold );
931  frobNormMat.getMatrix().truncateAccordingToSparsityPattern(this->getMatrix());
932  return mixed_norm_result;
933  }
934 
935 
936  /* B = alpha * A */
937  template<typename Treal, typename Tmatrix>
941 
942  if(this == &sm.B) // A = B
943  {
944  *this *= sm.A; // B *= alpha
945  return *this;
946  }
947  assert(!sm.tB);
948  /* Data structure set by assign - therefore set haveDataStructure to true */
949  this->matrixPtr.haveDataStructureSet(true);
950  this->matrixPtr->assign(sm.A, *sm.B.matrixPtr);
951  return *this;
952  }
953  /* C = alpha * A * A + beta * C : A and C are symmetric */
954  template<typename Treal, typename Tmatrix>
957  (const XYZpUV<Treal,
960  Treal,
962  assert(this != &sm2psm.B);
963  if (this == &sm2psm.E && &sm2psm.B == &sm2psm.C) {
964  /* Operation is C = alpha * A * A + beta * C */
965  Tmatrix::sysq('U',
966  sm2psm.A, *sm2psm.B.matrixPtr,
967  sm2psm.D, *this->matrixPtr);
968  return *this;
969  }
970  else /* this != &sm2psm.C || &sm2psm.B != &sm2psm.C */
971  throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator="
972  "(const XYZpUV<Treal, MatrixSymmetric"
973  "<Treal, Tmatrix> >& sm2psm) : "
974  "D = alpha * A * B + beta * C not supported for C != D"
975  " and for symmetric matrices not for A != B since this "
976  "generally will result in a nonsymmetric matrix");
977  }
978 
979  /* C = alpha * A * A : A and C are symmetric */
980  template<typename Treal, typename Tmatrix>
983  (const XYZ<Treal,
986  assert(this != &sm2.B);
987  if (&sm2.B == &sm2.C) {
988  this->matrixPtr.haveDataStructureSet(true);
989  Tmatrix::sysq('U', sm2.A, *sm2.B.matrixPtr, 0, *this->matrixPtr);
990  return *this;
991  }
992  else
993  throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator="
994  "(const XYZ<Treal, MatrixSymmetric<Treal, Tmatrix>,"
995  " MatrixSymmetric<Treal, Tmatrix> >& sm2) : "
996  "Operation C = alpha * A * B with only symmetric "
997  "matrices not supported for A != B");
998  }
999 
1000  /* C += alpha * A * A : A and C are symmetric */
1001  template<typename Treal, typename Tmatrix>
1004  (const XYZ<Treal,
1007  assert(this != &sm2.B);
1008  if (&sm2.B == &sm2.C) {
1009  Tmatrix::sysq('U', sm2.A, *sm2.B.matrixPtr, 1, *this->matrixPtr);
1010  return *this;
1011  }
1012  else
1013  throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator+="
1014  "(const XYZ<Treal, MatrixSymmetric<Treal, Tmatrix>,"
1015  " MatrixSymmetric<Treal, Tmatrix> >& sm2) : "
1016  "Operation C += alpha * A * B with only symmetric "
1017  "matrices not supported for A != B");
1018  }
1019 
1020  /* C = alpha * A * transpose(A) + beta * C : C is symmetric */
1021  template<typename Treal, typename Tmatrix>
1024  (const XYZpUV<Treal,
1027  Treal,
1028  MatrixSymmetric<Treal, Tmatrix> >& smmpsm) {
1029  if (this == &smmpsm.E)
1030  if (&smmpsm.B == &smmpsm.C)
1031  if (!smmpsm.tB && smmpsm.tC) {
1032  Tmatrix::syrk('U', false,
1033  smmpsm.A, *smmpsm.B.matrixPtr,
1034  smmpsm.D, *this->matrixPtr);
1035  }
1036  else
1037  throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator="
1038  "(const XYZpUV<Treal, MatrixGeneral"
1039  "<Treal, Tmatrix>, "
1040  "MatrixGeneral<Treal, Tmatrix>, Treal, "
1041  "MatrixSymmetric<Treal, Tmatrix> >&) : "
1042  "C = alpha * A' * A + beta * C, not implemented"
1043  " only C = alpha * A * A' + beta * C");
1044  else
1045  throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator="
1046  "(const XYZpUV<"
1047  "Treal, MatrixGeneral<Treal, Tmatrix>, "
1048  "MatrixGeneral<Treal, Tmatrix>, Treal, "
1049  "MatrixSymmetric<Treal, Tmatrix> >&) : "
1050  "You are trying to call C = alpha * A * A' + beta * C"
1051  " with A and A' being different objects");
1052  else
1053  throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator="
1054  "(const XYZpUV<"
1055  "Treal, MatrixGeneral<Treal, Tmatrix>, "
1056  "MatrixGeneral<Treal, Tmatrix>, Treal, "
1057  "MatrixSymmetric<Treal, Tmatrix> >&) : "
1058  "D = alpha * A * A' + beta * C not supported for C != D");
1059  return *this;
1060  }
1061 
1062  /* C = alpha * A * transpose(A) : C is symmetric */
1063  template<typename Treal, typename Tmatrix>
1066  (const XYZ<Treal,
1069  if (&smm.B == &smm.C)
1070  if (!smm.tB && smm.tC) {
1071  Tmatrix::syrk('U', false,
1072  smm.A, *smm.B.matrixPtr,
1073  0, *this->matrixPtr);
1074  }
1075  else
1076  throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator="
1077  "(const XYZ<"
1078  "Treal, MatrixGeneral<Treal, Tmatrix>, "
1079  "MatrixGeneral<Treal, Tmatrix> >&) : "
1080  "C = alpha * A' * A, not implemented "
1081  "only C = alpha * A * A'");
1082  else
1083  throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator="
1084  "(const XYZ<"
1085  "Treal, MatrixGeneral<Treal, Tmatrix>, "
1086  "MatrixGeneral<Treal, Tmatrix> >&) : "
1087  "You are trying to call C = alpha * A * A' "
1088  "with A and A' being different objects");
1089  return *this;
1090  }
1091 
1092  /* C += alpha * A * transpose(A) : C is symmetric */
1093  template<typename Treal, typename Tmatrix>
1096  (const XYZ<Treal,
1099  if (&smm.B == &smm.C)
1100  if (!smm.tB && smm.tC) {
1101  Tmatrix::syrk('U', false,
1102  smm.A, *smm.B.matrixPtr,
1103  1, *this->matrixPtr);
1104  }
1105  else
1106  throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator+="
1107  "(const XYZ<"
1108  "Treal, MatrixGeneral<Treal, Tmatrix>, "
1109  "MatrixGeneral<Treal, Tmatrix> >&) : "
1110  "C += alpha * A' * A, not implemented "
1111  "only C += alpha * A * A'");
1112  else
1113  throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator+="
1114  "(const XYZ<"
1115  "Treal, MatrixGeneral<Treal, Tmatrix>, "
1116  "MatrixGeneral<Treal, Tmatrix> >&) : "
1117  "You are trying to call C += alpha * A * A' "
1118  "with A and A' being different objects");
1119  return *this;
1120  }
1121 
1122 #if 1
1123  /* A = op1(Z) * A * op2(Z) : Z is upper triangular and A is symmetric */
1124  /* Either op1() or op2() is the transpose operation. */
1125  template<typename Treal, typename Tmatrix>
1131  if (this == &zaz.B) {
1132  if (&zaz.A == &zaz.C) {
1133  if (zaz.tA && !zaz.tC) {
1134  Tmatrix::trsytriplemm('R', *zaz.A.matrixPtr, *this->matrixPtr);
1135  }
1136  else if (!zaz.tA && zaz.tC) {
1137  Tmatrix::trsytriplemm('L', *zaz.A.matrixPtr, *this->matrixPtr);
1138  }
1139  else
1140  throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator="
1141  "(const XYZ<MatrixTriangular<Treal, Tmatrix>,"
1142  "MatrixSymmetric<Treal, Tmatrix>,"
1143  "MatrixTriangular<Treal, Tmatrix> >&) : "
1144  "A = op1(Z) * A * op2(Z) : Either op1 xor op2 must be "
1145  "the transpose operation.");
1146  }
1147  else
1148  throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator="
1149  "(const XYZ<MatrixTriangular<Treal, Tmatrix>,"
1150  "MatrixSymmetric<Treal, Tmatrix>,"
1151  "MatrixTriangular<Treal, Tmatrix> >&) : "
1152  "A = op1(Z1) * A * op2(Z2) : Z1 and Z2 must be the same "
1153  "object");
1154  }
1155  else
1156  throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator="
1157  "(const XYZ<MatrixTriangular<Treal, Tmatrix>,"
1158  "MatrixSymmetric<Treal, Tmatrix>,"
1159  "MatrixTriangular<Treal, Tmatrix> >&) : "
1160  "C = op1(Z) * A * op2(Z) : A and C must be the same "
1161  "object");
1162  return *this;
1163  }
1164 
1165 #endif
1166 
1167 
1172  template<typename Treal, typename Tmatrix>
1174  ssmmUpperTriangleOnly(const Treal alpha,
1177  const Treal beta,
1179  Tmatrix::ssmm_upper_tr_only(alpha, *A.matrixPtr, *B.matrixPtr,
1180  beta, *C.matrixPtr);
1181  }
1182 
1183 
1184 
1185 
1186 
1187  /* Addition */
1188  /* C = A + B */
1189  template<typename Treal, typename Tmatrix>
1194  assert(this != &mpm.A);
1195  (*this) = mpm.B;
1196  Tmatrix::add(1.0, *mpm.A.matrixPtr, *this->matrixPtr);
1197  return *this;
1198  }
1199  /* C = A - B */
1200  template<typename Treal, typename Tmatrix>
1205  assert(this != &mmm.B);
1206  (*this) = mmm.A;
1207  Tmatrix::add(-1.0, *mmm.B.matrixPtr, *this->matrixPtr);
1208  return *this;
1209  }
1210  /* B += A */
1211  template<typename Treal, typename Tmatrix>
1215  Tmatrix::add(1.0, *A.matrixPtr, *this->matrixPtr);
1216  return *this;
1217  }
1218  /* B -= A */
1219  template<typename Treal, typename Tmatrix>
1223  Tmatrix::add(-1.0, *A.matrixPtr, *this->matrixPtr);
1224  return *this;
1225  }
1226  /* B += alpha * A */
1227  template<typename Treal, typename Tmatrix>
1231  assert(!sm.tB);
1232  Tmatrix::add(sm.A, *sm.B.matrixPtr, *this->matrixPtr);
1233  return *this;
1234  }
1235 
1236  /* B -= alpha * A */
1237  template<typename Treal, typename Tmatrix>
1241  assert(!sm.tB);
1242  Tmatrix::add(-sm.A, *sm.B.matrixPtr, *this->matrixPtr);
1243  return *this;
1244  }
1245 
1250  template<typename Treal, typename Tmatrix, typename Top>
1252  Top & op) {
1253  return A.accumulateWith(op);
1254  }
1255 
1256 
1257 
1258 } /* end namespace mat */
1259 #endif
1260 
mat::euclNorm
@ euclNorm
Definition: matInclude.h:139
mat::MatrixSymmetric::randomZeroStructure
void randomZeroStructure(Treal probabilityBeingZero)
Definition: MatrixSymmetric.h:588
template_blas_sqrt
Treal template_blas_sqrt(Treal x)
mat::mixedNorm
@ mixedNorm
Definition: matInclude.h:139
mat::MatrixSymmetric::MatrixSymmetric
MatrixSymmetric(const XY< Treal, MatrixSymmetric< Treal, Tmatrix > > &sm)
Definition: MatrixSymmetric.h:83
mat::MatrixSymmetric::nnz
size_t nnz() const
Definition: MatrixSymmetric.h:484
mat::XpY
This proxy expresses the result of addition of two objects, of possibly different types,...
Definition: matrix_proxy.h:246
mat::MatrixSymmetric::assign_from_sparse
void assign_from_sparse(std::vector< int > const &rowind, std::vector< int > const &colind, std::vector< Treal > const &values, SizesAndBlocks const &newRows, SizesAndBlocks const &newCols)
Assign from sparse matrix given by three vectors.
Definition: MatrixSymmetric.h:193
mat::MatrixSymmetric::writeToFileProt
void writeToFileProt(std::ofstream &file) const
Definition: MatrixSymmetric.h:619
mat::MatrixSymmetric::readFromFileProt
void readFromFileProt(std::ifstream &file)
Definition: MatrixSymmetric.h:622
mat::Interval::midPoint
Treal midPoint() const
Definition: Interval.h:115
mat::MatrixSymmetric::frob_diff
static Treal frob_diff(const MatrixSymmetric< Treal, Tmatrix > &A, const MatrixSymmetric< Treal, Tmatrix > &B)
Returns the Frobenius norm of A - B ( || A - B ||_F )
Definition: MatrixSymmetric.h:403
mat::MatrixSymmetric::nvalues
size_t nvalues() const
Definition: MatrixSymmetric.h:487
truncation.h
mat_utils.h
Utilities used by other parts of the hierarchical matrix library.
mat::MatrixSymmetric::gershgorin
void gershgorin(Treal &lmin, Treal &lmax) const
Definition: MatrixSymmetric.h:476
MatrixBase.h
mat::MatrixSymmetric::assign_from_sparse
void assign_from_sparse(std::vector< int > const &rowind, std::vector< int > const &colind, std::vector< Treal > const &values)
Assign from sparse matrix given by three vectors.
Definition: MatrixSymmetric.h:176
mat::MatrixSymmetric::trace_ab
static Treal trace_ab(const MatrixSymmetric< Treal, Tmatrix > &A, const MatrixSymmetric< Treal, Tmatrix > &B)
Definition: MatrixSymmetric.h:480
mat::MatrixSymmetric::get_all_values
void get_all_values(std::vector< int > &rowind, std::vector< int > &colind, std::vector< Treal > &values) const
Get all values and corresponding row and column index lists, in matrix.
Definition: MatrixSymmetric.h:303
mat::MatrixSymmetric::assign_from_sparse
void assign_from_sparse(std::vector< int > const &rowind, std::vector< int > const &colind, std::vector< Treal > const &values, SizesAndBlocks const &newRows, SizesAndBlocks const &newCols, std::vector< int > const &rowPermutation, std::vector< int > const &colPermutation)
Same as above, except taking sizes and blocks arguments.
Definition: MatrixSymmetric.h:229
mat::MatrixSymmetric::transfer
void transfer(MatrixSymmetric< Treal, Tmatrix > &dest)
Transfer this matrix to dest, clearing previous content of dest if any.
Definition: MatrixSymmetric.h:604
mat::MatrixSymmetric::obj_type_id
std::string obj_type_id() const
Definition: MatrixSymmetric.h:617
mat::MatrixSymmetric::fullMatrix
void fullMatrix(std::vector< Treal > &fullMat, std::vector< int > const &rowInversePermutation, std::vector< int > const &colInversePermutation) const
Save matrix as full matrix.
Definition: MatrixSymmetric.h:148
mat::EuclTruncationSymm
Truncation of symmetric matrices.
Definition: truncation.h:154
mat::MatrixSymmetric::getSizesAndBlocksForFrobNormMat
void getSizesAndBlocksForFrobNormMat(SizesAndBlocks &rows_new, SizesAndBlocks &cols_new) const
Definition: MatrixSymmetric.h:872
mat::frobNorm
@ frobNorm
Definition: matInclude.h:139
mat::MatrixSymmetric::matVecProd
void matVecProd(Tvector &y, Tvector const &x) const
Definition: MatrixSymmetric.h:611
mat::MatrixSymmetric::write_to_buffer
void write_to_buffer(void *buffer, const int n_bytes) const
Definition: MatrixSymmetric.h:490
rows
mat::SizesAndBlocks rows
Definition: test.cc:51
mat::MatrixSymmetric::diff
static Interval< Treal > diff(const MatrixSymmetric< Treal, Tmatrix > &A, const MatrixSymmetric< Treal, Tmatrix > &B, normType const norm, Treal const requestedAccuracy)
Returns interval containing the Euclidean norm of A - B ( || A - B ||_2 )
Definition: MatrixSymmetric.h:710
mat::XYZ
This proxy expresses the result of multiplication of three objects, of possibly different types,...
Definition: matrix_proxy.h:67
mat::MatrixSymmetric::frob
Treal frob() const
Definition: MatrixSymmetric.h:360
mat::SizesAndBlocks::getBlockSizeVector
void getBlockSizeVector(std::vector< int > &blockSizesCopy) const
Definition: SizesAndBlocks.cc:87
template_blas_fabs
Treal template_blas_fabs(Treal x)
mat::MatrixSymmetric::simple_blockwise_frob_thresh
void simple_blockwise_frob_thresh(Treal const threshold)
Definition: MatrixSymmetric.h:472
mat::MatrixSymmetric::getPermutedAndSymmetrized
static void getPermutedAndSymmetrized(std::vector< int > const &rowind, std::vector< int > const &rowPermutation, std::vector< int > &newRowind, std::vector< int > const &colind, std::vector< int > const &colPermutation, std::vector< int > &newColind)
This function permutes row and column indices according to the specified permutation and gives the in...
Definition: MatrixSymmetric.h:632
B
#define B
mat::EuclTruncationSymmElementLevel
Truncation of symmetric matrices at the element level (used for mixed norm truncation)
Definition: truncation.h:246
LanczosLargestMagnitudeEig.h
mat::accumulate
Treal accumulate(MatrixSymmetric< Treal, Tmatrix > &A, Top &op)
Performs operation specified in 'op' on all nonzero matrix elements and sums up the result and return...
Definition: MatrixSymmetric.h:1251
mat::MatrixSymmetric::setElementsByRule
void setElementsByRule(TRule &rule)
Uses rule depending on the row and column indexes to set matrix elements The Trule class should have ...
Definition: MatrixSymmetric.h:597
mat::MatrixSymmetric
Symmetric matrix.
Definition: MatrixBase.h:51
mat::MatrixTriangular
Upper non-unit triangular matrix.
Definition: MatrixBase.h:53
mat::MatrixSymmetric::thresh
Treal thresh(Treal const threshold, normType const norm)
Does thresholding so that the error in the chosen norm is below the given threshold.
Definition: MatrixSymmetric.h:833
mat::VectorGeneral::resetSizesAndBlocks
void resetSizesAndBlocks(SizesAndBlocks const &newRows)
Definition: VectorGeneral.h:51
mat::Interval::length
Treal length() const
Returns the length of the interval.
Definition: Interval.h:109
mat::MatrixSymmetric::eucl
Treal eucl(Treal const requestedAccuracy, int maxIter=-1) const
Definition: MatrixSymmetric.h:673
mat::arn::LanczosLargestMagnitudeEig
Definition: LanczosLargestMagnitudeEig.h:47
mat::MatrixSymmetric::eucl_element_level_thresh
Treal eucl_element_level_thresh(Treal const threshold)
mat::MatrixSymmetric::eucl_thresh
Treal eucl_thresh(Treal const threshold, MatrixTriangular< Treal, Tmatrix > const *const Zptr=NULL)
Definition: MatrixSymmetric.h:857
mat::Failure
Definition: Failure.h:57
mat::arn::LanczosLargestMagnitudeEig::getLargestMagnitudeEig
void getLargestMagnitudeEig(Treal &ev, Treal &accuracy)
Definition: LanczosLargestMagnitudeEig.h:58
mat::MatrixSymmetric::assignFromFull
void assignFromFull(std::vector< Treal > const &fullMat, std::vector< int > const &rowPermutation, std::vector< int > const &colPermutation)
Definition: MatrixSymmetric.h:117
mat::MatrixSymmetric::mixed_norm_thresh
Treal mixed_norm_thresh(Treal const threshold)
Definition: MatrixSymmetric.h:913
mat::matrix_symm
@ matrix_symm
Definition: MatrixBase.h:56
mat::MatrixSymmetric::VectorType
VectorGeneral< Treal, typename Tmatrix::VectorType > VectorType
Definition: MatrixSymmetric.h:70
mat::ValidPtr::swap
static void swap(ValidPtr< Tobj > &ptrA, ValidPtr< Tobj > &ptrB)
Definition: ValidPtr.h:106
mat::XmY
This proxy expresses the result of substraction of two objects, of possibly different types,...
Definition: matrix_proxy.h:266
mat::MatrixSymmetric::real
Treal real
Definition: MatrixSymmetric.h:71
mat::MatrixSymmetric::operator=
MatrixSymmetric< Treal, Tmatrix > & operator=(const MatrixSymmetric< Treal, Tmatrix > &symm)
Definition: MatrixSymmetric.h:345
mat::EuclTruncationSymmWithZ
Truncation of symmetric matrices with Z.
Definition: truncation.h:203
mat::MatrixSymmetric::frob_thresh
Treal frob_thresh(Treal const threshold)
Does thresholding so that the error in the Frobenius norm is below the given threshold.
Definition: MatrixSymmetric.h:458
mat::VectorGeneral
Definition: MatrixBase.h:55
mat::euclIfSmall
Interval< Treal > euclIfSmall(Tmatrix const &M, Treal const requestedAbsAccuracy, Treal const requestedRelAccuracy, Treal const maxAbsVal, typename Tmatrix::VectorType *const eVecPtr=0, int maxIter=-1)
Returns interval containing the Euclidean norm of the matrix M.
Definition: LanczosLargestMagnitudeEig.h:260
mat::MatrixSymmetric::mixed_norm
Treal mixed_norm(Treal const requestedAccuracy, int maxIter=-1) const
Definition: MatrixSymmetric.h:656
mat::MatrixGeneral
Normal matrix.
Definition: MatrixBase.h:49
mat::MatrixSymmetric::random
void random()
Definition: MatrixSymmetric.h:584
mat::MatrixSymmetric::add_values
void add_values(std::vector< int > const &rowind, std::vector< int > const &colind, std::vector< Treal > const &values, std::vector< int > const &rowPermutation, std::vector< int > const &colPermutation)
Same as above, except taking two additional arguments specifying the permutation of rows and columns.
Definition: MatrixSymmetric.h:258
mat::MatrixSymmetric::add_values
void add_values(std::vector< int > const &rowind, std::vector< int > const &colind, std::vector< Treal > const &values)
Add given set of values to the matrix.
Definition: MatrixSymmetric.h:248
cols
mat::SizesAndBlocks cols
Definition: test.cc:52
mat
Definition: allocate.cc:39
A
#define A
mat::EuclTruncationBase::run
Treal run(Treal const threshold)
Definition: truncation.h:80
mat::MatrixSymmetric::accumulateWith
Treal accumulateWith(Top &op)
Definition: MatrixSymmetric.h:580
mat::normType
normType
Definition: matInclude.h:139
mat::XYZpUV
This proxy expresses the result of multiplication of three objects added to two other multiplied obje...
Definition: matrix_proxy.h:88
mat::Interval
Definition: Interval.h:46
mat::MatrixSymmetric::operator=
MatrixSymmetric< Treal, Tmatrix > & operator=(const MatrixGeneral< Treal, Tmatrix > &matr)
Definition: MatrixSymmetric.h:350
mat::symm
static void symm(const char *side, const char *uplo, const int *m, const int *n, const T *alpha, const T *A, const int *lda, const T *B, const int *ldb, const T *beta, T *C, const int *ldc)
Definition: mat_gblas.h:342
mat::VectorGeneral::rand
void rand()
Definition: VectorGeneral.h:108
mat::MatrixSymmetric::operator=
MatrixSymmetric< Treal, Tmatrix > & operator=(int const k)
Definition: MatrixSymmetric.h:355
mat::MatrixSymmetric::MatrixSymmetric
MatrixSymmetric(const MatrixGeneral< Treal, Tmatrix > &matr)
'Copy from normal matrix' - constructor
Definition: MatrixSymmetric.h:85
mat::MatrixBase
Base class for matrix API.
Definition: MatrixBase.h:69
mat::MatrixSymmetric::get_values
void get_values(std::vector< int > const &rowind, std::vector< int > const &colind, std::vector< Treal > &values) const
Get values given by row and column index lists.
Definition: MatrixSymmetric.h:273
mat::MatrixSymmetric::diffIfSmall
static Interval< Treal > diffIfSmall(const MatrixSymmetric< Treal, Tmatrix > &A, const MatrixSymmetric< Treal, Tmatrix > &B, normType const norm, Treal const requestedAccuracy, Treal const maxAbsVal)
Returns interval containing the Euclidean norm of A - B ( || A - B ||_2 ) based on the chosen norm.
Definition: MatrixSymmetric.h:738
mat::MatrixSymmetric::MatrixSymmetric
MatrixSymmetric(const MatrixSymmetric< Treal, Tmatrix > &symm)
Copy constructor
Definition: MatrixSymmetric.h:80
mat::XY
This proxy expresses the result of multiplication of two objects, of possibly different types,...
Definition: matrix_proxy.h:51
mat::MatrixSymmetric::quickEuclBounds
void quickEuclBounds(Treal &euclLowerBound, Treal &euclUpperBound) const
Definition: MatrixSymmetric.h:368
mat::SizesAndBlocks
Describes dimensions of matrix and its blocks on all levels.
Definition: SizesAndBlocks.h:45
mat::syrk
static void syrk(const char *uplo, const char *trans, const int *n, const int *k, const T *alpha, const T *A, const int *lda, const T *beta, T *C, const int *ldc)
Definition: mat_gblas.h:334
mat::MatrixSymmetric::euclDiffIfSmall
static Interval< Treal > euclDiffIfSmall(const MatrixSymmetric< Treal, Tmatrix > &A, const MatrixSymmetric< Treal, Tmatrix > &B, Treal const requestedAccuracy, Treal const maxAbsVal, VectorType *const eVecPtr=0)
Returns interval containing the Euclidean norm of A - B ( || A - B ||_2 ).
Definition: MatrixSymmetric.h:804
mat::MatrixBase::getPermutedIndexes
static void getPermutedIndexes(std::vector< int > const &index, std::vector< int > const &permutation, std::vector< int > &newIndex)
Definition: MatrixBase.h:205
mat::MatrixSymmetric::ssmmUpperTriangleOnly
static void ssmmUpperTriangleOnly(const Treal alpha, const MatrixSymmetric< Treal, Tmatrix > &A, const MatrixSymmetric< Treal, Tmatrix > &B, const Treal beta, MatrixSymmetric< Treal, Tmatrix > &C)
C = alpha * A * B + beta * C where A and B are symmetric and only the upper triangle of C is computed...
Definition: MatrixSymmetric.h:1174
mat::MatrixBase::operator=
MatrixBase< Treal, Tmatrix > & operator=(const MatrixBase< Treal, Tmatrix > &other)
Definition: MatrixBase.h:166
mat::MatrixSymmetric::assignFromFull
void assignFromFull(std::vector< Treal > const &fullMat)
Definition: MatrixSymmetric.h:110
mat::MatrixSymmetric::get_all_values
void get_all_values(std::vector< int > &rowind, std::vector< int > &colind, std::vector< Treal > &values, std::vector< int > const &rowInversePermutation, std::vector< int > const &colInversePermutation) const
Same as above, except taking two additional arguments specifying the permutation of rows and columns.
Definition: MatrixSymmetric.h:320
mat::MatrixSymmetric::mixed_diff
static Treal mixed_diff(const MatrixSymmetric< Treal, Tmatrix > &A, const MatrixSymmetric< Treal, Tmatrix > &B, Treal const requestedAccuracy)
Returns the 'mixed' norm of A - B ( || A - B ||_mixed )
Definition: MatrixSymmetric.h:787
mat::DiffMatrix
Definition: mat_utils.h:43
mat::SizesAndBlocks::getNTotalScalars
int getNTotalScalars() const
Definition: SizesAndBlocks.h:84
mat::MatrixSymmetric::read_from_buffer
void read_from_buffer(void *buffer, const int n_bytes)
Definition: MatrixSymmetric.h:493
mat::MatrixSymmetric::get_values
void get_values(std::vector< int > const &rowind, std::vector< int > const &colind, std::vector< Treal > &values, std::vector< int > const &rowPermutation, std::vector< int > const &colPermutation) const
Same as above, except taking two additional arguments specifying the permutation of rows and columns.
Definition: MatrixSymmetric.h:286
mat::MatrixSymmetric::fullMatrix
void fullMatrix(std::vector< Treal > &fullMat) const
Save matrix as full matrix.
Definition: MatrixSymmetric.h:138
mat::MatrixSymmetric::eucl_diff
static Treal eucl_diff(const MatrixSymmetric< Treal, Tmatrix > &A, const MatrixSymmetric< Treal, Tmatrix > &B, Treal const requestedAccuracy, int maxIter=-1)
Returns the Euclidean norm of A - B ( || A - B ||_2 )
Definition: MatrixSymmetric.h:769
mat::MatrixSymmetric::assign_from_sparse
void assign_from_sparse(std::vector< int > const &rowind, std::vector< int > const &colind, std::vector< Treal > const &values, std::vector< int > const &rowPermutation, std::vector< int > const &colPermutation)
Same as above, except taking two additional arguments specifying the permutation of rows and columns.
Definition: MatrixSymmetric.h:211
mat::MatrixSymmetric::MatrixSymmetric
MatrixSymmetric()
Default constructor
Definition: MatrixSymmetric.h:73
Interval.h