36 #ifndef MAT_MatrixSymmetric 37 #define MAT_MatrixSymmetric 65 template<
typename Treal,
typename Tmatrix>
66 class MatrixSymmetric :
public MatrixBase<Treal, Tmatrix> {
76 :
MatrixBase<Treal, Tmatrix>() { *
this = sm.A * sm.B; }
84 template<
typename Tfull>
85 inline void assign_from_full
86 (Tfull
const*
const fullmatrix,
87 int const totnrows,
int const totncols) {
88 assert(totnrows == totncols);
89 this->
matrixPtr->assign_from_full(fullmatrix, totnrows, totncols);
92 inline void assign_from_full
93 (Treal
const*
const fullmatrix,
94 int const totnrows,
int const totncols) {
95 assert(totnrows == totncols);
96 this->
matrixPtr->assign_from_full(fullmatrix, totnrows, totncols);
102 (std::vector<Treal>
const & fullMat) {
104 this->
matrixPtr->assignFromFull(fullMat);
109 (std::vector<Treal>
const & fullMat,
110 std::vector<int>
const & rowPermutation,
111 std::vector<int>
const & colPermutation) {
113 std::vector<int> rowind(fullMat.size());
114 std::vector<int> colind(fullMat.size());
116 for (
int col = 0; col < this->
get_ncols(); ++col)
117 for (
int row = 0; row < this->
get_nrows(); ++row) {
140 (std::vector<Treal> & fullMat,
141 std::vector<int>
const & rowInversePermutation,
142 std::vector<int>
const & colInversePermutation)
const {
143 std::vector<int> rowind;
144 std::vector<int> colind;
145 std::vector<Treal> values;
147 rowInversePermutation,
148 colInversePermutation);
150 assert(rowind.size() == values.size());
151 assert(rowind.size() == colind.size());
152 for (
unsigned int ind = 0; ind < values.size(); ++ind) {
153 assert(rowind[ind] + this->
get_nrows() * colind[ind] <
155 fullMat[rowind[ind] + this->
get_nrows() * colind[ind]] =
157 fullMat[colind[ind] + this->
get_nrows() * rowind[ind]] =
168 (std::vector<int>
const & rowind,
169 std::vector<int>
const & colind,
170 std::vector<Treal>
const & values) {
171 this->
matrixPtr->syAssignFromSparse(rowind, colind, values);
185 (std::vector<int>
const & rowind,
186 std::vector<int>
const & colind,
187 std::vector<Treal>
const & values,
191 this->
matrixPtr->syAssignFromSparse(rowind, colind, values);
203 (std::vector<int>
const & rowind,
204 std::vector<int>
const & colind,
205 std::vector<Treal>
const & values,
206 std::vector<int>
const & rowPermutation,
207 std::vector<int>
const & colPermutation) {
208 std::vector<int> newRowind;
209 std::vector<int> newColind;
211 colind, colPermutation, newColind);
213 this->
matrixPtr->syAssignFromSparse(newRowind, newColind, values);
221 (std::vector<int>
const & rowind,
222 std::vector<int>
const & colind,
223 std::vector<Treal>
const & values,
226 std::vector<int>
const & rowPermutation,
227 std::vector<int>
const & colPermutation) {
230 rowPermutation, colPermutation);
240 (std::vector<int>
const & rowind,
241 std::vector<int>
const & colind,
242 std::vector<Treal>
const & values) {
243 this->
matrixPtr->syAddValues(rowind, colind, values);
250 (std::vector<int>
const & rowind,
251 std::vector<int>
const & colind,
252 std::vector<Treal>
const & values,
253 std::vector<int>
const & rowPermutation,
254 std::vector<int>
const & colPermutation) {
255 std::vector<int> newRowind;
256 std::vector<int> newColind;
258 colind, colPermutation, newColind);
259 this->
matrixPtr->syAddValues(newRowind, newColind, values);
265 (std::vector<int>
const & rowind,
266 std::vector<int>
const & colind,
267 std::vector<Treal> & values)
const {
268 this->
matrixPtr->syGetValues(rowind, colind, values);
278 (std::vector<int>
const & rowind,
279 std::vector<int>
const & colind,
280 std::vector<Treal> & values,
281 std::vector<int>
const & rowPermutation,
282 std::vector<int>
const & colPermutation)
const {
283 std::vector<int> newRowind;
284 std::vector<int> newColind;
286 colind, colPermutation, newColind);
287 this->
matrixPtr->syGetValues(newRowind, newColind, values);
295 (std::vector<int> & rowind,
296 std::vector<int> & colind,
297 std::vector<Treal> & values)
const {
301 this->
matrixPtr->syGetAllValues(rowind, colind, values);
309 (std::vector<int> & rowind,
310 std::vector<int> & colind,
311 std::vector<Treal> & values,
312 std::vector<int>
const & rowInversePermutation,
313 std::vector<int>
const & colInversePermutation)
const {
314 std::vector<int> tmpRowind;
315 std::vector<int> tmpColind;
316 tmpRowind.reserve(rowind.capacity());
317 tmpColind.reserve(colind.capacity());
319 this->
matrixPtr->syGetAllValues(tmpRowind, tmpColind, values);
321 tmpColind, colInversePermutation, colind);
352 Treal
mixed_norm(Treal
const requestedAccuracy,
353 int maxIter = -1)
const;
354 Treal
eucl(Treal
const requestedAccuracy,
355 int maxIter = -1)
const;
358 Treal & euclUpperBound)
const {
359 Treal frobTmp =
frob();
361 euclUpperBound = frobTmp;
374 Treal
const requestedAccuracy);
386 Treal
const requestedAccuracy,
387 Treal
const maxAbsVal);
394 return Tmatrix::syFrobDiff(*A.matrixPtr, *B.matrixPtr);
403 Treal
const requestedAccuracy);
411 Treal
const requestedAccuracy);
422 Treal
const requestedAccuracy,
423 Treal
const maxAbsVal,
424 VectorType *
const eVecPtr = 0);
437 Treal
thresh(Treal
const threshold,
447 return 2.0 * this->
matrixPtr->frob_thresh(threshold / 2);
461 this->
matrixPtr->frobThreshLowestLevel(threshold*threshold, 0);
465 this->
matrixPtr->sy_gersgorin(lmin, lmax);
470 return Tmatrix::sy_trace_ab(*A.matrixPtr, *B.matrixPtr);
472 inline size_t nnz()
const {
493 MatrixSymmetric<Treal, Tmatrix>,
495 MatrixSymmetric<Treal, Tmatrix> >& sm2psm);
497 MatrixSymmetric<Treal, Tmatrix>&
operator=
499 MatrixSymmetric<Treal, Tmatrix>,
500 MatrixSymmetric<Treal, Tmatrix> >& sm2);
502 MatrixSymmetric<Treal, Tmatrix>&
operator+=
504 MatrixSymmetric<Treal, Tmatrix>,
505 MatrixSymmetric<Treal, Tmatrix> >& sm2);
507 MatrixSymmetric<Treal, Tmatrix>&
operator=
510 MatrixGeneral<Treal, Tmatrix>,
512 MatrixSymmetric<Treal, Tmatrix> >& smmpsm);
514 MatrixSymmetric<Treal, Tmatrix>&
operator=
516 MatrixGeneral<Treal, Tmatrix>,
517 MatrixGeneral<Treal, Tmatrix> >& smm);
519 MatrixSymmetric<Treal, Tmatrix>&
operator+=
521 MatrixGeneral<Treal, Tmatrix>,
522 MatrixGeneral<Treal, Tmatrix> >& smm);
527 MatrixSymmetric<Treal, Tmatrix>&
operator=
529 MatrixSymmetric<Treal, Tmatrix>,
537 const MatrixSymmetric<Treal, Tmatrix>& A,
538 const MatrixSymmetric<Treal, Tmatrix>& B,
540 MatrixSymmetric<Treal, Tmatrix>& C);
545 MatrixSymmetric<Treal, Tmatrix>&
operator=
547 MatrixSymmetric<Treal, Tmatrix> >
const & mpm);
549 MatrixSymmetric<Treal, Tmatrix>&
operator=
551 MatrixSymmetric<Treal, Tmatrix> >
const & mm);
553 MatrixSymmetric<Treal, Tmatrix>&
operator+=
554 (MatrixSymmetric<Treal, Tmatrix>
const &
A);
556 MatrixSymmetric<Treal, Tmatrix>&
operator-=
557 (MatrixSymmetric<Treal, Tmatrix>
const &
A);
560 MatrixSymmetric<Treal, Tmatrix>&
operator+=
564 MatrixSymmetric<Treal, Tmatrix>&
operator-=
567 template<
typename Top>
569 return this->
matrixPtr->syAccumulateWith(op);
577 this->
matrixPtr->syRandomZeroStructure(probabilityBeingZero);
584 template<
typename TRule>
586 this->
matrixPtr->sySetElementsByRule(rule);
592 void transfer( MatrixSymmetric<Treal, Tmatrix> & dest ) {
598 template<
typename Tvector>
600 Treal
const ONE = 1.0;
601 y = (ONE * (*this) * x);
620 (std::vector<int>
const & rowind,
621 std::vector<int>
const & rowPermutation,
622 std::vector<int> & newRowind,
623 std::vector<int>
const & colind,
624 std::vector<int>
const & colPermutation,
625 std::vector<int> & newColind) {
631 for (
unsigned int i = 0; i < newRowind.size(); ++i) {
632 if (newRowind[i] > newColind[i]) {
634 newRowind[i] = newColind[i];
642 template<
typename Treal,
typename Tmatrix>
653 frobNormMat.resetSizesAndBlocks(rows_new, cols_new);
654 frobNormMat.getMatrix().syAssignFrobNormsLowestLevel(this->
getMatrix());
655 return frobNormMat.eucl(requestedAccuracy, maxIter);
659 template<
typename Treal,
typename Tmatrix>
661 eucl(Treal
const requestedAccuracy,
663 assert(requestedAccuracy >= 0);
665 Treal frobTmp = this->
frob();
666 if (frobTmp < requestedAccuracy)
677 MatrixSymmetric<Treal, Tmatrix> Copy(*
this);
678 Copy.frob_thresh(requestedAccuracy / 2.0);
680 <Treal, MatrixSymmetric<Treal, Tmatrix>, VectorType>
681 lan(Copy, guess, maxIter);
682 lan.setAbsTol( requestedAccuracy / 2.0 );
685 <Treal, MatrixSymmetric<Treal, Tmatrix>, VectorType>
686 lan(*
this, guess, maxIter);
687 lan.setAbsTol( requestedAccuracy );
696 template<
typename Treal,
typename Tmatrix>
698 diff(
const MatrixSymmetric<Treal, Tmatrix>& A,
699 const MatrixSymmetric<Treal, Tmatrix>& B,
700 normType const norm, Treal
const requestedAccuracy) {
709 diff =
eucl_diff(A, B, requestedAccuracy);
710 eNMin = diff - requestedAccuracy;
711 eNMin = eNMin >= 0 ? eNMin : 0;
715 throw Failure(
"MatrixSymmetric<Treal, Tmatrix>::" 716 "diff(const MatrixSymmetric<Treal, Tmatrix>&, " 717 "const MatrixSymmetric<Treal, Tmatrix>&, " 718 "normType const, Treal): " 719 "Diff not implemented for selected norm");
724 template<
typename Treal,
typename Tmatrix>
727 const MatrixSymmetric<Treal, Tmatrix>& B,
729 Treal
const requestedAccuracy,
730 Treal
const maxAbsVal) {
743 throw Failure(
"MatrixSymmetric<Treal, Tmatrix>::" 745 "(const MatrixSymmetric<Treal, Tmatrix>&, " 746 "const MatrixSymmetric<Treal, Tmatrix>&, " 747 "normType const, Treal const, Treal const): " 748 "Diff not implemented for selected norm");
755 template<
typename Treal,
typename Tmatrix>
757 (
const MatrixSymmetric<Treal, Tmatrix>& A,
758 const MatrixSymmetric<Treal, Tmatrix>& B,
759 Treal
const requestedAccuracy) {
765 Treal relTol = sqrt(sqrt(std::numeric_limits<Treal>::epsilon()));
771 template<
typename Treal,
typename Tmatrix>
773 (
const MatrixSymmetric<Treal, Tmatrix>& A,
774 const MatrixSymmetric<Treal, Tmatrix>& B,
775 Treal
const requestedAccuracy) {
780 A.getSizesAndBlocksForFrobNormMat( rows_new, cols_new );
781 frobNormMat.resetSizesAndBlocks(rows_new, cols_new);
782 frobNormMat.getMatrix().syAssignDiffFrobNormsLowestLevel(A.getMatrix(),B.getMatrix());
784 return frobNormMat.eucl(requestedAccuracy);
788 template<
typename Treal,
typename Tmatrix>
790 (
const MatrixSymmetric<Treal, Tmatrix>& A,
791 const MatrixSymmetric<Treal, Tmatrix>& B,
792 Treal
const requestedAccuracy,
793 Treal
const maxAbsVal,
794 VectorType *
const eVecPtr) {
798 Treal relTol = sqrt(sqrt(std::numeric_limits<Treal>::epsilon()));
809 if ( tmpInterval.
length() < 2*requestedAccuracy )
817 template<
typename Treal,
typename Tmatrix>
832 throw Failure(
"MatrixSymmetric<Treal, Tmatrix>::" 833 "thresh(Treal const, " 835 "Thresholding not implemented for selected norm");
841 template<
typename Treal,
typename Tmatrix>
845 if ( Zptr == NULL ) {
847 return TruncObj.
run( threshold );
850 return TruncObj.
run( threshold );
856 template<
typename Treal,
typename Tmatrix>
859 std::vector<int> rows_block_sizes;
860 std::vector<int> cols_block_sizes;
870 rows_block_sizes.pop_back();
871 cols_block_sizes.pop_back();
874 int factor_rows = rows_block_sizes[rows_block_sizes.size()-1];
875 int factor_cols = cols_block_sizes[cols_block_sizes.size()-1];
876 for (
unsigned int ind = 0; ind < rows_block_sizes.size(); ++ind)
877 rows_block_sizes[ind] = rows_block_sizes[ind] / factor_rows;
878 for (
unsigned int ind = 0; ind < cols_block_sizes.size(); ++ind)
879 cols_block_sizes[ind] = cols_block_sizes[ind] / factor_cols;
883 if (n_rows % factor_rows)
884 n_rows = n_rows / factor_rows + 1;
886 n_rows = n_rows / factor_rows;
887 if (n_cols % factor_cols)
888 n_cols = n_cols / factor_cols + 1;
890 n_cols = n_cols / factor_cols;
897 template<
typename Treal,
typename Tmatrix>
900 assert(threshold >= (Treal)0.0);
901 if (threshold == (Treal)0.0)
910 frobNormMat.resetSizesAndBlocks(rows_new, cols_new);
914 frobNormMat.getMatrix().syAssignFrobNormsLowestLevel(this->
getMatrix());
916 Treal mixed_norm_result = TruncObj.
run( threshold );
917 frobNormMat.getMatrix().truncateAccordingToSparsityPattern(this->
getMatrix());
918 return mixed_norm_result;
923 template<
typename Treal,
typename Tmatrix>
924 inline MatrixSymmetric<Treal, Tmatrix>&
930 this->
matrixPtr->assign(sm.A, *sm.B.matrixPtr);
934 template<
typename Treal,
typename Tmatrix>
935 MatrixSymmetric<Treal, Tmatrix>&
936 MatrixSymmetric<Treal, Tmatrix>::operator=
938 MatrixSymmetric<Treal, Tmatrix>,
939 MatrixSymmetric<Treal, Tmatrix>,
941 MatrixSymmetric<Treal, Tmatrix> >& sm2psm) {
942 assert(
this != &sm2psm.B);
943 if (
this == &sm2psm.E && &sm2psm.B == &sm2psm.C) {
946 sm2psm.A, *sm2psm.B.matrixPtr,
947 sm2psm.D, *this->matrixPtr);
951 throw Failure(
"MatrixSymmetric<Treal, Tmatrix>::operator=" 952 "(const XYZpUV<Treal, MatrixSymmetric" 953 "<Treal, Tmatrix> >& sm2psm) : " 954 "D = alpha * A * B + beta * C not supported for C != D" 955 " and for symmetric matrices not for A != B since this " 956 "generally will result in a nonsymmetric matrix");
960 template<
typename Treal,
typename Tmatrix>
961 MatrixSymmetric<Treal, Tmatrix>&
962 MatrixSymmetric<Treal, Tmatrix>::operator=
964 MatrixSymmetric<Treal, Tmatrix>,
965 MatrixSymmetric<Treal, Tmatrix> >& sm2) {
966 assert(
this != &sm2.B);
967 if (&sm2.B == &sm2.C) {
969 Tmatrix::sysq(
'U', sm2.A, *sm2.B.matrixPtr, 0, *this->matrixPtr);
973 throw Failure(
"MatrixSymmetric<Treal, Tmatrix>::operator=" 974 "(const XYZ<Treal, MatrixSymmetric<Treal, Tmatrix>," 975 " MatrixSymmetric<Treal, Tmatrix> >& sm2) : " 976 "Operation C = alpha * A * B with only symmetric " 977 "matrices not supported for A != B");
981 template<
typename Treal,
typename Tmatrix>
982 MatrixSymmetric<Treal, Tmatrix>&
983 MatrixSymmetric<Treal, Tmatrix>::operator+=
985 MatrixSymmetric<Treal, Tmatrix>,
986 MatrixSymmetric<Treal, Tmatrix> >& sm2) {
987 assert(
this != &sm2.B);
988 if (&sm2.B == &sm2.C) {
989 Tmatrix::sysq(
'U', sm2.A, *sm2.B.matrixPtr, 1, *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>
1002 MatrixSymmetric<Treal, Tmatrix>&
1003 MatrixSymmetric<Treal, Tmatrix>::operator=
1005 MatrixGeneral<Treal, Tmatrix>,
1006 MatrixGeneral<Treal, Tmatrix>,
1008 MatrixSymmetric<Treal, Tmatrix> >& smmpsm) {
1009 if (
this == &smmpsm.E)
1010 if (&smmpsm.B == &smmpsm.C)
1011 if (!smmpsm.tB && smmpsm.tC) {
1013 smmpsm.A, *smmpsm.B.matrixPtr,
1014 smmpsm.D, *this->matrixPtr);
1017 throw Failure(
"MatrixSymmetric<Treal, Tmatrix>::operator=" 1018 "(const XYZpUV<Treal, MatrixGeneral" 1019 "<Treal, Tmatrix>, " 1020 "MatrixGeneral<Treal, Tmatrix>, Treal, " 1021 "MatrixSymmetric<Treal, Tmatrix> >&) : " 1022 "C = alpha * A' * A + beta * C, not implemented" 1023 " only C = alpha * A * A' + beta * C");
1025 throw Failure(
"MatrixSymmetric<Treal, Tmatrix>::operator=" 1027 "Treal, MatrixGeneral<Treal, Tmatrix>, " 1028 "MatrixGeneral<Treal, Tmatrix>, Treal, " 1029 "MatrixSymmetric<Treal, Tmatrix> >&) : " 1030 "You are trying to call C = alpha * A * A' + beta * C" 1031 " with A and A' being different objects");
1033 throw Failure(
"MatrixSymmetric<Treal, Tmatrix>::operator=" 1035 "Treal, MatrixGeneral<Treal, Tmatrix>, " 1036 "MatrixGeneral<Treal, Tmatrix>, Treal, " 1037 "MatrixSymmetric<Treal, Tmatrix> >&) : " 1038 "D = alpha * A * A' + beta * C not supported for C != D");
1043 template<
typename Treal,
typename Tmatrix>
1044 MatrixSymmetric<Treal, Tmatrix>&
1045 MatrixSymmetric<Treal, Tmatrix>::operator=
1047 MatrixGeneral<Treal, Tmatrix>,
1048 MatrixGeneral<Treal, Tmatrix> >& smm) {
1049 if (&smm.B == &smm.C)
1050 if (!smm.tB && smm.tC) {
1052 smm.A, *smm.B.matrixPtr,
1053 0, *this->matrixPtr);
1056 throw Failure(
"MatrixSymmetric<Treal, Tmatrix>::operator=" 1058 "Treal, MatrixGeneral<Treal, Tmatrix>, " 1059 "MatrixGeneral<Treal, Tmatrix> >&) : " 1060 "C = alpha * A' * A, not implemented " 1061 "only C = alpha * A * A'");
1063 throw Failure(
"MatrixSymmetric<Treal, Tmatrix>::operator=" 1065 "Treal, MatrixGeneral<Treal, Tmatrix>, " 1066 "MatrixGeneral<Treal, Tmatrix> >&) : " 1067 "You are trying to call C = alpha * A * A' " 1068 "with A and A' being different objects");
1073 template<
typename Treal,
typename Tmatrix>
1074 MatrixSymmetric<Treal, Tmatrix>&
1075 MatrixSymmetric<Treal, Tmatrix>::operator+=
1077 MatrixGeneral<Treal, Tmatrix>,
1078 MatrixGeneral<Treal, Tmatrix> >& smm) {
1079 if (&smm.B == &smm.C)
1080 if (!smm.tB && smm.tC) {
1082 smm.A, *smm.B.matrixPtr,
1083 1, *this->matrixPtr);
1086 throw Failure(
"MatrixSymmetric<Treal, Tmatrix>::operator+=" 1088 "Treal, MatrixGeneral<Treal, Tmatrix>, " 1089 "MatrixGeneral<Treal, Tmatrix> >&) : " 1090 "C += alpha * A' * A, not implemented " 1091 "only C += alpha * A * A'");
1093 throw Failure(
"MatrixSymmetric<Treal, Tmatrix>::operator+=" 1095 "Treal, MatrixGeneral<Treal, Tmatrix>, " 1096 "MatrixGeneral<Treal, Tmatrix> >&) : " 1097 "You are trying to call C += alpha * A * A' " 1098 "with A and A' being different objects");
1105 template<
typename Treal,
typename Tmatrix>
1106 MatrixSymmetric<Treal, Tmatrix>&
1107 MatrixSymmetric<Treal, Tmatrix>::operator=
1109 MatrixSymmetric<Treal, Tmatrix>,
1111 if (
this == &zaz.B) {
1112 if (&zaz.A == &zaz.C) {
1113 if (zaz.tA && !zaz.tC) {
1114 Tmatrix::trsytriplemm(
'R', *zaz.A.matrixPtr, *this->matrixPtr);
1116 else if (!zaz.tA && zaz.tC) {
1117 Tmatrix::trsytriplemm(
'L', *zaz.A.matrixPtr, *this->matrixPtr);
1120 throw Failure(
"MatrixSymmetric<Treal, Tmatrix>::operator=" 1121 "(const XYZ<MatrixTriangular<Treal, Tmatrix>," 1122 "MatrixSymmetric<Treal, Tmatrix>," 1123 "MatrixTriangular<Treal, Tmatrix> >&) : " 1124 "A = op1(Z) * A * op2(Z) : Either op1 xor op2 must be " 1125 "the transpose operation.");
1128 throw Failure(
"MatrixSymmetric<Treal, Tmatrix>::operator=" 1129 "(const XYZ<MatrixTriangular<Treal, Tmatrix>," 1130 "MatrixSymmetric<Treal, Tmatrix>," 1131 "MatrixTriangular<Treal, Tmatrix> >&) : " 1132 "A = op1(Z1) * A * op2(Z2) : Z1 and Z2 must be the same " 1136 throw Failure(
"MatrixSymmetric<Treal, Tmatrix>::operator=" 1137 "(const XYZ<MatrixTriangular<Treal, Tmatrix>," 1138 "MatrixSymmetric<Treal, Tmatrix>," 1139 "MatrixTriangular<Treal, Tmatrix> >&) : " 1140 "C = op1(Z) * A * op2(Z) : A and C must be the same " 1152 template<
typename Treal,
typename Tmatrix>
1155 const MatrixSymmetric<Treal, Tmatrix>& A,
1156 const MatrixSymmetric<Treal, Tmatrix>& B,
1158 MatrixSymmetric<Treal, Tmatrix>& C) {
1159 Tmatrix::ssmm_upper_tr_only(alpha, *A.matrixPtr, *B.matrixPtr,
1160 beta, *C.matrixPtr);
1167 template<
typename Treal,
typename Tmatrix>
1168 inline MatrixSymmetric<Treal, Tmatrix>&
1169 MatrixSymmetric<Treal, Tmatrix>::operator=
1171 MatrixSymmetric<Treal, Tmatrix> >& mpm) {
1172 assert(
this != &mpm.A);
1174 Tmatrix::add(1.0, *mpm.A.matrixPtr, *this->matrixPtr);
1178 template<
typename Treal,
typename Tmatrix>
1179 inline MatrixSymmetric<Treal, Tmatrix>&
1180 MatrixSymmetric<Treal, Tmatrix>::operator=
1182 MatrixSymmetric<Treal, Tmatrix> >& mmm) {
1183 assert(
this != &mmm.B);
1185 Tmatrix::add(-1.0, *mmm.B.matrixPtr, *this->matrixPtr);
1189 template<
typename Treal,
typename Tmatrix>
1190 inline MatrixSymmetric<Treal, Tmatrix>&
1191 MatrixSymmetric<Treal, Tmatrix>::operator+=
1192 (MatrixSymmetric<Treal, Tmatrix>
const &
A) {
1193 Tmatrix::add(1.0, *A.matrixPtr, *this->matrixPtr);
1197 template<
typename Treal,
typename Tmatrix>
1198 inline MatrixSymmetric<Treal, Tmatrix>&
1199 MatrixSymmetric<Treal, Tmatrix>::operator-=
1200 (MatrixSymmetric<Treal, Tmatrix>
const &
A) {
1201 Tmatrix::add(-1.0, *A.matrixPtr, *this->matrixPtr);
1205 template<
typename Treal,
typename Tmatrix>
1206 inline MatrixSymmetric<Treal, Tmatrix>&
1207 MatrixSymmetric<Treal, Tmatrix>::operator+=
1210 Tmatrix::add(sm.A, *sm.B.matrixPtr, *this->matrixPtr);
1215 template<
typename Treal,
typename Tmatrix>
1216 inline MatrixSymmetric<Treal, Tmatrix>&
1217 MatrixSymmetric<Treal, Tmatrix>::operator-=
1220 Tmatrix::add(-sm.A, *sm.B.matrixPtr, *this->matrixPtr);
1228 template<
typename Treal,
typename Tmatrix,
typename Top>
1231 return A.accumulateWith(op);
void random()
Definition: MatrixSymmetric.h:572
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:240
Treal mixed_norm(Treal const requestedAccuracy, int maxIter=-1) const
Definition: MatrixSymmetric.h:644
Treal mixed_norm_thresh(Treal const threshold)
Definition: MatrixSymmetric.h:899
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:265
void writeToFileProt(std::ofstream &file) const
Write object to file.
Definition: MatrixSymmetric.h:607
Treal frob_thresh(Treal const threshold)
Does thresholding so that the error in the Frobenius norm is below the given threshold.
Definition: MatrixSymmetric.h:446
Treal eucl(Treal const requestedAccuracy, int maxIter=-1) const
Definition: MatrixSymmetric.h:661
Normal matrix.
Definition: MatrixBase.h:47
Treal length() const
Returns the length of the interval.
Definition: Interval.h:107
MatrixSymmetric< Treal, Tmatrix > & operator=(const MatrixSymmetric< Treal, Tmatrix > &symm)
Definition: MatrixSymmetric.h:334
void getBlockSizeVector(std::vector< int > &blockSizesCopy) const
Definition: SizesAndBlocks.cc:79
MatrixBase< Treal, Tmatrix > & operator=(const MatrixBase< Treal, Tmatrix > &other)
Definition: MatrixBase.h:164
void quickEuclBounds(Treal &euclLowerBound, Treal &euclUpperBound) const
Definition: MatrixSymmetric.h:357
std::string obj_type_id() const
Definition: MatrixSymmetric.h:605
ValidPtr< Tmatrix > matrixPtr
Definition: MatrixBase.h:151
size_t nnz() const
Definition: MatrixSymmetric.h:472
int getNTotalScalars() const
Definition: SizesAndBlocks.h:76
MatrixSymmetric< Treal, Tmatrix > & operator=(const MatrixGeneral< Treal, Tmatrix > &matr)
Definition: MatrixSymmetric.h:339
Definition: LanczosLargestMagnitudeEig.h:44
void read_from_buffer_base(void *buffer, const int n_bytes, const matrix_type mattype)
Definition: MatrixBase.h:279
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:340
Truncation of symmetric matrices at the element level (used for mixed norm truncation) ...
Definition: truncation.h:244
MatrixSymmetric< Treal, Tmatrix > & operator=(int const k)
Definition: MatrixSymmetric.h:344
Treal accumulateWith(Top &op)
Definition: MatrixSymmetric.h:568
MatrixSymmetric(const MatrixSymmetric< Treal, Tmatrix > &symm)
Copy constructor.
Definition: MatrixSymmetric.h:73
This proxy expresses the result of addition of two objects, of possibly different types...
Definition: matrix_proxy.h:244
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:790
MatrixSymmetric(const MatrixGeneral< Treal, Tmatrix > &matr)
'Copy from normal matrix' - constructor
Definition: MatrixSymmetric.h:77
Definition: MatrixBase.h:53
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:332
void clear()
Release memory for the information written to file.
Definition: MatrixBase.h:116
Definition: allocate.cc:30
Describes dimensions of matrix and its blocks on all levels.
Definition: SizesAndBlocks.h:37
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:726
void fullMatrix(std::vector< Treal > &fullMat) const
Save matrix as full matrix.
Definition: MatrixSymmetric.h:130
Treal midPoint() const
Definition: Interval.h:113
Treal template_blas_fabs(Treal x)
Definition: Interval.h:44
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:773
void assignFromFull(std::vector< Treal > const &fullMat)
Definition: MatrixSymmetric.h:102
void readFromFileBase(std::ifstream &file, matrix_type const mattype)
Definition: MatrixBase.h:234
void read_from_buffer(void *buffer, const int n_bytes)
Definition: MatrixSymmetric.h:481
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:1154
static void swap(ValidPtr< Tobj > &ptrA, ValidPtr< Tobj > &ptrB)
Definition: ValidPtr.h:104
MatrixSymmetric(const XY< Treal, MatrixSymmetric< Treal, Tmatrix > > &sm)
Definition: MatrixSymmetric.h:75
Definition: matInclude.h:135
void getCols(SizesAndBlocks &colsCopy) const
Definition: MatrixBase.h:83
void matVecProd(Tvector &y, Tvector const &x) const
Definition: MatrixSymmetric.h:599
int get_nrows() const
Definition: MatrixBase.h:136
This proxy expresses the result of multiplication of three objects, of possibly different types...
Definition: matrix_proxy.h:65
Upper non-unit triangular matrix.
Definition: MatrixBase.h:51
MatrixSymmetric()
Default constructor.
Definition: MatrixSymmetric.h:71
void randomZeroStructure(Treal probabilityBeingZero)
Definition: MatrixSymmetric.h:576
Treal run(Treal const threshold)
Definition: truncation.h:78
void gersgorin(Treal &lmin, Treal &lmax) const
Definition: MatrixSymmetric.h:464
Interval< Treal > euclIfSmall(Tmatrix const &M, Treal const requestedAbsAccuracy, Treal const requestedRelAccuracy, Treal const maxAbsVal, typename Tmatrix::VectorType *const eVecPtr=0)
Returns interval containing the Euclidean norm of the matrix M.
Definition: LanczosLargestMagnitudeEig.h:257
size_t nvalues() const
Definition: MatrixSymmetric.h:475
int get_ncols() const
Definition: MatrixBase.h:139
static Treal eucl_diff(const MatrixSymmetric< Treal, Tmatrix > &A, const MatrixSymmetric< Treal, Tmatrix > &B, Treal const requestedAccuracy)
Returns the Euclidean norm of A - B ( || A - B ||_2 )
Definition: MatrixSymmetric.h:757
Definition: matInclude.h:135
void rand()
Definition: VectorGeneral.h:94
void resetSizesAndBlocks(SizesAndBlocks const &newRows, SizesAndBlocks const &newCols)
Definition: MatrixBase.h:74
void resetSizesAndBlocks(SizesAndBlocks const &newRows)
Definition: VectorGeneral.h:49
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:168
Base class for matrix API.
Definition: MatrixBase.h:67
void readFromFileProt(std::ifstream &file)
Read object from file.
Definition: MatrixSymmetric.h:610
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:392
Definition: mat_utils.h:34
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:295
void transfer(MatrixSymmetric< Treal, Tmatrix > &dest)
Transfer this matrix to dest, clearing previous content of dest if any.
Definition: MatrixSymmetric.h:592
void getLargestMagnitudeEig(Treal &ev, Treal &accuracy)
Definition: LanczosLargestMagnitudeEig.h:56
This proxy expresses the result of substraction of two objects, of possibly different types...
Definition: matrix_proxy.h:264
Tmatrix const & getMatrix() const
Definition: MatrixBase.h:144
static Treal trace_ab(const MatrixSymmetric< Treal, Tmatrix > &A, const MatrixSymmetric< Treal, Tmatrix > &B)
Definition: MatrixSymmetric.h:468
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:698
Base class for matrix API.
void getSizesAndBlocksForFrobNormMat(SizesAndBlocks &rows_new, SizesAndBlocks &cols_new) const
Definition: MatrixSymmetric.h:858
Definition: matInclude.h:135
Truncation of symmetric matrices.
Definition: truncation.h:152
Treal frob() const
Definition: MatrixSymmetric.h:349
void writeToFileBase(std::ofstream &file, matrix_type const mattype) const
Definition: MatrixBase.h:219
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:620
Truncation of symmetric matrices with Z.
Definition: truncation.h:201
This proxy expresses the result of multiplication of two objects, of possibly different types...
Definition: matrix_proxy.h:49
void write_to_buffer_base(void *buffer, const int n_bytes, const matrix_type mattype) const
Definition: MatrixBase.h:252
void simple_blockwise_frob_thresh(Treal const threshold)
Definition: MatrixSymmetric.h:460
void haveDataStructureSet(bool val)
Definition: ValidPtr.h:97
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:1229
void getRows(SizesAndBlocks &rowsCopy) const
Definition: MatrixBase.h:80
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:585
VectorGeneral< Treal, typename Tmatrix::VectorType > VectorType
Definition: MatrixSymmetric.h:68
Treal eucl_thresh(Treal const threshold, MatrixTriangular< Treal, Tmatrix > const *const Zptr=NULL)
Definition: MatrixSymmetric.h:843
Symmetric matrix.
Definition: MatrixBase.h:49
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:86
Definition: MatrixBase.h:54
Treal template_blas_sqrt(Treal x)
void write_to_buffer(void *buffer, const int n_bytes) const
Definition: MatrixSymmetric.h:478
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:819
normType
Definition: matInclude.h:135
Classes for truncation of small matrix elements.
Treal real
Definition: MatrixSymmetric.h:69
static void getPermutedIndexes(std::vector< int > const &index, std::vector< int > const &permutation, std::vector< int > &newIndex)
Definition: MatrixBase.h:203