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,
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 {
567 template<
typename Top>
569 return this->
matrixPtr->syAccumulateWith(op);
577 this->
matrixPtr->syRandomZeroStructure(probabilityBeingZero);
584 template<
typename TRule>
586 this->
matrixPtr->sySetElementsByRule(rule);
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>
649 this->getSizesAndBlocksForFrobNormMat( rows_new, cols_new );
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)
669 maxIter = this->get_nrows() * 100;
678 Copy.frob_thresh(requestedAccuracy / 2.0);
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>
700 normType const norm, Treal
const requestedAccuracy) {
705 diff = frob_diff(A, B);
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>
729 Treal
const requestedAccuracy,
730 Treal
const maxAbsVal) {
735 diff = frob_diff(A, B);
740 return euclDiffIfSmall(A, B, requestedAccuracy, 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>
759 Treal
const requestedAccuracy) {
762 Treal maxAbsVal = 2 * frob_diff(A,B);
765 Treal relTol = sqrt(sqrt(std::numeric_limits<Treal>::epsilon()));
771 template<
typename Treal,
typename Tmatrix>
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>
792 Treal
const requestedAccuracy,
793 Treal
const maxAbsVal,
798 Treal relTol = sqrt(sqrt(std::numeric_limits<Treal>::epsilon()));
809 if ( tmpInterval.
length() < 2*requestedAccuracy )
817 template<
typename Treal,
typename Tmatrix>
823 return this->frob_thresh(threshold);
826 return this->eucl_thresh(threshold);
829 return this->mixed_norm_thresh(threshold);
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)
906 this->getSizesAndBlocksForFrobNormMat( rows_new, cols_new );
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>
929 this->matrixPtr.haveDataStructureSet(
true);
930 this->matrixPtr->assign(sm.A, *sm.B.matrixPtr);
934 template<
typename Treal,
typename Tmatrix>
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>
966 assert(
this != &sm2.B);
967 if (&sm2.B == &sm2.C) {
968 this->matrixPtr.haveDataStructureSet(
true);
969 Tmatrix::sysq(
'U', sm2.A, *sm2.B.matrixPtr, 0, *this->matrixPtr);
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>
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>
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>
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>
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>
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>
1159 Tmatrix::ssmm_upper_tr_only(alpha, *A.matrixPtr, *B.matrixPtr,
1160 beta, *C.matrixPtr);
1167 template<
typename Treal,
typename Tmatrix>
1172 assert(
this != &mpm.A);
1174 Tmatrix::add(1.0, *mpm.A.matrixPtr, *this->matrixPtr);
1178 template<
typename Treal,
typename Tmatrix>
1183 assert(
this != &mmm.B);
1185 Tmatrix::add(-1.0, *mmm.B.matrixPtr, *this->matrixPtr);
1189 template<
typename Treal,
typename Tmatrix>
1193 Tmatrix::add(1.0, *A.matrixPtr, *this->matrixPtr);
1197 template<
typename Treal,
typename Tmatrix>
1201 Tmatrix::add(-1.0, *A.matrixPtr, *this->matrixPtr);
1205 template<
typename Treal,
typename Tmatrix>
1210 Tmatrix::add(sm.A, *sm.B.matrixPtr, *this->matrixPtr);
1215 template<
typename Treal,
typename Tmatrix>
1220 Tmatrix::add(-sm.A, *sm.B.matrixPtr, *this->matrixPtr);
1228 template<
typename Treal,
typename Tmatrix,
typename Top>
1231 return A.accumulateWith(op);