75 template<
class Treal,
class Telement>
77 template<
typename Tperm>
112 template<
class Treal,
class Telement = Treal>
126 #ifdef MAT_USE_ALLOC_TIMER
130 #ifdef MAT_USE_ALLOC_TIMER
148 void fullMatrix(std::vector<Treal> & fullMat)
const;
154 std::vector<int>
const & colind,
155 std::vector<Treal>
const & values);
157 std::vector<int>
const & colind,
158 std::vector<Treal>
const & values,
159 std::vector<int>
const & indexes);
161 void addValues(std::vector<int>
const & rowind,
162 std::vector<int>
const & colind,
163 std::vector<Treal>
const & values);
164 void addValues(std::vector<int>
const & rowind,
165 std::vector<int>
const & colind,
166 std::vector<Treal>
const & values,
167 std::vector<int>
const & indexes);
170 std::vector<int>
const & colind,
171 std::vector<Treal>
const & values);
174 std::vector<int>
const & colind,
175 std::vector<Treal>
const & values);
177 void getValues(std::vector<int>
const & rowind,
178 std::vector<int>
const & colind,
179 std::vector<Treal> & values)
const;
180 void getValues(std::vector<int>
const & rowind,
181 std::vector<int>
const & colind,
182 std::vector<Treal> &,
183 std::vector<int>
const & indexes)
const;
185 std::vector<int>
const & colind,
186 std::vector<Treal> & values)
const;
188 std::vector<int> & colind,
189 std::vector<Treal> &)
const;
191 std::vector<int> & colind,
192 std::vector<Treal> &)
const;
217 template<
typename TRule>
219 template<
typename TRule>
221 template<
typename TRule>
242 static void gemm(
const bool tA,
const bool tB,
const Treal alpha,
247 static void symm(
const char side,
const char uplo,
const Treal alpha,
252 static void syrk(
const char uplo,
const bool tA,
const Treal alpha,
257 static void sysq(
const char uplo,
const Treal alpha,
262 static void ssmm(
const Treal alpha,
276 static void trmm(
const char side,
const char uplo,
const bool tA,
321 static void add(
const Treal alpha,
324 void assign(Treal
const alpha,
341 (Treal
const threshold,
351 (
Matrix<Treal, Matrix<Treal, Telement> >
const & A );
354 (
Matrix<Treal, Matrix<Treal, Telement> >
const & A );
360 (
Matrix<Treal, Matrix<Treal, Telement> >
const & A,
361 Matrix<Treal, Matrix<Treal, Telement> >
const & B );
366 (
Matrix<Treal, Matrix<Treal, Telement> >
const & A,
367 Matrix<Treal, Matrix<Treal, Telement> >
const & B );
378 const Matrix<Treal, Telement>& A,
379 const Matrix<Treal, Telement>& B,
381 Matrix<Treal, Telement>& C);
383 Matrix<Treal, Telement>& A,
384 const Matrix<Treal, Telement>& Z);
386 const bool tA,
const Treal alpha,
387 const Matrix<Treal, Telement>& A,
388 Matrix<Treal, Telement>& B);
390 const Matrix<Treal, Telement>& Z,
391 Matrix<Treal, Telement>& A);
394 (Treal
const threshold,
405 (Treal
const threshold,
414 const Treal threshold = 0,
418 void gersgorin(Treal& lmin, Treal& lmax)
const;
454 template<
typename Top>
458 for (
int col = 0; col < this->
ncols(); col++) {
459 for (
int row = 0; row < col; row++)
467 template<
typename Top>
471 for (
int col = 0; col < this->
ncols(); col++)
472 for (
int row = 0; row < this->
nrows(); row++)
478 static inline unsigned int level() {
return Telement::level() + 1;}
484 Treal maxAbsGlobal = 0;
485 Treal maxAbsLocal = 0;
486 for (
int ind = 0; ind < this->
nElements(); ++ind) {
487 maxAbsLocal = this->
elements[ind].maxAbsValue();
488 maxAbsGlobal = maxAbsGlobal > maxAbsLocal ?
489 maxAbsGlobal : maxAbsLocal;
502 template<
class Treal,
class Telement>
505 int nTotalRows = this->rows.getNTotalScalars();
506 int nTotalCols = this->cols.getNTotalScalars();
507 assert((
int)fullMat.size() == nTotalRows * nTotalCols);
510 for (
int col = 0; col < this->ncols(); col++)
511 for (
int row = 0; row < this->nrows(); row++)
512 (*
this)(row, col).assignFromFull(fullMat);
515 template<
class Treal,
class Telement>
518 int nTotalRows = this->rows.getNTotalScalars();
519 int nTotalCols = this->cols.getNTotalScalars();
520 fullMat.resize(nTotalRows * nTotalCols);
521 if (this->is_zero()) {
522 int rowOffset = this->rows.getOffset();
523 int colOffset = this->cols.getOffset();
524 for (
int col = 0; col < this->nScalarsCols(); col++)
525 for (
int row = 0; row < this->nScalarsRows(); row++)
526 fullMat[(colOffset + col) * nTotalRows + (rowOffset + row)] = 0;
529 for (
int col = 0; col < this->ncols(); col++)
530 for (
int row = 0; row < this->nrows(); row++)
531 (*
this)(row, col).fullMatrix(fullMat);
535 template<
class Treal,
class Telement>
538 int nTotalRows = this->rows.getNTotalScalars();
539 int nTotalCols = this->cols.getNTotalScalars();
540 fullMat.resize(nTotalRows * nTotalCols);
541 if (this->is_zero()) {
542 int rowOffset = this->rows.getOffset();
543 int colOffset = this->cols.getOffset();
544 for (
int col = 0; col < this->nScalarsCols(); col++)
545 for (
int row = 0; row <= col; row++) {
546 fullMat[(colOffset + col) * nTotalRows + (rowOffset + row)] = 0;
547 fullMat[(rowOffset + row) * nTotalRows + (colOffset + col)] = 0;
551 for (
int col = 0; col < this->ncols(); col++) {
552 for (
int row = 0; row < col; row++)
553 (*
this)(row, col).syUpTriFullMatrix(fullMat);
554 (*this)(col, col).syFullMatrix(fullMat);
559 template<
class Treal,
class Telement>
562 int nTotalRows = this->rows.getNTotalScalars();
563 int nTotalCols = this->cols.getNTotalScalars();
564 fullMat.resize(nTotalRows * nTotalCols);
565 if (this->is_zero()) {
566 int rowOffset = this->rows.getOffset();
567 int colOffset = this->cols.getOffset();
568 for (
int col = 0; col < this->nScalarsCols(); col++)
569 for (
int row = 0; row <= this->nScalarsRows(); row++) {
570 fullMat[(colOffset + col) * nTotalRows + (rowOffset + row)] = 0;
571 fullMat[(rowOffset + row) * nTotalRows + (colOffset + col)] = 0;
575 for (
int col = 0; col < this->ncols(); col++)
576 for (
int row = 0; row < this->nrows(); row++)
577 (*
this)(row, col).syUpTriFullMatrix(fullMat);
584 template<
class Treal,
class Telement>
587 std::vector<int>
const & colind,
588 std::vector<Treal>
const & values) {
589 assert(rowind.size() == colind.size() &&
590 rowind.size() == values.size());
591 std::vector<int> indexes(values.size());
592 for (
unsigned int ind = 0; ind < values.size(); ++ind)
594 assignFromSparse(rowind, colind, values, indexes);
597 template<
class Treal,
class Telement>
600 std::vector<int>
const & colind,
601 std::vector<Treal>
const & values,
602 std::vector<int>
const & indexes) {
603 if (indexes.empty()) {
610 std::vector<std::vector<int> > columnBuckets (this->cols.getNBlocks());
611 std::vector<std::vector<int> > rowBuckets (this->rows.getNBlocks());
615 std::vector<int>::const_iterator it;
616 for ( it = indexes.begin(); it < indexes.end(); it++ )
617 columnBuckets[ this->cols.whichBlock( colind[*it] ) ].push_back (*it);
620 for (
int col = 0; col < this->ncols(); col++) {
622 while (!columnBuckets[col].empty()) {
623 currentInd = columnBuckets[col].back();
624 columnBuckets[col].pop_back();
625 rowBuckets[ this->rows.whichBlock
626 ( rowind[currentInd] ) ].push_back (currentInd);
629 for (
int row = 0; row < this->nrows(); row++) {
630 (*this)(row,col).assignFromSparse(rowind, colind, values, rowBuckets[row]);
631 rowBuckets[row].clear();
636 template<
class Treal,
class Telement>
639 std::vector<int>
const & colind,
640 std::vector<Treal>
const & values) {
641 assert(rowind.size() == colind.size() &&
642 rowind.size() == values.size());
643 std::vector<int> indexes(values.size());
644 for (
unsigned int ind = 0; ind < values.size(); ++ind)
646 addValues(rowind, colind, values, indexes);
649 template<
class Treal,
class Telement>
652 std::vector<int>
const & colind,
653 std::vector<Treal>
const & values,
654 std::vector<int>
const & indexes) {
660 std::vector<std::vector<int> > columnBuckets (this->cols.getNBlocks());
661 std::vector<std::vector<int> > rowBuckets (this->rows.getNBlocks());
664 std::vector<int>::const_iterator it;
665 for ( it = indexes.begin(); it < indexes.end(); it++ )
666 columnBuckets[ this->cols.whichBlock( colind[*it] ) ].push_back (*it);
669 for (
int col = 0; col < this->ncols(); col++) {
671 while (!columnBuckets[col].empty()) {
672 currentInd = columnBuckets[col].back();
673 columnBuckets[col].pop_back();
674 rowBuckets[ this->rows.whichBlock
675 ( rowind[currentInd] ) ].push_back (currentInd);
678 for (
int row = 0; row < this->nrows(); row++) {
679 (*this)(row,col).addValues(rowind, colind, values, rowBuckets[row]);
680 rowBuckets[row].clear();
685 template<
class Treal,
class Telement>
688 std::vector<int>
const & colind,
689 std::vector<Treal>
const & values) {
690 assert(rowind.size() == colind.size() &&
691 rowind.size() == values.size());
692 bool upperTriangleOnly =
true;
693 for (
unsigned int ind = 0; ind < values.size(); ++ind) {
695 upperTriangleOnly && colind[ind] >= rowind[ind];
697 if (!upperTriangleOnly)
698 throw Failure(
"Matrix<Treal, Telement>::"
699 "syAddValues(std::vector<int> const &, "
700 "std::vector<int> const &, "
701 "std::vector<Treal> const &, int const): "
702 "Only upper triangle can contain elements when assigning "
703 "symmetric or triangular matrix ");
704 assignFromSparse(rowind, colind, values);
707 template<
class Treal,
class Telement>
710 std::vector<int>
const & colind,
711 std::vector<Treal>
const & values) {
712 assert(rowind.size() == colind.size() &&
713 rowind.size() == values.size());
714 bool upperTriangleOnly =
true;
715 for (
unsigned int ind = 0; ind < values.size(); ++ind) {
717 upperTriangleOnly && colind[ind] >= rowind[ind];
719 if (!upperTriangleOnly)
720 throw Failure(
"Matrix<Treal, Telement>::"
721 "syAddValues(std::vector<int> const &, "
722 "std::vector<int> const &, "
723 "std::vector<Treal> const &, int const): "
724 "Only upper triangle can contain elements when assigning "
725 "symmetric or triangular matrix ");
726 addValues(rowind, colind, values);
729 template<
class Treal,
class Telement>
732 std::vector<int>
const & colind,
733 std::vector<Treal> & values)
const {
734 assert(rowind.size() == colind.size());
735 values.resize(rowind.size(),0);
736 std::vector<int> indexes(rowind.size());
737 for (
unsigned int ind = 0; ind < rowind.size(); ++ind)
739 getValues(rowind, colind, values, indexes);
742 template<
class Treal,
class Telement>
745 std::vector<int>
const & colind,
746 std::vector<Treal> & values,
747 std::vector<int>
const & indexes)
const {
748 assert(!this->is_empty());
751 std::vector<int>::const_iterator it;
752 if (this->is_zero()) {
753 for ( it = indexes.begin(); it < indexes.end(); it++ )
758 std::vector<std::vector<int> > columnBuckets (this->cols.getNBlocks());
759 std::vector<std::vector<int> > rowBuckets (this->rows.getNBlocks());
762 for ( it = indexes.begin(); it < indexes.end(); it++ )
763 columnBuckets[ this->cols.whichBlock( colind[*it] ) ].push_back (*it);
766 for (
int col = 0; col < this->ncols(); col++) {
768 while (!columnBuckets[col].empty()) {
769 currentInd = columnBuckets[col].back();
770 columnBuckets[col].pop_back();
771 rowBuckets[ this->rows.whichBlock
772 ( rowind[currentInd] ) ].push_back (currentInd);
775 for (
int row = 0; row < this->nrows(); row++) {
776 (*this)(row,col).getValues(rowind, colind, values, rowBuckets[row]);
777 rowBuckets[row].clear();
782 template<
class Treal,
class Telement>
785 std::vector<int>
const & colind,
786 std::vector<Treal> & values)
const {
787 assert(rowind.size() == colind.size());
788 bool upperTriangleOnly =
true;
789 for (
int unsigned ind = 0; ind < rowind.size(); ++ind) {
791 upperTriangleOnly && colind[ind] >= rowind[ind];
793 if (!upperTriangleOnly)
794 throw Failure(
"Matrix<Treal, Telement>::"
795 "syGetValues(std::vector<int> const &, "
796 "std::vector<int> const &, "
797 "std::vector<Treal> const &, int const): "
798 "Only upper triangle when retrieving elements from "
799 "symmetric or triangular matrix ");
800 getValues(rowind, colind, values);
804 template<
class Treal,
class Telement>
807 std::vector<int> & colind,
808 std::vector<Treal> & values)
const {
809 assert(rowind.size() == colind.size() &&
810 rowind.size() == values.size());
811 if (!this->is_zero())
812 for (
int ind = 0; ind < this->nElements(); ++ind)
813 this->elements[ind].getAllValues(rowind, colind, values);
816 template<
class Treal,
class Telement>
819 std::vector<int> & colind,
820 std::vector<Treal> & values)
const {
821 assert(rowind.size() == colind.size() &&
822 rowind.size() == values.size());
823 if (!this->is_zero())
824 for (
int col = 0; col < this->ncols(); ++col) {
825 for (
int row = 0; row < col; ++row)
826 (*
this)(row, col).getAllValues(rowind, colind, values);
827 (*this)(col, col).syGetAllValues(rowind, colind, values);
832 template<
class Treal,
class Telement>
834 if (!this->is_zero())
835 for (
int i = 0; i < this->nElements(); i++)
836 this->elements[i].clear();
841 template<
class Treal,
class Telement>
846 if (this->is_zero()) {
847 char * tmp = (
char*)&ZERO;
848 file.write(tmp,
sizeof(
int));
851 char * tmp = (
char*)&ONE;
852 file.write(tmp,
sizeof(
int));
853 for (
int i = 0; i < this->nElements(); i++)
854 this->elements[i].writeToFile(file);
857 template<
class Treal,
class Telement>
862 char tmp[
sizeof(int)];
863 file.read(tmp, (std::ifstream::pos_type)
sizeof(
int));
871 for (
int i = 0; i < this->nElements(); i++)
872 this->elements[i].readFromFile(file);
875 throw Failure(
"Matrix<Treal, Telement>::"
876 "readFromFile(std::ifstream & file):"
877 "File corruption int value not 0 or 1");
881 template<
class Treal,
class Telement>
885 for (
int ind = 0; ind < this->nElements(); ind++)
886 this->elements[ind].random();
889 template<
class Treal,
class Telement>
894 for (
int col = 1; col < this->ncols(); col++)
895 for (
int row = 0; row < col; row++)
896 (*
this)(row, col).random();
898 for (
int rc = 0; rc < this->nrows(); rc++)
899 (*
this)(rc,rc).syRandom();
902 template<
class Treal,
class Telement>
905 if (!this->highestLevel() &&
906 probabilityBeingZero > rand() / (Treal)RAND_MAX)
911 for (
int ind = 0; ind < this->nElements(); ind++)
912 this->elements[ind].randomZeroStructure(probabilityBeingZero);
916 template<
class Treal,
class Telement>
919 if (!this->highestLevel() &&
920 probabilityBeingZero > rand() / (Treal)RAND_MAX)
926 for (
int col = 1; col < this->ncols(); col++)
927 for (
int row = 0; row < col; row++)
928 (*
this)(row, col).randomZeroStructure(probabilityBeingZero);
930 for (
int rc = 0; rc < this->nrows(); rc++)
931 (*
this)(rc,rc).syRandomZeroStructure(probabilityBeingZero);
935 template<
class Treal,
class Telement>
936 template<
typename TRule>
941 for (
int ind = 0; ind < this->nElements(); ind++)
942 this->elements[ind].setElementsByRule(rule);
944 template<
class Treal,
class Telement>
945 template<
typename TRule>
951 for (
int col = 1; col < this->ncols(); col++)
952 for (
int row = 0; row < col; row++)
953 (*
this)(row, col).setElementsByRule(rule);
955 for (
int rc = 0; rc < this->nrows(); rc++)
956 (*
this)(rc,rc).sySetElementsByRule(rule);
960 template<
class Treal,
class Telement>
963 if (this->is_empty())
964 throw Failure(
"Matrix<Treal, Telement>::addIdentity(Treal): "
965 "Cannot add identity to empty matrix.");
966 if (this->ncols() != this->nrows())
967 throw Failure(
"Matrix<Treal, Telement>::addIdentity(Treal): "
968 "Matrix must be square to add identity");
971 for (
int ind = 0; ind < this->ncols(); ind++)
972 (*
this)(ind,ind).addIdentity(alpha);
975 template<
class Treal,
class Telement>
992 for (
int row = 0; row < AT.
nrows(); row++)
993 for (
int col = 0; col < AT.
ncols(); col++)
997 template<
class Treal,
class Telement>
1001 if (this->nrows() == this->ncols()) {
1002 if (!this->is_zero()) {
1004 for (
int rc = 0; rc < this->ncols(); rc++)
1005 (*
this)(rc, rc).symToNosym();
1007 for (
int row = 1; row < this->nrows(); row++)
1008 for (
int col = 0; col < row; col++)
1013 throw Failure(
"Matrix<Treal, Telement>::symToNosym(): "
1014 "Only quadratic matrices can be symmetric");
1017 std::cout<<
"Failure caught:Matrix<Treal, Telement>::symToNosym()"
1023 template<
class Treal,
class Telement>
1026 if (this->nrows() == this->ncols()) {
1027 if (!this->is_zero()) {
1029 for (
int rc = 0; rc < this->ncols(); rc++)
1030 (*
this)(rc, rc).nosymToSym();
1032 for (
int col = 0; col < this->ncols() - 1; col++)
1033 for (
int row = col + 1; row < this->nrows(); row++)
1034 (*
this)(row, col).clear();
1038 throw Failure(
"Matrix<Treal, Telement>::nosymToSym(): "
1039 "Only quadratic matrices can be symmetric");
1044 template<
class Treal,
class Telement>
1052 if (this->ncols() != this->nrows())
1054 (
"Matrix::operator=(int k = 1): "
1055 "Matrix must be quadratic to become identity matrix.");
1058 for (
int rc = 0; rc < this->ncols(); rc++)
1062 throw Failure(
"Matrix::operator=(int k) only "
1063 "implemented for k = 0, k = 1");
1069 template<
class Treal,
class Telement>
1072 if (!this->is_zero() && alpha != 1) {
1073 for (
int ind = 0; ind < this->nElements(); ind++)
1074 this->elements[ind] *= alpha;
1080 template<
class Treal,
class Telement>
1082 gemm(
const bool tA,
const bool tB,
const Treal alpha,
1105 Treal beta_tmp = beta;
1117 #pragma omp parallel for if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) schedule(dynamic)
1119 for (
int col = 0; col < C.
ncols(); col++) {
1121 for (
int row = 0; row < C.
nrows(); row++) {
1123 A(row, 0),
B(0, col),
1126 for(
int cola = 1; cola < A.
ncols(); cola++)
1128 A(row, cola),
B(cola, col),
1136 throw Failure(
"Matrix<class Treal, class Telement>::"
1137 "gemm(bool, bool, Treal, "
1138 "Matrix<Treal, Telement>, "
1139 "Matrix<Treal, Telement>, Treal, "
1140 "Matrix<Treal, Telement>): "
1141 "Incorrect matrixdimensions for multiplication");
1147 #pragma omp parallel for if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) schedule(dynamic)
1149 for (
int col = 0; col < C.
ncols(); col++) {
1151 for (
int row = 0; row < C.
nrows(); row++) {
1156 for(
int cola = 1; cola < A.
nrows(); cola++)
1158 A(cola, row),
B(cola, col),
1166 throw Failure(
"Matrix<class Treal, class Telement>::"
1167 "gemm(bool, bool, Treal, "
1168 "Matrix<Treal, Telement>, "
1169 "Matrix<Treal, Telement>, Treal, "
1170 "Matrix<Treal, Telement>): "
1171 "Incorrect matrixdimensions for multiplication");
1177 #pragma omp parallel for if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) schedule(dynamic)
1179 for (
int col = 0; col < C.
ncols(); col++) {
1181 for (
int row = 0; row < C.
nrows(); row++) {
1183 A(row, 0),
B(col, 0),
1186 for(
int cola = 1; cola < A.
ncols(); cola++)
1188 A(row, cola),
B(col, cola),
1196 throw Failure(
"Matrix<class Treal, class Telement>::"
1197 "gemm(bool, bool, Treal, "
1198 "Matrix<Treal, Telement>, "
1199 "Matrix<Treal, Telement>, Treal, "
1200 "Matrix<Treal, Telement>): "
1201 "Incorrect matrixdimensions for multiplication");
1207 #pragma omp parallel for if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) schedule(dynamic)
1209 for (
int col = 0; col < C.
ncols(); col++) {
1211 for (
int row = 0; row < C.
nrows(); row++) {
1213 A(0, row),
B(col, 0),
1216 for(
int cola = 1; cola < A.
nrows(); cola++)
1218 A(cola, row),
B(col, cola),
1226 throw Failure(
"Matrix<class Treal, class Telement>::"
1227 "gemm(bool, bool, Treal, "
1228 "Matrix<Treal, Telement>, "
1229 "Matrix<Treal, Telement>, "
1230 "Treal, Matrix<Treal, Telement>): "
1231 "Incorrect matrixdimensions for multiplication");
1232 else throw Failure(
"Matrix<class Treal, class Telement>::"
1233 "gemm(bool, bool, Treal, "
1234 "Matrix<Treal, Telement>, "
1235 "Matrix<Treal, Telement>, Treal, "
1236 "Matrix<Treal, Telement>):"
1237 "Very strange error!!");
1245 template<
class Treal,
class Telement>
1247 symm(
const char side,
const char uplo,
const Treal alpha,
1255 int dimA = A.
nrows();
1263 assert(side ==
'R');
1273 Treal beta_tmp = beta;
1281 if (dimA == B.
nrows() &&
1282 dimA == C.
nrows() &&
1285 #pragma omp parallel for if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) schedule(dynamic)
1287 for (
int col = 0; col < C.
ncols(); col++) {
1289 for (
int row = 0; row < C.
nrows(); row++) {
1292 B(row, col), beta_tmp, C(row, col));
1294 for(
int ind = 0; ind < row; ind++)
1296 A(ind, row),
B(ind, col),
1299 for(
int ind = row + 1; ind < dimA; ind++)
1301 A(row, ind),
B(ind, col),
1308 throw Failure(
"Matrix<class Treal, class Telement>"
1309 "::symm(bool, bool, Treal, Matrix<Treal, "
1310 "Telement>, Matrix<Treal, Telement>, "
1311 "Treal, Matrix<Treal, Telement>): "
1312 "Incorrect matrixdimensions for multiplication");
1314 assert(side ==
'R');
1315 if (B.
ncols() == dimA &&
1317 dimA == C.
ncols()) {
1319 #pragma omp parallel for if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) schedule(dynamic)
1321 for (
int col = 0; col < C.
ncols(); col++) {
1323 for (
int row = 0; row < C.
nrows(); row++) {
1326 B(row, col), beta_tmp, C(row, col));
1328 for(
int ind = col + 1; ind < dimA; ind++)
1330 B(row, ind),
A(col, ind),
1333 for(
int ind = 0; ind < col; ind++)
1335 B(row, ind),
A(ind, col),
1342 throw Failure(
"Matrix<class Treal, class Telement>"
1343 "::symm(bool, bool, Treal, Matrix<Treal, "
1344 "Telement>, Matrix<Treal, Telement>, "
1345 "Treal, Matrix<Treal, Telement>): "
1346 "Incorrect matrixdimensions for multiplication");
1354 throw Failure(
"Matrix<class Treal, class Telement>::"
1355 "symm only implemented for symmetric matrices in "
1356 "upper triangular storage");
1360 template<
class Treal,
class Telement>
1362 syrk(
const char uplo,
const bool tA,
const Treal alpha,
1381 if (alpha != 0 && !A.
is_zero()) {
1382 Treal beta_tmp = beta;
1388 if (!tA && uplo ==
'U') {
1390 #pragma omp parallel if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared)
1394 #pragma omp for schedule(dynamic) nowait
1396 for (
int rc = 0; rc < C.
ncols(); rc++) {
1399 for (
int cola = 1; cola < A.
ncols(); cola++)
1404 #pragma omp for schedule(dynamic) nowait
1406 for (
int row = 0; row < C.
nrows() - 1; row++) {
1408 for (
int col = row + 1; col < C.
ncols(); col++) {
1410 A(row, 0),
A(col,0), beta_tmp, C(row,col));
1411 for (
int ind = 1; ind < A.
ncols(); ind++)
1413 A(row, ind),
A(col,ind), 1.0, C(row,col));
1419 else if (tA && uplo ==
'U') {
1421 #pragma omp parallel if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared)
1425 #pragma omp for schedule(dynamic) nowait
1427 for (
int rc = 0; rc < C.
ncols(); rc++) {
1430 for (
int rowa = 1; rowa < A.
nrows(); rowa++)
1435 #pragma omp for schedule(dynamic) nowait
1437 for (
int row = 0; row < C.
nrows() - 1; row++) {
1439 for (
int col = row + 1; col < C.
ncols(); col++) {
1441 A(0, row),
A(0, col), beta_tmp, C(row,col));
1442 for (
int ind = 1; ind < A.
nrows(); ind++)
1444 A(ind, row),
A(ind, col), 1.0, C(row,col));
1451 throw Failure(
"Matrix<class Treal, class Telement>::"
1452 "syrk not implemented for lower triangular storage");
1459 throw Failure(
"Matrix<class Treal, class Telement>::syrk: "
1460 "Incorrect matrix dimensions for symmetric rank-k update");
1463 template<
class Treal,
class Telement>
1465 sysq(
const char uplo,
const Treal alpha,
1477 if (alpha != 0 && !A.
is_zero()) {
1479 Treal beta_tmp = beta;
1486 #pragma omp parallel if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared)
1490 #pragma omp for schedule(dynamic) nowait
1492 for (
int rc = 0; rc < C.
ncols(); rc++) {
1494 Telement::sysq(uplo, alpha,
A(rc, rc), beta_tmp, C(rc, rc));
1495 for (
int cola = 0; cola < rc; cola++)
1497 for (
int cola = rc + 1; cola < A.
ncols(); cola++)
1503 #pragma omp for schedule(dynamic) nowait
1505 for (
int row = 0; row < C.
nrows() - 1; row++) {
1507 for (
int col = row + 1; col < C.
ncols(); col++) {
1510 beta_tmp, C(row, col));
1514 for (
int ind = 0; ind < row; ind++)
1516 A(ind, row),
A(ind, col), 1.0, C(row, col));
1518 for (
int ind = row + 1; ind < col; ind++)
1520 A(row, ind),
A(ind, col), 1.0, C(row, col));
1522 for (
int ind = col + 1; ind < A.
ncols(); ind++)
1524 A(row, ind),
A(col, ind), 1.0, C(row, col));
1532 throw Failure(
"Matrix<class Treal, class Telement>::"
1533 "sysq only implemented for symmetric matrices in "
1534 "upper triangular storage");;
1540 throw Failure(
"Matrix<class Treal, class Telement>::sysq: "
1541 "Incorrect matrix dimensions for symmetric square "
1546 template<
class Treal,
class Telement>
1565 throw Failure(
"Matrix<class Treal, class Telement>::"
1567 "Matrix<Treal, Telement>, "
1568 "Matrix<Treal, Telement>, Treal, "
1569 "Matrix<Treal, Telement>): "
1570 "Incorrect matrixdimensions for ssmm multiplication");
1577 Treal beta_tmp = beta;
1588 #pragma omp parallel for if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) schedule(dynamic)
1590 for (
int col = 0; col < C.
ncols(); col++) {
1594 Telement::ssmm(alpha,
A(col,col),
B(col, col),
1595 beta_tmp, C(col,col));
1596 for (
int ind = 0; ind < col; ind++)
1599 alpha,
A(ind,col),
B(ind, col),
1601 for (
int ind = col + 1; ind < A.
ncols(); ind++)
1604 alpha,
A(col, ind),
B(col, ind),
1607 for (
int row = 0; row < col; row++) {
1612 alpha,
A(row, row),
B(row, col),
1613 beta_tmp, C(row, col));
1616 alpha,
B(col, col),
A(row, col),
1618 for (
int ind = 0; ind < row; ind++)
1621 alpha,
A(ind, row),
B(ind, col),
1623 for (
int ind = row + 1; ind < col; ind++)
1626 alpha,
A(row, ind),
B(ind, col),
1628 for (
int ind = col + 1; ind < A.
ncols(); ind++)
1631 alpha,
A(row, ind),
B(col, ind),
1636 for (
int row = col + 1; row < C.
nrows(); row++) {
1642 alpha,
B(col, col),
A(col, row),
1643 beta_tmp, tmpSubMat);
1646 alpha,
A(row, row),
B(col, row),
1648 for (
int ind = 0; ind < col; ind++)
1651 alpha,
B(ind, col),
A(ind, row),
1653 for (
int ind = col + 1; ind < row; ind++)
1656 alpha,
B(col, ind),
A(ind, row),
1658 for (
int ind = row + 1; ind < B.
nrows(); ind++)
1661 alpha,
B(col, ind),
A(row, ind),
1674 template<
class Treal,
class Telement>
1693 throw Failure(
"Matrix<class Treal, class Telement>::"
1694 "ssmm_upper_tr_only(Treal, "
1695 "Matrix<Treal, Telement>, "
1696 "Matrix<Treal, Telement>, Treal, "
1697 "Matrix<Treal, Telement>): "
1698 "Incorrect matrixdimensions for ssmm_upper_tr_only "
1706 Treal beta_tmp = beta;
1717 #pragma omp parallel for if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) schedule(dynamic)
1719 for (
int col = 0; col < C.
ncols(); col++) {
1723 Telement::ssmm_upper_tr_only(alpha,
A(col,col),
B(col, col),
1724 beta_tmp, C(col,col));
1725 for (
int ind = 0; ind < col; ind++)
1727 Telement::gemm_upper_tr_only(
true,
false,
1728 alpha,
A(ind,col),
B(ind, col),
1730 for (
int ind = col + 1; ind < A.
ncols(); ind++)
1732 Telement::gemm_upper_tr_only(
false,
true,
1733 alpha,
A(col, ind),
B(col, ind),
1736 for (
int row = 0; row < col; row++) {
1741 alpha,
A(row, row),
B(row, col),
1742 beta_tmp, C(row, col));
1745 alpha,
B(col, col),
A(row, col),
1747 for (
int ind = 0; ind < row; ind++)
1750 alpha,
A(ind, row),
B(ind, col),
1752 for (
int ind = row + 1; ind < col; ind++)
1755 alpha,
A(row, ind),
B(ind, col),
1757 for (
int ind = col + 1; ind < A.
ncols(); ind++)
1760 alpha,
A(row, ind),
B(col, ind),
1770 template<
class Treal,
class Telement>
1772 trmm(
const char side,
const char uplo,
const bool tA,
const Treal alpha,
1778 if (((side ==
'R' && B.
ncols() == A.
nrows()) ||
1785 for (
int col = B.
ncols() - 1; col >= 0; col--)
1786 for (
int row = 0; row < B.
nrows(); row++) {
1790 A(col, col),
B(row,col));
1792 for (
int ind = 0; ind < col; ind++)
1794 B(row,ind),
A(ind, col),
1799 assert(side ==
'L');
1801 for (
int row = 0; row < B.
nrows(); row++ )
1802 for (
int col = 0; col < B.
ncols(); col++) {
1804 A(row,row),
B(row,col));
1805 for (
int ind = row + 1 ; ind < B.
nrows(); ind++)
1807 A(row,ind),
B(ind,col),
1816 for (
int col = 0; col < B.
ncols(); col++)
1817 for (
int row = 0; row < B.
nrows(); row++) {
1819 A(col,col),
B(row,col));
1820 for (
int ind = col + 1; ind < A.
ncols(); ind++)
1822 B(row,ind),
A(col,ind),
1827 assert(side ==
'L');
1829 for (
int row = B.
nrows() - 1; row >= 0; row--)
1830 for (
int col = 0; col < B.
ncols(); col++) {
1832 A(row,row),
B(row,col));
1833 for (
int ind = 0; ind < row; ind++)
1835 A(ind,row),
B(ind,col),
1841 throw Failure(
"Matrix<class Treal, class Telement>::"
1842 "trmm not implemented for lower triangular matrices");
1844 throw Failure(
"Matrix<class Treal, class Telement>::trmm"
1845 ": Incorrect matrix dimensions for multiplication");
1851 template<
class Treal,
class Telement>
1853 assert(!this->is_empty());
1854 if (this->is_zero())
1858 for (
int i = 0; i < this->nElements(); i++)
1859 sum += this->elements[i].frobSquared();
1864 template<
class Treal,
class Telement>
1867 assert(!this->is_empty());
1868 if (this->nrows() != this->ncols())
1869 throw Failure(
"Matrix<Treal, Telement>::syFrobSquared(): "
1870 "Matrix must be have equal number of rows "
1871 "and cols to be symmetric");
1873 if (!this->is_zero()) {
1874 for (
int col = 1; col < this->ncols(); col++)
1875 for (
int row = 0; row < col; row++)
1876 sum += 2 * (*
this)(row, col).frobSquared();
1877 for (
int rc = 0; rc < this->ncols(); rc++)
1878 sum += (*
this)(rc, rc).syFrobSquared();
1883 template<
class Treal,
class Telement>
1891 throw Failure(
"Matrix<Treal, Telement>::"
1892 "frobSquaredDiff: Incorrect matrix dimensions");
1905 template<
class Treal,
class Telement>
1915 throw Failure(
"Matrix<Treal, Telement>::syFrobSquaredDiff: "
1916 "Incorrect matrix dimensions");
1919 for (
int col = 1; col < A.
ncols(); col++)
1920 for (
int row = 0; row < col; row++)
1921 sum += 2 * Telement::frobSquaredDiff(
A(row, col),
B(row, col));
1922 for (
int rc = 0; rc < A.
ncols(); rc++)
1923 sum += Telement::syFrobSquaredDiff(
A(rc, rc),
B(rc, rc));
1933 template<
class Treal,
class Telement>
1936 assert(!this->is_empty());
1937 if (this->nrows() != this->ncols())
1938 throw Failure(
"Matrix<Treal, Telement>::trace(): "
1939 "Matrix must be quadratic");
1940 if (this->is_zero())
1944 for (
int rc = 0; rc < this->ncols(); rc++)
1945 sum += (*
this)(rc,rc).
trace();
1950 template<
class Treal,
class Telement>
1957 throw Failure(
"Matrix<Treal, Telement>::trace_ab: "
1958 "Wrong matrix dimensions for trace(A * B)");
1961 for (
int rc = 0; rc < A.
nrows(); rc++)
1962 for (
int colA = 0; colA < A.
ncols(); colA++)
1963 tr += Telement::trace_ab(
A(rc,colA),
B(colA, rc));
1967 template<
class Treal,
class Telement>
1974 throw Failure(
"Matrix<Treal, Telement>::trace_aTb: "
1975 "Wrong matrix dimensions for trace(A' * B)");
1978 for (
int rc = 0; rc < A.
ncols(); rc++)
1979 for (
int rowB = 0; rowB < B.
nrows(); rowB++)
1980 tr += Telement::trace_aTb(
A(rowB,rc),
B(rowB, rc));
1984 template<
class Treal,
class Telement>
1992 throw Failure(
"Matrix<Treal, Telement>::sy_trace_ab: "
1993 "Wrong matrix dimensions for symmetric trace(A * B)");
1997 for (
int rc = 0; rc < A.
nrows(); rc++)
1998 tr += Telement::sy_trace_ab(
A(rc,rc),
B(rc, rc));
2000 for (
int colA = 1; colA < A.
ncols(); colA++)
2001 for (
int rowA = 0; rowA < colA; rowA++)
2002 tr += 2 * Telement::trace_aTb(
A(rowA, colA),
B(rowA, colA));
2007 template<
class Treal,
class Telement>
2015 throw Failure(
"Matrix<Treal, Telement>::add: "
2016 "Wrong matrix dimensions for addition");
2017 if (!A.
is_zero() && alpha != 0) {
2020 for (
int ind = 0; ind < A.
nElements(); ind++)
2025 template<
class Treal,
class Telement>
2030 if (this->is_empty()) {
2031 this->resetRows(A.
rows);
2032 this->resetCols(A.
cols);
2036 add(alpha, A, *
this);
2043 template<
class Treal,
class Telement>
2046 if (!this->is_zero())
2047 for (
int ind = 0; ind < this->nElements(); ind++)
2048 this->elements[ind].getFrobSqLowestLevel(frobsq);
2051 template<
class Treal,
class Telement>
2055 assert(!this->is_empty());
2056 bool thisMatIsZero =
true;
2059 bool EMatIsZero =
true;
2060 if (!ErrorMatrix->
is_zero() || !this->is_zero()) {
2063 if (this->is_zero())
2065 for (
int ind = 0; ind < this->nElements(); ind++) {
2066 this->elements[ind].
2067 frobThreshLowestLevel(threshold, &ErrorMatrix->
elements[ind]);
2068 thisMatIsZero = thisMatIsZero && this->elements[ind].is_zero();
2069 EMatIsZero = EMatIsZero && ErrorMatrix->
elements[ind].is_zero();
2074 ErrorMatrix->
clear();
2078 if (!this->is_zero()) {
2079 for (
int ind = 0; ind < this->nElements(); ind++) {
2080 this->elements[ind].frobThreshLowestLevel(threshold, 0);
2081 thisMatIsZero = thisMatIsZero && this->elements[ind].is_zero();
2088 template<
class Treal,
class Telement>
2091 if (!this->is_zero())
2092 for (
int ind = 0; ind < this->nElements(); ind++)
2093 this->elements[ind].getFrobSqElementLevel(frobsq);
2096 template<
class Treal,
class Telement>
2100 assert(!this->is_empty());
2101 bool thisMatIsZero =
true;
2104 bool EMatIsZero =
true;
2105 if (!ErrorMatrix->
is_zero() || !this->is_zero()) {
2108 if (this->is_zero())
2110 for (
int ind = 0; ind < this->nElements(); ind++) {
2111 this->elements[ind].
2112 frobThreshElementLevel(threshold, &ErrorMatrix->
elements[ind]);
2113 thisMatIsZero = thisMatIsZero && this->elements[ind].is_zero();
2114 EMatIsZero = EMatIsZero && ErrorMatrix->
elements[ind].is_zero();
2119 ErrorMatrix->
clear();
2123 if (!this->is_zero()) {
2124 for (
int ind = 0; ind < this->nElements(); ind++) {
2125 this->elements[ind].frobThreshElementLevel(threshold, 0);
2126 thisMatIsZero = thisMatIsZero && this->elements[ind].is_zero();
2135 template<
class Treal,
class Telement>
2139 if ( this->is_zero() )
2141 assert( this->nElements() == A.
nElements() );
2142 for (
int ind = 0; ind < this->nElements(); ind++)
2143 this->elements[ind].assignFrobNormsLowestLevel(A[ind]);
2149 template<
class Treal,
class Telement>
2154 throw Failure(
"Matrix<Treal, Telement>::syAssignFrobNormsLowestLevel(...): "
2155 "Matrix must be have equal number of rows "
2156 "and cols to be symmetric");
2158 if ( this->is_zero() )
2160 assert( this->nElements() == A.
nElements() );
2161 for (
int col = 1; col < this->ncols(); col++)
2162 for (
int row = 0; row < col; row++)
2163 (*
this)(row, col).assignFrobNormsLowestLevel(
A(row,col) );
2164 for (
int rc = 0; rc < this->ncols(); rc++)
2165 (*
this)(rc, rc).syAssignFrobNormsLowestLevel(
A(rc,rc) );
2171 template<
class Treal,
class Telement>
2181 if ( this->is_zero() )
2185 assert( this->nElements() == A.
nElements() );
2186 assert( this->nElements() == B.
nElements() );
2187 for (
int ind = 0; ind < this->nElements(); ind++)
2188 this->elements[ind].assignDiffFrobNormsLowestLevel( A[ind], B[ind] );
2193 this->assignFrobNormsLowestLevel( A );
2198 this->assignFrobNormsLowestLevel( B );
2202 template<
class Treal,
class Telement>
2212 if ( this->is_zero() )
2216 assert( this->nElements() == A.
nElements() );
2217 assert( this->nElements() == B.
nElements() );
2218 for (
int col = 1; col < this->ncols(); col++)
2219 for (
int row = 0; row < col; row++)
2220 (*
this)(row, col).assignDiffFrobNormsLowestLevel(
A(row,col),
B(row,col) );
2221 for (
int rc = 0; rc < this->ncols(); rc++)
2222 (*
this)(rc, rc).syAssignDiffFrobNormsLowestLevel(
A(rc,rc),
B(rc,rc) );
2227 this->syAssignFrobNormsLowestLevel( A );
2232 this->syAssignFrobNormsLowestLevel( B );
2239 template<
class Treal,
class Telement>
2242 if ( this->is_zero() ) {
2247 assert( this->nElements() == A.
nElements() );
2248 for (
int ind = 0; ind < this->nElements(); ind++)
2249 this->elements[ind].truncateAccordingToSparsityPattern( A[ind] );
2259 template<
class Treal,
class Telement>
2283 Treal beta_tmp = beta;
2293 for (
int col = 0; col < C.
ncols(); col++) {
2294 Telement::gemm_upper_tr_only(tA, tB, alpha,
2295 A(col, 0),
B(0, col),
2298 for(
int cola = 1; cola < A.
ncols(); cola++)
2299 Telement::gemm_upper_tr_only(tA, tB, alpha,
2300 A(col, cola),
B(cola, col),
2303 for (
int row = 0; row < col; row++) {
2305 A(row, 0),
B(0, col),
2308 for(
int cola = 1; cola < A.
ncols(); cola++)
2310 A(row, cola),
B(cola, col),
2317 throw Failure(
"Matrix<class Treal, class Telement>::"
2318 "gemm_upper_tr_only(bool, bool, Treal, "
2319 "Matrix<Treal, Telement>, "
2320 "Matrix<Treal, Telement>, "
2321 "Treal, Matrix<Treal, Telement>): "
2322 "Incorrect matrixdimensions for multiplication");
2327 for (
int col = 0; col < C.
ncols(); col++) {
2328 Telement::gemm_upper_tr_only(tA, tB, alpha,
2332 for(
int cola = 1; cola < A.
nrows(); cola++)
2333 Telement::gemm_upper_tr_only(tA, tB, alpha,
2334 A(cola, col),
B(cola, col),
2337 for (
int row = 0; row < col; row++) {
2342 for(
int cola = 1; cola < A.
nrows(); cola++)
2344 A(cola, row),
B(cola, col),
2351 throw Failure(
"Matrix<class Treal, class Telement>::"
2352 "gemm_upper_tr_only(bool, bool, Treal, "
2353 "Matrix<Treal, Telement>, "
2354 "Matrix<Treal, Telement>, "
2355 "Treal, Matrix<Treal, Telement>): "
2356 "Incorrect matrixdimensions for multiplication");
2361 for (
int col = 0; col < C.
ncols(); col++) {
2362 Telement::gemm_upper_tr_only(tA, tB, alpha,
2363 A(col, 0),
B(col, 0),
2366 for(
int cola = 1; cola < A.
ncols(); cola++)
2367 Telement::gemm_upper_tr_only(tA, tB, alpha,
2368 A(col, cola),
B(col, cola),
2371 for (
int row = 0; row < col; row++) {
2373 A(row, 0),
B(col, 0),
2376 for(
int cola = 1; cola < A.
ncols(); cola++)
2378 A(row, cola),
B(col, cola),
2385 throw Failure(
"Matrix<class Treal, class Telement>::"
2386 "gemm_upper_tr_only(bool, bool, Treal, "
2387 "Matrix<Treal, Telement>, "
2388 "Matrix<Treal, Telement>, "
2389 "Treal, Matrix<Treal, Telement>): "
2390 "Incorrect matrixdimensions for multiplication");
2395 for (
int col = 0; col < C.
ncols(); col++) {
2396 Telement::gemm_upper_tr_only(tA, tB, alpha,
2397 A(0, col),
B(col, 0),
2400 for(
int cola = 1; cola < A.
nrows(); cola++)
2401 Telement::gemm_upper_tr_only(tA, tB, alpha,
2402 A(cola, col),
B(col, cola),
2405 for (
int row = 0; row < col; row++) {
2407 A(0, row),
B(col, 0),
2410 for(
int cola = 1; cola < A.
nrows(); cola++)
2412 A(cola, row),
B(col, cola),
2419 throw Failure(
"Matrix<class Treal, class Telement>::"
2420 "gemm_upper_tr_only(bool, bool, Treal, "
2421 "Matrix<Treal, Telement>, "
2422 "Matrix<Treal, Telement>, Treal, "
2423 "Matrix<Treal, Telement>): "
2424 "Incorrect matrixdimensions for multiplication");
2425 else throw Failure(
"Matrix<class Treal, class Telement>::"
2426 "gemm_upper_tr_only(bool, bool, Treal, "
2427 "Matrix<Treal, Telement>, "
2428 "Matrix<Treal, Telement>, Treal, "
2429 "Matrix<Treal, Telement>):"
2430 "Very strange error!!");
2442 template<
class Treal,
class Telement>
2455 for (
int col = A.
ncols() - 1; col >= 0; col--) {
2457 Telement::sytr_upper_tr_only(side, alpha,
2458 A(col, col), Z(col, col));
2459 for (
int ind = 0; ind < col; ind++) {
2461 Telement::gemm_upper_tr_only(
true,
false, alpha,
A(ind, col),
2462 Z(ind, col), 1.0,
A(col, col));
2465 for (
int row = col - 1; row >= 0; row--) {
2468 alpha, Z(col, col),
A(row, col));
2472 for (
int ind = 0; ind < row; ind++) {
2477 for (
int ind = row + 1; ind < col; ind++) {
2486 assert(side ==
'L');
2488 for (
int col = 0; col < A.
ncols(); col++) {
2490 for (
int row = 0; row < col; row++) {
2493 Z(row, row),
A(row, col));
2497 for (
int ind = row + 1; ind < col; ind++)
2501 for (
int ind = col + 1; ind < A.
nrows(); ind++)
2506 Telement::sytr_upper_tr_only(side, alpha,
2507 A(col, col), Z(col, col));
2508 for (
int ind = col + 1; ind < A.
ncols(); ind++)
2509 Telement::gemm_upper_tr_only(
false,
true, alpha, Z(col, ind),
2510 A(col, ind), 1.0,
A(col, col));
2515 throw Failure(
"Matrix<class Treal, class Telement>::"
2516 "sytr_upper_tr_only: Incorrect matrix dimensions "
2517 "for symmetric-triangular multiplication");
2525 template<
class Treal,
class Telement>
2528 const bool tA,
const Treal alpha,
2534 if (((side ==
'R' && B.
ncols() == A.
nrows()) ||
2539 throw Failure(
"Matrix<Treal, class Telement>::"
2540 "trmm_upper_tr_only : "
2541 "not possible for upper triangular not transposed "
2542 "matrices A since lower triangle of B is needed");
2548 for (
int col = 0; col < B.
ncols(); col++) {
2549 Telement::trmm_upper_tr_only(side, uplo, tA, alpha,
2550 A(col,col),
B(col,col));
2551 for (
int ind = col + 1; ind < A.
ncols(); ind++)
2552 Telement::gemm_upper_tr_only(
false, tA, alpha,
2553 B(col,ind),
A(col,ind),
2555 for (
int row = 0; row < col; row++) {
2557 A(col,col),
B(row,col));
2558 for (
int ind = col + 1; ind < A.
ncols(); ind++)
2560 B(row,ind),
A(col,ind),
2566 assert(side ==
'L');
2568 for (
int row = B.
nrows() - 1; row >= 0; row--) {
2569 Telement::trmm_upper_tr_only(side, uplo, tA, alpha,
2570 A(row,row),
B(row,row));
2571 for (
int ind = 0; ind < row; ind++)
2572 Telement::gemm_upper_tr_only(tA,
false, alpha,
2573 A(ind,row),
B(ind,row),
2575 for (
int col = row + 1; col < B.
ncols(); col++) {
2577 A(row,row),
B(row,col));
2578 for (
int ind = 0; ind < row; ind++)
2580 A(ind,row),
B(ind,col),
2587 throw Failure(
"Matrix<class Treal, class Telement>::"
2588 "trmm_upper_tr_only not implemented for lower "
2589 "triangular matrices");
2591 throw Failure(
"Matrix<class Treal, class Telement>::"
2592 "trmm_upper_tr_only: Incorrect matrix dimensions "
2593 "for multiplication");
2602 template<
class Treal,
class Telement>
2614 assert(side ==
'L');
2622 template<
class Treal,
class Telement>
2626 assert(!this->is_empty());
2627 if (ErrorMatrix && ErrorMatrix->
is_empty()) {
2631 assert(threshold >= (Treal)0.0);
2632 if (threshold == (Treal)0.0)
2635 std::vector<Treal> frobsq_norms;
2636 getFrobSqLowestLevel(frobsq_norms);
2638 std::sort(frobsq_norms.begin(), frobsq_norms.end());
2639 int lastRemoved = -1;
2640 Treal frobsqSum = 0;
2641 int nnzBlocks = frobsq_norms.size();
2642 while(lastRemoved < nnzBlocks - 1 && frobsqSum < threshold) {
2644 frobsqSum += frobsq_norms[lastRemoved];
2648 if (lastRemoved == nnzBlocks - 1 && frobsqSum < threshold) {
2655 frobsqSum -= frobsq_norms[lastRemoved];
2657 while(lastRemoved > -1 && frobsq_norms[lastRemoved] ==
2658 frobsq_norms[lastRemoved + 1]) {
2659 frobsqSum -= frobsq_norms[lastRemoved];
2662 if (lastRemoved > -1) {
2663 Treal threshLowestLevel =
2664 (frobsq_norms[lastRemoved + 1] +
2665 frobsq_norms[lastRemoved]) / 2;
2666 this->frobThreshLowestLevel(threshLowestLevel, ErrorMatrix);
2672 template<
class Treal,
class Telement>
2676 const Treal threshold,
const side looking,
2680 throw Failure(
"Matrix<Treal, Telement>::sy_inch: "
2681 "Matrix must be quadratic!");
2683 throw Failure(
"Matrix<Treal>::sy_inch: Matrix is zero! "
2684 "Not possible to compute inverse cholesky.");
2686 throw Failure(
"Matrix<Treal>::sy_inch: Only unstable "
2687 "version of sy_inch implemented.");
2690 myThresh = threshold / (Z.
ncols() * Z.
nrows());
2695 if (looking ==
left) {
2697 throw Failure(
"Matrix<Treal, Telement>::syInch: "
2698 "Thresholding not implemented for left-looking inch.");
2700 Telement::syInch(
A(0,0), Z(0,0), threshold, looking, version);
2702 for (
int i = 1; i < Z.
ncols(); i++) {
2703 for (
int j = 0; j < i; j++) {
2707 for (
int k = 0; k < j; k++)
2709 A(k,j), Z(k,i), 1.0, Ptmp);
2712 for (
int k = 0; k < j; k++)
2714 Z(k,j), Ptmp, 1.0, Z(k,i));
2717 Telement::add(1.0, Ptmp, Z(j,i));
2720 for (
int k = 0; k < i; k++)
2721 Telement::gemm_upper_tr_only(
true,
false, 1.0,
2722 A(k,i), Z(k,i), 1.0, Ptmp);
2725 Telement::syInch(Ptmp, Z(i,i), threshold, looking, version);
2726 for (
int k = 0; k < i; k++) {
2733 assert(looking ==
right);
2735 Treal newThresh = 0;
2736 if (myThresh != 0 && Z.
ncols() > 1)
2737 newThresh = myThresh * 0.0001;
2739 newThresh = myThresh;
2741 for (
int i = 0; i < Z.
ncols(); i++) {
2744 for (
int k = 0; k < i; k++)
2745 Telement::gemm_upper_tr_only(
true,
false, 1.0,
2746 A(k,i), Z(k,i), 1.0, Ptmp);
2749 Telement::syInch(Ptmp, Z(i,i), newThresh, looking, version);
2750 for (
int k = 0; k < i; k++) {
2755 for (
int j = i + 1; j < Z.
ncols(); j++) {
2760 for (
int k = 0; k < i; k++)
2762 A(k,i), Z(k,j), 1.0, Ptmp);
2765 for (
int k = 0; k < i; k++)
2767 Z(k,i), Ptmp, 1.0, Z(k,j));
2770 Telement::add(1.0, Ptmp, Z(i,j));
2774 if (threshold != 0) {
2775 for (
int k = 0; k < i; k++)
2782 template<
class Treal,
class Telement>
2785 assert(!this->is_empty());
2786 if (this->nScalarsRows() == this->nScalarsCols()) {
2787 int size = this->nScalarsCols();
2788 Treal* colsums =
new Treal[size];
2789 Treal* diag =
new Treal[size];
2790 for (
int ind = 0; ind < size; ind++)
2792 this->add_abs_col_sums(colsums);
2793 this->get_diagonal(diag);
2796 lmin = diag[0] - tmp1;
2797 lmax = diag[0] + tmp1;
2798 for (
int col = 1; col < size; col++) {
2800 tmp2 = diag[col] - tmp1;
2801 lmin = (tmp2 < lmin ? tmp2 : lmin);
2802 tmp2 = diag[col] + tmp1;
2803 lmax = (tmp2 > lmax ? tmp2 : lmax);
2809 throw Failure(
"Matrix<Treal, Telement, int>::gersgorin(Treal, Treal): "
2810 "Matrix must be quadratic");
2814 template<
class Treal,
class Telement>
2818 if (!this->is_zero()) {
2820 for (
int col = 0; col < this->ncols(); col++) {
2821 for (
int row = 0; row < this->nrows(); row++) {
2822 (*this)(row,col).add_abs_col_sums(&sums[offset]);
2824 offset += (*this)(0,col).nScalarsCols();
2829 template<
class Treal,
class Telement>
2833 assert(this->nrows() == this->ncols());
2834 if (this->is_zero())
2835 for (
int rc = 0; rc < this->nScalarsCols(); rc++)
2839 for (
int rc = 0; rc < this->ncols(); rc++) {
2840 (*this)(rc,rc).get_diagonal(&diag[offset]);
2841 offset += (*this)(rc,rc).nScalarsCols();
2846 template<
class Treal,
class Telement>
2848 assert(!this->is_empty());
2850 if (this->is_zero())
2852 for (
int ind = 0; ind < this->nElements(); ind++)
2853 sum += this->elements[ind].memory_usage();
2857 template<
class Treal,
class Telement>
2860 if (!this->is_zero()) {
2861 for (
int ind = 0; ind < this->nElements(); ind++)
2862 sum += this->elements[ind].nnz();
2866 template<
class Treal,
class Telement>
2869 if (!this->is_zero()) {
2871 for (
int col = 1; col < this->ncols(); col++)
2872 for (
int row = 0; row < col; row++)
2873 sum += (*
this)(row, col).nnz();
2877 for (
int rc = 0; rc < this->nrows(); rc++)
2878 sum += (*
this)(rc,rc).sy_nnz();
2883 template<
class Treal,
class Telement>
2886 if (!this->is_zero()) {
2888 for (
int col = 1; col < this->ncols(); col++)
2889 for (
int row = 0; row < col; row++)
2890 sum += (*
this)(row, col).nvalues();
2892 for (
int rc = 0; rc < this->nrows(); rc++)
2893 sum += (*
this)(rc,rc).sy_nvalues();
2908 template<
class Treal>
2928 for (
int ind = 0; ind < this->
nElements(); ++ind)
2934 void fullMatrix(std::vector<Treal> & fullMat)
const;
2940 std::vector<int>
const & colind,
2941 std::vector<Treal>
const & values);
2943 std::vector<int>
const & colind,
2944 std::vector<Treal>
const & values,
2945 std::vector<int>
const & indexes);
2948 void addValues(std::vector<int>
const & rowind,
2949 std::vector<int>
const & colind,
2950 std::vector<Treal>
const & values);
2951 void addValues(std::vector<int>
const & rowind,
2952 std::vector<int>
const & colind,
2953 std::vector<Treal>
const & values,
2954 std::vector<int>
const & indexes);
2957 std::vector<int>
const & colind,
2958 std::vector<Treal>
const & values);
2961 std::vector<int>
const & colind,
2962 std::vector<Treal>
const & values);
2964 void getValues(std::vector<int>
const & rowind,
2965 std::vector<int>
const & colind,
2966 std::vector<Treal> & values)
const;
2967 void getValues(std::vector<int>
const & rowind,
2968 std::vector<int>
const & colind,
2969 std::vector<Treal> & values,
2970 std::vector<int>
const & indexes)
const;
2972 std::vector<int>
const & colind,
2973 std::vector<Treal> & values)
const;
2976 std::vector<int> & colind,
2977 std::vector<Treal> & values)
const;
2979 std::vector<int> & colind,
2980 std::vector<Treal> & values)
const;
3000 template<
typename TRule>
3002 template<
typename TRule>
3019 static void gemm(
const bool tA,
const bool tB,
const Treal alpha,
3024 static void symm(
const char side,
const char uplo,
const Treal alpha,
3029 static void syrk(
const char uplo,
const bool tA,
const Treal alpha,
3034 static void sysq(
const char uplo,
const Treal alpha,
3039 static void ssmm(
const Treal alpha,
3053 static void trmm(
const char side,
const char uplo,
const bool tA,
3084 Treal
trace()
const;
3092 static void add(
const Treal alpha,
3095 void assign(Treal
const alpha,
3112 (Treal
const threshold,
3125 for (
int ind = 0; ind < this->
nElements(); ind++)
3137 for (
int col = 1; col < this->
ncols(); col++)
3138 for (
int row = 0; row < col; row++)
3139 (*
this)(row,col) =
A(row,col).frob();
3140 for (
int rc = 0; rc < this->
nrows(); rc++)
3141 (*
this)(rc,rc) =
A(rc,rc).syFrob();
3161 for (
int ind = 0; ind < this->
nElements(); ind++) {
3193 for (
int col = 1; col < this->
ncols(); col++)
3194 for (
int row = 0; row < col; row++) {
3197 (*this)(row, col) = Diff.
frob();
3199 for (
int rc = 0; rc < this->
ncols(); rc++) {
3202 (*this)(rc, rc) = Diff.
syFrob();
3225 for (
int ind = 0; ind < this->
nElements(); ind++)
3244 const bool tA,
const Treal alpha,
3264 const Treal threshold = 0,
3269 const Treal threshold = 0,
3272 inch(A, Z, threshold, looking, version);
3275 void gersgorin(Treal& lmin, Treal& lmax)
const;
3290 return (
size_t)this->
nElements() *
sizeof(Treal);
3315 int n = this->
nrows();
3316 return this->
is_zero() ? 0 : n * (n + 1) / 2;
3328 for (
int col = 0; col < this->
ncols(); col++) {
3329 for (
int row = 0; row < col; row++) {
3330 res += 2*op.accumulate((*
this)(row, col),
3334 res += op.accumulate((*
this)(col, col),
3347 for (
int col = 0; col < this->
ncols(); col++)
3348 for (
int row = 0; row < this->
nrows(); row++)
3349 res += op.accumulate((*
this)(row, col),
3356 static inline unsigned int level() {
return 0;}
3362 Treal maxAbsGlobal = 0;
3363 Treal maxAbsLocal = 0;
3364 for (
int ind = 0; ind < this->
nElements(); ++ind) {
3366 maxAbsGlobal = maxAbsGlobal > maxAbsLocal ?
3367 maxAbsGlobal : maxAbsLocal;
3369 return maxAbsGlobal;
3381 std::cout<<
"operator="<<std::endl;
3394 void trim_memory_usage();
3396 void write_to_buffer_count(
int& zb_length,
int& vb_length)
const;
3397 void write_to_buffer(
int* zerobuf,
const int zb_length,
3398 Treal* valuebuf,
const int vb_length,
3399 int& zb_index,
int& vb_index)
const;
3400 void read_from_buffer(
int* zerobuf,
const int zb_length,
3401 Treal* valuebuf,
const int vb_length,
3402 int& zb_index,
int& vb_index);
3419 inline bool lowestLevel()
const {
return true;}
3429 template<
class Treal>
3431 template<
class Treal>
3436 template<
class Treal>
3441 assert((
int)fullMat.size() == nTotalRows * nTotalCols);
3446 for (
int col = 0; col < this->
ncols(); col++)
3447 for (
int row = 0; row < this->
nrows(); row++)
3449 fullMat[(colOffset + col) * nTotalRows + (rowOffset + row)];
3452 template<
class Treal>
3457 fullMat.resize(nTotalRows * nTotalCols);
3463 fullMat[(colOffset + col) * nTotalRows + (rowOffset + row)] = 0;
3466 for (
int col = 0; col < this->
ncols(); col++)
3467 for (
int row = 0; row < this->
nrows(); row++)
3468 fullMat[(colOffset + col) * nTotalRows + (rowOffset + row)] =
3473 template<
class Treal>
3478 fullMat.resize(nTotalRows * nTotalCols);
3483 for (
int row = 0; row <= col; row++) {
3484 fullMat[(colOffset + col) * nTotalRows + (rowOffset + row)] = 0;
3485 fullMat[(rowOffset + row) * nTotalRows + (colOffset + col)] = 0;
3489 for (
int col = 0; col < this->
ncols(); col++)
3490 for (
int row = 0; row <= col; row++) {
3491 fullMat[(colOffset + col) * nTotalRows + (rowOffset + row)] =
3493 fullMat[(rowOffset + row) * nTotalRows + (colOffset + col)] =
3499 template<
class Treal>
3504 fullMat.resize(nTotalRows * nTotalCols);
3509 for (
int row = 0; row <= this->
nScalarsRows(); row++) {
3510 fullMat[(colOffset + col) * nTotalRows + (rowOffset + row)] = 0;
3511 fullMat[(rowOffset + row) * nTotalRows + (colOffset + col)] = 0;
3515 for (
int col = 0; col < this->
ncols(); col++)
3516 for (
int row = 0; row < this->
nrows(); row++) {
3517 fullMat[(colOffset + col) * nTotalRows + (rowOffset + row)] =
3519 fullMat[(rowOffset + row) * nTotalRows + (colOffset + col)] =
3527 template<
class Treal>
3530 std::vector<int>
const & colind,
3531 std::vector<Treal>
const & values) {
3532 assert(rowind.size() == colind.size() &&
3533 rowind.size() == values.size());
3534 std::vector<int> indexes(values.size());
3535 for (
int ind = 0; ind < values.size(); ++ind) {
3541 template<
class Treal>
3544 std::vector<int>
const & colind,
3545 std::vector<Treal>
const & values,
3546 std::vector<int>
const & indexes) {
3547 if (indexes.empty()) {
3553 for (
int ind = 0; ind < this->
nElements(); ++ind)
3555 std::vector<int>::const_iterator it;
3556 for ( it = indexes.begin(); it < indexes.end(); it++ ) {
3563 template<
class Treal>
3566 std::vector<int>
const & colind,
3567 std::vector<Treal>
const & values) {
3568 assert(rowind.size() == colind.size() &&
3569 rowind.size() == values.size());
3570 std::vector<int> indexes(values.size());
3571 for (
int ind = 0; ind < values.size(); ++ind) {
3574 addValues(rowind, colind, values, indexes);
3577 template<
class Treal>
3580 std::vector<int>
const & colind,
3581 std::vector<Treal>
const & values,
3582 std::vector<int>
const & indexes) {
3583 if (indexes.empty())
3587 std::vector<int>::const_iterator it;
3588 for ( it = indexes.begin(); it < indexes.end(); it++ ) {
3594 template<
class Treal>
3597 std::vector<int>
const & colind,
3598 std::vector<Treal>
const & values) {
3599 assert(rowind.size() == colind.size() &&
3600 rowind.size() == values.size());
3601 bool upperTriangleOnly =
true;
3602 for (
int ind = 0; ind < values.size(); ++ind) {
3604 upperTriangleOnly && colind[ind] >= rowind[ind];
3606 if (!upperTriangleOnly)
3607 throw Failure(
"Matrix<Treal>::"
3608 "syAddValues(std::vector<int> const &, "
3609 "std::vector<int> const &, "
3610 "std::vector<Treal> const &, int const): "
3611 "Only upper triangle can contain elements when assigning "
3612 "symmetric or triangular matrix ");
3616 template<
class Treal>
3619 std::vector<int>
const & colind,
3620 std::vector<Treal>
const & values) {
3621 assert(rowind.size() == colind.size() &&
3622 rowind.size() == values.size());
3623 bool upperTriangleOnly =
true;
3624 for (
int ind = 0; ind < values.size(); ++ind) {
3626 upperTriangleOnly && colind[ind] >= rowind[ind];
3628 if (!upperTriangleOnly)
3629 throw Failure(
"Matrix<Treal>::"
3630 "syAddValues(std::vector<int> const &, "
3631 "std::vector<int> const &, "
3632 "std::vector<Treal> const &, int const): "
3633 "Only upper triangle can contain elements when assigning "
3634 "symmetric or triangular matrix ");
3638 template<
class Treal>
3641 std::vector<int>
const & colind,
3642 std::vector<Treal> & values)
const {
3643 assert(rowind.size() == colind.size());
3644 values.resize(rowind.size(),0);
3645 std::vector<int> indexes(rowind.size());
3646 for (
int ind = 0; ind < rowind.size(); ++ind) {
3649 getValues(rowind, colind, values, indexes);
3652 template<
class Treal>
3655 std::vector<int>
const & colind,
3656 std::vector<Treal> & values,
3657 std::vector<int>
const & indexes)
const {
3659 if (indexes.empty())
3661 std::vector<int>::const_iterator it;
3663 for ( it = indexes.begin(); it < indexes.end(); it++ )
3667 for ( it = indexes.begin(); it < indexes.end(); it++ )
3673 template<
class Treal>
3676 std::vector<int>
const & colind,
3677 std::vector<Treal> & values)
const {
3678 assert(rowind.size() == colind.size());
3679 bool upperTriangleOnly =
true;
3680 for (
int ind = 0; ind < rowind.size(); ++ind) {
3682 upperTriangleOnly && colind[ind] >= rowind[ind];
3684 if (!upperTriangleOnly)
3685 throw Failure(
"Matrix<Treal>::"
3686 "syGetValues(std::vector<int> const &, "
3687 "std::vector<int> const &, "
3688 "std::vector<Treal> const &, int const): "
3689 "Only upper triangle when retrieving elements from "
3690 "symmetric or triangular matrix ");
3694 template<
class Treal>
3697 std::vector<int> & colind,
3698 std::vector<Treal> & values)
const {
3699 assert(rowind.size() == colind.size() &&
3700 rowind.size() == values.size());
3702 for (
int col = 0; col < this->
ncols(); col++)
3703 for (
int row = 0; row < this->
nrows(); row++) {
3706 values.push_back((*
this)(row, col));
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());
3718 for (
int col = 0; col < this->
ncols(); ++col)
3719 for (
int row = 0; row <= col; ++row) {
3722 values.push_back((*
this)(row, col));
3727 template<
class Treal>
3733 template<
class Treal>
3739 char * tmp = (
char*)&ZERO;
3740 file.write(tmp,
sizeof(
int));
3743 char * tmp = (
char*)&ONE;
3744 file.write(tmp,
sizeof(
int));
3745 char * tmpel = (
char*)this->
elements;
3746 file.write(tmpel,
sizeof(Treal) * this->
nElements());
3750 template<
class Treal>
3755 char tmp[
sizeof(int)];
3756 file.read(tmp, (std::ifstream::pos_type)
sizeof(
int));
3757 switch ((
int)*tmp) {
3767 throw Failure(
"Matrix<Treal>::"
3768 "readFromFile(std::ifstream & file):"
3769 "File corruption, int value not 0 or 1");
3773 template<
class Treal>
3777 for (
int ind = 0; ind < this->
nElements(); ind++)
3778 this->
elements[ind] = rand() / (Treal)RAND_MAX;
3781 template<
class Treal>
3786 for (
int col = 1; col < this->
ncols(); col++)
3787 for (
int row = 0; row < col; row++)
3788 (*
this)(row, col) = rand() / (Treal)RAND_MAX;;
3790 for (
int rc = 0; rc < this->
nrows(); rc++)
3791 (*
this)(rc,rc) = rand() / (Treal)RAND_MAX;;
3794 template<
class Treal>
3798 probabilityBeingZero > rand() / (Treal)RAND_MAX)
3804 template<
class Treal>
3808 probabilityBeingZero > rand() / (Treal)RAND_MAX)
3814 template<
class Treal>
3815 template<
typename TRule>
3820 for (
int col = 0; col < this->
ncols(); col++)
3821 for (
int row = 0; row < this->
nrows(); row++)
3822 (*
this)(row,col) = rule.set(this->rows.getOffset() + row,
3825 template<
class Treal>
3826 template<
typename TRule>
3832 for (
int col = 0; col < this->
ncols(); col++)
3833 for (
int row = 0; row <= col; row++)
3834 (*
this)(row,col) = rule.set(this->rows.getOffset() + row,
3839 template<
class Treal>
3843 throw Failure(
"Matrix<Treal>::addIdentity(Treal): "
3844 "Cannot add identity to empty matrix.");
3846 throw Failure(
"Matrix<Treal, Telement>::addIdentity(Treal): "
3847 "Matrix must be square to add identity");
3850 for (
int ind = 0; ind < this->
ncols(); ind++)
3851 (*
this)(ind,ind) += alpha;
3854 template<
class Treal>
3870 for (
int row = 0; row < AT.
nrows(); row++)
3871 for (
int col = 0; col < AT.
ncols(); col++)
3872 AT(row,col) =
A(col,row);
3875 template<
class Treal>
3882 for (
int row = 1; row < this->
nrows(); row++)
3883 for (
int col = 0; col < row; col++)
3884 (*
this)(row, col) = (*
this)(col, row);
3888 throw Failure(
"Matrix<Treal>::symToNosym(): "
3889 "Only quadratic matrices can be symmetric");
3892 template<
class Treal>
3899 for (
int col = 0; col < this->
ncols() - 1; col++)
3900 for (
int row = col + 1; row < this->
nrows(); row++)
3901 (*
this)(row, col) = 0;
3905 throw Failure(
"Matrix<Treal>::nosymToSym(): "
3906 "Only quadratic matrices can be symmetric");
3909 template<
class Treal>
3918 throw Failure(
"Matrix<Treal>::operator=(int k = 1): "
3919 "Matrix must be quadratic to "
3920 "become identity matrix.");
3923 for (
int rc = 0; rc < this->
ncols(); rc++)
3928 (
"Matrix<Treal>::operator=(int k) only implemented for k = 0, k = 1");
3933 template<
class Treal>
3936 if (!this->
is_zero() && alpha != 1) {
3937 for (
int ind = 0; ind < this->
nElements(); ind++)
3943 template<
class Treal>
3945 gemm(
const bool tA,
const bool tB,
const Treal alpha,
3968 Treal beta_tmp = beta;
3983 throw Failure(
"Matrix<Treal>::"
3984 "gemm(bool, bool, Treal, Matrix<Treal>, "
3985 "Matrix<Treal>, Treal, "
3987 "Incorrect matrixdimensions for multiplication");
3996 throw Failure(
"Matrix<Treal>::"
3997 "gemm(bool, bool, Treal, Matrix<Treal>, "
3998 "Matrix<Treal>, Treal, "
4000 "Incorrect matrixdimensions for multiplication");
4009 throw Failure(
"Matrix<Treal>::"
4010 "gemm(bool, bool, Treal, Matrix<Treal>, "
4011 "Matrix<Treal>, Treal, "
4013 "Incorrect matrixdimensions for multiplication");
4022 throw Failure(
"Matrix<Treal>::"
4023 "gemm(bool, bool, Treal, Matrix<Treal>, "
4024 "Matrix<Treal>, Treal, "
4026 "Incorrect matrixdimensions for multiplication");
4027 else throw Failure(
"Matrix<Treal>::"
4028 "gemm(bool, bool, Treal, Matrix<Treal>, "
4029 "Matrix<Treal>, Treal, "
4030 "Matrix<Treal>):Very strange error!!");
4038 template<
class Treal>
4040 symm(
const char side,
const char uplo,
const Treal alpha,
4066 Treal beta_tmp = beta;
4080 throw Failure(
"Matrix<Treal>::symm: Incorrect matrix "
4081 "dimensions for symmetric multiply");
4088 template<
class Treal>
4090 syrk(
const char uplo,
const bool tA,
const Treal alpha,
4108 if (alpha != 0 && !A.
is_zero()) {
4109 Treal beta_tmp = beta;
4127 throw Failure(
"Matrix<Treal>::syrk: Incorrect matrix "
4128 "dimensions for symmetric rank-k update");
4131 template<
class Treal>
4133 sysq(
const char uplo,
const Treal alpha,
4150 template<
class Treal>
4173 template<
class Treal>
4181 ssmm(alpha, A, B, beta, C);
4186 template<
class Treal>
4195 if (((side ==
'R' && B.
ncols() == A.
nrows()) ||
4205 throw Failure(
"Matrix<Treal>::trmm: "
4206 "Incorrect matrix dimensions for multiplication");
4211 template<
class Treal>
4218 for (
int i = 0; i < this->
nElements(); i++)
4224 template<
class Treal>
4229 throw Failure(
"Matrix<Treal>::syFrobSquared(): "
4230 "Matrix must be have equal number of rows "
4231 "and cols to be symmetric");
4234 for (
int col = 1; col < this->
ncols(); col++)
4235 for (
int row = 0; row < col; row++)
4236 sum += 2 * (*
this)(row, col) * (*
this)(row, col);
4237 for (
int rc = 0; rc < this->
ncols(); rc++)
4238 sum += (*
this)(rc, rc) * (*
this)(rc, rc);
4243 template<
class Treal>
4251 throw Failure(
"Matrix<Treal>::frobSquaredDiff: "
4252 "Incorrect matrix dimensions");
4256 for (
int i = 0; i < A.
nElements(); i++) {
4270 template<
class Treal>
4280 throw Failure(
"Matrix<Treal>::syFrobSquaredDiff: "
4281 "Incorrect matrix dimensions");
4285 for (
int col = 1; col < A.
ncols(); col++)
4286 for (
int row = 0; row < col; row++) {
4287 diff =
A(row, col) -
B(row, col);
4288 sum += 2 * diff * diff;
4290 for (
int rc = 0; rc < A.
ncols(); rc++) {
4291 diff =
A(rc, rc) -
B(rc, rc);
4303 template<
class Treal>
4308 throw Failure(
"Matrix<Treal>::trace(): Matrix must be quadratic");
4313 for (
int rc = 0; rc < this->
ncols(); rc++)
4314 sum += (*
this)(rc,rc);
4319 template<
class Treal>
4326 throw Failure(
"Matrix<Treal>::trace_ab: "
4327 "Wrong matrix dimensions for trace(A * B)");
4330 for (
int rc = 0; rc < A.
nrows(); rc++)
4331 for (
int colA = 0; colA < A.
ncols(); colA++)
4332 tr +=
A(rc,colA) *
B(colA, rc);
4336 template<
class Treal>
4343 throw Failure(
"Matrix<Treal>::trace_aTb: "
4344 "Wrong matrix dimensions for trace(A' * B)");
4347 for (
int rc = 0; rc < A.
ncols(); rc++)
4348 for (
int rowB = 0; rowB < B.
nrows(); rowB++)
4349 tr +=
A(rowB,rc) *
B(rowB, rc);
4353 template<
class Treal>
4361 throw Failure(
"Matrix<Treal>::sy_trace_ab: "
4362 "Wrong matrix dimensions for symmetric trace(A * B)");
4368 for (
int rc = 0; rc < A.
nrows(); rc++)
4369 tr +=
A(rc,rc) *
B(rc, rc);
4371 for (
int colA = 1; colA < A.
ncols(); colA++)
4372 for (
int rowA = 0; rowA < colA; rowA++)
4373 tr += 2 *
A(rowA, colA) *
B(rowA, colA);
4377 template<
class Treal>
4385 throw Failure(
"Matrix<Treal>::add: "
4386 "Wrong matrix dimensions for addition");
4387 if (!A.
is_zero() && alpha != 0) {
4390 for (
int ind = 0; ind < A.
nElements(); ind++)
4395 template<
class Treal>
4410 template<
class Treal>
4417 template<
class Treal>
4421 for (
int ind = 0; ind < this->
nElements(); ind++)
4427 template<
class Treal>
4432 if ((!this->is_zero() && this->frobSquared() <= threshold) ||
4437 else if (!this->is_zero() && this->frobSquared() <= threshold)
4441 template<
class Treal>
4445 assert(!this->is_empty());
4446 bool thisMatIsZero =
true;
4449 bool EMatIsZero =
true;
4450 if (!ErrorMatrix->
is_zero() || !this->is_zero()) {
4453 if (this->is_zero())
4455 for (
int ind = 0; ind < this->nElements(); ind++) {
4456 if ( this->elements[ind] != 0 ) {
4457 assert(ErrorMatrix->
elements[ind] == 0);
4459 if ( this->elements[ind] * this->elements[ind] <= threshold ) {
4461 ErrorMatrix->
elements[ind] = this->elements[ind];
4462 this->elements[ind] = 0;
4466 thisMatIsZero =
false;
4469 if ( ErrorMatrix->
elements[ind] != 0 ) {
4470 assert(this->elements[ind] == 0);
4472 if ( ErrorMatrix->
elements[ind] * ErrorMatrix->
elements[ind] > threshold ) {
4474 this->elements[ind] = ErrorMatrix->
elements[ind];
4476 thisMatIsZero =
false;
4482 if (thisMatIsZero) {
4484 for (
int ind = 0; ind < this->nElements(); ind++)
4485 assert( this->elements[ind] == 0);
4491 for (
int ind = 0; ind < this->nElements(); ind++)
4492 assert( ErrorMatrix->
elements[ind] == 0);
4494 ErrorMatrix->
clear();
4499 if (!this->is_zero()) {
4500 for (
int ind = 0; ind < this->nElements(); ind++) {
4501 if ( this->elements[ind] * this->elements[ind] <= threshold )
4504 this->elements[ind] = 0;
4506 thisMatIsZero =
false;
4519 template<
class Treal>
4537 template<
class Treal>
4550 template<
class Treal>
4553 const bool tA,
const Treal alpha,
4567 template<
class Treal>
4588 template<
class Treal>
4591 assert(!this->is_empty());
4592 if (ErrorMatrix && ErrorMatrix->
is_empty()) {
4596 Treal fs = frobSquared();
4597 if (fs < threshold) {
4607 template<
class Treal>
4611 const Treal threshold,
const side looking,
4615 throw Failure(
"Matrix<Treal>::inch: Matrix must be quadratic!");
4617 throw Failure(
"Matrix<Treal>::inch: Matrix is zero! "
4618 "Not possible to compute inverse cholesky.");
4623 throw Failure(
"Matrix<Treal>::inch: potrf returned info > 0. The matrix is not positive definite.");
4625 throw Failure(
"Matrix<Treal>::inch: potrf returned info < 0");
4629 throw Failure(
"Matrix<Treal>::inch: trtri returned info > 0. The matrix is not positive definite.");
4631 throw Failure(
"Matrix<Treal>::inch: trtri returned info < 0");
4636 template<
class Treal>
4642 Treal* colsums =
new Treal[size];
4643 Treal* diag =
new Treal[size];
4644 for (
int ind = 0; ind < size; ind++)
4650 lmin = diag[0] - tmp1;
4651 lmax = diag[0] + tmp1;
4652 for (
int col = 1; col < size; col++) {
4654 tmp2 = diag[col] - tmp1;
4655 lmin = (tmp2 < lmin ? tmp2 : lmin);
4656 tmp2 = diag[col] + tmp1;
4657 lmax = (tmp2 > lmax ? tmp2 : lmax);
4663 throw Failure(
"Matrix<Treal>::gersgorin(Treal, Treal): Matrix must be quadratic");
4667 template<
class Treal>
4672 for (
int col = 0; col < this->
ncols(); col++)
4673 for (
int row = 0; row < this->
nrows(); row++)
4677 template<
class Treal>
4686 for (
int rc = 0; rc < this->
ncols(); rc++)
4687 diag[rc] = (*
this)(rc,rc);
4714 template<
class Treal>
4716 if (this->
is_zero() && this->cap > 0) {
4722 else if (this->cap > this->nel) {
4723 Treal* tmp =
new Treal[this->nel];
4724 for (
int i = 0; i < this->nel; i++)
4727 this->cap = this->nel;
4730 assert(this->cap == this->nel);
4737 template<
class Treal>
4738 void Matrix<Treal>::
4739 write_to_buffer_count(
int& zb_length,
int& vb_length)
const {
4745 vb_length += this->nel;
4749 template<
class Treal>
4750 void Matrix<Treal>::
4751 write_to_buffer(
int* zerobuf,
const int zb_length,
4752 Treal* valuebuf,
const int vb_length,
4753 int& zb_index,
int& vb_index)
const {
4754 if (zb_index < zb_length) {
4756 zerobuf[zb_index] = 0;
4760 if (vb_index + this->nel < vb_length + 1) {
4761 zerobuf[zb_index] = 1;
4763 for (
int i = 0; i < this->nel; i++)
4764 valuebuf[vb_index + i] = this->
elements[i];
4765 vb_index += this->nel;
4768 throw Failure(
"Matrix<Treal, Telement>::write_to_buffer: "
4769 "Insufficient space in buffers");
4773 throw Failure(
"Matrix<Treal, Telement>::write_to_buffer: "
4774 "Insufficient space in buffers");
4777 template<
class Treal>
4778 void Matrix<Treal>::
4779 read_from_buffer(
int* zerobuf,
const int zb_length,
4780 Treal* valuebuf,
const int vb_length,
4781 int& zb_index,
int& vb_index) {
4782 if (zb_index < zb_length) {
4783 if (zerobuf[zb_index] == 0) {
4788 this->content =
ful;
4790 this->assert_alloc();
4791 if (vb_index + this->nel < vb_length + 1) {
4792 assert(zerobuf[zb_index] == 1);
4794 for (
int i = 0; i < this->nel; i++)
4795 this->
elements[i] = valuebuf[vb_index + i];
4796 vb_index += this->nel;
4799 throw Failure(
"Matrix<Treal, Telement>::read_from_buffer: "
4800 "Mismatch, buffers does not match matrix");
4804 throw Failure(
"Matrix<Treal, Telement>::read_from_buffer: "
4805 "Mismatch, buffers does not match matrix");
void syFullMatrix(std::vector< Treal > &fullMat) const
Definition: Matrix.h:537
Vector< Treal, typename ElementType::VectorType > VectorType
Definition: Matrix.h:116
static Treal sy_trace_ab(const Matrix< Treal, Telement > &A, const Matrix< Treal, Telement > &B)
Definition: Matrix.h:1986
static void trifulltofull(Treal *trifull, const int size)
Definition: gblas.h:1191
void syAssignFrobNormsLowestLevel(Matrix< Treal, Matrix< Treal > > const &A)
Definition: Matrix.h:3132
size_t nnz() const
Returns number of nonzeros in matrix.
Definition: Matrix.h:2858
Definition: Matrix.h:2909
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:71
void assign(Treal const alpha, Matrix< Treal, Telement > const &A)
Definition: Matrix.h:2027
void frobThreshElementLevel(Treal const threshold, Matrix< Treal, Telement > *ErrorMatrix)
Definition: Matrix.h:2099
SizesAndBlocks getSizesAndBlocksForLowerLevel(int const blockNumber) const
Definition: SizesAndBlocks.cc:63
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:2624
void getValues(std::vector< int > const &rowind, std::vector< int > const &colind, std::vector< Treal > &values) const
Definition: Matrix.h:731
size_t nvalues() const
Returns number of stored values in matrix.
Definition: Matrix.h:442
Vector< Treal, Treal > VectorType
Definition: Matrix.h:2912
SizesAndBlocks cols
Definition: MatrixHierarchicBase.h:163
Treal frob() const
Definition: Matrix.h:284
double accumulatedTime
Definition: Matrix.h:82
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:2173
Treal syFrobSquared() const
Definition: Matrix.h:1866
#define MAT_OMP_FINALIZE
Definition: matInclude.h:88
inchversion
Definition: Matrix.h:74
void writeToFile(std::ofstream &file) const
Definition: Matrix.h:843
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:2674
void symToNosym()
Definition: Matrix.h:3877
void tic()
Definition: matInclude.cc:77
void getAllValues(std::vector< int > &rowind, std::vector< int > &colind, std::vector< Treal > &) const
Definition: Matrix.h:806
static void add(const Treal alpha, const Matrix< Treal, Telement > &A, Matrix< Treal, Telement > &B)
Definition: Matrix.h:2009
void getFrobSqLowestLevel(std::vector< Treal > &frobsq) const
Definition: Matrix.h:2045
bool highestLevel() const
Definition: MatrixHierarchicBase.h:131
int getOffset() const
Definition: SizesAndBlocks.h:75
void getFrobSqElementLevel(std::vector< Treal > &frobsq) const
Definition: Matrix.h:2090
int getNTotalScalars() const
Definition: SizesAndBlocks.h:76
Proxy structs used by the matrix API.
Treal syFrobSquared() const
Definition: Matrix.h:4226
void fullMatrix(std::vector< Treal > &fullMat) const
Definition: Matrix.h:517
size_t sy_nvalues() const
Returns number of stored values in matrix.
Definition: Matrix.h:2884
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: gblas.h:340
static void transpose(Matrix< Treal, Telement > const &A, Matrix< Treal, Telement > &AT)
Definition: Matrix.h:977
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: gblas.h:230
Matrix()
Definition: Matrix.h:2915
void setElementsByRule(TRule &rule)
Definition: Matrix.h:938
Vector class.
Definition: Matrix.h:76
MatrixHierarchicBase< Treal, Telement > & operator=(const MatrixHierarchicBase< Treal, Telement > &mat)
Definition: MatrixHierarchicBase.h:200
return elements[row+col *nrows()]
Definition: MatrixHierarchicBase.h:79
void syRandomZeroStructure(Treal probabilityBeingZero)
Definition: Matrix.h:918
size_t sy_nvalues() const
Returns number of stored values in matrix.
Definition: Matrix.h:3313
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:1548
void reset()
Definition: Matrix.h:84
Treal frob_thresh(Treal const threshold, Matrix< Treal > *ErrorMatrix=0)
Definition: Matrix.h:3251
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:2241
Matrix< Treal, Telement > & operator=(const Matrix< Treal, Telement > &mat)
Definition: Matrix.h:196
size_t memory_usage() const
Definition: Matrix.h:3285
void syAssignDiffFrobNormsLowestLevel(Matrix< Treal, Matrix< Treal > > const &A, Matrix< Treal, Matrix< Treal > > const &B)
Definition: Matrix.h:3179
const int & ncols() const
Definition: MatrixHierarchicBase.h:69
const int & nScalarsRows() const
Definition: MatrixHierarchicBase.h:61
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: gblas.h:332
Treal frob() const
Definition: Matrix.h:3059
const int & nrows() const
Definition: MatrixHierarchicBase.h:67
Describes dimensions of matrix and its blocks on all levels.
Definition: SizesAndBlocks.h:37
void get_diagonal(Treal *diag) const
Definition: Matrix.h:2831
void addValues(std::vector< int > const &rowind, std::vector< int > const &colind, std::vector< Treal > const &values)
Definition: Matrix.h:638
static Treal syFrobSquaredDiff(const Matrix< Treal, Telement > &A, const Matrix< Treal, Telement > &B)
Definition: Matrix.h:1908
void syAddValues(std::vector< int > const &rowind, std::vector< int > const &colind, std::vector< Treal > const &values)
Definition: Matrix.h:709
void clear()
Set matrix to zero and delete all arrays.
Definition: Matrix.h:3728
Definition: matInclude.h:158
Treal template_blas_fabs(Treal x)
float toc()
Definition: matInclude.cc:81
Treal syFrob() const
Definition: Matrix.h:289
Definition: matInclude.h:134
void allocate()
Definition: Matrix.h:122
void symToNosym()
Definition: Matrix.h:999
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:2527
Matrix< Treal, Telement > & operator*=(const Treal alpha)
Definition: Matrix.h:1071
static Treal trace_ab(const Matrix< Treal, Telement > &A, const Matrix< Treal, Telement > &B)
Definition: Matrix.h:1952
C++ interface to a subset of BLAS and LAPACK.
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:2261
Base class for Matrix and Matrix specialization.
Definition: MatrixHierarchicBase.h:50
void sy_gersgorin(Treal &lmin, Treal &lmax) const
Definition: Matrix.h:420
static void trtri(const char *uplo, const char *diag, const int *n, T *A, const int *lda, int *info)
Definition: gblas.h:319
void addTime(double timeToAdd)
Definition: Matrix.h:86
static Treal trace_aTb(const Matrix< Treal, Telement > &A, const Matrix< Treal, Telement > &B)
Definition: Matrix.h:1969
void sy_gersgorin(Treal &lmin, Treal &lmax) const
Definition: Matrix.h:3276
void syAssignFrobNormsLowestLevel(Matrix< Treal, Matrix< Treal, Telement > > const &A)
Version of assignFrobNormsLowestLevelToMatrix for symmetric matrices.
Definition: Matrix.h:2151
void random()
Definition: Matrix.h:882
static unsigned int level()
Definition: Matrix.h:3356
Treal trace() const
Definition: Matrix.h:1935
static Treal syFrobDiff(const Matrix< Treal, Telement > &A, const Matrix< Treal, Telement > &B)
Definition: Matrix.h:304
void syRandom()
Definition: Matrix.h:890
void truncateAccordingToSparsityPattern(Matrix< Treal, Matrix< Treal > > &A) const
Definition: Matrix.h:3219
Treal geAccumulateWith(Top &op)
Definition: Matrix.h:3342
static Treal frobDiff(const Matrix< Treal, Telement > &A, const Matrix< Treal, Telement > &B)
Definition: Matrix.h:295
static void potrf(const char *uplo, const int *n, T *A, const int *lda, int *info)
Definition: gblas.h:312
void allocate()
Definition: Matrix.h:2924
static const Treal ONE
Definition: Matrix.h:3426
void freeElements(float *ptr)
Definition: allocate.cc:40
void gersgorin(Treal &lmin, Treal &lmax) const
Definition: Matrix.h:4638
int nElements() const
Definition: MatrixHierarchicBase.h:108
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:1676
size_t sy_nnz() const
Returns number of nonzeros in matrix including lower triangle elements.
Definition: Matrix.h:2867
#define MAT_OMP_START
Definition: matInclude.h:86
static unsigned int level()
Definition: Matrix.h:478
void syGetAllValues(std::vector< int > &rowind, std::vector< int > &colind, std::vector< Treal > &) const
Definition: Matrix.h:818
Telement ElementType
Definition: Matrix.h:115
Matrix class and heart of the matrix library.
Definition: Matrix.h:113
Matrix< Treal > & operator=(const Matrix< Treal > &mat)
Definition: Matrix.h:2983
void randomZeroStructure(Treal probabilityBeingZero)
Get a random zero structure with a specified probability that each submatrix is zero.
Definition: Matrix.h:904
double getAccumulatedTime()
Definition: Matrix.h:85
static Treal frobSquaredDiff(const Matrix< Treal, Telement > &A, const Matrix< Treal, Telement > &B)
Definition: Matrix.h:1886
void frobThreshLowestLevel(Treal const threshold, Matrix< Treal, Telement > *ErrorMatrix)
Definition: Matrix.h:2054
void assignFromSparse(std::vector< int > const &rowind, std::vector< int > const &colind, std::vector< Treal > const &values)
Definition: Matrix.h:586
#define MAT_OMP_END
Definition: matInclude.h:87
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:3267
~Matrix()
Definition: Matrix.h:203
size_t nnz() const
Returns number of nonzeros in matrix.
Definition: Matrix.h:3293
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: gblas.h:279
size_t nvalues() const
Returns number of stored values in matrix.
Definition: Matrix.h:3308
Treal frobSquared() const
Definition: Matrix.h:1852
void assignFromFull(std::vector< Treal > const &fullMat)
Definition: Matrix.h:504
SingletonForTimings()
Definition: Matrix.h:100
size_t sy_nnz() const
Returns number of nonzeros in matrix including lower triangle elements.
Definition: Matrix.h:3299
bool is_empty() const
Check if matrix is empty Empty is different from zero, a zero matrix contains information about block...
Definition: MatrixHierarchicBase.h:141
void syAssignFromSparse(std::vector< int > const &rowind, std::vector< int > const &colind, std::vector< Treal > const &values)
Definition: Matrix.h:687
void syGetValues(std::vector< int > const &rowind, std::vector< int > const &colind, std::vector< Treal > &values) const
Definition: Matrix.h:784
void sySetElementsByRule(TRule &rule)
Definition: Matrix.h:947
static void trsytriplemm(char const side, const Matrix< Treal, Telement > &Z, Matrix< Treal, Telement > &A)
Definition: Matrix.h:2604
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:1362
side
Definition: Matrix.h:73
static SingletonForTimings & instance()
Definition: Matrix.h:94
void syUpTriFullMatrix(std::vector< Treal > &fullMat) const
Definition: Matrix.h:561
void resetCols(SizesAndBlocks const &newCols)
Definition: MatrixHierarchicBase.h:117
int const & getNBlocks() const
Definition: SizesAndBlocks.h:64
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:1772
size_t memory_usage() const
Definition: Matrix.h:2847
void assignDiffFrobNormsLowestLevel(Matrix< Treal, Matrix< Treal > > const &A, Matrix< Treal, Matrix< Treal > > const &B)
Definition: Matrix.h:3147
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:2204
void add_abs_col_sums(Treal *abscolsums) const
Definition: Matrix.h:2816
~Matrix()
Definition: Matrix.h:2989
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:2137
void trSetElementsByRule(TRule &rule)
Definition: Matrix.h:222
Treal geAccumulateWith(Top &op)
Accumulation algorithm for general matrices.
Definition: Matrix.h:468
#define MAT_OMP_INIT
Definition: matInclude.h:85
Treal frobSquared() const
Definition: Matrix.h:4212
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:1247
void readFromFile(std::ifstream &file)
Definition: Matrix.h:859
void nosymToSym()
Definition: Matrix.h:3894
Matrix()
Definition: Matrix.h:119
bool is_zero() const
Definition: MatrixHierarchicBase.h:106
static void sysq(const char uplo, const Treal alpha, const Matrix< Treal, Telement > &A, const Treal beta, Matrix< Treal, Telement > &C)
Definition: Matrix.h:1465
static const Treal ZERO
Definition: Matrix.h:3425
SizesAndBlocks rows
Definition: MatrixHierarchicBase.h:162
Treal syAccumulateWith(Top &op)
Definition: Matrix.h:3323
Xtrans< TX > transpose(TX const &A)
Transposition.
Definition: matrix_proxy.h:129
const int & nScalarsCols() const
Definition: MatrixHierarchicBase.h:63
Treal syFrob() const
Definition: Matrix.h:3061
Treal trace(const XYZ< Treal, MatrixGeneral< Treal, Tmatrix >, MatrixGeneral< Treal, Tmatrix > > &smm)
Definition: MatrixGeneral.h:902
void resetRows(SizesAndBlocks const &newRows)
Definition: MatrixHierarchicBase.h:112
void addIdentity(Treal alpha)
Definition: Matrix.h:962
static void sytr_upper_tr_only(char const side, const Treal alpha, Matrix< Treal, Telement > &A, const Matrix< Treal, Telement > &Z)
Definition: Matrix.h:2444
Treal template_blas_sqrt(Treal x)
Treal maxAbsValue() const
Definition: Matrix.h:480
void nosymToSym()
Definition: Matrix.h:1025
void clear()
Definition: Matrix.h:833
Treal maxAbsValue() const
Definition: Matrix.h:3358
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:1082
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:394
void gersgorin(Treal &lmin, Treal &lmax) const
Definition: Matrix.h:2784
Treal syAccumulateWith(Top &op)
Definition: Matrix.h:455
static void swap(MatrixHierarchicBase< Treal, Telement > &A, MatrixHierarchicBase< Treal, Telement > &B)
Definition: MatrixHierarchicBase.h:221