ergo
MatrixSymmetric.h
Go to the documentation of this file.
1 /* Ergo, version 3.2, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2012 Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek.
4  *
5  * This program is free software: you can redistribute it and/or modify
6  * it under the terms of the GNU General Public License as published by
7  * the Free Software Foundation, either version 3 of the License, or
8  * (at your option) any later version.
9  *
10  * This program is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  * GNU General Public License for more details.
14  *
15  * You should have received a copy of the GNU General Public License
16  * along with this program. If not, see <http://www.gnu.org/licenses/>.
17  *
18  * Primary academic reference:
19  * Kohn−Sham Density Functional Theory Electronic Structure Calculations
20  * with Linearly Scaling Computational Time and Memory Usage,
21  * Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek,
22  * J. Chem. Theory Comput. 7, 340 (2011),
23  * <http://dx.doi.org/10.1021/ct100611z>
24  *
25  * For further information about Ergo, see <http://www.ergoscf.org>.
26  */
27 
36 #ifndef MAT_MatrixSymmetric
37 #define MAT_MatrixSymmetric
38 
39 #include <algorithm>
40 
41 #include "MatrixBase.h"
42 #include "Interval.h"
44 #include "mat_utils.h"
45 #include "truncation.h"
46 
47 namespace mat {
65  template<typename Treal, typename Tmatrix>
66  class MatrixSymmetric : public MatrixBase<Treal, Tmatrix> {
67  public:
69  typedef Treal real;
70 
72  :MatrixBase<Treal, Tmatrix>() {}
74  :MatrixBase<Treal, Tmatrix>(symm) {}
75  explicit MatrixSymmetric(const XY<Treal, MatrixSymmetric<Treal, Tmatrix> >& sm)
76  :MatrixBase<Treal, Tmatrix>() { *this = sm.A * sm.B; }
78  :MatrixBase<Treal, Tmatrix>(matr) {
79  this->matrixPtr->nosymToSym();
80  }
83 #if 0
84  template<typename Tfull>
85  inline void assign_from_full
86  (Tfull const* const fullmatrix,
87  int const totnrows, int const totncols) {
88  assert(totnrows == totncols);
89  this->matrixPtr->assign_from_full(fullmatrix, totnrows, totncols);
90  this->matrixPtr->nosym_to_sym();
91  }
92  inline void assign_from_full
93  (Treal const* const fullmatrix,
94  int const totnrows, int const totncols) {
95  assert(totnrows == totncols);
96  this->matrixPtr->assign_from_full(fullmatrix, totnrows, totncols);
97  this->matrixPtr->nosym_to_sym();
98  }
99 #endif
100 
101  inline void assignFromFull
102  (std::vector<Treal> const & fullMat) {
103  assert((int)fullMat.size() == this->get_nrows() * this->get_ncols());
104  this->matrixPtr->assignFromFull(fullMat);
105  this->matrixPtr->nosymToSym();
106  }
107 
108  inline void assignFromFull
109  (std::vector<Treal> const & fullMat,
110  std::vector<int> const & rowPermutation,
111  std::vector<int> const & colPermutation) {
112  assert((int)fullMat.size() == this->get_nrows() * this->get_ncols());
113  std::vector<int> rowind(fullMat.size());
114  std::vector<int> colind(fullMat.size());
115  int ind = 0;
116  for (int col = 0; col < this->get_ncols(); ++col)
117  for (int row = 0; row < this->get_nrows(); ++row) {
118  rowind[ind] = row;
119  colind[ind] = col;
120  ++ind;
121  }
122  this->assign_from_sparse(rowind,
123  colind,
124  fullMat,
125  rowPermutation,
126  colPermutation);
127  this->matrixPtr->nosymToSym();
128  }
129 
130  inline void fullMatrix(std::vector<Treal> & fullMat) const {
131  this->matrixPtr->syFullMatrix(fullMat);
132  }
139  inline void fullMatrix
140  (std::vector<Treal> & fullMat,
141  std::vector<int> const & rowInversePermutation,
142  std::vector<int> const & colInversePermutation) const {
143  std::vector<int> rowind;
144  std::vector<int> colind;
145  std::vector<Treal> values;
146  get_all_values(rowind, colind, values,
147  rowInversePermutation,
148  colInversePermutation);
149  fullMat.assign(this->get_nrows()*this->get_ncols(),0);
150  assert(rowind.size() == values.size());
151  assert(rowind.size() == colind.size());
152  for (unsigned int ind = 0; ind < values.size(); ++ind) {
153  assert(rowind[ind] + this->get_nrows() * colind[ind] <
154  this->get_nrows()*this->get_ncols());
155  fullMat[rowind[ind] + this->get_nrows() * colind[ind]] =
156  values[ind];
157  fullMat[colind[ind] + this->get_nrows() * rowind[ind]] =
158  values[ind];
159  }
160  }
167  inline void assign_from_sparse
168  (std::vector<int> const & rowind,
169  std::vector<int> const & colind,
170  std::vector<Treal> const & values) {
171  this->matrixPtr->syAssignFromSparse(rowind, colind, values);
172  }
184  inline void assign_from_sparse
185  (std::vector<int> const & rowind,
186  std::vector<int> const & colind,
187  std::vector<Treal> const & values,
188  SizesAndBlocks const & newRows,
189  SizesAndBlocks const & newCols) {
190  this->resetSizesAndBlocks(newRows, newCols);
191  this->matrixPtr->syAssignFromSparse(rowind, colind, values);
192  }
202  inline void assign_from_sparse
203  (std::vector<int> const & rowind,
204  std::vector<int> const & colind,
205  std::vector<Treal> const & values,
206  std::vector<int> const & rowPermutation,
207  std::vector<int> const & colPermutation) {
208  std::vector<int> newRowind;
209  std::vector<int> newColind;
210  this->getPermutedAndSymmetrized(rowind, rowPermutation, newRowind,
211  colind, colPermutation, newColind);
212 
213  this->matrixPtr->syAssignFromSparse(newRowind, newColind, values);
214  }
220  inline void assign_from_sparse
221  (std::vector<int> const & rowind,
222  std::vector<int> const & colind,
223  std::vector<Treal> const & values,
224  SizesAndBlocks const & newRows,
225  SizesAndBlocks const & newCols,
226  std::vector<int> const & rowPermutation,
227  std::vector<int> const & colPermutation) {
228  this->resetSizesAndBlocks(newRows, newCols);
229  assign_from_sparse(rowind, colind, values,
230  rowPermutation, colPermutation);
231  }
239  inline void add_values
240  (std::vector<int> const & rowind,
241  std::vector<int> const & colind,
242  std::vector<Treal> const & values) {
243  this->matrixPtr->syAddValues(rowind, colind, values);
244  }
245 
249  inline void add_values
250  (std::vector<int> const & rowind,
251  std::vector<int> const & colind,
252  std::vector<Treal> const & values,
253  std::vector<int> const & rowPermutation,
254  std::vector<int> const & colPermutation) {
255  std::vector<int> newRowind;
256  std::vector<int> newColind;
257  this->getPermutedAndSymmetrized(rowind, rowPermutation, newRowind,
258  colind, colPermutation, newColind);
259  this->matrixPtr->syAddValues(newRowind, newColind, values);
260  }
261 
262 
263 
264  inline void get_values
265  (std::vector<int> const & rowind,
266  std::vector<int> const & colind,
267  std::vector<Treal> & values) const {
268  this->matrixPtr->syGetValues(rowind, colind, values);
269  }
277  inline void get_values
278  (std::vector<int> const & rowind,
279  std::vector<int> const & colind,
280  std::vector<Treal> & values,
281  std::vector<int> const & rowPermutation,
282  std::vector<int> const & colPermutation) const {
283  std::vector<int> newRowind;
284  std::vector<int> newColind;
285  this->getPermutedAndSymmetrized(rowind, rowPermutation, newRowind,
286  colind, colPermutation, newColind);
287  this->matrixPtr->syGetValues(newRowind, newColind, values);
288  }
294  inline void get_all_values
295  (std::vector<int> & rowind,
296  std::vector<int> & colind,
297  std::vector<Treal> & values) const {
298  rowind.resize(0);
299  colind.resize(0);
300  values.resize(0);
301  this->matrixPtr->syGetAllValues(rowind, colind, values);
302  }
308  inline void get_all_values
309  (std::vector<int> & rowind,
310  std::vector<int> & colind,
311  std::vector<Treal> & values,
312  std::vector<int> const & rowInversePermutation,
313  std::vector<int> const & colInversePermutation) const {
314  std::vector<int> tmpRowind;
315  std::vector<int> tmpColind;
316  tmpRowind.reserve(rowind.capacity());
317  tmpColind.reserve(colind.capacity());
318  values.resize(0);
319  this->matrixPtr->syGetAllValues(tmpRowind, tmpColind, values);
320  this->getPermutedAndSymmetrized(tmpRowind, rowInversePermutation, rowind,
321  tmpColind, colInversePermutation, colind);
322  }
336  return *this;
337  }
341  this->matrixPtr->nosymToSym();
342  return *this;
343  }
345  *this->matrixPtr = k;
346  return *this;
347  }
348 
349  inline Treal frob() const {
350  return this->matrixPtr->syFrob();
351  }
352  Treal mixed_norm(Treal const requestedAccuracy,
353  int maxIter = -1) const;
354  Treal eucl(Treal const requestedAccuracy,
355  int maxIter = -1) const;
356 
357  void quickEuclBounds(Treal & euclLowerBound,
358  Treal & euclUpperBound) const {
359  Treal frobTmp = frob();
360  euclLowerBound = frobTmp / template_blas_sqrt( (Treal)this->get_nrows() );
361  euclUpperBound = frobTmp;
362  }
363 
364 
370  static Interval<Treal>
373  normType const norm,
374  Treal const requestedAccuracy);
382  static Interval<Treal>
385  normType const norm,
386  Treal const requestedAccuracy,
387  Treal const maxAbsVal);
391  static inline Treal frob_diff
394  return Tmatrix::syFrobDiff(*A.matrixPtr, *B.matrixPtr);
395  }
396 
400  static Treal eucl_diff
403  Treal const requestedAccuracy);
404 
408  static Treal mixed_diff
411  Treal const requestedAccuracy);
412 
422  Treal const requestedAccuracy,
423  Treal const maxAbsVal,
424  VectorType * const eVecPtr = 0);
425 
426 
437  Treal thresh(Treal const threshold,
438  normType const norm);
439 
446  inline Treal frob_thresh(Treal const threshold) {
447  return 2.0 * this->matrixPtr->frob_thresh(threshold / 2);
448  }
449 
450  Treal eucl_thresh(Treal const threshold,
451  MatrixTriangular<Treal, Tmatrix> const * const Zptr = NULL);
452 
453  Treal eucl_element_level_thresh(Treal const threshold);
454 
456  ( SizesAndBlocks & rows_new, SizesAndBlocks & cols_new ) const;
457 
458  Treal mixed_norm_thresh(Treal const threshold);
459 
460  void simple_blockwise_frob_thresh(Treal const threshold) {
461  this->matrixPtr->frobThreshLowestLevel(threshold*threshold, 0);
462  }
463 
464  inline void gersgorin(Treal& lmin, Treal& lmax) const {
465  this->matrixPtr->sy_gersgorin(lmin, lmax);
466  }
467  static inline Treal trace_ab
470  return Tmatrix::sy_trace_ab(*A.matrixPtr, *B.matrixPtr);
471  }
472  inline size_t nnz() const { /* Note: size_t instead of int here to avoid integer overflow. */
473  return this->matrixPtr->sy_nnz();
474  }
475  inline size_t nvalues() const { /* Note: size_t instead of int here to avoid integer overflow. */
476  return this->matrixPtr->sy_nvalues();
477  }
478  inline void write_to_buffer(void* buffer, const int n_bytes) const {
479  this->write_to_buffer_base(buffer, n_bytes, matrix_symm);
480  }
481  inline void read_from_buffer(void* buffer, const int n_bytes) {
482  this->read_from_buffer_base(buffer, n_bytes, matrix_symm);
483  }
484 
485 
491  (const XYZpUV<Treal,
494  Treal,
498  (const XYZ<Treal,
503  (const XYZ<Treal,
508  (const XYZpUV<Treal,
511  Treal,
515  (const XYZ<Treal,
520  (const XYZ<Treal,
523 
531 
536  static void ssmmUpperTriangleOnly(const Treal alpha,
539  const Treal beta,
541 
542 
543  /* Addition */
547  MatrixSymmetric<Treal, Tmatrix> > const & mpm);
551  MatrixSymmetric<Treal, Tmatrix> > const & mm);
558 
562 
566 
567  template<typename Top>
568  Treal accumulateWith(Top & op) {
569  return this->matrixPtr->syAccumulateWith(op);
570  }
571 
572  void random() {
573  this->matrixPtr->syRandom();
574  }
575 
576  void randomZeroStructure(Treal probabilityBeingZero) {
577  this->matrixPtr->syRandomZeroStructure(probabilityBeingZero);
578  }
579 
584  template<typename TRule>
585  void setElementsByRule(TRule & rule) {
586  this->matrixPtr->sySetElementsByRule(rule);
587  return;
588  }
589 
593  ValidPtr<Tmatrix>::swap( this->matrixPtr, dest.matrixPtr );
594  // *this now contains previous content of dest
595  this->clear();
596  }
597 
598  template<typename Tvector>
599  void matVecProd(Tvector & y, Tvector const & x) const {
600  Treal const ONE = 1.0;
601  y = (ONE * (*this) * x);
602  }
603 
604 
605  std::string obj_type_id() const {return "MatrixSymmetric";}
606  protected:
607  inline void writeToFileProt(std::ofstream & file) const {
608  this->writeToFileBase(file, matrix_symm);
609  }
610  inline void readFromFileProt(std::ifstream & file) {
611  this->readFromFileBase(file, matrix_symm);
612  }
613 
619  static void getPermutedAndSymmetrized
620  (std::vector<int> const & rowind,
621  std::vector<int> const & rowPermutation,
622  std::vector<int> & newRowind,
623  std::vector<int> const & colind,
624  std::vector<int> const & colPermutation,
625  std::vector<int> & newColind) {
627  getPermutedIndexes(rowind, rowPermutation, newRowind);
629  getPermutedIndexes(colind, colPermutation, newColind);
630  int tmp;
631  for (unsigned int i = 0; i < newRowind.size(); ++i) {
632  if (newRowind[i] > newColind[i]) {
633  tmp = newRowind[i];
634  newRowind[i] = newColind[i];
635  newColind[i] = tmp;
636  }
637  }
638  }
639  private:
640  };
641 
642  template<typename Treal, typename Tmatrix>
644  mixed_norm(Treal const requestedAccuracy,
645  int maxIter) const {
646  // Construct SizesAndBlocks for frobNormMat
647  SizesAndBlocks rows_new;
648  SizesAndBlocks cols_new;
649  this->getSizesAndBlocksForFrobNormMat( rows_new, cols_new );
650  // Now we can construct an empty matrix where the Frobenius norms
651  // of lowest level nonzero submatrices will be stored
653  frobNormMat.resetSizesAndBlocks(rows_new, cols_new);
654  frobNormMat.getMatrix().syAssignFrobNormsLowestLevel(this->getMatrix());
655  return frobNormMat.eucl(requestedAccuracy, maxIter);
656  }
657 
658 
659  template<typename Treal, typename Tmatrix>
661  eucl(Treal const requestedAccuracy,
662  int maxIter) const {
663  assert(requestedAccuracy >= 0);
664  /* Check if norm is really small, in that case quick return */
665  Treal frobTmp = this->frob();
666  if (frobTmp < requestedAccuracy)
667  return (Treal)0.0;
668  if (maxIter < 0)
669  maxIter = this->get_nrows() * 100;
670  VectorType guess;
671  SizesAndBlocks cols;
672  this->getCols(cols);
673  guess.resetSizesAndBlocks(cols);
674  guess.rand();
675  // Elias note 2010-03-26: changed this back from "new code" to "old code" to reduce memory usage.
676 #if 0 // "new code"
678  Copy.frob_thresh(requestedAccuracy / 2.0);
681  lan(Copy, guess, maxIter);
682  lan.setAbsTol( requestedAccuracy / 2.0 );
683 #else // "old code"
685  <Treal, MatrixSymmetric<Treal, Tmatrix>, VectorType>
686  lan(*this, guess, maxIter);
687  lan.setAbsTol( requestedAccuracy );
688 #endif
689  lan.run();
690  Treal eVal = 0;
691  Treal acc = 0;
692  lan.getLargestMagnitudeEig(eVal, acc);
693  return template_blas_fabs(eVal);
694  }
695 
696  template<typename Treal, typename Tmatrix>
700  normType const norm, Treal const requestedAccuracy) {
701  Treal diff;
702  Treal eNMin;
703  switch (norm) {
704  case frobNorm:
705  diff = frob_diff(A, B);
706  return Interval<Treal>(diff / template_blas_sqrt((Treal)A.get_nrows()), diff);
707  break;
708  case euclNorm:
709  diff = eucl_diff(A, B, requestedAccuracy);
710  eNMin = diff - requestedAccuracy;
711  eNMin = eNMin >= 0 ? eNMin : 0;
712  return Interval<Treal>(eNMin, diff + requestedAccuracy);
713  break;
714  default:
715  throw Failure("MatrixSymmetric<Treal, Tmatrix>::"
716  "diff(const MatrixSymmetric<Treal, Tmatrix>&, "
717  "const MatrixSymmetric<Treal, Tmatrix>&, "
718  "normType const, Treal): "
719  "Diff not implemented for selected norm");
720  }
721  }
722 
723 #if 1
724  template<typename Treal, typename Tmatrix>
728  normType const norm,
729  Treal const requestedAccuracy,
730  Treal const maxAbsVal) {
731  Treal diff;
732  switch (norm) {
733  case frobNorm:
734  {
735  diff = frob_diff(A, B);
736  return Interval<Treal>(diff / template_blas_sqrt((Treal)A.get_nrows()), diff);
737  }
738  break;
739  case euclNorm:
740  return euclDiffIfSmall(A, B, requestedAccuracy, maxAbsVal);
741  break;
742  default:
743  throw Failure("MatrixSymmetric<Treal, Tmatrix>::"
744  "diffIfSmall"
745  "(const MatrixSymmetric<Treal, Tmatrix>&, "
746  "const MatrixSymmetric<Treal, Tmatrix>&, "
747  "normType const, Treal const, Treal const): "
748  "Diff not implemented for selected norm");
749  }
750  }
751 
752 #endif
753 
754 
755  template<typename Treal, typename Tmatrix>
759  Treal const requestedAccuracy) {
760  // DiffMatrix is a lightweight proxy object:
762  Treal maxAbsVal = 2 * frob_diff(A,B);
763  // Now, maxAbsVal should be larger than the Eucl norm
764  // Note that mat::euclIfSmall lies outside this class
765  Treal relTol = sqrt(sqrt(std::numeric_limits<Treal>::epsilon()));
766  Interval<Treal> euclInt =
767  mat::euclIfSmall(Diff, requestedAccuracy, relTol, maxAbsVal);
768  return euclInt.midPoint();
769  }
770 
771  template<typename Treal, typename Tmatrix>
775  Treal const requestedAccuracy) {
777  {
778  SizesAndBlocks rows_new;
779  SizesAndBlocks cols_new;
780  A.getSizesAndBlocksForFrobNormMat( rows_new, cols_new );
781  frobNormMat.resetSizesAndBlocks(rows_new, cols_new);
782  frobNormMat.getMatrix().syAssignDiffFrobNormsLowestLevel(A.getMatrix(),B.getMatrix());
783  }
784  return frobNormMat.eucl(requestedAccuracy);
785  }
786 
787 #if 1
788  template<typename Treal, typename Tmatrix>
792  Treal const requestedAccuracy,
793  Treal const maxAbsVal,
794  VectorType * const eVecPtr) {
795  // DiffMatrix is a lightweight proxy object:
797  // Note that this function lies outside this class
798  Treal relTol = sqrt(sqrt(std::numeric_limits<Treal>::epsilon()));
799  Interval<Treal> tmpInterval = mat::euclIfSmall(Diff, requestedAccuracy, relTol, maxAbsVal, eVecPtr);
800  // Emanuel note: Ugly fix to make certain tests pass, we expand
801  // the interval up to the requested accuracy. Note that larger
802  // intervals may occur if the norm is not 'small'. It happens that
803  // Lanczos misconverges to for example the second largest
804  // eigenvalue. This happens in particular when the first and second
805  // eigenvalues are very close (of the order of the requested
806  // accuracy). Expanding the interval makes the largest eigenvalue
807  // (at least for certain cases) end up inside the interval even
808  // though Lanczos has misconverged.
809  if ( tmpInterval.length() < 2*requestedAccuracy )
810  return Interval<Treal>( tmpInterval.midPoint()-requestedAccuracy, tmpInterval.midPoint()+requestedAccuracy );
811  return tmpInterval;
812  }
813 
814 #endif
815 
816 
817  template<typename Treal, typename Tmatrix>
819  thresh(Treal const threshold,
820  normType const norm) {
821  switch (norm) {
822  case frobNorm:
823  return this->frob_thresh(threshold);
824  break;
825  case euclNorm:
826  return this->eucl_thresh(threshold);
827  break;
828  case mixedNorm:
829  return this->mixed_norm_thresh(threshold);
830  break;
831  default:
832  throw Failure("MatrixSymmetric<Treal, Tmatrix>::"
833  "thresh(Treal const, "
834  "normType const): "
835  "Thresholding not implemented for selected norm");
836  }
837  }
838 
839 #if 1
840 
841  template<typename Treal, typename Tmatrix>
843  eucl_thresh(Treal const threshold,
844  MatrixTriangular<Treal, Tmatrix> const * const Zptr) {
845  if ( Zptr == NULL ) {
847  return TruncObj.run( threshold );
848  }
850  return TruncObj.run( threshold );
851  }
852 
853 #endif
854 
855 
856  template<typename Treal, typename Tmatrix>
858  ( SizesAndBlocks & rows_new, SizesAndBlocks & cols_new ) const {
859  std::vector<int> rows_block_sizes;
860  std::vector<int> cols_block_sizes;
861  int n_rows;
862  int n_cols;
863  {
864  SizesAndBlocks rows;
865  SizesAndBlocks cols;
866  this->getRows(rows);
867  this->getCols(cols);
868  rows.getBlockSizeVector( rows_block_sizes );
869  cols.getBlockSizeVector( cols_block_sizes );
870  rows_block_sizes.pop_back(); // Remove the '1' at the end
871  cols_block_sizes.pop_back(); // Remove the '1' at the end
872  n_rows = rows.getNTotalScalars();
873  n_cols = cols.getNTotalScalars();
874  int factor_rows = rows_block_sizes[rows_block_sizes.size()-1];
875  int factor_cols = cols_block_sizes[cols_block_sizes.size()-1];
876  for (unsigned int ind = 0; ind < rows_block_sizes.size(); ++ind)
877  rows_block_sizes[ind] = rows_block_sizes[ind] / factor_rows;
878  for (unsigned int ind = 0; ind < cols_block_sizes.size(); ++ind)
879  cols_block_sizes[ind] = cols_block_sizes[ind] / factor_cols;
880  // Now set the number of (scalar) rows and cols, should be equal
881  // to the number of blocks at the lowest level of the original
882  // matrix
883  if (n_rows % factor_rows)
884  n_rows = n_rows / factor_rows + 1;
885  else
886  n_rows = n_rows / factor_rows;
887  if (n_cols % factor_cols)
888  n_cols = n_cols / factor_cols + 1;
889  else
890  n_cols = n_cols / factor_cols;
891  }
892  rows_new = SizesAndBlocks( rows_block_sizes, n_rows );
893  cols_new = SizesAndBlocks( cols_block_sizes, n_cols );
894  }
895 
896 
897  template<typename Treal, typename Tmatrix>
899  mixed_norm_thresh(Treal const threshold) {
900  assert(threshold >= (Treal)0.0);
901  if (threshold == (Treal)0.0)
902  return (Treal)0;
903  // Construct SizesAndBlocks for frobNormMat
904  SizesAndBlocks rows_new;
905  SizesAndBlocks cols_new;
906  this->getSizesAndBlocksForFrobNormMat( rows_new, cols_new );
907  // Now we can construct an empty matrix where the Frobenius norms
908  // of lowest level nonzero submatrices will be stored
910  frobNormMat.resetSizesAndBlocks(rows_new, cols_new);
911  // We want the following step to dominate the mixed_norm truncation (this
912  // is where Frobenius norms of submatrices are computed, i.e. it
913  // is here we loop over all matrix elements.)
914  frobNormMat.getMatrix().syAssignFrobNormsLowestLevel(this->getMatrix());
916  Treal mixed_norm_result = TruncObj.run( threshold );
917  frobNormMat.getMatrix().truncateAccordingToSparsityPattern(this->getMatrix());
918  return mixed_norm_result;
919  }
920 
921 
922  /* B = alpha * A */
923  template<typename Treal, typename Tmatrix>
927  assert(!sm.tB);
928  /* Data structure set by assign - therefore set haveDataStructure to true */
929  this->matrixPtr.haveDataStructureSet(true);
930  this->matrixPtr->assign(sm.A, *sm.B.matrixPtr);
931  return *this;
932  }
933  /* C = alpha * A * A + beta * C : A and C are symmetric */
934  template<typename Treal, typename Tmatrix>
937  (const XYZpUV<Treal,
940  Treal,
942  assert(this != &sm2psm.B);
943  if (this == &sm2psm.E && &sm2psm.B == &sm2psm.C) {
944  /* Operation is C = alpha * A * A + beta * C */
945  Tmatrix::sysq('U',
946  sm2psm.A, *sm2psm.B.matrixPtr,
947  sm2psm.D, *this->matrixPtr);
948  return *this;
949  }
950  else /* this != &sm2psm.C || &sm2psm.B != &sm2psm.C */
951  throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator="
952  "(const XYZpUV<Treal, MatrixSymmetric"
953  "<Treal, Tmatrix> >& sm2psm) : "
954  "D = alpha * A * B + beta * C not supported for C != D"
955  " and for symmetric matrices not for A != B since this "
956  "generally will result in a nonsymmetric matrix");
957  }
958 
959  /* C = alpha * A * A : A and C are symmetric */
960  template<typename Treal, typename Tmatrix>
963  (const XYZ<Treal,
966  assert(this != &sm2.B);
967  if (&sm2.B == &sm2.C) {
968  this->matrixPtr.haveDataStructureSet(true);
969  Tmatrix::sysq('U', sm2.A, *sm2.B.matrixPtr, 0, *this->matrixPtr);
970  return *this;
971  }
972  else
973  throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator="
974  "(const XYZ<Treal, MatrixSymmetric<Treal, Tmatrix>,"
975  " MatrixSymmetric<Treal, Tmatrix> >& sm2) : "
976  "Operation C = alpha * A * B with only symmetric "
977  "matrices not supported for A != B");
978  }
979 
980  /* C += alpha * A * A : A and C are symmetric */
981  template<typename Treal, typename Tmatrix>
984  (const XYZ<Treal,
987  assert(this != &sm2.B);
988  if (&sm2.B == &sm2.C) {
989  Tmatrix::sysq('U', sm2.A, *sm2.B.matrixPtr, 1, *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 * transpose(A) + beta * C : C is symmetric */
1001  template<typename Treal, typename Tmatrix>
1004  (const XYZpUV<Treal,
1007  Treal,
1008  MatrixSymmetric<Treal, Tmatrix> >& smmpsm) {
1009  if (this == &smmpsm.E)
1010  if (&smmpsm.B == &smmpsm.C)
1011  if (!smmpsm.tB && smmpsm.tC) {
1012  Tmatrix::syrk('U', false,
1013  smmpsm.A, *smmpsm.B.matrixPtr,
1014  smmpsm.D, *this->matrixPtr);
1015  }
1016  else
1017  throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator="
1018  "(const XYZpUV<Treal, MatrixGeneral"
1019  "<Treal, Tmatrix>, "
1020  "MatrixGeneral<Treal, Tmatrix>, Treal, "
1021  "MatrixSymmetric<Treal, Tmatrix> >&) : "
1022  "C = alpha * A' * A + beta * C, not implemented"
1023  " only C = alpha * A * A' + beta * C");
1024  else
1025  throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator="
1026  "(const XYZpUV<"
1027  "Treal, MatrixGeneral<Treal, Tmatrix>, "
1028  "MatrixGeneral<Treal, Tmatrix>, Treal, "
1029  "MatrixSymmetric<Treal, Tmatrix> >&) : "
1030  "You are trying to call C = alpha * A * A' + beta * C"
1031  " with A and A' being different objects");
1032  else
1033  throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator="
1034  "(const XYZpUV<"
1035  "Treal, MatrixGeneral<Treal, Tmatrix>, "
1036  "MatrixGeneral<Treal, Tmatrix>, Treal, "
1037  "MatrixSymmetric<Treal, Tmatrix> >&) : "
1038  "D = alpha * A * A' + beta * C not supported for C != D");
1039  return *this;
1040  }
1041 
1042  /* C = alpha * A * transpose(A) : C is symmetric */
1043  template<typename Treal, typename Tmatrix>
1046  (const XYZ<Treal,
1049  if (&smm.B == &smm.C)
1050  if (!smm.tB && smm.tC) {
1051  Tmatrix::syrk('U', false,
1052  smm.A, *smm.B.matrixPtr,
1053  0, *this->matrixPtr);
1054  }
1055  else
1056  throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator="
1057  "(const XYZ<"
1058  "Treal, MatrixGeneral<Treal, Tmatrix>, "
1059  "MatrixGeneral<Treal, Tmatrix> >&) : "
1060  "C = alpha * A' * A, not implemented "
1061  "only C = alpha * A * A'");
1062  else
1063  throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator="
1064  "(const XYZ<"
1065  "Treal, MatrixGeneral<Treal, Tmatrix>, "
1066  "MatrixGeneral<Treal, Tmatrix> >&) : "
1067  "You are trying to call C = alpha * A * A' "
1068  "with A and A' being different objects");
1069  return *this;
1070  }
1071 
1072  /* C += alpha * A * transpose(A) : C is symmetric */
1073  template<typename Treal, typename Tmatrix>
1076  (const XYZ<Treal,
1079  if (&smm.B == &smm.C)
1080  if (!smm.tB && smm.tC) {
1081  Tmatrix::syrk('U', false,
1082  smm.A, *smm.B.matrixPtr,
1083  1, *this->matrixPtr);
1084  }
1085  else
1086  throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator+="
1087  "(const XYZ<"
1088  "Treal, MatrixGeneral<Treal, Tmatrix>, "
1089  "MatrixGeneral<Treal, Tmatrix> >&) : "
1090  "C += alpha * A' * A, not implemented "
1091  "only C += alpha * A * A'");
1092  else
1093  throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator+="
1094  "(const XYZ<"
1095  "Treal, MatrixGeneral<Treal, Tmatrix>, "
1096  "MatrixGeneral<Treal, Tmatrix> >&) : "
1097  "You are trying to call C += alpha * A * A' "
1098  "with A and A' being different objects");
1099  return *this;
1100  }
1101 
1102 #if 1
1103  /* A = op1(Z) * A * op2(Z) : Z is upper triangular and A is symmetric */
1104  /* Either op1() or op2() is the transpose operation. */
1105  template<typename Treal, typename Tmatrix>
1111  if (this == &zaz.B) {
1112  if (&zaz.A == &zaz.C) {
1113  if (zaz.tA && !zaz.tC) {
1114  Tmatrix::trsytriplemm('R', *zaz.A.matrixPtr, *this->matrixPtr);
1115  }
1116  else if (!zaz.tA && zaz.tC) {
1117  Tmatrix::trsytriplemm('L', *zaz.A.matrixPtr, *this->matrixPtr);
1118  }
1119  else
1120  throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator="
1121  "(const XYZ<MatrixTriangular<Treal, Tmatrix>,"
1122  "MatrixSymmetric<Treal, Tmatrix>,"
1123  "MatrixTriangular<Treal, Tmatrix> >&) : "
1124  "A = op1(Z) * A * op2(Z) : Either op1 xor op2 must be "
1125  "the transpose operation.");
1126  }
1127  else
1128  throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator="
1129  "(const XYZ<MatrixTriangular<Treal, Tmatrix>,"
1130  "MatrixSymmetric<Treal, Tmatrix>,"
1131  "MatrixTriangular<Treal, Tmatrix> >&) : "
1132  "A = op1(Z1) * A * op2(Z2) : Z1 and Z2 must be the same "
1133  "object");
1134  }
1135  else
1136  throw Failure("MatrixSymmetric<Treal, Tmatrix>::operator="
1137  "(const XYZ<MatrixTriangular<Treal, Tmatrix>,"
1138  "MatrixSymmetric<Treal, Tmatrix>,"
1139  "MatrixTriangular<Treal, Tmatrix> >&) : "
1140  "C = op1(Z) * A * op2(Z) : A and C must be the same "
1141  "object");
1142  return *this;
1143  }
1144 
1145 #endif
1146 
1147 
1152  template<typename Treal, typename Tmatrix>
1154  ssmmUpperTriangleOnly(const Treal alpha,
1157  const Treal beta,
1159  Tmatrix::ssmm_upper_tr_only(alpha, *A.matrixPtr, *B.matrixPtr,
1160  beta, *C.matrixPtr);
1161  }
1162 
1163 
1164 
1165  /* Addition */
1166  /* C = A + B */
1167  template<typename Treal, typename Tmatrix>
1172  assert(this != &mpm.A);
1173  (*this) = mpm.B;
1174  Tmatrix::add(1.0, *mpm.A.matrixPtr, *this->matrixPtr);
1175  return *this;
1176  }
1177  /* C = A - B */
1178  template<typename Treal, typename Tmatrix>
1183  assert(this != &mmm.B);
1184  (*this) = mmm.A;
1185  Tmatrix::add(-1.0, *mmm.B.matrixPtr, *this->matrixPtr);
1186  return *this;
1187  }
1188  /* B += A */
1189  template<typename Treal, typename Tmatrix>
1193  Tmatrix::add(1.0, *A.matrixPtr, *this->matrixPtr);
1194  return *this;
1195  }
1196  /* B -= A */
1197  template<typename Treal, typename Tmatrix>
1201  Tmatrix::add(-1.0, *A.matrixPtr, *this->matrixPtr);
1202  return *this;
1203  }
1204  /* B += alpha * A */
1205  template<typename Treal, typename Tmatrix>
1209  assert(!sm.tB);
1210  Tmatrix::add(sm.A, *sm.B.matrixPtr, *this->matrixPtr);
1211  return *this;
1212  }
1213 
1214  /* B -= alpha * A */
1215  template<typename Treal, typename Tmatrix>
1219  assert(!sm.tB);
1220  Tmatrix::add(-sm.A, *sm.B.matrixPtr, *this->matrixPtr);
1221  return *this;
1222  }
1223 
1228  template<typename Treal, typename Tmatrix, typename Top>
1230  Top & op) {
1231  return A.accumulateWith(op);
1232  }
1233 
1234 
1235 
1236 } /* end namespace mat */
1237 #endif
1238