ergo
MatrixSymmetric.h
Go to the documentation of this file.
1 /* Ergo, version 3.7, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2018 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>() {}
76  :MatrixBase<Treal, Tmatrix>(symm) {}
77  explicit MatrixSymmetric(const XY<Treal, MatrixSymmetric<Treal, Tmatrix> >& sm)
78  :MatrixBase<Treal, Tmatrix>() { *this = sm.A * sm.B; }
80  :MatrixBase<Treal, Tmatrix>(matr) {
81  this->matrixPtr->nosymToSym();
82  }
85 #if 0
86  template<typename Tfull>
87  inline void assign_from_full
88  (Tfull const* const fullmatrix,
89  int const totnrows, int const totncols) {
90  assert(totnrows == totncols);
91  this->matrixPtr->assign_from_full(fullmatrix, totnrows, totncols);
92  this->matrixPtr->nosym_to_sym();
93  }
94  inline void assign_from_full
95  (Treal const* const fullmatrix,
96  int const totnrows, int const totncols) {
97  assert(totnrows == totncols);
98  this->matrixPtr->assign_from_full(fullmatrix, totnrows, totncols);
99  this->matrixPtr->nosym_to_sym();
100  }
101 #endif
102 
103  inline void assignFromFull
104  (std::vector<Treal> const & fullMat) {
105  assert((int)fullMat.size() == this->get_nrows() * this->get_ncols());
106  this->matrixPtr->assignFromFull(fullMat);
107  this->matrixPtr->nosymToSym();
108  }
109 
110  inline void assignFromFull
111  (std::vector<Treal> const & fullMat,
112  std::vector<int> const & rowPermutation,
113  std::vector<int> const & colPermutation) {
114  assert((int)fullMat.size() == this->get_nrows() * this->get_ncols());
115  std::vector<int> rowind(fullMat.size());
116  std::vector<int> colind(fullMat.size());
117  int ind = 0;
118  for (int col = 0; col < this->get_ncols(); ++col)
119  for (int row = 0; row < this->get_nrows(); ++row) {
120  rowind[ind] = row;
121  colind[ind] = col;
122  ++ind;
123  }
124  this->assign_from_sparse(rowind,
125  colind,
126  fullMat,
127  rowPermutation,
128  colPermutation);
129  this->matrixPtr->nosymToSym();
130  }
131 
132  inline void fullMatrix(std::vector<Treal> & fullMat) const {
133  this->matrixPtr->syFullMatrix(fullMat);
134  }
141  inline void fullMatrix
142  (std::vector<Treal> & fullMat,
143  std::vector<int> const & rowInversePermutation,
144  std::vector<int> const & colInversePermutation) const {
145  std::vector<int> rowind;
146  std::vector<int> colind;
147  std::vector<Treal> values;
148  get_all_values(rowind, colind, values,
149  rowInversePermutation,
150  colInversePermutation);
151  fullMat.assign(this->get_nrows()*this->get_ncols(),0);
152  assert(rowind.size() == values.size());
153  assert(rowind.size() == colind.size());
154  for (unsigned int ind = 0; ind < values.size(); ++ind) {
155  assert(rowind[ind] + this->get_nrows() * colind[ind] <
156  this->get_nrows()*this->get_ncols());
157  fullMat[rowind[ind] + this->get_nrows() * colind[ind]] =
158  values[ind];
159  fullMat[colind[ind] + this->get_nrows() * rowind[ind]] =
160  values[ind];
161  }
162  }
169  inline void assign_from_sparse
170  (std::vector<int> const & rowind,
171  std::vector<int> const & colind,
172  std::vector<Treal> const & values) {
173  this->matrixPtr->syAssignFromSparse(rowind, colind, values);
174  }
186  inline void assign_from_sparse
187  (std::vector<int> const & rowind,
188  std::vector<int> const & colind,
189  std::vector<Treal> const & values,
190  SizesAndBlocks const & newRows,
191  SizesAndBlocks const & newCols) {
192  this->resetSizesAndBlocks(newRows, newCols);
193  this->matrixPtr->syAssignFromSparse(rowind, colind, values);
194  }
204  inline void assign_from_sparse
205  (std::vector<int> const & rowind,
206  std::vector<int> const & colind,
207  std::vector<Treal> const & values,
208  std::vector<int> const & rowPermutation,
209  std::vector<int> const & colPermutation) {
210  std::vector<int> newRowind;
211  std::vector<int> newColind;
212  this->getPermutedAndSymmetrized(rowind, rowPermutation, newRowind,
213  colind, colPermutation, newColind);
214 
215  this->matrixPtr->syAssignFromSparse(newRowind, newColind, values);
216  }
222  inline void assign_from_sparse
223  (std::vector<int> const & rowind,
224  std::vector<int> const & colind,
225  std::vector<Treal> const & values,
226  SizesAndBlocks const & newRows,
227  SizesAndBlocks const & newCols,
228  std::vector<int> const & rowPermutation,
229  std::vector<int> const & colPermutation) {
230  this->resetSizesAndBlocks(newRows, newCols);
231  assign_from_sparse(rowind, colind, values,
232  rowPermutation, colPermutation);
233  }
241  inline void add_values
242  (std::vector<int> const & rowind,
243  std::vector<int> const & colind,
244  std::vector<Treal> const & values) {
245  this->matrixPtr->syAddValues(rowind, colind, values);
246  }
247 
251  inline void add_values
252  (std::vector<int> const & rowind,
253  std::vector<int> const & colind,
254  std::vector<Treal> const & values,
255  std::vector<int> const & rowPermutation,
256  std::vector<int> const & colPermutation) {
257  std::vector<int> newRowind;
258  std::vector<int> newColind;
259  this->getPermutedAndSymmetrized(rowind, rowPermutation, newRowind,
260  colind, colPermutation, newColind);
261  this->matrixPtr->syAddValues(newRowind, newColind, values);
262  }
263 
264 
265 
266  inline void get_values
267  (std::vector<int> const & rowind,
268  std::vector<int> const & colind,
269  std::vector<Treal> & values) const {
270  this->matrixPtr->syGetValues(rowind, colind, values);
271  }
279  inline void get_values
280  (std::vector<int> const & rowind,
281  std::vector<int> const & colind,
282  std::vector<Treal> & values,
283  std::vector<int> const & rowPermutation,
284  std::vector<int> const & colPermutation) const {
285  std::vector<int> newRowind;
286  std::vector<int> newColind;
287  this->getPermutedAndSymmetrized(rowind, rowPermutation, newRowind,
288  colind, colPermutation, newColind);
289  this->matrixPtr->syGetValues(newRowind, newColind, values);
290  }
296  inline void get_all_values
297  (std::vector<int> & rowind,
298  std::vector<int> & colind,
299  std::vector<Treal> & values) const {
300  rowind.resize(0);
301  colind.resize(0);
302  values.resize(0);
303  rowind.reserve(nnz());
304  colind.reserve(nnz());
305  values.reserve(nnz());
306  this->matrixPtr->syGetAllValues(rowind, colind, values);
307  }
313  inline void get_all_values
314  (std::vector<int> & rowind,
315  std::vector<int> & colind,
316  std::vector<Treal> & values,
317  std::vector<int> const & rowInversePermutation,
318  std::vector<int> const & colInversePermutation) const {
319  std::vector<int> tmpRowind;
320  std::vector<int> tmpColind;
321  tmpRowind.reserve(rowind.capacity());
322  tmpColind.reserve(colind.capacity());
323  values.resize(0);
324  this->matrixPtr->syGetAllValues(tmpRowind, tmpColind, values);
325  this->getPermutedAndSymmetrized(tmpRowind, rowInversePermutation, rowind,
326  tmpColind, colInversePermutation, colind);
327  }
341  return *this;
342  }
346  this->matrixPtr->nosymToSym();
347  return *this;
348  }
350  *this->matrixPtr = k;
351  return *this;
352  }
353 
354  inline Treal frob() const {
355  return this->matrixPtr->syFrob();
356  }
357  Treal mixed_norm(Treal const requestedAccuracy,
358  int maxIter = -1) const;
359  Treal eucl(Treal const requestedAccuracy,
360  int maxIter = -1) const;
361 
362  void quickEuclBounds(Treal & euclLowerBound,
363  Treal & euclUpperBound) const {
364  Treal frobTmp = frob();
365  euclLowerBound = frobTmp / template_blas_sqrt( (Treal)this->get_nrows() );
366  euclUpperBound = frobTmp;
367  }
368 
369 
375  static Interval<Treal>
378  normType const norm,
379  Treal const requestedAccuracy);
387  static Interval<Treal>
390  normType const norm,
391  Treal const requestedAccuracy,
392  Treal const maxAbsVal);
396  static inline Treal frob_diff
399  return Tmatrix::syFrobDiff(*A.matrixPtr, *B.matrixPtr);
400  }
401 
405  static Treal eucl_diff
408  Treal const requestedAccuracy,
409  int maxIter = -1);
410 
414  static Treal mixed_diff
417  Treal const requestedAccuracy);
418 
428  Treal const requestedAccuracy,
429  Treal const maxAbsVal,
430  VectorType * const eVecPtr = 0);
431 
432 
443  Treal thresh(Treal const threshold,
444  normType const norm);
445 
452  inline Treal frob_thresh(Treal const threshold) {
453  return 2.0 * this->matrixPtr->frob_thresh(threshold / 2);
454  }
455 
456  Treal eucl_thresh(Treal const threshold,
457  MatrixTriangular<Treal, Tmatrix> const * const Zptr = NULL);
458 
459  Treal eucl_element_level_thresh(Treal const threshold);
460 
462  ( SizesAndBlocks & rows_new, SizesAndBlocks & cols_new ) const;
463 
464  Treal mixed_norm_thresh(Treal const threshold);
465 
466  void simple_blockwise_frob_thresh(Treal const threshold) {
467  this->matrixPtr->frobThreshLowestLevel(threshold*threshold, 0);
468  }
469 
470  inline void gershgorin(Treal& lmin, Treal& lmax) const {
471  this->matrixPtr->sy_gershgorin(lmin, lmax);
472  }
473  static inline Treal trace_ab
476  return Tmatrix::sy_trace_ab(*A.matrixPtr, *B.matrixPtr);
477  }
478  inline size_t nnz() const { /* Note: size_t instead of int here to avoid integer overflow. */
479  return this->matrixPtr->sy_nnz();
480  }
481  inline size_t nvalues() const { /* Note: size_t instead of int here to avoid integer overflow. */
482  return this->matrixPtr->sy_nvalues();
483  }
484  inline void write_to_buffer(void* buffer, const int n_bytes) const {
485  this->write_to_buffer_base(buffer, n_bytes, matrix_symm);
486  }
487  inline void read_from_buffer(void* buffer, const int n_bytes) {
488  this->read_from_buffer_base(buffer, n_bytes, matrix_symm);
489  }
490 
491 
497  (const XYZpUV<Treal,
500  Treal,
504  (const XYZ<Treal,
509  (const XYZ<Treal,
514  (const XYZpUV<Treal,
517  Treal,
521  (const XYZ<Treal,
526  (const XYZ<Treal,
529 
537 
542  static void ssmmUpperTriangleOnly(const Treal alpha,
545  const Treal beta,
547 
548 
549  /* Addition */
553  MatrixSymmetric<Treal, Tmatrix> > const & mpm);
557  MatrixSymmetric<Treal, Tmatrix> > const & mm);
564 
568 
572 
573  template<typename Top>
574  Treal accumulateWith(Top & op) {
575  return this->matrixPtr->syAccumulateWith(op);
576  }
577 
578  void random() {
579  this->matrixPtr->syRandom();
580  }
581 
582  void randomZeroStructure(Treal probabilityBeingZero) {
583  this->matrixPtr->syRandomZeroStructure(probabilityBeingZero);
584  }
585 
590  template<typename TRule>
591  void setElementsByRule(TRule & rule) {
592  this->matrixPtr->sySetElementsByRule(rule);
593  return;
594  }
595 
599  ValidPtr<Tmatrix>::swap( this->matrixPtr, dest.matrixPtr );
600  // *this now contains previous content of dest
601  this->clear();
602  }
603 
604  template<typename Tvector>
605  void matVecProd(Tvector & y, Tvector const & x) const {
606  Treal const ONE = 1.0;
607  y = (ONE * (*this) * x);
608  }
609 
610 
611  std::string obj_type_id() const {return "MatrixSymmetric";}
612  protected:
613  inline void writeToFileProt(std::ofstream & file) const {
614  this->writeToFileBase(file, matrix_symm);
615  }
616  inline void readFromFileProt(std::ifstream & file) {
617  this->readFromFileBase(file, matrix_symm);
618  }
619 
625  static void getPermutedAndSymmetrized
626  (std::vector<int> const & rowind,
627  std::vector<int> const & rowPermutation,
628  std::vector<int> & newRowind,
629  std::vector<int> const & colind,
630  std::vector<int> const & colPermutation,
631  std::vector<int> & newColind) {
633  getPermutedIndexes(rowind, rowPermutation, newRowind);
635  getPermutedIndexes(colind, colPermutation, newColind);
636  int tmp;
637  for (unsigned int i = 0; i < newRowind.size(); ++i) {
638  if (newRowind[i] > newColind[i]) {
639  tmp = newRowind[i];
640  newRowind[i] = newColind[i];
641  newColind[i] = tmp;
642  }
643  }
644  }
645  private:
646  };
647 
648  template<typename Treal, typename Tmatrix>
650  mixed_norm(Treal const requestedAccuracy,
651  int maxIter) const {
652  // Construct SizesAndBlocks for frobNormMat
653  SizesAndBlocks rows_new;
654  SizesAndBlocks cols_new;
655  this->getSizesAndBlocksForFrobNormMat( rows_new, cols_new );
656  // Now we can construct an empty matrix where the Frobenius norms
657  // of lowest level nonzero submatrices will be stored
659  frobNormMat.resetSizesAndBlocks(rows_new, cols_new);
660  frobNormMat.getMatrix().syAssignFrobNormsLowestLevel(this->getMatrix());
661  return frobNormMat.eucl(requestedAccuracy, maxIter);
662  }
663 
664 
665  template<typename Treal, typename Tmatrix>
667  eucl(Treal const requestedAccuracy,
668  int maxIter) const {
669  assert(requestedAccuracy >= 0);
670  /* Check if norm is really small, in that case quick return */
671  Treal frobTmp = this->frob();
672  if (frobTmp < requestedAccuracy)
673  return (Treal)0.0;
674  if (maxIter < 0)
675  maxIter = this->get_nrows() * 100;
676  VectorType guess;
678  this->getCols(cols);
679  guess.resetSizesAndBlocks(cols);
680  guess.rand();
681  // Elias note 2010-03-26: changed this back from "new code" to "old code" to reduce memory usage.
682 #if 0 // "new code"
684  Copy.frob_thresh(requestedAccuracy / 2.0);
687  lan(Copy, guess, maxIter);
688  lan.setAbsTol( requestedAccuracy / 2.0 );
689 #else // "old code"
692  lan(*this, guess, maxIter);
693  lan.setAbsTol( requestedAccuracy );
694 #endif
695  lan.run();
696  Treal eVal = 0;
697  Treal acc = 0;
698  lan.getLargestMagnitudeEig(eVal, acc);
699  return template_blas_fabs(eVal);
700  }
701 
702  template<typename Treal, typename Tmatrix>
706  normType const norm, Treal const requestedAccuracy) {
707  Treal diff;
708  Treal eNMin;
709  switch (norm) {
710  case frobNorm:
711  diff = frob_diff(A, B);
712  return Interval<Treal>(diff / template_blas_sqrt((Treal)A.get_nrows()), diff);
713  break;
714  case euclNorm:
715  diff = eucl_diff(A, B, requestedAccuracy);
716  eNMin = diff - requestedAccuracy;
717  eNMin = eNMin >= 0 ? eNMin : 0;
718  return Interval<Treal>(eNMin, diff + requestedAccuracy);
719  break;
720  default:
721  throw Failure("MatrixSymmetric<Treal, Tmatrix>::"
722  "diff(const MatrixSymmetric<Treal, Tmatrix>&, "
723  "const MatrixSymmetric<Treal, Tmatrix>&, "
724  "normType const, Treal): "
725  "Diff not implemented for selected norm");
726  }
727  }
728 
729 #if 1
730  template<typename Treal, typename Tmatrix>
734  normType const norm,
735  Treal const requestedAccuracy,
736  Treal const maxAbsVal) {
737  Treal diff;
738  switch (norm) {
739  case frobNorm:
740  {
741  diff = frob_diff(A, B);
742  return Interval<Treal>(diff / template_blas_sqrt((Treal)A.get_nrows()), diff);
743  }
744  break;
745  case euclNorm:
746  return euclDiffIfSmall(A, B, requestedAccuracy, maxAbsVal);
747  break;
748  default:
749  throw Failure("MatrixSymmetric<Treal, Tmatrix>::"
750  "diffIfSmall"
751  "(const MatrixSymmetric<Treal, Tmatrix>&, "
752  "const MatrixSymmetric<Treal, Tmatrix>&, "
753  "normType const, Treal const, Treal const): "
754  "Diff not implemented for selected norm");
755  }
756  }
757 
758 #endif
759 
760 
761  template<typename Treal, typename Tmatrix>
765  Treal const requestedAccuracy,
766  int maxIter) {
767  // DiffMatrix is a lightweight proxy object:
769  Treal maxAbsVal = 2 * frob_diff(A,B);
770  // Now, maxAbsVal should be larger than the Eucl norm
771  // Note that mat::euclIfSmall lies outside this class
772  Treal relTol = template_blas_sqrt(template_blas_sqrt(mat::getMachineEpsilon<Treal>()));
773  VectorType * const eVecPtrNotUsed = 0;
774  Interval<Treal> euclInt =
775  mat::euclIfSmall(Diff, requestedAccuracy, relTol, maxAbsVal, eVecPtrNotUsed, maxIter);
776  return euclInt.midPoint();
777  }
778 
779  template<typename Treal, typename Tmatrix>
783  Treal const requestedAccuracy) {
785  {
786  SizesAndBlocks rows_new;
787  SizesAndBlocks cols_new;
788  A.getSizesAndBlocksForFrobNormMat( rows_new, cols_new );
789  frobNormMat.resetSizesAndBlocks(rows_new, cols_new);
790  frobNormMat.getMatrix().syAssignDiffFrobNormsLowestLevel(A.getMatrix(),B.getMatrix());
791  }
792  return frobNormMat.eucl(requestedAccuracy);
793  }
794 
795 #if 1
796  template<typename Treal, typename Tmatrix>
800  Treal const requestedAccuracy,
801  Treal const maxAbsVal,
802  VectorType * const eVecPtr) {
803  // DiffMatrix is a lightweight proxy object:
805  // Note that this function lies outside this class
806  Treal relTol = template_blas_sqrt(template_blas_sqrt(mat::getMachineEpsilon<Treal>()));
807  Interval<Treal> tmpInterval = mat::euclIfSmall(Diff, requestedAccuracy, relTol, maxAbsVal, eVecPtr);
808  // Emanuel note: Ugly fix to make certain tests pass, we expand
809  // the interval up to the requested accuracy. Note that larger
810  // intervals may occur if the norm is not 'small'. It happens that
811  // Lanczos misconverges to for example the second largest
812  // eigenvalue. This happens in particular when the first and second
813  // eigenvalues are very close (of the order of the requested
814  // accuracy). Expanding the interval makes the largest eigenvalue
815  // (at least for certain cases) end up inside the interval even
816  // though Lanczos has misconverged.
817  if ( tmpInterval.length() < 2*requestedAccuracy )
818  return Interval<Treal>( tmpInterval.midPoint()-requestedAccuracy, tmpInterval.midPoint()+requestedAccuracy );
819  return tmpInterval;
820  }
821 
822 #endif
823 
824 
825  template<typename Treal, typename Tmatrix>
827  thresh(Treal const threshold,
828  normType const norm) {
829  switch (norm) {
830  case frobNorm:
831  return this->frob_thresh(threshold);
832  break;
833  case euclNorm:
834  return this->eucl_thresh(threshold);
835  break;
836  case mixedNorm:
837  return this->mixed_norm_thresh(threshold);
838  break;
839  default:
840  throw Failure("MatrixSymmetric<Treal, Tmatrix>::"
841  "thresh(Treal const, "
842  "normType const): "
843  "Thresholding not implemented for selected norm");
844  }
845  }
846 
847 #if 1
848 
849  template<typename Treal, typename Tmatrix>
851  eucl_thresh(Treal const threshold,
852  MatrixTriangular<Treal, Tmatrix> const * const Zptr) {
853  if ( Zptr == NULL ) {
855  return TruncObj.run( threshold );
856  }
858  return TruncObj.run( threshold );
859  }
860 
861 #endif
862 
863 
864  template<typename Treal, typename Tmatrix>
866  ( SizesAndBlocks & rows_new, SizesAndBlocks & cols_new ) const {
867  std::vector<int> rows_block_sizes;
868  std::vector<int> cols_block_sizes;
869  int n_rows;
870  int n_cols;
871  {
874  this->getRows(rows);
875  this->getCols(cols);
876  rows.getBlockSizeVector( rows_block_sizes );
877  cols.getBlockSizeVector( cols_block_sizes );
878  rows_block_sizes.pop_back(); // Remove the '1' at the end
879  cols_block_sizes.pop_back(); // Remove the '1' at the end
880  n_rows = rows.getNTotalScalars();
881  n_cols = cols.getNTotalScalars();
882  int factor_rows = rows_block_sizes[rows_block_sizes.size()-1];
883  int factor_cols = cols_block_sizes[cols_block_sizes.size()-1];
884  for (unsigned int ind = 0; ind < rows_block_sizes.size(); ++ind)
885  rows_block_sizes[ind] = rows_block_sizes[ind] / factor_rows;
886  for (unsigned int ind = 0; ind < cols_block_sizes.size(); ++ind)
887  cols_block_sizes[ind] = cols_block_sizes[ind] / factor_cols;
888  // Now set the number of (scalar) rows and cols, should be equal
889  // to the number of blocks at the lowest level of the original
890  // matrix
891  if (n_rows % factor_rows)
892  n_rows = n_rows / factor_rows + 1;
893  else
894  n_rows = n_rows / factor_rows;
895  if (n_cols % factor_cols)
896  n_cols = n_cols / factor_cols + 1;
897  else
898  n_cols = n_cols / factor_cols;
899  }
900  rows_new = SizesAndBlocks( rows_block_sizes, n_rows );
901  cols_new = SizesAndBlocks( cols_block_sizes, n_cols );
902  }
903 
904 
905  template<typename Treal, typename Tmatrix>
907  mixed_norm_thresh(Treal const threshold) {
908  assert(threshold >= (Treal)0.0);
909  if (threshold == (Treal)0.0)
910  return (Treal)0;
911  // Construct SizesAndBlocks for frobNormMat
912  SizesAndBlocks rows_new;
913  SizesAndBlocks cols_new;
914  this->getSizesAndBlocksForFrobNormMat( rows_new, cols_new );
915  // Now we can construct an empty matrix where the Frobenius norms
916  // of lowest level nonzero submatrices will be stored
918  frobNormMat.resetSizesAndBlocks(rows_new, cols_new);
919  // We want the following step to dominate the mixed_norm truncation (this
920  // is where Frobenius norms of submatrices are computed, i.e. it
921  // is here we loop over all matrix elements.)
922  frobNormMat.getMatrix().syAssignFrobNormsLowestLevel(this->getMatrix());
924  Treal mixed_norm_result = TruncObj.run( threshold );
925  frobNormMat.getMatrix().truncateAccordingToSparsityPattern(this->getMatrix());
926  return mixed_norm_result;
927  }
928 
929 
930  /* B = alpha * A */
931  template<typename Treal, typename Tmatrix>
935 
936  if(this == &sm.B) // A = B
937  {
938  *this *= sm.A; // B *= alpha
939  return *this;
940  }
941  assert(!sm.tB);
942  /* Data structure set by assign - therefore set haveDataStructure to true */
943  this->matrixPtr.haveDataStructureSet(true);
944  this->matrixPtr->assign(sm.A, *sm.B.matrixPtr);
945  return *this;
946  }
947  /* C = alpha * A * A + beta * C : A and C are symmetric */
948  template<typename Treal, typename Tmatrix>
951  (const XYZpUV<Treal,
954  Treal,
956  assert(this != &sm2psm.B);
957  if (this == &sm2psm.E && &sm2psm.B == &sm2psm.C) {
958  /* Operation is C = alpha * A * A + beta * C */
959  Tmatrix::sysq('U',
960  sm2psm.A, *sm2psm.B.matrixPtr,
961  sm2psm.D, *this->matrixPtr);
962  return *this;
963  }
964  else /* this != &sm2psm.C || &sm2psm.B != &sm2psm.C */
965  throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator="
966  "(const XYZpUV<Treal, MatrixSymmetric"
967  "<Treal, Tmatrix> >& sm2psm) : "
968  "D = alpha * A * B + beta * C not supported for C != D"
969  " and for symmetric matrices not for A != B since this "
970  "generally will result in a nonsymmetric matrix");
971  }
972 
973  /* C = alpha * A * A : A and C are symmetric */
974  template<typename Treal, typename Tmatrix>
977  (const XYZ<Treal,
980  assert(this != &sm2.B);
981  if (&sm2.B == &sm2.C) {
982  this->matrixPtr.haveDataStructureSet(true);
983  Tmatrix::sysq('U', sm2.A, *sm2.B.matrixPtr, 0, *this->matrixPtr);
984  return *this;
985  }
986  else
987  throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator="
988  "(const XYZ<Treal, MatrixSymmetric<Treal, Tmatrix>,"
989  " MatrixSymmetric<Treal, Tmatrix> >& sm2) : "
990  "Operation C = alpha * A * B with only symmetric "
991  "matrices not supported for A != B");
992  }
993 
994  /* C += alpha * A * A : A and C are symmetric */
995  template<typename Treal, typename Tmatrix>
998  (const XYZ<Treal,
1001  assert(this != &sm2.B);
1002  if (&sm2.B == &sm2.C) {
1003  Tmatrix::sysq('U', sm2.A, *sm2.B.matrixPtr, 1, *this->matrixPtr);
1004  return *this;
1005  }
1006  else
1007  throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator+="
1008  "(const XYZ<Treal, MatrixSymmetric<Treal, Tmatrix>,"
1009  " MatrixSymmetric<Treal, Tmatrix> >& sm2) : "
1010  "Operation C += alpha * A * B with only symmetric "
1011  "matrices not supported for A != B");
1012  }
1013 
1014  /* C = alpha * A * transpose(A) + beta * C : C is symmetric */
1015  template<typename Treal, typename Tmatrix>
1018  (const XYZpUV<Treal,
1021  Treal,
1022  MatrixSymmetric<Treal, Tmatrix> >& smmpsm) {
1023  if (this == &smmpsm.E)
1024  if (&smmpsm.B == &smmpsm.C)
1025  if (!smmpsm.tB && smmpsm.tC) {
1026  Tmatrix::syrk('U', false,
1027  smmpsm.A, *smmpsm.B.matrixPtr,
1028  smmpsm.D, *this->matrixPtr);
1029  }
1030  else
1031  throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator="
1032  "(const XYZpUV<Treal, MatrixGeneral"
1033  "<Treal, Tmatrix>, "
1034  "MatrixGeneral<Treal, Tmatrix>, Treal, "
1035  "MatrixSymmetric<Treal, Tmatrix> >&) : "
1036  "C = alpha * A' * A + beta * C, not implemented"
1037  " only C = alpha * A * A' + beta * C");
1038  else
1039  throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator="
1040  "(const XYZpUV<"
1041  "Treal, MatrixGeneral<Treal, Tmatrix>, "
1042  "MatrixGeneral<Treal, Tmatrix>, Treal, "
1043  "MatrixSymmetric<Treal, Tmatrix> >&) : "
1044  "You are trying to call C = alpha * A * A' + beta * C"
1045  " with A and A' being different objects");
1046  else
1047  throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator="
1048  "(const XYZpUV<"
1049  "Treal, MatrixGeneral<Treal, Tmatrix>, "
1050  "MatrixGeneral<Treal, Tmatrix>, Treal, "
1051  "MatrixSymmetric<Treal, Tmatrix> >&) : "
1052  "D = alpha * A * A' + beta * C not supported for C != D");
1053  return *this;
1054  }
1055 
1056  /* C = alpha * A * transpose(A) : C is symmetric */
1057  template<typename Treal, typename Tmatrix>
1060  (const XYZ<Treal,
1063  if (&smm.B == &smm.C)
1064  if (!smm.tB && smm.tC) {
1065  Tmatrix::syrk('U', false,
1066  smm.A, *smm.B.matrixPtr,
1067  0, *this->matrixPtr);
1068  }
1069  else
1070  throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator="
1071  "(const XYZ<"
1072  "Treal, MatrixGeneral<Treal, Tmatrix>, "
1073  "MatrixGeneral<Treal, Tmatrix> >&) : "
1074  "C = alpha * A' * A, not implemented "
1075  "only C = alpha * A * A'");
1076  else
1077  throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator="
1078  "(const XYZ<"
1079  "Treal, MatrixGeneral<Treal, Tmatrix>, "
1080  "MatrixGeneral<Treal, Tmatrix> >&) : "
1081  "You are trying to call C = alpha * A * A' "
1082  "with A and A' being different objects");
1083  return *this;
1084  }
1085 
1086  /* C += alpha * A * transpose(A) : C is symmetric */
1087  template<typename Treal, typename Tmatrix>
1090  (const XYZ<Treal,
1093  if (&smm.B == &smm.C)
1094  if (!smm.tB && smm.tC) {
1095  Tmatrix::syrk('U', false,
1096  smm.A, *smm.B.matrixPtr,
1097  1, *this->matrixPtr);
1098  }
1099  else
1100  throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator+="
1101  "(const XYZ<"
1102  "Treal, MatrixGeneral<Treal, Tmatrix>, "
1103  "MatrixGeneral<Treal, Tmatrix> >&) : "
1104  "C += alpha * A' * A, not implemented "
1105  "only C += alpha * A * A'");
1106  else
1107  throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator+="
1108  "(const XYZ<"
1109  "Treal, MatrixGeneral<Treal, Tmatrix>, "
1110  "MatrixGeneral<Treal, Tmatrix> >&) : "
1111  "You are trying to call C += alpha * A * A' "
1112  "with A and A' being different objects");
1113  return *this;
1114  }
1115 
1116 #if 1
1117  /* A = op1(Z) * A * op2(Z) : Z is upper triangular and A is symmetric */
1118  /* Either op1() or op2() is the transpose operation. */
1119  template<typename Treal, typename Tmatrix>
1125  if (this == &zaz.B) {
1126  if (&zaz.A == &zaz.C) {
1127  if (zaz.tA && !zaz.tC) {
1128  Tmatrix::trsytriplemm('R', *zaz.A.matrixPtr, *this->matrixPtr);
1129  }
1130  else if (!zaz.tA && zaz.tC) {
1131  Tmatrix::trsytriplemm('L', *zaz.A.matrixPtr, *this->matrixPtr);
1132  }
1133  else
1134  throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator="
1135  "(const XYZ<MatrixTriangular<Treal, Tmatrix>,"
1136  "MatrixSymmetric<Treal, Tmatrix>,"
1137  "MatrixTriangular<Treal, Tmatrix> >&) : "
1138  "A = op1(Z) * A * op2(Z) : Either op1 xor op2 must be "
1139  "the transpose operation.");
1140  }
1141  else
1142  throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator="
1143  "(const XYZ<MatrixTriangular<Treal, Tmatrix>,"
1144  "MatrixSymmetric<Treal, Tmatrix>,"
1145  "MatrixTriangular<Treal, Tmatrix> >&) : "
1146  "A = op1(Z1) * A * op2(Z2) : Z1 and Z2 must be the same "
1147  "object");
1148  }
1149  else
1150  throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator="
1151  "(const XYZ<MatrixTriangular<Treal, Tmatrix>,"
1152  "MatrixSymmetric<Treal, Tmatrix>,"
1153  "MatrixTriangular<Treal, Tmatrix> >&) : "
1154  "C = op1(Z) * A * op2(Z) : A and C must be the same "
1155  "object");
1156  return *this;
1157  }
1158 
1159 #endif
1160 
1161 
1166  template<typename Treal, typename Tmatrix>
1168  ssmmUpperTriangleOnly(const Treal alpha,
1171  const Treal beta,
1173  Tmatrix::ssmm_upper_tr_only(alpha, *A.matrixPtr, *B.matrixPtr,
1174  beta, *C.matrixPtr);
1175  }
1176 
1177 
1178 
1179 
1180 
1181  /* Addition */
1182  /* C = A + B */
1183  template<typename Treal, typename Tmatrix>
1188  assert(this != &mpm.A);
1189  (*this) = mpm.B;
1190  Tmatrix::add(1.0, *mpm.A.matrixPtr, *this->matrixPtr);
1191  return *this;
1192  }
1193  /* C = A - B */
1194  template<typename Treal, typename Tmatrix>
1199  assert(this != &mmm.B);
1200  (*this) = mmm.A;
1201  Tmatrix::add(-1.0, *mmm.B.matrixPtr, *this->matrixPtr);
1202  return *this;
1203  }
1204  /* B += A */
1205  template<typename Treal, typename Tmatrix>
1209  Tmatrix::add(1.0, *A.matrixPtr, *this->matrixPtr);
1210  return *this;
1211  }
1212  /* B -= A */
1213  template<typename Treal, typename Tmatrix>
1217  Tmatrix::add(-1.0, *A.matrixPtr, *this->matrixPtr);
1218  return *this;
1219  }
1220  /* B += alpha * A */
1221  template<typename Treal, typename Tmatrix>
1225  assert(!sm.tB);
1226  Tmatrix::add(sm.A, *sm.B.matrixPtr, *this->matrixPtr);
1227  return *this;
1228  }
1229 
1230  /* B -= alpha * A */
1231  template<typename Treal, typename Tmatrix>
1235  assert(!sm.tB);
1236  Tmatrix::add(-sm.A, *sm.B.matrixPtr, *this->matrixPtr);
1237  return *this;
1238  }
1239 
1244  template<typename Treal, typename Tmatrix, typename Top>
1246  Top & op) {
1247  return A.accumulateWith(op);
1248  }
1249 
1250 
1251 
1252 } /* end namespace mat */
1253 #endif
1254 
void random()
Definition: MatrixSymmetric.h:578
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:242
#define A
Treal mixed_norm_thresh(Treal const threshold)
Definition: MatrixSymmetric.h:907
Treal frob_thresh(Treal const threshold)
Does thresholding so that the error in the Frobenius norm is below the given threshold.
Definition: MatrixSymmetric.h:452
mat::SizesAndBlocks rows
Definition: test.cc:51
Normal matrix.
Definition: MatrixBase.h:49
mat::SizesAndBlocks cols
Definition: test.cc:52
MatrixSymmetric< Treal, Tmatrix > & operator=(const MatrixSymmetric< Treal, Tmatrix > &symm)
Definition: MatrixSymmetric.h:339
MatrixBase< Treal, Tmatrix > & operator=(const MatrixBase< Treal, Tmatrix > &other)
Definition: MatrixBase.h:166
ValidPtr< Tmatrix > matrixPtr
Definition: MatrixBase.h:153
MatrixSymmetric< Treal, Tmatrix > & operator=(const MatrixGeneral< Treal, Tmatrix > &matr)
Definition: MatrixSymmetric.h:344
Definition: LanczosLargestMagnitudeEig.h:46
void read_from_buffer_base(void *buffer, const int n_bytes, const matrix_type mattype)
Definition: MatrixBase.h:281
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
Truncation of symmetric matrices at the element level (used for mixed norm truncation) ...
Definition: truncation.h:246
MatrixSymmetric< Treal, Tmatrix > & operator=(int const k)
Definition: MatrixSymmetric.h:349
Treal accumulateWith(Top &op)
Definition: MatrixSymmetric.h:574
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
MatrixSymmetric(const MatrixSymmetric< Treal, Tmatrix > &symm)
Copy constructor.
Definition: MatrixSymmetric.h:75
This proxy expresses the result of addition of two objects, of possibly different types...
Definition: matrix_proxy.h:246
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:798
MatrixSymmetric(const MatrixGeneral< Treal, Tmatrix > &matr)
&#39;Copy from normal matrix&#39; - constructor
Definition: MatrixSymmetric.h:79
void quickEuclBounds(Treal &euclLowerBound, Treal &euclUpperBound) const
Definition: MatrixSymmetric.h:362
Definition: MatrixBase.h:55
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:267
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
void clear()
Release memory for the information written to file.
Definition: MatrixBase.h:118
Definition: allocate.cc:39
Describes dimensions of matrix and its blocks on all levels.
Definition: SizesAndBlocks.h:45
int getNTotalScalars() const
Definition: SizesAndBlocks.h:84
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:732
Treal template_blas_fabs(Treal x)
Definition: Interval.h:46
static Treal mixed_diff(const MatrixSymmetric< Treal, Tmatrix > &A, const MatrixSymmetric< Treal, Tmatrix > &B, Treal const requestedAccuracy)
Returns the &#39;mixed&#39; norm of A - B ( || A - B ||_mixed )
Definition: MatrixSymmetric.h:781
void assignFromFull(std::vector< Treal > const &fullMat)
Definition: MatrixSymmetric.h:104
Utilities used by other parts of the hierarchical matrix library.
void readFromFileBase(std::ifstream &file, matrix_type const mattype)
Definition: MatrixBase.h:236
void read_from_buffer(void *buffer, const int n_bytes)
Definition: MatrixSymmetric.h:487
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:1168
static void swap(ValidPtr< Tobj > &ptrA, ValidPtr< Tobj > &ptrB)
Definition: ValidPtr.h:106
MatrixSymmetric(const XY< Treal, MatrixSymmetric< Treal, Tmatrix > > &sm)
Definition: MatrixSymmetric.h:77
Definition: matInclude.h:139
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:763
int get_ncols() const
Definition: MatrixBase.h:141
void fullMatrix(std::vector< Treal > &fullMat) const
Save matrix as full matrix.
Definition: MatrixSymmetric.h:132
This proxy expresses the result of multiplication of three objects, of possibly different types...
Definition: matrix_proxy.h:67
void write_to_buffer_base(void *buffer, const int n_bytes, const matrix_type mattype) const
Definition: MatrixBase.h:254
Upper non-unit triangular matrix.
Definition: MatrixBase.h:53
MatrixSymmetric()
Default constructor.
Definition: MatrixSymmetric.h:73
void randomZeroStructure(Treal probabilityBeingZero)
Definition: MatrixSymmetric.h:582
Treal run(Treal const threshold)
Definition: truncation.h:80
void writeToFileBase(std::ofstream &file, matrix_type const mattype) const
Definition: MatrixBase.h:221
void gershgorin(Treal &lmin, Treal &lmax) const
Definition: MatrixSymmetric.h:470
Treal midPoint() const
Definition: Interval.h:115
Definition: matInclude.h:139
void rand()
Definition: VectorGeneral.h:102
void resetSizesAndBlocks(SizesAndBlocks const &newRows, SizesAndBlocks const &newCols)
Definition: MatrixBase.h:76
void writeToFileProt(std::ofstream &file) const
Write object to file.
Definition: MatrixSymmetric.h:613
void getBlockSizeVector(std::vector< int > &blockSizesCopy) const
Definition: SizesAndBlocks.cc:87
void resetSizesAndBlocks(SizesAndBlocks const &newRows)
Definition: VectorGeneral.h:51
size_t nvalues() const
Definition: MatrixSymmetric.h:481
size_t nnz() const
Definition: MatrixSymmetric.h:478
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:170
Base class for matrix API.
Definition: MatrixBase.h:69
void readFromFileProt(std::ifstream &file)
Read object from file.
Definition: MatrixSymmetric.h:616
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:397
Definition: mat_utils.h:43
void transfer(MatrixSymmetric< Treal, Tmatrix > &dest)
Transfer this matrix to dest, clearing previous content of dest if any.
Definition: MatrixSymmetric.h:598
void getLargestMagnitudeEig(Treal &ev, Treal &accuracy)
Definition: LanczosLargestMagnitudeEig.h:58
This proxy expresses the result of substraction of two objects, of possibly different types...
Definition: matrix_proxy.h:266
std::string obj_type_id() const
Definition: MatrixSymmetric.h:611
static Treal trace_ab(const MatrixSymmetric< Treal, Tmatrix > &A, const MatrixSymmetric< Treal, Tmatrix > &B)
Definition: MatrixSymmetric.h:474
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:297
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:704
Base class for matrix API.
void matVecProd(Tvector &y, Tvector const &x) const
Definition: MatrixSymmetric.h:605
Treal length() const
Returns the length of the interval.
Definition: Interval.h:109
Definition: matInclude.h:139
Truncation of symmetric matrices.
Definition: truncation.h:154
void getSizesAndBlocksForFrobNormMat(SizesAndBlocks &rows_new, SizesAndBlocks &cols_new) const
Definition: MatrixSymmetric.h:866
Definition: Failure.h:57
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:626
Truncation of symmetric matrices with Z.
Definition: truncation.h:203
This proxy expresses the result of multiplication of two objects, of possibly different types...
Definition: matrix_proxy.h:51
void write_to_buffer(void *buffer, const int n_bytes) const
Definition: MatrixSymmetric.h:484
void simple_blockwise_frob_thresh(Treal const threshold)
Definition: MatrixSymmetric.h:466
Treal eucl(Treal const requestedAccuracy, int maxIter=-1) const
Definition: MatrixSymmetric.h:667
Treal frob() const
Definition: MatrixSymmetric.h:354
Treal eucl_element_level_thresh(Treal const threshold)
Treal accumulate(MatrixSymmetric< Treal, Tmatrix > &A, Top &op)
Performs operation specified in &#39;op&#39; on all nonzero matrix elements and sums up the result and return...
Definition: MatrixSymmetric.h:1245
#define B
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:591
Interval class.
VectorGeneral< Treal, typename Tmatrix::VectorType > VectorType
Definition: MatrixSymmetric.h:70
int get_nrows() const
Definition: MatrixBase.h:138
Treal eucl_thresh(Treal const threshold, MatrixTriangular< Treal, Tmatrix > const *const Zptr=NULL)
Definition: MatrixSymmetric.h:851
Symmetric matrix.
Definition: MatrixBase.h:51
Class for computing the largest magnitude eigenvalue of a symmetric matrix with the Lanczos method...
This proxy expresses the result of multiplication of three objects added to two other multiplied obje...
Definition: matrix_proxy.h:88
Definition: MatrixBase.h:56
Treal template_blas_sqrt(Treal x)
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:827
normType
Definition: matInclude.h:139
Classes for truncation of small matrix elements.
Treal real
Definition: MatrixSymmetric.h:71
Treal mixed_norm(Treal const requestedAccuracy, int maxIter=-1) const
Definition: MatrixSymmetric.h:650
static void getPermutedIndexes(std::vector< int > const &index, std::vector< int > const &permutation, std::vector< int > &newIndex)
Definition: MatrixBase.h:205