25 #include <queso/GslMatrix.h>
26 #include <queso/GslVector.h>
27 #include <queso/Defines.h>
28 #include <gsl/gsl_linalg.h>
29 #include <gsl/gsl_eigen.h>
41 m_mat (gsl_matrix_calloc(map.NumGlobalElements(),nCols)),
49 m_determinant (-INFINITY),
50 m_lnDeterminant(-INFINITY),
65 m_mat (gsl_matrix_calloc(map.NumGlobalElements(),map.NumGlobalElements())),
73 m_determinant (-INFINITY),
74 m_lnDeterminant(-INFINITY),
81 for (
unsigned int i = 0; i <
m_mat->size1; ++i) {
82 (*this)(i,i) = diagValue;
91 m_mat (gsl_matrix_calloc(v.sizeLocal(),v.sizeLocal())),
99 m_determinant (-INFINITY),
100 m_lnDeterminant(-INFINITY),
101 m_permutation (NULL),
107 for (
unsigned int i = 0; i <
m_mat->size1; ++i) {
108 (*this)(i,i) = diagValue;
115 m_mat (gsl_matrix_calloc(v.sizeLocal(),v.sizeLocal())),
123 m_determinant (-INFINITY),
124 m_lnDeterminant(-INFINITY),
125 m_permutation (NULL),
132 for (
unsigned int i = 0; i <
dim; ++i) {
141 m_mat (gsl_matrix_calloc(B.numRowsLocal(),B.numCols())),
149 m_determinant (-INFINITY),
150 m_lnDeterminant(-INFINITY),
151 m_permutation (NULL),
181 iRC = gsl_matrix_scale(
m_mat,a);
225 iRC = gsl_matrix_memcpy(this->
m_mat, src.
m_mat);
235 gsl_matrix_free(
m_LU);
298 unsigned int nCols = this->
numCols();
300 for (
unsigned int i = 0; i < nRows; i++) {
301 for (
unsigned int j = 0; j < nCols; j++) {
316 unsigned int nCols = this->
numCols();
318 for (
unsigned int i = 0; i < nRows; i++) {
319 for (
unsigned int j = 0; j < nCols; j++) {
320 aux = fabs((*
this)(i,j));
321 if (aux > value) value = aux;
331 double value = -INFINITY;
334 unsigned int nCols = this->
numCols();
336 for (
unsigned int i = 0; i < nRows; i++) {
337 for (
unsigned int j = 0; j < nCols; j++) {
339 if (aux > value) value = aux;
350 unsigned int nCols = this->
numCols();
352 for (
unsigned int row = 0; row < nRows; ++row) {
353 for (
unsigned int col = 0; col < nCols; ++col) {
354 *gsl_matrix_ptr(
m_mat,row,col) = value;
363 unsigned int initialTargetRowId,
364 unsigned int initialTargetColId,
376 for (
unsigned int j = 0; j < mat.
numCols(); ++j) {
377 (*this)(initialTargetRowId+i,initialTargetColId+j) = mat(i,j);
386 unsigned int initialTargetRowId,
387 unsigned int initialTargetColId,
399 for (
unsigned int j = 0; j < mat.
numCols(); ++j) {
400 mat(i,j) = (*this)(initialTargetRowId+i,initialTargetColId+j) ;
412 gsl_error_handler_t* oldHandler;
413 oldHandler = gsl_set_error_handler_off();
414 iRC = gsl_linalg_cholesky_decomp(
m_mat);
416 std::cerr <<
"In GslMatrix::chol()"
418 <<
", gsl error message = " << gsl_strerror(iRC)
420 std::cerr <<
"Here is the offending matrix: " << std::endl;
421 std::cerr << *
this << std::endl;
423 gsl_set_error_handler(oldHandler);
428 "matrix is not positive definite",
438 unsigned int nCols = this->
numCols();
459 unsigned int nCols = this->
numCols();
470 <<
", this->numCols() = " << this->
numCols()
476 <<
"\n rhsVec.sizeLocal() = " << rhsVec.
sizeLocal()
477 <<
"\n solVec.sizeLocal() = " << solVec.
sizeLocal()
490 unsigned int nCols = this->
numCols();
501 for (
unsigned int j = 0; j < rhsMat.
numCols(); ++j) {
503 iRC = this->
svdSolve(rhsVec, solVec);
538 unsigned int nCols = this->
numCols();
552 struct timeval timevalBegin;
553 gettimeofday(&timevalBegin, NULL);
554 gsl_error_handler_t* oldHandler;
555 oldHandler = gsl_set_error_handler_off();
563 std::cerr <<
"In GslMatrix::internalSvd()"
565 <<
", gsl error message = " << gsl_strerror(iRC)
568 gsl_set_error_handler(oldHandler);
570 struct timeval timevalNow;
571 gettimeofday(&timevalNow, NULL);
579 "GslMatrix::internalSvd()",
594 unsigned int nCols = this->
numCols();
600 if (includeDiagonal) {
601 for (
unsigned int i = 0; i < nRows; i++) {
602 for (
unsigned int j = 0; j <= i; j++) {
608 for (
unsigned int i = 0; i < nRows; i++) {
609 for (
unsigned int j = 0; j < i; j++) {
622 unsigned int nCols = this->
numCols();
628 if (includeDiagonal) {
629 for (
unsigned int i = 0; i < nRows; i++) {
630 for (
unsigned int j = i; j < nCols; j++) {
636 for (
unsigned int i = 0; i < nRows; i++) {
637 for (
unsigned int j = (i+1); j < nCols; j++) {
650 unsigned int nCols = this->
numCols();
651 for (
unsigned int i = 0; i < nRows; ++i) {
652 for (
unsigned int j = 0; j < nCols; ++j) {
653 double aux = (*this)(i,j);
656 (-thresholdValue < aux)) {
660 (thresholdValue > aux)) {
673 unsigned int nCols = this->
numCols();
674 for (
unsigned int i = 0; i < nRows; ++i) {
675 for (
unsigned int j = 0; j < nCols; ++j) {
676 double aux = (*this)(i,j);
679 (-thresholdValue > aux)) {
683 (thresholdValue < aux)) {
696 unsigned int nCols = this->
numCols();
701 for (
unsigned int row = 0; row < nRows; ++row) {
702 for (
unsigned int col = 0; col < nCols; ++col) {
703 mat(row,col) = (*this)(col,row);
714 unsigned int nCols = this->
numCols();
721 unitVector.
cwSet(0.);
723 for (
unsigned int j = 0; j < nCols; ++j) {
724 if (j > 0) unitVector[j-1] = 0.;
727 for (
unsigned int i = 0; i < nRows; ++i) {
728 (*m_inverse)(i,j) = multVector[i];
740 <<
"\n M = " << *
this
742 <<
"\n M*M^{-1} = " << (*this)*(*m_inverse)
743 <<
"\n M^{-1}*M = " << (*m_inverse)*(*this)
752 unsigned int initialTargetRowId,
753 unsigned int initialTargetColId,
754 const std::vector<const GslMatrix* >& matrices,
755 bool checkForExactNumRowsMatching,
756 bool checkForExactNumColsMatching)
758 unsigned int sumNumRowsLocals = 0;
759 unsigned int sumNumCols = 0;
760 for (
unsigned int i = 0; i < matrices.size(); ++i) {
761 sumNumRowsLocals += matrices[i]->numRowsLocal();
762 sumNumCols += matrices[i]->numCols();
769 unsigned int cumulativeRowId = 0;
770 unsigned int cumulativeColId = 0;
771 for (
unsigned int i = 0; i < matrices.size(); ++i) {
772 unsigned int nRows = matrices[i]->numRowsLocal();
773 unsigned int nCols = matrices[i]->numCols();
774 for (
unsigned int rowId = 0; rowId < nRows; ++rowId) {
775 for (
unsigned int colId = 0; colId < nCols; ++colId) {
776 (*this)(initialTargetRowId + cumulativeRowId + rowId, initialTargetColId + cumulativeColId + colId) = (*(matrices[i]))(rowId,colId);
779 cumulativeRowId += nRows;
780 cumulativeColId += nCols;
788 unsigned int initialTargetRowId,
789 unsigned int initialTargetColId,
790 const std::vector<GslMatrix* >& matrices,
791 bool checkForExactNumRowsMatching,
792 bool checkForExactNumColsMatching)
794 unsigned int sumNumRowsLocals = 0;
795 unsigned int sumNumCols = 0;
796 for (
unsigned int i = 0; i < matrices.size(); ++i) {
797 sumNumRowsLocals += matrices[i]->numRowsLocal();
798 sumNumCols += matrices[i]->numCols();
805 unsigned int cumulativeRowId = 0;
806 unsigned int cumulativeColId = 0;
807 for (
unsigned int i = 0; i < matrices.size(); ++i) {
808 unsigned int nRows = matrices[i]->numRowsLocal();
809 unsigned int nCols = matrices[i]->numCols();
810 for (
unsigned int rowId = 0; rowId < nRows; ++rowId) {
811 for (
unsigned int colId = 0; colId < nCols; ++colId) {
812 (*this)(initialTargetRowId + cumulativeRowId + rowId, initialTargetColId + cumulativeColId + colId) = (*(matrices[i]))(rowId,colId);
815 cumulativeRowId += nRows;
816 cumulativeColId += nCols;
824 unsigned int initialTargetRowId,
825 unsigned int initialTargetColId,
826 const std::vector<const GslMatrix* >& matrices,
827 bool checkForExactNumRowsMatching,
828 bool checkForExactNumColsMatching)
830 unsigned int sumNumCols = 0;
831 for (
unsigned int i = 0; i < matrices.size(); ++i) {
834 sumNumCols += matrices[i]->numCols();
839 unsigned int cumulativeColId = 0;
840 for (
unsigned int i = 0; i < matrices.size(); ++i) {
841 unsigned int nRows = matrices[i]->numRowsLocal();
842 unsigned int nCols = matrices[i]->numCols();
843 for (
unsigned int rowId = 0; rowId < nRows; ++rowId) {
844 for (
unsigned int colId = 0; colId < nCols; ++colId) {
845 (*this)(initialTargetRowId + rowId, initialTargetColId + cumulativeColId + colId) = (*(matrices[i]))(rowId,colId);
848 cumulativeColId += nCols;
856 unsigned int initialTargetRowId,
857 unsigned int initialTargetColId,
858 const std::vector<GslMatrix* >& matrices,
859 bool checkForExactNumRowsMatching,
860 bool checkForExactNumColsMatching)
862 unsigned int sumNumCols = 0;
863 for (
unsigned int i = 0; i < matrices.size(); ++i) {
866 sumNumCols += matrices[i]->numCols();
871 unsigned int cumulativeColId = 0;
872 for (
unsigned int i = 0; i < matrices.size(); ++i) {
873 unsigned int nRows = matrices[i]->numRowsLocal();
874 unsigned int nCols = matrices[i]->numCols();
875 for (
unsigned int rowId = 0; rowId < nRows; ++rowId) {
876 for (
unsigned int colId = 0; colId < nCols; ++colId) {
877 (*this)(initialTargetRowId + rowId, initialTargetColId + cumulativeColId + colId) = (*(matrices[i]))(rowId,colId);
880 cumulativeColId += nCols;
888 unsigned int initialTargetRowId,
889 unsigned int initialTargetColId,
890 const std::vector<const GslMatrix* >& matrices,
891 bool checkForExactNumRowsMatching,
892 bool checkForExactNumColsMatching)
894 unsigned int sumNumRows = 0;
895 for (
unsigned int i = 0; i < matrices.size(); ++i) {
898 sumNumRows += matrices[i]->numRowsLocal();
903 unsigned int cumulativeRowId = 0;
904 for (
unsigned int i = 0; i < matrices.size(); ++i) {
905 unsigned int nRows = matrices[i]->numRowsLocal();
906 unsigned int nCols = matrices[i]->numCols();
907 for (
unsigned int rowId = 0; rowId < nRows; ++rowId) {
908 for (
unsigned int colId = 0; colId < nCols; ++colId) {
909 (*this)(initialTargetRowId + cumulativeRowId + rowId, initialTargetColId + colId) = (*(matrices[i]))(rowId,colId);
912 cumulativeRowId += nRows;
920 unsigned int initialTargetRowId,
921 unsigned int initialTargetColId,
922 const std::vector<GslMatrix* >& matrices,
923 bool checkForExactNumRowsMatching,
924 bool checkForExactNumColsMatching)
926 unsigned int sumNumRows = 0;
927 for (
unsigned int i = 0; i < matrices.size(); ++i) {
930 sumNumRows += matrices[i]->numRowsLocal();
935 unsigned int cumulativeRowId = 0;
936 for (
unsigned int i = 0; i < matrices.size(); ++i) {
937 unsigned int nRows = matrices[i]->numRowsLocal();
938 unsigned int nCols = matrices[i]->numCols();
939 for (
unsigned int rowId = 0; rowId < nRows; ++rowId) {
940 for (
unsigned int colId = 0; colId < nCols; ++colId) {
941 (*this)(initialTargetRowId + cumulativeRowId + rowId, initialTargetColId + colId) = (*(matrices[i]))(rowId,colId);
944 cumulativeRowId += nRows;
952 unsigned int initialTargetRowId,
953 unsigned int initialTargetColId,
956 bool checkForExactNumRowsMatching,
957 bool checkForExactNumColsMatching)
964 for (
unsigned int rowId1 = 0; rowId1 < mat1.
numRowsLocal(); ++rowId1) {
965 for (
unsigned int colId1 = 0; colId1 < mat1.
numCols(); ++colId1) {
966 double multiplicativeFactor = mat1(rowId1,colId1);
967 unsigned int targetRowId = rowId1 * mat2.
numRowsLocal();
968 unsigned int targetColId = colId1 * mat2.
numCols();
969 for (
unsigned int rowId2 = 0; rowId2 < mat2.
numRowsLocal(); ++rowId2) {
970 for (
unsigned int colId2 = 0; colId2 < mat2.
numCols(); ++colId2) {
971 (*this)(initialTargetRowId + targetRowId + rowId2, initialTargetColId + targetColId + colId2) = multiplicativeFactor * mat2(rowId2,colId2);
982 unsigned int initialTargetRowId,
983 unsigned int initialTargetColId,
986 bool checkForExactNumRowsMatching,
987 bool checkForExactNumColsMatching)
994 for (
unsigned int rowId1 = 0; rowId1 < mat1.
numRowsLocal(); ++rowId1) {
995 for (
unsigned int colId1 = 0; colId1 < mat1.
numCols(); ++colId1) {
996 double multiplicativeFactor = mat1(rowId1,colId1);
997 unsigned int targetRowId = rowId1 * vec2.
sizeLocal();
998 unsigned int targetColId = colId1 * 1;
999 for (
unsigned int rowId2 = 0; rowId2 < vec2.
sizeLocal(); ++rowId2) {
1000 for (
unsigned int colId2 = 0; colId2 < 1; ++colId2) {
1001 (*this)(initialTargetRowId + targetRowId + rowId2, initialTargetColId + targetColId + colId2) = multiplicativeFactor * vec2[rowId2];
1013 unsigned int initialTargetRowId,
1014 unsigned int initialTargetColId,
1016 bool checkForExactNumRowsMatching,
1017 bool checkForExactNumColsMatching)
1020 unsigned int nCols = mat.
numCols();
1026 for (
unsigned int row = 0; row < nRows; ++row) {
1027 for (
unsigned int col = 0; col < nCols; ++col) {
1028 (*this)(initialTargetRowId + col, initialTargetColId + row) = mat(row,col);
1044 <<
": before 'this->invertMultiply()'"
1050 <<
": after 'this->invertMultiply()'"
1056 <<
": before 'gsl_linalg_LU_det()'"
1062 <<
": after 'gsl_linalg_LU_det()'"
1068 <<
": after 'gsl_linalg_LU_lndet()'"
1085 <<
": before 'this->invertMultiply()'"
1091 <<
": after 'this->invertMultiply()'"
1097 <<
": before 'gsl_linalg_LU_det()'"
1103 <<
": after 'gsl_linalg_LU_det()'"
1109 <<
": before 'gsl_linalg_LU_lndet()'"
1125 if (relativeVec[0] > 0.) {
1126 relativeVec = (1./relativeVec[0])*relativeVec;
1129 unsigned int rankValue = 0;
1130 for (
unsigned int i = 0; i < relativeVec.
sizeLocal(); ++i) {
1131 if (( (*
m_svdSvec)[i] >= absoluteZeroThreshold ) &&
1132 ( relativeVec [i] >= relativeZeroThreshold )) {
1140 <<
", this->numCols() = " << this->
numCols()
1141 <<
", absoluteZeroThreshold = " << absoluteZeroThreshold
1142 <<
", relativeZeroThreshold = " << relativeZeroThreshold
1143 <<
", rankValue = " << rankValue
1145 <<
", relativeVec = " << relativeVec
1173 unsigned int sizeX = this->
numCols();
1175 for (
unsigned int i = 0; i < sizeY; ++i) {
1177 for (
unsigned int j = 0; j < sizeX; ++j) {
1178 value += (*this)(i,j)*x[j];
1221 std::cout <<
"In GslMatrix::invertMultiply()"
1222 <<
": before LU decomposition, m_LU = ";
1223 gsl_matrix_fprintf(stdout,
m_LU,
"%f");
1224 std::cout << std::endl;
1227 gsl_error_handler_t* oldHandler;
1228 oldHandler = gsl_set_error_handler_off();
1231 <<
": before 'gsl_linalg_LU_decomp()'"
1236 std::cerr <<
"In GslMatrix::invertMultiply()"
1237 <<
", after gsl_linalg_LU_decomp()"
1238 <<
": iRC = " << iRC
1239 <<
", gsl error message = " << gsl_strerror(iRC)
1242 gsl_set_error_handler(oldHandler);
1245 <<
": after 'gsl_linalg_LU_decomp()'"
1246 <<
", IRC = " << iRC
1252 std::cout <<
"In GslMatrix::invertMultiply()"
1253 <<
": after LU decomposition, m_LU = ";
1254 gsl_matrix_fprintf(stdout,
m_LU,
"%f");
1255 std::cout << std::endl;
1259 gsl_error_handler_t* oldHandler;
1260 oldHandler = gsl_set_error_handler_off();
1263 <<
": before 'gsl_linalg_LU_solve()'"
1269 std::cerr <<
"In GslMatrix::invertMultiply()"
1270 <<
", after gsl_linalg_LU_solve()"
1271 <<
": iRC = " << iRC
1272 <<
", gsl error message = " << gsl_strerror(iRC)
1275 gsl_set_error_handler(oldHandler);
1278 <<
": after 'gsl_linalg_LU_solve()'"
1279 <<
", IRC = " << iRC
1290 std::cout <<
"In GslMatrix::invertMultiply()"
1291 <<
": ||b - Ax||_2 = " << tmpVec.
norm2()
1292 <<
": ||b - Ax||_2/||b||_2 = " << tmpVec.
norm2()/b.
norm2()
1313 "Matrices B and X are incompatible");
1315 "Matrices B and X are incompatible");
1317 "This and X matrices are incompatible");
1323 for(
unsigned int j = 0; j < B.
numCols(); ++j )
1360 if (
m_LU == NULL ) {
1387 unsigned int n = eigenValues.
sizeLocal();
1395 if (eigenVectors == NULL) {
1396 gsl_eigen_symm_workspace* w = gsl_eigen_symm_alloc((
size_t) n);
1397 gsl_eigen_symm(
m_mat,eigenValues.
data(),w);
1398 gsl_eigen_symm_free(w);
1401 gsl_eigen_symmv_workspace* w = gsl_eigen_symmv_alloc((
size_t) n);
1403 gsl_eigen_symmv_sort(eigenValues.
data(),eigenVectors->
m_mat,GSL_EIGEN_SORT_VAL_ASC);
1404 gsl_eigen_symmv_free(w);
1415 unsigned int n = eigenVector.
sizeLocal();
1426 const unsigned int max_num_iterations = 10000;
1427 const double tolerance = 1.0e-13;
1441 for(
unsigned int k = 0;
k < max_num_iterations; ++
k )
1448 index = (w.
abs()).getMaxValueIndex();
1452 z = (1.0/lambda) * w;
1456 residual = ( (*this)*z - lambda*z ).norm2();
1458 if( residual < tolerance )
1460 eigenValue = lambda;
1483 unsigned int n = eigenVector.
sizeLocal();
1494 const unsigned int max_num_iterations = 1000;
1495 const double tolerance = 1.0e-13;
1507 double one_over_lambda;
1510 for(
unsigned int k = 0;
k < max_num_iterations; ++
k )
1512 w = (*this).invertMultiplyForceLU( z );
1517 index = (w.
abs()).getMaxValueIndex();
1520 one_over_lambda = w[index];
1522 lambda = 1.0/one_over_lambda;
1528 residual = ( (*this)*z - lambda*z ).norm2();
1530 if( residual < tolerance )
1532 eigenValue = lambda;
1559 gsl_vector* gsl_column = gsl_vector_alloc( column.
sizeLocal() );
1561 int error_code = gsl_matrix_get_col( gsl_column,
m_mat, column_num );
1565 for(
unsigned int i = 0; i < column.
sizeLocal(); ++i )
1567 column[i] = gsl_vector_get( gsl_column, i );
1571 gsl_vector_free( gsl_column );
1585 gsl_vector* gsl_row = gsl_vector_alloc( row.
sizeLocal() );
1587 int error_code = gsl_matrix_get_row( gsl_row,
m_mat, row_num );
1591 for(
unsigned int i = 0; i < row.
sizeLocal(); ++i )
1593 row[i] = gsl_vector_get( gsl_row, i );
1597 gsl_vector_free( gsl_row );
1610 this->
getRow( row_num, row );
1637 gsl_vector* gsl_row = row.
data();
1639 int error_code = gsl_matrix_set_row(
m_mat, row_num, gsl_row );
1654 gsl_vector* gsl_column = column.
data();
1656 int error_code = gsl_matrix_set_col(
m_mat, column_num, gsl_column );
1669 "GslMatrix::mpiSum()",
1670 "local and global matrices incompatible");
1674 std::vector<double> local( size, 0.0 );
1675 std::vector<double> global( size, 0.0 );
1678 for(
unsigned int i = 0; i < this->
numRowsLocal(); i++ )
1680 for(
unsigned int j = 0; j < this->
numCols(); j++ )
1684 local[
k] = (*this)(i,j);
1689 "GslMatrix::mpiSum()",
1690 "failed MPI.Allreduce()");
1692 for(
unsigned int i = 0; i < this->
numRowsLocal(); i++ )
1694 for(
unsigned int j = 0; j < this->
numCols(); j++ )
1698 M_global(i,j) = global[
k];
1721 for (
unsigned int colId = 0; colId < y1Mat.
numCols(); ++colId) {
1740 unsigned int nCols = this->
numCols();
1743 for (
unsigned int i = 0; i < nRows; ++i) {
1744 for (
unsigned int j = 0; j < nCols; ++j) {
1748 if (i != (nRows-1)) os <<
"; ";
1753 for (
unsigned int i = 0; i < nRows; ++i) {
1754 for (
unsigned int j = 0; j < nCols; ++j) {
1767 const std::string& varNamePrefix,
1768 const std::string& fileName,
1769 const std::string& fileType,
1770 const std::set<unsigned int>& allowedSubEnvIds)
const
1783 unsigned int nCols = this->
numCols();
1790 for (
unsigned int i = 0; i < nRows; ++i) {
1791 for (
unsigned int j = 0; j < nCols; ++j) {
1792 *filePtrSet.
ofsVar << (*this)(i,j)
1795 *filePtrSet.
ofsVar <<
"\n";
1797 *filePtrSet.
ofsVar <<
"];\n";
1807 const std::string& fileName,
1808 const std::string& fileType,
1809 const std::set<unsigned int>& allowedSubEnvIds)
1825 unsigned int idOfMyFirstLine = 1;
1826 unsigned int idOfMyLastLine = nRowsLocal;
1827 unsigned int nCols = this->
numCols();
1831 std::string tmpString;
1834 *filePtrSet.
ifsVar >> tmpString;
1838 *filePtrSet.
ifsVar >> tmpString;
1843 *filePtrSet.
ifsVar >> tmpString;
1845 unsigned int posInTmpString = 6;
1848 char nRowsString[tmpString.size()-posInTmpString+1];
1849 unsigned int posInRowsString = 0;
1852 nRowsString[posInRowsString++] = tmpString[posInTmpString++];
1853 }
while (tmpString[posInTmpString] !=
',');
1854 nRowsString[posInRowsString] =
'\0';
1858 char nColsString[tmpString.size()-posInTmpString+1];
1859 unsigned int posInColsString = 0;
1862 nColsString[posInColsString++] = tmpString[posInTmpString++];
1863 }
while (tmpString[posInTmpString] !=
')');
1864 nColsString[posInColsString] =
'\0';
1867 unsigned int numRowsInFile = (
unsigned int) strtod(nRowsString,NULL);
1868 unsigned int numColsInFile = (
unsigned int) strtod(nColsString,NULL);
1872 <<
", numRowsInFile = " << numRowsInFile
1873 <<
", numColsInFile = " << numColsInFile
1874 <<
", nRowsLocal = " << nRowsLocal
1875 <<
", nCols = " << nCols
1883 queso_require_equal_to_msg(numColsInFile, nCols,
"number of parameters of vec in file is different than number of parameters in this vec object");
1886 unsigned int maxCharsPerLine = 64*nCols;
1888 unsigned int lineId = 0;
1889 while (lineId < idOfMyFirstLine) {
1890 filePtrSet.
ifsVar->ignore(maxCharsPerLine,
'\n');
1896 <<
": beginning to read input actual data"
1903 *filePtrSet.
ifsVar >> tmpString;
1907 *filePtrSet.
ifsVar >> tmpString;
1912 std::streampos tmpPos = filePtrSet.
ifsVar->tellg();
1913 filePtrSet.
ifsVar->seekg(tmpPos+(std::streampos)2);
1917 <<
": beginning to read lines with numbers only"
1918 <<
", lineId = " << lineId
1919 <<
", idOfMyFirstLine = " << idOfMyFirstLine
1920 <<
", idOfMyLastLine = " << idOfMyLastLine
1925 while (lineId <= idOfMyLastLine) {
1926 for (
unsigned int j = 0; j < nCols; ++j) {
1927 *filePtrSet.
ifsVar >> tmpRead;
1928 (*this)(lineId-idOfMyFirstLine,j) = tmpRead;
1962 unsigned int m1Cols = m1.
numCols();
1964 unsigned int m2Cols = m2.
numCols();
1972 unsigned int commonSize = m1Cols;
1973 for (
unsigned int row1 = 0; row1 < m1Rows; ++row1) {
1974 for (
unsigned int col2 = 0; col2 < m2Cols; ++col2) {
1976 for (
unsigned int k = 0;
k < commonSize; ++
k) {
1977 result += m1(row1,
k)*m2(
k,col2);
1979 mat(row1,col2) = result;
2007 for (
unsigned int i = 0; i < nRows; ++i) {
2008 double value1 = v1[i];
2009 for (
unsigned int j = 0; j < nCols; ++j) {
2010 answer(i,j) = value1*v2[j];
2021 unsigned int mCols = mat.
numCols();
2028 for (
unsigned int i = 0; i < mRows; ++i) {
2029 double vecValue = vec[i];
2030 for (
unsigned int j = 0; j < mCols; ++j) {
2031 answer(i,j) *= vecValue;
2042 unsigned int mCols = mat.
numCols();
2049 for (
unsigned int j = 0; j < mCols; ++j) {
2050 double vecValue = vec[j];
2051 for (
unsigned int i = 0; i < mRows; ++i) {
2052 answer(i,j) *= vecValue;
GslMatrix transpose() const
This function calculated the transpose of this matrix (square).
unsigned int displayVerbosity() const
void cwSet(double value)
Component-wise set all values to this with value.
int chol()
Computes Cholesky factorization of a real symmetric positive definite matrix this.
bool m_printHorizontally
Flag for either or not print this matrix.
const GslMatrix & svdMatU() const
This function calls private member GslMatrix::internalSvd() to set a M-by-N orthogonal matrix U of th...
double normFrob() const
Returns the Frobenius norm of this matrix.
void mpiSum(const MpiComm &comm, GslMatrix &M_global) const
void setRow(const unsigned int row_num, const GslVector &row)
This function copies vector row into the row_num-th column of this matrix.
std::ofstream * ofsVar
Provides a stream interface to write data to files.
GslMatrix & operator-=(const GslMatrix &rhs)
Stores in this the coordinate-wise subtraction of this by rhs.
int subRank() const
Access function for sub-rank.
GslMatrix & operator*=(double a)
Stores in this the coordinate-wise multiplication of this and a.
GslMatrix()
Default Constructor.
void setColumn(const unsigned int column_num, const GslVector &column)
This function copies vector column into the column_num-th column of this matrix.
const int UQ_MATRIX_IS_NOT_POS_DEFINITE_RC
Map * m_svdColMap
Mapping for matrices involved in the singular value decomposition (svd) routine.
#define UQ_FATAL_RC_MACRO(macroIRc, givenRank, where, what)
bool m_inDebugMode
Flag for either or not QUESO is in debug mode.
void zeroUpper(bool includeDiagonal=false)
This function sets all the entries above the main diagonal of this matrix to zero.
double m_determinant
The determinant of this matrix.
const Map & m_map
Mapping variable.
Class for matrix operations (virtual).
void subWriteContents(const std::string &varNamePrefix, const std::string &fileName, const std::string &fileType, const std::set< unsigned int > &allowedSubEnvIds) const
Write contents of subenvironment in file fileName.
bool m_isSingular
Indicates whether or not this matrix is singular.
A class for partitioning vectors and matrices.
void cwExtract(unsigned int rowId, unsigned int colId, GslMatrix &mat) const
#define queso_require_less_msg(expr1, expr2, msg)
unsigned int numCols() const
Returns the column dimension of this matrix.
int fullRank() const
Returns the process full rank.
int m_signum
m_signum stores the sign of the permutation of the LU decomposition PA = LU.
void matlabLinearInterpExtrap(const GslVector &x1Vec, const GslVector &y1Vec, const GslVector &x2Vec)
void zeroLower(bool includeDiagonal=false)
This function sets all the entries bellow the main diagonal of this matrix to zero.
GslMatrix matrixProduct(const GslVector &v1, const GslVector &v2)
unsigned int numRowsLocal() const
Returns the local row dimension of this matrix.
#define queso_require_not_equal_to_msg(expr1, expr2, msg)
void fillWithBlocksVertically(unsigned int rowId, unsigned int colId, const std::vector< const GslMatrix * > &matrices, bool checkForExactNumRowsMatching, bool checkForExactNumColsMatching)
This function fills this matrix vertically with const block matrices.
GslMatrix operator*(double a, const GslMatrix &mat)
const int UQ_MATRIX_SVD_FAILED_RC
void print(std::ostream &os) const
Print method. Defines the behavior of the ostream << operator inherited from the Object class...
#define queso_require_less_equal_msg(expr1, expr2, msg)
double m_lnDeterminant
The natural logarithm of the determinant of this matrix.
int svd(GslMatrix &matU, GslVector &vecS, GslMatrix &matVt) const
Checks for the dimension of this matrix, matU, VecS and matVt, and calls the protected routine intern...
#define queso_require_msg(asserted, msg)
const std::string & subIdString() const
Access to the attribute m_subIdString; which stores the string for the sub-environment, and it will be used, for instance, to create the output files for each sub-environment.
#define queso_require_equal_to_msg(expr1, expr2, msg)
double lnDeterminant() const
Calculates the ln(determinant) of this matrix.
int internalSvd() const
This function factorizes the M-by-N matrix A into the singular value decomposition A = U S V^T for M ...
#define queso_require_greater_equal_msg(expr1, expr2, msg)
void cwSet(double value)
Component-wise sets all values to this with value.
void fillWithTensorProduct(unsigned int rowId, unsigned int colId, const GslMatrix &mat1, const GslMatrix &mat2, bool checkForExactNumRowsMatching, bool checkForExactNumColsMatching)
This function calculates the tensor product of matrices mat1 and mat2 and stores it in this matrix...
void closeFile(FilePtrSetStruct &filePtrSet, const std::string &fileType) const
Closes the file.
GslMatrix * m_svdVmat
m_svdVmat stores the N-by-N orthogonal square matrix V after the singular value decomposition of a ma...
void matlabLinearInterpExtrap(const GslVector &x1Vec, const GslMatrix &y1Mat, const GslVector &x2Vec)
GslVector * m_svdSvec
m_svdSvec stores the diagonal of the N-by-N diagonal matrix of singular values S after the singular v...
GslMatrix leftDiagScaling(const GslVector &vec, const GslMatrix &mat)
std::ostream & operator<<(std::ostream &os, const BaseEnvironment &obj)
GslVector abs() const
This function returns absolute value of elements in this.
void fillWithTranspose(unsigned int rowId, unsigned int colId, const GslMatrix &mat, bool checkForExactNumRowsMatching, bool checkForExactNumColsMatching)
This function stores the transpose of this matrix into this matrix.
GslMatrix rightDiagScaling(const GslMatrix &mat, const GslVector &vec)
gsl_matrix * m_mat
GSL matrix, also referred to as this matrix.
Class for matrix operations using GSL library.
Struct for handling data input and output from files.
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
gsl_matrix * m_LU
GSL matrix for the LU decomposition of m_mat.
const BaseEnvironment & env() const
double max() const
Returns the maximum element value of the matrix.
This (virtual) class sets up the environment underlying the use of the QUESO library by an executable...
void filterLargeValues(double thresholdValue)
This function sets to zero (filters) all entries of this matrix which are greater than thresholdValue...
gsl_permutation * m_permutation
The permutation matrix of a LU decomposition.
GslMatrix * m_inverse
Inverse matrix of this.
unsigned int numOfProcsForStorage() const
void smallestEigen(double &eigenValue, GslVector &eigenVector) const
This function finds smallest eigenvalue, namely eigenValue, of this matrix and its corresponding eige...
GslVector invertMultiplyForceLU(const GslVector &b) const
This function calculates the inverse of this matrix and multiplies it with vector b...
unsigned int numRowsGlobal() const
Returns the global row dimension of this matrix.
gsl_vector * data() const
double normMax() const
Returns the Frobenius norm of this matrix.
void fillWithBlocksDiagonally(unsigned int rowId, unsigned int colId, const std::vector< const GslMatrix * > &matrices, bool checkForExactNumRowsMatching, bool checkForExactNumColsMatching)
This function fills this matrix diagonally with const block matrices.
const BaseEnvironment & m_env
QUESO environment variable.
unsigned int sizeLocal() const
Returns the length of this vector.
std::ifstream * ifsVar
Provides a stream interface to read data from files.
gsl_matrix * data()
Returns this matrix.
GslMatrix operator-(const GslMatrix &m1, const GslMatrix &m2)
GslMatrix inverse() const
This function calculated the inverse of this matrix (square).
void subReadContents(const std::string &fileName, const std::string &fileType, const std::set< unsigned int > &allowedSubEnvIds)
Read contents of subenvironment from file fileName.
void filterSmallValues(double thresholdValue)
This function sets to zero (filters) all entries of this matrix which are smaller than thresholdValue...
GslMatrix & operator+=(const GslMatrix &rhs)
Stores in this the coordinate-wise addition of this and rhs.
void getRow(const unsigned int row_num, GslVector &row) const
This function gets the row_num-th column of this matrix and stores it into vector row...
GslVector multiply(const GslVector &x) const
This function multiplies this matrix by vector x and returns the resulting vector.
GslMatrix operator+(const GslMatrix &m1, const GslMatrix &m2)
void Allreduce(void *sendbuf, void *recvbuf, int count, RawType_MPI_Datatype datatype, RawType_MPI_Op op, const char *whereMsg, const char *whatMsg) const
Combines values from all processes and distributes the result back to all processes.
void fillWithBlocksHorizontally(unsigned int rowId, unsigned int colId, const std::vector< const GslMatrix * > &matrices, bool checkForExactNumRowsMatching, bool checkForExactNumColsMatching)
This function fills this matrix horizontally with const block matrices.
bool openOutputFile(const std::string &fileName, const std::string &fileType, const std::set< unsigned int > &allowedSubEnvIds, bool writeOver, FilePtrSetStruct &filePtrSet) const
Opens an output file for each sub-environment that was chosen to send data to the file...
GslMatrix * m_svdUmat
m_svdUmat stores the M-by-N orthogonal matrix U after the singular value decomposition of a matrix...
Class for vector operations using GSL library.
unsigned int checkingLevel() const
Access function to private attribute m_checkingLevel.
const GslMatrix & svdMatV() const
This function calls private member GslMatrix::internalSvd() to set a N-by-N orthogonal square matrix ...
The QUESO MPI Communicator Class.
void copy(const GslMatrix &src)
In this function this matrix receives a copy of matrix src.
GslMatrix & operator=(const GslMatrix &rhs)
Copies values from matrix rhs to this.
int worldRank() const
Returns the process world rank.
const BaseEnvironment & env() const
void resetLU()
In this function resets the LU decomposition of this matrix, as well as deletes the private member po...
unsigned int rank(double absoluteZeroThreshold, double relativeZeroThreshold) const
This function returns the number of singular values of this matrix (rank).
void getColumn(const unsigned int column_num, GslVector &column) const
This function gets the column_num-th column of this matrix and stores it into vector column...
#define queso_require_greater_msg(expr1, expr2, msg)
GslVector invertMultiply(const GslVector &b) const
This function calculates the inverse of this matrix and multiplies it with vector b...
GslMatrix & operator/=(double a)
Stores in this the coordinate-wise division of this by a.
GslMatrix * m_svdVTmat
m_svdVmatT stores the transpose of N-by-N orthogonal square matrix V, namely V^T, after the singular ...
virtual void base_copy(const Matrix &src)
Copies base data from matrix src to this matrix.
void largestEigen(double &eigenValue, GslVector &eigenVector) const
This function finds largest eigenvalue, namely eigenValue, of this matrix and its corresponding eigen...
void eigen(GslVector &eigenValues, GslMatrix *eigenVectors) const
This function computes the eigenvalues of a real symmetric matrix.
#define UQ_RC_MACRO(macroIRc, givenRank, where, what, retValue)
double norm2() const
Returns the 2-norm (Euclidean norm) of the vector.
double determinant() const
Calculates the determinant of this matrix.
int svdSolve(const GslVector &rhsVec, GslVector &solVec) const
This function solves the system A x = b using the singular value decomposition (U, S, V) of A which must have been computed previously with GslMatrix::svd (x=solVec, b=rhsVec).
bool openInputFile(const std::string &fileName, const std::string &fileType, const std::set< unsigned int > &allowedSubEnvIds, FilePtrSetStruct &filePtrSet) const
Opens an input file.