Go to the documentation of this file.
38 #ifndef MAT_MatrixSymmetric
39 #define MAT_MatrixSymmetric
67 template<
typename Treal,
typename Tmatrix>
68 class MatrixSymmetric :
public MatrixBase<Treal, Tmatrix> {
76 #if __cplusplus >= 201103L
84 :
MatrixBase<Treal, Tmatrix>() { *
this = sm.A * sm.B; }
87 this->matrixPtr->nosymToSym();
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();
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();
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();
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());
124 for (
int col = 0; col < this->get_ncols(); ++col)
125 for (
int row = 0; row < this->get_nrows(); ++row) {
135 this->matrixPtr->nosymToSym();
139 this->matrixPtr->syFullMatrix(fullMat);
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;
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]] =
165 fullMat[colind[ind] + this->get_nrows() * rowind[ind]] =
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);
193 (std::vector<int>
const & rowind,
194 std::vector<int>
const & colind,
195 std::vector<Treal>
const & values,
198 this->resetSizesAndBlocks(newRows, newCols);
199 this->matrixPtr->syAssignFromSparse(rowind, colind, values);
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;
219 colind, colPermutation, newColind);
221 this->matrixPtr->syAssignFromSparse(newRowind, newColind, values);
229 (std::vector<int>
const & rowind,
230 std::vector<int>
const & colind,
231 std::vector<Treal>
const & values,
234 std::vector<int>
const & rowPermutation,
235 std::vector<int>
const & colPermutation) {
236 this->resetSizesAndBlocks(newRows, newCols);
238 rowPermutation, colPermutation);
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);
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;
266 colind, colPermutation, newColind);
267 this->matrixPtr->syAddValues(newRowind, newColind, 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);
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;
294 colind, colPermutation, newColind);
295 this->matrixPtr->syGetValues(newRowind, newColind, values);
303 (std::vector<int> & rowind,
304 std::vector<int> & colind,
305 std::vector<Treal> & values)
const {
309 rowind.reserve(
nnz());
310 colind.reserve(
nnz());
311 values.reserve(
nnz());
312 this->matrixPtr->syGetAllValues(rowind, colind, 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());
330 this->matrixPtr->syGetAllValues(tmpRowind, tmpColind, values);
332 tmpColind, colInversePermutation, colind);
352 this->matrixPtr->nosymToSym();
356 *this->matrixPtr = k;
361 return this->matrixPtr->syFrob();
363 Treal
mixed_norm(Treal
const requestedAccuracy,
364 int maxIter = -1)
const;
365 Treal
eucl(Treal
const requestedAccuracy,
366 int maxIter = -1)
const;
369 Treal & euclUpperBound)
const {
370 Treal frobTmp =
frob();
372 euclUpperBound = frobTmp;
385 Treal
const requestedAccuracy);
397 Treal
const requestedAccuracy,
398 Treal
const maxAbsVal);
405 return Tmatrix::syFrobDiff(*
A.matrixPtr, *
B.matrixPtr);
414 Treal
const requestedAccuracy,
423 Treal
const requestedAccuracy);
434 Treal
const requestedAccuracy,
435 Treal
const maxAbsVal,
449 Treal
thresh(Treal
const threshold,
459 return 2.0 * this->matrixPtr->frob_thresh(threshold / 2);
473 this->matrixPtr->frobThreshLowestLevel(threshold*threshold, 0);
477 this->matrixPtr->sy_gershgorin(lmin, lmax);
482 return Tmatrix::sy_trace_ab(*
A.matrixPtr, *
B.matrixPtr);
484 inline size_t nnz()
const {
485 return this->matrixPtr->sy_nnz();
488 return this->matrixPtr->sy_nvalues();
491 this->write_to_buffer_base(buffer, n_bytes,
matrix_symm);
494 this->read_from_buffer_base(buffer, n_bytes,
matrix_symm);
579 template<
typename Top>
581 return this->matrixPtr->syAccumulateWith(op);
585 this->matrixPtr->syRandom();
589 this->matrixPtr->syRandomZeroStructure(probabilityBeingZero);
596 template<
typename TRule>
598 this->matrixPtr->sySetElementsByRule(rule);
610 template<
typename Tvector>
612 Treal
const ONE = 1.0;
613 y = (ONE * (*this) * x);
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) {
643 for (
unsigned int i = 0; i < newRowind.size(); ++i) {
644 if (newRowind[i] > newColind[i]) {
646 newRowind[i] = newColind[i];
654 template<
typename Treal,
typename Tmatrix>
661 this->getSizesAndBlocksForFrobNormMat( rows_new, cols_new );
665 frobNormMat.resetSizesAndBlocks(rows_new, cols_new);
666 frobNormMat.getMatrix().syAssignFrobNormsLowestLevel(this->getMatrix());
667 return frobNormMat.eucl(requestedAccuracy, maxIter);
671 template<
typename Treal,
typename Tmatrix>
673 eucl(Treal
const requestedAccuracy,
675 assert(requestedAccuracy >= 0);
677 Treal frobTmp = this->frob();
678 if (frobTmp < requestedAccuracy)
681 maxIter = this->get_nrows() * 100;
690 Copy.frob_thresh(requestedAccuracy / 2.0);
693 lan(Copy, guess, maxIter);
694 lan.setAbsTol( requestedAccuracy / 2.0 );
698 lan(*
this, guess, maxIter);
699 lan.setAbsTol( requestedAccuracy );
708 template<
typename Treal,
typename Tmatrix>
712 normType const norm, Treal
const requestedAccuracy) {
717 diff = frob_diff(
A,
B);
721 diff = eucl_diff(
A,
B, requestedAccuracy);
722 eNMin = diff - requestedAccuracy;
723 eNMin = eNMin >= 0 ? eNMin : 0;
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");
736 template<
typename Treal,
typename Tmatrix>
741 Treal
const requestedAccuracy,
742 Treal
const maxAbsVal) {
747 diff = frob_diff(
A,
B);
752 return euclDiffIfSmall(
A,
B, requestedAccuracy, maxAbsVal);
755 throw Failure(
"MatrixSymmetric<Treal, Tmatrix>::"
757 "(const MatrixSymmetric<Treal, Tmatrix>&, "
758 "const MatrixSymmetric<Treal, Tmatrix>&, "
759 "normType const, Treal const, Treal const): "
760 "Diff not implemented for selected norm");
767 template<
typename Treal,
typename Tmatrix>
771 Treal
const requestedAccuracy,
775 Treal maxAbsVal = 2 * frob_diff(
A,
B);
781 mat::euclIfSmall(Diff, requestedAccuracy, relTol, maxAbsVal, eVecPtrNotUsed, maxIter);
785 template<
typename Treal,
typename Tmatrix>
789 Treal
const requestedAccuracy) {
794 A.getSizesAndBlocksForFrobNormMat( rows_new, cols_new );
795 frobNormMat.resetSizesAndBlocks(rows_new, cols_new);
796 frobNormMat.getMatrix().syAssignDiffFrobNormsLowestLevel(
A.getMatrix(),
B.getMatrix());
798 return frobNormMat.eucl(requestedAccuracy);
802 template<
typename Treal,
typename Tmatrix>
806 Treal
const requestedAccuracy,
807 Treal
const maxAbsVal,
823 if ( tmpInterval.
length() < 2*requestedAccuracy )
831 template<
typename Treal,
typename Tmatrix>
833 thresh(Treal
const threshold,
837 return this->frob_thresh(threshold);
840 return this->eucl_thresh(threshold);
843 return this->mixed_norm_thresh(threshold);
846 throw Failure(
"MatrixSymmetric<Treal, Tmatrix>::"
847 "thresh(Treal const, "
849 "Thresholding not implemented for selected norm");
855 template<
typename Treal,
typename Tmatrix>
859 if ( Zptr == NULL ) {
861 return TruncObj.
run( threshold );
864 return TruncObj.
run( threshold );
870 template<
typename Treal,
typename Tmatrix>
873 std::vector<int> rows_block_sizes;
874 std::vector<int> cols_block_sizes;
884 rows_block_sizes.pop_back();
885 cols_block_sizes.pop_back();
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;
897 if (n_rows % factor_rows)
898 n_rows = n_rows / factor_rows + 1;
900 n_rows = n_rows / factor_rows;
901 if (n_cols % factor_cols)
902 n_cols = n_cols / factor_cols + 1;
904 n_cols = n_cols / factor_cols;
911 template<
typename Treal,
typename Tmatrix>
914 assert(threshold >= (Treal)0.0);
915 if (threshold == (Treal)0.0)
920 this->getSizesAndBlocksForFrobNormMat( rows_new, cols_new );
924 frobNormMat.resetSizesAndBlocks(rows_new, cols_new);
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;
937 template<
typename Treal,
typename Tmatrix>
949 this->matrixPtr.haveDataStructureSet(
true);
950 this->matrixPtr->assign(sm.A, *sm.B.matrixPtr);
954 template<
typename Treal,
typename Tmatrix>
962 assert(
this != &sm2psm.B);
963 if (
this == &sm2psm.E && &sm2psm.B == &sm2psm.C) {
966 sm2psm.A, *sm2psm.B.matrixPtr,
967 sm2psm.D, *this->matrixPtr);
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");
980 template<
typename Treal,
typename Tmatrix>
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);
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");
1001 template<
typename Treal,
typename Tmatrix>
1007 assert(
this != &sm2.B);
1008 if (&sm2.B == &sm2.C) {
1009 Tmatrix::sysq(
'U', sm2.A, *sm2.B.matrixPtr, 1, *this->matrixPtr);
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");
1021 template<
typename Treal,
typename Tmatrix>
1029 if (
this == &smmpsm.E)
1030 if (&smmpsm.B == &smmpsm.C)
1031 if (!smmpsm.tB && smmpsm.tC) {
1033 smmpsm.A, *smmpsm.B.matrixPtr,
1034 smmpsm.D, *this->matrixPtr);
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");
1045 throw Failure(
"MatrixSymmetric<Treal, Tmatrix>::operator="
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");
1053 throw Failure(
"MatrixSymmetric<Treal, Tmatrix>::operator="
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");
1063 template<
typename Treal,
typename Tmatrix>
1069 if (&smm.B == &smm.C)
1070 if (!smm.tB && smm.tC) {
1072 smm.A, *smm.B.matrixPtr,
1073 0, *this->matrixPtr);
1076 throw Failure(
"MatrixSymmetric<Treal, Tmatrix>::operator="
1078 "Treal, MatrixGeneral<Treal, Tmatrix>, "
1079 "MatrixGeneral<Treal, Tmatrix> >&) : "
1080 "C = alpha * A' * A, not implemented "
1081 "only C = alpha * A * A'");
1083 throw Failure(
"MatrixSymmetric<Treal, Tmatrix>::operator="
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");
1093 template<
typename Treal,
typename Tmatrix>
1099 if (&smm.B == &smm.C)
1100 if (!smm.tB && smm.tC) {
1102 smm.A, *smm.B.matrixPtr,
1103 1, *this->matrixPtr);
1106 throw Failure(
"MatrixSymmetric<Treal, Tmatrix>::operator+="
1108 "Treal, MatrixGeneral<Treal, Tmatrix>, "
1109 "MatrixGeneral<Treal, Tmatrix> >&) : "
1110 "C += alpha * A' * A, not implemented "
1111 "only C += alpha * A * A'");
1113 throw Failure(
"MatrixSymmetric<Treal, Tmatrix>::operator+="
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");
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);
1136 else if (!zaz.tA && zaz.tC) {
1137 Tmatrix::trsytriplemm(
'L', *zaz.A.matrixPtr, *this->matrixPtr);
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.");
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 "
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 "
1172 template<
typename Treal,
typename Tmatrix>
1179 Tmatrix::ssmm_upper_tr_only(alpha, *
A.matrixPtr, *
B.matrixPtr,
1180 beta, *C.matrixPtr);
1189 template<
typename Treal,
typename Tmatrix>
1194 assert(
this != &mpm.A);
1196 Tmatrix::add(1.0, *mpm.A.matrixPtr, *this->matrixPtr);
1200 template<
typename Treal,
typename Tmatrix>
1205 assert(
this != &mmm.B);
1207 Tmatrix::add(-1.0, *mmm.B.matrixPtr, *this->matrixPtr);
1211 template<
typename Treal,
typename Tmatrix>
1215 Tmatrix::add(1.0, *
A.matrixPtr, *this->matrixPtr);
1219 template<
typename Treal,
typename Tmatrix>
1223 Tmatrix::add(-1.0, *
A.matrixPtr, *this->matrixPtr);
1227 template<
typename Treal,
typename Tmatrix>
1232 Tmatrix::add(sm.A, *sm.B.matrixPtr, *this->matrixPtr);
1237 template<
typename Treal,
typename Tmatrix>
1242 Tmatrix::add(-sm.A, *sm.B.matrixPtr, *this->matrixPtr);
1250 template<
typename Treal,
typename Tmatrix,
typename Top>
1253 return A.accumulateWith(op);
@ euclNorm
Definition: matInclude.h:139
void randomZeroStructure(Treal probabilityBeingZero)
Definition: MatrixSymmetric.h:588
Treal template_blas_sqrt(Treal x)
@ mixedNorm
Definition: matInclude.h:139
MatrixSymmetric(const XY< Treal, MatrixSymmetric< Treal, Tmatrix > > &sm)
Definition: MatrixSymmetric.h:83
size_t nnz() const
Definition: MatrixSymmetric.h:484
This proxy expresses the result of addition of two objects, of possibly different types,...
Definition: matrix_proxy.h:246
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
void writeToFileProt(std::ofstream &file) const
Definition: MatrixSymmetric.h:619
void readFromFileProt(std::ifstream &file)
Definition: MatrixSymmetric.h:622
Treal midPoint() const
Definition: Interval.h:115
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
size_t nvalues() const
Definition: MatrixSymmetric.h:487
Utilities used by other parts of the hierarchical matrix library.
void gershgorin(Treal &lmin, Treal &lmax) const
Definition: MatrixSymmetric.h:476
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
static Treal trace_ab(const MatrixSymmetric< Treal, Tmatrix > &A, const MatrixSymmetric< Treal, Tmatrix > &B)
Definition: MatrixSymmetric.h:480
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
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
void transfer(MatrixSymmetric< Treal, Tmatrix > &dest)
Transfer this matrix to dest, clearing previous content of dest if any.
Definition: MatrixSymmetric.h:604
std::string obj_type_id() const
Definition: MatrixSymmetric.h:617
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
Truncation of symmetric matrices.
Definition: truncation.h:154
void getSizesAndBlocksForFrobNormMat(SizesAndBlocks &rows_new, SizesAndBlocks &cols_new) const
Definition: MatrixSymmetric.h:872
@ frobNorm
Definition: matInclude.h:139
void matVecProd(Tvector &y, Tvector const &x) const
Definition: MatrixSymmetric.h:611
void write_to_buffer(void *buffer, const int n_bytes) const
Definition: MatrixSymmetric.h:490
mat::SizesAndBlocks rows
Definition: test.cc:51
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
This proxy expresses the result of multiplication of three objects, of possibly different types,...
Definition: matrix_proxy.h:67
Treal frob() const
Definition: MatrixSymmetric.h:360
void getBlockSizeVector(std::vector< int > &blockSizesCopy) const
Definition: SizesAndBlocks.cc:87
Treal template_blas_fabs(Treal x)
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
Truncation of symmetric matrices at the element level (used for mixed norm truncation)
Definition: truncation.h:246
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
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
Symmetric matrix.
Definition: MatrixBase.h:51
Upper non-unit triangular matrix.
Definition: MatrixBase.h:53
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 resetSizesAndBlocks(SizesAndBlocks const &newRows)
Definition: VectorGeneral.h:51
Treal length() const
Returns the length of the interval.
Definition: Interval.h:109
Treal eucl(Treal const requestedAccuracy, int maxIter=-1) const
Definition: MatrixSymmetric.h:673
Definition: LanczosLargestMagnitudeEig.h:47
Treal eucl_element_level_thresh(Treal const threshold)
Treal eucl_thresh(Treal const threshold, MatrixTriangular< Treal, Tmatrix > const *const Zptr=NULL)
Definition: MatrixSymmetric.h:857
void getLargestMagnitudeEig(Treal &ev, Treal &accuracy)
Definition: LanczosLargestMagnitudeEig.h:58
void assignFromFull(std::vector< Treal > const &fullMat, std::vector< int > const &rowPermutation, std::vector< int > const &colPermutation)
Definition: MatrixSymmetric.h:117
Treal mixed_norm_thresh(Treal const threshold)
Definition: MatrixSymmetric.h:913
@ matrix_symm
Definition: MatrixBase.h:56
VectorGeneral< Treal, typename Tmatrix::VectorType > VectorType
Definition: MatrixSymmetric.h:70
static void swap(ValidPtr< Tobj > &ptrA, ValidPtr< Tobj > &ptrB)
Definition: ValidPtr.h:106
This proxy expresses the result of substraction of two objects, of possibly different types,...
Definition: matrix_proxy.h:266
Treal real
Definition: MatrixSymmetric.h:71
MatrixSymmetric< Treal, Tmatrix > & operator=(const MatrixSymmetric< Treal, Tmatrix > &symm)
Definition: MatrixSymmetric.h:345
Truncation of symmetric matrices with Z.
Definition: truncation.h:203
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
Definition: MatrixBase.h:55
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
Treal mixed_norm(Treal const requestedAccuracy, int maxIter=-1) const
Definition: MatrixSymmetric.h:656
Normal matrix.
Definition: MatrixBase.h:49
void random()
Definition: MatrixSymmetric.h:584
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
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
mat::SizesAndBlocks cols
Definition: test.cc:52
Definition: allocate.cc:39
Treal run(Treal const threshold)
Definition: truncation.h:80
Treal accumulateWith(Top &op)
Definition: MatrixSymmetric.h:580
normType
Definition: matInclude.h:139
This proxy expresses the result of multiplication of three objects added to two other multiplied obje...
Definition: matrix_proxy.h:88
Definition: Interval.h:46
MatrixSymmetric< Treal, Tmatrix > & operator=(const MatrixGeneral< Treal, Tmatrix > &matr)
Definition: MatrixSymmetric.h:350
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
void rand()
Definition: VectorGeneral.h:108
MatrixSymmetric< Treal, Tmatrix > & operator=(int const k)
Definition: MatrixSymmetric.h:355
MatrixSymmetric(const MatrixGeneral< Treal, Tmatrix > &matr)
'Copy from normal matrix' - constructor
Definition: MatrixSymmetric.h:85
Base class for matrix API.
Definition: MatrixBase.h:69
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 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
MatrixSymmetric(const MatrixSymmetric< Treal, Tmatrix > &symm)
Copy constructor
Definition: MatrixSymmetric.h:80
This proxy expresses the result of multiplication of two objects, of possibly different types,...
Definition: matrix_proxy.h:51
void quickEuclBounds(Treal &euclLowerBound, Treal &euclUpperBound) const
Definition: MatrixSymmetric.h:368
Describes dimensions of matrix and its blocks on all levels.
Definition: SizesAndBlocks.h:45
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
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
static void getPermutedIndexes(std::vector< int > const &index, std::vector< int > const &permutation, std::vector< int > &newIndex)
Definition: MatrixBase.h:205
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
MatrixBase< Treal, Tmatrix > & operator=(const MatrixBase< Treal, Tmatrix > &other)
Definition: MatrixBase.h:166
void assignFromFull(std::vector< Treal > const &fullMat)
Definition: MatrixSymmetric.h:110
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
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
Definition: mat_utils.h:43
int getNTotalScalars() const
Definition: SizesAndBlocks.h:84
void read_from_buffer(void *buffer, const int n_bytes)
Definition: MatrixSymmetric.h:493
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
void fullMatrix(std::vector< Treal > &fullMat) const
Save matrix as full matrix.
Definition: MatrixSymmetric.h:138
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
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()
Default constructor
Definition: MatrixSymmetric.h:73