77 template<
class Treal,
class Telement>
79 template<
typename Tperm>
114 template<
class Treal,
class Telement = Treal>
115 class Matrix:
public MatrixHierarchicBase<Treal, Telement> {
120 friend class Vector<Treal, Telement>;
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,
388 const bool tA,
const Treal alpha,
396 (Treal
const threshold,
407 (Treal
const threshold,
416 const Treal threshold = 0,
420 void gershgorin(Treal& lmin, Treal& lmax)
const;
425 tmp.gershgorin(lmin, lmax);
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>
519 fullMatrix(std::vector<Treal> & fullMat)
const {
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>
640 addValues(std::vector<int>
const & rowind,
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>
653 addValues(std::vector<int>
const & rowind,
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>
733 getValues(std::vector<int>
const & rowind,
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>
746 getValues(std::vector<int>
const & rowind,
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>
988 if (AT.is_zero() || (AT.nElements() !=
A.nElements())) {
994 for (
int row = 0; row < AT.nrows(); row++)
995 for (
int col = 0; col < AT.ncols(); col++)
996 Telement::transpose(
A(col,row), AT(row,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++)
1011 Telement::transpose((*
this)(col, row), (*
this)(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());
1094 C.resetRows(
A.cols);
1096 C.resetRows(
A.rows);
1098 C.resetCols(
B.rows);
1100 C.resetCols(
B.cols);
1103 if ((
A.is_zero() ||
B.is_zero() || alpha == 0) &&
1104 (C.is_zero() || beta == 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++) {
1125 Telement::gemm(tA, tB, alpha,
1126 A(row, 0),
B(0, col),
1129 for(
int cola = 1; cola <
A.ncols(); cola++)
1130 Telement::gemm(tA, tB, alpha,
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++) {
1156 Telement::gemm(tA, tB, alpha,
1160 for(
int cola = 1; cola <
A.nrows(); cola++)
1161 Telement::gemm(tA, tB, alpha,
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++) {
1187 Telement::gemm(tA, tB, alpha,
1188 A(row, 0),
B(col, 0),
1191 for(
int cola = 1; cola <
A.ncols(); cola++)
1192 Telement::gemm(tA, tB, alpha,
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++) {
1218 Telement::gemm(tA, tB, alpha,
1219 A(0, row),
B(col, 0),
1222 for(
int cola = 1; cola <
A.nrows(); cola++)
1223 Telement::gemm(tA, tB, alpha,
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();
1265 C.resetRows(
A.rows);
1266 C.resetCols(
B.cols);
1269 assert(
side ==
'R');
1270 C.resetRows(
B.rows);
1271 C.resetCols(
A.cols);
1274 if ((
A.is_zero() ||
B.is_zero() || alpha == 0) &&
1275 (C.is_zero() || beta == 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++) {
1298 Telement::symm(
side, uplo, alpha,
A(row, row),
1299 B(row, col), beta_tmp, C(row, col));
1301 for(
int ind = 0; ind < row; ind++)
1302 Telement::gemm(
true,
false, alpha,
1303 A(ind, row),
B(ind, col),
1306 for(
int ind = row + 1; ind < dimA; ind++)
1307 Telement::gemm(
false,
false, alpha,
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++) {
1333 Telement::symm(
side, uplo, alpha,
A(col, col),
1334 B(row, col), beta_tmp, C(row, col));
1336 for(
int ind = col + 1; ind < dimA; ind++)
1337 Telement::gemm(
false,
true, alpha,
1338 B(row, ind),
A(col, ind),
1341 for(
int ind = 0; ind < col; ind++)
1342 Telement::gemm(
false,
false, alpha,
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());
1378 C.resetRows(
A.cols);
1379 C.resetCols(
A.cols);
1382 C.resetRows(
A.rows);
1383 C.resetCols(
A.rows);
1387 if (C.nrows() == C.ncols() &&
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++) {
1407 Telement::syrk(uplo, tA, alpha,
A(rc, 0), beta_tmp, C(rc, rc));
1408 for (
int cola = 1; cola <
A.ncols(); cola++)
1409 Telement::syrk(uplo, tA, alpha,
A(rc, cola), 1.0, C(rc, rc));
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++) {
1419 Telement::gemm(tA, !tA, alpha,
1420 A(row, 0),
A(col,0), beta_tmp, C(row,col));
1421 for (
int ind = 1; ind <
A.ncols(); ind++)
1422 Telement::gemm(tA, !tA, alpha,
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++) {
1440 Telement::syrk(uplo, tA, alpha,
A(0, rc), beta_tmp, C(rc, rc));
1441 for (
int rowa = 1; rowa <
A.nrows(); rowa++)
1442 Telement::syrk(uplo, tA, alpha,
A(rowa, rc), 1.0, C(rc, rc));
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++) {
1452 Telement::gemm(tA, !tA, alpha,
1453 A(0, row),
A(0, col), beta_tmp, C(row,col));
1454 for (
int ind = 1; ind <
A.nrows(); ind++)
1455 Telement::gemm(tA, !tA, alpha,
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());
1484 C.resetRows(
A.rows);
1485 C.resetCols(
A.cols);
1487 if (C.nrows() == C.ncols() &&
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++)
1509 Telement::syrk(uplo,
true, alpha,
A(cola, rc), 1.0, C(rc, rc));
1510 for (
int cola = rc + 1; cola <
A.ncols(); cola++)
1511 Telement::syrk(uplo,
false, alpha,
A(rc, cola), 1.0, C(rc, rc));
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++) {
1523 Telement::symm(
'L',
'U', alpha,
A(row, row),
A(row, col),
1524 beta_tmp, C(row, col));
1525 Telement::symm(
'R',
'U', alpha,
A(col, col),
A(row, col),
1528 for (
int ind = 0; ind < row; ind++)
1529 Telement::gemm(
true,
false, alpha,
1530 A(ind, row),
A(ind, col), 1.0, C(row, col));
1532 for (
int ind = row + 1; ind < col; ind++)
1533 Telement::gemm(
false,
false, alpha,
1534 A(row, ind),
A(ind, col), 1.0, C(row, col));
1536 for (
int ind = col + 1; ind <
A.ncols(); ind++)
1537 Telement::gemm(
false,
true, alpha,
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>
1562 ssmm(
const Treal alpha,
1567 assert(!
A.is_empty());
1568 assert(!
B.is_empty());
1571 C.resetRows(
A.rows);
1572 C.resetCols(
B.cols);
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) &&
1587 (C.is_zero() || beta == 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++)
1613 Telement::gemm(
true,
false,
1614 alpha,
A(ind,col),
B(ind, col),
1616 for (
int ind = col + 1; ind <
A.ncols(); ind++)
1618 Telement::gemm(
false,
true,
1619 alpha,
A(col, ind),
B(col, ind),
1622 for (
int row = 0; row < col; row++) {
1626 Telement::symm(
'L',
'U',
1627 alpha,
A(row, row),
B(row, col),
1628 beta_tmp, C(row, col));
1630 Telement::symm(
'R',
'U',
1631 alpha,
B(col, col),
A(row, col),
1633 for (
int ind = 0; ind < row; ind++)
1635 Telement::gemm(
true,
false,
1636 alpha,
A(ind, row),
B(ind, col),
1638 for (
int ind = row + 1; ind < col; ind++)
1640 Telement::gemm(
false,
false,
1641 alpha,
A(row, ind),
B(ind, col),
1643 for (
int ind = col + 1; ind <
A.ncols(); ind++)
1645 Telement::gemm(
false,
true,
1646 alpha,
A(row, ind),
B(col, ind),
1651 for (
int row = col + 1; row < C.nrows(); row++) {
1652 Telement::transpose(C(row, col), tmpSubMat);
1656 Telement::symm(
'L',
'U',
1657 alpha,
B(col, col),
A(col, row),
1658 beta_tmp, tmpSubMat);
1660 Telement::symm(
'R',
'U',
1661 alpha,
A(row, row),
B(col, row),
1663 for (
int ind = 0; ind < col; ind++)
1665 Telement::gemm(
true,
false,
1666 alpha,
B(ind, col),
A(ind, row),
1668 for (
int ind = col + 1; ind < row; ind++)
1670 Telement::gemm(
false,
false,
1671 alpha,
B(col, ind),
A(ind, row),
1673 for (
int ind = row + 1; ind <
B.nrows(); ind++)
1675 Telement::gemm(
false,
true,
1676 alpha,
B(col, ind),
A(row, ind),
1678 Telement::transpose(tmpSubMat, C(row, col));
1689 template<
class Treal,
class Telement>
1696 assert(!
A.is_empty());
1697 assert(!
B.is_empty());
1700 C.resetRows(
A.rows);
1701 C.resetCols(
B.cols);
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) &&
1717 (C.is_zero() || beta == 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++) {
1756 Telement::symm(
'L',
'U',
1757 alpha,
A(row, row),
B(row, col),
1758 beta_tmp, C(row, col));
1760 Telement::symm(
'R',
'U',
1761 alpha,
B(col, col),
A(row, col),
1763 for (
int ind = 0; ind < row; ind++)
1765 Telement::gemm(
true,
false,
1766 alpha,
A(ind, row),
B(ind, col),
1768 for (
int ind = row + 1; ind < col; ind++)
1770 Telement::gemm(
false,
false,
1771 alpha,
A(row, ind),
B(ind, col),
1773 for (
int ind = col + 1; ind <
A.ncols(); ind++)
1775 Telement::gemm(
false,
true,
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++) {
1805 Telement::trmm(
side, uplo, tA, alpha,
1806 A(col, col),
B(row,col));
1808 for (
int ind = 0; ind < col; ind++)
1809 Telement::gemm(
false, tA, alpha,
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++) {
1819 Telement::trmm(
side, uplo, tA, alpha,
1820 A(row,row),
B(row,col));
1821 for (
int ind = row + 1 ; ind <
B.nrows(); ind++)
1822 Telement::gemm(tA,
false, alpha,
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++) {
1834 Telement::trmm(
side, uplo, tA, alpha,
1835 A(col,col),
B(row,col));
1836 for (
int ind = col + 1; ind <
A.ncols(); ind++)
1837 Telement::gemm(
false, tA, alpha,
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++) {
1847 Telement::trmm(
side, uplo, tA, alpha,
1848 A(row,row),
B(row,col));
1849 for (
int ind = 0; ind < row; ind++)
1850 Telement::gemm(tA,
false, alpha,
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>
2025 add(
const Treal alpha,
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>
2043 assign(Treal
const alpha,
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;
2074 assert(!ErrorMatrix->is_empty());
2075 bool EMatIsZero =
true;
2076 if (!ErrorMatrix->is_zero() || !this->is_zero()) {
2077 if (ErrorMatrix->is_zero())
2078 ErrorMatrix->allocate();
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;
2119 assert(!ErrorMatrix->is_empty());
2120 bool EMatIsZero =
true;
2121 if (!ErrorMatrix->is_zero() || !this->is_zero()) {
2122 if (ErrorMatrix->is_zero())
2123 ErrorMatrix->allocate();
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());
2287 C.resetRows(
A.cols);
2289 C.resetRows(
A.rows);
2291 C.resetCols(
B.rows);
2293 C.resetCols(
B.cols);
2295 if ((
A.is_zero() ||
B.is_zero() || alpha == 0) &&
2296 (C.is_zero() || beta == 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++) {
2320 Telement::gemm(tA, tB, alpha,
2321 A(row, 0),
B(0, col),
2324 for(
int cola = 1; cola <
A.ncols(); cola++)
2325 Telement::gemm(tA, tB, alpha,
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++) {
2354 Telement::gemm(tA, tB, alpha,
2358 for(
int cola = 1; cola <
A.nrows(); cola++)
2359 Telement::gemm(tA, tB, alpha,
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++) {
2388 Telement::gemm(tA, tB, alpha,
2389 A(row, 0),
B(col, 0),
2392 for(
int cola = 1; cola <
A.ncols(); cola++)
2393 Telement::gemm(tA, tB, alpha,
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++) {
2422 Telement::gemm(tA, tB, alpha,
2423 A(0, row),
B(col, 0),
2426 for(
int cola = 1; cola <
A.nrows(); cola++)
2427 Telement::gemm(tA, tB, alpha,
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());
2464 assert(!Z.is_empty());
2465 if (alpha != 0 && !
A.is_zero() && !Z.is_zero()) {
2466 if (
A.nrows() ==
A.ncols() &&
2467 Z.nrows() == Z.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--) {
2483 Telement::trmm(
side,
'U',
false,
2484 alpha, Z(col, col),
A(row, col));
2486 Telement::symm(
'L',
'U', alpha,
A(row, row), Z(row, col),
2488 for (
int ind = 0; ind < row; ind++) {
2490 Telement::gemm(
true,
false, alpha,
A(ind, row), Z(ind, col),
2493 for (
int ind = row + 1; ind < col; ind++) {
2495 Telement::gemm(
false,
false, alpha,
A(row, ind), Z(ind, col),
2502 assert(
side ==
'L');
2504 for (
int col = 0; col <
A.ncols(); col++) {
2506 for (
int row = 0; row < col; row++) {
2508 Telement::trmm(
side,
'U',
false, alpha,
2509 Z(row, row),
A(row, col));
2511 Telement::symm(
'R',
'U', alpha,
A(col, col), Z(row, col),
2513 for (
int ind = row + 1; ind < col; ind++)
2515 Telement::gemm(
false,
false, alpha, Z(row, ind),
A(ind, col),
2517 for (
int ind = col + 1; ind <
A.nrows(); ind++)
2519 Telement::gemm(
false,
true, alpha, Z(row, ind),
A(col, 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++) {
2572 Telement::trmm(
side, uplo, tA, alpha,
2573 A(col,col),
B(row,col));
2574 for (
int ind = col + 1; ind <
A.ncols(); ind++)
2575 Telement::gemm(
false, tA, alpha,
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++) {
2592 Telement::trmm(
side, uplo, tA, alpha,
2593 A(row,row),
B(row,col));
2594 for (
int ind = 0; ind < row; ind++)
2595 Telement::gemm(tA,
false, alpha,
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()) {
2644 ErrorMatrix->resetRows(this->
rows);
2645 ErrorMatrix->resetCols(this->
cols);
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());
2707 Z.resetRows(
A.rows);
2708 Z.resetCols(
A.cols);
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++)
2724 Telement::gemm(
true,
false, 1.0,
2725 A(k,j), Z(k,i), 1.0, Ptmp);
2726 Telement::trmm(
'L',
'U',
true, 1.0, Z(j,j), Ptmp);
2728 for (
int k = 0; k < j; k++)
2729 Telement::gemm(
false,
false, -1.0,
2730 Z(k,j), Ptmp, 1.0, Z(k,i));
2732 Telement::trmm(
'L',
'U',
false, -1.0, Z(j,j), Ptmp);
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++) {
2743 Telement::trmm(
'R',
'U',
false, 1.0, Z(i,i), Z(k,i));
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++) {
2767 Telement::trmm(
'R',
'U',
false, 1.0, Z(i,i), Z(k,i));
2771 for (
int j = i + 1; j < Z.ncols(); j++) {
2776 for (
int k = 0; k < i; k++)
2777 Telement::gemm(
true,
false, 1.0,
2778 A(k,i), Z(k,j), 1.0, Ptmp);
2779 Telement::trmm(
'L',
'U',
true, 1.0, Z(i,i), Ptmp);
2781 for (
int k = 0; k < i; k++)
2782 Telement::gemm(
false,
false, -1.0,
2783 Z(k,i), Ptmp, 1.0, Z(k,j));
2785 Telement::trmm(
'L',
'U',
false, -1.0, Z(i,i), Ptmp);
2786 Telement::add(1.0, Ptmp, Z(i,j));
2790 if (threshold != 0) {
2791 for (
int k = 0; k < i; k++)
2792 Z(k,i).frob_thresh(myThresh);
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>
2929 friend class Vector<Treal, 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;
3295 tmp.gershgorin(lmin, lmax);
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);
3460 if (this->is_zero())
3462 for (
int col = 0; col < this->ncols(); col++)
3463 for (
int row = 0; row < this->nrows(); row++)
3465 fullMat[(colOffset + col) * nTotalRows + (rowOffset + row)];
3468 template<
class Treal>
3470 fullMatrix(std::vector<Treal> & fullMat)
const {
3473 fullMat.resize(nTotalRows * nTotalCols);
3476 if (this->is_zero()) {
3477 for (
int col = 0; col < this->nScalarsCols(); col++)
3478 for (
int row = 0; row < this->nScalarsRows(); row++)
3479 fullMat[(colOffset + col) * nTotalRows + (rowOffset + row)] = 0;
3482 for (
int col = 0; col < this->ncols(); col++)
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);
3497 if (this->is_zero()) {
3498 for (
int col = 0; col < this->nScalarsCols(); col++)
3499 for (
int row = 0; row <= col; row++) {
3500 fullMat[(colOffset + col) * nTotalRows + (rowOffset + row)] = 0;
3501 fullMat[(rowOffset + row) * nTotalRows + (colOffset + col)] = 0;
3505 for (
int col = 0; col < this->ncols(); col++)
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);
3523 if (this->is_zero()) {
3524 for (
int col = 0; col < this->nScalarsCols(); col++)
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;
3531 for (
int col = 0; col < this->ncols(); col++)
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) {
3554 assignFromSparse(rowind, colind, values, indexes);
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()) {
3567 if (this->is_zero())
3569 for (
int ind = 0; ind < this->nElements(); ++ind)
3570 this->elements[ind] = 0;
3571 std::vector<int>::const_iterator it;
3572 for ( it = indexes.begin(); it < indexes.end(); it++ ) {
3579 template<
class Treal>
3581 addValues(std::vector<int>
const & rowind,
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>
3595 addValues(std::vector<int>
const & rowind,
3596 std::vector<int>
const & colind,
3597 std::vector<Treal>
const & values,
3598 std::vector<int>
const & indexes) {
3599 if (indexes.empty())
3601 if (this->is_zero())
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 ");
3629 assignFromSparse(rowind, colind, values);
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 ");
3651 addValues(rowind, colind, values);
3654 template<
class Treal>
3656 getValues(std::vector<int>
const & rowind,
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>
3670 getValues(std::vector<int>
const & rowind,
3671 std::vector<int>
const & colind,
3672 std::vector<Treal> & values,
3673 std::vector<int>
const & indexes)
const {
3674 assert(!this->is_empty());
3675 if (indexes.empty())
3677 std::vector<int>::const_iterator it;
3678 if (this->is_zero()) {
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 ");
3707 getValues(rowind, colind, values);
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());
3717 if (!this->is_zero())
3718 for (
int col = 0; col < this->ncols(); col++)
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());
3733 if (!this->is_zero())
3734 for (
int col = 0; col < this->ncols(); ++col)
3735 for (
int row = 0; row <= col; ++row) {
3738 values.push_back((*
this)(row, col));
3743 template<
class Treal>
3749 template<
class Treal>
3754 if (this->is_zero()) {
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) {
3778 if (this->is_zero())
3780 file.read((
char*)this->elements,
sizeof(Treal) * this->nElements());
3783 throw Failure(
"Matrix<Treal>::"
3784 "readFromFile(std::ifstream & file):"
3785 "File corruption, int value not 0 or 1");
3789 template<
class Treal>
3791 if (this->is_zero())
3793 for (
int ind = 0; ind < this->nElements(); ind++)
3794 this->elements[ind] = rand() / (Treal)RAND_MAX;
3797 template<
class Treal>
3799 if (this->is_zero())
3802 for (
int col = 1; col < this->ncols(); col++)
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>
3813 if (!this->highestLevel() &&
3814 probabilityBeingZero > rand() / (Treal)RAND_MAX)
3820 template<
class Treal>
3823 if (!this->highestLevel() &&
3824 probabilityBeingZero > rand() / (Treal)RAND_MAX)
3830 template<
class Treal>
3831 template<
typename TRule>
3834 if (this->is_zero())
3836 for (
int col = 0; col < this->ncols(); col++)
3837 for (
int row = 0; row < this->nrows(); row++)
3838 (*
this)(row,col) = rule.set(this->rows.getOffset() + row,
3839 this->cols.getOffset() + col);
3841 template<
class Treal>
3842 template<
typename TRule>
3845 if (this->is_zero())
3848 for (
int col = 0; col < this->ncols(); col++)
3849 for (
int row = 0; row <= col; row++)
3850 (*
this)(row,col) = rule.set(this->rows.getOffset() + row,
3851 this->cols.getOffset() + col);
3855 template<
class Treal>
3858 if (this->is_empty())
3859 throw Failure(
"Matrix<Treal>::addIdentity(Treal): "
3860 "Cannot add identity to empty matrix.");
3861 if (this->ncols() != this->nrows())
3862 throw Failure(
"Matrix<Treal, Telement>::addIdentity(Treal): "
3863 "Matrix must be square to add identity");
3864 if (this->is_zero())
3866 for (
int ind = 0; ind < this->ncols(); ind++)
3867 (*
this)(ind,ind) += alpha;
3870 template<
class Treal>
3880 if (AT.is_zero() || (AT.nElements() !=
A.nElements())) {
3886 for (
int row = 0; row < AT.nrows(); row++)
3887 for (
int col = 0; col < AT.ncols(); col++)
3888 AT(row,col) =
A(col,row);
3891 template<
class Treal>
3894 if (this->nrows() == this->ncols()) {
3895 if (!this->is_zero()) {
3898 for (
int row = 1; row < this->nrows(); row++)
3899 for (
int col = 0; col < row; col++)
3900 (*
this)(row, col) = (*
this)(col, row);
3904 throw Failure(
"Matrix<Treal>::symToNosym(): "
3905 "Only quadratic matrices can be symmetric");
3908 template<
class Treal>
3911 if (this->nrows() == this->ncols()) {
3912 if (!this->is_zero()) {
3915 for (
int col = 0; col < this->ncols() - 1; col++)
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>
3933 if (this->ncols() != this->nrows())
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++)
3954 this->elements[ind] *= alpha;
3959 template<
class Treal>
3961 gemm(
const bool tA,
const bool tB,
const Treal alpha,
3966 assert(!
A.is_empty());
3967 assert(!
B.is_empty());
3971 C.resetRows(
A.cols);
3973 C.resetRows(
A.rows);
3975 C.resetCols(
B.rows);
3977 C.resetCols(
B.cols);
3980 if ((
A.is_zero() ||
B.is_zero() || alpha == 0) &&
3981 (C.is_zero() || beta == 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() &&
3994 B.ncols() == C.ncols())
3995 mat::gemm(
"N",
"N", &
A.nrows(), &
B.ncols(), &
A.ncols(), &alpha,
3996 A.elements, &
A.nrows(),
B.elements, &
B.nrows(),
3997 &beta_tmp, C.elements, &C.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() &&
4007 B.ncols() == C.ncols())
4008 mat::gemm(
"T",
"N", &
A.ncols(), &
B.ncols(), &
A.nrows(), &alpha,
4009 A.elements, &
A.nrows(),
B.elements, &
B.nrows(),
4010 &beta_tmp, C.elements, &C.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() &&
4020 B.nrows() == C.ncols())
4021 mat::gemm(
"N",
"T", &
A.nrows(), &
B.nrows(), &
A.ncols(), &alpha,
4022 A.elements, &
A.nrows(),
B.elements, &
B.nrows(),
4023 &beta_tmp, C.elements, &C.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() &&
4033 B.nrows() == C.ncols())
4034 mat::gemm(
"T",
"T", &
A.ncols(), &
B.nrows(), &
A.nrows(), &alpha,
4035 A.elements, &
A.nrows(),
B.elements, &
B.nrows(),
4036 &beta_tmp, C.elements, &C.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,
4061 assert(!
A.is_empty());
4062 assert(!
B.is_empty());
4063 assert(
A.nrows() ==
A.ncols());
4068 C.resetRows(
A.rows);
4069 C.resetCols(
B.cols);
4072 assert(
side ==
'R');
4073 C.resetRows(
B.rows);
4074 C.resetCols(
A.cols);
4078 if ((
A.is_zero() ||
B.is_zero() || alpha == 0) &&
4079 (C.is_zero() || beta == 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(),
4094 &beta_tmp, C.elements, &C.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,
4110 assert(!
A.is_empty());
4114 C.resetRows(
A.cols);
4115 C.resetCols(
A.cols);
4118 C.resetRows(
A.rows);
4119 C.resetCols(
A.rows);
4122 if (C.nrows() == C.ncols() &&
4123 ((!tA &&
A.nrows() == C.nrows()) || (tA &&
A.ncols() == C.nrows())))
4124 if (alpha != 0 && !
A.is_zero()) {
4125 Treal beta_tmp = beta;
4131 mat::syrk(&uplo,
"N", &C.nrows(), &
A.ncols(), &alpha,
4132 A.elements, &
A.nrows(), &beta_tmp,
4133 C.elements, &C.nrows());
4136 mat::syrk(&uplo,
"T", &C.nrows(), &
A.nrows(), &alpha,
4137 A.elements, &
A.nrows(), &beta_tmp,
4138 C.elements, &C.nrows());
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,
4153 assert(!
A.is_empty());
4156 C.resetRows(
A.rows);
4157 C.resetCols(
A.cols);
4166 template<
class Treal>
4168 ssmm(
const Treal alpha,
4173 assert(!
A.is_empty());
4174 assert(!
B.is_empty());
4177 C.resetRows(
A.rows);
4178 C.resetCols(
B.cols);
4189 template<
class Treal>
4197 ssmm(alpha,
A,
B, beta, C);
4202 template<
class Treal>
4204 trmm(
const char side,
const char uplo,
const bool tA,
4208 assert(!
A.is_empty());
4209 assert(!
B.is_empty());
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>
4229 assert(!this->is_empty());
4230 if (this->is_zero())
4234 for (
int i = 0; i < this->nElements(); i++)
4235 sum += this->elements[i] * this->elements[i];
4240 template<
class Treal>
4243 assert(!this->is_empty());
4244 if (this->nrows() != this->ncols())
4245 throw Failure(
"Matrix<Treal>::syFrobSquared(): "
4246 "Matrix must be have equal number of rows "
4247 "and cols to be symmetric");
4249 if (!this->is_zero()) {
4250 for (
int col = 1; col < this->ncols(); col++)
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>
4322 assert(!this->is_empty());
4323 if (this->nrows() != this->ncols())
4324 throw Failure(
"Matrix<Treal>::trace(): Matrix must be quadratic");
4325 if (this->is_zero())
4329 for (
int rc = 0; rc < this->ncols(); rc++)
4330 sum += (*
this)(rc,rc);
4335 template<
class Treal>
4339 assert(!
A.is_empty());
4340 assert(!
B.is_empty());
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>
4356 assert(!
A.is_empty());
4357 assert(!
B.is_empty());
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>
4373 assert(!
A.is_empty());
4374 assert(!
B.is_empty());
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>
4395 add(
const Treal alpha,
4398 assert(!
A.is_empty());
4399 assert(!
B.is_empty());
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>
4413 assign(Treal
const alpha,
4415 assert(!
A.is_empty());
4416 if (this->is_empty()) {
4417 this->resetRows(
A.rows);
4418 this->resetCols(
A.cols);
4426 template<
class Treal>
4429 if (!this->is_zero())
4430 frobsq.push_back(this->frobSquared());
4433 template<
class Treal>
4436 if (!this->is_zero())
4437 for (
int ind = 0; ind < this->nElements(); ind++)
4438 if ( this->elements[ind] != 0 )
4439 frobsq.push_back(this->elements[ind] * this->elements[ind]);
4443 template<
class Treal>
4448 if ((!this->is_zero() && this->frobSquared() <= threshold) ||
4449 (!ErrorMatrix->is_zero() &&
4450 ErrorMatrix->frobSquared() > threshold))
4453 else if (!this->is_zero() && this->frobSquared() <= threshold)
4457 template<
class Treal>
4461 assert(!this->is_empty());
4462 bool thisMatIsZero =
true;
4464 assert(!ErrorMatrix->is_empty());
4465 bool EMatIsZero =
true;
4466 if (!ErrorMatrix->is_zero() || !this->is_zero()) {
4467 if (ErrorMatrix->is_zero())
4468 ErrorMatrix->allocate();
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];
4491 ErrorMatrix->elements[ind] = 0;
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>
4595 assert(
side ==
'L');
4604 template<
class Treal>
4607 assert(!this->is_empty());
4608 if (ErrorMatrix && ErrorMatrix->is_empty()) {
4609 ErrorMatrix->resetRows(this->
rows);
4610 ErrorMatrix->resetCols(this->
cols);
4612 Treal fs = frobSquared();
4613 if (fs < threshold) {
4623 template<
class Treal>
4627 const Treal threshold,
const side looking,
4629 assert(!
A.is_empty());
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.");
4637 potrf(
"U", &
A.nrows(), Z.elements, &
A.nrows(), &info);
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");
4643 trtri(
"U",
"N", &
A.nrows(), Z.elements, &
A.nrows(), &info);
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>
4655 assert(!this->is_empty());
4656 if (this->nScalarsRows() == this->nScalarsCols()) {
4657 int size = this->nScalarsCols();
4658 Treal* colsums =
new Treal[size];
4659 Treal* diag =
new Treal[size];
4660 for (
int ind = 0; ind < size; ind++)
4662 this->add_abs_col_sums(colsums);
4663 this->get_diagonal(diag);
4666 lmin = diag[0] - tmp1;
4667 lmax = diag[0] + tmp1;
4668 for (
int col = 1; col < size; col++) {
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>
4687 if (!this->is_zero())
4688 for (
int col = 0; col < this->ncols(); col++)
4689 for (
int row = 0; row < this->nrows(); row++)
4693 template<
class Treal>
4697 assert(this->nrows() == this->ncols());
4698 if (this->is_zero())
4699 for (
int rc = 0; rc < this->nScalarsCols(); rc++)
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) {
4734 this->elements = NULL;
4738 else if (this->cap > this->nel) {
4739 Treal* tmp =
new Treal[this->nel];
4740 for (
int i = 0; i < this->nel; i++)
4741 tmp[i] = this->elements[i];
4743 this->cap = this->nel;
4744 this->elements = tmp;
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 {
4756 if (this->is_zero()) {
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) {
4771 if (this->is_zero()) {
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;
4805 this->nel = this->nrows() * this->ncols();
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");
Code for memory allocation/deallocation routines used by matrix library.
Base class for Matrix and Matrix specialization.
Definition MatrixHierarchicBase.h:52
const int & ncols() const
Definition MatrixHierarchicBase.h:71
static void swap(MatrixHierarchicBase< Treal, Telement > &A, MatrixHierarchicBase< Treal, Telement > &B)
Definition MatrixHierarchicBase.h:223
const int & nrows() const
Definition MatrixHierarchicBase.h:69
SizesAndBlocks cols
Definition MatrixHierarchicBase.h:165
SizesAndBlocks rows
Definition MatrixHierarchicBase.h:164
void resetRows(SizesAndBlocks const &newRows)
Definition MatrixHierarchicBase.h:114
const int & nScalarsCols() const
Definition MatrixHierarchicBase.h:65
Telement int col
Definition MatrixHierarchicBase.h:75
const int & nScalarsRows() const
Definition MatrixHierarchicBase.h:63
return elements[row+col *nrows()]
Definition MatrixHierarchicBase.h:81
int nElements() const
Definition MatrixHierarchicBase.h:110
bool is_empty() const
Definition MatrixHierarchicBase.h:143
MatrixHierarchicBase< Treal, Telement > & operator=(const MatrixHierarchicBase< Treal, Telement > &mat)
Definition MatrixHierarchicBase.h:202
bool is_zero() const
Definition MatrixHierarchicBase.h:108
void resetCols(SizesAndBlocks const &newCols)
Definition MatrixHierarchicBase.h:119
size_t sy_nnz() const
Returns number of nonzeros in matrix including lower triangle elements.
Definition Matrix.h:3315
Treal maxAbsValue() const
Definition Matrix.h:3374
void assignFrobNormsLowestLevel(Matrix< Treal, Matrix< Treal > > const &A)
Definition Matrix.h:3136
static const Treal ZERO
Definition Matrix.h:3441
Treal frob() const
Definition Matrix.h:3075
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:2931
void allocate()
Definition Matrix.h:2940
Matrix< Treal > & operator=(const Matrix< Treal > &mat)
Definition Matrix.h:2999
size_t memory_usage() const
Definition Matrix.h:3301
Treal syAccumulateWith(Top &op)
Definition Matrix.h:3339
static Treal frobDiff(const Matrix< Treal > &A, const Matrix< Treal > &B)
Definition Matrix.h:3083
static unsigned int level()
Definition Matrix.h:3372
Treal frob_thresh(Treal const threshold, Matrix< Treal > *ErrorMatrix=0)
Definition Matrix.h:3267
size_t sy_nvalues() const
Returns number of stored values in matrix.
Definition Matrix.h:3329
static const Treal ONE
Definition Matrix.h:3442
size_t nnz() const
Returns number of nonzeros in matrix.
Definition Matrix.h:3309
void assignDiffFrobNormsLowestLevel(Matrix< Treal, Matrix< Treal > > const &A, Matrix< Treal, Matrix< Treal > > const &B)
Definition Matrix.h:3163
~Matrix()
Definition Matrix.h:3005
size_t nvalues() const
Returns number of stored values in matrix.
Definition Matrix.h:3324
Treal syFrob() const
Definition Matrix.h:3077
Treal geAccumulateWith(Top &op)
Definition Matrix.h:3358
static Treal syFrobDiff(const Matrix< Treal > &A, const Matrix< Treal > &B)
Definition Matrix.h:3092
void sy_gershgorin(Treal &lmin, Treal &lmax) const
Definition Matrix.h:3292
void syAssignFrobNormsLowestLevel(Matrix< Treal, Matrix< Treal > > const &A)
Definition Matrix.h:3148
void syAssignDiffFrobNormsLowestLevel(Matrix< Treal, Matrix< Treal > > const &A, Matrix< Treal, Matrix< Treal > > const &B)
Definition Matrix.h:3195
void truncateAccordingToSparsityPattern(Matrix< Treal, Matrix< Treal > > &A) const
Definition Matrix.h:3235
Matrix class and heart of the matrix library.
Definition Vector.h:44
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 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
size_t memory_usage() const
Definition Matrix.h:2863
void writeToFile(std::ofstream &file) const
Definition Matrix.h:845
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
static void transpose(Matrix< Treal, Telement > const &A, Matrix< Treal, Telement > &AT)
Definition Matrix.h:979
void sySetElementsByRule(TRule &rule)
Definition Matrix.h:949
void syRandom()
Definition Matrix.h:892
size_t nnz() const
Returns number of nonzeros in matrix.
Definition Matrix.h:2874
Treal syFrobSquared() const
Definition Matrix.h:1882
static Treal frobDiff(const Matrix< Treal, Telement > &A, const Matrix< Treal, Telement > &B)
Definition Matrix.h:297
void syGetValues(std::vector< int > const &rowind, std::vector< int > const &colind, std::vector< Treal > &values) const
Definition Matrix.h:786
void syFullMatrix(std::vector< Treal > &fullMat) const
Definition Matrix.h:539
Matrix< Treal, Telement > & operator=(const Matrix< Treal, Telement > &mat)
Definition Matrix.h:198
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
Treal syFrob() const
Definition Matrix.h:291
void trSetElementsByRule(TRule &rule)
Definition Matrix.h:224
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
void getFrobSqLowestLevel(std::vector< Treal > &frobsq) const
Definition Matrix.h:2061
static Treal syFrobSquaredDiff(const Matrix< Treal, Telement > &A, const Matrix< Treal, Telement > &B)
Definition Matrix.h:1924
static Treal syFrobDiff(const Matrix< Treal, Telement > &A, const Matrix< Treal, Telement > &B)
Definition Matrix.h:306
void fullMatrix(std::vector< Treal > &fullMat) const
Definition Matrix.h:519
void sy_gershgorin(Treal &lmin, Treal &lmax) const
Definition Matrix.h:422
void frobThreshElementLevel(Treal const threshold, Matrix< Treal, Telement > *ErrorMatrix)
Definition Matrix.h:2115
void setElementsByRule(TRule &rule)
Definition Matrix.h:940
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
static void add(const Treal alpha, const Matrix< Treal, Telement > &A, Matrix< Treal, Telement > &B)
Definition Matrix.h:2025
void getValues(std::vector< int > const &rowind, std::vector< int > const &colind, std::vector< Treal > &values) const
Definition Matrix.h:733
void addIdentity(Treal alpha)
Definition Matrix.h:964
static void trsytriplemm(char const side, const Matrix< Treal, Telement > &Z, Matrix< Treal, Telement > &A)
Definition Matrix.h:2620
size_t nvalues() const
Returns number of stored values in matrix.
Definition Matrix.h:444
void assignFromFull(std::vector< Treal > const &fullMat)
Definition Matrix.h:506
Treal trace() const
Definition Matrix.h:1951
void randomZeroStructure(Treal probabilityBeingZero)
Get a random zero structure with a specified probability that each submatrix is zero.
Definition Matrix.h:906
void gershgorin(Treal &lmin, Treal &lmax) const
Definition Matrix.h:2800
static Treal sy_trace_ab(const Matrix< Treal, Telement > &A, const Matrix< Treal, Telement > &B)
Definition Matrix.h:2002
size_t sy_nvalues() const
Returns number of stored values in matrix.
Definition Matrix.h:2900
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
void frobThreshLowestLevel(Treal const threshold, Matrix< Treal, Telement > *ErrorMatrix)
Definition Matrix.h:2070
void getAllValues(std::vector< int > &rowind, std::vector< int > &colind, std::vector< Treal > &) const
Definition Matrix.h:808
Treal syAccumulateWith(Top &op)
Definition Matrix.h:457
void random()
Definition Matrix.h:884
void symToNosym()
Definition Matrix.h:1001
void add_abs_col_sums(Treal *abscolsums) const
Definition Matrix.h:2832
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 get_diagonal(Treal *diag) const
Definition Matrix.h:2847
void getFrobSqElementLevel(std::vector< Treal > &frobsq) const
Definition Matrix.h:2106
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
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
static Treal trace_ab(const Matrix< Treal, Telement > &A, const Matrix< Treal, Telement > &B)
Definition Matrix.h:1968
Treal maxAbsValue() const
Definition Matrix.h:482
void nosymToSym()
Definition Matrix.h:1027
static unsigned int level()
Definition Matrix.h:480
void syAssignFrobNormsLowestLevel(Matrix< Treal, Matrix< Treal, Telement > > const &A)
Version of assignFrobNormsLowestLevelToMatrix for symmetric matrices.
Definition Matrix.h:2167
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 clear()
Definition Matrix.h:835
void syGetAllValues(std::vector< int > &rowind, std::vector< int > &colind, std::vector< Treal > &) const
Definition Matrix.h:820
Treal geAccumulateWith(Top &op)
Accumulation algorithm for general matrices.
Definition Matrix.h:470
static Treal trace_aTb(const Matrix< Treal, Telement > &A, const Matrix< Treal, Telement > &B)
Definition Matrix.h:1985
static Treal frobSquaredDiff(const Matrix< Treal, Telement > &A, const Matrix< Treal, Telement > &B)
Definition Matrix.h:1902
Telement ElementType
Definition Matrix.h:117
Vector< Treal, typename ElementType::VectorType > VectorType
Definition Matrix.h:118
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
Matrix()
Definition Matrix.h:121
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
size_t sy_nnz() const
Returns number of nonzeros in matrix including lower triangle elements.
Definition Matrix.h:2883
void assign(Treal const alpha, Matrix< Treal, Telement > const &A)
Definition Matrix.h:2043
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
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
void readFromFile(std::ifstream &file)
Definition Matrix.h:861
~Matrix()
Definition Matrix.h:205
void syRandomZeroStructure(Treal probabilityBeingZero)
Definition Matrix.h:920
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
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
Matrix< Treal, Telement > & operator*=(const Treal alpha)
Definition Matrix.h:1073
Treal frob() const
Definition Matrix.h:286
void allocate()
Definition Matrix.h:124
void syAddValues(std::vector< int > const &rowind, std::vector< int > const &colind, std::vector< Treal > const &values)
Definition Matrix.h:711
Treal frobSquared() const
Definition Matrix.h:1868
void syAssignFromSparse(std::vector< int > const &rowind, std::vector< int > const &colind, std::vector< Treal > const &values)
Definition Matrix.h:689
void addValues(std::vector< int > const &rowind, std::vector< int > const &colind, std::vector< Treal > const &values)
Definition Matrix.h:640
void syUpTriFullMatrix(std::vector< Treal > &fullMat) const
Definition Matrix.h:563
void assignFromSparse(std::vector< int > const &rowind, std::vector< int > const &colind, std::vector< Treal > const &values)
Definition Matrix.h:588
static SingletonForTimings & instance()
Definition Matrix.h:96
double getAccumulatedTime()
Definition Matrix.h:87
double accumulatedTime
Definition Matrix.h:84
SingletonForTimings(const SingletonForTimings &other)
void addTime(double timeToAdd)
Definition Matrix.h:88
void reset()
Definition Matrix.h:86
SingletonForTimings()
Definition Matrix.h:102
Describes dimensions of matrix and its blocks on all levels.
Definition SizesAndBlocks.h:45
int const & getNBlocks() const
Definition SizesAndBlocks.h:72
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
int getNTotalScalars() const
Definition SizesAndBlocks.h:84
int getOffset() const
Definition SizesAndBlocks.h:83
SizesAndBlocks getSizesAndBlocksForLowerLevel(int const blockNumber) const
Definition SizesAndBlocks.cc:71
Definition matInclude.h:156
void tic()
Definition matInclude.cc:84
float toc()
Definition matInclude.cc:88
Vector class.
Definition Vector.h:54
mat::SizesAndBlocks rows
Definition test.cc:51
mat::SizesAndBlocks cols
Definition test.cc:52
#define MAT_OMP_FINALIZE
Definition matInclude.h:92
#define MAT_OMP_INIT
Definition matInclude.h:89
#define MAT_OMP_END
Definition matInclude.h:91
#define MAT_OMP_START
Definition matInclude.h:90
C++ interface to a subset of BLAS and LAPACK.
Proxy structs used by the matrix API.
Definition allocate.cc:39
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
static void trtri(const char *uplo, const char *diag, const int *n, T *A, const int *lda, int *info)
Definition mat_gblas.h:321
side
Definition Matrix.h:75
@ right
Definition Matrix.h:75
@ left
Definition Matrix.h:75
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
void freeElements(float *ptr)
Definition allocate.cc:49
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
@ ful
Definition matInclude.h:138
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
T * allocateElements(int n)
Definition allocate.h:42
static void potrf(const char *uplo, const int *n, T *A, const int *lda, int *info)
Definition mat_gblas.h:314
Treal trace(const XYZ< Treal, MatrixGeneral< Treal, Tmatrix >, MatrixGeneral< Treal, Tmatrix > > &smm)
Definition MatrixGeneral.h:904
inchversion
Definition Matrix.h:76
@ unstable
Definition Matrix.h:76
@ stable
Definition Matrix.h:76
static void trifulltofull(Treal *trifull, const int size)
Definition mat_gblas.h:1193
A quicksort implementation.
Treal template_blas_sqrt(Treal x)
Treal template_blas_fabs(Treal x)