ergo
MatrixSymmetric.h
Go to the documentation of this file.
1/* Ergo, version 3.8.2, a program for linear scaling electronic structure
2 * calculations.
3 * Copyright (C) 2023 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
49namespace 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
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 }
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 }
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 }
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 }
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 }
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);
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
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);
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
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
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,
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
Interval class.
Class for computing the largest magnitude eigenvalue of a symmetric matrix with the Lanczos method.
Base class for matrix API.
Treal run(Treal const threshold)
Definition truncation.h:80
Truncation of symmetric matrices at the element level (used for mixed norm truncation)
Definition truncation.h:246
Truncation of symmetric matrices with Z.
Definition truncation.h:203
Truncation of symmetric matrices.
Definition truncation.h:154
Definition Failure.h:57
Definition Interval.h:46
Treal midPoint() const
Definition Interval.h:115
Treal length() const
Returns the length of the interval.
Definition Interval.h:109
Base class for matrix API.
Definition MatrixBase.h:69
int get_ncols() const
Definition MatrixBase.h:141
Tmatrix const & getMatrix() const
Definition MatrixBase.h:146
static void getPermutedIndexes(std::vector< int > const &index, std::vector< int > const &permutation, std::vector< int > &newIndex)
Definition MatrixBase.h:205
void read_from_buffer_base(void *buffer, const int n_bytes, const matrix_type mattype)
Definition MatrixBase.h:281
int get_nrows() const
Definition MatrixBase.h:138
void writeToFileBase(std::ofstream &file, matrix_type const mattype) const
Definition MatrixBase.h:221
MatrixBase< Treal, Tmatrix > & operator=(const MatrixBase< Treal, Tmatrix > &other)
Definition MatrixBase.h:166
void resetSizesAndBlocks(SizesAndBlocks const &newRows, SizesAndBlocks const &newCols)
Definition MatrixBase.h:76
void clear()
Release memory for the information written to file.
Definition MatrixBase.h:118
void write_to_buffer_base(void *buffer, const int n_bytes, const matrix_type mattype) const
Definition MatrixBase.h:254
ValidPtr< Tmatrix > matrixPtr
Definition MatrixBase.h:153
void readFromFileBase(std::ifstream &file, matrix_type const mattype)
Definition MatrixBase.h:236
Normal matrix.
Definition MatrixGeneral.h:59
Symmetric matrix.
Definition MatrixSymmetric.h:68
void gershgorin(Treal &lmin, Treal &lmax) const
Definition MatrixSymmetric.h:476
Treal eucl_element_level_thresh(Treal const threshold)
void getSizesAndBlocksForFrobNormMat(SizesAndBlocks &rows_new, SizesAndBlocks &cols_new) const
Definition MatrixSymmetric.h:872
MatrixSymmetric(const MatrixSymmetric< Treal, Tmatrix > &symm)
Copy constructor
Definition MatrixSymmetric.h:80
Treal accumulateWith(Top &op)
Definition MatrixSymmetric.h:580
void write_to_buffer(void *buffer, const int n_bytes) const
Definition MatrixSymmetric.h:490
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
VectorGeneral< Treal, typename Tmatrix::VectorType > VectorType
Definition MatrixSymmetric.h:70
size_t nnz() const
Definition MatrixSymmetric.h:484
void writeToFileProt(std::ofstream &file) const
Write object to file.
Definition MatrixSymmetric.h:619
Treal mixed_norm_thresh(Treal const threshold)
Definition MatrixSymmetric.h:913
void readFromFileProt(std::ifstream &file)
Read object from file.
Definition MatrixSymmetric.h:622
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
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
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
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
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
std::string obj_type_id() const
Definition MatrixSymmetric.h:617
Treal mixed_norm(Treal const requestedAccuracy, int maxIter=-1) const
Definition MatrixSymmetric.h:656
Treal eucl(Treal const requestedAccuracy, int maxIter=-1) const
Definition MatrixSymmetric.h:673
MatrixSymmetric< Treal, Tmatrix > & operator=(const MatrixSymmetric< Treal, Tmatrix > &symm)
Definition MatrixSymmetric.h:345
void randomZeroStructure(Treal probabilityBeingZero)
Definition MatrixSymmetric.h:588
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
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
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
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
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
void assignFromFull(std::vector< Treal > const &fullMat)
Definition MatrixSymmetric.h:110
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
static Treal trace_ab(const MatrixSymmetric< Treal, Tmatrix > &A, const MatrixSymmetric< Treal, Tmatrix > &B)
Definition MatrixSymmetric.h:480
Treal frob() const
Definition MatrixSymmetric.h:360
void read_from_buffer(void *buffer, const int n_bytes)
Definition MatrixSymmetric.h:493
void matVecProd(Tvector &y, Tvector const &x) const
Definition MatrixSymmetric.h:611
void simple_blockwise_frob_thresh(Treal const threshold)
Definition MatrixSymmetric.h:472
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
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
void quickEuclBounds(Treal &euclLowerBound, Treal &euclUpperBound) const
Definition MatrixSymmetric.h:368
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
Treal real
Definition MatrixSymmetric.h:71
void assignFromFull(std::vector< Treal > const &fullMat, std::vector< int > const &rowPermutation, std::vector< int > const &colPermutation)
Definition MatrixSymmetric.h:117
MatrixSymmetric(const MatrixGeneral< Treal, Tmatrix > &matr)
'Copy from normal matrix' - constructor
Definition MatrixSymmetric.h:85
MatrixSymmetric(const XY< Treal, MatrixSymmetric< Treal, Tmatrix > > &sm)
Definition MatrixSymmetric.h:83
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
Treal eucl_thresh(Treal const threshold, MatrixTriangular< Treal, Tmatrix > const *const Zptr=NULL)
Definition MatrixSymmetric.h:857
void random()
Definition MatrixSymmetric.h:584
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
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
MatrixSymmetric()
Default constructor
Definition MatrixSymmetric.h:73
void transfer(MatrixSymmetric< Treal, Tmatrix > &dest)
Transfer this matrix to dest, clearing previous content of dest if any.
Definition MatrixSymmetric.h:604
MatrixSymmetric< Treal, Tmatrix > & operator=(int const k)
Definition MatrixSymmetric.h:355
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
size_t nvalues() const
Definition MatrixSymmetric.h:487
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
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
MatrixSymmetric< Treal, Tmatrix > & operator=(const MatrixGeneral< Treal, Tmatrix > &matr)
Definition MatrixSymmetric.h:350
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
void fullMatrix(std::vector< Treal > &fullMat) const
Save matrix as full matrix.
Definition MatrixSymmetric.h:138
Upper non-unit triangular matrix.
Definition MatrixTriangular.h:59
Describes dimensions of matrix and its blocks on all levels.
Definition SizesAndBlocks.h:45
int getNTotalScalars() const
Definition SizesAndBlocks.h:84
void getBlockSizeVector(std::vector< int > &blockSizesCopy) const
Definition SizesAndBlocks.cc:87
static void swap(ValidPtr< Tobj > &ptrA, ValidPtr< Tobj > &ptrB)
Definition ValidPtr.h:106
Definition VectorGeneral.h:48
void resetSizesAndBlocks(SizesAndBlocks const &newRows)
Definition VectorGeneral.h:51
void rand()
Definition VectorGeneral.h:108
Definition LanczosLargestMagnitudeEig.h:47
void getLargestMagnitudeEig(Treal &ev, Treal &accuracy)
Definition LanczosLargestMagnitudeEig.h:58
mat::SizesAndBlocks rows
Definition test.cc:51
mat::SizesAndBlocks cols
Definition test.cc:52
#define B
#define A
Utilities used by other parts of the hierarchical matrix library.
Definition allocate.cc:39
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
@ matrix_symm
Definition MatrixBase.h:56
normType
Definition matInclude.h:139
@ euclNorm
Definition matInclude.h:139
@ frobNorm
Definition matInclude.h:139
@ mixedNorm
Definition matInclude.h:139
static Treal getMachineEpsilon()
Definition matInclude.h:147
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
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
Definition mat_utils.h:60
This proxy expresses the result of multiplication of three objects, of possibly different types,...
Definition matrix_proxy.h:67
This proxy expresses the result of multiplication of three objects added to two other multiplied obje...
Definition matrix_proxy.h:88
This proxy expresses the result of multiplication of two objects, of possibly different types,...
Definition matrix_proxy.h:51
This proxy expresses the result of substraction of two objects, of possibly different types,...
Definition matrix_proxy.h:266
This proxy expresses the result of addition of two objects, of possibly different types,...
Definition matrix_proxy.h:246
Treal template_blas_sqrt(Treal x)
Treal template_blas_fabs(Treal x)
Classes for truncation of small matrix elements.