38 #ifndef MAT_MatrixSymmetric 39 #define MAT_MatrixSymmetric 67 template<
typename Treal,
typename Tmatrix>
68 class MatrixSymmetric :
public MatrixBase<Treal, Tmatrix> {
78 :
MatrixBase<Treal, Tmatrix>() { *
this = sm.A * sm.B; }
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);
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);
104 (std::vector<Treal>
const & fullMat) {
106 this->
matrixPtr->assignFromFull(fullMat);
111 (std::vector<Treal>
const & fullMat,
112 std::vector<int>
const & rowPermutation,
113 std::vector<int>
const & colPermutation) {
115 std::vector<int> rowind(fullMat.size());
116 std::vector<int> colind(fullMat.size());
118 for (
int col = 0; col < this->
get_ncols(); ++col)
119 for (
int row = 0; row < this->
get_nrows(); ++row) {
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;
149 rowInversePermutation,
150 colInversePermutation);
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] <
157 fullMat[rowind[ind] + this->
get_nrows() * colind[ind]] =
159 fullMat[colind[ind] + this->
get_nrows() * rowind[ind]] =
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);
187 (std::vector<int>
const & rowind,
188 std::vector<int>
const & colind,
189 std::vector<Treal>
const & values,
193 this->
matrixPtr->syAssignFromSparse(rowind, colind, values);
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;
213 colind, colPermutation, newColind);
215 this->
matrixPtr->syAssignFromSparse(newRowind, newColind, values);
223 (std::vector<int>
const & rowind,
224 std::vector<int>
const & colind,
225 std::vector<Treal>
const & values,
228 std::vector<int>
const & rowPermutation,
229 std::vector<int>
const & colPermutation) {
232 rowPermutation, colPermutation);
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);
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;
260 colind, colPermutation, newColind);
261 this->
matrixPtr->syAddValues(newRowind, newColind, 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);
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;
288 colind, colPermutation, newColind);
289 this->
matrixPtr->syGetValues(newRowind, newColind, values);
297 (std::vector<int> & rowind,
298 std::vector<int> & colind,
299 std::vector<Treal> & values)
const {
303 rowind.reserve(
nnz());
304 colind.reserve(
nnz());
305 values.reserve(
nnz());
306 this->
matrixPtr->syGetAllValues(rowind, colind, 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());
324 this->
matrixPtr->syGetAllValues(tmpRowind, tmpColind, values);
326 tmpColind, colInversePermutation, colind);
357 Treal
mixed_norm(Treal
const requestedAccuracy,
358 int maxIter = -1)
const;
359 Treal
eucl(Treal
const requestedAccuracy,
360 int maxIter = -1)
const;
363 Treal & euclUpperBound)
const {
364 Treal frobTmp =
frob();
366 euclUpperBound = frobTmp;
379 Treal
const requestedAccuracy);
391 Treal
const requestedAccuracy,
392 Treal
const maxAbsVal);
399 return Tmatrix::syFrobDiff(*
A.matrixPtr, *
B.matrixPtr);
408 Treal
const requestedAccuracy,
417 Treal
const requestedAccuracy);
428 Treal
const requestedAccuracy,
429 Treal
const maxAbsVal,
443 Treal
thresh(Treal
const threshold,
453 return 2.0 * this->
matrixPtr->frob_thresh(threshold / 2);
467 this->
matrixPtr->frobThreshLowestLevel(threshold*threshold, 0);
471 this->
matrixPtr->sy_gershgorin(lmin, lmax);
476 return Tmatrix::sy_trace_ab(*
A.matrixPtr, *
B.matrixPtr);
478 inline size_t nnz()
const {
573 template<
typename Top>
575 return this->
matrixPtr->syAccumulateWith(op);
583 this->
matrixPtr->syRandomZeroStructure(probabilityBeingZero);
590 template<
typename TRule>
592 this->
matrixPtr->sySetElementsByRule(rule);
604 template<
typename Tvector>
606 Treal
const ONE = 1.0;
607 y = (ONE * (*this) * x);
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) {
637 for (
unsigned int i = 0; i < newRowind.size(); ++i) {
638 if (newRowind[i] > newColind[i]) {
640 newRowind[i] = newColind[i];
648 template<
typename Treal,
typename Tmatrix>
655 this->getSizesAndBlocksForFrobNormMat( rows_new, cols_new );
659 frobNormMat.resetSizesAndBlocks(rows_new, cols_new);
660 frobNormMat.getMatrix().syAssignFrobNormsLowestLevel(this->getMatrix());
661 return frobNormMat.eucl(requestedAccuracy, maxIter);
665 template<
typename Treal,
typename Tmatrix>
667 eucl(Treal
const requestedAccuracy,
669 assert(requestedAccuracy >= 0);
671 Treal frobTmp = this->frob();
672 if (frobTmp < requestedAccuracy)
675 maxIter = this->get_nrows() * 100;
684 Copy.frob_thresh(requestedAccuracy / 2.0);
687 lan(Copy, guess, maxIter);
688 lan.setAbsTol( requestedAccuracy / 2.0 );
692 lan(*
this, guess, maxIter);
693 lan.setAbsTol( requestedAccuracy );
702 template<
typename Treal,
typename Tmatrix>
706 normType const norm, Treal
const requestedAccuracy) {
711 diff = frob_diff(
A,
B);
715 diff = eucl_diff(
A,
B, requestedAccuracy);
716 eNMin = diff - requestedAccuracy;
717 eNMin = eNMin >= 0 ? eNMin : 0;
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");
730 template<
typename Treal,
typename Tmatrix>
735 Treal
const requestedAccuracy,
736 Treal
const maxAbsVal) {
741 diff = frob_diff(
A,
B);
746 return euclDiffIfSmall(
A,
B, requestedAccuracy, maxAbsVal);
749 throw Failure(
"MatrixSymmetric<Treal, Tmatrix>::" 751 "(const MatrixSymmetric<Treal, Tmatrix>&, " 752 "const MatrixSymmetric<Treal, Tmatrix>&, " 753 "normType const, Treal const, Treal const): " 754 "Diff not implemented for selected norm");
761 template<
typename Treal,
typename Tmatrix>
765 Treal
const requestedAccuracy,
769 Treal maxAbsVal = 2 * frob_diff(
A,
B);
775 mat::euclIfSmall(Diff, requestedAccuracy, relTol, maxAbsVal, eVecPtrNotUsed, maxIter);
779 template<
typename Treal,
typename Tmatrix>
783 Treal
const requestedAccuracy) {
788 A.getSizesAndBlocksForFrobNormMat( rows_new, cols_new );
789 frobNormMat.resetSizesAndBlocks(rows_new, cols_new);
790 frobNormMat.getMatrix().syAssignDiffFrobNormsLowestLevel(
A.getMatrix(),
B.getMatrix());
792 return frobNormMat.eucl(requestedAccuracy);
796 template<
typename Treal,
typename Tmatrix>
800 Treal
const requestedAccuracy,
801 Treal
const maxAbsVal,
817 if ( tmpInterval.
length() < 2*requestedAccuracy )
825 template<
typename Treal,
typename Tmatrix>
831 return this->frob_thresh(threshold);
834 return this->eucl_thresh(threshold);
837 return this->mixed_norm_thresh(threshold);
840 throw Failure(
"MatrixSymmetric<Treal, Tmatrix>::" 841 "thresh(Treal const, " 843 "Thresholding not implemented for selected norm");
849 template<
typename Treal,
typename Tmatrix>
853 if ( Zptr == NULL ) {
855 return TruncObj.
run( threshold );
858 return TruncObj.
run( threshold );
864 template<
typename Treal,
typename Tmatrix>
867 std::vector<int> rows_block_sizes;
868 std::vector<int> cols_block_sizes;
878 rows_block_sizes.pop_back();
879 cols_block_sizes.pop_back();
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;
891 if (n_rows % factor_rows)
892 n_rows = n_rows / factor_rows + 1;
894 n_rows = n_rows / factor_rows;
895 if (n_cols % factor_cols)
896 n_cols = n_cols / factor_cols + 1;
898 n_cols = n_cols / factor_cols;
905 template<
typename Treal,
typename Tmatrix>
908 assert(threshold >= (Treal)0.0);
909 if (threshold == (Treal)0.0)
914 this->getSizesAndBlocksForFrobNormMat( rows_new, cols_new );
918 frobNormMat.resetSizesAndBlocks(rows_new, cols_new);
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;
931 template<
typename Treal,
typename Tmatrix>
943 this->matrixPtr.haveDataStructureSet(
true);
944 this->matrixPtr->assign(sm.A, *sm.B.matrixPtr);
948 template<
typename Treal,
typename Tmatrix>
956 assert(
this != &sm2psm.B);
957 if (
this == &sm2psm.E && &sm2psm.B == &sm2psm.C) {
960 sm2psm.A, *sm2psm.B.matrixPtr,
961 sm2psm.D, *this->matrixPtr);
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");
974 template<
typename Treal,
typename Tmatrix>
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);
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");
995 template<
typename Treal,
typename Tmatrix>
1001 assert(
this != &sm2.B);
1002 if (&sm2.B == &sm2.C) {
1003 Tmatrix::sysq(
'U', sm2.A, *sm2.B.matrixPtr, 1, *this->matrixPtr);
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");
1015 template<
typename Treal,
typename Tmatrix>
1023 if (
this == &smmpsm.E)
1024 if (&smmpsm.B == &smmpsm.C)
1025 if (!smmpsm.tB && smmpsm.tC) {
1027 smmpsm.A, *smmpsm.B.matrixPtr,
1028 smmpsm.D, *this->matrixPtr);
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");
1039 throw Failure(
"MatrixSymmetric<Treal, Tmatrix>::operator=" 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");
1047 throw Failure(
"MatrixSymmetric<Treal, Tmatrix>::operator=" 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");
1057 template<
typename Treal,
typename Tmatrix>
1063 if (&smm.B == &smm.C)
1064 if (!smm.tB && smm.tC) {
1066 smm.A, *smm.B.matrixPtr,
1067 0, *this->matrixPtr);
1070 throw Failure(
"MatrixSymmetric<Treal, Tmatrix>::operator=" 1072 "Treal, MatrixGeneral<Treal, Tmatrix>, " 1073 "MatrixGeneral<Treal, Tmatrix> >&) : " 1074 "C = alpha * A' * A, not implemented " 1075 "only C = alpha * A * A'");
1077 throw Failure(
"MatrixSymmetric<Treal, Tmatrix>::operator=" 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");
1087 template<
typename Treal,
typename Tmatrix>
1093 if (&smm.B == &smm.C)
1094 if (!smm.tB && smm.tC) {
1096 smm.A, *smm.B.matrixPtr,
1097 1, *this->matrixPtr);
1100 throw Failure(
"MatrixSymmetric<Treal, Tmatrix>::operator+=" 1102 "Treal, MatrixGeneral<Treal, Tmatrix>, " 1103 "MatrixGeneral<Treal, Tmatrix> >&) : " 1104 "C += alpha * A' * A, not implemented " 1105 "only C += alpha * A * A'");
1107 throw Failure(
"MatrixSymmetric<Treal, Tmatrix>::operator+=" 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");
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);
1130 else if (!zaz.tA && zaz.tC) {
1131 Tmatrix::trsytriplemm(
'L', *zaz.A.matrixPtr, *this->matrixPtr);
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.");
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 " 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 " 1166 template<
typename Treal,
typename Tmatrix>
1173 Tmatrix::ssmm_upper_tr_only(alpha, *
A.matrixPtr, *
B.matrixPtr,
1174 beta, *C.matrixPtr);
1183 template<
typename Treal,
typename Tmatrix>
1188 assert(
this != &mpm.A);
1190 Tmatrix::add(1.0, *mpm.A.matrixPtr, *this->matrixPtr);
1194 template<
typename Treal,
typename Tmatrix>
1199 assert(
this != &mmm.B);
1201 Tmatrix::add(-1.0, *mmm.B.matrixPtr, *this->matrixPtr);
1205 template<
typename Treal,
typename Tmatrix>
1209 Tmatrix::add(1.0, *
A.matrixPtr, *this->matrixPtr);
1213 template<
typename Treal,
typename Tmatrix>
1217 Tmatrix::add(-1.0, *
A.matrixPtr, *this->matrixPtr);
1221 template<
typename Treal,
typename Tmatrix>
1226 Tmatrix::add(sm.A, *sm.B.matrixPtr, *this->matrixPtr);
1231 template<
typename Treal,
typename Tmatrix>
1236 Tmatrix::add(-sm.A, *sm.B.matrixPtr, *this->matrixPtr);
1244 template<
typename Treal,
typename Tmatrix,
typename Top>
1247 return A.accumulateWith(op);
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
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)
'Copy from normal matrix' - 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 'mixed' 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
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 'op' on all nonzero matrix elements and sums up the result and return...
Definition: MatrixSymmetric.h:1245
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
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