74 template<
class Treal,
class Telement>
76 template<
typename Tperm>
111 template<
class Treal,
class Telement = Treal>
122 assert(!this->is_empty());
125 #ifdef MAT_USE_ALLOC_TIMER
128 this->
elements =
new Telement[this->nElements()];
129 #ifdef MAT_USE_ALLOC_TIMER
139 (*this)(row,col).resetRows(rowSAB);
140 (*this)(row,col).resetCols(colSAB);
147 void fullMatrix(std::vector<Treal> & fullMat)
const;
153 std::vector<int>
const & colind,
154 std::vector<Treal>
const & values);
156 std::vector<int>
const & colind,
157 std::vector<Treal>
const & values,
158 std::vector<int>
const & indexes);
160 void addValues(std::vector<int>
const & rowind,
161 std::vector<int>
const & colind,
162 std::vector<Treal>
const & values);
163 void addValues(std::vector<int>
const & rowind,
164 std::vector<int>
const & colind,
165 std::vector<Treal>
const & values,
166 std::vector<int>
const & indexes);
169 std::vector<int>
const & colind,
170 std::vector<Treal>
const & values);
173 std::vector<int>
const & colind,
174 std::vector<Treal>
const & values);
176 void getValues(std::vector<int>
const & rowind,
177 std::vector<int>
const & colind,
178 std::vector<Treal> & values)
const;
179 void getValues(std::vector<int>
const & rowind,
180 std::vector<int>
const & colind,
181 std::vector<Treal> &,
182 std::vector<int>
const & indexes)
const;
184 std::vector<int>
const & colind,
185 std::vector<Treal> & values)
const;
187 std::vector<int> & colind,
188 std::vector<Treal> &)
const;
190 std::vector<int> & colind,
191 std::vector<Treal> &)
const;
216 template<
typename TRule>
218 template<
typename TRule>
220 template<
typename TRule>
241 static void gemm(
const bool tA,
const bool tB,
const Treal alpha,
246 static void symm(
const char side,
const char uplo,
const Treal alpha,
251 static void syrk(
const char uplo,
const bool tA,
const Treal alpha,
256 static void sysq(
const char uplo,
const Treal alpha,
261 static void ssmm(
const Treal alpha,
275 static void trmm(
const char side,
const char uplo,
const bool tA,
320 static void add(
const Treal alpha,
323 void assign(Treal
const alpha,
340 (Treal
const threshold,
350 (
Matrix<Treal, Matrix<Treal, Telement> >
const & A );
353 (
Matrix<Treal, Matrix<Treal, Telement> >
const & A );
359 (
Matrix<Treal, Matrix<Treal, Telement> >
const & A,
360 Matrix<Treal, Matrix<Treal, Telement> >
const & B );
365 (
Matrix<Treal, Matrix<Treal, Telement> >
const & A,
366 Matrix<Treal, Matrix<Treal, Telement> >
const & B );
377 const Matrix<Treal, Telement>& A,
378 const Matrix<Treal, Telement>& B,
380 Matrix<Treal, Telement>& C);
382 Matrix<Treal, Telement>& A,
383 const Matrix<Treal, Telement>& Z);
385 const bool tA,
const Treal alpha,
386 const Matrix<Treal, Telement>& A,
387 Matrix<Treal, Telement>& B);
389 const Matrix<Treal, Telement>& Z,
390 Matrix<Treal, Telement>& A);
393 (Treal
const threshold,
404 (Treal
const threshold,
413 const Treal threshold = 0,
417 void gersgorin(Treal& lmin, Treal& lmax)
const;
453 template<
typename Top>
456 if (!this->is_zero()) {
457 for (
int col = 0; col < this->
ncols(); col++) {
458 for (
int row = 0; row < col; row++)
466 template<
typename Top>
469 if (!this->is_zero()) {
470 for (
int col = 0; col < this->
ncols(); col++)
471 for (
int row = 0; row < this->
nrows(); row++)
483 Treal maxAbsGlobal = 0;
484 Treal maxAbsLocal = 0;
485 for (
int ind = 0; ind < this->nElements(); ++ind) {
486 maxAbsLocal = this->
elements[ind].maxAbsValue();
487 maxAbsGlobal = maxAbsGlobal > maxAbsLocal ?
488 maxAbsGlobal : maxAbsLocal;
501 template<
class Treal,
class Telement>
504 int nTotalRows = this->rows.getNTotalScalars();
505 int nTotalCols = this->cols.getNTotalScalars();
506 assert((
int)fullMat.size() == nTotalRows * nTotalCols);
509 for (
int col = 0; col < this->ncols(); col++)
510 for (
int row = 0; row < this->nrows(); row++)
511 (*
this)(row, col).assignFromFull(fullMat);
514 template<
class Treal,
class Telement>
517 int nTotalRows = this->rows.getNTotalScalars();
518 int nTotalCols = this->cols.getNTotalScalars();
519 fullMat.resize(nTotalRows * nTotalCols);
520 if (this->is_zero()) {
521 int rowOffset = this->rows.getOffset();
522 int colOffset = this->cols.getOffset();
523 for (
int col = 0; col < this->nScalarsCols(); col++)
524 for (
int row = 0; row < this->nScalarsRows(); row++)
525 fullMat[(colOffset + col) * nTotalRows + (rowOffset + row)] = 0;
528 for (
int col = 0; col < this->ncols(); col++)
529 for (
int row = 0; row < this->nrows(); row++)
530 (*
this)(row, col).fullMatrix(fullMat);
534 template<
class Treal,
class Telement>
537 int nTotalRows = this->rows.getNTotalScalars();
538 int nTotalCols = this->cols.getNTotalScalars();
539 fullMat.resize(nTotalRows * nTotalCols);
540 if (this->is_zero()) {
541 int rowOffset = this->rows.getOffset();
542 int colOffset = this->cols.getOffset();
543 for (
int col = 0; col < this->nScalarsCols(); col++)
544 for (
int row = 0; row <= col; row++) {
545 fullMat[(colOffset + col) * nTotalRows + (rowOffset + row)] = 0;
546 fullMat[(rowOffset + row) * nTotalRows + (colOffset + col)] = 0;
550 for (
int col = 0; col < this->ncols(); col++) {
551 for (
int row = 0; row < col; row++)
552 (*
this)(row, col).syUpTriFullMatrix(fullMat);
553 (*this)(col, col).syFullMatrix(fullMat);
558 template<
class Treal,
class Telement>
561 int nTotalRows = this->rows.getNTotalScalars();
562 int nTotalCols = this->cols.getNTotalScalars();
563 fullMat.resize(nTotalRows * nTotalCols);
564 if (this->is_zero()) {
565 int rowOffset = this->rows.getOffset();
566 int colOffset = this->cols.getOffset();
567 for (
int col = 0; col < this->nScalarsCols(); col++)
568 for (
int row = 0; row <= this->nScalarsRows(); row++) {
569 fullMat[(colOffset + col) * nTotalRows + (rowOffset + row)] = 0;
570 fullMat[(rowOffset + row) * nTotalRows + (colOffset + col)] = 0;
574 for (
int col = 0; col < this->ncols(); col++)
575 for (
int row = 0; row < this->nrows(); row++)
576 (*
this)(row, col).syUpTriFullMatrix(fullMat);
583 template<
class Treal,
class Telement>
586 std::vector<int>
const & colind,
587 std::vector<Treal>
const & values) {
588 assert(rowind.size() == colind.size() &&
589 rowind.size() == values.size());
590 std::vector<int> indexes(values.size());
591 for (
unsigned int ind = 0; ind < values.size(); ++ind)
593 assignFromSparse(rowind, colind, values, indexes);
596 template<
class Treal,
class Telement>
599 std::vector<int>
const & colind,
600 std::vector<Treal>
const & values,
601 std::vector<int>
const & indexes) {
602 if (indexes.empty()) {
609 std::vector<std::vector<int> > columnBuckets (this->cols.getNBlocks());
610 std::vector<std::vector<int> > rowBuckets (this->rows.getNBlocks());
614 std::vector<int>::const_iterator it;
615 for ( it = indexes.begin(); it < indexes.end(); it++ )
616 columnBuckets[ this->cols.whichBlock( colind[*it] ) ].push_back (*it);
619 for (
int col = 0; col < this->ncols(); col++) {
621 while (!columnBuckets[col].empty()) {
622 currentInd = columnBuckets[col].back();
623 columnBuckets[col].pop_back();
624 rowBuckets[ this->rows.whichBlock
625 ( rowind[currentInd] ) ].push_back (currentInd);
628 for (
int row = 0; row < this->nrows(); row++) {
629 (*this)(row,col).assignFromSparse(rowind, colind, values, rowBuckets[row]);
630 rowBuckets[row].clear();
635 template<
class Treal,
class Telement>
638 std::vector<int>
const & colind,
639 std::vector<Treal>
const & values) {
640 assert(rowind.size() == colind.size() &&
641 rowind.size() == values.size());
642 std::vector<int> indexes(values.size());
643 for (
unsigned int ind = 0; ind < values.size(); ++ind)
645 addValues(rowind, colind, values, indexes);
648 template<
class Treal,
class Telement>
651 std::vector<int>
const & colind,
652 std::vector<Treal>
const & values,
653 std::vector<int>
const & indexes) {
659 std::vector<std::vector<int> > columnBuckets (this->cols.getNBlocks());
660 std::vector<std::vector<int> > rowBuckets (this->rows.getNBlocks());
663 std::vector<int>::const_iterator it;
664 for ( it = indexes.begin(); it < indexes.end(); it++ )
665 columnBuckets[ this->cols.whichBlock( colind[*it] ) ].push_back (*it);
668 for (
int col = 0; col < this->ncols(); col++) {
670 while (!columnBuckets[col].empty()) {
671 currentInd = columnBuckets[col].back();
672 columnBuckets[col].pop_back();
673 rowBuckets[ this->rows.whichBlock
674 ( rowind[currentInd] ) ].push_back (currentInd);
677 for (
int row = 0; row < this->nrows(); row++) {
678 (*this)(row,col).addValues(rowind, colind, values, rowBuckets[row]);
679 rowBuckets[row].clear();
684 template<
class Treal,
class Telement>
687 std::vector<int>
const & colind,
688 std::vector<Treal>
const & values) {
689 assert(rowind.size() == colind.size() &&
690 rowind.size() == values.size());
691 bool upperTriangleOnly =
true;
692 for (
unsigned int ind = 0; ind < values.size(); ++ind) {
694 upperTriangleOnly && colind[ind] >= rowind[ind];
696 if (!upperTriangleOnly)
697 throw Failure(
"Matrix<Treal, Telement>::"
698 "syAddValues(std::vector<int> const &, "
699 "std::vector<int> const &, "
700 "std::vector<Treal> const &, int const): "
701 "Only upper triangle can contain elements when assigning "
702 "symmetric or triangular matrix ");
703 assignFromSparse(rowind, colind, values);
706 template<
class Treal,
class Telement>
709 std::vector<int>
const & colind,
710 std::vector<Treal>
const & values) {
711 assert(rowind.size() == colind.size() &&
712 rowind.size() == values.size());
713 bool upperTriangleOnly =
true;
714 for (
unsigned int ind = 0; ind < values.size(); ++ind) {
716 upperTriangleOnly && colind[ind] >= rowind[ind];
718 if (!upperTriangleOnly)
719 throw Failure(
"Matrix<Treal, Telement>::"
720 "syAddValues(std::vector<int> const &, "
721 "std::vector<int> const &, "
722 "std::vector<Treal> const &, int const): "
723 "Only upper triangle can contain elements when assigning "
724 "symmetric or triangular matrix ");
725 addValues(rowind, colind, values);
728 template<
class Treal,
class Telement>
731 std::vector<int>
const & colind,
732 std::vector<Treal> & values)
const {
733 assert(rowind.size() == colind.size());
734 values.resize(rowind.size(),0);
735 std::vector<int> indexes(rowind.size());
736 for (
unsigned int ind = 0; ind < rowind.size(); ++ind)
738 getValues(rowind, colind, values, indexes);
741 template<
class Treal,
class Telement>
744 std::vector<int>
const & colind,
745 std::vector<Treal> & values,
746 std::vector<int>
const & indexes)
const {
747 assert(!this->is_empty());
750 std::vector<int>::const_iterator it;
751 if (this->is_zero()) {
752 for ( it = indexes.begin(); it < indexes.end(); it++ )
757 std::vector<std::vector<int> > columnBuckets (this->cols.getNBlocks());
758 std::vector<std::vector<int> > rowBuckets (this->rows.getNBlocks());
761 for ( it = indexes.begin(); it < indexes.end(); it++ )
762 columnBuckets[ this->cols.whichBlock( colind[*it] ) ].push_back (*it);
765 for (
int col = 0; col < this->ncols(); col++) {
767 while (!columnBuckets[col].empty()) {
768 currentInd = columnBuckets[col].back();
769 columnBuckets[col].pop_back();
770 rowBuckets[ this->rows.whichBlock
771 ( rowind[currentInd] ) ].push_back (currentInd);
774 for (
int row = 0; row < this->nrows(); row++) {
775 (*this)(row,col).getValues(rowind, colind, values, rowBuckets[row]);
776 rowBuckets[row].clear();
781 template<
class Treal,
class Telement>
784 std::vector<int>
const & colind,
785 std::vector<Treal> & values)
const {
786 assert(rowind.size() == colind.size());
787 bool upperTriangleOnly =
true;
788 for (
int unsigned ind = 0; ind < rowind.size(); ++ind) {
790 upperTriangleOnly && colind[ind] >= rowind[ind];
792 if (!upperTriangleOnly)
793 throw Failure(
"Matrix<Treal, Telement>::"
794 "syGetValues(std::vector<int> const &, "
795 "std::vector<int> const &, "
796 "std::vector<Treal> const &, int const): "
797 "Only upper triangle when retrieving elements from "
798 "symmetric or triangular matrix ");
799 getValues(rowind, colind, values);
803 template<
class Treal,
class Telement>
806 std::vector<int> & colind,
807 std::vector<Treal> & values)
const {
808 assert(rowind.size() == colind.size() &&
809 rowind.size() == values.size());
810 if (!this->is_zero())
811 for (
int ind = 0; ind < this->nElements(); ++ind)
812 this->elements[ind].getAllValues(rowind, colind, values);
815 template<
class Treal,
class Telement>
818 std::vector<int> & colind,
819 std::vector<Treal> & values)
const {
820 assert(rowind.size() == colind.size() &&
821 rowind.size() == values.size());
822 if (!this->is_zero())
823 for (
int col = 0; col < this->ncols(); ++col) {
824 for (
int row = 0; row < col; ++row)
825 (*
this)(row, col).getAllValues(rowind, colind, values);
826 (*this)(col, col).syGetAllValues(rowind, colind, values);
831 template<
class Treal,
class Telement>
833 if (!this->is_zero())
834 for (
int i = 0; i < this->nElements(); i++)
835 this->elements[i].clear();
836 delete[] this->elements;
840 template<
class Treal,
class Telement>
845 if (this->is_zero()) {
846 char * tmp = (
char*)&ZERO;
847 file.write(tmp,
sizeof(
int));
850 char * tmp = (
char*)&ONE;
851 file.write(tmp,
sizeof(
int));
852 for (
int i = 0; i < this->nElements(); i++)
853 this->elements[i].writeToFile(file);
856 template<
class Treal,
class Telement>
861 char tmp[
sizeof(int)];
862 file.read(tmp, (std::ifstream::pos_type)
sizeof(
int));
870 for (
int i = 0; i < this->nElements(); i++)
871 this->elements[i].readFromFile(file);
874 throw Failure(
"Matrix<Treal, Telement>::"
875 "readFromFile(std::ifstream & file):"
876 "File corruption int value not 0 or 1");
880 template<
class Treal,
class Telement>
884 for (
int ind = 0; ind < this->nElements(); ind++)
885 this->elements[ind].random();
888 template<
class Treal,
class Telement>
893 for (
int col = 1; col < this->ncols(); col++)
894 for (
int row = 0; row < col; row++)
895 (*
this)(row, col).random();
897 for (
int rc = 0; rc < this->nrows(); rc++)
898 (*
this)(rc,rc).syRandom();
901 template<
class Treal,
class Telement>
904 if (!this->highestLevel() &&
905 probabilityBeingZero > rand() / (Treal)RAND_MAX)
910 for (
int ind = 0; ind < this->nElements(); ind++)
911 this->elements[ind].randomZeroStructure(probabilityBeingZero);
915 template<
class Treal,
class Telement>
918 if (!this->highestLevel() &&
919 probabilityBeingZero > rand() / (Treal)RAND_MAX)
925 for (
int col = 1; col < this->ncols(); col++)
926 for (
int row = 0; row < col; row++)
927 (*
this)(row, col).randomZeroStructure(probabilityBeingZero);
929 for (
int rc = 0; rc < this->nrows(); rc++)
930 (*
this)(rc,rc).syRandomZeroStructure(probabilityBeingZero);
934 template<
class Treal,
class Telement>
935 template<
typename TRule>
940 for (
int ind = 0; ind < this->nElements(); ind++)
941 this->elements[ind].setElementsByRule(rule);
943 template<
class Treal,
class Telement>
944 template<
typename TRule>
950 for (
int col = 1; col < this->ncols(); col++)
951 for (
int row = 0; row < col; row++)
952 (*
this)(row, col).setElementsByRule(rule);
954 for (
int rc = 0; rc < this->nrows(); rc++)
955 (*
this)(rc,rc).sySetElementsByRule(rule);
959 template<
class Treal,
class Telement>
962 if (this->is_empty())
963 throw Failure(
"Matrix<Treal, Telement>::addIdentity(Treal): "
964 "Cannot add identity to empty matrix.");
965 if (this->ncols() != this->nrows())
966 throw Failure(
"Matrix<Treal, Telement>::addIdentity(Treal): "
967 "Matrix must be square to add identity");
970 for (
int ind = 0; ind < this->ncols(); ind++)
971 (*
this)(ind,ind).addIdentity(alpha);
974 template<
class Treal,
class Telement>
985 if (AT.is_zero() || (AT.nElements() != A.nElements())) {
987 AT.
elements =
new Telement[A.nElements()];
991 for (
int row = 0; row < AT.
nrows(); row++)
992 for (
int col = 0; col < AT.
ncols(); col++)
996 template<
class Treal,
class Telement>
1000 if (this->nrows() == this->ncols()) {
1001 if (!this->is_zero()) {
1003 for (
int rc = 0; rc < this->ncols(); rc++)
1004 (*
this)(rc, rc).symToNosym();
1006 for (
int row = 1; row < this->nrows(); row++)
1007 for (
int col = 0; col < row; col++)
1012 throw Failure(
"Matrix<Treal, Telement>::symToNosym(): "
1013 "Only quadratic matrices can be symmetric");
1016 std::cout<<
"Failure caught:Matrix<Treal, Telement>::symToNosym()"
1022 template<
class Treal,
class Telement>
1025 if (this->nrows() == this->ncols()) {
1026 if (!this->is_zero()) {
1028 for (
int rc = 0; rc < this->ncols(); rc++)
1029 (*
this)(rc, rc).nosymToSym();
1031 for (
int col = 0; col < this->ncols() - 1; col++)
1032 for (
int row = col + 1; row < this->nrows(); row++)
1033 (*
this)(row, col).clear();
1037 throw Failure(
"Matrix<Treal, Telement>::nosymToSym(): "
1038 "Only quadratic matrices can be symmetric");
1043 template<
class Treal,
class Telement>
1051 if (this->ncols() != this->nrows())
1053 (
"Matrix::operator=(int k = 1): "
1054 "Matrix must be quadratic to become identity matrix.");
1057 for (
int rc = 0; rc < this->ncols(); rc++)
1061 throw Failure(
"Matrix::operator=(int k) only "
1062 "implemented for k = 0, k = 1");
1068 template<
class Treal,
class Telement>
1071 if (!this->is_zero() && alpha != 1) {
1072 for (
int ind = 0; ind < this->nElements(); ind++)
1073 this->elements[ind] *= alpha;
1079 template<
class Treal,
class Telement>
1081 gemm(
const bool tA,
const bool tB,
const Treal alpha,
1086 assert(!A.is_empty());
1087 assert(!B.is_empty());
1091 C.resetRows(A.
cols);
1093 C.resetRows(A.
rows);
1095 C.resetCols(B.
rows);
1097 C.resetCols(B.
cols);
1100 if ((A.is_zero() || B.is_zero() || alpha == 0) &&
1101 (C.is_zero() || beta == 0))
1104 Treal beta_tmp = beta;
1109 if (!A.is_zero() && !B.is_zero() && alpha != 0) {
1116 #pragma omp parallel for if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) schedule(dynamic)
1118 for (
int col = 0; col < C.
ncols(); col++) {
1120 for (
int row = 0; row < C.
nrows(); row++) {
1122 A(row, 0), B(0, col),
1125 for(
int cola = 1; cola < A.
ncols(); cola++)
1127 A(row, cola), B(cola, col),
1135 throw Failure(
"Matrix<class Treal, class Telement>::"
1136 "gemm(bool, bool, Treal, "
1137 "Matrix<Treal, Telement>, "
1138 "Matrix<Treal, Telement>, Treal, "
1139 "Matrix<Treal, Telement>): "
1140 "Incorrect matrixdimensions for multiplication");
1146 #pragma omp parallel for if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) schedule(dynamic)
1148 for (
int col = 0; col < C.
ncols(); col++) {
1150 for (
int row = 0; row < C.
nrows(); row++) {
1155 for(
int cola = 1; cola < A.
nrows(); cola++)
1157 A(cola, row), B(cola, col),
1165 throw Failure(
"Matrix<class Treal, class Telement>::"
1166 "gemm(bool, bool, Treal, "
1167 "Matrix<Treal, Telement>, "
1168 "Matrix<Treal, Telement>, Treal, "
1169 "Matrix<Treal, Telement>): "
1170 "Incorrect matrixdimensions for multiplication");
1176 #pragma omp parallel for if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) schedule(dynamic)
1178 for (
int col = 0; col < C.
ncols(); col++) {
1180 for (
int row = 0; row < C.
nrows(); row++) {
1182 A(row, 0), B(col, 0),
1185 for(
int cola = 1; cola < A.
ncols(); cola++)
1187 A(row, cola), B(col, cola),
1195 throw Failure(
"Matrix<class Treal, class Telement>::"
1196 "gemm(bool, bool, Treal, "
1197 "Matrix<Treal, Telement>, "
1198 "Matrix<Treal, Telement>, Treal, "
1199 "Matrix<Treal, Telement>): "
1200 "Incorrect matrixdimensions for multiplication");
1206 #pragma omp parallel for if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) schedule(dynamic)
1208 for (
int col = 0; col < C.
ncols(); col++) {
1210 for (
int row = 0; row < C.
nrows(); row++) {
1212 A(0, row), B(col, 0),
1215 for(
int cola = 1; cola < A.
nrows(); cola++)
1217 A(cola, row), B(col, cola),
1225 throw Failure(
"Matrix<class Treal, class Telement>::"
1226 "gemm(bool, bool, Treal, "
1227 "Matrix<Treal, Telement>, "
1228 "Matrix<Treal, Telement>, "
1229 "Treal, Matrix<Treal, Telement>): "
1230 "Incorrect matrixdimensions for multiplication");
1231 else throw Failure(
"Matrix<class Treal, class Telement>::"
1232 "gemm(bool, bool, Treal, "
1233 "Matrix<Treal, Telement>, "
1234 "Matrix<Treal, Telement>, Treal, "
1235 "Matrix<Treal, Telement>):"
1236 "Very strange error!!");
1244 template<
class Treal,
class Telement>
1246 symm(
const char side,
const char uplo,
const Treal alpha,
1251 assert(!A.is_empty());
1252 assert(!B.is_empty());
1254 int dimA = A.
nrows();
1258 C.resetRows(A.
rows);
1259 C.resetCols(B.
cols);
1262 assert(side ==
'R');
1263 C.resetRows(B.
rows);
1264 C.resetCols(A.
cols);
1267 if ((A.is_zero() || B.is_zero() || alpha == 0) &&
1268 (C.is_zero() || beta == 0))
1272 Treal beta_tmp = beta;
1277 if (!A.is_zero() && !B.is_zero() && alpha != 0) {
1280 if (dimA == B.
nrows() &&
1281 dimA == C.
nrows() &&
1284 #pragma omp parallel for if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) schedule(dynamic)
1286 for (
int col = 0; col < C.
ncols(); col++) {
1288 for (
int row = 0; row < C.
nrows(); row++) {
1291 B(row, col), beta_tmp, C(row, col));
1293 for(
int ind = 0; ind < row; ind++)
1295 A(ind, row), B(ind, col),
1298 for(
int ind = row + 1; ind < dimA; ind++)
1300 A(row, ind), B(ind, col),
1307 throw Failure(
"Matrix<class Treal, class Telement>"
1308 "::symm(bool, bool, Treal, Matrix<Treal, "
1309 "Telement>, Matrix<Treal, Telement>, "
1310 "Treal, Matrix<Treal, Telement>): "
1311 "Incorrect matrixdimensions for multiplication");
1313 assert(side ==
'R');
1314 if (B.
ncols() == dimA &&
1316 dimA == C.
ncols()) {
1318 #pragma omp parallel for if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) schedule(dynamic)
1320 for (
int col = 0; col < C.
ncols(); col++) {
1322 for (
int row = 0; row < C.
nrows(); row++) {
1325 B(row, col), beta_tmp, C(row, col));
1327 for(
int ind = col + 1; ind < dimA; ind++)
1329 B(row, ind), A(col, ind),
1332 for(
int ind = 0; ind < col; ind++)
1334 B(row, ind), A(ind, col),
1341 throw Failure(
"Matrix<class Treal, class Telement>"
1342 "::symm(bool, bool, Treal, Matrix<Treal, "
1343 "Telement>, Matrix<Treal, Telement>, "
1344 "Treal, Matrix<Treal, Telement>): "
1345 "Incorrect matrixdimensions for multiplication");
1353 throw Failure(
"Matrix<class Treal, class Telement>::"
1354 "symm only implemented for symmetric matrices in "
1355 "upper triangular storage");
1359 template<
class Treal,
class Telement>
1361 syrk(
const char uplo,
const bool tA,
const Treal alpha,
1365 assert(!A.is_empty());
1369 C.resetRows(A.
cols);
1370 C.resetCols(A.
cols);
1373 C.resetRows(A.
rows);
1374 C.resetCols(A.
rows);
1380 if (alpha != 0 && !A.is_zero()) {
1381 Treal beta_tmp = beta;
1387 if (!tA && uplo ==
'U') {
1389 #pragma omp parallel if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared)
1393 #pragma omp for schedule(dynamic) nowait
1395 for (
int rc = 0; rc < C.
ncols(); rc++) {
1398 for (
int cola = 1; cola < A.
ncols(); cola++)
1403 #pragma omp for schedule(dynamic) nowait
1405 for (
int row = 0; row < C.
nrows() - 1; row++) {
1407 for (
int col = row + 1; col < C.
ncols(); col++) {
1409 A(row, 0), A(col,0), beta_tmp, C(row,col));
1410 for (
int ind = 1; ind < A.
ncols(); ind++)
1412 A(row, ind), A(col,ind), 1.0, C(row,col));
1418 else if (tA && uplo ==
'U') {
1420 #pragma omp parallel if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared)
1424 #pragma omp for schedule(dynamic) nowait
1426 for (
int rc = 0; rc < C.
ncols(); rc++) {
1429 for (
int rowa = 1; rowa < A.
nrows(); rowa++)
1434 #pragma omp for schedule(dynamic) nowait
1436 for (
int row = 0; row < C.
nrows() - 1; row++) {
1438 for (
int col = row + 1; col < C.
ncols(); col++) {
1440 A(0, row), A(0, col), beta_tmp, C(row,col));
1441 for (
int ind = 1; ind < A.
nrows(); ind++)
1443 A(ind, row), A(ind, col), 1.0, C(row,col));
1450 throw Failure(
"Matrix<class Treal, class Telement>::"
1451 "syrk not implemented for lower triangular storage");
1458 throw Failure(
"Matrix<class Treal, class Telement>::syrk: "
1459 "Incorrect matrix dimensions for symmetric rank-k update");
1462 template<
class Treal,
class Telement>
1464 sysq(
const char uplo,
const Treal alpha,
1468 assert(!A.is_empty());
1471 C.resetRows(A.
rows);
1472 C.resetCols(A.
cols);
1476 if (alpha != 0 && !A.is_zero()) {
1478 Treal beta_tmp = beta;
1485 #pragma omp parallel if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared)
1489 #pragma omp for schedule(dynamic) nowait
1491 for (
int rc = 0; rc < C.
ncols(); rc++) {
1493 Telement::sysq(uplo, alpha, A(rc, rc), beta_tmp, C(rc, rc));
1494 for (
int cola = 0; cola < rc; cola++)
1496 for (
int cola = rc + 1; cola < A.
ncols(); cola++)
1502 #pragma omp for schedule(dynamic) nowait
1504 for (
int row = 0; row < C.
nrows() - 1; row++) {
1506 for (
int col = row + 1; col < C.
ncols(); col++) {
1509 beta_tmp, C(row, col));
1513 for (
int ind = 0; ind < row; ind++)
1515 A(ind, row), A(ind, col), 1.0, C(row, col));
1517 for (
int ind = row + 1; ind < col; ind++)
1519 A(row, ind), A(ind, col), 1.0, C(row, col));
1521 for (
int ind = col + 1; ind < A.
ncols(); ind++)
1523 A(row, ind), A(col, ind), 1.0, C(row, col));
1531 throw Failure(
"Matrix<class Treal, class Telement>::"
1532 "sysq only implemented for symmetric matrices in "
1533 "upper triangular storage");;
1539 throw Failure(
"Matrix<class Treal, class Telement>::sysq: "
1540 "Incorrect matrix dimensions for symmetric square "
1545 template<
class Treal,
class Telement>
1552 assert(!A.is_empty());
1553 assert(!B.is_empty());
1556 C.resetRows(A.
rows);
1557 C.resetCols(B.
cols);
1564 throw Failure(
"Matrix<class Treal, class Telement>::"
1566 "Matrix<Treal, Telement>, "
1567 "Matrix<Treal, Telement>, Treal, "
1568 "Matrix<Treal, Telement>): "
1569 "Incorrect matrixdimensions for ssmm multiplication");
1571 if ((A.is_zero() || B.is_zero() || alpha == 0) &&
1572 (C.is_zero() || beta == 0)) {
1576 Treal beta_tmp = beta;
1581 if (A.is_zero() || B.is_zero() || alpha == 0) {
1587 #pragma omp parallel for if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) schedule(dynamic)
1589 for (
int col = 0; col < C.
ncols(); col++) {
1593 Telement::ssmm(alpha, A(col,col), B(col, col),
1594 beta_tmp, C(col,col));
1595 for (
int ind = 0; ind < col; ind++)
1598 alpha, A(ind,col), B(ind, col),
1600 for (
int ind = col + 1; ind < A.
ncols(); ind++)
1603 alpha, A(col, ind), B(col, ind),
1606 for (
int row = 0; row < col; row++) {
1611 alpha, A(row, row), B(row, col),
1612 beta_tmp, C(row, col));
1615 alpha, B(col, col), A(row, col),
1617 for (
int ind = 0; ind < row; ind++)
1620 alpha, A(ind, row), B(ind, col),
1622 for (
int ind = row + 1; ind < col; ind++)
1625 alpha, A(row, ind), B(ind, col),
1627 for (
int ind = col + 1; ind < A.
ncols(); ind++)
1630 alpha, A(row, ind), B(col, ind),
1635 for (
int row = col + 1; row < C.
nrows(); row++) {
1641 alpha, B(col, col), A(col, row),
1642 beta_tmp, tmpSubMat);
1645 alpha, A(row, row), B(col, row),
1647 for (
int ind = 0; ind < col; ind++)
1650 alpha, B(ind, col), A(ind, row),
1652 for (
int ind = col + 1; ind < row; ind++)
1655 alpha, B(col, ind), A(ind, row),
1657 for (
int ind = row + 1; ind < B.
nrows(); ind++)
1660 alpha, B(col, ind), A(row, ind),
1673 template<
class Treal,
class Telement>
1680 assert(!A.is_empty());
1681 assert(!B.is_empty());
1684 C.resetRows(A.
rows);
1685 C.resetCols(B.
cols);
1692 throw Failure(
"Matrix<class Treal, class Telement>::"
1693 "ssmm_upper_tr_only(Treal, "
1694 "Matrix<Treal, Telement>, "
1695 "Matrix<Treal, Telement>, Treal, "
1696 "Matrix<Treal, Telement>): "
1697 "Incorrect matrixdimensions for ssmm_upper_tr_only "
1700 if ((A.is_zero() || B.is_zero() || alpha == 0) &&
1701 (C.is_zero() || beta == 0)) {
1705 Treal beta_tmp = beta;
1710 if (A.is_zero() || B.is_zero() || alpha == 0) {
1716 #pragma omp parallel for if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) schedule(dynamic)
1718 for (
int col = 0; col < C.
ncols(); col++) {
1722 Telement::ssmm_upper_tr_only(alpha, A(col,col), B(col, col),
1723 beta_tmp, C(col,col));
1724 for (
int ind = 0; ind < col; ind++)
1726 Telement::gemm_upper_tr_only(
true,
false,
1727 alpha, A(ind,col), B(ind, col),
1729 for (
int ind = col + 1; ind < A.
ncols(); ind++)
1731 Telement::gemm_upper_tr_only(
false,
true,
1732 alpha, A(col, ind), B(col, ind),
1735 for (
int row = 0; row < col; row++) {
1740 alpha, A(row, row), B(row, col),
1741 beta_tmp, C(row, col));
1744 alpha, B(col, col), A(row, col),
1746 for (
int ind = 0; ind < row; ind++)
1749 alpha, A(ind, row), B(ind, col),
1751 for (
int ind = row + 1; ind < col; ind++)
1754 alpha, A(row, ind), B(ind, col),
1756 for (
int ind = col + 1; ind < A.
ncols(); ind++)
1759 alpha, A(row, ind), B(col, ind),
1769 template<
class Treal,
class Telement>
1771 trmm(
const char side,
const char uplo,
const bool tA,
const Treal alpha,
1774 assert(!A.is_empty());
1775 assert(!B.is_empty());
1776 if (alpha != 0 && !A.is_zero() && !B.is_zero())
1777 if (((side ==
'R' && B.
ncols() == A.
nrows()) ||
1784 for (
int col = B.
ncols() - 1; col >= 0; col--)
1785 for (
int row = 0; row < B.
nrows(); row++) {
1789 A(col, col), B(row,col));
1791 for (
int ind = 0; ind < col; ind++)
1793 B(row,ind), A(ind, col),
1798 assert(side ==
'L');
1800 for (
int row = 0; row < B.
nrows(); row++ )
1801 for (
int col = 0; col < B.
ncols(); col++) {
1803 A(row,row), B(row,col));
1804 for (
int ind = row + 1 ; ind < B.
nrows(); ind++)
1806 A(row,ind), B(ind,col),
1815 for (
int col = 0; col < B.
ncols(); col++)
1816 for (
int row = 0; row < B.
nrows(); row++) {
1818 A(col,col), B(row,col));
1819 for (
int ind = col + 1; ind < A.
ncols(); ind++)
1821 B(row,ind), A(col,ind),
1826 assert(side ==
'L');
1828 for (
int row = B.
nrows() - 1; row >= 0; row--)
1829 for (
int col = 0; col < B.
ncols(); col++) {
1831 A(row,row), B(row,col));
1832 for (
int ind = 0; ind < row; ind++)
1834 A(ind,row), B(ind,col),
1840 throw Failure(
"Matrix<class Treal, class Telement>::"
1841 "trmm not implemented for lower triangular matrices");
1843 throw Failure(
"Matrix<class Treal, class Telement>::trmm"
1844 ": Incorrect matrix dimensions for multiplication");
1850 template<
class Treal,
class Telement>
1852 assert(!this->is_empty());
1853 if (this->is_zero())
1857 for (
int i = 0; i < this->nElements(); i++)
1858 sum += this->elements[i].frobSquared();
1863 template<
class Treal,
class Telement>
1866 assert(!this->is_empty());
1867 if (this->nrows() != this->ncols())
1868 throw Failure(
"Matrix<Treal, Telement>::syFrobSquared(): "
1869 "Matrix must be have equal number of rows "
1870 "and cols to be symmetric");
1872 if (!this->is_zero()) {
1873 for (
int col = 1; col < this->ncols(); col++)
1874 for (
int row = 0; row < col; row++)
1875 sum += 2 * (*
this)(row, col).frobSquared();
1876 for (
int rc = 0; rc < this->ncols(); rc++)
1877 sum += (*
this)(rc, rc).syFrobSquared();
1882 template<
class Treal,
class Telement>
1887 assert(!A.is_empty());
1888 assert(!B.is_empty());
1890 throw Failure(
"Matrix<Treal, Telement>::"
1891 "frobSquaredDiff: Incorrect matrix dimensions");
1893 if (!A.is_zero() && !B.is_zero())
1894 for (
int i = 0; i < A.nElements(); i++)
1896 else if (A.is_zero() && !B.is_zero())
1898 else if (!A.is_zero() && B.is_zero())
1904 template<
class Treal,
class Telement>
1909 assert(!A.is_empty());
1910 assert(!B.is_empty());
1914 throw Failure(
"Matrix<Treal, Telement>::syFrobSquaredDiff: "
1915 "Incorrect matrix dimensions");
1917 if (!A.is_zero() && !B.is_zero()) {
1918 for (
int col = 1; col < A.
ncols(); col++)
1919 for (
int row = 0; row < col; row++)
1920 sum += 2 * Telement::frobSquaredDiff(A(row, col), B(row, col));
1921 for (
int rc = 0; rc < A.
ncols(); rc++)
1922 sum += Telement::syFrobSquaredDiff(A(rc, rc),B(rc, rc));
1924 else if (A.is_zero() && !B.is_zero())
1926 else if (!A.is_zero() && B.is_zero())
1932 template<
class Treal,
class Telement>
1935 assert(!this->is_empty());
1936 if (this->nrows() != this->ncols())
1937 throw Failure(
"Matrix<Treal, Telement>::trace(): "
1938 "Matrix must be quadratic");
1939 if (this->is_zero())
1943 for (
int rc = 0; rc < this->ncols(); rc++)
1944 sum += (*
this)(rc,rc).
trace();
1949 template<
class Treal,
class Telement>
1953 assert(!A.is_empty());
1954 assert(!B.is_empty());
1956 throw Failure(
"Matrix<Treal, Telement>::trace_ab: "
1957 "Wrong matrix dimensions for trace(A * B)");
1959 if (!A.is_zero() && !B.is_zero())
1960 for (
int rc = 0; rc < A.
nrows(); rc++)
1961 for (
int colA = 0; colA < A.
ncols(); colA++)
1962 tr += Telement::trace_ab(A(rc,colA), B(colA, rc));
1966 template<
class Treal,
class Telement>
1970 assert(!A.is_empty());
1971 assert(!B.is_empty());
1973 throw Failure(
"Matrix<Treal, Telement>::trace_aTb: "
1974 "Wrong matrix dimensions for trace(A' * B)");
1976 if (!A.is_zero() && !B.is_zero())
1977 for (
int rc = 0; rc < A.
ncols(); rc++)
1978 for (
int rowB = 0; rowB < B.
nrows(); rowB++)
1979 tr += Telement::trace_aTb(A(rowB,rc), B(rowB, rc));
1983 template<
class Treal,
class Telement>
1987 assert(!A.is_empty());
1988 assert(!B.is_empty());
1991 throw Failure(
"Matrix<Treal, Telement>::sy_trace_ab: "
1992 "Wrong matrix dimensions for symmetric trace(A * B)");
1994 if (!A.is_zero() && !B.is_zero()) {
1996 for (
int rc = 0; rc < A.
nrows(); rc++)
1997 tr += Telement::sy_trace_ab(A(rc,rc), B(rc, rc));
1999 for (
int colA = 1; colA < A.
ncols(); colA++)
2000 for (
int rowA = 0; rowA < colA; rowA++)
2001 tr += 2 * Telement::trace_aTb(A(rowA, colA), B(rowA, colA));
2006 template<
class Treal,
class Telement>
2011 assert(!A.is_empty());
2012 assert(!B.is_empty());
2014 throw Failure(
"Matrix<Treal, Telement>::add: "
2015 "Wrong matrix dimensions for addition");
2016 if (!A.is_zero() && alpha != 0) {
2019 for (
int ind = 0; ind < A.nElements(); ind++)
2024 template<
class Treal,
class Telement>
2028 assert(!A.is_empty());
2029 if (this->is_empty()) {
2030 this->resetRows(A.
rows);
2031 this->resetCols(A.
cols);
2035 add(alpha, A, *
this);
2042 template<
class Treal,
class Telement>
2045 if (!this->is_zero())
2046 for (
int ind = 0; ind < this->nElements(); ind++)
2047 this->elements[ind].getFrobSqLowestLevel(frobsq);
2050 template<
class Treal,
class Telement>
2054 assert(!this->is_empty());
2055 bool thisMatIsZero =
true;
2057 assert(!ErrorMatrix->is_empty());
2058 bool EMatIsZero =
true;
2059 if (!ErrorMatrix->is_zero() || !this->is_zero()) {
2060 if (ErrorMatrix->is_zero())
2062 if (this->is_zero())
2064 for (
int ind = 0; ind < this->nElements(); ind++) {
2065 this->elements[ind].
2066 frobThreshLowestLevel(threshold, &ErrorMatrix->
elements[ind]);
2067 thisMatIsZero = thisMatIsZero && this->elements[ind].is_zero();
2068 EMatIsZero = EMatIsZero && ErrorMatrix->
elements[ind].is_zero();
2073 ErrorMatrix->
clear();
2077 if (!this->is_zero()) {
2078 for (
int ind = 0; ind < this->nElements(); ind++) {
2079 this->elements[ind].frobThreshLowestLevel(threshold, 0);
2080 thisMatIsZero = thisMatIsZero && this->elements[ind].is_zero();
2087 template<
class Treal,
class Telement>
2090 if (!this->is_zero())
2091 for (
int ind = 0; ind < this->nElements(); ind++)
2092 this->elements[ind].getFrobSqElementLevel(frobsq);
2095 template<
class Treal,
class Telement>
2099 assert(!this->is_empty());
2100 bool thisMatIsZero =
true;
2102 assert(!ErrorMatrix->is_empty());
2103 bool EMatIsZero =
true;
2104 if (!ErrorMatrix->is_zero() || !this->is_zero()) {
2105 if (ErrorMatrix->is_zero())
2107 if (this->is_zero())
2109 for (
int ind = 0; ind < this->nElements(); ind++) {
2110 this->elements[ind].
2111 frobThreshElementLevel(threshold, &ErrorMatrix->
elements[ind]);
2112 thisMatIsZero = thisMatIsZero && this->elements[ind].is_zero();
2113 EMatIsZero = EMatIsZero && ErrorMatrix->
elements[ind].is_zero();
2118 ErrorMatrix->
clear();
2122 if (!this->is_zero()) {
2123 for (
int ind = 0; ind < this->nElements(); ind++) {
2124 this->elements[ind].frobThreshElementLevel(threshold, 0);
2125 thisMatIsZero = thisMatIsZero && this->elements[ind].is_zero();
2134 template<
class Treal,
class Telement>
2138 if ( this->is_zero() )
2140 assert( this->nElements() == A.nElements() );
2141 for (
int ind = 0; ind < this->nElements(); ind++)
2142 this->elements[ind].assignFrobNormsLowestLevel(A[ind]);
2148 template<
class Treal,
class Telement>
2151 assert(!A.is_empty());
2153 throw Failure(
"Matrix<Treal, Telement>::syAssignFrobNormsLowestLevel(...): "
2154 "Matrix must be have equal number of rows "
2155 "and cols to be symmetric");
2157 if ( this->is_zero() )
2159 assert( this->nElements() == A.nElements() );
2160 for (
int col = 1; col < this->ncols(); col++)
2161 for (
int row = 0; row < col; row++)
2162 (*
this)(row, col).assignFrobNormsLowestLevel( A(row,col) );
2163 for (
int rc = 0; rc < this->ncols(); rc++)
2164 (*
this)(rc, rc).syAssignFrobNormsLowestLevel( A(rc,rc) );
2170 template<
class Treal,
class Telement>
2174 if ( A.is_zero() && B.is_zero() ) {
2180 if ( this->is_zero() )
2182 if ( !A.is_zero() && !B.is_zero() ) {
2184 assert( this->nElements() == A.nElements() );
2185 assert( this->nElements() == B.nElements() );
2186 for (
int ind = 0; ind < this->nElements(); ind++)
2187 this->elements[ind].assignDiffFrobNormsLowestLevel( A[ind], B[ind] );
2190 if ( !A.is_zero() ) {
2192 this->assignFrobNormsLowestLevel( A );
2195 if ( !B.is_zero() ) {
2197 this->assignFrobNormsLowestLevel( B );
2201 template<
class Treal,
class Telement>
2205 if ( A.is_zero() && B.is_zero() ) {
2211 if ( this->is_zero() )
2213 if ( !A.is_zero() && !B.is_zero() ) {
2215 assert( this->nElements() == A.nElements() );
2216 assert( this->nElements() == B.nElements() );
2217 for (
int col = 1; col < this->ncols(); col++)
2218 for (
int row = 0; row < col; row++)
2219 (*
this)(row, col).assignDiffFrobNormsLowestLevel( A(row,col), B(row,col) );
2220 for (
int rc = 0; rc < this->ncols(); rc++)
2221 (*
this)(rc, rc).syAssignDiffFrobNormsLowestLevel( A(rc,rc), B(rc,rc) );
2224 if ( !A.is_zero() ) {
2226 this->syAssignFrobNormsLowestLevel( A );
2229 if ( !B.is_zero() ) {
2231 this->syAssignFrobNormsLowestLevel( B );
2238 template<
class Treal,
class Telement>
2241 if ( this->is_zero() ) {
2245 assert( !A.is_zero() );
2246 assert( this->nElements() == A.nElements() );
2247 for (
int ind = 0; ind < this->nElements(); ind++)
2248 this->elements[ind].truncateAccordingToSparsityPattern( A[ind] );
2258 template<
class Treal,
class Telement>
2265 assert(!A.is_empty());
2266 assert(!B.is_empty());
2270 C.resetRows(A.
cols);
2272 C.resetRows(A.
rows);
2274 C.resetCols(B.
rows);
2276 C.resetCols(B.
cols);
2278 if ((A.is_zero() || B.is_zero() || alpha == 0) &&
2279 (C.is_zero() || beta == 0))
2282 Treal beta_tmp = beta;
2287 if (!A.is_zero() && !B.is_zero() && alpha != 0) {
2292 for (
int col = 0; col < C.
ncols(); col++) {
2293 Telement::gemm_upper_tr_only(tA, tB, alpha,
2294 A(col, 0), B(0, col),
2297 for(
int cola = 1; cola < A.
ncols(); cola++)
2298 Telement::gemm_upper_tr_only(tA, tB, alpha,
2299 A(col, cola), B(cola, col),
2302 for (
int row = 0; row < col; row++) {
2304 A(row, 0), B(0, col),
2307 for(
int cola = 1; cola < A.
ncols(); cola++)
2309 A(row, cola), B(cola, col),
2316 throw Failure(
"Matrix<class Treal, class Telement>::"
2317 "gemm_upper_tr_only(bool, bool, Treal, "
2318 "Matrix<Treal, Telement>, "
2319 "Matrix<Treal, Telement>, "
2320 "Treal, Matrix<Treal, Telement>): "
2321 "Incorrect matrixdimensions for multiplication");
2326 for (
int col = 0; col < C.
ncols(); col++) {
2327 Telement::gemm_upper_tr_only(tA, tB, alpha,
2331 for(
int cola = 1; cola < A.
nrows(); cola++)
2332 Telement::gemm_upper_tr_only(tA, tB, alpha,
2333 A(cola, col), B(cola, col),
2336 for (
int row = 0; row < col; row++) {
2341 for(
int cola = 1; cola < A.
nrows(); cola++)
2343 A(cola, row), B(cola, col),
2350 throw Failure(
"Matrix<class Treal, class Telement>::"
2351 "gemm_upper_tr_only(bool, bool, Treal, "
2352 "Matrix<Treal, Telement>, "
2353 "Matrix<Treal, Telement>, "
2354 "Treal, Matrix<Treal, Telement>): "
2355 "Incorrect matrixdimensions for multiplication");
2360 for (
int col = 0; col < C.
ncols(); col++) {
2361 Telement::gemm_upper_tr_only(tA, tB, alpha,
2362 A(col, 0), B(col, 0),
2365 for(
int cola = 1; cola < A.
ncols(); cola++)
2366 Telement::gemm_upper_tr_only(tA, tB, alpha,
2367 A(col, cola), B(col, cola),
2370 for (
int row = 0; row < col; row++) {
2372 A(row, 0), B(col, 0),
2375 for(
int cola = 1; cola < A.
ncols(); cola++)
2377 A(row, cola), B(col, cola),
2384 throw Failure(
"Matrix<class Treal, class Telement>::"
2385 "gemm_upper_tr_only(bool, bool, Treal, "
2386 "Matrix<Treal, Telement>, "
2387 "Matrix<Treal, Telement>, "
2388 "Treal, Matrix<Treal, Telement>): "
2389 "Incorrect matrixdimensions for multiplication");
2394 for (
int col = 0; col < C.
ncols(); col++) {
2395 Telement::gemm_upper_tr_only(tA, tB, alpha,
2396 A(0, col), B(col, 0),
2399 for(
int cola = 1; cola < A.
nrows(); cola++)
2400 Telement::gemm_upper_tr_only(tA, tB, alpha,
2401 A(cola, col), B(col, cola),
2404 for (
int row = 0; row < col; row++) {
2406 A(0, row), B(col, 0),
2409 for(
int cola = 1; cola < A.
nrows(); cola++)
2411 A(cola, row), B(col, cola),
2418 throw Failure(
"Matrix<class Treal, class Telement>::"
2419 "gemm_upper_tr_only(bool, bool, Treal, "
2420 "Matrix<Treal, Telement>, "
2421 "Matrix<Treal, Telement>, Treal, "
2422 "Matrix<Treal, Telement>): "
2423 "Incorrect matrixdimensions for multiplication");
2424 else throw Failure(
"Matrix<class Treal, class Telement>::"
2425 "gemm_upper_tr_only(bool, bool, Treal, "
2426 "Matrix<Treal, Telement>, "
2427 "Matrix<Treal, Telement>, Treal, "
2428 "Matrix<Treal, Telement>):"
2429 "Very strange error!!");
2441 template<
class Treal,
class Telement>
2446 assert(!A.is_empty());
2447 assert(!Z.is_empty());
2448 if (alpha != 0 && !A.is_zero() && !Z.is_zero()) {
2454 for (
int col = A.
ncols() - 1; col >= 0; col--) {
2456 Telement::sytr_upper_tr_only(side, alpha,
2457 A(col, col), Z(col, col));
2458 for (
int ind = 0; ind < col; ind++) {
2460 Telement::gemm_upper_tr_only(
true,
false, alpha, A(ind, col),
2461 Z(ind, col), 1.0, A(col, col));
2464 for (
int row = col - 1; row >= 0; row--) {
2467 alpha, Z(col, col), A(row, col));
2471 for (
int ind = 0; ind < row; ind++) {
2476 for (
int ind = row + 1; ind < col; ind++) {
2485 assert(side ==
'L');
2487 for (
int col = 0; col < A.
ncols(); col++) {
2489 for (
int row = 0; row < col; row++) {
2492 Z(row, row), A(row, col));
2496 for (
int ind = row + 1; ind < col; ind++)
2500 for (
int ind = col + 1; ind < A.
nrows(); ind++)
2505 Telement::sytr_upper_tr_only(side, alpha,
2506 A(col, col), Z(col, col));
2507 for (
int ind = col + 1; ind < A.
ncols(); ind++)
2508 Telement::gemm_upper_tr_only(
false,
true, alpha, Z(col, ind),
2509 A(col, ind), 1.0, A(col, col));
2514 throw Failure(
"Matrix<class Treal, class Telement>::"
2515 "sytr_upper_tr_only: Incorrect matrix dimensions "
2516 "for symmetric-triangular multiplication");
2524 template<
class Treal,
class Telement>
2527 const bool tA,
const Treal alpha,
2530 assert(!A.is_empty());
2531 assert(!B.is_empty());
2532 if (alpha != 0 && !A.is_zero() && !B.is_zero())
2533 if (((side ==
'R' && B.
ncols() == A.
nrows()) ||
2538 throw Failure(
"Matrix<Treal, class Telement>::"
2539 "trmm_upper_tr_only : "
2540 "not possible for upper triangular not transposed "
2541 "matrices A since lower triangle of B is needed");
2547 for (
int col = 0; col < B.
ncols(); col++) {
2548 Telement::trmm_upper_tr_only(side, uplo, tA, alpha,
2549 A(col,col), B(col,col));
2550 for (
int ind = col + 1; ind < A.
ncols(); ind++)
2551 Telement::gemm_upper_tr_only(
false, tA, alpha,
2552 B(col,ind), A(col,ind),
2554 for (
int row = 0; row < col; row++) {
2556 A(col,col), B(row,col));
2557 for (
int ind = col + 1; ind < A.
ncols(); ind++)
2559 B(row,ind), A(col,ind),
2565 assert(side ==
'L');
2567 for (
int row = B.
nrows() - 1; row >= 0; row--) {
2568 Telement::trmm_upper_tr_only(side, uplo, tA, alpha,
2569 A(row,row), B(row,row));
2570 for (
int ind = 0; ind < row; ind++)
2571 Telement::gemm_upper_tr_only(tA,
false, alpha,
2572 A(ind,row), B(ind,row),
2574 for (
int col = row + 1; col < B.
ncols(); col++) {
2576 A(row,row), B(row,col));
2577 for (
int ind = 0; ind < row; ind++)
2579 A(ind,row), B(ind,col),
2586 throw Failure(
"Matrix<class Treal, class Telement>::"
2587 "trmm_upper_tr_only not implemented for lower "
2588 "triangular matrices");
2590 throw Failure(
"Matrix<class Treal, class Telement>::"
2591 "trmm_upper_tr_only: Incorrect matrix dimensions "
2592 "for multiplication");
2601 template<
class Treal,
class Telement>
2613 assert(side ==
'L');
2621 template<
class Treal,
class Telement>
2625 assert(!this->is_empty());
2626 if (ErrorMatrix && ErrorMatrix->is_empty()) {
2627 ErrorMatrix->resetRows(this->rows);
2628 ErrorMatrix->resetCols(this->cols);
2630 assert(threshold >= (Treal)0.0);
2631 if (threshold == (Treal)0.0)
2634 std::vector<Treal> frobsq_norms;
2635 getFrobSqLowestLevel(frobsq_norms);
2637 std::sort(frobsq_norms.begin(), frobsq_norms.end());
2638 int lastRemoved = -1;
2639 Treal frobsqSum = 0;
2640 int nnzBlocks = frobsq_norms.size();
2641 while(lastRemoved < nnzBlocks - 1 && frobsqSum < threshold) {
2643 frobsqSum += frobsq_norms[lastRemoved];
2647 if (lastRemoved == nnzBlocks - 1 && frobsqSum < threshold) {
2654 frobsqSum -= frobsq_norms[lastRemoved];
2656 while(lastRemoved > -1 && frobsq_norms[lastRemoved] ==
2657 frobsq_norms[lastRemoved + 1]) {
2658 frobsqSum -= frobsq_norms[lastRemoved];
2661 if (lastRemoved > -1) {
2662 Treal threshLowestLevel =
2663 (frobsq_norms[lastRemoved + 1] +
2664 frobsq_norms[lastRemoved]) / 2;
2665 this->frobThreshLowestLevel(threshLowestLevel, ErrorMatrix);
2671 template<
class Treal,
class Telement>
2675 const Treal threshold,
const side looking,
2677 assert(!A.is_empty());
2679 throw Failure(
"Matrix<Treal, Telement>::sy_inch: "
2680 "Matrix must be quadratic!");
2682 throw Failure(
"Matrix<Treal>::sy_inch: Matrix is zero! "
2683 "Not possible to compute inverse cholesky.");
2685 throw Failure(
"Matrix<Treal>::sy_inch: Only unstable "
2686 "version of sy_inch implemented.");
2689 myThresh = threshold / (Z.
ncols() * Z.
nrows());
2690 Z.resetRows(A.
rows);
2691 Z.resetCols(A.
cols);
2694 if (looking ==
left) {
2696 throw Failure(
"Matrix<Treal, Telement>::syInch: "
2697 "Thresholding not implemented for left-looking inch.");
2699 Telement::syInch(A(0,0), Z(0,0), threshold, looking, version);
2701 for (
int i = 1; i < Z.
ncols(); i++) {
2702 for (
int j = 0; j < i; j++) {
2706 for (
int k = 0; k < j; k++)
2708 A(k,j), Z(k,i), 1.0, Ptmp);
2711 for (
int k = 0; k < j; k++)
2713 Z(k,j), Ptmp, 1.0, Z(k,i));
2716 Telement::add(1.0, Ptmp, Z(j,i));
2719 for (
int k = 0; k < i; k++)
2720 Telement::gemm_upper_tr_only(
true,
false, 1.0,
2721 A(k,i), Z(k,i), 1.0, Ptmp);
2724 Telement::syInch(Ptmp, Z(i,i), threshold, looking, version);
2725 for (
int k = 0; k < i; k++) {
2732 assert(looking ==
right);
2734 Treal newThresh = 0;
2735 if (myThresh != 0 && Z.
ncols() > 1)
2736 newThresh = myThresh * 0.0001;
2738 newThresh = myThresh;
2740 for (
int i = 0; i < Z.
ncols(); i++) {
2743 for (
int k = 0; k < i; k++)
2744 Telement::gemm_upper_tr_only(
true,
false, 1.0,
2745 A(k,i), Z(k,i), 1.0, Ptmp);
2748 Telement::syInch(Ptmp, Z(i,i), newThresh, looking, version);
2749 for (
int k = 0; k < i; k++) {
2754 for (
int j = i + 1; j < Z.
ncols(); j++) {
2759 for (
int k = 0; k < i; k++)
2761 A(k,i), Z(k,j), 1.0, Ptmp);
2764 for (
int k = 0; k < i; k++)
2766 Z(k,i), Ptmp, 1.0, Z(k,j));
2769 Telement::add(1.0, Ptmp, Z(i,j));
2773 if (threshold != 0) {
2774 for (
int k = 0; k < i; k++)
2781 template<
class Treal,
class Telement>
2784 assert(!this->is_empty());
2785 if (this->nScalarsRows() == this->nScalarsCols()) {
2786 int size = this->nScalarsCols();
2787 Treal* colsums =
new Treal[size];
2788 Treal* diag =
new Treal[size];
2789 for (
int ind = 0; ind < size; ind++)
2791 this->add_abs_col_sums(colsums);
2792 this->get_diagonal(diag);
2795 lmin = diag[0] - tmp1;
2796 lmax = diag[0] + tmp1;
2797 for (
int col = 1; col < size; col++) {
2799 tmp2 = diag[col] - tmp1;
2800 lmin = (tmp2 < lmin ? tmp2 : lmin);
2801 tmp2 = diag[col] + tmp1;
2802 lmax = (tmp2 > lmax ? tmp2 : lmax);
2808 throw Failure(
"Matrix<Treal, Telement, int>::gersgorin(Treal, Treal): "
2809 "Matrix must be quadratic");
2813 template<
class Treal,
class Telement>
2817 if (!this->is_zero()) {
2819 for (
int col = 0; col < this->ncols(); col++) {
2820 for (
int row = 0; row < this->nrows(); row++) {
2821 (*this)(row,col).add_abs_col_sums(&sums[offset]);
2823 offset += (*this)(0,col).nScalarsCols();
2828 template<
class Treal,
class Telement>
2832 assert(this->nrows() == this->ncols());
2833 if (this->is_zero())
2834 for (
int rc = 0; rc < this->nScalarsCols(); rc++)
2838 for (
int rc = 0; rc < this->ncols(); rc++) {
2839 (*this)(rc,rc).get_diagonal(&diag[offset]);
2840 offset += (*this)(rc,rc).nScalarsCols();
2845 template<
class Treal,
class Telement>
2847 assert(!this->is_empty());
2849 if (this->is_zero())
2851 for (
int ind = 0; ind < this->nElements(); ind++)
2852 sum += this->elements[ind].memory_usage();
2856 template<
class Treal,
class Telement>
2859 if (!this->is_zero()) {
2860 for (
int ind = 0; ind < this->nElements(); ind++)
2861 sum += this->elements[ind].nnz();
2865 template<
class Treal,
class Telement>
2868 if (!this->is_zero()) {
2870 for (
int col = 1; col < this->ncols(); col++)
2871 for (
int row = 0; row < col; row++)
2872 sum += (*
this)(row, col).nnz();
2876 for (
int rc = 0; rc < this->nrows(); rc++)
2877 sum += (*
this)(rc,rc).sy_nnz();
2882 template<
class Treal,
class Telement>
2885 if (!this->is_zero()) {
2887 for (
int col = 1; col < this->ncols(); col++)
2888 for (
int row = 0; row < col; row++)
2889 sum += (*
this)(row, col).nvalues();
2891 for (
int rc = 0; rc < this->nrows(); rc++)
2892 sum += (*
this)(rc,rc).sy_nvalues();
2907 template<
class Treal>
2924 assert(!this->is_empty());
2926 this->
elements =
new Treal[this->nElements()];
2927 for (
int ind = 0; ind < this->nElements(); ++ind)
2933 void fullMatrix(std::vector<Treal> & fullMat)
const;
2939 std::vector<int>
const & colind,
2940 std::vector<Treal>
const & values);
2942 std::vector<int>
const & colind,
2943 std::vector<Treal>
const & values,
2944 std::vector<int>
const & indexes);
2947 void addValues(std::vector<int>
const & rowind,
2948 std::vector<int>
const & colind,
2949 std::vector<Treal>
const & values);
2950 void addValues(std::vector<int>
const & rowind,
2951 std::vector<int>
const & colind,
2952 std::vector<Treal>
const & values,
2953 std::vector<int>
const & indexes);
2956 std::vector<int>
const & colind,
2957 std::vector<Treal>
const & values);
2960 std::vector<int>
const & colind,
2961 std::vector<Treal>
const & values);
2963 void getValues(std::vector<int>
const & rowind,
2964 std::vector<int>
const & colind,
2965 std::vector<Treal> & values)
const;
2966 void getValues(std::vector<int>
const & rowind,
2967 std::vector<int>
const & colind,
2968 std::vector<Treal> & values,
2969 std::vector<int>
const & indexes)
const;
2971 std::vector<int>
const & colind,
2972 std::vector<Treal> & values)
const;
2975 std::vector<int> & colind,
2976 std::vector<Treal> & values)
const;
2978 std::vector<int> & colind,
2979 std::vector<Treal> & values)
const;
2999 template<
typename TRule>
3001 template<
typename TRule>
3018 static void gemm(
const bool tA,
const bool tB,
const Treal alpha,
3023 static void symm(
const char side,
const char uplo,
const Treal alpha,
3028 static void syrk(
const char uplo,
const bool tA,
const Treal alpha,
3033 static void sysq(
const char uplo,
const Treal alpha,
3038 static void ssmm(
const Treal alpha,
3052 static void trmm(
const char side,
const char uplo,
const bool tA,
3083 Treal
trace()
const;
3091 static void add(
const Treal alpha,
3094 void assign(Treal
const alpha,
3111 (Treal
const threshold,
3121 if ( this->is_zero() )
3123 assert( this->nElements() == A.nElements() );
3124 for (
int ind = 0; ind < this->nElements(); ind++)
3133 if ( this->is_zero() )
3135 assert( this->nElements() == A.nElements() );
3136 for (
int col = 1; col < this->
ncols(); col++)
3137 for (
int row = 0; row < col; row++)
3138 (*
this)(row,col) = A(row,col).
frob();
3139 for (
int rc = 0; rc < this->
nrows(); rc++)
3140 (*
this)(rc,rc) = A(rc,rc).
syFrob();
3148 if ( A.is_zero() && B.is_zero() ) {
3154 if ( this->is_zero() )
3156 if ( !A.is_zero() && !B.is_zero() ) {
3158 assert( this->nElements() == A.nElements() );
3159 assert( this->nElements() == B.nElements() );
3160 for (
int ind = 0; ind < this->nElements(); ind++) {
3167 if ( !A.is_zero() ) {
3172 if ( !B.is_zero() ) {
3180 if ( A.is_zero() && B.is_zero() ) {
3186 if ( this->is_zero() )
3188 if ( !A.is_zero() && !B.is_zero() ) {
3190 assert( this->nElements() == A.nElements() );
3191 assert( this->nElements() == B.nElements() );
3192 for (
int col = 1; col < this->
ncols(); col++)
3193 for (
int row = 0; row < col; row++) {
3196 (*this)(row, col) = Diff.
frob();
3198 for (
int rc = 0; rc < this->
ncols(); rc++) {
3201 (*this)(rc, rc) = Diff.
syFrob();
3205 if ( !A.is_zero() ) {
3210 if ( !B.is_zero() ) {
3219 if ( this->is_zero() )
3223 assert( this->nElements() == A.nElements() );
3224 for (
int ind = 0; ind < this->nElements(); ind++)
3243 const bool tA,
const Treal alpha,
3263 const Treal threshold = 0,
3268 const Treal threshold = 0,
3271 inch(A, Z, threshold, looking, version);
3274 void gersgorin(Treal& lmin, Treal& lmax)
const;
3285 assert(!this->is_empty());
3286 if (this->is_zero())
3289 return (
size_t)this->nElements() *
sizeof(Treal);
3293 if (this->is_zero())
3296 return this->nElements();
3299 if (this->is_zero())
3302 return this->nElements();
3314 int n = this->
nrows();
3315 return this->is_zero() ? 0 : n * (n + 1) / 2;
3324 if (!this->is_zero()) {
3327 for (
int col = 0; col < this->
ncols(); col++) {
3328 for (
int row = 0; row < col; row++) {
3329 res += 2*op.accumulate((*
this)(row, col),
3333 res += op.accumulate((*
this)(col, col),
3343 if (!this->is_zero()) {
3346 for (
int col = 0; col < this->
ncols(); col++)
3347 for (
int row = 0; row < this->
nrows(); row++)
3348 res += op.accumulate((*
this)(row, col),
3355 static inline unsigned int level() {
return 0;}
3358 if (this->is_zero())
3361 Treal maxAbsGlobal = 0;
3362 Treal maxAbsLocal = 0;
3363 for (
int ind = 0; ind < this->nElements(); ++ind) {
3365 maxAbsGlobal = maxAbsGlobal > maxAbsLocal ?
3366 maxAbsGlobal : maxAbsLocal;
3368 return maxAbsGlobal;
3380 std::cout<<
"operator="<<std::endl;
3393 void trim_memory_usage();
3395 void write_to_buffer_count(
int& zb_length,
int& vb_length)
const;
3396 void write_to_buffer(
int* zerobuf,
const int zb_length,
3397 Treal* valuebuf,
const int vb_length,
3398 int& zb_index,
int& vb_index)
const;
3399 void read_from_buffer(
int* zerobuf,
const int zb_length,
3400 Treal* valuebuf,
const int vb_length,
3401 int& zb_index,
int& vb_index);
3418 inline bool lowestLevel()
const {
return true;}
3428 template<
class Treal>
3430 template<
class Treal>
3435 template<
class Treal>
3440 assert((
int)fullMat.size() == nTotalRows * nTotalCols);
3443 if (this->is_zero())
3445 for (
int col = 0; col < this->
ncols(); col++)
3446 for (
int row = 0; row < this->
nrows(); row++)
3448 fullMat[(colOffset + col) * nTotalRows + (rowOffset + row)];
3451 template<
class Treal>
3456 fullMat.resize(nTotalRows * nTotalCols);
3459 if (this->is_zero()) {
3462 fullMat[(colOffset + col) * nTotalRows + (rowOffset + row)] = 0;
3465 for (
int col = 0; col < this->
ncols(); col++)
3466 for (
int row = 0; row < this->
nrows(); row++)
3467 fullMat[(colOffset + col) * nTotalRows + (rowOffset + row)] =
3472 template<
class Treal>
3477 fullMat.resize(nTotalRows * nTotalCols);
3480 if (this->is_zero()) {
3482 for (
int row = 0; row <= col; row++) {
3483 fullMat[(colOffset + col) * nTotalRows + (rowOffset + row)] = 0;
3484 fullMat[(rowOffset + row) * nTotalRows + (colOffset + col)] = 0;
3488 for (
int col = 0; col < this->
ncols(); col++)
3489 for (
int row = 0; row <= col; row++) {
3490 fullMat[(colOffset + col) * nTotalRows + (rowOffset + row)] =
3492 fullMat[(rowOffset + row) * nTotalRows + (colOffset + col)] =
3498 template<
class Treal>
3503 fullMat.resize(nTotalRows * nTotalCols);
3506 if (this->is_zero()) {
3508 for (
int row = 0; row <= this->
nScalarsRows(); row++) {
3509 fullMat[(colOffset + col) * nTotalRows + (rowOffset + row)] = 0;
3510 fullMat[(rowOffset + row) * nTotalRows + (colOffset + col)] = 0;
3514 for (
int col = 0; col < this->
ncols(); col++)
3515 for (
int row = 0; row < this->
nrows(); row++) {
3516 fullMat[(colOffset + col) * nTotalRows + (rowOffset + row)] =
3518 fullMat[(rowOffset + row) * nTotalRows + (colOffset + col)] =
3526 template<
class Treal>
3529 std::vector<int>
const & colind,
3530 std::vector<Treal>
const & values) {
3531 assert(rowind.size() == colind.size() &&
3532 rowind.size() == values.size());
3533 std::vector<int> indexes(values.size());
3534 for (
int ind = 0; ind < values.size(); ++ind) {
3540 template<
class Treal>
3543 std::vector<int>
const & colind,
3544 std::vector<Treal>
const & values,
3545 std::vector<int>
const & indexes) {
3546 if (indexes.empty()) {
3550 if (this->is_zero())
3552 for (
int ind = 0; ind < this->nElements(); ++ind)
3554 std::vector<int>::const_iterator it;
3555 for ( it = indexes.begin(); it < indexes.end(); it++ ) {
3562 template<
class Treal>
3565 std::vector<int>
const & colind,
3566 std::vector<Treal>
const & values) {
3567 assert(rowind.size() == colind.size() &&
3568 rowind.size() == values.size());
3569 std::vector<int> indexes(values.size());
3570 for (
int ind = 0; ind < values.size(); ++ind) {
3573 addValues(rowind, colind, values, indexes);
3576 template<
class Treal>
3579 std::vector<int>
const & colind,
3580 std::vector<Treal>
const & values,
3581 std::vector<int>
const & indexes) {
3582 if (indexes.empty())
3584 if (this->is_zero())
3586 std::vector<int>::const_iterator it;
3587 for ( it = indexes.begin(); it < indexes.end(); it++ ) {
3593 template<
class Treal>
3596 std::vector<int>
const & colind,
3597 std::vector<Treal>
const & values) {
3598 assert(rowind.size() == colind.size() &&
3599 rowind.size() == values.size());
3600 bool upperTriangleOnly =
true;
3601 for (
int ind = 0; ind < values.size(); ++ind) {
3603 upperTriangleOnly && colind[ind] >= rowind[ind];
3605 if (!upperTriangleOnly)
3606 throw Failure(
"Matrix<Treal>::"
3607 "syAddValues(std::vector<int> const &, "
3608 "std::vector<int> const &, "
3609 "std::vector<Treal> const &, int const): "
3610 "Only upper triangle can contain elements when assigning "
3611 "symmetric or triangular matrix ");
3615 template<
class Treal>
3618 std::vector<int>
const & colind,
3619 std::vector<Treal>
const & values) {
3620 assert(rowind.size() == colind.size() &&
3621 rowind.size() == values.size());
3622 bool upperTriangleOnly =
true;
3623 for (
int ind = 0; ind < values.size(); ++ind) {
3625 upperTriangleOnly && colind[ind] >= rowind[ind];
3627 if (!upperTriangleOnly)
3628 throw Failure(
"Matrix<Treal>::"
3629 "syAddValues(std::vector<int> const &, "
3630 "std::vector<int> const &, "
3631 "std::vector<Treal> const &, int const): "
3632 "Only upper triangle can contain elements when assigning "
3633 "symmetric or triangular matrix ");
3637 template<
class Treal>
3640 std::vector<int>
const & colind,
3641 std::vector<Treal> & values)
const {
3642 assert(rowind.size() == colind.size());
3643 values.resize(rowind.size(),0);
3644 std::vector<int> indexes(rowind.size());
3645 for (
int ind = 0; ind < rowind.size(); ++ind) {
3648 getValues(rowind, colind, values, indexes);
3651 template<
class Treal>
3654 std::vector<int>
const & colind,
3655 std::vector<Treal> & values,
3656 std::vector<int>
const & indexes)
const {
3657 assert(!this->is_empty());
3658 if (indexes.empty())
3660 std::vector<int>::const_iterator it;
3661 if (this->is_zero()) {
3662 for ( it = indexes.begin(); it < indexes.end(); it++ )
3666 for ( it = indexes.begin(); it < indexes.end(); it++ )
3672 template<
class Treal>
3675 std::vector<int>
const & colind,
3676 std::vector<Treal> & values)
const {
3677 assert(rowind.size() == colind.size());
3678 bool upperTriangleOnly =
true;
3679 for (
int ind = 0; ind < rowind.size(); ++ind) {
3681 upperTriangleOnly && colind[ind] >= rowind[ind];
3683 if (!upperTriangleOnly)
3684 throw Failure(
"Matrix<Treal>::"
3685 "syGetValues(std::vector<int> const &, "
3686 "std::vector<int> const &, "
3687 "std::vector<Treal> const &, int const): "
3688 "Only upper triangle when retrieving elements from "
3689 "symmetric or triangular matrix ");
3693 template<
class Treal>
3696 std::vector<int> & colind,
3697 std::vector<Treal> & values)
const {
3698 assert(rowind.size() == colind.size() &&
3699 rowind.size() == values.size());
3700 if (!this->is_zero())
3701 for (
int col = 0; col < this->
ncols(); col++)
3702 for (
int row = 0; row < this->
nrows(); row++) {
3705 values.push_back((*
this)(row, col));
3709 template<
class Treal>
3712 std::vector<int> & colind,
3713 std::vector<Treal> & values)
const {
3714 assert(rowind.size() == colind.size() &&
3715 rowind.size() == values.size());
3716 if (!this->is_zero())
3717 for (
int col = 0; col < this->
ncols(); ++col)
3718 for (
int row = 0; row <= col; ++row) {
3721 values.push_back((*
this)(row, col));
3726 template<
class Treal>
3732 template<
class Treal>
3737 if (this->is_zero()) {
3738 char * tmp = (
char*)&ZERO;
3739 file.write(tmp,
sizeof(
int));
3742 char * tmp = (
char*)&ONE;
3743 file.write(tmp,
sizeof(
int));
3744 char * tmpel = (
char*)this->
elements;
3745 file.write(tmpel,
sizeof(Treal) * this->nElements());
3749 template<
class Treal>
3754 char tmp[
sizeof(int)];
3755 file.read(tmp, (std::ifstream::pos_type)
sizeof(
int));
3756 switch ((
int)*tmp) {
3761 if (this->is_zero())
3763 file.read((
char*)this->
elements,
sizeof(Treal) * this->nElements());
3766 throw Failure(
"Matrix<Treal>::"
3767 "readFromFile(std::ifstream & file):"
3768 "File corruption, int value not 0 or 1");
3772 template<
class Treal>
3774 if (this->is_zero())
3776 for (
int ind = 0; ind < this->nElements(); ind++)
3777 this->
elements[ind] = rand() / (Treal)RAND_MAX;
3780 template<
class Treal>
3782 if (this->is_zero())
3785 for (
int col = 1; col < this->
ncols(); col++)
3786 for (
int row = 0; row < col; row++)
3787 (*
this)(row, col) = rand() / (Treal)RAND_MAX;;
3789 for (
int rc = 0; rc < this->
nrows(); rc++)
3790 (*
this)(rc,rc) = rand() / (Treal)RAND_MAX;;
3793 template<
class Treal>
3796 if (!this->highestLevel() &&
3797 probabilityBeingZero > rand() / (Treal)RAND_MAX)
3803 template<
class Treal>
3806 if (!this->highestLevel() &&
3807 probabilityBeingZero > rand() / (Treal)RAND_MAX)
3813 template<
class Treal>
3814 template<
typename TRule>
3817 if (this->is_zero())
3819 for (
int col = 0; col < this->
ncols(); col++)
3820 for (
int row = 0; row < this->
nrows(); row++)
3821 (*
this)(row,col) = rule.set(this->rows.getOffset() + row,
3824 template<
class Treal>
3825 template<
typename TRule>
3828 if (this->is_zero())
3831 for (
int col = 0; col < this->
ncols(); col++)
3832 for (
int row = 0; row <= col; row++)
3833 (*
this)(row,col) = rule.set(this->rows.getOffset() + row,
3838 template<
class Treal>
3841 if (this->is_empty())
3842 throw Failure(
"Matrix<Treal>::addIdentity(Treal): "
3843 "Cannot add identity to empty matrix.");
3845 throw Failure(
"Matrix<Treal, Telement>::addIdentity(Treal): "
3846 "Matrix must be square to add identity");
3847 if (this->is_zero())
3849 for (
int ind = 0; ind < this->
ncols(); ind++)
3850 (*
this)(ind,ind) += alpha;
3853 template<
class Treal>
3863 if (AT.is_zero() || (AT.nElements() != A.nElements())) {
3865 AT.
elements =
new Treal[A.nElements()];
3869 for (
int row = 0; row < AT.
nrows(); row++)
3870 for (
int col = 0; col < AT.
ncols(); col++)
3871 AT(row,col) = A(col,row);
3874 template<
class Treal>
3878 if (!this->is_zero()) {
3881 for (
int row = 1; row < this->
nrows(); row++)
3882 for (
int col = 0; col < row; col++)
3883 (*
this)(row, col) = (*
this)(col, row);
3887 throw Failure(
"Matrix<Treal>::symToNosym(): "
3888 "Only quadratic matrices can be symmetric");
3891 template<
class Treal>
3895 if (!this->is_zero()) {
3898 for (
int col = 0; col < this->
ncols() - 1; col++)
3899 for (
int row = col + 1; row < this->
nrows(); row++)
3900 (*
this)(row, col) = 0;
3904 throw Failure(
"Matrix<Treal>::nosymToSym(): "
3905 "Only quadratic matrices can be symmetric");
3908 template<
class Treal>
3917 throw Failure(
"Matrix<Treal>::operator=(int k = 1): "
3918 "Matrix must be quadratic to "
3919 "become identity matrix.");
3922 for (
int rc = 0; rc < this->
ncols(); rc++)
3927 (
"Matrix<Treal>::operator=(int k) only implemented for k = 0, k = 1");
3932 template<
class Treal>
3935 if (!this->is_zero() && alpha != 1) {
3936 for (
int ind = 0; ind < this->nElements(); ind++)
3942 template<
class Treal>
3944 gemm(
const bool tA,
const bool tB,
const Treal alpha,
3954 C.resetRows(A.
cols);
3956 C.resetRows(A.
rows);
3958 C.resetCols(B.
rows);
3960 C.resetCols(B.
cols);
3963 if ((A.is_zero() || B.is_zero() || alpha == 0) &&
3964 (C.is_zero() || beta == 0))
3967 Treal beta_tmp = beta;
3973 if (!A.is_zero() && !B.is_zero() && alpha != 0) {
3982 throw Failure(
"Matrix<Treal>::"
3983 "gemm(bool, bool, Treal, Matrix<Treal>, "
3984 "Matrix<Treal>, Treal, "
3986 "Incorrect matrixdimensions for multiplication");
3995 throw Failure(
"Matrix<Treal>::"
3996 "gemm(bool, bool, Treal, Matrix<Treal>, "
3997 "Matrix<Treal>, Treal, "
3999 "Incorrect matrixdimensions for multiplication");
4008 throw Failure(
"Matrix<Treal>::"
4009 "gemm(bool, bool, Treal, Matrix<Treal>, "
4010 "Matrix<Treal>, Treal, "
4012 "Incorrect matrixdimensions for multiplication");
4021 throw Failure(
"Matrix<Treal>::"
4022 "gemm(bool, bool, Treal, Matrix<Treal>, "
4023 "Matrix<Treal>, Treal, "
4025 "Incorrect matrixdimensions for multiplication");
4026 else throw Failure(
"Matrix<Treal>::"
4027 "gemm(bool, bool, Treal, Matrix<Treal>, "
4028 "Matrix<Treal>, Treal, "
4029 "Matrix<Treal>):Very strange error!!");
4037 template<
class Treal>
4039 symm(
const char side,
const char uplo,
const Treal alpha,
4051 C.resetRows(A.
rows);
4052 C.resetCols(B.
cols);
4056 C.resetRows(B.
rows);
4057 C.resetCols(A.
cols);
4061 if ((A.is_zero() || B.is_zero() || alpha == 0) &&
4062 (C.is_zero() || beta == 0))
4065 Treal beta_tmp = beta;
4070 if (!A.is_zero() && !B.is_zero() && alpha != 0) {
4079 throw Failure(
"Matrix<Treal>::symm: Incorrect matrix "
4080 "dimensions for symmetric multiply");
4087 template<
class Treal>
4089 syrk(
const char uplo,
const bool tA,
const Treal alpha,
4097 C.resetRows(A.
cols);
4098 C.resetCols(A.
cols);
4101 C.resetRows(A.
rows);
4102 C.resetCols(A.
rows);
4107 if (alpha != 0 && !A.is_zero()) {
4108 Treal beta_tmp = beta;
4126 throw Failure(
"Matrix<Treal>::syrk: Incorrect matrix "
4127 "dimensions for symmetric rank-k update");
4130 template<
class Treal>
4132 sysq(
const char uplo,
const Treal alpha,
4139 C.resetRows(A.
rows);
4140 C.resetCols(A.
cols);
4149 template<
class Treal>
4160 C.resetRows(A.
rows);
4161 C.resetCols(B.
cols);
4172 template<
class Treal>
4180 ssmm(alpha, A, B, beta, C);
4185 template<
class Treal>
4193 if (alpha != 0 && !A.is_zero() && !B.is_zero())
4194 if (((side ==
'R' && B.
ncols() == A.
nrows()) ||
4204 throw Failure(
"Matrix<Treal>::trmm: "
4205 "Incorrect matrix dimensions for multiplication");
4210 template<
class Treal>
4212 assert(!this->is_empty());
4213 if (this->is_zero())
4217 for (
int i = 0; i < this->nElements(); i++)
4223 template<
class Treal>
4226 assert(!this->is_empty());
4228 throw Failure(
"Matrix<Treal>::syFrobSquared(): "
4229 "Matrix must be have equal number of rows "
4230 "and cols to be symmetric");
4232 if (!this->is_zero()) {
4233 for (
int col = 1; col < this->
ncols(); col++)
4234 for (
int row = 0; row < col; row++)
4235 sum += 2 * (*
this)(row, col) * (*
this)(row, col);
4236 for (
int rc = 0; rc < this->
ncols(); rc++)
4237 sum += (*
this)(rc, rc) * (*
this)(rc, rc);
4242 template<
class Treal>
4247 assert(!A.is_empty());
4248 assert(!B.is_empty());
4250 throw Failure(
"Matrix<Treal>::frobSquaredDiff: "
4251 "Incorrect matrix dimensions");
4253 if (!A.is_zero() && !B.is_zero()) {
4255 for (
int i = 0; i < A.nElements(); i++) {
4260 else if (A.is_zero() && !B.is_zero())
4262 else if (!A.is_zero() && B.is_zero())
4269 template<
class Treal>
4274 assert(!A.is_empty());
4275 assert(!B.is_empty());
4279 throw Failure(
"Matrix<Treal>::syFrobSquaredDiff: "
4280 "Incorrect matrix dimensions");
4282 if (!A.is_zero() && !B.is_zero()) {
4284 for (
int col = 1; col < A.
ncols(); col++)
4285 for (
int row = 0; row < col; row++) {
4286 diff = A(row, col) - B(row, col);
4287 sum += 2 * diff * diff;
4289 for (
int rc = 0; rc < A.
ncols(); rc++) {
4290 diff = A(rc, rc) - B(rc, rc);
4294 else if (A.is_zero() && !B.is_zero())
4296 else if (!A.is_zero() && B.is_zero())
4302 template<
class Treal>
4305 assert(!this->is_empty());
4307 throw Failure(
"Matrix<Treal>::trace(): Matrix must be quadratic");
4308 if (this->is_zero())
4312 for (
int rc = 0; rc < this->
ncols(); rc++)
4313 sum += (*
this)(rc,rc);
4318 template<
class Treal>
4325 throw Failure(
"Matrix<Treal>::trace_ab: "
4326 "Wrong matrix dimensions for trace(A * B)");
4328 if (!A.is_zero() && !B.is_zero())
4329 for (
int rc = 0; rc < A.
nrows(); rc++)
4330 for (
int colA = 0; colA < A.
ncols(); colA++)
4331 tr += A(rc,colA) * B(colA, rc);
4335 template<
class Treal>
4342 throw Failure(
"Matrix<Treal>::trace_aTb: "
4343 "Wrong matrix dimensions for trace(A' * B)");
4345 if (!A.is_zero() && !B.is_zero())
4346 for (
int rc = 0; rc < A.
ncols(); rc++)
4347 for (
int rowB = 0; rowB < B.
nrows(); rowB++)
4348 tr += A(rowB,rc) * B(rowB, rc);
4352 template<
class Treal>
4360 throw Failure(
"Matrix<Treal>::sy_trace_ab: "
4361 "Wrong matrix dimensions for symmetric trace(A * B)");
4362 if (A.is_zero() || B.is_zero())
4367 for (
int rc = 0; rc < A.
nrows(); rc++)
4368 tr += A(rc,rc) * B(rc, rc);
4370 for (
int colA = 1; colA < A.
ncols(); colA++)
4371 for (
int rowA = 0; rowA < colA; rowA++)
4372 tr += 2 * A(rowA, colA) * B(rowA, colA);
4376 template<
class Treal>
4384 throw Failure(
"Matrix<Treal>::add: "
4385 "Wrong matrix dimensions for addition");
4386 if (!A.is_zero() && alpha != 0) {
4389 for (
int ind = 0; ind < A.nElements(); ind++)
4394 template<
class Treal>
4399 if (this->is_empty()) {
4400 this->resetRows(A.
rows);
4401 this->resetCols(A.
cols);
4409 template<
class Treal>
4412 if (!this->is_zero())
4416 template<
class Treal>
4419 if (!this->is_zero())
4420 for (
int ind = 0; ind < this->nElements(); ind++)
4426 template<
class Treal>
4431 if ((!this->is_zero() && this->frobSquared() <= threshold) ||
4432 (!ErrorMatrix->is_zero() &&
4436 else if (!this->is_zero() && this->frobSquared() <= threshold)
4440 template<
class Treal>
4444 assert(!this->is_empty());
4445 bool thisMatIsZero =
true;
4447 assert(!ErrorMatrix->is_empty());
4448 bool EMatIsZero =
true;
4449 if (!ErrorMatrix->is_zero() || !this->is_zero()) {
4450 if (ErrorMatrix->is_zero())
4452 if (this->is_zero())
4454 for (
int ind = 0; ind < this->nElements(); ind++) {
4455 if ( this->elements[ind] != 0 ) {
4456 assert(ErrorMatrix->
elements[ind] == 0);
4458 if ( this->elements[ind] * this->elements[ind] <= threshold ) {
4460 ErrorMatrix->
elements[ind] = this->elements[ind];
4461 this->elements[ind] = 0;
4465 thisMatIsZero =
false;
4468 if ( ErrorMatrix->
elements[ind] != 0 ) {
4469 assert(this->elements[ind] == 0);
4471 if ( ErrorMatrix->
elements[ind] * ErrorMatrix->
elements[ind] > threshold ) {
4473 this->elements[ind] = ErrorMatrix->
elements[ind];
4475 thisMatIsZero =
false;
4481 if (thisMatIsZero) {
4483 for (
int ind = 0; ind < this->nElements(); ind++)
4484 assert( this->elements[ind] == 0);
4490 for (
int ind = 0; ind < this->nElements(); ind++)
4491 assert( ErrorMatrix->
elements[ind] == 0);
4493 ErrorMatrix->
clear();
4498 if (!this->is_zero()) {
4499 for (
int ind = 0; ind < this->nElements(); ind++) {
4500 if ( this->elements[ind] * this->elements[ind] <= threshold )
4503 this->elements[ind] = 0;
4505 thisMatIsZero =
false;
4518 template<
class Treal>
4536 template<
class Treal>
4549 template<
class Treal>
4552 const bool tA,
const Treal alpha,
4566 template<
class Treal>
4587 template<
class Treal>
4590 assert(!this->is_empty());
4591 if (ErrorMatrix && ErrorMatrix->is_empty()) {
4592 ErrorMatrix->resetRows(this->rows);
4593 ErrorMatrix->resetCols(this->cols);
4595 Treal fs = frobSquared();
4596 if (fs < threshold) {
4606 template<
class Treal>
4610 const Treal threshold,
const side looking,
4614 throw Failure(
"Matrix<Treal>::inch: Matrix must be quadratic!");
4616 throw Failure(
"Matrix<Treal>::inch: Matrix is zero! "
4617 "Not possible to compute inverse cholesky.");
4622 throw Failure(
"Matrix<Treal>::inch: potrf returned info > 0. The matrix is not positive definite.");
4624 throw Failure(
"Matrix<Treal>::inch: potrf returned info < 0");
4628 throw Failure(
"Matrix<Treal>::inch: trtri returned info > 0. The matrix is not positive definite.");
4630 throw Failure(
"Matrix<Treal>::inch: trtri returned info < 0");
4635 template<
class Treal>
4638 assert(!this->is_empty());
4641 Treal* colsums =
new Treal[size];
4642 Treal* diag =
new Treal[size];
4643 for (
int ind = 0; ind < size; ind++)
4649 lmin = diag[0] - tmp1;
4650 lmax = diag[0] + tmp1;
4651 for (
int col = 1; col < size; col++) {
4653 tmp2 = diag[col] - tmp1;
4654 lmin = (tmp2 < lmin ? tmp2 : lmin);
4655 tmp2 = diag[col] + tmp1;
4656 lmax = (tmp2 > lmax ? tmp2 : lmax);
4662 throw Failure(
"Matrix<Treal>::gersgorin(Treal, Treal): Matrix must be quadratic");
4666 template<
class Treal>
4670 if (!this->is_zero())
4671 for (
int col = 0; col < this->
ncols(); col++)
4672 for (
int row = 0; row < this->
nrows(); row++)
4676 template<
class Treal>
4681 if (this->is_zero())
4685 for (
int rc = 0; rc < this->
ncols(); rc++)
4686 diag[rc] = (*
this)(rc,rc);
4713 template<
class Treal>
4715 if (this->is_zero() && this->cap > 0) {
4721 else if (this->cap > this->nel) {
4722 Treal* tmp =
new Treal[this->nel];
4723 for (
int i = 0; i < this->nel; i++)
4726 this->cap = this->nel;
4729 assert(this->cap == this->nel);
4736 template<
class Treal>
4737 void Matrix<Treal>::
4738 write_to_buffer_count(
int& zb_length,
int& vb_length)
const {
4739 if (this->is_zero()) {
4744 vb_length += this->nel;
4748 template<
class Treal>
4749 void Matrix<Treal>::
4750 write_to_buffer(
int* zerobuf,
const int zb_length,
4751 Treal* valuebuf,
const int vb_length,
4752 int& zb_index,
int& vb_index)
const {
4753 if (zb_index < zb_length) {
4754 if (this->is_zero()) {
4755 zerobuf[zb_index] = 0;
4759 if (vb_index + this->nel < vb_length + 1) {
4760 zerobuf[zb_index] = 1;
4762 for (
int i = 0; i < this->nel; i++)
4763 valuebuf[vb_index + i] = this->
elements[i];
4764 vb_index += this->nel;
4767 throw Failure(
"Matrix<Treal, Telement>::write_to_buffer: "
4768 "Insufficient space in buffers");
4772 throw Failure(
"Matrix<Treal, Telement>::write_to_buffer: "
4773 "Insufficient space in buffers");
4776 template<
class Treal>
4777 void Matrix<Treal>::
4778 read_from_buffer(
int* zerobuf,
const int zb_length,
4779 Treal* valuebuf,
const int vb_length,
4780 int& zb_index,
int& vb_index) {
4781 if (zb_index < zb_length) {
4782 if (zerobuf[zb_index] == 0) {
4787 this->content =
ful;
4789 this->assert_alloc();
4790 if (vb_index + this->nel < vb_length + 1) {
4791 assert(zerobuf[zb_index] == 1);
4793 for (
int i = 0; i < this->nel; i++)
4794 this->
elements[i] = valuebuf[vb_index + i];
4795 vb_index += this->nel;
4798 throw Failure(
"Matrix<Treal, Telement>::read_from_buffer: "
4799 "Mismatch, buffers does not match matrix");
4803 throw Failure(
"Matrix<Treal, Telement>::read_from_buffer: "
4804 "Mismatch, buffers does not match matrix");