36 #ifndef MAT_MATRIXGENERAL
37 #define MAT_MATRIXGENERAL
56 template<
typename Treal,
typename Tmatrix>
57 class MatrixGeneral :
public MatrixBase<Treal, Tmatrix> {
74 template<
typename Tfull>
75 inline void assign_from_full
76 (Tfull
const*
const fullmatrix,
77 const int totnrows,
const int totncols) {
78 this->
matrixPtr->assign_from_full(fullmatrix, totnrows, totncols);
80 inline void assign_from_full
81 (Treal
const*
const fullmatrix,
82 const int totnrows,
const int totncols) {
83 this->
matrixPtr->assign_from_full(fullmatrix, totnrows, totncols);
88 (std::vector<Treal>
const & fullMat) {
93 inline void fullMatrix(std::vector<Treal> & fullMat)
const {
98 (std::vector<Treal> & fullMat,
99 std::vector<int>
const & rowInversePermutation,
100 std::vector<int>
const & colInversePermutation)
const {
101 std::vector<int> rowind;
102 std::vector<int> colind;
103 std::vector<Treal> values;
105 rowInversePermutation,
106 colInversePermutation);
108 for (
unsigned int ind = 0; ind < values.size(); ++ind)
109 fullMat[rowind[ind] + this->
get_nrows() * colind[ind]] =
120 (std::vector<int>
const & rowind,
121 std::vector<int>
const & colind,
122 std::vector<Treal>
const & values,
126 this->
matrixPtr->assignFromSparse(rowind, colind, values);
136 (std::vector<int>
const & rowind,
137 std::vector<int>
const & colind,
138 std::vector<Treal>
const & values,
139 std::vector<int>
const & rowPermutation,
140 std::vector<int>
const & colPermutation) {
141 std::vector<int> newRowind;
142 std::vector<int> newColind;
147 this->
matrixPtr->assignFromSparse(newRowind, newColind, values);
155 (std::vector<int>
const & rowind,
156 std::vector<int>
const & colind,
157 std::vector<Treal>
const & values,
160 std::vector<int>
const & rowPermutation,
161 std::vector<int>
const & colPermutation) {
164 rowPermutation, colPermutation);
171 (std::vector<int>
const & rowind,
172 std::vector<int>
const & colind,
173 std::vector<Treal> & values)
const {
174 this->
matrixPtr->getValues(rowind, colind, values);
184 (std::vector<int>
const & rowind,
185 std::vector<int>
const & colind,
186 std::vector<Treal> & values,
187 std::vector<int>
const & rowPermutation,
188 std::vector<int>
const & colPermutation)
const {
189 std::vector<int> newRowind;
190 std::vector<int> newColind;
195 this->
matrixPtr->getValues(newRowind, newColind, values);
202 (std::vector<int> & rowind,
203 std::vector<int> & colind,
204 std::vector<Treal> & values)
const {
208 this->
matrixPtr->getAllValues(rowind, colind, values);
220 (std::vector<int> & rowind,
221 std::vector<int> & colind,
222 std::vector<Treal> & values,
223 std::vector<int>
const & rowInversePermutation,
224 std::vector<int>
const & colInversePermutation)
const {
225 std::vector<int> tmpRowind;
226 std::vector<int> tmpColind;
227 tmpRowind.reserve(rowind.capacity());
228 tmpColind.reserve(colind.capacity());
230 this->
matrixPtr->getAllValues(tmpRowind, tmpColind, values);
247 inline void fullmatrix(Treal*
const full,
248 const int totnrows,
const int totncols)
const {
249 this->
matrixPtr->fullmatrix(full, totnrows, totncols);
288 return Tmatrix::frobDiff(*A.matrixPtr, *B.matrixPtr);
290 Treal
eucl(Treal
const requestedAccuracy,
291 int maxIter = -1)
const;
294 void thresh(Treal
const threshold,
308 return Tmatrix::trace_ab(*A.matrixPtr, *B.matrixPtr);
313 return Tmatrix::trace_aTb(*A.matrixPtr, *B.matrixPtr);
315 inline size_t nnz()
const {
451 this->
matrixPtr->randomZeroStructure(probabilityBeingZero);
454 template<
typename TRule>
456 this->
matrixPtr->setElementsByRule(rule);
472 template<
typename Treal,
typename Tmatrix>
474 eucl(Treal
const requestedAccuracy,
483 maxIter = this->get_nrows() * 100;
486 lan(ata, guess, maxIter);
487 lan.setRelTol( 100 * std::numeric_limits<Treal>::epsilon() );
494 if ( euclInt.
low() < 0 )
496 if ( euclInt.
length() / 2.0 > requestedAccuracy ) {
497 std::cout <<
"req: " << requestedAccuracy
498 <<
" obt: " << euclInt.
length() / 2.0 << std::endl;
499 throw std::runtime_error(
"Desired accuracy not obtained in Lanczos.");
505 template<
typename Treal,
typename Tmatrix>
511 this->frob_thresh(threshold);
514 throw Failure(
"MatrixGeneral<Treal, Tmatrix>::"
515 "thresh(Treal const, "
517 "Thresholding not imlpemented for selected norm");
521 template<
typename Treal,
typename Tmatrix>
525 return TruncObj.
run( threshold );
531 template<
typename Treal,
typename Tmatrix>
537 assert(
this != &smm.B &&
this != &smm.C);
538 this->matrixPtr.haveDataStructureSet(
true);
540 *smm.B.matrixPtr, *smm.C.matrixPtr,
541 0, *this->matrixPtr);
546 template<
typename Treal,
typename Tmatrix>
551 assert(
this != &mm.A &&
this != &mm.B);
553 *mm.A.matrixPtr, *mm.B.matrixPtr,
554 0, *this->matrixPtr);
559 template<
typename Treal,
typename Tmatrix>
565 assert(
this != &smm.B &&
this != &smm.C);
567 *smm.B.matrixPtr, *smm.C.matrixPtr,
568 1, *this->matrixPtr);
573 template<
typename Treal,
typename Tmatrix>
581 assert(
this != &smmpsm.B &&
this != &smmpsm.C);
583 if (
this == &smmpsm.E)
585 *smmpsm.B.matrixPtr, *smmpsm.C.matrixPtr,
586 smmpsm.D, *this->matrixPtr);
588 throw Failure(
"MatrixGeneral<Treal, Tmatrix>::operator="
589 "(const XYZpUV<Treal, MatrixGeneral"
590 "<Treal, Tmatrix> >&) : D = alpha "
591 "* op(A) * op(B) + beta * C not supported for C != D");
597 template<
typename Treal,
typename Tmatrix>
602 assert(
this != &mpm.A);
604 Tmatrix::add(1.0, *mpm.A.matrixPtr, *this->matrixPtr);
608 template<
typename Treal,
typename Tmatrix>
612 Tmatrix::add(1.0, *A.matrixPtr, *this->matrixPtr);
616 template<
typename Treal,
typename Tmatrix>
620 Tmatrix::add(-1.0, *A.matrixPtr, *this->matrixPtr);
624 template<
typename Treal,
typename Tmatrix>
629 Tmatrix::add(sm.A, *sm.B.matrixPtr, *this->matrixPtr);
636 template<
typename Treal,
typename Tmatrix>
642 assert(
this != &smm.C);
643 assert(!smm.tB && !smm.tC);
644 this->matrixPtr.haveDataStructureSet(
true);
646 *smm.B.matrixPtr, *smm.C.matrixPtr,
647 0, *this->matrixPtr);
652 template<
typename Treal,
typename Tmatrix>
657 assert(
this != &mm.B);
660 *mm.A.matrixPtr, *mm.B.matrixPtr,
661 0, *this->matrixPtr);
666 template<
typename Treal,
typename Tmatrix>
672 assert(
this != &smm.C);
673 assert(!smm.tB && !smm.tC);
675 *smm.B.matrixPtr, *smm.C.matrixPtr,
676 1, *this->matrixPtr);
681 template<
typename Treal,
typename Tmatrix>
689 assert(
this != &smmpsm.C);
690 assert(!smmpsm.tB && !smmpsm.tC && !smmpsm.tE);
691 if (
this == &smmpsm.E)
693 *smmpsm.B.matrixPtr, *smmpsm.C.matrixPtr,
694 smmpsm.D, *this->matrixPtr);
696 throw Failure(
"MatrixGeneral<Treal, Tmatrix>::operator="
697 "(const XYZpUV<Treal, MatrixGeneral"
698 "<Treal, Tmatrix>, MatrixSymmetric<Treal, "
699 "Tmatrix>, Treal, MatrixGeneral"
700 "<Treal, Tmatrix> >&) "
701 ": D = alpha * A * B + beta * C (with A symmetric)"
702 " not supported for C != D");
707 template<
typename Treal,
typename Tmatrix>
713 assert(
this != &smm.B);
714 assert(!smm.tB && !smm.tC);
715 this->matrixPtr.haveDataStructureSet(
true);
717 *smm.C.matrixPtr, *smm.B.matrixPtr,
718 0, *this->matrixPtr);
723 template<
typename Treal,
typename Tmatrix>
728 assert(
this != &mm.A);
729 assert(!mm.tA && !mm.tB);
731 *mm.B.matrixPtr, *mm.A.matrixPtr,
732 0, *this->matrixPtr);
737 template<
typename Treal,
typename Tmatrix>
743 assert(
this != &smm.B);
744 assert(!smm.tB && !smm.tC);
746 *smm.C.matrixPtr, *smm.B.matrixPtr,
747 1, *this->matrixPtr);
752 template<
typename Treal,
typename Tmatrix>
760 assert(
this != &smmpsm.B);
761 assert(!smmpsm.tB && !smmpsm.tC && !smmpsm.tE);
762 if (
this == &smmpsm.E)
764 *smmpsm.C.matrixPtr, *smmpsm.B.matrixPtr,
765 smmpsm.D, *this->matrixPtr);
767 throw Failure(
"MatrixGeneral<Treal, Tmatrix>::operator="
768 "(const XYZpUV<Treal, MatrixSymmetric"
769 "<Treal, Tmatrix>, MatrixGeneral<Treal, "
770 "Tmatrix>, Treal, MatrixGeneral"
771 "<Treal, Tmatrix> >&) "
772 ": D = alpha * B * A + beta * C (with A symmetric)"
773 " not supported for C != D");
779 template<
typename Treal,
typename Tmatrix>
785 assert(!smm.tB && !smm.tC);
786 this->matrixPtr.haveDataStructureSet(
true);
796 template<
typename Treal,
typename Tmatrix>
811 template<
typename Treal,
typename Tmatrix>
817 assert(!smm.tB && !smm.tC);
828 template<
typename Treal,
typename Tmatrix>
836 assert(!smmpsm.tB && !smmpsm.tC && !smmpsm.tE);
837 if (
this == &smmpsm.E)
838 Tmatrix::ssmm(smmpsm.A,
844 throw Failure(
"MatrixGeneral<Treal, Tmatrix>::"
845 "operator=(const XYZpUV<"
846 "Treal, MatrixSymmetric<Treal, Tmatrix>, "
847 "MatrixSymmetric<Treal, Tmatrix>, Treal, "
848 "MatrixGeneral<Treal, Tmatrix> >&) "
849 ": D = alpha * A * B + beta * C (with A and B symmetric)"
850 " not supported for C != D");
859 template<
typename Treal,
typename Tmatrix>
868 *smm.B.matrixPtr, *this->matrixPtr);
870 throw Failure(
"MatrixGeneral<Treal, Tmatrix>::operator="
871 "(const XYZ<Treal, MatrixTriangular"
872 "<Treal, Tmatrix>, MatrixGeneral<Treal,"
873 " Tmatrix> >& : D = alpha * op(A) * B (with"
874 " A upper triangular) not supported for B != D");
880 template<
typename Treal,
typename Tmatrix>
889 *smm.C.matrixPtr, *this->matrixPtr);
891 throw Failure(
"MatrixGeneral<Treal, Tmatrix>::operator="
892 "(const XYZ<Treal, MatrixGeneral"
893 "<Treal, Tmatrix>, MatrixTriangular<Treal,"
894 " Tmatrix> >& : D = alpha * A * op(B) (with"
895 " B upper triangular) not supported for A != D");
901 template<
typename Treal,
typename Tmatrix>
905 if ((!smm.tB && !smm.tC) || (smm.tB && smm.tC))