ergo
MatrixGeneral.h
Go to the documentation of this file.
1/* Ergo, version 3.8.2, a program for linear scaling electronic structure
2 * calculations.
3 * Copyright (C) 2023 Elias Rudberg, Emanuel H. Rubensson, Pawel Salek,
4 * and Anastasia Kruchinina.
5 *
6 * This program is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation, either version 3 of the License, or
9 * (at your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program. If not, see <http://www.gnu.org/licenses/>.
18 *
19 * Primary academic reference:
20 * Ergo: An open-source program for linear-scaling electronic structure
21 * calculations,
22 * Elias Rudberg, Emanuel H. Rubensson, Pawel Salek, and Anastasia
23 * Kruchinina,
24 * SoftwareX 7, 107 (2018),
25 * <http://dx.doi.org/10.1016/j.softx.2018.03.005>
26 *
27 * For further information about Ergo, see <http://www.ergoscf.org>.
28 */
29
38#ifndef MAT_MATRIXGENERAL
39#define MAT_MATRIXGENERAL
40#include "MatrixBase.h"
41namespace mat {
58 template<typename Treal, typename Tmatrix>
59 class MatrixGeneral : public MatrixBase<Treal, Tmatrix> {
60 public:
62
64 :MatrixBase<Treal, Tmatrix>() {}
66 :MatrixBase<Treal, Tmatrix>(matr) {}
68 :MatrixBase<Treal, Tmatrix>(symm) {
69 this->matrixPtr->symToNosym();
70 }
72 :MatrixBase<Treal, Tmatrix>(triang) {}
75#if 0
76 template<typename Tfull>
77 inline void assign_from_full
78 (Tfull const* const fullmatrix,
79 const int totnrows, const int totncols) {
80 this->matrixPtr->assign_from_full(fullmatrix, totnrows, totncols);
81 }
82 inline void assign_from_full
83 (Treal const* const fullmatrix,
84 const int totnrows, const int totncols) {
85 this->matrixPtr->assign_from_full(fullmatrix, totnrows, totncols);
86 }
87#endif
88
89 inline void assignFromFull
90 (std::vector<Treal> const & fullMat) {
91 assert((int)fullMat.size() == this->get_nrows() * this->get_ncols());
92 this->matrixPtr->assignFromFull(fullMat);
93 }
94
95 inline void fullMatrix(std::vector<Treal> & fullMat) const {
96 this->matrixPtr->fullMatrix(fullMat);
97 }
98
99 inline void fullMatrix
100 (std::vector<Treal> & fullMat,
101 std::vector<int> const & rowInversePermutation,
102 std::vector<int> const & colInversePermutation) const {
103 std::vector<int> rowind;
104 std::vector<int> colind;
105 std::vector<Treal> values;
106 get_all_values(rowind, colind, values,
107 rowInversePermutation,
108 colInversePermutation);
109 fullMat.assign(this->get_nrows()*this->get_ncols(),0);
110 for (unsigned int ind = 0; ind < values.size(); ++ind)
111 fullMat[rowind[ind] + this->get_nrows() * colind[ind]] =
112 values[ind];
113 }
122 (std::vector<int> const & rowind,
123 std::vector<int> const & colind,
124 std::vector<Treal> const & values,
125 SizesAndBlocks const & newRows,
126 SizesAndBlocks const & newCols) {
127 this->resetSizesAndBlocks(newRows, newCols);
128 this->matrixPtr->assignFromSparse(rowind, colind, values);
129 }
138 (std::vector<int> const & rowind,
139 std::vector<int> const & colind,
140 std::vector<Treal> const & values,
141 std::vector<int> const & rowPermutation,
142 std::vector<int> const & colPermutation) {
143 std::vector<int> newRowind;
144 std::vector<int> newColind;
146 getPermutedIndexes(rowind, rowPermutation, newRowind);
148 getPermutedIndexes(colind, colPermutation, newColind);
149 this->matrixPtr->assignFromSparse(newRowind, newColind, values);
150 }
157 (std::vector<int> const & rowind,
158 std::vector<int> const & colind,
159 std::vector<Treal> const & values,
160 SizesAndBlocks const & newRows,
161 SizesAndBlocks const & newCols,
162 std::vector<int> const & rowPermutation,
163 std::vector<int> const & colPermutation) {
164 this->resetSizesAndBlocks(newRows, newCols);
165 assign_from_sparse(rowind, colind, values,
166 rowPermutation, colPermutation);
167 }
172 inline void get_values
173 (std::vector<int> const & rowind,
174 std::vector<int> const & colind,
175 std::vector<Treal> & values) const {
176 this->matrixPtr->getValues(rowind, colind, values);
177 }
185 inline void get_values
186 (std::vector<int> const & rowind,
187 std::vector<int> const & colind,
188 std::vector<Treal> & values,
189 std::vector<int> const & rowPermutation,
190 std::vector<int> const & colPermutation) const {
191 std::vector<int> newRowind;
192 std::vector<int> newColind;
194 getPermutedIndexes(rowind, rowPermutation, newRowind);
196 getPermutedIndexes(colind, colPermutation, newColind);
197 this->matrixPtr->getValues(newRowind, newColind, values);
198 }
203 inline void get_all_values
204 (std::vector<int> & rowind,
205 std::vector<int> & colind,
206 std::vector<Treal> & values) const {
207 rowind.resize(0);
208 colind.resize(0);
209 values.resize(0);
210 this->matrixPtr->getAllValues(rowind, colind, values);
211 }
221 inline void get_all_values
222 (std::vector<int> & rowind,
223 std::vector<int> & colind,
224 std::vector<Treal> & values,
225 std::vector<int> const & rowInversePermutation,
226 std::vector<int> const & colInversePermutation) const {
227 std::vector<int> tmpRowind;
228 std::vector<int> tmpColind;
229 tmpRowind.reserve(rowind.capacity());
230 tmpColind.reserve(colind.capacity());
231 values.resize(0);
232 this->matrixPtr->getAllValues(tmpRowind, tmpColind, values);
234 getPermutedIndexes(tmpRowind, rowInversePermutation, rowind);
236 getPermutedIndexes(tmpColind, colInversePermutation, colind);
237
238 }
248#if 0
249 inline void fullmatrix(Treal* const full,
250 const int totnrows, const int totncols) const {
251 this->matrixPtr->fullmatrix(full, totnrows, totncols);
252 }
253#endif
257 return *this;
258 }
261 if (mt.tA)
263 else
265 return *this;
266 }
267
271 this->matrixPtr->symToNosym();
272 return *this;
273 }
277 return *this;
278 }
279
281 *this->matrixPtr = k;
282 return *this;
283 }
284 inline Treal frob() const {
285 return this->matrixPtr->frob();
286 }
287 static inline Treal frob_diff
290 return Tmatrix::frobDiff(*A.matrixPtr, *B.matrixPtr);
291 }
292 Treal eucl(Treal const requestedAccuracy,
293 int maxIter = -1) const;
294
295
296 void thresh(Treal const threshold,
297 normType const norm);
298
299 inline void frob_thresh(Treal threshold) {
300 this->matrixPtr->frob_thresh(threshold);
301 }
302 Treal eucl_thresh(Treal const threshold);
303
304 inline void gershgorin(Treal& lmin, Treal& lmax) {
305 this->matrixPtr->gershgorin(lmin, lmax);
306 }
307 static inline Treal trace_ab
310 return Tmatrix::trace_ab(*A.matrixPtr, *B.matrixPtr);
311 }
312 static inline Treal trace_aTb
315 return Tmatrix::trace_aTb(*A.matrixPtr, *B.matrixPtr);
316 }
317 inline size_t nnz() const { /* Note: size_t instead of int here to avoid integer overflow. */
318 return this->matrixPtr->nnz();
319 }
320 inline size_t nvalues() const { /* Note: size_t instead of int here to avoid integer overflow. */
321 return this->matrixPtr->nvalues();
322 }
323
324 inline void write_to_buffer(void* buffer, const int n_bytes) const {
325 this->write_to_buffer_base(buffer, n_bytes, matrix_matr);
326 }
327 inline void read_from_buffer(void* buffer, const int n_bytes) {
328 this->read_from_buffer_base(buffer, n_bytes, matrix_matr);
329 }
330
331 /* OPERATIONS ONLY INVOLVING ORDINARY MATRICES */
334 (const XYZ<Treal,
337
342
345 (const XYZ<Treal,
348
351 (const XYZpUV<Treal,
354 Treal,
356
360 MatrixGeneral<Treal, Tmatrix> > const & mpm);
369
370
371 /* OPERATIONS INVOLVING SYMMETRIC MATRICES */
374 (const XYZ<Treal,
383 (const XYZ<Treal,
388 (const XYZpUV<Treal,
391 Treal,
395 (const XYZ<Treal,
404 (const XYZ<Treal,
409 (const XYZpUV<Treal,
412 Treal,
416 (const XYZ<Treal,
425 (const XYZ<Treal,
430 (const XYZpUV<Treal,
433 Treal,
435
436 /* OPERATIONS INVOLVING UPPER TRIANGULAR MATRICES */
439 (const XYZ<Treal,
444 (const XYZ<Treal,
447
448 void random() {
449 this->matrixPtr->random();
450 }
451
452 void randomZeroStructure(Treal probabilityBeingZero) {
453 this->matrixPtr->randomZeroStructure(probabilityBeingZero);
454 }
455
456 template<typename TRule>
457 void setElementsByRule(TRule & rule) {
458 this->matrixPtr->setElementsByRule(rule);
459 return;
460 }
461
462 std::string obj_type_id() const {return "MatrixGeneral";}
463 protected:
464 inline void writeToFileProt(std::ofstream & file) const {
465 this->writeToFileBase(file, matrix_matr);
466 }
467 inline void readFromFileProt(std::ifstream & file) {
468 this->readFromFileBase(file, matrix_matr);
469 }
470 private:
471
472 };
473
474 template<typename Treal, typename Tmatrix>
476 eucl(Treal const requestedAccuracy,
477 int maxIter) const {
478 VectorType guess;
480 this->getCols(cols);
482 guess.rand();
484 if (maxIter < 0)
485 maxIter = this->get_nrows() * 100;
488 lan(ata, guess, maxIter);
489 lan.setRelTol( 100 * mat::getMachineEpsilon<Treal>() );
490 lan.run();
491 Treal eVal = 0;
492 Treal acc = 0;
493 lan.getLargestMagnitudeEig(eVal, acc);
494 Interval<Treal> euclInt( sqrt(eVal-acc),
495 sqrt(eVal+acc) );
496 if ( euclInt.low() < 0 )
497 euclInt = Interval<Treal>( 0, sqrt(eVal+acc) );
498 if ( euclInt.length() / 2.0 > requestedAccuracy ) {
499 std::cout << "req: " << requestedAccuracy
500 << " obt: " << euclInt.length() / 2.0 << std::endl;
501 throw std::runtime_error("Desired accuracy not obtained in Lanczos.");
502 }
503 return euclInt.midPoint();
504 }
505
506
507 template<typename Treal, typename Tmatrix>
509 thresh(Treal const threshold,
510 normType const norm) {
511 switch (norm) {
512 case frobNorm:
513 this->frob_thresh(threshold);
514 break;
515 default:
516 throw Failure("MatrixGeneral<Treal, Tmatrix>::"
517 "thresh(Treal const, "
518 "normType const): "
519 "Thresholding not imlpemented for selected norm");
520 }
521 }
522
523 template<typename Treal, typename Tmatrix>
525 eucl_thresh(Treal const threshold) {
527 return TruncObj.run( threshold );
528 }
529
530
531 /* OPERATIONS ONLY INVOLVING ORDINARY MATRICES */
532 /* C = alpha * op(A) * op(B) */
533 template<typename Treal, typename Tmatrix>
536 (const XYZ<Treal,
539 assert(this != &smm.B && this != &smm.C);
540 this->matrixPtr.haveDataStructureSet(true);
541 Tmatrix::gemm(smm.tB, smm.tC, smm.A,
542 *smm.B.matrixPtr, *smm.C.matrixPtr,
543 0, *this->matrixPtr);
544 return *this;
545 }
546
547 /* C = op(A) * op(B) */
548 template<typename Treal, typename Tmatrix>
553 assert(this != &mm.A && this != &mm.B);
554 Tmatrix::gemm(mm.tA, mm.tB, 1.0,
555 *mm.A.matrixPtr, *mm.B.matrixPtr,
556 0, *this->matrixPtr);
557 return *this;
558 }
559
560 /* C += alpha * op(A) * op(B) */
561 template<typename Treal, typename Tmatrix>
564 (const XYZ<Treal,
567 assert(this != &smm.B && this != &smm.C);
568 Tmatrix::gemm(smm.tB, smm.tC, smm.A,
569 *smm.B.matrixPtr, *smm.C.matrixPtr,
570 1, *this->matrixPtr);
571 return *this;
572 }
573
574 /* C = alpha * op(A) * op(B) + beta * C */
575 template<typename Treal, typename Tmatrix>
578 (const XYZpUV<Treal,
581 Treal,
583 assert(this != &smmpsm.B && this != &smmpsm.C);
584 assert(!smmpsm.tE);
585 if (this == &smmpsm.E)
586 Tmatrix::gemm(smmpsm.tB, smmpsm.tC, smmpsm.A,
587 *smmpsm.B.matrixPtr, *smmpsm.C.matrixPtr,
588 smmpsm.D, *this->matrixPtr);
589 else
590 throw Failure("MatrixGeneral<Treal, Tmatrix>::operator="
591 "(const XYZpUV<Treal, MatrixGeneral"
592 "<Treal, Tmatrix> >&) : D = alpha "
593 "* op(A) * op(B) + beta * C not supported for C != D");
594 return *this;
595 }
596
597
598 /* C = A + B */
599 template<typename Treal, typename Tmatrix>
604 assert(this != &mpm.A);
605 (*this) = mpm.B;
606 Tmatrix::add(1.0, *mpm.A.matrixPtr, *this->matrixPtr);
607 return *this;
608 }
609 /* B += A */
610 template<typename Treal, typename Tmatrix>
614 Tmatrix::add(1.0, *A.matrixPtr, *this->matrixPtr);
615 return *this;
616 }
617 /* B -= A */
618 template<typename Treal, typename Tmatrix>
622 Tmatrix::add(-1.0, *A.matrixPtr, *this->matrixPtr);
623 return *this;
624 }
625 /* B += alpha * A */
626 template<typename Treal, typename Tmatrix>
630 assert(!sm.tB);
631 Tmatrix::add(sm.A, *sm.B.matrixPtr, *this->matrixPtr);
632 return *this;
633 }
634
635
636 /* OPERATIONS INVOLVING SYMMETRIC MATRICES */
637 /* C = alpha * A * B : A is symmetric */
638 template<typename Treal, typename Tmatrix>
641 (const XYZ<Treal,
644 assert(this != &smm.C);
645 assert(!smm.tB && !smm.tC);
646 this->matrixPtr.haveDataStructureSet(true);
647 Tmatrix::symm('L', 'U', smm.A,
648 *smm.B.matrixPtr, *smm.C.matrixPtr,
649 0, *this->matrixPtr);
650 return *this;
651 }
652
653 /* C = A * B : A is symmetric */
654 template<typename Treal, typename Tmatrix>
659 assert(this != &mm.B);
660 assert(!mm.tB);
661 Tmatrix::symm('L', 'U', 1.0,
662 *mm.A.matrixPtr, *mm.B.matrixPtr,
663 0, *this->matrixPtr);
664 return *this;
665 }
666
667 /* C += alpha * A * B : A is symmetric */
668 template<typename Treal, typename Tmatrix>
671 (const XYZ<Treal,
674 assert(this != &smm.C);
675 assert(!smm.tB && !smm.tC);
676 Tmatrix::symm('L', 'U', smm.A,
677 *smm.B.matrixPtr, *smm.C.matrixPtr,
678 1, *this->matrixPtr);
679 return *this;
680 }
681
682 /* C = alpha * A * B + beta * C : A is symmetric */
683 template<typename Treal, typename Tmatrix>
686 (const XYZpUV<Treal,
689 Treal,
691 assert(this != &smmpsm.C);
692 assert(!smmpsm.tB && !smmpsm.tC && !smmpsm.tE);
693 if (this == &smmpsm.E)
694 Tmatrix::symm('L', 'U', smmpsm.A,
695 *smmpsm.B.matrixPtr, *smmpsm.C.matrixPtr,
696 smmpsm.D, *this->matrixPtr);
697 else
698 throw Failure("MatrixGeneral<Treal, Tmatrix>::operator="
699 "(const XYZpUV<Treal, MatrixGeneral"
700 "<Treal, Tmatrix>, MatrixSymmetric<Treal, "
701 "Tmatrix>, Treal, MatrixGeneral"
702 "<Treal, Tmatrix> >&) "
703 ": D = alpha * A * B + beta * C (with A symmetric)"
704 " not supported for C != D");
705 return *this;
706 }
707
708 /* C = alpha * B * A : A is symmetric */
709 template<typename Treal, typename Tmatrix>
712 (const XYZ<Treal,
715 assert(this != &smm.B);
716 assert(!smm.tB && !smm.tC);
717 this->matrixPtr.haveDataStructureSet(true);
718 Tmatrix::symm('R', 'U', smm.A,
719 *smm.C.matrixPtr, *smm.B.matrixPtr,
720 0, *this->matrixPtr);
721 return *this;
722 }
723
724 /* C = B * A : A is symmetric */
725 template<typename Treal, typename Tmatrix>
730 assert(this != &mm.A);
731 assert(!mm.tA && !mm.tB);
732 Tmatrix::symm('R', 'U', 1.0,
733 *mm.B.matrixPtr, *mm.A.matrixPtr,
734 0, *this->matrixPtr);
735 return *this;
736 }
737
738 /* C += alpha * B * A : A is symmetric */
739 template<typename Treal, typename Tmatrix>
742 (const XYZ<Treal,
745 assert(this != &smm.B);
746 assert(!smm.tB && !smm.tC);
747 Tmatrix::symm('R', 'U', smm.A,
748 *smm.C.matrixPtr, *smm.B.matrixPtr,
749 1, *this->matrixPtr);
750 return *this;
751 }
752
753 /* C = alpha * B * A + beta * C : A is symmetric */
754 template<typename Treal, typename Tmatrix>
757 (const XYZpUV<Treal,
760 Treal,
762 assert(this != &smmpsm.B);
763 assert(!smmpsm.tB && !smmpsm.tC && !smmpsm.tE);
764 if (this == &smmpsm.E)
765 Tmatrix::symm('R', 'U', smmpsm.A,
766 *smmpsm.C.matrixPtr, *smmpsm.B.matrixPtr,
767 smmpsm.D, *this->matrixPtr);
768 else
769 throw Failure("MatrixGeneral<Treal, Tmatrix>::operator="
770 "(const XYZpUV<Treal, MatrixSymmetric"
771 "<Treal, Tmatrix>, MatrixGeneral<Treal, "
772 "Tmatrix>, Treal, MatrixGeneral"
773 "<Treal, Tmatrix> >&) "
774 ": D = alpha * B * A + beta * C (with A symmetric)"
775 " not supported for C != D");
776 return *this;
777 }
778
779
781 template<typename Treal, typename Tmatrix>
784 (const XYZ<Treal,
787 assert(!smm.tB && !smm.tC);
788 this->matrixPtr.haveDataStructureSet(true);
789 Tmatrix::ssmm(smm.A,
790 *smm.B.matrixPtr,
791 *smm.C.matrixPtr,
792 0,
793 *this->matrixPtr);
794 return *this;
795 }
796
798 template<typename Treal, typename Tmatrix>
803 assert(!mm.tB);
804 Tmatrix::ssmm(1.0,
805 *mm.A.matrixPtr,
806 *mm.B.matrixPtr,
807 0,
808 *this->matrixPtr);
809 return *this;
810 }
811
813 template<typename Treal, typename Tmatrix>
816 (const XYZ<Treal,
819 assert(!smm.tB && !smm.tC);
820 Tmatrix::ssmm(smm.A,
821 *smm.B.matrixPtr,
822 *smm.C.matrixPtr,
823 1,
824 *this->matrixPtr);
825 return *this;
826 }
827
828
830 template<typename Treal, typename Tmatrix>
833 (const XYZpUV<Treal,
836 Treal,
838 assert(!smmpsm.tB && !smmpsm.tC && !smmpsm.tE);
839 if (this == &smmpsm.E)
840 Tmatrix::ssmm(smmpsm.A,
841 *smmpsm.B.matrixPtr,
842 *smmpsm.C.matrixPtr,
843 smmpsm.D,
844 *this->matrixPtr);
845 else
846 throw Failure("MatrixGeneral<Treal, Tmatrix>::"
847 "operator=(const XYZpUV<"
848 "Treal, MatrixSymmetric<Treal, Tmatrix>, "
849 "MatrixSymmetric<Treal, Tmatrix>, Treal, "
850 "MatrixGeneral<Treal, Tmatrix> >&) "
851 ": D = alpha * A * B + beta * C (with A and B symmetric)"
852 " not supported for C != D");
853 return *this;
854 }
855
856
857
858 /* OPERATIONS INVOLVING UPPER TRIANGULAR MATRICES */
859
860 /* B = alpha * op(A) * B : A is upper triangular */
861 template<typename Treal, typename Tmatrix>
864 (const XYZ<Treal,
867 assert(!smm.tC);
868 if (this == &smm.C)
869 Tmatrix::trmm('L', 'U', smm.tB, smm.A,
870 *smm.B.matrixPtr, *this->matrixPtr);
871 else
872 throw Failure("MatrixGeneral<Treal, Tmatrix>::operator="
873 "(const XYZ<Treal, MatrixTriangular"
874 "<Treal, Tmatrix>, MatrixGeneral<Treal,"
875 " Tmatrix> >& : D = alpha * op(A) * B (with"
876 " A upper triangular) not supported for B != D");
877 return *this;
878 }
879
880
881 /* A = alpha * A * op(B) : B is upper triangular */
882 template<typename Treal, typename Tmatrix>
885 (const XYZ<Treal,
888 assert(!smm.tB);
889 if (this == &smm.B)
890 Tmatrix::trmm('R', 'U', smm.tC, smm.A,
891 *smm.C.matrixPtr, *this->matrixPtr);
892 else
893 throw Failure("MatrixGeneral<Treal, Tmatrix>::operator="
894 "(const XYZ<Treal, MatrixGeneral"
895 "<Treal, Tmatrix>, MatrixTriangular<Treal,"
896 " Tmatrix> >& : D = alpha * A * op(B) (with"
897 " B upper triangular) not supported for A != D");
898 return *this;
899 }
900
901
902 /******* FUNCTIONS DECLARED OUTSIDE CLASS */
903 template<typename Treal, typename Tmatrix>
904 Treal trace(const XYZ<Treal,
907 if ((!smm.tB && !smm.tC) || (smm.tB && smm.tC))
909 trace_ab(smm.B, smm.C);
910 else if (smm.tB)
912 trace_aTb(smm.B, smm.C);
913 else
915 trace_aTb(smm.C, smm.B);
916 }
917
918
919} /* end namespace mat */
920#endif
921
922
Base class for matrix API.
Treal run(Treal const threshold)
Definition truncation.h:80
Truncation of general matrices.
Definition truncation.h:272
Definition Failure.h:57
Definition Interval.h:46
Treal low() const
Definition Interval.h:144
Treal midPoint() const
Definition Interval.h:115
Treal length() const
Returns the length of the interval.
Definition Interval.h:109
Base class for matrix API.
Definition MatrixBase.h:69
int get_ncols() const
Definition MatrixBase.h:141
static void getPermutedIndexes(std::vector< int > const &index, std::vector< int > const &permutation, std::vector< int > &newIndex)
Definition MatrixBase.h:205
void read_from_buffer_base(void *buffer, const int n_bytes, const matrix_type mattype)
Definition MatrixBase.h:281
int get_nrows() const
Definition MatrixBase.h:138
void writeToFileBase(std::ofstream &file, matrix_type const mattype) const
Definition MatrixBase.h:221
MatrixBase< Treal, Tmatrix > & operator=(const MatrixBase< Treal, Tmatrix > &other)
Definition MatrixBase.h:166
void resetSizesAndBlocks(SizesAndBlocks const &newRows, SizesAndBlocks const &newCols)
Definition MatrixBase.h:76
void write_to_buffer_base(void *buffer, const int n_bytes, const matrix_type mattype) const
Definition MatrixBase.h:254
ValidPtr< Tmatrix > matrixPtr
Definition MatrixBase.h:153
void readFromFileBase(std::ifstream &file, matrix_type const mattype)
Definition MatrixBase.h:236
Normal matrix.
Definition MatrixGeneral.h:59
MatrixGeneral< Treal, Tmatrix > & operator=(const MatrixGeneral< Treal, Tmatrix > &mat)
Definition MatrixGeneral.h:255
MatrixGeneral(const MatrixTriangular< Treal, Tmatrix > &triang)
Copy from triangular matrix constructor
Definition MatrixGeneral.h:71
Treal eucl_thresh(Treal const threshold)
Definition MatrixGeneral.h:525
void thresh(Treal const threshold, normType const norm)
Definition MatrixGeneral.h:509
void fullMatrix(std::vector< Treal > &fullMat) const
Definition MatrixGeneral.h:95
MatrixGeneral(const MatrixGeneral< Treal, Tmatrix > &matr)
Copy constructor
Definition MatrixGeneral.h:65
MatrixGeneral()
Default constructor
Definition MatrixGeneral.h:63
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 MatrixGeneral.h:204
Treal frob() const
Definition MatrixGeneral.h:284
std::string obj_type_id() const
Definition MatrixGeneral.h:462
MatrixGeneral< Treal, Tmatrix > & operator=(const MatrixTriangular< Treal, Tmatrix > &triang)
Definition MatrixGeneral.h:275
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 not assuming that sizes and blocks are set.
Definition MatrixGeneral.h:157
void frob_thresh(Treal threshold)
Definition MatrixGeneral.h:299
void write_to_buffer(void *buffer, const int n_bytes) const
Definition MatrixGeneral.h:324
void setElementsByRule(TRule &rule)
Definition MatrixGeneral.h:457
void assignFromFull(std::vector< Treal > const &fullMat)
Definition MatrixGeneral.h:90
void random()
Definition MatrixGeneral.h:448
MatrixGeneral< Treal, Tmatrix > & operator=(const MatrixSymmetric< Treal, Tmatrix > &symm)
Definition MatrixGeneral.h:269
void readFromFileProt(std::ifstream &file)
Read object from file.
Definition MatrixGeneral.h:467
static Treal trace_ab(const MatrixGeneral< Treal, Tmatrix > &A, const MatrixGeneral< Treal, Tmatrix > &B)
Definition MatrixGeneral.h:308
size_t nnz() const
Definition MatrixGeneral.h:317
void read_from_buffer(void *buffer, const int n_bytes)
Definition MatrixGeneral.h:327
MatrixGeneral< Treal, Tmatrix > & operator=(int const k)
Definition MatrixGeneral.h:280
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 MatrixGeneral.h:186
static Treal frob_diff(const MatrixGeneral< Treal, Tmatrix > &A, const MatrixGeneral< Treal, Tmatrix > &B)
Definition MatrixGeneral.h:288
static Treal trace_aTb(const MatrixGeneral< Treal, Tmatrix > &A, const MatrixGeneral< Treal, Tmatrix > &B)
Definition MatrixGeneral.h:313
VectorGeneral< Treal, typename Tmatrix::VectorType > VectorType
Definition MatrixGeneral.h:61
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 MatrixGeneral.h:222
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 arrays.
Definition MatrixGeneral.h:122
void writeToFileProt(std::ofstream &file) const
Write object to file.
Definition MatrixGeneral.h:464
MatrixGeneral< Treal, Tmatrix > & operator=(const Xtrans< MatrixGeneral< Treal, Tmatrix > > &mt)
Definition MatrixGeneral.h:260
MatrixGeneral(const MatrixSymmetric< Treal, Tmatrix > &symm)
Copy from symmetric matrix constructor
Definition MatrixGeneral.h:67
void gershgorin(Treal &lmin, Treal &lmax)
Definition MatrixGeneral.h:304
size_t nvalues() const
Definition MatrixGeneral.h:320
Treal eucl(Treal const requestedAccuracy, int maxIter=-1) const
Definition MatrixGeneral.h:476
void fullMatrix(std::vector< Treal > &fullMat, std::vector< int > const &rowInversePermutation, std::vector< int > const &colInversePermutation) const
Save matrix as full matrix.
Definition MatrixGeneral.h:100
void randomZeroStructure(Treal probabilityBeingZero)
Definition MatrixGeneral.h:452
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 MatrixGeneral.h:173
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 MatrixGeneral.h:138
Symmetric matrix.
Definition MatrixSymmetric.h:68
Upper non-unit triangular matrix.
Definition MatrixTriangular.h:59
Describes dimensions of matrix and its blocks on all levels.
Definition SizesAndBlocks.h:45
Definition VectorGeneral.h:48
void resetSizesAndBlocks(SizesAndBlocks const &newRows)
Definition VectorGeneral.h:51
void rand()
Definition VectorGeneral.h:108
Definition LanczosLargestMagnitudeEig.h:47
void getLargestMagnitudeEig(Treal &ev, Treal &accuracy)
Definition LanczosLargestMagnitudeEig.h:58
mat::SizesAndBlocks cols
Definition test.cc:52
#define B
#define A
Definition allocate.cc:39
@ matrix_matr
Definition MatrixBase.h:56
normType
Definition matInclude.h:139
@ frobNorm
Definition matInclude.h:139
Xtrans< TX > transpose(TX const &A)
Transposition.
Definition matrix_proxy.h:131
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
Treal trace(const XYZ< Treal, MatrixGeneral< Treal, Tmatrix >, MatrixGeneral< Treal, Tmatrix > > &smm)
Definition MatrixGeneral.h:904
normalMatrix MatrixGeneral
Definition random_matrices.h:74
Definition mat_utils.h:95
This proxy expresses the result of multiplication of three objects, of possibly different types,...
Definition matrix_proxy.h:67
This proxy expresses the result of multiplication of three objects added to two other multiplied obje...
Definition matrix_proxy.h:88
This proxy expresses the result of multiplication of two objects, of possibly different types,...
Definition matrix_proxy.h:51
This proxy expresses the result of addition of two objects, of possibly different types,...
Definition matrix_proxy.h:246
This proxy expresses the result of transposition of an object of type TX.
Definition matrix_proxy.h:118