77 template<
class Treal,
class Telement>
79 template<
typename Tperm>
114 template<
class Treal,
class Telement = Treal>
128 #ifdef MAT_USE_ALLOC_TIMER 132 #ifdef MAT_USE_ALLOC_TIMER 150 void fullMatrix(std::vector<Treal> & fullMat)
const;
156 std::vector<int>
const & colind,
157 std::vector<Treal>
const & values);
159 std::vector<int>
const & colind,
160 std::vector<Treal>
const & values,
161 std::vector<int>
const & indexes);
163 void addValues(std::vector<int>
const & rowind,
164 std::vector<int>
const & colind,
165 std::vector<Treal>
const & values);
166 void addValues(std::vector<int>
const & rowind,
167 std::vector<int>
const & colind,
168 std::vector<Treal>
const & values,
169 std::vector<int>
const & indexes);
172 std::vector<int>
const & colind,
173 std::vector<Treal>
const & values);
176 std::vector<int>
const & colind,
177 std::vector<Treal>
const & values);
179 void getValues(std::vector<int>
const & rowind,
180 std::vector<int>
const & colind,
181 std::vector<Treal> & values)
const;
182 void getValues(std::vector<int>
const & rowind,
183 std::vector<int>
const & colind,
184 std::vector<Treal> &,
185 std::vector<int>
const & indexes)
const;
187 std::vector<int>
const & colind,
188 std::vector<Treal> & values)
const;
190 std::vector<int> & colind,
191 std::vector<Treal> &)
const;
193 std::vector<int> & colind,
194 std::vector<Treal> &)
const;
219 template<
typename TRule>
221 template<
typename TRule>
223 template<
typename TRule>
244 static void gemm(
const bool tA,
const bool tB,
const Treal alpha,
249 static void symm(
const char side,
const char uplo,
const Treal alpha,
254 static void syrk(
const char uplo,
const bool tA,
const Treal alpha,
259 static void sysq(
const char uplo,
const Treal alpha,
264 static void ssmm(
const Treal alpha,
278 static void trmm(
const char side,
const char uplo,
const bool tA,
323 static void add(
const Treal alpha,
326 void assign(Treal
const alpha,
343 (Treal
const threshold,
353 (
Matrix<Treal, Matrix<Treal, Telement> >
const &
A );
356 (
Matrix<Treal, Matrix<Treal, Telement> >
const &
A );
362 (
Matrix<Treal, Matrix<Treal, Telement> >
const &
A,
363 Matrix<Treal, Matrix<Treal, Telement> >
const &
B );
368 (
Matrix<Treal, Matrix<Treal, Telement> >
const &
A,
369 Matrix<Treal, Matrix<Treal, Telement> >
const &
B );
380 const Matrix<Treal, Telement>&
A,
381 const Matrix<Treal, Telement>&
B,
383 Matrix<Treal, Telement>& C);
385 Matrix<Treal, Telement>&
A,
386 const Matrix<Treal, Telement>& Z);
388 const bool tA,
const Treal alpha,
389 const Matrix<Treal, Telement>&
A,
390 Matrix<Treal, Telement>&
B);
392 const Matrix<Treal, Telement>& Z,
393 Matrix<Treal, Telement>&
A);
396 (Treal
const threshold,
407 (Treal
const threshold,
416 const Treal threshold = 0,
420 void gershgorin(Treal& lmin, Treal& lmax)
const;
456 template<
typename Top>
461 for (
int row = 0; row <
col; row++)
469 template<
typename Top>
474 for (
int row = 0; row < this->
nrows(); row++)
480 static inline unsigned int level() {
return Telement::level() + 1;}
486 Treal maxAbsGlobal = 0;
487 Treal maxAbsLocal = 0;
488 for (
int ind = 0; ind < this->
nElements(); ++ind) {
489 maxAbsLocal = this->
elements[ind].maxAbsValue();
490 maxAbsGlobal = maxAbsGlobal > maxAbsLocal ?
491 maxAbsGlobal : maxAbsLocal;
504 template<
class Treal,
class Telement>
509 assert((
int)fullMat.size() == nTotalRows * nTotalCols);
512 for (
int col = 0; col < this->ncols(); col++)
513 for (
int row = 0; row < this->nrows(); row++)
514 (*
this)(row, col).assignFromFull(fullMat);
517 template<
class Treal,
class Telement>
522 fullMat.resize(nTotalRows * nTotalCols);
523 if (this->is_zero()) {
526 for (
int col = 0; col < this->nScalarsCols(); col++)
527 for (
int row = 0; row < this->nScalarsRows(); row++)
528 fullMat[(colOffset + col) * nTotalRows + (rowOffset + row)] = 0;
531 for (
int col = 0; col < this->ncols(); col++)
532 for (
int row = 0; row < this->nrows(); row++)
533 (*
this)(row, col).fullMatrix(fullMat);
537 template<
class Treal,
class Telement>
542 fullMat.resize(nTotalRows * nTotalCols);
543 if (this->is_zero()) {
546 for (
int col = 0; col < this->nScalarsCols(); col++)
547 for (
int row = 0; row <= col; row++) {
548 fullMat[(colOffset + col) * nTotalRows + (rowOffset + row)] = 0;
549 fullMat[(rowOffset + row) * nTotalRows + (colOffset + col)] = 0;
553 for (
int col = 0; col < this->ncols(); col++) {
554 for (
int row = 0; row < col; row++)
555 (*
this)(row, col).syUpTriFullMatrix(fullMat);
556 (*this)(col, col).syFullMatrix(fullMat);
561 template<
class Treal,
class Telement>
566 fullMat.resize(nTotalRows * nTotalCols);
567 if (this->is_zero()) {
570 for (
int col = 0; col < this->nScalarsCols(); col++)
571 for (
int row = 0; row <= this->nScalarsRows(); row++) {
572 fullMat[(colOffset + col) * nTotalRows + (rowOffset + row)] = 0;
573 fullMat[(rowOffset + row) * nTotalRows + (colOffset + col)] = 0;
577 for (
int col = 0; col < this->ncols(); col++)
578 for (
int row = 0; row < this->nrows(); row++)
579 (*
this)(row, col).syUpTriFullMatrix(fullMat);
586 template<
class Treal,
class Telement>
589 std::vector<int>
const & colind,
590 std::vector<Treal>
const & values) {
591 assert(rowind.size() == colind.size() &&
592 rowind.size() == values.size());
593 std::vector<int> indexes(values.size());
594 for (
unsigned int ind = 0; ind < values.size(); ++ind)
596 assignFromSparse(rowind, colind, values, indexes);
599 template<
class Treal,
class Telement>
602 std::vector<int>
const & colind,
603 std::vector<Treal>
const & values,
604 std::vector<int>
const & indexes) {
605 if (indexes.empty()) {
612 std::vector<std::vector<int> > columnBuckets (this->
cols.
getNBlocks());
613 std::vector<std::vector<int> > rowBuckets (this->
rows.
getNBlocks());
617 std::vector<int>::const_iterator it;
618 for ( it = indexes.begin(); it < indexes.end(); it++ )
619 columnBuckets[ this->
cols.
whichBlock( colind[*it] ) ].push_back (*it);
622 for (
int col = 0; col < this->ncols(); col++) {
624 while (!columnBuckets[col].empty()) {
625 currentInd = columnBuckets[col].back();
626 columnBuckets[col].pop_back();
628 ( rowind[currentInd] ) ].push_back (currentInd);
631 for (
int row = 0; row < this->nrows(); row++) {
632 (*this)(row,col).assignFromSparse(rowind, colind, values, rowBuckets[row]);
633 rowBuckets[row].clear();
638 template<
class Treal,
class Telement>
641 std::vector<int>
const & colind,
642 std::vector<Treal>
const & values) {
643 assert(rowind.size() == colind.size() &&
644 rowind.size() == values.size());
645 std::vector<int> indexes(values.size());
646 for (
unsigned int ind = 0; ind < values.size(); ++ind)
648 addValues(rowind, colind, values, indexes);
651 template<
class Treal,
class Telement>
654 std::vector<int>
const & colind,
655 std::vector<Treal>
const & values,
656 std::vector<int>
const & indexes) {
662 std::vector<std::vector<int> > columnBuckets (this->
cols.
getNBlocks());
663 std::vector<std::vector<int> > rowBuckets (this->
rows.
getNBlocks());
666 std::vector<int>::const_iterator it;
667 for ( it = indexes.begin(); it < indexes.end(); it++ )
668 columnBuckets[ this->
cols.
whichBlock( colind[*it] ) ].push_back (*it);
671 for (
int col = 0; col < this->ncols(); col++) {
673 while (!columnBuckets[col].empty()) {
674 currentInd = columnBuckets[col].back();
675 columnBuckets[col].pop_back();
677 ( rowind[currentInd] ) ].push_back (currentInd);
680 for (
int row = 0; row < this->nrows(); row++) {
681 (*this)(row,col).addValues(rowind, colind, values, rowBuckets[row]);
682 rowBuckets[row].clear();
687 template<
class Treal,
class Telement>
690 std::vector<int>
const & colind,
691 std::vector<Treal>
const & values) {
692 assert(rowind.size() == colind.size() &&
693 rowind.size() == values.size());
694 bool upperTriangleOnly =
true;
695 for (
unsigned int ind = 0; ind < values.size(); ++ind) {
697 upperTriangleOnly && colind[ind] >= rowind[ind];
699 if (!upperTriangleOnly)
700 throw Failure(
"Matrix<Treal, Telement>::" 701 "syAddValues(std::vector<int> const &, " 702 "std::vector<int> const &, " 703 "std::vector<Treal> const &, int const): " 704 "Only upper triangle can contain elements when assigning " 705 "symmetric or triangular matrix ");
706 assignFromSparse(rowind, colind, values);
709 template<
class Treal,
class Telement>
712 std::vector<int>
const & colind,
713 std::vector<Treal>
const & values) {
714 assert(rowind.size() == colind.size() &&
715 rowind.size() == values.size());
716 bool upperTriangleOnly =
true;
717 for (
unsigned int ind = 0; ind < values.size(); ++ind) {
719 upperTriangleOnly && colind[ind] >= rowind[ind];
721 if (!upperTriangleOnly)
722 throw Failure(
"Matrix<Treal, Telement>::" 723 "syAddValues(std::vector<int> const &, " 724 "std::vector<int> const &, " 725 "std::vector<Treal> const &, int const): " 726 "Only upper triangle can contain elements when assigning " 727 "symmetric or triangular matrix ");
728 addValues(rowind, colind, values);
731 template<
class Treal,
class Telement>
734 std::vector<int>
const & colind,
735 std::vector<Treal> & values)
const {
736 assert(rowind.size() == colind.size());
737 values.resize(rowind.size(),0);
738 std::vector<int> indexes(rowind.size());
739 for (
unsigned int ind = 0; ind < rowind.size(); ++ind)
741 getValues(rowind, colind, values, indexes);
744 template<
class Treal,
class Telement>
747 std::vector<int>
const & colind,
748 std::vector<Treal> & values,
749 std::vector<int>
const & indexes)
const {
750 assert(!this->is_empty());
753 std::vector<int>::const_iterator it;
754 if (this->is_zero()) {
755 for ( it = indexes.begin(); it < indexes.end(); it++ )
760 std::vector<std::vector<int> > columnBuckets (this->
cols.
getNBlocks());
761 std::vector<std::vector<int> > rowBuckets (this->
rows.
getNBlocks());
764 for ( it = indexes.begin(); it < indexes.end(); it++ )
765 columnBuckets[ this->
cols.
whichBlock( colind[*it] ) ].push_back (*it);
768 for (
int col = 0; col < this->ncols(); col++) {
770 while (!columnBuckets[col].empty()) {
771 currentInd = columnBuckets[col].back();
772 columnBuckets[col].pop_back();
774 ( rowind[currentInd] ) ].push_back (currentInd);
777 for (
int row = 0; row < this->nrows(); row++) {
778 (*this)(row,col).getValues(rowind, colind, values, rowBuckets[row]);
779 rowBuckets[row].clear();
784 template<
class Treal,
class Telement>
787 std::vector<int>
const & colind,
788 std::vector<Treal> & values)
const {
789 assert(rowind.size() == colind.size());
790 bool upperTriangleOnly =
true;
791 for (
int unsigned ind = 0; ind < rowind.size(); ++ind) {
793 upperTriangleOnly && colind[ind] >= rowind[ind];
795 if (!upperTriangleOnly)
796 throw Failure(
"Matrix<Treal, Telement>::" 797 "syGetValues(std::vector<int> const &, " 798 "std::vector<int> const &, " 799 "std::vector<Treal> const &, int const): " 800 "Only upper triangle when retrieving elements from " 801 "symmetric or triangular matrix ");
802 getValues(rowind, colind, values);
806 template<
class Treal,
class Telement>
809 std::vector<int> & colind,
810 std::vector<Treal> & values)
const {
811 assert(rowind.size() == colind.size() &&
812 rowind.size() == values.size());
813 if (!this->is_zero())
814 for (
int ind = 0; ind < this->nElements(); ++ind)
815 this->elements[ind].getAllValues(rowind, colind, values);
818 template<
class Treal,
class Telement>
821 std::vector<int> & colind,
822 std::vector<Treal> & values)
const {
823 assert(rowind.size() == colind.size() &&
824 rowind.size() == values.size());
825 if (!this->is_zero())
826 for (
int col = 0; col < this->ncols(); ++col) {
827 for (
int row = 0; row < col; ++row)
828 (*
this)(row, col).getAllValues(rowind, colind, values);
829 (*this)(col, col).syGetAllValues(rowind, colind, values);
834 template<
class Treal,
class Telement>
836 if (!this->is_zero())
837 for (
int i = 0; i < this->nElements(); i++)
838 this->elements[i].clear();
843 template<
class Treal,
class Telement>
848 if (this->is_zero()) {
849 char * tmp = (
char*)&ZERO;
850 file.write(tmp,
sizeof(
int));
853 char * tmp = (
char*)&ONE;
854 file.write(tmp,
sizeof(
int));
855 for (
int i = 0; i < this->nElements(); i++)
856 this->elements[i].writeToFile(file);
859 template<
class Treal,
class Telement>
864 char tmp[
sizeof(int)];
865 file.read(tmp, (std::ifstream::pos_type)
sizeof(
int));
873 for (
int i = 0; i < this->nElements(); i++)
874 this->elements[i].readFromFile(file);
877 throw Failure(
"Matrix<Treal, Telement>::" 878 "readFromFile(std::ifstream & file):" 879 "File corruption int value not 0 or 1");
883 template<
class Treal,
class Telement>
887 for (
int ind = 0; ind < this->nElements(); ind++)
888 this->elements[ind].random();
891 template<
class Treal,
class Telement>
896 for (
int col = 1; col < this->ncols(); col++)
897 for (
int row = 0; row < col; row++)
898 (*
this)(row, col).random();
900 for (
int rc = 0; rc < this->nrows(); rc++)
901 (*
this)(rc,rc).syRandom();
904 template<
class Treal,
class Telement>
907 if (!this->highestLevel() &&
908 probabilityBeingZero > rand() / (Treal)RAND_MAX)
913 for (
int ind = 0; ind < this->nElements(); ind++)
914 this->elements[ind].randomZeroStructure(probabilityBeingZero);
918 template<
class Treal,
class Telement>
921 if (!this->highestLevel() &&
922 probabilityBeingZero > rand() / (Treal)RAND_MAX)
928 for (
int col = 1; col < this->ncols(); col++)
929 for (
int row = 0; row < col; row++)
930 (*
this)(row, col).randomZeroStructure(probabilityBeingZero);
932 for (
int rc = 0; rc < this->nrows(); rc++)
933 (*
this)(rc,rc).syRandomZeroStructure(probabilityBeingZero);
937 template<
class Treal,
class Telement>
938 template<
typename TRule>
943 for (
int ind = 0; ind < this->nElements(); ind++)
944 this->elements[ind].setElementsByRule(rule);
946 template<
class Treal,
class Telement>
947 template<
typename TRule>
953 for (
int col = 1; col < this->ncols(); col++)
954 for (
int row = 0; row < col; row++)
955 (*
this)(row, col).setElementsByRule(rule);
957 for (
int rc = 0; rc < this->nrows(); rc++)
958 (*
this)(rc,rc).sySetElementsByRule(rule);
962 template<
class Treal,
class Telement>
965 if (this->is_empty())
966 throw Failure(
"Matrix<Treal, Telement>::addIdentity(Treal): " 967 "Cannot add identity to empty matrix.");
968 if (this->ncols() != this->nrows())
969 throw Failure(
"Matrix<Treal, Telement>::addIdentity(Treal): " 970 "Matrix must be square to add identity");
973 for (
int ind = 0; ind < this->ncols(); ind++)
974 (*
this)(ind,ind).addIdentity(alpha);
977 template<
class Treal,
class Telement>
990 AT.
elements = allocateElements<Telement>(
A.nElements());
994 for (
int row = 0; row < AT.
nrows(); row++)
995 for (
int col = 0; col < AT.
ncols(); col++)
999 template<
class Treal,
class Telement>
1003 if (this->nrows() == this->ncols()) {
1004 if (!this->is_zero()) {
1006 for (
int rc = 0; rc < this->ncols(); rc++)
1007 (*
this)(rc, rc).symToNosym();
1009 for (
int row = 1; row < this->nrows(); row++)
1010 for (
int col = 0; col < row; col++)
1015 throw Failure(
"Matrix<Treal, Telement>::symToNosym(): " 1016 "Only quadratic matrices can be symmetric");
1019 std::cout<<
"Failure caught:Matrix<Treal, Telement>::symToNosym()" 1025 template<
class Treal,
class Telement>
1028 if (this->nrows() == this->ncols()) {
1029 if (!this->is_zero()) {
1031 for (
int rc = 0; rc < this->ncols(); rc++)
1032 (*
this)(rc, rc).nosymToSym();
1034 for (
int col = 0; col < this->ncols() - 1; col++)
1035 for (
int row = col + 1; row < this->nrows(); row++)
1036 (*
this)(row, col).clear();
1040 throw Failure(
"Matrix<Treal, Telement>::nosymToSym(): " 1041 "Only quadratic matrices can be symmetric");
1046 template<
class Treal,
class Telement>
1054 if (this->ncols() != this->nrows())
1056 (
"Matrix::operator=(int k = 1): " 1057 "Matrix must be quadratic to become identity matrix.");
1060 for (
int rc = 0; rc < this->ncols(); rc++)
1064 throw Failure(
"Matrix::operator=(int k) only " 1065 "implemented for k = 0, k = 1");
1071 template<
class Treal,
class Telement>
1074 if (!this->is_zero() && alpha != 1) {
1075 for (
int ind = 0; ind < this->nElements(); ind++)
1076 this->elements[ind] *= alpha;
1082 template<
class Treal,
class Telement>
1084 gemm(
const bool tA,
const bool tB,
const Treal alpha,
1089 assert(!
A.is_empty());
1090 assert(!
B.is_empty());
1103 if ((
A.is_zero() ||
B.is_zero() || alpha == 0) &&
1107 Treal beta_tmp = beta;
1112 if (!
A.is_zero() && !
B.is_zero() && alpha != 0) {
1115 if (
A.ncols() ==
B.nrows() &&
1116 A.nrows() == C.
nrows() &&
1117 B.ncols() == C.
ncols()) {
1118 int C_ncols = C.
ncols();
1120 #pragma omp parallel for if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) schedule(dynamic) 1122 for (
int col = 0; col < C_ncols; col++) {
1124 for (
int row = 0; row < C.
nrows(); row++) {
1126 A(row, 0),
B(0, col),
1129 for(
int cola = 1; cola <
A.ncols(); cola++)
1131 A(row, cola),
B(cola, col),
1139 throw Failure(
"Matrix<class Treal, class Telement>::" 1140 "gemm(bool, bool, Treal, " 1141 "Matrix<Treal, Telement>, " 1142 "Matrix<Treal, Telement>, Treal, " 1143 "Matrix<Treal, Telement>): " 1144 "Incorrect matrixdimensions for multiplication");
1146 if (
A.nrows() ==
B.nrows() &&
1147 A.ncols() == C.
nrows() &&
1148 B.ncols() == C.
ncols()) {
1149 int C_ncols = C.
ncols();
1151 #pragma omp parallel for if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) schedule(dynamic) 1153 for (
int col = 0; col < C_ncols; col++) {
1155 for (
int row = 0; row < C.
nrows(); row++) {
1160 for(
int cola = 1; cola <
A.nrows(); cola++)
1162 A(cola, row),
B(cola, col),
1170 throw Failure(
"Matrix<class Treal, class Telement>::" 1171 "gemm(bool, bool, Treal, " 1172 "Matrix<Treal, Telement>, " 1173 "Matrix<Treal, Telement>, Treal, " 1174 "Matrix<Treal, Telement>): " 1175 "Incorrect matrixdimensions for multiplication");
1177 if (
A.ncols() ==
B.ncols() &&
1178 A.nrows() == C.
nrows() &&
1179 B.nrows() == C.
ncols()) {
1180 int C_ncols = C.
ncols();
1182 #pragma omp parallel for if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) schedule(dynamic) 1184 for (
int col = 0; col < C_ncols; col++) {
1186 for (
int row = 0; row < C.
nrows(); row++) {
1188 A(row, 0),
B(col, 0),
1191 for(
int cola = 1; cola <
A.ncols(); cola++)
1193 A(row, cola),
B(col, cola),
1201 throw Failure(
"Matrix<class Treal, class Telement>::" 1202 "gemm(bool, bool, Treal, " 1203 "Matrix<Treal, Telement>, " 1204 "Matrix<Treal, Telement>, Treal, " 1205 "Matrix<Treal, Telement>): " 1206 "Incorrect matrixdimensions for multiplication");
1208 if (
A.nrows() ==
B.ncols() &&
1209 A.ncols() == C.
nrows() &&
1210 B.nrows() == C.
ncols()) {
1211 int C_ncols = C.
ncols();
1213 #pragma omp parallel for if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) schedule(dynamic) 1215 for (
int col = 0; col < C_ncols; col++) {
1217 for (
int row = 0; row < C.
nrows(); row++) {
1219 A(0, row),
B(col, 0),
1222 for(
int cola = 1; cola <
A.nrows(); cola++)
1224 A(cola, row),
B(col, cola),
1232 throw Failure(
"Matrix<class Treal, class Telement>::" 1233 "gemm(bool, bool, Treal, " 1234 "Matrix<Treal, Telement>, " 1235 "Matrix<Treal, Telement>, " 1236 "Treal, Matrix<Treal, Telement>): " 1237 "Incorrect matrixdimensions for multiplication");
1238 else throw Failure(
"Matrix<class Treal, class Telement>::" 1239 "gemm(bool, bool, Treal, " 1240 "Matrix<Treal, Telement>, " 1241 "Matrix<Treal, Telement>, Treal, " 1242 "Matrix<Treal, Telement>):" 1243 "Very strange error!!");
1251 template<
class Treal,
class Telement>
1253 symm(
const char side,
const char uplo,
const Treal alpha,
1258 assert(!
A.is_empty());
1259 assert(!
B.is_empty());
1260 assert(
A.nrows() ==
A.ncols());
1261 int dimA =
A.nrows();
1269 assert(
side ==
'R');
1274 if ((
A.is_zero() ||
B.is_zero() || alpha == 0) &&
1279 Treal beta_tmp = beta;
1284 if (!
A.is_zero() && !
B.is_zero() && alpha != 0) {
1287 if (dimA ==
B.nrows() &&
1288 dimA == C.
nrows() &&
1289 B.ncols() == C.
ncols()) {
1290 int C_ncols = C.
ncols();
1292 #pragma omp parallel for if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) schedule(dynamic) 1294 for (
int col = 0; col < C_ncols; col++) {
1296 for (
int row = 0; row < C.
nrows(); row++) {
1299 B(row, col), beta_tmp, C(row, col));
1301 for(
int ind = 0; ind < row; ind++)
1303 A(ind, row),
B(ind, col),
1306 for(
int ind = row + 1; ind < dimA; ind++)
1308 A(row, ind),
B(ind, col),
1315 throw Failure(
"Matrix<class Treal, class Telement>" 1316 "::symm(bool, bool, Treal, Matrix<Treal, " 1317 "Telement>, Matrix<Treal, Telement>, " 1318 "Treal, Matrix<Treal, Telement>): " 1319 "Incorrect matrixdimensions for multiplication");
1321 assert(
side ==
'R');
1322 if (
B.ncols() == dimA &&
1323 B.nrows() == C.
nrows() &&
1324 dimA == C.
ncols()) {
1325 int C_ncols = C.
ncols();
1327 #pragma omp parallel for if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) schedule(dynamic) 1329 for (
int col = 0; col < C_ncols; col++) {
1331 for (
int row = 0; row < C.
nrows(); row++) {
1334 B(row, col), beta_tmp, C(row, col));
1336 for(
int ind = col + 1; ind < dimA; ind++)
1338 B(row, ind),
A(col, ind),
1341 for(
int ind = 0; ind < col; ind++)
1343 B(row, ind),
A(ind, col),
1350 throw Failure(
"Matrix<class Treal, class Telement>" 1351 "::symm(bool, bool, Treal, Matrix<Treal, " 1352 "Telement>, Matrix<Treal, Telement>, " 1353 "Treal, Matrix<Treal, Telement>): " 1354 "Incorrect matrixdimensions for multiplication");
1362 throw Failure(
"Matrix<class Treal, class Telement>::" 1363 "symm only implemented for symmetric matrices in " 1364 "upper triangular storage");
1368 template<
class Treal,
class Telement>
1370 syrk(
const char uplo,
const bool tA,
const Treal alpha,
1374 assert(!
A.is_empty());
1388 ((!tA &&
A.nrows() == C.
nrows()) || (tA &&
A.ncols() == C.
nrows())))
1389 if (alpha != 0 && !
A.is_zero()) {
1390 Treal beta_tmp = beta;
1396 if (!tA && uplo ==
'U') {
1398 #pragma omp parallel if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) 1401 int C_ncols = C.
ncols();
1403 #pragma omp for schedule(dynamic) nowait 1405 for (
int rc = 0; rc < C_ncols; rc++) {
1408 for (
int cola = 1; cola <
A.ncols(); cola++)
1412 int C_nrows = C.
nrows();
1414 #pragma omp for schedule(dynamic) nowait 1416 for (
int row = 0; row < C_nrows - 1; row++) {
1418 for (
int col = row + 1; col < C.
ncols(); col++) {
1420 A(row, 0),
A(col,0), beta_tmp, C(row,col));
1421 for (
int ind = 1; ind <
A.ncols(); ind++)
1423 A(row, ind),
A(col,ind), 1.0, C(row,col));
1429 else if (tA && uplo ==
'U') {
1431 #pragma omp parallel if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) 1434 int C_ncols = C.
ncols();
1436 #pragma omp for schedule(dynamic) nowait 1438 for (
int rc = 0; rc < C_ncols; rc++) {
1441 for (
int rowa = 1; rowa <
A.nrows(); rowa++)
1445 int C_nrows = C.
nrows();
1447 #pragma omp for schedule(dynamic) nowait 1449 for (
int row = 0; row < C_nrows - 1; row++) {
1451 for (
int col = row + 1; col < C.
ncols(); col++) {
1453 A(0, row),
A(0, col), beta_tmp, C(row,col));
1454 for (
int ind = 1; ind <
A.nrows(); ind++)
1456 A(ind, row),
A(ind, col), 1.0, C(row,col));
1463 throw Failure(
"Matrix<class Treal, class Telement>::" 1464 "syrk not implemented for lower triangular storage");
1471 throw Failure(
"Matrix<class Treal, class Telement>::syrk: " 1472 "Incorrect matrix dimensions for symmetric rank-k update");
1475 template<
class Treal,
class Telement>
1477 sysq(
const char uplo,
const Treal alpha,
1481 assert(!
A.is_empty());
1488 A.nrows() == C.
nrows() &&
A.nrows() ==
A.ncols())
1489 if (alpha != 0 && !
A.is_zero()) {
1491 Treal beta_tmp = beta;
1498 #pragma omp parallel if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) 1501 int C_ncols = C.
ncols();
1503 #pragma omp for schedule(dynamic) nowait 1505 for (
int rc = 0; rc < C_ncols; rc++) {
1507 Telement::sysq(uplo, alpha,
A(rc, rc), beta_tmp, C(rc, rc));
1508 for (
int cola = 0; cola < rc; cola++)
1510 for (
int cola = rc + 1; cola <
A.ncols(); cola++)
1515 int C_nrows = C.
nrows();
1517 #pragma omp for schedule(dynamic) nowait 1519 for (
int row = 0; row < C_nrows - 1; row++) {
1521 for (
int col = row + 1; col < C.
ncols(); col++) {
1524 beta_tmp, C(row, col));
1528 for (
int ind = 0; ind < row; ind++)
1530 A(ind, row),
A(ind, col), 1.0, C(row, col));
1532 for (
int ind = row + 1; ind < col; ind++)
1534 A(row, ind),
A(ind, col), 1.0, C(row, col));
1536 for (
int ind = col + 1; ind <
A.ncols(); ind++)
1538 A(row, ind),
A(col, ind), 1.0, C(row, col));
1546 throw Failure(
"Matrix<class Treal, class Telement>::" 1547 "sysq only implemented for symmetric matrices in " 1548 "upper triangular storage");;
1554 throw Failure(
"Matrix<class Treal, class Telement>::sysq: " 1555 "Incorrect matrix dimensions for symmetric square " 1560 template<
class Treal,
class Telement>
1567 assert(!
A.is_empty());
1568 assert(!
B.is_empty());
1574 if (
A.ncols() !=
B.nrows() ||
1575 A.nrows() != C.
nrows() ||
1576 B.ncols() != C.
ncols() ||
1577 A.nrows() !=
A.ncols() ||
1578 B.nrows() !=
B.ncols()) {
1579 throw Failure(
"Matrix<class Treal, class Telement>::" 1581 "Matrix<Treal, Telement>, " 1582 "Matrix<Treal, Telement>, Treal, " 1583 "Matrix<Treal, Telement>): " 1584 "Incorrect matrixdimensions for ssmm multiplication");
1586 if ((
A.is_zero() ||
B.is_zero() || alpha == 0) &&
1591 Treal beta_tmp = beta;
1596 if (
A.is_zero() ||
B.is_zero() || alpha == 0) {
1601 int C_ncols = C.
ncols();
1603 #pragma omp parallel for if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) schedule(dynamic) 1605 for (
int col = 0; col < C_ncols; col++) {
1609 Telement::ssmm(alpha,
A(col,col),
B(col, col),
1610 beta_tmp, C(col,col));
1611 for (
int ind = 0; ind < col; ind++)
1614 alpha,
A(ind,col),
B(ind, col),
1616 for (
int ind = col + 1; ind <
A.ncols(); ind++)
1619 alpha,
A(col, ind),
B(col, ind),
1622 for (
int row = 0; row < col; row++) {
1627 alpha,
A(row, row),
B(row, col),
1628 beta_tmp, C(row, col));
1631 alpha,
B(col, col),
A(row, col),
1633 for (
int ind = 0; ind < row; ind++)
1636 alpha,
A(ind, row),
B(ind, col),
1638 for (
int ind = row + 1; ind < col; ind++)
1641 alpha,
A(row, ind),
B(ind, col),
1643 for (
int ind = col + 1; ind <
A.ncols(); ind++)
1646 alpha,
A(row, ind),
B(col, ind),
1651 for (
int row = col + 1; row < C.
nrows(); row++) {
1657 alpha,
B(col, col),
A(col, row),
1658 beta_tmp, tmpSubMat);
1661 alpha,
A(row, row),
B(col, row),
1663 for (
int ind = 0; ind < col; ind++)
1666 alpha,
B(ind, col),
A(ind, row),
1668 for (
int ind = col + 1; ind < row; ind++)
1671 alpha,
B(col, ind),
A(ind, row),
1673 for (
int ind = row + 1; ind <
B.nrows(); ind++)
1676 alpha,
B(col, ind),
A(row, ind),
1689 template<
class Treal,
class Telement>
1696 assert(!
A.is_empty());
1697 assert(!
B.is_empty());
1703 if (
A.ncols() !=
B.nrows() ||
1704 A.nrows() != C.
nrows() ||
1705 B.ncols() != C.
ncols() ||
1706 A.nrows() !=
A.ncols() ||
1707 B.nrows() !=
B.ncols()) {
1708 throw Failure(
"Matrix<class Treal, class Telement>::" 1709 "ssmm_upper_tr_only(Treal, " 1710 "Matrix<Treal, Telement>, " 1711 "Matrix<Treal, Telement>, Treal, " 1712 "Matrix<Treal, Telement>): " 1713 "Incorrect matrixdimensions for ssmm_upper_tr_only " 1716 if ((
A.is_zero() ||
B.is_zero() || alpha == 0) &&
1721 Treal beta_tmp = beta;
1726 if (
A.is_zero() ||
B.is_zero() || alpha == 0) {
1731 int C_ncols = C.
ncols();
1733 #pragma omp parallel for if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) schedule(dynamic) 1735 for (
int col = 0; col < C_ncols; col++) {
1739 Telement::ssmm_upper_tr_only(alpha,
A(col,col),
B(col, col),
1740 beta_tmp, C(col,col));
1741 for (
int ind = 0; ind < col; ind++)
1743 Telement::gemm_upper_tr_only(
true,
false,
1744 alpha,
A(ind,col),
B(ind, col),
1746 for (
int ind = col + 1; ind <
A.ncols(); ind++)
1748 Telement::gemm_upper_tr_only(
false,
true,
1749 alpha,
A(col, ind),
B(col, ind),
1752 for (
int row = 0; row < col; row++) {
1757 alpha,
A(row, row),
B(row, col),
1758 beta_tmp, C(row, col));
1761 alpha,
B(col, col),
A(row, col),
1763 for (
int ind = 0; ind < row; ind++)
1766 alpha,
A(ind, row),
B(ind, col),
1768 for (
int ind = row + 1; ind < col; ind++)
1771 alpha,
A(row, ind),
B(ind, col),
1773 for (
int ind = col + 1; ind <
A.ncols(); ind++)
1776 alpha,
A(row, ind),
B(col, ind),
1786 template<
class Treal,
class Telement>
1788 trmm(
const char side,
const char uplo,
const bool tA,
const Treal alpha,
1791 assert(!
A.is_empty());
1792 assert(!
B.is_empty());
1793 if (alpha != 0 && !
A.is_zero() && !
B.is_zero())
1794 if (((
side ==
'R' &&
B.ncols() ==
A.nrows()) ||
1795 (
side ==
'L' &&
A.ncols() ==
B.nrows())) &&
1796 A.nrows() ==
A.ncols())
1801 for (
int col =
B.ncols() - 1; col >= 0; col--)
1802 for (
int row = 0; row <
B.nrows(); row++) {
1806 A(col, col),
B(row,col));
1808 for (
int ind = 0; ind < col; ind++)
1810 B(row,ind),
A(ind, col),
1815 assert(
side ==
'L');
1817 for (
int row = 0; row <
B.nrows(); row++ )
1818 for (
int col = 0; col <
B.ncols(); col++) {
1820 A(row,row),
B(row,col));
1821 for (
int ind = row + 1 ; ind <
B.nrows(); ind++)
1823 A(row,ind),
B(ind,col),
1832 for (
int col = 0; col <
B.ncols(); col++)
1833 for (
int row = 0; row <
B.nrows(); row++) {
1835 A(col,col),
B(row,col));
1836 for (
int ind = col + 1; ind <
A.ncols(); ind++)
1838 B(row,ind),
A(col,ind),
1843 assert(
side ==
'L');
1845 for (
int row =
B.nrows() - 1; row >= 0; row--)
1846 for (
int col = 0; col <
B.ncols(); col++) {
1848 A(row,row),
B(row,col));
1849 for (
int ind = 0; ind < row; ind++)
1851 A(ind,row),
B(ind,col),
1857 throw Failure(
"Matrix<class Treal, class Telement>::" 1858 "trmm not implemented for lower triangular matrices");
1860 throw Failure(
"Matrix<class Treal, class Telement>::trmm" 1861 ": Incorrect matrix dimensions for multiplication");
1867 template<
class Treal,
class Telement>
1869 assert(!this->is_empty());
1870 if (this->is_zero())
1874 for (
int i = 0; i < this->nElements(); i++)
1875 sum += this->elements[i].frobSquared();
1880 template<
class Treal,
class Telement>
1883 assert(!this->is_empty());
1884 if (this->nrows() != this->ncols())
1885 throw Failure(
"Matrix<Treal, Telement>::syFrobSquared(): " 1886 "Matrix must be have equal number of rows " 1887 "and cols to be symmetric");
1889 if (!this->is_zero()) {
1890 for (
int col = 1; col < this->ncols(); col++)
1891 for (
int row = 0; row < col; row++)
1892 sum += 2 * (*
this)(row, col).frobSquared();
1893 for (
int rc = 0; rc < this->ncols(); rc++)
1894 sum += (*
this)(rc, rc).syFrobSquared();
1899 template<
class Treal,
class Telement>
1904 assert(!
A.is_empty());
1905 assert(!
B.is_empty());
1906 if (
A.nrows() !=
B.nrows() ||
A.ncols() !=
B.ncols())
1907 throw Failure(
"Matrix<Treal, Telement>::" 1908 "frobSquaredDiff: Incorrect matrix dimensions");
1910 if (!
A.is_zero() && !
B.is_zero())
1911 for (
int i = 0; i <
A.nElements(); i++)
1912 sum += Telement::frobSquaredDiff(
A.elements[i],
B.elements[i]);
1913 else if (
A.is_zero() && !
B.is_zero())
1914 sum =
B.frobSquared();
1915 else if (!
A.is_zero() &&
B.is_zero())
1916 sum =
A.frobSquared();
1921 template<
class Treal,
class Telement>
1926 assert(!
A.is_empty());
1927 assert(!
B.is_empty());
1928 if (
A.nrows() !=
B.nrows() ||
1929 A.ncols() !=
B.ncols() ||
1930 A.nrows() !=
A.ncols())
1931 throw Failure(
"Matrix<Treal, Telement>::syFrobSquaredDiff: " 1932 "Incorrect matrix dimensions");
1934 if (!
A.is_zero() && !
B.is_zero()) {
1935 for (
int col = 1; col <
A.ncols(); col++)
1936 for (
int row = 0; row < col; row++)
1937 sum += 2 * Telement::frobSquaredDiff(
A(row, col),
B(row, col));
1938 for (
int rc = 0; rc <
A.ncols(); rc++)
1939 sum += Telement::syFrobSquaredDiff(
A(rc, rc),
B(rc, rc));
1941 else if (
A.is_zero() && !
B.is_zero())
1942 sum =
B.syFrobSquared();
1943 else if (!
A.is_zero() &&
B.is_zero())
1944 sum =
A.syFrobSquared();
1949 template<
class Treal,
class Telement>
1952 assert(!this->is_empty());
1953 if (this->nrows() != this->ncols())
1954 throw Failure(
"Matrix<Treal, Telement>::trace(): " 1955 "Matrix must be quadratic");
1956 if (this->is_zero())
1960 for (
int rc = 0; rc < this->ncols(); rc++)
1961 sum += (*
this)(rc,rc).
trace();
1966 template<
class Treal,
class Telement>
1970 assert(!
A.is_empty());
1971 assert(!
B.is_empty());
1972 if (
A.nrows() !=
B.ncols() ||
A.ncols() !=
B.nrows())
1973 throw Failure(
"Matrix<Treal, Telement>::trace_ab: " 1974 "Wrong matrix dimensions for trace(A * B)");
1976 if (!
A.is_zero() && !
B.is_zero())
1977 for (
int rc = 0; rc <
A.nrows(); rc++)
1978 for (
int colA = 0; colA <
A.ncols(); colA++)
1979 tr += Telement::trace_ab(
A(rc,colA),
B(colA, rc));
1983 template<
class Treal,
class Telement>
1987 assert(!
A.is_empty());
1988 assert(!
B.is_empty());
1989 if (
A.ncols() !=
B.ncols() ||
A.nrows() !=
B.nrows())
1990 throw Failure(
"Matrix<Treal, Telement>::trace_aTb: " 1991 "Wrong matrix dimensions for trace(A' * B)");
1993 if (!
A.is_zero() && !
B.is_zero())
1994 for (
int rc = 0; rc <
A.ncols(); rc++)
1995 for (
int rowB = 0; rowB <
B.nrows(); rowB++)
1996 tr += Telement::trace_aTb(
A(rowB,rc),
B(rowB, rc));
2000 template<
class Treal,
class Telement>
2004 assert(!
A.is_empty());
2005 assert(!
B.is_empty());
2006 if (
A.nrows() !=
B.ncols() ||
A.ncols() !=
B.nrows() ||
2007 A.nrows() !=
A.ncols())
2008 throw Failure(
"Matrix<Treal, Telement>::sy_trace_ab: " 2009 "Wrong matrix dimensions for symmetric trace(A * B)");
2011 if (!
A.is_zero() && !
B.is_zero()) {
2013 for (
int rc = 0; rc <
A.nrows(); rc++)
2014 tr += Telement::sy_trace_ab(
A(rc,rc),
B(rc, rc));
2016 for (
int colA = 1; colA <
A.ncols(); colA++)
2017 for (
int rowA = 0; rowA < colA; rowA++)
2018 tr += 2 * Telement::trace_aTb(
A(rowA, colA),
B(rowA, colA));
2023 template<
class Treal,
class Telement>
2028 assert(!
A.is_empty());
2029 assert(!
B.is_empty());
2030 if (
A.nrows() !=
B.nrows() ||
A.ncols() !=
B.ncols())
2031 throw Failure(
"Matrix<Treal, Telement>::add: " 2032 "Wrong matrix dimensions for addition");
2033 if (!
A.is_zero() && alpha != 0) {
2036 for (
int ind = 0; ind <
A.nElements(); ind++)
2037 Telement::add(alpha,
A.elements[ind],
B.elements[ind]);
2041 template<
class Treal,
class Telement>
2045 assert(!
A.is_empty());
2046 if (this->is_empty()) {
2047 this->resetRows(
A.rows);
2048 this->resetCols(
A.cols);
2052 add(alpha,
A, *
this);
2059 template<
class Treal,
class Telement>
2062 if (!this->is_zero())
2063 for (
int ind = 0; ind < this->nElements(); ind++)
2064 this->elements[ind].getFrobSqLowestLevel(frobsq);
2067 template<
class Treal,
class Telement>
2071 assert(!this->is_empty());
2072 bool thisMatIsZero =
true;
2075 bool EMatIsZero =
true;
2076 if (!ErrorMatrix->
is_zero() || !this->is_zero()) {
2079 if (this->is_zero())
2081 for (
int ind = 0; ind < this->nElements(); ind++) {
2082 this->elements[ind].
2083 frobThreshLowestLevel(threshold, &ErrorMatrix->
elements[ind]);
2084 thisMatIsZero = thisMatIsZero && this->elements[ind].is_zero();
2085 EMatIsZero = EMatIsZero && ErrorMatrix->
elements[ind].is_zero();
2090 ErrorMatrix->
clear();
2094 if (!this->is_zero()) {
2095 for (
int ind = 0; ind < this->nElements(); ind++) {
2096 this->elements[ind].frobThreshLowestLevel(threshold, 0);
2097 thisMatIsZero = thisMatIsZero && this->elements[ind].is_zero();
2104 template<
class Treal,
class Telement>
2107 if (!this->is_zero())
2108 for (
int ind = 0; ind < this->nElements(); ind++)
2109 this->elements[ind].getFrobSqElementLevel(frobsq);
2112 template<
class Treal,
class Telement>
2116 assert(!this->is_empty());
2117 bool thisMatIsZero =
true;
2120 bool EMatIsZero =
true;
2121 if (!ErrorMatrix->
is_zero() || !this->is_zero()) {
2124 if (this->is_zero())
2126 for (
int ind = 0; ind < this->nElements(); ind++) {
2127 this->elements[ind].
2128 frobThreshElementLevel(threshold, &ErrorMatrix->
elements[ind]);
2129 thisMatIsZero = thisMatIsZero && this->elements[ind].is_zero();
2130 EMatIsZero = EMatIsZero && ErrorMatrix->
elements[ind].is_zero();
2135 ErrorMatrix->
clear();
2139 if (!this->is_zero()) {
2140 for (
int ind = 0; ind < this->nElements(); ind++) {
2141 this->elements[ind].frobThreshElementLevel(threshold, 0);
2142 thisMatIsZero = thisMatIsZero && this->elements[ind].is_zero();
2151 template<
class Treal,
class Telement>
2155 if ( this->is_zero() )
2157 assert( this->nElements() ==
A.nElements() );
2158 for (
int ind = 0; ind < this->nElements(); ind++)
2159 this->elements[ind].assignFrobNormsLowestLevel(
A[ind]);
2165 template<
class Treal,
class Telement>
2168 assert(!
A.is_empty());
2169 if (
A.nrows() !=
A.ncols())
2170 throw Failure(
"Matrix<Treal, Telement>::syAssignFrobNormsLowestLevel(...): " 2171 "Matrix must be have equal number of rows " 2172 "and cols to be symmetric");
2174 if ( this->is_zero() )
2176 assert( this->nElements() ==
A.nElements() );
2177 for (
int col = 1; col < this->ncols(); col++)
2178 for (
int row = 0; row < col; row++)
2179 (*
this)(row, col).assignFrobNormsLowestLevel(
A(row,col) );
2180 for (
int rc = 0; rc < this->ncols(); rc++)
2181 (*
this)(rc, rc).syAssignFrobNormsLowestLevel(
A(rc,rc) );
2187 template<
class Treal,
class Telement>
2191 if (
A.is_zero() &&
B.is_zero() ) {
2197 if ( this->is_zero() )
2199 if ( !
A.is_zero() && !
B.is_zero() ) {
2201 assert( this->nElements() ==
A.nElements() );
2202 assert( this->nElements() ==
B.nElements() );
2203 for (
int ind = 0; ind < this->nElements(); ind++)
2204 this->elements[ind].assignDiffFrobNormsLowestLevel(
A[ind],
B[ind] );
2207 if ( !
A.is_zero() ) {
2209 this->assignFrobNormsLowestLevel(
A );
2212 if ( !
B.is_zero() ) {
2214 this->assignFrobNormsLowestLevel(
B );
2218 template<
class Treal,
class Telement>
2222 if (
A.is_zero() &&
B.is_zero() ) {
2228 if ( this->is_zero() )
2230 if ( !
A.is_zero() && !
B.is_zero() ) {
2232 assert( this->nElements() ==
A.nElements() );
2233 assert( this->nElements() ==
B.nElements() );
2234 for (
int col = 1; col < this->ncols(); col++)
2235 for (
int row = 0; row < col; row++)
2236 (*
this)(row, col).assignDiffFrobNormsLowestLevel(
A(row,col),
B(row,col) );
2237 for (
int rc = 0; rc < this->ncols(); rc++)
2238 (*
this)(rc, rc).syAssignDiffFrobNormsLowestLevel(
A(rc,rc),
B(rc,rc) );
2241 if ( !
A.is_zero() ) {
2243 this->syAssignFrobNormsLowestLevel(
A );
2246 if ( !
B.is_zero() ) {
2248 this->syAssignFrobNormsLowestLevel(
B );
2255 template<
class Treal,
class Telement>
2258 if ( this->is_zero() ) {
2262 assert( !
A.is_zero() );
2263 assert( this->nElements() ==
A.nElements() );
2264 for (
int ind = 0; ind < this->nElements(); ind++)
2265 this->elements[ind].truncateAccordingToSparsityPattern(
A[ind] );
2275 template<
class Treal,
class Telement>
2282 assert(!
A.is_empty());
2283 assert(!
B.is_empty());
2295 if ((
A.is_zero() ||
B.is_zero() || alpha == 0) &&
2299 Treal beta_tmp = beta;
2304 if (!
A.is_zero() && !
B.is_zero() && alpha != 0) {
2306 if (
A.ncols() ==
B.nrows() &&
2307 A.nrows() == C.
nrows() &&
2308 B.ncols() == C.
ncols()) {
2309 for (
int col = 0; col < C.
ncols(); col++) {
2310 Telement::gemm_upper_tr_only(tA, tB, alpha,
2311 A(col, 0),
B(0, col),
2314 for(
int cola = 1; cola <
A.ncols(); cola++)
2315 Telement::gemm_upper_tr_only(tA, tB, alpha,
2316 A(col, cola),
B(cola, col),
2319 for (
int row = 0; row < col; row++) {
2321 A(row, 0),
B(0, col),
2324 for(
int cola = 1; cola <
A.ncols(); cola++)
2326 A(row, cola),
B(cola, col),
2333 throw Failure(
"Matrix<class Treal, class Telement>::" 2334 "gemm_upper_tr_only(bool, bool, Treal, " 2335 "Matrix<Treal, Telement>, " 2336 "Matrix<Treal, Telement>, " 2337 "Treal, Matrix<Treal, Telement>): " 2338 "Incorrect matrixdimensions for multiplication");
2340 if (
A.nrows() ==
B.nrows() &&
2341 A.ncols() == C.
nrows() &&
2342 B.ncols() == C.
ncols()) {
2343 for (
int col = 0; col < C.
ncols(); col++) {
2344 Telement::gemm_upper_tr_only(tA, tB, alpha,
2348 for(
int cola = 1; cola <
A.nrows(); cola++)
2349 Telement::gemm_upper_tr_only(tA, tB, alpha,
2350 A(cola, col),
B(cola, col),
2353 for (
int row = 0; row < col; row++) {
2358 for(
int cola = 1; cola <
A.nrows(); cola++)
2360 A(cola, row),
B(cola, col),
2367 throw Failure(
"Matrix<class Treal, class Telement>::" 2368 "gemm_upper_tr_only(bool, bool, Treal, " 2369 "Matrix<Treal, Telement>, " 2370 "Matrix<Treal, Telement>, " 2371 "Treal, Matrix<Treal, Telement>): " 2372 "Incorrect matrixdimensions for multiplication");
2374 if (
A.ncols() ==
B.ncols() &&
2375 A.nrows() == C.
nrows() &&
2376 B.nrows() == C.
ncols()) {
2377 for (
int col = 0; col < C.
ncols(); col++) {
2378 Telement::gemm_upper_tr_only(tA, tB, alpha,
2379 A(col, 0),
B(col, 0),
2382 for(
int cola = 1; cola <
A.ncols(); cola++)
2383 Telement::gemm_upper_tr_only(tA, tB, alpha,
2384 A(col, cola),
B(col, cola),
2387 for (
int row = 0; row < col; row++) {
2389 A(row, 0),
B(col, 0),
2392 for(
int cola = 1; cola <
A.ncols(); cola++)
2394 A(row, cola),
B(col, cola),
2401 throw Failure(
"Matrix<class Treal, class Telement>::" 2402 "gemm_upper_tr_only(bool, bool, Treal, " 2403 "Matrix<Treal, Telement>, " 2404 "Matrix<Treal, Telement>, " 2405 "Treal, Matrix<Treal, Telement>): " 2406 "Incorrect matrixdimensions for multiplication");
2408 if (
A.nrows() ==
B.ncols() &&
2409 A.ncols() == C.
nrows() &&
2410 B.nrows() == C.
ncols()) {
2411 for (
int col = 0; col < C.
ncols(); col++) {
2412 Telement::gemm_upper_tr_only(tA, tB, alpha,
2413 A(0, col),
B(col, 0),
2416 for(
int cola = 1; cola <
A.nrows(); cola++)
2417 Telement::gemm_upper_tr_only(tA, tB, alpha,
2418 A(cola, col),
B(col, cola),
2421 for (
int row = 0; row < col; row++) {
2423 A(0, row),
B(col, 0),
2426 for(
int cola = 1; cola <
A.nrows(); cola++)
2428 A(cola, row),
B(col, cola),
2435 throw Failure(
"Matrix<class Treal, class Telement>::" 2436 "gemm_upper_tr_only(bool, bool, Treal, " 2437 "Matrix<Treal, Telement>, " 2438 "Matrix<Treal, Telement>, Treal, " 2439 "Matrix<Treal, Telement>): " 2440 "Incorrect matrixdimensions for multiplication");
2441 else throw Failure(
"Matrix<class Treal, class Telement>::" 2442 "gemm_upper_tr_only(bool, bool, Treal, " 2443 "Matrix<Treal, Telement>, " 2444 "Matrix<Treal, Telement>, Treal, " 2445 "Matrix<Treal, Telement>):" 2446 "Very strange error!!");
2458 template<
class Treal,
class Telement>
2463 assert(!
A.is_empty());
2465 if (alpha != 0 && !
A.is_zero() && !Z.
is_zero()) {
2466 if (
A.nrows() ==
A.ncols() &&
2468 A.nrows() == Z.
nrows()) {
2471 for (
int col =
A.ncols() - 1; col >= 0; col--) {
2473 Telement::sytr_upper_tr_only(
side, alpha,
2474 A(col, col), Z(col, col));
2475 for (
int ind = 0; ind < col; ind++) {
2477 Telement::gemm_upper_tr_only(
true,
false, alpha,
A(ind, col),
2478 Z(ind, col), 1.0,
A(col, col));
2481 for (
int row = col - 1; row >= 0; row--) {
2484 alpha, Z(col, col),
A(row, col));
2488 for (
int ind = 0; ind < row; ind++) {
2493 for (
int ind = row + 1; ind < col; ind++) {
2502 assert(
side ==
'L');
2504 for (
int col = 0; col <
A.ncols(); col++) {
2506 for (
int row = 0; row < col; row++) {
2509 Z(row, row),
A(row, col));
2513 for (
int ind = row + 1; ind < col; ind++)
2517 for (
int ind = col + 1; ind <
A.nrows(); ind++)
2522 Telement::sytr_upper_tr_only(
side, alpha,
2523 A(col, col), Z(col, col));
2524 for (
int ind = col + 1; ind <
A.ncols(); ind++)
2525 Telement::gemm_upper_tr_only(
false,
true, alpha, Z(col, ind),
2526 A(col, ind), 1.0,
A(col, col));
2531 throw Failure(
"Matrix<class Treal, class Telement>::" 2532 "sytr_upper_tr_only: Incorrect matrix dimensions " 2533 "for symmetric-triangular multiplication");
2541 template<
class Treal,
class Telement>
2544 const bool tA,
const Treal alpha,
2547 assert(!
A.is_empty());
2548 assert(!
B.is_empty());
2549 if (alpha != 0 && !
A.is_zero() && !
B.is_zero())
2550 if (((
side ==
'R' &&
B.ncols() ==
A.nrows()) ||
2551 (
side ==
'L' &&
A.ncols() ==
B.nrows())) &&
2552 A.nrows() ==
A.ncols())
2555 throw Failure(
"Matrix<Treal, class Telement>::" 2556 "trmm_upper_tr_only : " 2557 "not possible for upper triangular not transposed " 2558 "matrices A since lower triangle of B is needed");
2564 for (
int col = 0; col <
B.ncols(); col++) {
2565 Telement::trmm_upper_tr_only(
side, uplo, tA, alpha,
2566 A(col,col),
B(col,col));
2567 for (
int ind = col + 1; ind <
A.ncols(); ind++)
2568 Telement::gemm_upper_tr_only(
false, tA, alpha,
2569 B(col,ind),
A(col,ind),
2571 for (
int row = 0; row < col; row++) {
2573 A(col,col),
B(row,col));
2574 for (
int ind = col + 1; ind <
A.ncols(); ind++)
2576 B(row,ind),
A(col,ind),
2582 assert(
side ==
'L');
2584 for (
int row =
B.nrows() - 1; row >= 0; row--) {
2585 Telement::trmm_upper_tr_only(
side, uplo, tA, alpha,
2586 A(row,row),
B(row,row));
2587 for (
int ind = 0; ind < row; ind++)
2588 Telement::gemm_upper_tr_only(tA,
false, alpha,
2589 A(ind,row),
B(ind,row),
2591 for (
int col = row + 1; col <
B.ncols(); col++) {
2593 A(row,row),
B(row,col));
2594 for (
int ind = 0; ind < row; ind++)
2596 A(ind,row),
B(ind,col),
2603 throw Failure(
"Matrix<class Treal, class Telement>::" 2604 "trmm_upper_tr_only not implemented for lower " 2605 "triangular matrices");
2607 throw Failure(
"Matrix<class Treal, class Telement>::" 2608 "trmm_upper_tr_only: Incorrect matrix dimensions " 2609 "for multiplication");
2618 template<
class Treal,
class Telement>
2630 assert(
side ==
'L');
2638 template<
class Treal,
class Telement>
2642 assert(!this->is_empty());
2643 if (ErrorMatrix && ErrorMatrix->
is_empty()) {
2647 assert(threshold >= (Treal)0.0);
2648 if (threshold == (Treal)0.0)
2651 std::vector<Treal> frobsq_norms;
2652 getFrobSqLowestLevel(frobsq_norms);
2654 std::sort(frobsq_norms.begin(), frobsq_norms.end());
2655 int lastRemoved = -1;
2656 Treal frobsqSum = 0;
2657 int nnzBlocks = frobsq_norms.size();
2658 while(lastRemoved < nnzBlocks - 1 && frobsqSum < threshold) {
2660 frobsqSum += frobsq_norms[lastRemoved];
2664 if (lastRemoved == nnzBlocks - 1 && frobsqSum < threshold) {
2671 frobsqSum -= frobsq_norms[lastRemoved];
2673 while(lastRemoved > -1 && frobsq_norms[lastRemoved] ==
2674 frobsq_norms[lastRemoved + 1]) {
2675 frobsqSum -= frobsq_norms[lastRemoved];
2678 if (lastRemoved > -1) {
2679 Treal threshLowestLevel =
2680 (frobsq_norms[lastRemoved + 1] +
2681 frobsq_norms[lastRemoved]) / 2;
2682 this->frobThreshLowestLevel(threshLowestLevel, ErrorMatrix);
2688 template<
class Treal,
class Telement>
2692 const Treal threshold,
const side looking,
2694 assert(!
A.is_empty());
2695 if (
A.nrows() !=
A.ncols())
2696 throw Failure(
"Matrix<Treal, Telement>::sy_inch: " 2697 "Matrix must be quadratic!");
2699 throw Failure(
"Matrix<Treal>::sy_inch: Matrix is zero! " 2700 "Not possible to compute inverse cholesky.");
2702 throw Failure(
"Matrix<Treal>::sy_inch: Only unstable " 2703 "version of sy_inch implemented.");
2706 myThresh = threshold / (Z.
ncols() * Z.
nrows());
2711 if (looking ==
left) {
2713 throw Failure(
"Matrix<Treal, Telement>::syInch: " 2714 "Thresholding not implemented for left-looking inch.");
2716 Telement::syInch(
A(0,0), Z(0,0), threshold, looking, version);
2718 for (
int i = 1; i < Z.
ncols(); i++) {
2719 for (
int j = 0; j < i; j++) {
2723 for (
int k = 0; k < j; k++)
2725 A(k,j), Z(k,i), 1.0, Ptmp);
2728 for (
int k = 0; k < j; k++)
2730 Z(k,j), Ptmp, 1.0, Z(k,i));
2733 Telement::add(1.0, Ptmp, Z(j,i));
2736 for (
int k = 0; k < i; k++)
2737 Telement::gemm_upper_tr_only(
true,
false, 1.0,
2738 A(k,i), Z(k,i), 1.0, Ptmp);
2741 Telement::syInch(Ptmp, Z(i,i), threshold, looking, version);
2742 for (
int k = 0; k < i; k++) {
2749 assert(looking ==
right);
2751 Treal newThresh = 0;
2752 if (myThresh != 0 && Z.
ncols() > 1)
2753 newThresh = myThresh * 0.0001;
2755 newThresh = myThresh;
2757 for (
int i = 0; i < Z.
ncols(); i++) {
2760 for (
int k = 0; k < i; k++)
2761 Telement::gemm_upper_tr_only(
true,
false, 1.0,
2762 A(k,i), Z(k,i), 1.0, Ptmp);
2765 Telement::syInch(Ptmp, Z(i,i), newThresh, looking, version);
2766 for (
int k = 0; k < i; k++) {
2771 for (
int j = i + 1; j < Z.
ncols(); j++) {
2776 for (
int k = 0; k < i; k++)
2778 A(k,i), Z(k,j), 1.0, Ptmp);
2781 for (
int k = 0; k < i; k++)
2783 Z(k,i), Ptmp, 1.0, Z(k,j));
2786 Telement::add(1.0, Ptmp, Z(i,j));
2790 if (threshold != 0) {
2791 for (
int k = 0; k < i; k++)
2798 template<
class Treal,
class Telement>
2801 assert(!this->is_empty());
2802 if (this->nScalarsRows() == this->nScalarsCols()) {
2803 int size = this->nScalarsCols();
2804 Treal* colsums =
new Treal[size];
2805 Treal* diag =
new Treal[size];
2806 for (
int ind = 0; ind < size; ind++)
2808 this->add_abs_col_sums(colsums);
2809 this->get_diagonal(diag);
2812 lmin = diag[0] - tmp1;
2813 lmax = diag[0] + tmp1;
2814 for (
int col = 1; col < size; col++) {
2816 tmp2 = diag[col] - tmp1;
2817 lmin = (tmp2 < lmin ? tmp2 : lmin);
2818 tmp2 = diag[col] + tmp1;
2819 lmax = (tmp2 > lmax ? tmp2 : lmax);
2825 throw Failure(
"Matrix<Treal, Telement, int>::gershgorin(Treal, Treal): " 2826 "Matrix must be quadratic");
2830 template<
class Treal,
class Telement>
2834 if (!this->is_zero()) {
2836 for (
int col = 0; col < this->ncols(); col++) {
2837 for (
int row = 0; row < this->nrows(); row++) {
2838 (*this)(row,col).add_abs_col_sums(&sums[offset]);
2840 offset += (*this)(0,col).nScalarsCols();
2845 template<
class Treal,
class Telement>
2849 assert(this->nrows() == this->ncols());
2850 if (this->is_zero())
2851 for (
int rc = 0; rc < this->nScalarsCols(); rc++)
2855 for (
int rc = 0; rc < this->ncols(); rc++) {
2856 (*this)(rc,rc).get_diagonal(&diag[offset]);
2857 offset += (*this)(rc,rc).nScalarsCols();
2862 template<
class Treal,
class Telement>
2864 assert(!this->is_empty());
2866 if (this->is_zero())
2868 for (
int ind = 0; ind < this->nElements(); ind++)
2869 sum += this->elements[ind].memory_usage();
2873 template<
class Treal,
class Telement>
2876 if (!this->is_zero()) {
2877 for (
int ind = 0; ind < this->nElements(); ind++)
2878 sum += this->elements[ind].nnz();
2882 template<
class Treal,
class Telement>
2885 if (!this->is_zero()) {
2887 for (
int col = 1; col < this->ncols(); col++)
2888 for (
int row = 0; row < col; row++)
2889 sum += (*
this)(row, col).nnz();
2893 for (
int rc = 0; rc < this->nrows(); rc++)
2894 sum += (*
this)(rc,rc).sy_nnz();
2899 template<
class Treal,
class Telement>
2902 if (!this->is_zero()) {
2904 for (
int col = 1; col < this->ncols(); col++)
2905 for (
int row = 0; row < col; row++)
2906 sum += (*
this)(row, col).nvalues();
2908 for (
int rc = 0; rc < this->nrows(); rc++)
2909 sum += (*
this)(rc,rc).sy_nvalues();
2924 template<
class Treal>
2944 for (
int ind = 0; ind < this->
nElements(); ++ind)
2950 void fullMatrix(std::vector<Treal> & fullMat)
const;
2956 std::vector<int>
const & colind,
2957 std::vector<Treal>
const & values);
2959 std::vector<int>
const & colind,
2960 std::vector<Treal>
const & values,
2961 std::vector<int>
const & indexes);
2964 void addValues(std::vector<int>
const & rowind,
2965 std::vector<int>
const & colind,
2966 std::vector<Treal>
const & values);
2967 void addValues(std::vector<int>
const & rowind,
2968 std::vector<int>
const & colind,
2969 std::vector<Treal>
const & values,
2970 std::vector<int>
const & indexes);
2973 std::vector<int>
const & colind,
2974 std::vector<Treal>
const & values);
2977 std::vector<int>
const & colind,
2978 std::vector<Treal>
const & values);
2980 void getValues(std::vector<int>
const & rowind,
2981 std::vector<int>
const & colind,
2982 std::vector<Treal> & values)
const;
2983 void getValues(std::vector<int>
const & rowind,
2984 std::vector<int>
const & colind,
2985 std::vector<Treal> & values,
2986 std::vector<int>
const & indexes)
const;
2988 std::vector<int>
const & colind,
2989 std::vector<Treal> & values)
const;
2992 std::vector<int> & colind,
2993 std::vector<Treal> & values)
const;
2995 std::vector<int> & colind,
2996 std::vector<Treal> & values)
const;
3016 template<
typename TRule>
3018 template<
typename TRule>
3035 static void gemm(
const bool tA,
const bool tB,
const Treal alpha,
3040 static void symm(
const char side,
const char uplo,
const Treal alpha,
3045 static void syrk(
const char uplo,
const bool tA,
const Treal alpha,
3050 static void sysq(
const char uplo,
const Treal alpha,
3055 static void ssmm(
const Treal alpha,
3069 static void trmm(
const char side,
const char uplo,
const bool tA,
3100 Treal
trace()
const;
3108 static void add(
const Treal alpha,
3111 void assign(Treal
const alpha,
3128 (Treal
const threshold,
3141 for (
int ind = 0; ind < this->
nElements(); ind++)
3154 for (
int row = 0; row <
col; row++)
3155 (*
this)(row,
col) =
A(row,
col).frob();
3156 for (
int rc = 0; rc < this->
nrows(); rc++)
3157 (*
this)(rc,rc) =
A(rc,rc).syFrob();
3165 if (
A.is_zero() &&
B.is_zero() ) {
3173 if ( !A.
is_zero() && !
B.is_zero() ) {
3177 for (
int ind = 0; ind < this->
nElements(); ind++) {
3184 if ( !
A.is_zero() ) {
3189 if ( !
B.is_zero() ) {
3197 if (
A.is_zero() &&
B.is_zero() ) {
3205 if ( !A.
is_zero() && !
B.is_zero() ) {
3210 for (
int row = 0; row <
col; row++) {
3213 (*this)(row,
col) = Diff.
frob();
3215 for (
int rc = 0; rc < this->
ncols(); rc++) {
3218 (*this)(rc, rc) = Diff.
syFrob();
3222 if ( !
A.is_zero() ) {
3227 if ( !
B.is_zero() ) {
3241 for (
int ind = 0; ind < this->
nElements(); ind++)
3260 const bool tA,
const Treal alpha,
3280 const Treal threshold = 0,
3285 const Treal threshold = 0,
3288 inch(
A, Z, threshold, looking, version);
3291 void gershgorin(Treal& lmin, Treal& lmax)
const;
3306 return (
size_t)this->
nElements() *
sizeof(Treal);
3331 int n = this->
nrows();
3332 return this->
is_zero() ? 0 : n * (n + 1) / 2;
3345 for (
int row = 0; row <
col; row++) {
3346 res += 2*op.accumulate((*
this)(row,
col),
3350 res += op.accumulate((*
this)(
col,
col),
3364 for (
int row = 0; row < this->
nrows(); row++)
3365 res += op.accumulate((*
this)(row,
col),
3372 static inline unsigned int level() {
return 0;}
3378 Treal maxAbsGlobal = 0;
3379 Treal maxAbsLocal = 0;
3380 for (
int ind = 0; ind < this->
nElements(); ++ind) {
3382 maxAbsGlobal = maxAbsGlobal > maxAbsLocal ?
3383 maxAbsGlobal : maxAbsLocal;
3385 return maxAbsGlobal;
3397 std::cout<<
"operator="<<std::endl;
3410 void trim_memory_usage();
3412 void write_to_buffer_count(
int& zb_length,
int& vb_length)
const;
3413 void write_to_buffer(
int* zerobuf,
const int zb_length,
3414 Treal* valuebuf,
const int vb_length,
3415 int& zb_index,
int& vb_index)
const;
3416 void read_from_buffer(
int* zerobuf,
const int zb_length,
3417 Treal* valuebuf,
const int vb_length,
3418 int& zb_index,
int& vb_index);
3435 inline bool lowestLevel()
const {
return true;}
3445 template<
class Treal>
3447 template<
class Treal>
3452 template<
class Treal>
3457 assert((
int)fullMat.size() == nTotalRows * nTotalCols);
3463 for (
int row = 0; row < this->
nrows(); row++)
3465 fullMat[(colOffset +
col) * nTotalRows + (rowOffset + row)];
3468 template<
class Treal>
3473 fullMat.resize(nTotalRows * nTotalCols);
3479 fullMat[(colOffset +
col) * nTotalRows + (rowOffset + row)] = 0;
3483 for (
int row = 0; row < this->
nrows(); row++)
3484 fullMat[(colOffset +
col) * nTotalRows + (rowOffset + row)] =
3489 template<
class Treal>
3494 fullMat.resize(nTotalRows * nTotalCols);
3499 for (
int row = 0; row <=
col; row++) {
3500 fullMat[(colOffset +
col) * nTotalRows + (rowOffset + row)] = 0;
3501 fullMat[(rowOffset + row) * nTotalRows + (colOffset +
col)] = 0;
3506 for (
int row = 0; row <=
col; row++) {
3507 fullMat[(colOffset +
col) * nTotalRows + (rowOffset + row)] =
3509 fullMat[(rowOffset + row) * nTotalRows + (colOffset +
col)] =
3515 template<
class Treal>
3520 fullMat.resize(nTotalRows * nTotalCols);
3525 for (
int row = 0; row <= this->
nScalarsRows(); row++) {
3526 fullMat[(colOffset +
col) * nTotalRows + (rowOffset + row)] = 0;
3527 fullMat[(rowOffset + row) * nTotalRows + (colOffset +
col)] = 0;
3532 for (
int row = 0; row < this->
nrows(); row++) {
3533 fullMat[(colOffset +
col) * nTotalRows + (rowOffset + row)] =
3535 fullMat[(rowOffset + row) * nTotalRows + (colOffset +
col)] =
3543 template<
class Treal>
3546 std::vector<int>
const & colind,
3547 std::vector<Treal>
const & values) {
3548 assert(rowind.size() == colind.size() &&
3549 rowind.size() == values.size());
3550 std::vector<int> indexes(values.size());
3551 for (
int ind = 0; ind < values.size(); ++ind) {
3557 template<
class Treal>
3560 std::vector<int>
const & colind,
3561 std::vector<Treal>
const & values,
3562 std::vector<int>
const & indexes) {
3563 if (indexes.empty()) {
3569 for (
int ind = 0; ind < this->
nElements(); ++ind)
3571 std::vector<int>::const_iterator it;
3572 for ( it = indexes.begin(); it < indexes.end(); it++ ) {
3579 template<
class Treal>
3582 std::vector<int>
const & colind,
3583 std::vector<Treal>
const & values) {
3584 assert(rowind.size() == colind.size() &&
3585 rowind.size() == values.size());
3586 std::vector<int> indexes(values.size());
3587 for (
int ind = 0; ind < values.size(); ++ind) {
3590 addValues(rowind, colind, values, indexes);
3593 template<
class Treal>
3596 std::vector<int>
const & colind,
3597 std::vector<Treal>
const & values,
3598 std::vector<int>
const & indexes) {
3599 if (indexes.empty())
3603 std::vector<int>::const_iterator it;
3604 for ( it = indexes.begin(); it < indexes.end(); it++ ) {
3610 template<
class Treal>
3613 std::vector<int>
const & colind,
3614 std::vector<Treal>
const & values) {
3615 assert(rowind.size() == colind.size() &&
3616 rowind.size() == values.size());
3617 bool upperTriangleOnly =
true;
3618 for (
int ind = 0; ind < values.size(); ++ind) {
3620 upperTriangleOnly && colind[ind] >= rowind[ind];
3622 if (!upperTriangleOnly)
3623 throw Failure(
"Matrix<Treal>::" 3624 "syAddValues(std::vector<int> const &, " 3625 "std::vector<int> const &, " 3626 "std::vector<Treal> const &, int const): " 3627 "Only upper triangle can contain elements when assigning " 3628 "symmetric or triangular matrix ");
3632 template<
class Treal>
3635 std::vector<int>
const & colind,
3636 std::vector<Treal>
const & values) {
3637 assert(rowind.size() == colind.size() &&
3638 rowind.size() == values.size());
3639 bool upperTriangleOnly =
true;
3640 for (
int ind = 0; ind < values.size(); ++ind) {
3642 upperTriangleOnly && colind[ind] >= rowind[ind];
3644 if (!upperTriangleOnly)
3645 throw Failure(
"Matrix<Treal>::" 3646 "syAddValues(std::vector<int> const &, " 3647 "std::vector<int> const &, " 3648 "std::vector<Treal> const &, int const): " 3649 "Only upper triangle can contain elements when assigning " 3650 "symmetric or triangular matrix ");
3654 template<
class Treal>
3657 std::vector<int>
const & colind,
3658 std::vector<Treal> & values)
const {
3659 assert(rowind.size() == colind.size());
3660 values.resize(rowind.size(),0);
3661 std::vector<int> indexes(rowind.size());
3662 for (
int ind = 0; ind < rowind.size(); ++ind) {
3665 getValues(rowind, colind, values, indexes);
3668 template<
class Treal>
3671 std::vector<int>
const & colind,
3672 std::vector<Treal> & values,
3673 std::vector<int>
const & indexes)
const {
3675 if (indexes.empty())
3677 std::vector<int>::const_iterator it;
3679 for ( it = indexes.begin(); it < indexes.end(); it++ )
3683 for ( it = indexes.begin(); it < indexes.end(); it++ )
3689 template<
class Treal>
3692 std::vector<int>
const & colind,
3693 std::vector<Treal> & values)
const {
3694 assert(rowind.size() == colind.size());
3695 bool upperTriangleOnly =
true;
3696 for (
int ind = 0; ind < rowind.size(); ++ind) {
3698 upperTriangleOnly && colind[ind] >= rowind[ind];
3700 if (!upperTriangleOnly)
3701 throw Failure(
"Matrix<Treal>::" 3702 "syGetValues(std::vector<int> const &, " 3703 "std::vector<int> const &, " 3704 "std::vector<Treal> const &, int const): " 3705 "Only upper triangle when retrieving elements from " 3706 "symmetric or triangular matrix ");
3710 template<
class Treal>
3713 std::vector<int> & colind,
3714 std::vector<Treal> & values)
const {
3715 assert(rowind.size() == colind.size() &&
3716 rowind.size() == values.size());
3719 for (
int row = 0; row < this->
nrows(); row++) {
3722 values.push_back((*
this)(row,
col));
3726 template<
class Treal>
3729 std::vector<int> & colind,
3730 std::vector<Treal> & values)
const {
3731 assert(rowind.size() == colind.size() &&
3732 rowind.size() == values.size());
3735 for (
int row = 0; row <=
col; ++row) {
3738 values.push_back((*
this)(row,
col));
3743 template<
class Treal>
3749 template<
class Treal>
3755 char * tmp = (
char*)&ZERO;
3756 file.write(tmp,
sizeof(
int));
3759 char * tmp = (
char*)&ONE;
3760 file.write(tmp,
sizeof(
int));
3761 char * tmpel = (
char*)this->
elements;
3762 file.write(tmpel,
sizeof(Treal) * this->
nElements());
3766 template<
class Treal>
3771 char tmp[
sizeof(int)];
3772 file.read(tmp, (std::ifstream::pos_type)
sizeof(
int));
3773 switch ((
int)*tmp) {
3783 throw Failure(
"Matrix<Treal>::" 3784 "readFromFile(std::ifstream & file):" 3785 "File corruption, int value not 0 or 1");
3789 template<
class Treal>
3793 for (
int ind = 0; ind < this->
nElements(); ind++)
3794 this->
elements[ind] = rand() / (Treal)RAND_MAX;
3797 template<
class Treal>
3803 for (
int row = 0; row <
col; row++)
3804 (*
this)(row,
col) = rand() / (Treal)RAND_MAX;;
3806 for (
int rc = 0; rc < this->
nrows(); rc++)
3807 (*
this)(rc,rc) = rand() / (Treal)RAND_MAX;;
3810 template<
class Treal>
3814 probabilityBeingZero > rand() / (Treal)RAND_MAX)
3820 template<
class Treal>
3824 probabilityBeingZero > rand() / (Treal)RAND_MAX)
3830 template<
class Treal>
3831 template<
typename TRule>
3837 for (
int row = 0; row < this->
nrows(); row++)
3838 (*
this)(row,
col) = rule.set(this->rows.getOffset() + row,
3841 template<
class Treal>
3842 template<
typename TRule>
3849 for (
int row = 0; row <=
col; row++)
3850 (*
this)(row,
col) = rule.set(this->rows.getOffset() + row,
3855 template<
class Treal>
3859 throw Failure(
"Matrix<Treal>::addIdentity(Treal): " 3860 "Cannot add identity to empty matrix.");
3862 throw Failure(
"Matrix<Treal, Telement>::addIdentity(Treal): " 3863 "Matrix must be square to add identity");
3866 for (
int ind = 0; ind < this->
ncols(); ind++)
3867 (*
this)(ind,ind) += alpha;
3870 template<
class Treal>
3882 AT.
elements = allocateElements<Treal>(
A.nElements());
3886 for (
int row = 0; row < AT.
nrows(); row++)
3891 template<
class Treal>
3898 for (
int row = 1; row < this->
nrows(); row++)
3900 (*
this)(row,
col) = (*
this)(
col, row);
3904 throw Failure(
"Matrix<Treal>::symToNosym(): " 3905 "Only quadratic matrices can be symmetric");
3908 template<
class Treal>
3916 for (
int row =
col + 1; row < this->
nrows(); row++)
3917 (*
this)(row,
col) = 0;
3921 throw Failure(
"Matrix<Treal>::nosymToSym(): " 3922 "Only quadratic matrices can be symmetric");
3925 template<
class Treal>
3934 throw Failure(
"Matrix<Treal>::operator=(int k = 1): " 3935 "Matrix must be quadratic to " 3936 "become identity matrix.");
3939 for (
int rc = 0; rc < this->
ncols(); rc++)
3944 (
"Matrix<Treal>::operator=(int k) only implemented for k = 0, k = 1");
3949 template<
class Treal>
3952 if (!this->
is_zero() && alpha != 1) {
3953 for (
int ind = 0; ind < this->
nElements(); ind++)
3959 template<
class Treal>
3961 gemm(
const bool tA,
const bool tB,
const Treal alpha,
3980 if ((
A.is_zero() ||
B.is_zero() || alpha == 0) &&
3984 Treal beta_tmp = beta;
3990 if (!
A.is_zero() && !
B.is_zero() && alpha != 0) {
3992 if (
A.ncols() ==
B.nrows() &&
3993 A.nrows() == C.
nrows() &&
3995 mat::gemm(
"N",
"N", &
A.nrows(), &
B.ncols(), &
A.ncols(), &alpha,
3996 A.elements, &
A.nrows(),
B.elements, &
B.nrows(),
3999 throw Failure(
"Matrix<Treal>::" 4000 "gemm(bool, bool, Treal, Matrix<Treal>, " 4001 "Matrix<Treal>, Treal, " 4003 "Incorrect matrixdimensions for multiplication");
4005 if (
A.nrows() ==
B.nrows() &&
4006 A.ncols() == C.
nrows() &&
4008 mat::gemm(
"T",
"N", &
A.ncols(), &
B.ncols(), &
A.nrows(), &alpha,
4009 A.elements, &
A.nrows(),
B.elements, &
B.nrows(),
4012 throw Failure(
"Matrix<Treal>::" 4013 "gemm(bool, bool, Treal, Matrix<Treal>, " 4014 "Matrix<Treal>, Treal, " 4016 "Incorrect matrixdimensions for multiplication");
4018 if (
A.ncols() ==
B.ncols() &&
4019 A.nrows() == C.
nrows() &&
4021 mat::gemm(
"N",
"T", &
A.nrows(), &
B.nrows(), &
A.ncols(), &alpha,
4022 A.elements, &
A.nrows(),
B.elements, &
B.nrows(),
4025 throw Failure(
"Matrix<Treal>::" 4026 "gemm(bool, bool, Treal, Matrix<Treal>, " 4027 "Matrix<Treal>, Treal, " 4029 "Incorrect matrixdimensions for multiplication");
4031 if (
A.nrows() ==
B.ncols() &&
4032 A.ncols() == C.
nrows() &&
4034 mat::gemm(
"T",
"T", &
A.ncols(), &
B.nrows(), &
A.nrows(), &alpha,
4035 A.elements, &
A.nrows(),
B.elements, &
B.nrows(),
4038 throw Failure(
"Matrix<Treal>::" 4039 "gemm(bool, bool, Treal, Matrix<Treal>, " 4040 "Matrix<Treal>, Treal, " 4042 "Incorrect matrixdimensions for multiplication");
4043 else throw Failure(
"Matrix<Treal>::" 4044 "gemm(bool, bool, Treal, Matrix<Treal>, " 4045 "Matrix<Treal>, Treal, " 4046 "Matrix<Treal>):Very strange error!!");
4054 template<
class Treal>
4056 symm(
const char side,
const char uplo,
const Treal alpha,
4078 if ((
A.is_zero() ||
B.is_zero() || alpha == 0) &&
4082 Treal beta_tmp = beta;
4087 if (!
A.is_zero() && !
B.is_zero() && alpha != 0) {
4088 if (
A.nrows() ==
A.ncols() && C.
nrows() ==
B.nrows() &&
4089 C.
ncols() ==
B.ncols() &&
4090 ((
side ==
'L' &&
A.ncols() ==
B.nrows()) ||
4091 (
side ==
'R' &&
B.ncols() ==
A.nrows())))
4093 A.elements, &
A.nrows(),
B.elements, &
B.nrows(),
4096 throw Failure(
"Matrix<Treal>::symm: Incorrect matrix " 4097 "dimensions for symmetric multiply");
4104 template<
class Treal>
4106 syrk(
const char uplo,
const bool tA,
const Treal alpha,
4123 ((!tA &&
A.nrows() == C.
nrows()) || (tA &&
A.ncols() == C.
nrows())))
4124 if (alpha != 0 && !
A.is_zero()) {
4125 Treal beta_tmp = beta;
4132 A.elements, &
A.nrows(), &beta_tmp,
4137 A.elements, &
A.nrows(), &beta_tmp,
4143 throw Failure(
"Matrix<Treal>::syrk: Incorrect matrix " 4144 "dimensions for symmetric rank-k update");
4147 template<
class Treal>
4149 sysq(
const char uplo,
const Treal alpha,
4166 template<
class Treal>
4189 template<
class Treal>
4197 ssmm(alpha,
A,
B, beta, C);
4202 template<
class Treal>
4210 if (alpha != 0 && !
A.is_zero() && !
B.is_zero())
4211 if (((
side ==
'R' &&
B.ncols() ==
A.nrows()) ||
4212 (
side ==
'L' &&
A.ncols() ==
B.nrows())) &&
4213 A.nrows() ==
A.ncols())
4216 &alpha,
A.elements, &
A.nrows(),
B.elements, &
B.nrows());
4219 &alpha,
A.elements, &
A.nrows(),
B.elements, &
B.nrows());
4221 throw Failure(
"Matrix<Treal>::trmm: " 4222 "Incorrect matrix dimensions for multiplication");
4227 template<
class Treal>
4234 for (
int i = 0; i < this->
nElements(); i++)
4240 template<
class Treal>
4245 throw Failure(
"Matrix<Treal>::syFrobSquared(): " 4246 "Matrix must be have equal number of rows " 4247 "and cols to be symmetric");
4251 for (
int row = 0; row <
col; row++)
4252 sum += 2 * (*
this)(row,
col) * (*
this)(row,
col);
4253 for (
int rc = 0; rc < this->
ncols(); rc++)
4254 sum += (*
this)(rc, rc) * (*
this)(rc, rc);
4259 template<
class Treal>
4264 assert(!
A.is_empty());
4265 assert(!
B.is_empty());
4266 if (
A.nrows() !=
B.nrows() ||
A.ncols() !=
B.ncols())
4267 throw Failure(
"Matrix<Treal>::frobSquaredDiff: " 4268 "Incorrect matrix dimensions");
4270 if (!
A.is_zero() && !
B.is_zero()) {
4272 for (
int i = 0; i <
A.nElements(); i++) {
4273 diff =
A.elements[i] -
B.elements[i];
4277 else if (
A.is_zero() && !
B.is_zero())
4278 sum =
B.frobSquared();
4279 else if (!
A.is_zero() &&
B.is_zero())
4280 sum =
A.frobSquared();
4286 template<
class Treal>
4291 assert(!
A.is_empty());
4292 assert(!
B.is_empty());
4293 if (
A.nrows() !=
B.nrows() ||
4294 A.ncols() !=
B.ncols() ||
4295 A.nrows() !=
A.ncols())
4296 throw Failure(
"Matrix<Treal>::syFrobSquaredDiff: " 4297 "Incorrect matrix dimensions");
4299 if (!
A.is_zero() && !
B.is_zero()) {
4301 for (
int col = 1; col <
A.ncols(); col++)
4302 for (
int row = 0; row < col; row++) {
4303 diff =
A(row, col) -
B(row, col);
4304 sum += 2 * diff * diff;
4306 for (
int rc = 0; rc <
A.ncols(); rc++) {
4307 diff =
A(rc, rc) -
B(rc, rc);
4311 else if (
A.is_zero() && !
B.is_zero())
4312 sum =
B.syFrobSquared();
4313 else if (!
A.is_zero() &&
B.is_zero())
4314 sum =
A.syFrobSquared();
4319 template<
class Treal>
4324 throw Failure(
"Matrix<Treal>::trace(): Matrix must be quadratic");
4329 for (
int rc = 0; rc < this->
ncols(); rc++)
4330 sum += (*
this)(rc,rc);
4335 template<
class Treal>
4341 if (
A.nrows() !=
B.ncols() ||
A.ncols() !=
B.nrows())
4342 throw Failure(
"Matrix<Treal>::trace_ab: " 4343 "Wrong matrix dimensions for trace(A * B)");
4345 if (!
A.is_zero() && !
B.is_zero())
4346 for (
int rc = 0; rc <
A.nrows(); rc++)
4347 for (
int colA = 0; colA <
A.ncols(); colA++)
4348 tr +=
A(rc,colA) *
B(colA, rc);
4352 template<
class Treal>
4358 if (
A.ncols() !=
B.ncols() ||
A.nrows() !=
B.nrows())
4359 throw Failure(
"Matrix<Treal>::trace_aTb: " 4360 "Wrong matrix dimensions for trace(A' * B)");
4362 if (!
A.is_zero() && !
B.is_zero())
4363 for (
int rc = 0; rc <
A.ncols(); rc++)
4364 for (
int rowB = 0; rowB <
B.nrows(); rowB++)
4365 tr +=
A(rowB,rc) *
B(rowB, rc);
4369 template<
class Treal>
4375 if (
A.nrows() !=
B.ncols() ||
A.ncols() !=
B.nrows() ||
4376 A.nrows() !=
A.ncols())
4377 throw Failure(
"Matrix<Treal>::sy_trace_ab: " 4378 "Wrong matrix dimensions for symmetric trace(A * B)");
4379 if (
A.is_zero() ||
B.is_zero())
4384 for (
int rc = 0; rc <
A.nrows(); rc++)
4385 tr +=
A(rc,rc) *
B(rc, rc);
4387 for (
int colA = 1; colA <
A.ncols(); colA++)
4388 for (
int rowA = 0; rowA < colA; rowA++)
4389 tr += 2 *
A(rowA, colA) *
B(rowA, colA);
4393 template<
class Treal>
4400 if (
A.nrows() !=
B.nrows() ||
A.ncols() !=
B.ncols())
4401 throw Failure(
"Matrix<Treal>::add: " 4402 "Wrong matrix dimensions for addition");
4403 if (!
A.is_zero() && alpha != 0) {
4406 for (
int ind = 0; ind <
A.nElements(); ind++)
4407 B.elements[ind] += alpha *
A.elements[ind];
4411 template<
class Treal>
4426 template<
class Treal>
4433 template<
class Treal>
4437 for (
int ind = 0; ind < this->
nElements(); ind++)
4443 template<
class Treal>
4448 if ((!this->is_zero() && this->frobSquared() <= threshold) ||
4453 else if (!this->is_zero() && this->frobSquared() <= threshold)
4457 template<
class Treal>
4461 assert(!this->is_empty());
4462 bool thisMatIsZero =
true;
4465 bool EMatIsZero =
true;
4466 if (!ErrorMatrix->
is_zero() || !this->is_zero()) {
4469 if (this->is_zero())
4471 for (
int ind = 0; ind < this->nElements(); ind++) {
4472 if ( this->elements[ind] != 0 ) {
4473 assert(ErrorMatrix->
elements[ind] == 0);
4475 if ( this->elements[ind] * this->elements[ind] <= threshold ) {
4477 ErrorMatrix->
elements[ind] = this->elements[ind];
4478 this->elements[ind] = 0;
4482 thisMatIsZero =
false;
4485 if ( ErrorMatrix->
elements[ind] != 0 ) {
4486 assert(this->elements[ind] == 0);
4488 if ( ErrorMatrix->
elements[ind] * ErrorMatrix->
elements[ind] > threshold ) {
4490 this->elements[ind] = ErrorMatrix->
elements[ind];
4492 thisMatIsZero =
false;
4498 if (thisMatIsZero) {
4500 for (
int ind = 0; ind < this->nElements(); ind++)
4501 assert( this->elements[ind] == 0);
4507 for (
int ind = 0; ind < this->nElements(); ind++)
4508 assert( ErrorMatrix->
elements[ind] == 0);
4510 ErrorMatrix->
clear();
4515 if (!this->is_zero()) {
4516 for (
int ind = 0; ind < this->nElements(); ind++) {
4517 if ( this->elements[ind] * this->elements[ind] <= threshold )
4520 this->elements[ind] = 0;
4522 thisMatIsZero =
false;
4535 template<
class Treal>
4553 template<
class Treal>
4566 template<
class Treal>
4569 const bool tA,
const Treal alpha,
4583 template<
class Treal>
4604 template<
class Treal>
4607 assert(!this->is_empty());
4608 if (ErrorMatrix && ErrorMatrix->
is_empty()) {
4612 Treal fs = frobSquared();
4613 if (fs < threshold) {
4623 template<
class Treal>
4627 const Treal threshold,
const side looking,
4630 if (
A.nrows() !=
A.ncols())
4631 throw Failure(
"Matrix<Treal>::inch: Matrix must be quadratic!");
4633 throw Failure(
"Matrix<Treal>::inch: Matrix is zero! " 4634 "Not possible to compute inverse cholesky.");
4639 throw Failure(
"Matrix<Treal>::inch: potrf returned info > 0. The matrix is not positive definite.");
4641 throw Failure(
"Matrix<Treal>::inch: potrf returned info < 0");
4645 throw Failure(
"Matrix<Treal>::inch: trtri returned info > 0. The matrix is not positive definite.");
4647 throw Failure(
"Matrix<Treal>::inch: trtri returned info < 0");
4652 template<
class Treal>
4658 Treal* colsums =
new Treal[size];
4659 Treal* diag =
new Treal[size];
4660 for (
int ind = 0; ind < size; ind++)
4666 lmin = diag[0] - tmp1;
4667 lmax = diag[0] + tmp1;
4670 tmp2 = diag[
col] - tmp1;
4671 lmin = (tmp2 < lmin ? tmp2 : lmin);
4672 tmp2 = diag[
col] + tmp1;
4673 lmax = (tmp2 > lmax ? tmp2 : lmax);
4679 throw Failure(
"Matrix<Treal>::gershgorin(Treal, Treal): Matrix must be quadratic");
4683 template<
class Treal>
4689 for (
int row = 0; row < this->
nrows(); row++)
4693 template<
class Treal>
4702 for (
int rc = 0; rc < this->
ncols(); rc++)
4703 diag[rc] = (*
this)(rc,rc);
4730 template<
class Treal>
4732 if (this->
is_zero() && this->cap > 0) {
4738 else if (this->cap > this->nel) {
4739 Treal* tmp =
new Treal[this->nel];
4740 for (
int i = 0; i < this->nel; i++)
4743 this->cap = this->nel;
4746 assert(this->cap == this->nel);
4753 template<
class Treal>
4754 void Matrix<Treal>::
4755 write_to_buffer_count(
int& zb_length,
int& vb_length)
const {
4761 vb_length += this->nel;
4765 template<
class Treal>
4766 void Matrix<Treal>::
4767 write_to_buffer(
int* zerobuf,
const int zb_length,
4768 Treal* valuebuf,
const int vb_length,
4769 int& zb_index,
int& vb_index)
const {
4770 if (zb_index < zb_length) {
4772 zerobuf[zb_index] = 0;
4776 if (vb_index + this->nel < vb_length + 1) {
4777 zerobuf[zb_index] = 1;
4779 for (
int i = 0; i < this->nel; i++)
4780 valuebuf[vb_index + i] = this->
elements[i];
4781 vb_index += this->nel;
4784 throw Failure(
"Matrix<Treal, Telement>::write_to_buffer: " 4785 "Insufficient space in buffers");
4789 throw Failure(
"Matrix<Treal, Telement>::write_to_buffer: " 4790 "Insufficient space in buffers");
4793 template<
class Treal>
4794 void Matrix<Treal>::
4795 read_from_buffer(
int* zerobuf,
const int zb_length,
4796 Treal* valuebuf,
const int vb_length,
4797 int& zb_index,
int& vb_index) {
4798 if (zb_index < zb_length) {
4799 if (zerobuf[zb_index] == 0) {
4804 this->content =
ful;
4806 this->assert_alloc();
4807 if (vb_index + this->nel < vb_length + 1) {
4808 assert(zerobuf[zb_index] == 1);
4810 for (
int i = 0; i < this->nel; i++)
4811 this->
elements[i] = valuebuf[vb_index + i];
4812 vb_index += this->nel;
4815 throw Failure(
"Matrix<Treal, Telement>::read_from_buffer: " 4816 "Mismatch, buffers does not match matrix");
4820 throw Failure(
"Matrix<Treal, Telement>::read_from_buffer: " 4821 "Mismatch, buffers does not match matrix");
int getOffset() const
Definition: SizesAndBlocks.h:83
Vector< Treal, typename ElementType::VectorType > VectorType
Definition: Matrix.h:118
void syUpTriFullMatrix(std::vector< Treal > &fullMat) const
Definition: Matrix.h:563
static Treal sy_trace_ab(const Matrix< Treal, Telement > &A, const Matrix< Treal, Telement > &B)
Definition: Matrix.h:2002
static void trifulltofull(Treal *trifull, const int size)
Definition: mat_gblas.h:1193
void syAssignFrobNormsLowestLevel(Matrix< Treal, Matrix< Treal > > const &A)
Definition: Matrix.h:3148
Definition: Matrix.h:2925
void assign(Treal const alpha, Matrix< Treal, Telement > const &A)
Definition: Matrix.h:2043
void frobThreshElementLevel(Treal const threshold, Matrix< Treal, Telement > *ErrorMatrix)
Definition: Matrix.h:2115
Treal frob_squared_thresh(Treal const threshold, Matrix< Treal, Telement > *ErrorMatrix=0)
Removes small elements so that the introduced error is smaller than threshold in the squared Frobeniu...
Definition: Matrix.h:2640
mat::SizesAndBlocks rows
Definition: test.cc:51
void fullMatrix(std::vector< Treal > &fullMat) const
Definition: Matrix.h:519
size_t sy_nnz() const
Returns number of nonzeros in matrix including lower triangle elements.
Definition: Matrix.h:2883
void gershgorin(Treal &lmin, Treal &lmax) const
Definition: Matrix.h:4654
Vector< Treal, Treal > VectorType
Definition: Matrix.h:2928
SizesAndBlocks cols
Definition: MatrixHierarchicBase.h:165
mat::SizesAndBlocks cols
Definition: test.cc:52
double accumulatedTime
Definition: Matrix.h:84
void assignDiffFrobNormsLowestLevel(Matrix< Treal, Matrix< Treal, Telement > > const &A, Matrix< Treal, Matrix< Treal, Telement > > const &B)
Same as assignFrobNormsLowestLevel except that the Frobenius norms of the differences between submatr...
Definition: Matrix.h:2189
#define MAT_OMP_FINALIZE
Definition: matInclude.h:92
void writeToFile(std::ofstream &file) const
Definition: Matrix.h:845
inchversion
Definition: Matrix.h:76
static void syInch(const Matrix< Treal, Telement > &A, Matrix< Treal, Telement > &Z, const Treal threshold=0, const side looking=left, const inchversion version=unstable)
Definition: Matrix.h:2690
Treal syFrobSquared() const
Definition: Matrix.h:1882
const int & nScalarsRows() const
Definition: MatrixHierarchicBase.h:63
void symToNosym()
Definition: Matrix.h:3893
void tic()
Definition: matInclude.cc:84
int const & getNBlocks() const
Definition: SizesAndBlocks.h:72
static void add(const Treal alpha, const Matrix< Treal, Telement > &A, Matrix< Treal, Telement > &B)
Definition: Matrix.h:2025
SizesAndBlocks getSizesAndBlocksForLowerLevel(int const blockNumber) const
Definition: SizesAndBlocks.cc:71
Proxy structs used by the matrix API.
Treal frobSquared() const
Definition: Matrix.h:4228
int whichBlock(int const globalIndex) const
Returns the blocknumber (between 0 and nBlocks-1) that contains elements with the given global index...
Definition: SizesAndBlocks.h:79
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
bool highestLevel() const
Definition: MatrixHierarchicBase.h:133
static void transpose(Matrix< Treal, Telement > const &A, Matrix< Treal, Telement > &AT)
Definition: Matrix.h:979
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
Matrix()
Definition: Matrix.h:2931
void setElementsByRule(TRule &rule)
Definition: Matrix.h:940
Vector class.
Definition: Matrix.h:78
MatrixHierarchicBase< Treal, Telement > & operator=(const MatrixHierarchicBase< Treal, Telement > &mat)
Definition: MatrixHierarchicBase.h:202
Treal frob() const
Definition: Matrix.h:286
return elements[row+col *nrows()]
Definition: MatrixHierarchicBase.h:81
void add_abs_col_sums(Treal *abscolsums) const
Definition: Matrix.h:2832
void syRandomZeroStructure(Treal probabilityBeingZero)
Definition: Matrix.h:920
static void ssmm(const Treal alpha, const Matrix< Treal, Telement > &A, const Matrix< Treal, Telement > &B, const Treal beta, Matrix< Treal, Telement > &C)
Definition: Matrix.h:1562
void reset()
Definition: Matrix.h:86
void getValues(std::vector< int > const &rowind, std::vector< int > const &colind, std::vector< Treal > &values) const
Definition: Matrix.h:733
Treal frob_thresh(Treal const threshold, Matrix< Treal > *ErrorMatrix=0)
Definition: Matrix.h:3267
Treal frobSquared() const
Definition: Matrix.h:1868
void gershgorin(Treal &lmin, Treal &lmax) const
Definition: Matrix.h:2800
Matrix< Treal, Telement > & operator=(const Matrix< Treal, Telement > &mat)
Definition: Matrix.h:198
void syAssignDiffFrobNormsLowestLevel(Matrix< Treal, Matrix< Treal > > const &A, Matrix< Treal, Matrix< Treal > > const &B)
Definition: Matrix.h:3195
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:334
size_t nnz() const
Returns number of nonzeros in matrix.
Definition: Matrix.h:2874
Definition: allocate.cc:39
Describes dimensions of matrix and its blocks on all levels.
Definition: SizesAndBlocks.h:45
int getNTotalScalars() const
Definition: SizesAndBlocks.h:84
const int & ncols() const
Definition: MatrixHierarchicBase.h:71
void syGetAllValues(std::vector< int > &rowind, std::vector< int > &colind, std::vector< Treal > &) const
Definition: Matrix.h:820
void addValues(std::vector< int > const &rowind, std::vector< int > const &colind, std::vector< Treal > const &values)
Definition: Matrix.h:640
void truncateAccordingToSparsityPattern(Matrix< Treal, Matrix< Treal > > &A) const
Definition: Matrix.h:3235
static Treal syFrobSquaredDiff(const Matrix< Treal, Telement > &A, const Matrix< Treal, Telement > &B)
Definition: Matrix.h:1924
void truncateAccordingToSparsityPattern(Matrix< Treal, Matrix< Treal, Telement > > &A) const
Truncate matrix A according to the sparsity pattern of the this matrix (frobNormMat).
Definition: Matrix.h:2257
void syAddValues(std::vector< int > const &rowind, std::vector< int > const &colind, std::vector< Treal > const &values)
Definition: Matrix.h:711
void getFrobSqLowestLevel(std::vector< Treal > &frobsq) const
Definition: Matrix.h:2061
A quicksort implementation.
Code for memory allocation/deallocation routines used by matrix library.
void clear()
Set matrix to zero and delete all arrays.
Definition: Matrix.h:3744
Definition: matInclude.h:156
Treal template_blas_fabs(Treal x)
float toc()
Definition: matInclude.cc:88
Treal frob() const
Definition: Matrix.h:3075
Definition: matInclude.h:138
void allocate()
Definition: Matrix.h:124
void symToNosym()
Definition: Matrix.h:1001
static void trmm_upper_tr_only(const char side, const char uplo, const bool tA, const Treal alpha, const Matrix< Treal, Telement > &A, Matrix< Treal, Telement > &B)
Definition: Matrix.h:2543
Matrix< Treal, Telement > & operator*=(const Treal alpha)
Definition: Matrix.h:1073
size_t memory_usage() const
Definition: Matrix.h:2863
static Treal trace_ab(const Matrix< Treal, Telement > &A, const Matrix< Treal, Telement > &B)
Definition: Matrix.h:1968
Telement int col
Definition: MatrixHierarchicBase.h:75
static void gemm_upper_tr_only(const bool tA, const bool tB, const Treal alpha, const Matrix< Treal, Telement > &A, const Matrix< Treal, Telement > &B, const Treal beta, Matrix< Treal, Telement > &C)
Definition: Matrix.h:2277
Base class for Matrix and Matrix specialization.
Definition: MatrixHierarchicBase.h:52
static void trtri(const char *uplo, const char *diag, const int *n, T *A, const int *lda, int *info)
Definition: mat_gblas.h:321
size_t nvalues() const
Returns number of stored values in matrix.
Definition: Matrix.h:444
void addTime(double timeToAdd)
Definition: Matrix.h:88
static Treal trace_aTb(const Matrix< Treal, Telement > &A, const Matrix< Treal, Telement > &B)
Definition: Matrix.h:1985
C++ interface to a subset of BLAS and LAPACK.
Treal syFrob() const
Definition: Matrix.h:3077
void syAssignFrobNormsLowestLevel(Matrix< Treal, Matrix< Treal, Telement > > const &A)
Version of assignFrobNormsLowestLevelToMatrix for symmetric matrices.
Definition: Matrix.h:2167
bool is_zero() const
Definition: MatrixHierarchicBase.h:108
size_t nvalues() const
Returns number of stored values in matrix.
Definition: Matrix.h:3324
void random()
Definition: Matrix.h:884
int nElements() const
Definition: MatrixHierarchicBase.h:110
static unsigned int level()
Definition: Matrix.h:3372
static Treal syFrobDiff(const Matrix< Treal, Telement > &A, const Matrix< Treal, Telement > &B)
Definition: Matrix.h:306
void syRandom()
Definition: Matrix.h:892
size_t sy_nvalues() const
Returns number of stored values in matrix.
Definition: Matrix.h:3329
Treal geAccumulateWith(Top &op)
Definition: Matrix.h:3358
static Treal frobDiff(const Matrix< Treal, Telement > &A, const Matrix< Treal, Telement > &B)
Definition: Matrix.h:297
static void potrf(const char *uplo, const int *n, T *A, const int *lda, int *info)
Definition: mat_gblas.h:314
void allocate()
Definition: Matrix.h:2940
static const Treal ONE
Definition: Matrix.h:3442
void get_diagonal(Treal *diag) const
Definition: Matrix.h:2847
const int & nrows() const
Definition: MatrixHierarchicBase.h:69
Treal maxAbsValue() const
Definition: Matrix.h:3374
void sy_gershgorin(Treal &lmin, Treal &lmax) const
Definition: Matrix.h:3292
void freeElements(float *ptr)
Definition: allocate.cc:49
Treal maxAbsValue() const
Definition: Matrix.h:482
static void ssmm_upper_tr_only(const Treal alpha, const Matrix< Treal, Telement > &A, const Matrix< Treal, Telement > &B, const Treal beta, Matrix< Treal, Telement > &C)
Definition: Matrix.h:1691
void syFullMatrix(std::vector< Treal > &fullMat) const
Definition: Matrix.h:539
#define MAT_OMP_START
Definition: matInclude.h:90
static unsigned int level()
Definition: Matrix.h:480
Telement ElementType
Definition: Matrix.h:117
Matrix class and heart of the matrix library.
Definition: Matrix.h:115
Matrix< Treal > & operator=(const Matrix< Treal > &mat)
Definition: Matrix.h:2999
void randomZeroStructure(Treal probabilityBeingZero)
Get a random zero structure with a specified probability that each submatrix is zero.
Definition: Matrix.h:906
double getAccumulatedTime()
Definition: Matrix.h:87
void getAllValues(std::vector< int > &rowind, std::vector< int > &colind, std::vector< Treal > &) const
Definition: Matrix.h:808
static Treal frobSquaredDiff(const Matrix< Treal, Telement > &A, const Matrix< Treal, Telement > &B)
Definition: Matrix.h:1902
void frobThreshLowestLevel(Treal const threshold, Matrix< Treal, Telement > *ErrorMatrix)
Definition: Matrix.h:2070
void assignFromSparse(std::vector< int > const &rowind, std::vector< int > const &colind, std::vector< Treal > const &values)
Definition: Matrix.h:588
#define MAT_OMP_END
Definition: matInclude.h:91
static void syInch(const Matrix< Treal > &A, Matrix< Treal > &Z, const Treal threshold=0, const side looking=left, const inchversion version=unstable)
Definition: Matrix.h:3283
~Matrix()
Definition: Matrix.h:205
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
const int & nScalarsCols() const
Definition: MatrixHierarchicBase.h:65
void assignFromFull(std::vector< Treal > const &fullMat)
Definition: Matrix.h:506
size_t memory_usage() const
Definition: Matrix.h:3301
void getFrobSqElementLevel(std::vector< Treal > &frobsq) const
Definition: Matrix.h:2106
SingletonForTimings()
Definition: Matrix.h:102
void syAssignFromSparse(std::vector< int > const &rowind, std::vector< int > const &colind, std::vector< Treal > const &values)
Definition: Matrix.h:689
void sySetElementsByRule(TRule &rule)
Definition: Matrix.h:949
size_t sy_nvalues() const
Returns number of stored values in matrix.
Definition: Matrix.h:2900
static void trsytriplemm(char const side, const Matrix< Treal, Telement > &Z, Matrix< Treal, Telement > &A)
Definition: Matrix.h:2620
static void syrk(const char uplo, const bool tA, const Treal alpha, const Matrix< Treal, Telement > &A, const Treal beta, Matrix< Treal, Telement > &C)
Definition: Matrix.h:1370
void sy_gershgorin(Treal &lmin, Treal &lmax) const
Definition: Matrix.h:422
side
Definition: Matrix.h:75
static SingletonForTimings & instance()
Definition: Matrix.h:96
void resetCols(SizesAndBlocks const &newCols)
Definition: MatrixHierarchicBase.h:119
Treal syFrob() const
Definition: Matrix.h:291
bool is_empty() const
Check if matrix is empty Empty is different from zero, a zero matrix contains information about block...
Definition: MatrixHierarchicBase.h:143
static void trmm(const char side, const char uplo, const bool tA, const Treal alpha, const Matrix< Treal, Telement > &A, Matrix< Treal, Telement > &B)
Definition: Matrix.h:1788
void assignDiffFrobNormsLowestLevel(Matrix< Treal, Matrix< Treal > > const &A, Matrix< Treal, Matrix< Treal > > const &B)
Definition: Matrix.h:3163
void syAssignDiffFrobNormsLowestLevel(Matrix< Treal, Matrix< Treal, Telement > > const &A, Matrix< Treal, Matrix< Treal, Telement > > const &B)
Same as syAssignFrobNormsLowestLevel except that the Frobenius norms of the differences between subma...
Definition: Matrix.h:2220
~Matrix()
Definition: Matrix.h:3005
void assignFrobNormsLowestLevel(Matrix< Treal, Matrix< Treal, Telement > > const &A)
Build a matrix with single matrix elements at the lowest level containing the Frobenius norms of the ...
Definition: Matrix.h:2153
void trSetElementsByRule(TRule &rule)
Definition: Matrix.h:224
Treal geAccumulateWith(Top &op)
Accumulation algorithm for general matrices.
Definition: Matrix.h:470
#define MAT_OMP_INIT
Definition: matInclude.h:89
size_t sy_nnz() const
Returns number of nonzeros in matrix including lower triangle elements.
Definition: Matrix.h:3315
static void symm(const char side, const char uplo, const Treal alpha, const Matrix< Treal, Telement > &A, const Matrix< Treal, Telement > &B, const Treal beta, Matrix< Treal, Telement > &C)
Definition: Matrix.h:1253
void readFromFile(std::ifstream &file)
Definition: Matrix.h:861
void nosymToSym()
Definition: Matrix.h:3910
Matrix()
Definition: Matrix.h:121
size_t nnz() const
Returns number of nonzeros in matrix.
Definition: Matrix.h:3309
static void sysq(const char uplo, const Treal alpha, const Matrix< Treal, Telement > &A, const Treal beta, Matrix< Treal, Telement > &C)
Definition: Matrix.h:1477
static const Treal ZERO
Definition: Matrix.h:3441
SizesAndBlocks rows
Definition: MatrixHierarchicBase.h:164
Treal syAccumulateWith(Top &op)
Definition: Matrix.h:3339
Xtrans< TX > transpose(TX const &A)
Transposition.
Definition: matrix_proxy.h:131
void syGetValues(std::vector< int > const &rowind, std::vector< int > const &colind, std::vector< Treal > &values) const
Definition: Matrix.h:786
Treal trace(const XYZ< Treal, MatrixGeneral< Treal, Tmatrix >, MatrixGeneral< Treal, Tmatrix > > &smm)
Definition: MatrixGeneral.h:904
void resetRows(SizesAndBlocks const &newRows)
Definition: MatrixHierarchicBase.h:114
void addIdentity(Treal alpha)
Definition: Matrix.h:964
static void sytr_upper_tr_only(char const side, const Treal alpha, Matrix< Treal, Telement > &A, const Matrix< Treal, Telement > &Z)
Definition: Matrix.h:2460
Treal template_blas_sqrt(Treal x)
void nosymToSym()
Definition: Matrix.h:1027
void clear()
Definition: Matrix.h:835
static void gemm(const bool tA, const bool tB, const Treal alpha, const Matrix< Treal, Telement > &A, const Matrix< Treal, Telement > &B, const Treal beta, Matrix< Treal, Telement > &C)
Definition: Matrix.h:1084
Treal frob_thresh(Treal const threshold, Matrix< Treal, Telement > *ErrorMatrix=0)
Removes small elements so that the introduced error is smaller than the threshold in the Frobenius no...
Definition: Matrix.h:396
Treal trace() const
Definition: Matrix.h:1951
Treal syAccumulateWith(Top &op)
Definition: Matrix.h:457
static void swap(MatrixHierarchicBase< Treal, Telement > &A, MatrixHierarchicBase< Treal, Telement > &B)
Definition: MatrixHierarchicBase.h:223