38 #ifndef MAT_MATRIXGENERAL 39 #define MAT_MATRIXGENERAL 58 template<
typename Treal,
typename Tmatrix>
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);
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);
90 (std::vector<Treal>
const & fullMat) {
95 inline void fullMatrix(std::vector<Treal> & fullMat)
const {
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;
107 rowInversePermutation,
108 colInversePermutation);
110 for (
unsigned int ind = 0; ind < values.size(); ++ind)
111 fullMat[rowind[ind] + this->
get_nrows() * colind[ind]] =
122 (std::vector<int>
const & rowind,
123 std::vector<int>
const & colind,
124 std::vector<Treal>
const & values,
128 this->
matrixPtr->assignFromSparse(rowind, colind, values);
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;
149 this->
matrixPtr->assignFromSparse(newRowind, newColind, values);
157 (std::vector<int>
const & rowind,
158 std::vector<int>
const & colind,
159 std::vector<Treal>
const & values,
162 std::vector<int>
const & rowPermutation,
163 std::vector<int>
const & colPermutation) {
166 rowPermutation, colPermutation);
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);
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;
197 this->
matrixPtr->getValues(newRowind, newColind, values);
204 (std::vector<int> & rowind,
205 std::vector<int> & colind,
206 std::vector<Treal> & values)
const {
210 this->
matrixPtr->getAllValues(rowind, colind, 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());
232 this->
matrixPtr->getAllValues(tmpRowind, tmpColind, values);
249 inline void fullmatrix(Treal*
const full,
250 const int totnrows,
const int totncols)
const {
251 this->
matrixPtr->fullmatrix(full, totnrows, totncols);
290 return Tmatrix::frobDiff(*
A.matrixPtr, *
B.matrixPtr);
292 Treal
eucl(Treal
const requestedAccuracy,
293 int maxIter = -1)
const;
296 void thresh(Treal
const threshold,
310 return Tmatrix::trace_ab(*
A.matrixPtr, *
B.matrixPtr);
315 return Tmatrix::trace_aTb(*
A.matrixPtr, *
B.matrixPtr);
317 inline size_t nnz()
const {
453 this->
matrixPtr->randomZeroStructure(probabilityBeingZero);
456 template<
typename TRule>
458 this->
matrixPtr->setElementsByRule(rule);
474 template<
typename Treal,
typename Tmatrix>
476 eucl(Treal
const requestedAccuracy,
485 maxIter = this->get_nrows() * 100;
488 lan(ata, guess, maxIter);
489 lan.setRelTol( 100 * mat::getMachineEpsilon<Treal>() );
496 if ( euclInt.
low() < 0 )
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.");
507 template<
typename Treal,
typename Tmatrix>
513 this->frob_thresh(threshold);
516 throw Failure(
"MatrixGeneral<Treal, Tmatrix>::" 517 "thresh(Treal const, " 519 "Thresholding not imlpemented for selected norm");
523 template<
typename Treal,
typename Tmatrix>
527 return TruncObj.
run( threshold );
533 template<
typename Treal,
typename Tmatrix>
539 assert(
this != &smm.B &&
this != &smm.C);
540 this->matrixPtr.haveDataStructureSet(
true);
542 *smm.B.matrixPtr, *smm.C.matrixPtr,
543 0, *this->matrixPtr);
548 template<
typename Treal,
typename Tmatrix>
553 assert(
this != &mm.A &&
this != &mm.B);
555 *mm.A.matrixPtr, *mm.B.matrixPtr,
556 0, *this->matrixPtr);
561 template<
typename Treal,
typename Tmatrix>
567 assert(
this != &smm.B &&
this != &smm.C);
569 *smm.B.matrixPtr, *smm.C.matrixPtr,
570 1, *this->matrixPtr);
575 template<
typename Treal,
typename Tmatrix>
583 assert(
this != &smmpsm.B &&
this != &smmpsm.C);
585 if (
this == &smmpsm.E)
587 *smmpsm.B.matrixPtr, *smmpsm.C.matrixPtr,
588 smmpsm.D, *this->matrixPtr);
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");
599 template<
typename Treal,
typename Tmatrix>
604 assert(
this != &mpm.A);
606 Tmatrix::add(1.0, *mpm.A.matrixPtr, *this->matrixPtr);
610 template<
typename Treal,
typename Tmatrix>
614 Tmatrix::add(1.0, *
A.matrixPtr, *this->matrixPtr);
618 template<
typename Treal,
typename Tmatrix>
622 Tmatrix::add(-1.0, *
A.matrixPtr, *this->matrixPtr);
626 template<
typename Treal,
typename Tmatrix>
631 Tmatrix::add(sm.A, *sm.B.matrixPtr, *this->matrixPtr);
638 template<
typename Treal,
typename Tmatrix>
644 assert(
this != &smm.C);
645 assert(!smm.tB && !smm.tC);
646 this->matrixPtr.haveDataStructureSet(
true);
648 *smm.B.matrixPtr, *smm.C.matrixPtr,
649 0, *this->matrixPtr);
654 template<
typename Treal,
typename Tmatrix>
659 assert(
this != &mm.B);
662 *mm.A.matrixPtr, *mm.B.matrixPtr,
663 0, *this->matrixPtr);
668 template<
typename Treal,
typename Tmatrix>
674 assert(
this != &smm.C);
675 assert(!smm.tB && !smm.tC);
677 *smm.B.matrixPtr, *smm.C.matrixPtr,
678 1, *this->matrixPtr);
683 template<
typename Treal,
typename Tmatrix>
691 assert(
this != &smmpsm.C);
692 assert(!smmpsm.tB && !smmpsm.tC && !smmpsm.tE);
693 if (
this == &smmpsm.E)
695 *smmpsm.B.matrixPtr, *smmpsm.C.matrixPtr,
696 smmpsm.D, *this->matrixPtr);
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");
709 template<
typename Treal,
typename Tmatrix>
715 assert(
this != &smm.B);
716 assert(!smm.tB && !smm.tC);
717 this->matrixPtr.haveDataStructureSet(
true);
719 *smm.C.matrixPtr, *smm.B.matrixPtr,
720 0, *this->matrixPtr);
725 template<
typename Treal,
typename Tmatrix>
730 assert(
this != &mm.A);
731 assert(!mm.tA && !mm.tB);
733 *mm.B.matrixPtr, *mm.A.matrixPtr,
734 0, *this->matrixPtr);
739 template<
typename Treal,
typename Tmatrix>
745 assert(
this != &smm.B);
746 assert(!smm.tB && !smm.tC);
748 *smm.C.matrixPtr, *smm.B.matrixPtr,
749 1, *this->matrixPtr);
754 template<
typename Treal,
typename Tmatrix>
762 assert(
this != &smmpsm.B);
763 assert(!smmpsm.tB && !smmpsm.tC && !smmpsm.tE);
764 if (
this == &smmpsm.E)
766 *smmpsm.C.matrixPtr, *smmpsm.B.matrixPtr,
767 smmpsm.D, *this->matrixPtr);
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");
781 template<
typename Treal,
typename Tmatrix>
787 assert(!smm.tB && !smm.tC);
788 this->matrixPtr.haveDataStructureSet(
true);
798 template<
typename Treal,
typename Tmatrix>
813 template<
typename Treal,
typename Tmatrix>
819 assert(!smm.tB && !smm.tC);
830 template<
typename Treal,
typename Tmatrix>
838 assert(!smmpsm.tB && !smmpsm.tC && !smmpsm.tE);
839 if (
this == &smmpsm.E)
840 Tmatrix::ssmm(smmpsm.A,
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");
861 template<
typename Treal,
typename Tmatrix>
870 *smm.B.matrixPtr, *this->matrixPtr);
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");
882 template<
typename Treal,
typename Tmatrix>
891 *smm.C.matrixPtr, *this->matrixPtr);
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");
903 template<
typename Treal,
typename Tmatrix>
907 if ((!smm.tB && !smm.tC) || (smm.tB && smm.tC))
Normal matrix.
Definition: MatrixBase.h:49
normalMatrix MatrixGeneral
Definition: random_matrices.h:70
mat::SizesAndBlocks cols
Definition: test.cc:52
MatrixBase< Treal, Tmatrix > & operator=(const MatrixBase< Treal, Tmatrix > &other)
Definition: MatrixBase.h:166
Treal low() const
Definition: Interval.h:144
VectorGeneral< Treal, typename Tmatrix::VectorType > VectorType
Definition: MatrixGeneral.h:61
static Treal trace_ab(const MatrixGeneral< Treal, Tmatrix > &A, const MatrixGeneral< Treal, Tmatrix > &B)
Definition: MatrixGeneral.h:308
ValidPtr< Tmatrix > matrixPtr
Definition: MatrixBase.h:153
Treal eucl(Treal const requestedAccuracy, int maxIter=-1) const
Definition: MatrixGeneral.h:476
Definition: LanczosLargestMagnitudeEig.h:46
void read_from_buffer_base(void *buffer, const int n_bytes, const matrix_type mattype)
Definition: MatrixBase.h:281
MatrixGeneral(const MatrixGeneral< Treal, Tmatrix > &matr)
Copy constructor.
Definition: MatrixGeneral.h:65
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
static void gemm(const char *ta, const char *tb, const int *n, const int *k, const int *l, 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:232
MatrixGeneral< Treal, Tmatrix > & operator=(const Xtrans< MatrixGeneral< Treal, Tmatrix > > &mt)
Definition: MatrixGeneral.h:260
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
MatrixGeneral< Treal, Tmatrix > & operator=(const MatrixGeneral< Treal, Tmatrix > &mat)
Definition: MatrixGeneral.h:255
static Treal frob_diff(const MatrixGeneral< Treal, Tmatrix > &A, const MatrixGeneral< Treal, Tmatrix > &B)
Definition: MatrixGeneral.h:288
void gershgorin(Treal &lmin, Treal &lmax)
Definition: MatrixGeneral.h:304
Truncation of general matrices.
Definition: truncation.h:272
This proxy expresses the result of addition of two objects, of possibly different types...
Definition: matrix_proxy.h:246
MatrixGeneral< Treal, Tmatrix > & operator=(const MatrixSymmetric< Treal, Tmatrix > &symm)
Definition: MatrixGeneral.h:269
Definition: MatrixBase.h:55
Definition: allocate.cc:39
Describes dimensions of matrix and its blocks on all levels.
Definition: SizesAndBlocks.h:45
Treal frob() const
Definition: MatrixGeneral.h:284
MatrixGeneral< Treal, Tmatrix > & operator=(const MatrixTriangular< Treal, Tmatrix > &triang)
Definition: MatrixGeneral.h:275
Definition: Interval.h:46
void assignFromFull(std::vector< Treal > const &fullMat)
Definition: MatrixGeneral.h:90
size_t nnz() const
Definition: MatrixGeneral.h:317
Definition: MatrixBase.h:56
void readFromFileBase(std::ifstream &file, matrix_type const mattype)
Definition: MatrixBase.h:236
MatrixGeneral(const MatrixSymmetric< Treal, Tmatrix > &symm)
Copy from symmetric matrix constructor.
Definition: MatrixGeneral.h:67
int get_ncols() const
Definition: MatrixBase.h:141
This proxy expresses the result of multiplication of three objects, of possibly different types...
Definition: matrix_proxy.h:67
void write_to_buffer_base(void *buffer, const int n_bytes, const matrix_type mattype) const
Definition: MatrixBase.h:254
Upper non-unit triangular matrix.
Definition: MatrixBase.h:53
This proxy expresses the result of transposition of an object of type TX.
Definition: matrix_proxy.h:118
Treal run(Treal const threshold)
Definition: truncation.h:80
void writeToFileBase(std::ofstream &file, matrix_type const mattype) const
Definition: MatrixBase.h:221
Treal midPoint() const
Definition: Interval.h:115
void setElementsByRule(TRule &rule)
Definition: MatrixGeneral.h:457
void thresh(Treal const threshold, normType const norm)
Definition: MatrixGeneral.h:509
void rand()
Definition: VectorGeneral.h:102
size_t nvalues() const
Definition: MatrixGeneral.h:320
void resetSizesAndBlocks(SizesAndBlocks const &newRows, SizesAndBlocks const &newCols)
Definition: MatrixBase.h:76
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 resetSizesAndBlocks(SizesAndBlocks const &newRows)
Definition: VectorGeneral.h:51
std::string obj_type_id() const
Definition: MatrixGeneral.h:462
Definition: mat_utils.h:78
static void trmm(const char *side, const char *uplo, const char *transa, const char *diag, const int *m, const int *n, const T *alpha, const T *A, const int *lda, T *B, const int *ldb)
Definition: mat_gblas.h:281
Treal eucl_thresh(Treal const threshold)
Definition: MatrixGeneral.h:525
void random()
Definition: MatrixGeneral.h:448
Base class for matrix API.
Definition: MatrixBase.h:69
void getLargestMagnitudeEig(Treal &ev, Treal &accuracy)
Definition: LanczosLargestMagnitudeEig.h:58
void writeToFileProt(std::ofstream &file) const
Write object to file.
Definition: MatrixGeneral.h:464
Base class for matrix API.
MatrixGeneral()
Default constructor.
Definition: MatrixGeneral.h:63
Treal length() const
Returns the length of the interval.
Definition: Interval.h:109
Definition: matInclude.h:139
This proxy expresses the result of multiplication of two objects, of possibly different types...
Definition: matrix_proxy.h:51
void read_from_buffer(void *buffer, const int n_bytes)
Definition: MatrixGeneral.h:327
MatrixGeneral(const MatrixTriangular< Treal, Tmatrix > &triang)
Copy from triangular matrix constructor.
Definition: MatrixGeneral.h:71
void frob_thresh(Treal threshold)
Definition: MatrixGeneral.h:299
int get_nrows() const
Definition: MatrixBase.h:138
Xtrans< TX > transpose(TX const &A)
Transposition.
Definition: matrix_proxy.h:131
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
Symmetric matrix.
Definition: MatrixBase.h:51
Treal trace(const XYZ< Treal, MatrixGeneral< Treal, Tmatrix >, MatrixGeneral< Treal, Tmatrix > > &smm)
Definition: MatrixGeneral.h:904
This proxy expresses the result of multiplication of three objects added to two other multiplied obje...
Definition: matrix_proxy.h:88
MatrixGeneral< Treal, Tmatrix > & operator=(int const k)
Definition: MatrixGeneral.h:280
void readFromFileProt(std::ifstream &file)
Read object from file.
Definition: MatrixGeneral.h:467
static Treal trace_aTb(const MatrixGeneral< Treal, Tmatrix > &A, const MatrixGeneral< Treal, Tmatrix > &B)
Definition: MatrixGeneral.h:313
void fullMatrix(std::vector< Treal > &fullMat) const
Definition: MatrixGeneral.h:95
normType
Definition: matInclude.h:139
void write_to_buffer(void *buffer, const int n_bytes) const
Definition: MatrixGeneral.h:324
static void getPermutedIndexes(std::vector< int > const &index, std::vector< int > const &permutation, std::vector< int > &newIndex)
Definition: MatrixBase.h:205