77 template<
class Treal,
class Telement>
79 template<
typename Tperm>
114 template<
class Treal,
class Telement = Treal>
120 friend class Vector<Treal, Telement>;
128#ifdef MAT_USE_ALLOC_TIMER
132#ifdef MAT_USE_ALLOC_TIMER
138 colSAB = this->
cols.getSizesAndBlocksForLowerLevel(
col);
139 for (
int row = 0; row < this->
rows.getNBlocks(); row++) {
141 rowSAB = this->
rows.getSizesAndBlocksForLowerLevel(row);
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);
164 std::vector<int>
const & colind,
165 std::vector<Treal>
const & values);
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);
180 std::vector<int>
const & colind,
181 std::vector<Treal> & values)
const;
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,
343 (Treal
const threshold,
388 const bool tA,
const Treal alpha,
396 (Treal
const threshold,
407 (Treal
const threshold,
416 const Treal threshold = 0,
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>
507 int nTotalRows = this->
rows.getNTotalScalars();
508 int nTotalCols = this->
cols.getNTotalScalars();
509 assert((
int)fullMat.size() == nTotalRows * nTotalCols);
513 for (
int row = 0; row < this->
nrows(); row++)
517 template<
class Treal,
class Telement>
519 fullMatrix(std::vector<Treal> & fullMat)
const {
520 int nTotalRows = this->
rows.getNTotalScalars();
521 int nTotalCols = this->
cols.getNTotalScalars();
522 fullMat.resize(nTotalRows * nTotalCols);
524 int rowOffset = this->
rows.getOffset();
525 int colOffset = this->
cols.getOffset();
528 fullMat[(colOffset +
col) * nTotalRows + (rowOffset + row)] = 0;
532 for (
int row = 0; row < this->
nrows(); row++)
537 template<
class Treal,
class Telement>
540 int nTotalRows = this->
rows.getNTotalScalars();
541 int nTotalCols = this->
cols.getNTotalScalars();
542 fullMat.resize(nTotalRows * nTotalCols);
544 int rowOffset = this->
rows.getOffset();
545 int colOffset = this->
cols.getOffset();
547 for (
int row = 0; row <=
col; row++) {
548 fullMat[(colOffset +
col) * nTotalRows + (rowOffset + row)] = 0;
549 fullMat[(rowOffset + row) * nTotalRows + (colOffset +
col)] = 0;
554 for (
int row = 0; row <
col; row++)
561 template<
class Treal,
class Telement>
564 int nTotalRows = this->
rows.getNTotalScalars();
565 int nTotalCols = this->
cols.getNTotalScalars();
566 fullMat.resize(nTotalRows * nTotalCols);
568 int rowOffset = this->
rows.getOffset();
569 int colOffset = this->
cols.getOffset();
572 fullMat[(colOffset +
col) * nTotalRows + (rowOffset + row)] = 0;
573 fullMat[(rowOffset + row) * nTotalRows + (colOffset +
col)] = 0;
578 for (
int row = 0; row < this->
nrows(); row++)
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)
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);
624 while (!columnBuckets[
col].empty()) {
625 currentInd = columnBuckets[
col].back();
626 columnBuckets[
col].pop_back();
627 rowBuckets[ this->
rows.whichBlock
628 ( rowind[currentInd] ) ].push_back (currentInd);
631 for (
int row = 0; row < this->
nrows(); 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);
673 while (!columnBuckets[
col].empty()) {
674 currentInd = columnBuckets[
col].back();
675 columnBuckets[
col].pop_back();
676 rowBuckets[ this->
rows.whichBlock
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 ");
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 ");
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 {
753 std::vector<int>::const_iterator it;
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);
770 while (!columnBuckets[
col].empty()) {
771 currentInd = columnBuckets[
col].back();
772 columnBuckets[
col].pop_back();
773 rowBuckets[ this->
rows.whichBlock
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 ");
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());
814 for (
int ind = 0; ind < this->
nElements(); ++ind)
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());
827 for (
int row = 0; row <
col; ++row)
834 template<
class Treal,
class Telement>
837 for (
int i = 0; i < this->
nElements(); i++)
843 template<
class Treal,
class Telement>
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++)
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++)
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++)
891 template<
class Treal,
class Telement>
897 for (
int row = 0; row <
col; row++)
900 for (
int rc = 0; rc < this->
nrows(); rc++)
904 template<
class Treal,
class Telement>
908 probabilityBeingZero > rand() / (Treal)RAND_MAX)
913 for (
int ind = 0; ind < this->
nElements(); ind++)
918 template<
class Treal,
class Telement>
922 probabilityBeingZero > rand() / (Treal)RAND_MAX)
929 for (
int row = 0; row <
col; row++)
932 for (
int rc = 0; rc < this->
nrows(); rc++)
937 template<
class Treal,
class Telement>
938 template<
typename TRule>
943 for (
int ind = 0; ind < this->
nElements(); ind++)
946 template<
class Treal,
class Telement>
947 template<
typename TRule>
954 for (
int row = 0; row <
col; row++)
957 for (
int rc = 0; rc < this->
nrows(); rc++)
962 template<
class Treal,
class Telement>
966 throw Failure(
"Matrix<Treal, Telement>::addIdentity(Treal): "
967 "Cannot add identity to empty matrix.");
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++)
977 template<
class Treal,
class Telement>
994 for (
int row = 0; row < AT.
nrows(); row++)
996 Telement::transpose(
A(
col,row), AT(row,
col));
999 template<
class Treal,
class Telement>
1006 for (
int rc = 0; rc < this->
ncols(); rc++)
1009 for (
int row = 1; row < this->
nrows(); row++)
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>
1031 for (
int rc = 0; rc < this->
ncols(); rc++)
1034 for (
int col = 0;
col < this->ncols() - 1;
col++)
1035 for (
int row =
col + 1; row < this->
nrows(); row++)
1040 throw Failure(
"Matrix<Treal, Telement>::nosymToSym(): "
1041 "Only quadratic matrices can be symmetric");
1046 template<
class Treal,
class Telement>
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++)
1082 template<
class Treal,
class Telement>
1084 gemm(
const bool tA,
const bool tB,
const Treal alpha,
1103 if ((
A.is_zero() ||
B.is_zero() || alpha == 0) &&
1107 Treal beta_tmp = beta;
1112 if (!
A.is_zero() && !
B.is_zero() && alpha != 0) {
1115 if (
A.ncols() ==
B.nrows() &&
1116 A.nrows() == C.
nrows() &&
1117 B.ncols() == C.
ncols()) {
1118 int C_ncols = C.
ncols();
1120#pragma omp parallel for if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) schedule(dynamic)
1122 for (
int col = 0;
col < C_ncols;
col++) {
1124 for (
int row = 0; row < C.
nrows(); row++) {
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,
1261 int dimA =
A.nrows();
1274 if ((
A.is_zero() ||
B.is_zero() || alpha == 0) &&
1279 Treal beta_tmp = beta;
1284 if (!
A.is_zero() && !
B.is_zero() && alpha != 0) {
1287 if (dimA ==
B.nrows() &&
1288 dimA == C.
nrows() &&
1289 B.ncols() == C.
ncols()) {
1290 int C_ncols = C.
ncols();
1292#pragma omp parallel for if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) schedule(dynamic)
1294 for (
int col = 0;
col < C_ncols;
col++) {
1296 for (
int row = 0; row < C.
nrows(); row++) {
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");
1322 if (
B.ncols() == dimA &&
1323 B.nrows() == C.
nrows() &&
1324 dimA == C.
ncols()) {
1325 int C_ncols = C.
ncols();
1327#pragma omp parallel for if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) schedule(dynamic)
1329 for (
int col = 0;
col < C_ncols;
col++) {
1331 for (
int row = 0; row < C.
nrows(); row++) {
1334 B(row,
col), beta_tmp, C(row,
col));
1336 for(
int ind =
col + 1; ind < dimA; ind++)
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,
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++) {
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++) {
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,
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++) {
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,
1574 if (
A.ncols() !=
B.nrows() ||
1575 A.nrows() != C.
nrows() ||
1576 B.ncols() != C.
ncols() ||
1577 A.nrows() !=
A.ncols() ||
1578 B.nrows() !=
B.ncols()) {
1579 throw Failure(
"Matrix<class Treal, class Telement>::"
1581 "Matrix<Treal, Telement>, "
1582 "Matrix<Treal, Telement>, Treal, "
1583 "Matrix<Treal, Telement>): "
1584 "Incorrect matrixdimensions for ssmm multiplication");
1586 if ((
A.is_zero() ||
B.is_zero() || alpha == 0) &&
1591 Treal beta_tmp = beta;
1596 if (
A.is_zero() ||
B.is_zero() || alpha == 0) {
1601 int C_ncols = C.
ncols();
1603#pragma omp parallel for if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) schedule(dynamic)
1605 for (
int col = 0;
col < C_ncols;
col++) {
1611 for (
int ind = 0; ind <
col; ind++)
1613 Telement::gemm(
true,
false,
1616 for (
int ind =
col + 1; ind <
A.ncols(); ind++)
1618 Telement::gemm(
false,
true,
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',
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',
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>
1703 if (
A.ncols() !=
B.nrows() ||
1704 A.nrows() != C.
nrows() ||
1705 B.ncols() != C.
ncols() ||
1706 A.nrows() !=
A.ncols() ||
1707 B.nrows() !=
B.ncols()) {
1708 throw Failure(
"Matrix<class Treal, class Telement>::"
1709 "ssmm_upper_tr_only(Treal, "
1710 "Matrix<Treal, Telement>, "
1711 "Matrix<Treal, Telement>, Treal, "
1712 "Matrix<Treal, Telement>): "
1713 "Incorrect matrixdimensions for ssmm_upper_tr_only "
1716 if ((
A.is_zero() ||
B.is_zero() || alpha == 0) &&
1721 Treal beta_tmp = beta;
1726 if (
A.is_zero() ||
B.is_zero() || alpha == 0) {
1731 int C_ncols = C.
ncols();
1733#pragma omp parallel for if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) schedule(dynamic)
1735 for (
int col = 0;
col < C_ncols;
col++) {
1741 for (
int ind = 0; ind <
col; ind++)
1743 Telement::gemm_upper_tr_only(
true,
false,
1746 for (
int ind =
col + 1; ind <
A.ncols(); ind++)
1748 Telement::gemm_upper_tr_only(
false,
true,
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',
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,
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,
1808 for (
int ind = 0; ind <
col; ind++)
1809 Telement::gemm(
false, tA, alpha,
1810 B(row,ind),
A(ind,
col),
1817 for (
int row = 0; row <
B.nrows(); row++ )
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),
1833 for (
int row = 0; row <
B.nrows(); row++) {
1834 Telement::trmm(
side, uplo, tA, alpha,
1836 for (
int ind =
col + 1; ind <
A.ncols(); ind++)
1837 Telement::gemm(
false, tA, alpha,
1838 B(row,ind),
A(
col,ind),
1845 for (
int row =
B.nrows() - 1; row >= 0; row--)
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>
1874 for (
int i = 0; i < this->
nElements(); i++)
1880 template<
class Treal,
class Telement>
1885 throw Failure(
"Matrix<Treal, Telement>::syFrobSquared(): "
1886 "Matrix must be have equal number of rows "
1887 "and cols to be symmetric");
1891 for (
int row = 0; row <
col; row++)
1893 for (
int rc = 0; rc < this->ncols(); rc++)
1899 template<
class Treal,
class Telement>
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>
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()) {
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>
1954 throw Failure(
"Matrix<Treal, Telement>::trace(): "
1955 "Matrix must be quadratic");
1960 for (
int rc = 0; rc < this->
ncols(); rc++)
1961 sum += (*
this)(rc,rc).
trace();
1966 template<
class Treal,
class Telement>
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>
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>
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,
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,
2052 add(alpha,
A, *
this);
2059 template<
class Treal,
class Telement>
2063 for (
int ind = 0; ind < this->
nElements(); ind++)
2067 template<
class Treal,
class Telement>
2072 bool thisMatIsZero =
true;
2075 bool EMatIsZero =
true;
2076 if (!ErrorMatrix->
is_zero() || !this->is_zero()) {
2081 for (
int ind = 0; ind < this->
nElements(); ind++) {
2084 thisMatIsZero = thisMatIsZero && this->
elements[ind].is_zero();
2085 EMatIsZero = EMatIsZero && ErrorMatrix->
elements[ind].is_zero();
2090 ErrorMatrix->
clear();
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>
2108 for (
int ind = 0; ind < this->
nElements(); ind++)
2112 template<
class Treal,
class Telement>
2117 bool thisMatIsZero =
true;
2120 bool EMatIsZero =
true;
2121 if (!ErrorMatrix->
is_zero() || !this->is_zero()) {
2126 for (
int ind = 0; ind < this->
nElements(); ind++) {
2129 thisMatIsZero = thisMatIsZero && this->
elements[ind].is_zero();
2130 EMatIsZero = EMatIsZero && ErrorMatrix->
elements[ind].is_zero();
2135 ErrorMatrix->
clear();
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>
2158 for (
int ind = 0; ind < this->
nElements(); ind++)
2165 template<
class Treal,
class Telement>
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");
2178 for (
int row = 0; row <
col; row++)
2180 for (
int rc = 0; rc < this->ncols(); rc++)
2187 template<
class Treal,
class Telement>
2191 if (
A.is_zero() &&
B.is_zero() ) {
2199 if ( !
A.is_zero() && !
B.is_zero() ) {
2203 for (
int ind = 0; ind < this->
nElements(); ind++)
2207 if ( !
A.is_zero() ) {
2212 if ( !
B.is_zero() ) {
2218 template<
class Treal,
class Telement>
2222 if (
A.is_zero() &&
B.is_zero() ) {
2230 if ( !
A.is_zero() && !
B.is_zero() ) {
2235 for (
int row = 0; row <
col; row++)
2237 for (
int rc = 0; rc < this->ncols(); rc++)
2241 if ( !
A.is_zero() ) {
2246 if ( !
B.is_zero() ) {
2255 template<
class Treal,
class Telement>
2264 for (
int ind = 0; ind < this->
nElements(); ind++)
2275 template<
class Treal,
class Telement>
2295 if ((
A.is_zero() ||
B.is_zero() || alpha == 0) &&
2299 Treal beta_tmp = beta;
2304 if (!
A.is_zero() && !
B.is_zero() && alpha != 0) {
2306 if (
A.ncols() ==
B.nrows() &&
2307 A.nrows() == C.
nrows() &&
2308 B.ncols() == C.
ncols()) {
2310 Telement::gemm_upper_tr_only(tA, tB, alpha,
2314 for(
int cola = 1; cola <
A.ncols(); cola++)
2315 Telement::gemm_upper_tr_only(tA, tB, alpha,
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()) {
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,
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()) {
2378 Telement::gemm_upper_tr_only(tA, tB, alpha,
2382 for(
int cola = 1; cola <
A.ncols(); cola++)
2383 Telement::gemm_upper_tr_only(tA, tB, alpha,
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()) {
2412 Telement::gemm_upper_tr_only(tA, tB, alpha,
2416 for(
int cola = 1; cola <
A.nrows(); cola++)
2417 Telement::gemm_upper_tr_only(tA, tB, alpha,
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>
2465 if (alpha != 0 && !
A.is_zero() && !Z.
is_zero()) {
2466 if (
A.nrows() ==
A.ncols() &&
2468 A.nrows() == Z.
nrows()) {
2471 for (
int col =
A.ncols() - 1;
col >= 0;
col--) {
2473 Telement::sytr_upper_tr_only(
side, alpha,
2475 for (
int ind = 0; ind <
col; ind++) {
2477 Telement::gemm_upper_tr_only(
true,
false, alpha,
A(ind,
col),
2481 for (
int row =
col - 1; row >= 0; row--) {
2483 Telement::trmm(
side,
'U',
false,
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),
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,
2524 for (
int ind =
col + 1; ind <
A.ncols(); ind++)
2525 Telement::gemm_upper_tr_only(
false,
true, alpha, Z(
col, ind),
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,
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");
2565 Telement::trmm_upper_tr_only(
side, uplo, tA, alpha,
2567 for (
int ind =
col + 1; ind <
A.ncols(); ind++)
2568 Telement::gemm_upper_tr_only(
false, tA, alpha,
2571 for (
int row = 0; row <
col; row++) {
2572 Telement::trmm(
side, uplo, tA, alpha,
2574 for (
int ind =
col + 1; ind <
A.ncols(); ind++)
2575 Telement::gemm(
false, tA, alpha,
2576 B(row,ind),
A(
col,ind),
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>
2638 template<
class Treal,
class Telement>
2643 if (ErrorMatrix && ErrorMatrix->
is_empty()) {
2647 assert(threshold >= (Treal)0.0);
2648 if (threshold == (Treal)0.0)
2651 std::vector<Treal> 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;
2688 template<
class Treal,
class Telement>
2692 const Treal threshold,
const side looking,
2695 if (
A.nrows() !=
A.ncols())
2696 throw Failure(
"Matrix<Treal, Telement>::sy_inch: "
2697 "Matrix must be quadratic!");
2699 throw Failure(
"Matrix<Treal>::sy_inch: Matrix is zero! "
2700 "Not possible to compute inverse cholesky.");
2702 throw Failure(
"Matrix<Treal>::sy_inch: Only unstable "
2703 "version of sy_inch implemented.");
2706 myThresh = threshold / (Z.
ncols() * Z.
nrows());
2711 if (looking ==
left) {
2713 throw Failure(
"Matrix<Treal, Telement>::syInch: "
2714 "Thresholding not implemented for left-looking inch.");
2716 Telement::syInch(
A(0,0), Z(0,0), threshold, looking, version);
2718 for (
int i = 1; i < Z.
ncols(); i++) {
2719 for (
int j = 0; j < i; j++) {
2723 for (
int k = 0; k < j; k++)
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));
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++)
2798 template<
class Treal,
class Telement>
2804 Treal* colsums =
new Treal[size];
2805 Treal* diag =
new Treal[size];
2806 for (
int ind = 0; ind < size; ind++)
2812 lmin = diag[0] - tmp1;
2813 lmax = diag[0] + tmp1;
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>
2837 for (
int row = 0; row < this->
nrows(); row++) {
2845 template<
class Treal,
class Telement>
2855 for (
int rc = 0; rc < this->
ncols(); rc++) {
2862 template<
class Treal,
class Telement>
2868 for (
int ind = 0; ind < this->
nElements(); ind++)
2873 template<
class Treal,
class Telement>
2877 for (
int ind = 0; ind < this->
nElements(); ind++)
2882 template<
class Treal,
class Telement>
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>
2905 for (
int row = 0; row <
col; row++)
2908 for (
int rc = 0; rc < this->
nrows(); rc++)
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;
3306 return (
size_t)this->
nElements() *
sizeof(Treal);
3331 int n = this->
nrows();
3332 return this->
is_zero() ? 0 : n * (n + 1) / 2;
3342 int rowOffset = this->
rows.getOffset();
3343 int colOffset = this->
cols.getOffset();
3345 for (
int row = 0; row <
col; row++) {
3346 res += 2*op.accumulate((*
this)(row,
col),
3350 res += op.accumulate((*
this)(
col,
col),
3361 int rowOffset = this->
rows.getOffset();
3362 int colOffset = this->
cols.getOffset();
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>
3455 int nTotalRows = this->
rows.getNTotalScalars();
3456 int nTotalCols = this->
cols.getNTotalScalars();
3457 assert((
int)fullMat.size() == nTotalRows * nTotalCols);
3458 int rowOffset = this->
rows.getOffset();
3459 int colOffset = this->
cols.getOffset();
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 {
3471 int nTotalRows = this->
rows.getNTotalScalars();
3472 int nTotalCols = this->
cols.getNTotalScalars();
3473 fullMat.resize(nTotalRows * nTotalCols);
3474 int rowOffset = this->
rows.getOffset();
3475 int colOffset = this->
cols.getOffset();
3479 fullMat[(colOffset +
col) * nTotalRows + (rowOffset + row)] = 0;
3483 for (
int row = 0; row < this->
nrows(); row++)
3484 fullMat[(colOffset +
col) * nTotalRows + (rowOffset + row)] =
3489 template<
class Treal>
3492 int nTotalRows = this->
rows.getNTotalScalars();
3493 int nTotalCols = this->
cols.getNTotalScalars();
3494 fullMat.resize(nTotalRows * nTotalCols);
3495 int rowOffset = this->
rows.getOffset();
3496 int colOffset = this->
cols.getOffset();
3499 for (
int row = 0; row <=
col; row++) {
3500 fullMat[(colOffset +
col) * nTotalRows + (rowOffset + row)] = 0;
3501 fullMat[(rowOffset + row) * nTotalRows + (colOffset +
col)] = 0;
3506 for (
int row = 0; row <=
col; row++) {
3507 fullMat[(colOffset +
col) * nTotalRows + (rowOffset + row)] =
3509 fullMat[(rowOffset + row) * nTotalRows + (colOffset +
col)] =
3515 template<
class Treal>
3518 int nTotalRows = this->
rows.getNTotalScalars();
3519 int nTotalCols = this->
cols.getNTotalScalars();
3520 fullMat.resize(nTotalRows * nTotalCols);
3521 int rowOffset = this->
rows.getOffset();
3522 int colOffset = this->
cols.getOffset();
3525 for (
int row = 0; row <= this->
nScalarsRows(); row++) {
3526 fullMat[(colOffset +
col) * nTotalRows + (rowOffset + row)] = 0;
3527 fullMat[(rowOffset + row) * nTotalRows + (colOffset +
col)] = 0;
3532 for (
int row = 0; row < this->
nrows(); row++) {
3533 fullMat[(colOffset +
col) * nTotalRows + (rowOffset + row)] =
3535 fullMat[(rowOffset + row) * nTotalRows + (colOffset +
col)] =
3543 template<
class Treal>
3546 std::vector<int>
const & colind,
3547 std::vector<Treal>
const & values) {
3548 assert(rowind.size() == colind.size() &&
3549 rowind.size() == values.size());
3550 std::vector<int> indexes(values.size());
3551 for (
int ind = 0; ind < values.size(); ++ind) {
3557 template<
class Treal>
3560 std::vector<int>
const & colind,
3561 std::vector<Treal>
const & values,
3562 std::vector<int>
const & indexes) {
3563 if (indexes.empty()) {
3569 for (
int ind = 0; ind < this->
nElements(); ++ind)
3571 std::vector<int>::const_iterator it;
3572 for ( it = indexes.begin(); it < indexes.end(); it++ ) {
3573 (*this)(this->
rows.whichBlock( rowind[*it] ),
3574 this->
cols.whichBlock( colind[*it] )) = values[*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())
3603 std::vector<int>::const_iterator it;
3604 for ( it = indexes.begin(); it < indexes.end(); it++ ) {
3605 (*this)(this->
rows.whichBlock( rowind[*it] ),
3606 this->
cols.whichBlock( colind[*it] )) += values[*it];
3610 template<
class Treal>
3613 std::vector<int>
const & colind,
3614 std::vector<Treal>
const & values) {
3615 assert(rowind.size() == colind.size() &&
3616 rowind.size() == values.size());
3617 bool upperTriangleOnly =
true;
3618 for (
int ind = 0; ind < values.size(); ++ind) {
3620 upperTriangleOnly && colind[ind] >= rowind[ind];
3622 if (!upperTriangleOnly)
3623 throw Failure(
"Matrix<Treal>::"
3624 "syAddValues(std::vector<int> const &, "
3625 "std::vector<int> const &, "
3626 "std::vector<Treal> const &, int const): "
3627 "Only upper triangle can contain elements when assigning "
3628 "symmetric or triangular matrix ");
3632 template<
class Treal>
3635 std::vector<int>
const & colind,
3636 std::vector<Treal>
const & values) {
3637 assert(rowind.size() == colind.size() &&
3638 rowind.size() == values.size());
3639 bool upperTriangleOnly =
true;
3640 for (
int ind = 0; ind < values.size(); ++ind) {
3642 upperTriangleOnly && colind[ind] >= rowind[ind];
3644 if (!upperTriangleOnly)
3645 throw Failure(
"Matrix<Treal>::"
3646 "syAddValues(std::vector<int> const &, "
3647 "std::vector<int> const &, "
3648 "std::vector<Treal> const &, int const): "
3649 "Only upper triangle can contain elements when assigning "
3650 "symmetric or triangular matrix ");
3654 template<
class Treal>
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 {
3675 if (indexes.empty())
3677 std::vector<int>::const_iterator it;
3679 for ( it = indexes.begin(); it < indexes.end(); it++ )
3683 for ( it = indexes.begin(); it < indexes.end(); it++ )
3684 values[*it] = (*this)(this->
rows.whichBlock( rowind[*it] ),
3685 this->
cols.whichBlock( colind[*it] ));
3689 template<
class Treal>
3692 std::vector<int>
const & colind,
3693 std::vector<Treal> & values)
const {
3694 assert(rowind.size() == colind.size());
3695 bool upperTriangleOnly =
true;
3696 for (
int ind = 0; ind < rowind.size(); ++ind) {
3698 upperTriangleOnly && colind[ind] >= rowind[ind];
3700 if (!upperTriangleOnly)
3701 throw Failure(
"Matrix<Treal>::"
3702 "syGetValues(std::vector<int> const &, "
3703 "std::vector<int> const &, "
3704 "std::vector<Treal> const &, int const): "
3705 "Only upper triangle when retrieving elements from "
3706 "symmetric or triangular matrix ");
3710 template<
class Treal>
3713 std::vector<int> & colind,
3714 std::vector<Treal> & values)
const {
3715 assert(rowind.size() == colind.size() &&
3716 rowind.size() == values.size());
3719 for (
int row = 0; row < this->
nrows(); row++) {
3720 rowind.push_back(this->
rows.getOffset() + row);
3721 colind.push_back(this->
cols.getOffset() +
col);
3722 values.push_back((*
this)(row,
col));
3726 template<
class Treal>
3729 std::vector<int> & colind,
3730 std::vector<Treal> & values)
const {
3731 assert(rowind.size() == colind.size() &&
3732 rowind.size() == values.size());
3735 for (
int row = 0; row <=
col; ++row) {
3736 rowind.push_back(this->
rows.getOffset() + row);
3737 colind.push_back(this->
cols.getOffset() +
col);
3738 values.push_back((*
this)(row,
col));
3743 template<
class Treal>
3749 template<
class Treal>
3755 char * tmp = (
char*)&
ZERO;
3756 file.write(tmp,
sizeof(
int));
3759 char * tmp = (
char*)&
ONE;
3760 file.write(tmp,
sizeof(
int));
3761 char * tmpel = (
char*)this->
elements;
3762 file.write(tmpel,
sizeof(Treal) * this->
nElements());
3766 template<
class Treal>
3771 char tmp[
sizeof(int)];
3772 file.read(tmp, (std::ifstream::pos_type)
sizeof(
int));
3773 switch ((
int)*tmp) {
3783 throw Failure(
"Matrix<Treal>::"
3784 "readFromFile(std::ifstream & file):"
3785 "File corruption, int value not 0 or 1");
3789 template<
class Treal>
3793 for (
int ind = 0; ind < this->
nElements(); ind++)
3794 this->
elements[ind] = rand() / (Treal)RAND_MAX;
3797 template<
class Treal>
3803 for (
int row = 0; row <
col; row++)
3804 (*
this)(row,
col) = rand() / (Treal)RAND_MAX;;
3806 for (
int rc = 0; rc < this->
nrows(); rc++)
3807 (*
this)(rc,rc) = rand() / (Treal)RAND_MAX;;
3810 template<
class Treal>
3814 probabilityBeingZero > rand() / (Treal)RAND_MAX)
3820 template<
class Treal>
3824 probabilityBeingZero > rand() / (Treal)RAND_MAX)
3830 template<
class Treal>
3831 template<
typename TRule>
3837 for (
int row = 0; row < this->
nrows(); row++)
3838 (*
this)(row,
col) = rule.set(this->rows.getOffset() + row,
3839 this->cols.getOffset() +
col);
3841 template<
class Treal>
3842 template<
typename TRule>
3849 for (
int row = 0; row <=
col; row++)
3850 (*
this)(row,
col) = rule.set(this->rows.getOffset() + row,
3851 this->cols.getOffset() +
col);
3855 template<
class Treal>
3859 throw Failure(
"Matrix<Treal>::addIdentity(Treal): "
3860 "Cannot add identity to empty matrix.");
3862 throw Failure(
"Matrix<Treal, Telement>::addIdentity(Treal): "
3863 "Matrix must be square to add identity");
3866 for (
int ind = 0; ind < this->
ncols(); ind++)
3867 (*
this)(ind,ind) += alpha;
3870 template<
class Treal>
3886 for (
int row = 0; row < AT.
nrows(); row++)
3891 template<
class Treal>
3898 for (
int row = 1; row < this->
nrows(); row++)
3900 (*
this)(row,
col) = (*
this)(
col, row);
3904 throw Failure(
"Matrix<Treal>::symToNosym(): "
3905 "Only quadratic matrices can be symmetric");
3908 template<
class Treal>
3916 for (
int row =
col + 1; row < this->
nrows(); row++)
3917 (*
this)(row,
col) = 0;
3921 throw Failure(
"Matrix<Treal>::nosymToSym(): "
3922 "Only quadratic matrices can be symmetric");
3925 template<
class Treal>
3934 throw Failure(
"Matrix<Treal>::operator=(int k = 1): "
3935 "Matrix must be quadratic to "
3936 "become identity matrix.");
3939 for (
int rc = 0; rc < this->
ncols(); rc++)
3944 (
"Matrix<Treal>::operator=(int k) only implemented for k = 0, k = 1");
3949 template<
class Treal>
3952 if (!this->
is_zero() && alpha != 1) {
3953 for (
int ind = 0; ind < this->
nElements(); ind++)
3959 template<
class Treal>
3961 gemm(
const bool tA,
const bool tB,
const Treal alpha,
3980 if ((
A.is_zero() ||
B.is_zero() || alpha == 0) &&
3984 Treal beta_tmp = beta;
3990 if (!
A.is_zero() && !
B.is_zero() && alpha != 0) {
3992 if (
A.ncols() ==
B.nrows() &&
3993 A.nrows() == C.
nrows() &&
3995 mat::gemm(
"N",
"N", &
A.nrows(), &
B.ncols(), &
A.ncols(), &alpha,
3996 A.elements, &
A.nrows(),
B.elements, &
B.nrows(),
3999 throw Failure(
"Matrix<Treal>::"
4000 "gemm(bool, bool, Treal, Matrix<Treal>, "
4001 "Matrix<Treal>, Treal, "
4003 "Incorrect matrixdimensions for multiplication");
4005 if (
A.nrows() ==
B.nrows() &&
4006 A.ncols() == C.
nrows() &&
4008 mat::gemm(
"T",
"N", &
A.ncols(), &
B.ncols(), &
A.nrows(), &alpha,
4009 A.elements, &
A.nrows(),
B.elements, &
B.nrows(),
4012 throw Failure(
"Matrix<Treal>::"
4013 "gemm(bool, bool, Treal, Matrix<Treal>, "
4014 "Matrix<Treal>, Treal, "
4016 "Incorrect matrixdimensions for multiplication");
4018 if (
A.ncols() ==
B.ncols() &&
4019 A.nrows() == C.
nrows() &&
4021 mat::gemm(
"N",
"T", &
A.nrows(), &
B.nrows(), &
A.ncols(), &alpha,
4022 A.elements, &
A.nrows(),
B.elements, &
B.nrows(),
4025 throw Failure(
"Matrix<Treal>::"
4026 "gemm(bool, bool, Treal, Matrix<Treal>, "
4027 "Matrix<Treal>, Treal, "
4029 "Incorrect matrixdimensions for multiplication");
4031 if (
A.nrows() ==
B.ncols() &&
4032 A.ncols() == C.
nrows() &&
4034 mat::gemm(
"T",
"T", &
A.ncols(), &
B.nrows(), &
A.nrows(), &alpha,
4035 A.elements, &
A.nrows(),
B.elements, &
B.nrows(),
4038 throw Failure(
"Matrix<Treal>::"
4039 "gemm(bool, bool, Treal, Matrix<Treal>, "
4040 "Matrix<Treal>, Treal, "
4042 "Incorrect matrixdimensions for multiplication");
4043 else throw Failure(
"Matrix<Treal>::"
4044 "gemm(bool, bool, Treal, Matrix<Treal>, "
4045 "Matrix<Treal>, Treal, "
4046 "Matrix<Treal>):Very strange error!!");
4054 template<
class Treal>
4056 symm(
const char side,
const char uplo,
const Treal alpha,
4078 if ((
A.is_zero() ||
B.is_zero() || alpha == 0) &&
4082 Treal beta_tmp = beta;
4087 if (!
A.is_zero() && !
B.is_zero() && alpha != 0) {
4088 if (
A.nrows() ==
A.ncols() && C.
nrows() ==
B.nrows() &&
4089 C.
ncols() ==
B.ncols() &&
4090 ((
side ==
'L' &&
A.ncols() ==
B.nrows()) ||
4091 (
side ==
'R' &&
B.ncols() ==
A.nrows())))
4093 A.elements, &
A.nrows(),
B.elements, &
B.nrows(),
4096 throw Failure(
"Matrix<Treal>::symm: Incorrect matrix "
4097 "dimensions for symmetric multiply");
4104 template<
class Treal>
4106 syrk(
const char uplo,
const bool tA,
const Treal alpha,
4123 ((!tA &&
A.nrows() == C.
nrows()) || (tA &&
A.ncols() == C.
nrows())))
4124 if (alpha != 0 && !
A.is_zero()) {
4125 Treal beta_tmp = beta;
4132 A.elements, &
A.nrows(), &beta_tmp,
4137 A.elements, &
A.nrows(), &beta_tmp,
4143 throw Failure(
"Matrix<Treal>::syrk: Incorrect matrix "
4144 "dimensions for symmetric rank-k update");
4147 template<
class Treal>
4149 sysq(
const char uplo,
const Treal alpha,
4166 template<
class Treal>
4168 ssmm(
const Treal alpha,
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,
4210 if (alpha != 0 && !
A.is_zero() && !
B.is_zero())
4211 if (((
side ==
'R' &&
B.ncols() ==
A.nrows()) ||
4212 (
side ==
'L' &&
A.ncols() ==
B.nrows())) &&
4213 A.nrows() ==
A.ncols())
4216 &alpha,
A.elements, &
A.nrows(),
B.elements, &
B.nrows());
4219 &alpha,
A.elements, &
A.nrows(),
B.elements, &
B.nrows());
4221 throw Failure(
"Matrix<Treal>::trmm: "
4222 "Incorrect matrix dimensions for multiplication");
4227 template<
class Treal>
4234 for (
int i = 0; i < this->
nElements(); i++)
4240 template<
class Treal>
4245 throw Failure(
"Matrix<Treal>::syFrobSquared(): "
4246 "Matrix must be have equal number of rows "
4247 "and cols to be symmetric");
4251 for (
int row = 0; row <
col; row++)
4252 sum += 2 * (*
this)(row,
col) * (*
this)(row,
col);
4253 for (
int rc = 0; rc < this->ncols(); rc++)
4254 sum += (*
this)(rc, rc) * (*
this)(rc, rc);
4259 template<
class Treal>
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>
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()) {
4302 for (
int row = 0; row <
col; row++) {
4304 sum += 2 * diff * diff;
4306 for (
int rc = 0; rc <
A.ncols(); rc++) {
4307 diff =
A(rc, rc) -
B(rc, rc);
4311 else if (
A.is_zero() && !
B.is_zero())
4312 sum =
B.syFrobSquared();
4313 else if (!
A.is_zero() &&
B.is_zero())
4314 sum =
A.syFrobSquared();
4319 template<
class Treal>
4324 throw Failure(
"Matrix<Treal>::trace(): Matrix must be quadratic");
4329 for (
int rc = 0; rc < this->
ncols(); rc++)
4330 sum += (*
this)(rc,rc);
4335 template<
class Treal>
4341 if (
A.nrows() !=
B.ncols() ||
A.ncols() !=
B.nrows())
4342 throw Failure(
"Matrix<Treal>::trace_ab: "
4343 "Wrong matrix dimensions for trace(A * B)");
4345 if (!
A.is_zero() && !
B.is_zero())
4346 for (
int rc = 0; rc <
A.nrows(); rc++)
4347 for (
int colA = 0; colA <
A.ncols(); colA++)
4348 tr +=
A(rc,colA) *
B(colA, rc);
4352 template<
class Treal>
4358 if (
A.ncols() !=
B.ncols() ||
A.nrows() !=
B.nrows())
4359 throw Failure(
"Matrix<Treal>::trace_aTb: "
4360 "Wrong matrix dimensions for trace(A' * B)");
4362 if (!
A.is_zero() && !
B.is_zero())
4363 for (
int rc = 0; rc <
A.ncols(); rc++)
4364 for (
int rowB = 0; rowB <
B.nrows(); rowB++)
4365 tr +=
A(rowB,rc) *
B(rowB, rc);
4369 template<
class Treal>
4375 if (
A.nrows() !=
B.ncols() ||
A.ncols() !=
B.nrows() ||
4376 A.nrows() !=
A.ncols())
4377 throw Failure(
"Matrix<Treal>::sy_trace_ab: "
4378 "Wrong matrix dimensions for symmetric trace(A * B)");
4379 if (
A.is_zero() ||
B.is_zero())
4384 for (
int rc = 0; rc <
A.nrows(); rc++)
4385 tr +=
A(rc,rc) *
B(rc, rc);
4387 for (
int colA = 1; colA <
A.ncols(); colA++)
4388 for (
int rowA = 0; rowA < colA; rowA++)
4389 tr += 2 *
A(rowA, colA) *
B(rowA, colA);
4393 template<
class Treal>
4395 add(
const Treal alpha,
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,
4426 template<
class Treal>
4433 template<
class Treal>
4437 for (
int ind = 0; ind < this->
nElements(); ind++)
4443 template<
class Treal>
4457 template<
class Treal>
4462 bool thisMatIsZero =
true;
4465 bool EMatIsZero =
true;
4466 if (!ErrorMatrix->
is_zero() || !this->is_zero()) {
4471 for (
int ind = 0; ind < this->
nElements(); ind++) {
4482 thisMatIsZero =
false;
4485 if ( ErrorMatrix->
elements[ind] != 0 ) {
4488 if ( ErrorMatrix->
elements[ind] * ErrorMatrix->
elements[ind] > threshold ) {
4492 thisMatIsZero =
false;
4498 if (thisMatIsZero) {
4500 for (
int ind = 0; ind < this->nElements(); ind++)
4507 for (
int ind = 0; ind < this->nElements(); ind++)
4510 ErrorMatrix->
clear();
4516 for (
int ind = 0; ind < this->
nElements(); ind++) {
4522 thisMatIsZero =
false;
4535 template<
class Treal>
4553 template<
class Treal>
4566 template<
class Treal>
4569 const bool tA,
const Treal alpha,
4583 template<
class Treal>
4604 template<
class Treal>
4608 if (ErrorMatrix && ErrorMatrix->
is_empty()) {
4613 if (fs < threshold) {
4623 template<
class Treal>
4627 const Treal threshold,
const side looking,
4630 if (
A.nrows() !=
A.ncols())
4631 throw Failure(
"Matrix<Treal>::inch: Matrix must be quadratic!");
4633 throw Failure(
"Matrix<Treal>::inch: Matrix is zero! "
4634 "Not possible to compute inverse cholesky.");
4639 throw Failure(
"Matrix<Treal>::inch: potrf returned info > 0. The matrix is not positive definite.");
4641 throw Failure(
"Matrix<Treal>::inch: potrf returned info < 0");
4645 throw Failure(
"Matrix<Treal>::inch: trtri returned info > 0. The matrix is not positive definite.");
4647 throw Failure(
"Matrix<Treal>::inch: trtri returned info < 0");
4652 template<
class Treal>
4658 Treal* colsums =
new Treal[size];
4659 Treal* diag =
new Treal[size];
4660 for (
int ind = 0; ind < size; ind++)
4666 lmin = diag[0] - tmp1;
4667 lmax = diag[0] + tmp1;
4670 tmp2 = diag[
col] - tmp1;
4671 lmin = (tmp2 < lmin ? tmp2 : lmin);
4672 tmp2 = diag[
col] + tmp1;
4673 lmax = (tmp2 > lmax ? tmp2 : lmax);
4679 throw Failure(
"Matrix<Treal>::gershgorin(Treal, Treal): Matrix must be quadratic");
4683 template<
class Treal>
4689 for (
int row = 0; row < this->
nrows(); row++)
4693 template<
class Treal>
4702 for (
int rc = 0; rc < this->
ncols(); rc++)
4703 diag[rc] = (*
this)(rc,rc);
4730 template<
class Treal>
4732 if (this->is_zero() && this->cap > 0) {
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>
4756 if (this->is_zero()) {
4761 vb_length += this->nel;
4765 template<
class Treal>
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>
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.
bool highestLevel() const
Definition MatrixHierarchicBase.h:133
const int & ncols() const
Definition MatrixHierarchicBase.h:71
static void swap(MatrixHierarchicBase< Treal, Treal > &A, MatrixHierarchicBase< Treal, Treal > &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
MatrixHierarchicBase()
Definition MatrixHierarchicBase.h:149
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
Treal frobSquared() const
Definition Matrix.h:4228
void syRandom()
Definition Matrix.h:3798
size_t sy_nnz() const
Returns number of nonzeros in matrix including lower triangle elements.
Definition Matrix.h:3315
void add_abs_col_sums(Treal *abscolsums) const
Definition Matrix.h:4685
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
Vector< Treal, Treal > VectorType
Definition Matrix.h:2928
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
void clear()
Set matrix to zero and delete all arrays.
Definition Matrix.h:3744
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
void addValues(std::vector< int > const &rowind, std::vector< int > const &colind, std::vector< Treal > const &values)
Definition Matrix.h:3581
void random()
Definition Matrix.h:3790
void getValues(std::vector< int > const &rowind, std::vector< int > const &colind, std::vector< Treal > &values) const
Definition Matrix.h:3656
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
static void ssmm(const Treal alpha, const Matrix< Treal > &A, const Matrix< Treal > &B, const Treal beta, Matrix< Treal > &C)
Definition Matrix.h:4168
void sy_gershgorin(Treal &lmin, Treal &lmax) const
Definition Matrix.h:3292
void assignFromSparse(std::vector< int > const &rowind, std::vector< int > const &colind, std::vector< Treal > const &values)
Definition Matrix.h:3545
void syAssignFrobNormsLowestLevel(Matrix< Treal, Matrix< Treal > > const &A)
Definition Matrix.h:3148
static void inch(const Matrix< Treal > &A, Matrix< Treal > &Z, const Treal threshold=0, const side looking=left, const inchversion version=unstable)
Definition Matrix.h:4625
void syAssignDiffFrobNormsLowestLevel(Matrix< Treal, Matrix< Treal > > const &A, Matrix< Treal, Matrix< Treal > > const &B)
Definition Matrix.h:3195
void get_diagonal(Treal *diag) const
Definition Matrix.h:4695
void truncateAccordingToSparsityPattern(Matrix< Treal, Matrix< Treal > > &A) const
Definition Matrix.h:3235
Matrix class and heart of the matrix library.
Definition Matrix.h:115
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 addValues(std::vector< int > const &rowind, std::vector< int > const &colind, std::vector< Treal > const &values, std::vector< int > const &indexes)
Definition Matrix.h:653
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
Telement ElementType
Definition Matrix.h:117
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
Matrix< Treal, Telement > & operator=(int const k)
Definition Matrix.h:1048
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 assignFromSparse(std::vector< int > const &rowind, std::vector< int > const &colind, std::vector< Treal > const &values, std::vector< int > const &indexes)
Definition Matrix.h:601
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
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 getValues(std::vector< int > const &rowind, std::vector< int > const &colind, std::vector< Treal > &, std::vector< int > const &indexes) const
Definition Matrix.h:746
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
Vector< real, typename ElementType::VectorType > VectorType
Definition Matrix.h:118
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
Definition matInclude.h:156
void tic()
Definition matInclude.cc:84
float toc()
Definition matInclude.cc:88
Vector class.
Definition Vector.h:54
#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
const Treal Matrix< Treal >::ONE
Definition Matrix.h:3448
@ 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
const Treal Matrix< Treal >::ZERO
Definition Matrix.h:3446
A quicksort implementation.
Treal template_blas_sqrt(Treal x)
Treal template_blas_fabs(Treal x)