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();
699 Map serial_map(nCols, 0, comm);
703 for (
unsigned int row = 0; row < nRows; ++row) {
704 for (
unsigned int col = 0; col < nCols; ++col) {
705 mat(col,row) = (*this)(row,col);
716 unsigned int nCols = this->
numCols();
723 unitVector.
cwSet(0.);
725 for (
unsigned int j = 0; j < nCols; ++j) {
726 if (j > 0) unitVector[j-1] = 0.;
729 for (
unsigned int i = 0; i < nRows; ++i) {
730 (*m_inverse)(i,j) = multVector[i];
742 <<
"\n M = " << *
this
744 <<
"\n M*M^{-1} = " << (*this)*(*m_inverse)
745 <<
"\n M^{-1}*M = " << (*m_inverse)*(*this)
754 unsigned int initialTargetRowId,
755 unsigned int initialTargetColId,
756 const std::vector<const GslMatrix* >& matrices,
757 bool checkForExactNumRowsMatching,
758 bool checkForExactNumColsMatching)
760 unsigned int sumNumRowsLocals = 0;
761 unsigned int sumNumCols = 0;
762 for (
unsigned int i = 0; i < matrices.size(); ++i) {
763 sumNumRowsLocals += matrices[i]->numRowsLocal();
764 sumNumCols += matrices[i]->numCols();
771 unsigned int cumulativeRowId = 0;
772 unsigned int cumulativeColId = 0;
773 for (
unsigned int i = 0; i < matrices.size(); ++i) {
774 unsigned int nRows = matrices[i]->numRowsLocal();
775 unsigned int nCols = matrices[i]->numCols();
776 for (
unsigned int rowId = 0; rowId < nRows; ++rowId) {
777 for (
unsigned int colId = 0; colId < nCols; ++colId) {
778 (*this)(initialTargetRowId + cumulativeRowId + rowId, initialTargetColId + cumulativeColId + colId) = (*(matrices[i]))(rowId,colId);
781 cumulativeRowId += nRows;
782 cumulativeColId += nCols;
790 unsigned int initialTargetRowId,
791 unsigned int initialTargetColId,
792 const std::vector<GslMatrix* >& matrices,
793 bool checkForExactNumRowsMatching,
794 bool checkForExactNumColsMatching)
796 unsigned int sumNumRowsLocals = 0;
797 unsigned int sumNumCols = 0;
798 for (
unsigned int i = 0; i < matrices.size(); ++i) {
799 sumNumRowsLocals += matrices[i]->numRowsLocal();
800 sumNumCols += matrices[i]->numCols();
807 unsigned int cumulativeRowId = 0;
808 unsigned int cumulativeColId = 0;
809 for (
unsigned int i = 0; i < matrices.size(); ++i) {
810 unsigned int nRows = matrices[i]->numRowsLocal();
811 unsigned int nCols = matrices[i]->numCols();
812 for (
unsigned int rowId = 0; rowId < nRows; ++rowId) {
813 for (
unsigned int colId = 0; colId < nCols; ++colId) {
814 (*this)(initialTargetRowId + cumulativeRowId + rowId, initialTargetColId + cumulativeColId + colId) = (*(matrices[i]))(rowId,colId);
817 cumulativeRowId += nRows;
818 cumulativeColId += nCols;
826 unsigned int initialTargetRowId,
827 unsigned int initialTargetColId,
828 const std::vector<const GslMatrix* >& matrices,
829 bool checkForExactNumRowsMatching,
830 bool checkForExactNumColsMatching)
832 unsigned int sumNumCols = 0;
833 for (
unsigned int i = 0; i < matrices.size(); ++i) {
836 sumNumCols += matrices[i]->numCols();
841 unsigned int cumulativeColId = 0;
842 for (
unsigned int i = 0; i < matrices.size(); ++i) {
843 unsigned int nRows = matrices[i]->numRowsLocal();
844 unsigned int nCols = matrices[i]->numCols();
845 for (
unsigned int rowId = 0; rowId < nRows; ++rowId) {
846 for (
unsigned int colId = 0; colId < nCols; ++colId) {
847 (*this)(initialTargetRowId + rowId, initialTargetColId + cumulativeColId + colId) = (*(matrices[i]))(rowId,colId);
850 cumulativeColId += nCols;
858 unsigned int initialTargetRowId,
859 unsigned int initialTargetColId,
860 const std::vector<GslMatrix* >& matrices,
861 bool checkForExactNumRowsMatching,
862 bool checkForExactNumColsMatching)
864 unsigned int sumNumCols = 0;
865 for (
unsigned int i = 0; i < matrices.size(); ++i) {
868 sumNumCols += matrices[i]->numCols();
873 unsigned int cumulativeColId = 0;
874 for (
unsigned int i = 0; i < matrices.size(); ++i) {
875 unsigned int nRows = matrices[i]->numRowsLocal();
876 unsigned int nCols = matrices[i]->numCols();
877 for (
unsigned int rowId = 0; rowId < nRows; ++rowId) {
878 for (
unsigned int colId = 0; colId < nCols; ++colId) {
879 (*this)(initialTargetRowId + rowId, initialTargetColId + cumulativeColId + colId) = (*(matrices[i]))(rowId,colId);
882 cumulativeColId += nCols;
890 unsigned int initialTargetRowId,
891 unsigned int initialTargetColId,
892 const std::vector<const GslMatrix* >& matrices,
893 bool checkForExactNumRowsMatching,
894 bool checkForExactNumColsMatching)
896 unsigned int sumNumRows = 0;
897 for (
unsigned int i = 0; i < matrices.size(); ++i) {
900 sumNumRows += matrices[i]->numRowsLocal();
905 unsigned int cumulativeRowId = 0;
906 for (
unsigned int i = 0; i < matrices.size(); ++i) {
907 unsigned int nRows = matrices[i]->numRowsLocal();
908 unsigned int nCols = matrices[i]->numCols();
909 for (
unsigned int rowId = 0; rowId < nRows; ++rowId) {
910 for (
unsigned int colId = 0; colId < nCols; ++colId) {
911 (*this)(initialTargetRowId + cumulativeRowId + rowId, initialTargetColId + colId) = (*(matrices[i]))(rowId,colId);
914 cumulativeRowId += nRows;
922 unsigned int initialTargetRowId,
923 unsigned int initialTargetColId,
924 const std::vector<GslMatrix* >& matrices,
925 bool checkForExactNumRowsMatching,
926 bool checkForExactNumColsMatching)
928 unsigned int sumNumRows = 0;
929 for (
unsigned int i = 0; i < matrices.size(); ++i) {
932 sumNumRows += matrices[i]->numRowsLocal();
937 unsigned int cumulativeRowId = 0;
938 for (
unsigned int i = 0; i < matrices.size(); ++i) {
939 unsigned int nRows = matrices[i]->numRowsLocal();
940 unsigned int nCols = matrices[i]->numCols();
941 for (
unsigned int rowId = 0; rowId < nRows; ++rowId) {
942 for (
unsigned int colId = 0; colId < nCols; ++colId) {
943 (*this)(initialTargetRowId + cumulativeRowId + rowId, initialTargetColId + colId) = (*(matrices[i]))(rowId,colId);
946 cumulativeRowId += nRows;
954 unsigned int initialTargetRowId,
955 unsigned int initialTargetColId,
958 bool checkForExactNumRowsMatching,
959 bool checkForExactNumColsMatching)
966 for (
unsigned int rowId1 = 0; rowId1 < mat1.
numRowsLocal(); ++rowId1) {
967 for (
unsigned int colId1 = 0; colId1 < mat1.
numCols(); ++colId1) {
968 double multiplicativeFactor = mat1(rowId1,colId1);
969 unsigned int targetRowId = rowId1 * mat2.
numRowsLocal();
970 unsigned int targetColId = colId1 * mat2.
numCols();
971 for (
unsigned int rowId2 = 0; rowId2 < mat2.
numRowsLocal(); ++rowId2) {
972 for (
unsigned int colId2 = 0; colId2 < mat2.
numCols(); ++colId2) {
973 (*this)(initialTargetRowId + targetRowId + rowId2, initialTargetColId + targetColId + colId2) = multiplicativeFactor * mat2(rowId2,colId2);
984 unsigned int initialTargetRowId,
985 unsigned int initialTargetColId,
988 bool checkForExactNumRowsMatching,
989 bool checkForExactNumColsMatching)
996 for (
unsigned int rowId1 = 0; rowId1 < mat1.
numRowsLocal(); ++rowId1) {
997 for (
unsigned int colId1 = 0; colId1 < mat1.
numCols(); ++colId1) {
998 double multiplicativeFactor = mat1(rowId1,colId1);
999 unsigned int targetRowId = rowId1 * vec2.
sizeLocal();
1000 unsigned int targetColId = colId1 * 1;
1001 for (
unsigned int rowId2 = 0; rowId2 < vec2.
sizeLocal(); ++rowId2) {
1002 for (
unsigned int colId2 = 0; colId2 < 1; ++colId2) {
1003 (*this)(initialTargetRowId + targetRowId + rowId2, initialTargetColId + targetColId + colId2) = multiplicativeFactor * vec2[rowId2];
1015 unsigned int initialTargetRowId,
1016 unsigned int initialTargetColId,
1018 bool checkForExactNumRowsMatching,
1019 bool checkForExactNumColsMatching)
1022 unsigned int nCols = mat.
numCols();
1028 for (
unsigned int row = 0; row < nRows; ++row) {
1029 for (
unsigned int col = 0; col < nCols; ++col) {
1030 (*this)(initialTargetRowId + col, initialTargetColId + row) = mat(row,col);
1046 <<
": before 'this->invertMultiply()'"
1052 <<
": after 'this->invertMultiply()'"
1058 <<
": before 'gsl_linalg_LU_det()'"
1064 <<
": after 'gsl_linalg_LU_det()'"
1070 <<
": after 'gsl_linalg_LU_lndet()'"
1087 <<
": before 'this->invertMultiply()'"
1093 <<
": after 'this->invertMultiply()'"
1099 <<
": before 'gsl_linalg_LU_det()'"
1105 <<
": after 'gsl_linalg_LU_det()'"
1111 <<
": before 'gsl_linalg_LU_lndet()'"
1127 if (relativeVec[0] > 0.) {
1128 relativeVec = (1./relativeVec[0])*relativeVec;
1131 unsigned int rankValue = 0;
1132 for (
unsigned int i = 0; i < relativeVec.
sizeLocal(); ++i) {
1133 if (( (*
m_svdSvec)[i] >= absoluteZeroThreshold ) &&
1134 ( relativeVec [i] >= relativeZeroThreshold )) {
1142 <<
", this->numCols() = " << this->
numCols()
1143 <<
", absoluteZeroThreshold = " << absoluteZeroThreshold
1144 <<
", relativeZeroThreshold = " << relativeZeroThreshold
1145 <<
", rankValue = " << rankValue
1147 <<
", relativeVec = " << relativeVec
1175 unsigned int sizeX = this->
numCols();
1177 for (
unsigned int i = 0; i < sizeY; ++i) {
1179 for (
unsigned int j = 0; j < sizeX; ++j) {
1180 value += (*this)(i,j)*x[j];
1210 const unsigned int p_s = this->
numCols();
1211 const unsigned int n_s = X.
numCols();
1213 for (
unsigned int k=0;
k<p_s;
k++)
1214 for (
unsigned int j=0; j<n_s; j++)
1216 for (
unsigned int i=0; i<m_s; i++)
1217 Y(i,j) += (*this)(i,
k) * X(
k,j);
1256 std::cout <<
"In GslMatrix::invertMultiply()"
1257 <<
": before LU decomposition, m_LU = ";
1258 gsl_matrix_fprintf(stdout,
m_LU,
"%f");
1259 std::cout << std::endl;
1262 gsl_error_handler_t* oldHandler;
1263 oldHandler = gsl_set_error_handler_off();
1266 <<
": before 'gsl_linalg_LU_decomp()'"
1271 std::cerr <<
"In GslMatrix::invertMultiply()"
1272 <<
", after gsl_linalg_LU_decomp()"
1273 <<
": iRC = " << iRC
1274 <<
", gsl error message = " << gsl_strerror(iRC)
1277 gsl_set_error_handler(oldHandler);
1280 <<
": after 'gsl_linalg_LU_decomp()'"
1281 <<
", IRC = " << iRC
1287 std::cout <<
"In GslMatrix::invertMultiply()"
1288 <<
": after LU decomposition, m_LU = ";
1289 gsl_matrix_fprintf(stdout,
m_LU,
"%f");
1290 std::cout << std::endl;
1294 gsl_error_handler_t* oldHandler;
1295 oldHandler = gsl_set_error_handler_off();
1298 <<
": before 'gsl_linalg_LU_solve()'"
1304 std::cerr <<
"In GslMatrix::invertMultiply()"
1305 <<
", after gsl_linalg_LU_solve()"
1306 <<
": iRC = " << iRC
1307 <<
", gsl error message = " << gsl_strerror(iRC)
1310 gsl_set_error_handler(oldHandler);
1313 <<
": after 'gsl_linalg_LU_solve()'"
1314 <<
", IRC = " << iRC
1325 std::cout <<
"In GslMatrix::invertMultiply()"
1326 <<
": ||b - Ax||_2 = " << tmpVec.
norm2()
1327 <<
": ||b - Ax||_2/||b||_2 = " << tmpVec.
norm2()/b.
norm2()
1348 "Matrices B and X are incompatible");
1350 "Matrices B and X are incompatible");
1352 "This and X matrices are incompatible");
1358 for(
unsigned int j = 0; j < B.
numCols(); ++j )
1395 if (
m_LU == NULL ) {
1422 unsigned int n = eigenValues.
sizeLocal();
1430 if (eigenVectors == NULL) {
1431 gsl_eigen_symm_workspace* w = gsl_eigen_symm_alloc((
size_t) n);
1432 gsl_eigen_symm(
m_mat,eigenValues.
data(),w);
1433 gsl_eigen_symm_free(w);
1436 gsl_eigen_symmv_workspace* w = gsl_eigen_symmv_alloc((
size_t) n);
1438 gsl_eigen_symmv_sort(eigenValues.
data(),eigenVectors->
m_mat,GSL_EIGEN_SORT_VAL_ASC);
1439 gsl_eigen_symmv_free(w);
1450 unsigned int n = eigenVector.
sizeLocal();
1461 const unsigned int max_num_iterations = 10000;
1462 const double tolerance = 1.0e-13;
1476 for(
unsigned int k = 0;
k < max_num_iterations; ++
k )
1483 index = (w.
abs()).getMaxValueIndex();
1487 z = (1.0/lambda) * w;
1491 residual = ( (*this)*z - lambda*z ).norm2();
1493 if( residual < tolerance )
1495 eigenValue = lambda;
1518 unsigned int n = eigenVector.
sizeLocal();
1529 const unsigned int max_num_iterations = 1000;
1530 const double tolerance = 1.0e-13;
1542 double one_over_lambda;
1545 for(
unsigned int k = 0;
k < max_num_iterations; ++
k )
1547 w = (*this).invertMultiplyForceLU( z );
1552 index = (w.
abs()).getMaxValueIndex();
1555 one_over_lambda = w[index];
1557 lambda = 1.0/one_over_lambda;
1563 residual = ( (*this)*z - lambda*z ).norm2();
1565 if( residual < tolerance )
1567 eigenValue = lambda;
1594 gsl_vector* gsl_column = gsl_vector_alloc( column.
sizeLocal() );
1596 int error_code = gsl_matrix_get_col( gsl_column,
m_mat, column_num );
1600 for(
unsigned int i = 0; i < column.
sizeLocal(); ++i )
1602 column[i] = gsl_vector_get( gsl_column, i );
1606 gsl_vector_free( gsl_column );
1620 gsl_vector* gsl_row = gsl_vector_alloc( row.
sizeLocal() );
1622 int error_code = gsl_matrix_get_row( gsl_row,
m_mat, row_num );
1626 for(
unsigned int i = 0; i < row.
sizeLocal(); ++i )
1628 row[i] = gsl_vector_get( gsl_row, i );
1632 gsl_vector_free( gsl_row );
1645 this->
getRow( row_num, row );
1672 gsl_vector* gsl_row = row.
data();
1674 int error_code = gsl_matrix_set_row(
m_mat, row_num, gsl_row );
1689 gsl_vector* gsl_column = column.
data();
1691 int error_code = gsl_matrix_set_col(
m_mat, column_num, gsl_column );
1704 "GslMatrix::mpiSum()",
1705 "local and global matrices incompatible");
1709 std::vector<double> local( size, 0.0 );
1710 std::vector<double> global( size, 0.0 );
1713 for(
unsigned int i = 0; i < this->
numRowsLocal(); i++ )
1715 for(
unsigned int j = 0; j < this->
numCols(); j++ )
1719 local[
k] = (*this)(i,j);
1724 "GslMatrix::mpiSum()",
1725 "failed MPI.Allreduce()");
1727 for(
unsigned int i = 0; i < this->
numRowsLocal(); i++ )
1729 for(
unsigned int j = 0; j < this->
numCols(); j++ )
1733 M_global(i,j) = global[
k];
1756 for (
unsigned int colId = 0; colId < y1Mat.
numCols(); ++colId) {
1775 unsigned int nCols = this->
numCols();
1778 for (
unsigned int i = 0; i < nRows; ++i) {
1779 for (
unsigned int j = 0; j < nCols; ++j) {
1783 if (i != (nRows-1)) os <<
"; ";
1788 for (
unsigned int i = 0; i < nRows; ++i) {
1789 for (
unsigned int j = 0; j < nCols; ++j) {
1802 const std::string& varNamePrefix,
1803 const std::string& fileName,
1804 const std::string& fileType,
1805 const std::set<unsigned int>& allowedSubEnvIds)
const
1818 unsigned int nCols = this->
numCols();
1825 for (
unsigned int i = 0; i < nRows; ++i) {
1826 for (
unsigned int j = 0; j < nCols; ++j) {
1827 *filePtrSet.
ofsVar << (*this)(i,j)
1830 *filePtrSet.
ofsVar <<
"\n";
1832 *filePtrSet.
ofsVar <<
"];\n";
1842 const std::string& fileName,
1843 const std::string& fileType,
1844 const std::set<unsigned int>& allowedSubEnvIds)
1860 unsigned int idOfMyFirstLine = 1;
1861 unsigned int idOfMyLastLine = nRowsLocal;
1862 unsigned int nCols = this->
numCols();
1866 std::string tmpString;
1869 *filePtrSet.
ifsVar >> tmpString;
1873 *filePtrSet.
ifsVar >> tmpString;
1878 *filePtrSet.
ifsVar >> tmpString;
1880 unsigned int posInTmpString = 6;
1883 char nRowsString[tmpString.size()-posInTmpString+1];
1884 unsigned int posInRowsString = 0;
1887 nRowsString[posInRowsString++] = tmpString[posInTmpString++];
1888 }
while (tmpString[posInTmpString] !=
',');
1889 nRowsString[posInRowsString] =
'\0';
1893 char nColsString[tmpString.size()-posInTmpString+1];
1894 unsigned int posInColsString = 0;
1897 nColsString[posInColsString++] = tmpString[posInTmpString++];
1898 }
while (tmpString[posInTmpString] !=
')');
1899 nColsString[posInColsString] =
'\0';
1902 unsigned int numRowsInFile = (
unsigned int) strtod(nRowsString,NULL);
1903 unsigned int numColsInFile = (
unsigned int) strtod(nColsString,NULL);
1907 <<
", numRowsInFile = " << numRowsInFile
1908 <<
", numColsInFile = " << numColsInFile
1909 <<
", nRowsLocal = " << nRowsLocal
1910 <<
", nCols = " << nCols
1918 queso_require_equal_to_msg(numColsInFile, nCols,
"number of parameters of vec in file is different than number of parameters in this vec object");
1921 unsigned int maxCharsPerLine = 64*nCols;
1923 unsigned int lineId = 0;
1924 while (lineId < idOfMyFirstLine) {
1925 filePtrSet.
ifsVar->ignore(maxCharsPerLine,
'\n');
1931 <<
": beginning to read input actual data"
1938 *filePtrSet.
ifsVar >> tmpString;
1942 *filePtrSet.
ifsVar >> tmpString;
1947 std::streampos tmpPos = filePtrSet.
ifsVar->tellg();
1948 filePtrSet.
ifsVar->seekg(tmpPos+(std::streampos)2);
1952 <<
": beginning to read lines with numbers only"
1953 <<
", lineId = " << lineId
1954 <<
", idOfMyFirstLine = " << idOfMyFirstLine
1955 <<
", idOfMyLastLine = " << idOfMyLastLine
1960 while (lineId <= idOfMyLastLine) {
1961 for (
unsigned int j = 0; j < nCols; ++j) {
1962 *filePtrSet.
ifsVar >> tmpRead;
1963 (*this)(lineId-idOfMyFirstLine,j) = tmpRead;
1997 unsigned int m1Cols = m1.
numCols();
1999 unsigned int m2Cols = m2.
numCols();
2007 unsigned int commonSize = m1Cols;
2008 for (
unsigned int row1 = 0; row1 < m1Rows; ++row1) {
2009 for (
unsigned int col2 = 0; col2 < m2Cols; ++col2) {
2011 for (
unsigned int k = 0;
k < commonSize; ++
k) {
2012 result += m1(row1,
k)*m2(
k,col2);
2014 mat(row1,col2) = result;
2042 for (
unsigned int i = 0; i < nRows; ++i) {
2043 double value1 = v1[i];
2044 for (
unsigned int j = 0; j < nCols; ++j) {
2045 answer(i,j) = value1*v2[j];
2056 unsigned int mCols = mat.
numCols();
2063 for (
unsigned int i = 0; i < mRows; ++i) {
2064 double vecValue = vec[i];
2065 for (
unsigned int j = 0; j < mCols; ++j) {
2066 answer(i,j) *= vecValue;
2077 unsigned int mCols = mat.
numCols();
2084 for (
unsigned int j = 0; j < mCols; ++j) {
2085 double vecValue = vec[j];
2086 for (
unsigned int i = 0; i < mRows; ++i) {
2087 answer(i,j) *= vecValue;
double norm2() const
Returns the 2-norm (Euclidean norm) of the vector.
GslMatrix operator+(const GslMatrix &m1, const GslMatrix &m2)
void zeroUpper(bool includeDiagonal=false)
This function sets all the entries above the main diagonal of this matrix to zero.
double normMax() const
Returns the Frobenius norm of this matrix.
void zeroLower(bool includeDiagonal=false)
This function sets all the entries bellow the main diagonal of this matrix to zero.
void eigen(GslVector &eigenValues, GslMatrix *eigenVectors) const
This function computes the eigenvalues of a real symmetric matrix.
GslMatrix operator-(const GslMatrix &m1, const GslMatrix &m2)
gsl_matrix * m_mat
GSL matrix, also referred to as this matrix.
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.
bool m_isSingular
Indicates whether or not this matrix is singular.
int worldRank() const
Returns the process world rank.
const int UQ_MATRIX_IS_NOT_POS_DEFINITE_RC
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
#define queso_require_greater_equal_msg(expr1, expr2, msg)
const GslMatrix & svdMatV() const
This function calls private member GslMatrix::internalSvd() to set a N-by-N orthogonal square matrix ...
void filterLargeValues(double thresholdValue)
This function sets to zero (filters) all entries of this matrix which are greater than thresholdValue...
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.
int subRank() const
Access function for sub-rank.
std::ifstream * ifsVar
Provides a stream interface to read data from files.
void largestEigen(double &eigenValue, GslVector &eigenVector) const
This function finds largest eigenvalue, namely eigenValue, of this matrix and its corresponding eigen...
int internalSvd() const
This function factorizes the M-by-N matrix A into the singular value decomposition A = U S V^T for M ...
This (virtual) class sets up the environment underlying the use of the QUESO library by an executable...
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.
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...
gsl_matrix * data()
Returns this matrix.
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_printHorizontally
Flag for either or not print this matrix.
GslMatrix & operator/=(double a)
Stores in this the coordinate-wise division of this by a.
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).
#define UQ_RC_MACRO(macroIRc, givenRank, where, what, retValue)
double lnDeterminant() const
Calculates the ln(determinant) of this matrix.
GslVector * m_svdSvec
m_svdSvec stores the diagonal of the N-by-N diagonal matrix of singular values S after the singular v...
double max() const
Returns the maximum element value of the matrix.
GslVector multiply(const GslVector &x) const
This function multiplies this matrix by vector x and returns the resulting vector.
GslMatrix operator*(double a, const GslMatrix &mat)
#define queso_require_equal_to_msg(expr1, expr2, msg)
bool openInputFile(const std::string &fileName, const std::string &fileType, const std::set< unsigned int > &allowedSubEnvIds, FilePtrSetStruct &filePtrSet) const
Opens an input file.
#define queso_require_less_msg(expr1, expr2, msg)
void mpiSum(const MpiComm &comm, GslMatrix &M_global) const
const BaseEnvironment & env() const
bool m_inDebugMode
Flag for either or not QUESO is in debug mode.
A class for partitioning vectors and matrices.
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...
const GslMatrix & svdMatU() const
This function calls private member GslMatrix::internalSvd() to set a M-by-N orthogonal matrix U of th...
GslMatrix * m_svdUmat
m_svdUmat stores the M-by-N orthogonal matrix U after the singular value decomposition of a matrix...
GslVector invertMultiplyForceLU(const GslVector &b) const
This function calculates the inverse of this matrix and multiplies it with vector b...
Struct for handling data input and output from files.
void setColumn(const unsigned int column_num, const GslVector &column)
This function copies vector column into the column_num-th column of this matrix.
unsigned int numRowsGlobal() const
Returns the global row dimension of this matrix.
const BaseEnvironment & env() const
const Map & m_map
Mapping variable.
Class for vector operations using GSL library.
void resetLU()
In this function resets the LU decomposition of this matrix, as well as deletes the private member po...
unsigned int numOfProcsForStorage() const
#define UQ_FATAL_RC_MACRO(macroIRc, givenRank, where, what)
std::ostream & operator<<(std::ostream &os, const BaseEnvironment &obj)
unsigned int numRowsLocal() const
Returns the local row dimension of this matrix.
int m_signum
m_signum stores the sign of the permutation of the LU decomposition PA = LU.
void cwSet(double value)
Component-wise sets all values to this with value.
const BaseEnvironment & m_env
QUESO environment variable.
GslVector invertMultiply(const GslVector &b) const
This function calculates the inverse of this matrix and multiplies it with vector b...
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.
GslMatrix * m_svdVmat
m_svdVmat stores the N-by-N orthogonal square matrix V after the singular value decomposition of a ma...
GslMatrix transpose() const
This function calculated the transpose of this matrix (square).
int chol()
Computes Cholesky factorization of a real symmetric positive definite matrix this.
GslMatrix & operator*=(double a)
Stores in this the coordinate-wise multiplication of this and a.
gsl_matrix * m_LU
GSL matrix for the LU decomposition of m_mat.
#define queso_require_greater_msg(expr1, expr2, msg)
double m_lnDeterminant
The natural logarithm of the determinant of this matrix.
void closeFile(FilePtrSetStruct &filePtrSet, const std::string &fileType) const
Closes the file.
int fullRank() const
Returns the process full rank.
void matlabLinearInterpExtrap(const GslVector &x1Vec, const GslVector &y1Vec, const GslVector &x2Vec)
GslMatrix & operator-=(const GslMatrix &rhs)
Stores in this the coordinate-wise subtraction of this by rhs.
std::ofstream * ofsVar
Provides a stream interface to write data to files.
GslMatrix inverse() const
This function calculated the inverse of this matrix (square).
void smallestEigen(double &eigenValue, GslVector &eigenVector) const
This function finds smallest eigenvalue, namely eigenValue, of this matrix and its corresponding eige...
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...
#define queso_require_not_equal_to_msg(expr1, expr2, msg)
void copy(const GslMatrix &src)
In this function this matrix receives a copy of matrix src.
unsigned int checkingLevel() const
Access function to private attribute m_checkingLevel.
GslMatrix & operator+=(const GslMatrix &rhs)
Stores in this the coordinate-wise addition of this and rhs.
#define queso_require_less_equal_msg(expr1, expr2, msg)
const int UQ_MATRIX_SVD_FAILED_RC
unsigned int displayVerbosity() const
GslMatrix leftDiagScaling(const GslVector &vec, const GslMatrix &mat)
unsigned int sizeLocal() const
Returns the length of this vector.
void cwExtract(unsigned int rowId, unsigned int colId, GslMatrix &mat) const
#define queso_require_msg(asserted, msg)
The QUESO MPI Communicator Class.
void cwSet(double value)
Component-wise set all values to this with value.
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.
GslMatrix * m_inverse
Inverse matrix of this.
Class for matrix operations (virtual).
void setRow(const unsigned int row_num, const GslVector &row)
This function copies vector row into the row_num-th column of this matrix.
void matlabLinearInterpExtrap(const GslVector &x1Vec, const GslMatrix &y1Mat, const GslVector &x2Vec)
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...
double determinant() const
Calculates the determinant of this matrix.
const MpiComm & Comm() const
Access function for MpiComm communicator.
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.
GslMatrix rightDiagScaling(const GslMatrix &mat, const GslVector &vec)
double normFrob() const
Returns the Frobenius norm of this matrix.
void subReadContents(const std::string &fileName, const std::string &fileType, const std::set< unsigned int > &allowedSubEnvIds)
Read contents of subenvironment from file fileName.
double m_determinant
The determinant of this matrix.
void filterSmallValues(double thresholdValue)
This function sets to zero (filters) all entries of this matrix which are smaller than thresholdValue...
GslMatrix()
Default Constructor.
void print(std::ostream &os) const
Print method. Defines the behavior of the ostream << operator inherited from the Object class...
gsl_permutation * m_permutation
The permutation matrix of a LU decomposition.
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.
GslMatrix matrixProduct(const GslVector &v1, const GslVector &v2)
GslMatrix & operator=(const GslMatrix &rhs)
Copies values from matrix rhs to this.
Map * m_svdColMap
Mapping for matrices involved in the singular value decomposition (svd) routine.
gsl_vector * data() const
unsigned int rank(double absoluteZeroThreshold, double relativeZeroThreshold) const
This function returns the number of singular values of this matrix (rank).
unsigned int numCols() const
Returns the column dimension 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...
GslVector abs() const
This function returns absolute value of elements in this.
Class for matrix operations using GSL library.