queso-0.57.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
QUESO::GslMatrix Class Reference

Class for matrix operations using GSL library. More...

#include <GslMatrix.h>

Inheritance diagram for QUESO::GslMatrix:
QUESO::Matrix

Public Member Functions

Constructor/Destructor methods
 GslMatrix (const BaseEnvironment &env, const Map &map, unsigned int numCols)
 Shaped Constructor: creates a shaped matrix with numCols columns. More...
 
 GslMatrix (const BaseEnvironment &env, const Map &map, double diagValue)
 Shaped Constructor: creates a square matrix with size map.NumGlobalElements() and diagonal values all equal to diagValue. More...
 
 GslMatrix (const GslVector &v, double diagValue)
 Shaped Constructor: creates a square matrix with size v.sizeLocal() and diagonal values all equal to diagValue. More...
 
 GslMatrix (const GslVector &v)
 Shaped Constructor: creates a square matrix with size v.sizeLocal(). More...
 
 GslMatrix (const GslMatrix &B)
 Shaped Constructor: creates a matrix with B.numCols() columns and B.numRowsLocal() rows. More...
 
 ~GslMatrix ()
 Destructor. More...
 
Set methods
GslMatrixoperator= (const GslMatrix &rhs)
 Copies values from matrix rhs to this. More...
 
GslMatrixoperator*= (double a)
 Stores in this the coordinate-wise multiplication of this and a. More...
 
GslMatrixoperator/= (double a)
 Stores in this the coordinate-wise division of this by a. More...
 
GslMatrixoperator+= (const GslMatrix &rhs)
 Stores in this the coordinate-wise addition of this and rhs. More...
 
GslMatrixoperator-= (const GslMatrix &rhs)
 Stores in this the coordinate-wise subtraction of this by rhs. More...
 
Accessor methods
double & operator() (unsigned int i, unsigned int j)
 Element access method (non-const). More...
 
const double & operator() (unsigned int i, unsigned int j) const
 Element access method (const). More...
 
Attribute methods
unsigned int numRowsLocal () const
 Returns the local row dimension of this matrix. More...
 
unsigned int numRowsGlobal () const
 Returns the global row dimension of this matrix. More...
 
unsigned int numCols () const
 Returns the column dimension of this matrix. More...
 
double max () const
 Returns the maximum element value of the matrix. More...
 
unsigned int rank (double absoluteZeroThreshold, double relativeZeroThreshold) const
 This function returns the number of singular values of this matrix (rank). More...
 
GslMatrix transpose () const
 This function calculated the transpose of this matrix (square). More...
 
GslMatrix inverse () const
 This function calculated the inverse of this matrix (square). More...
 
double determinant () const
 Calculates the determinant of this matrix. More...
 
double lnDeterminant () const
 Calculates the ln(determinant) of this matrix. More...
 
Norm methods
double normFrob () const
 Returns the Frobenius norm of this matrix. More...
 
double normMax () const
 Returns the Frobenius norm of this matrix. More...
 
Mathematical methods
int chol ()
 Computes Cholesky factorization of a real symmetric positive definite matrix this. More...
 
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 internalSvd to compute the singular values of this. More...
 
const GslMatrixsvdMatU () const
 This function calls private member GslMatrix::internalSvd() to set a M-by-N orthogonal matrix U of the singular value decomposition (svd) of a general rectangular M-by-N matrix A. More...
 
const GslMatrixsvdMatV () const
 This function calls private member GslMatrix::internalSvd() to set a N-by-N orthogonal square matrix V of the singular value decomposition (svd) of a general rectangular M-by-N matrix A. More...
 
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). More...
 
int svdSolve (const GslMatrix &rhsMat, GslMatrix &solMat) 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=solMat, b=rhsMat). More...
 
GslVector multiply (const GslVector &x) const
 This function multiplies this matrix by vector x and returns the resulting vector. More...
 
void multiply (const GslVector &x, GslVector &y) const
 This function multiplies this matrix by vector x and fills. More...
 
GslMatrix multiply (const GslMatrix &X) const
 This function multiplies this matrix by matrix X and returns the resulting matrix. More...
 
void multiply (const GslMatrix &X, GslMatrix &Y) const
 This function multiplies this matrix by matrix X and fills. More...
 
GslVector invertMultiply (const GslVector &b) const
 This function calculates the inverse of this matrix and multiplies it with vector b. More...
 
void invertMultiply (const GslVector &b, GslVector &x) const
 This function calculates the inverse of this matrix, multiplies it with vector b and stores the result in vector x. More...
 
GslMatrix invertMultiply (const GslMatrix &B) const
 This function calculates the inverse of this matrix and multiplies it with matrix B. More...
 
void invertMultiply (const GslMatrix &B, GslMatrix &X) const
 This function calculates the inverse of this matrix, multiplies it with matrix B and stores the result in matrix X. More...
 
GslVector invertMultiplyForceLU (const GslVector &b) const
 This function calculates the inverse of this matrix and multiplies it with vector b. More...
 
void invertMultiplyForceLU (const GslVector &b, GslVector &x) const
 This function calculates the inverse of this matrix, multiplies it with vector b and stores the result in vector x. More...
 
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. More...
 
GslVector getColumn (const unsigned int column_num) const
 This function gets the column_num-th column of this matrix. More...
 
void setColumn (const unsigned int column_num, const GslVector &column)
 This function copies vector column into the column_num-th column of this matrix. More...
 
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. More...
 
GslVector getRow (const unsigned int row_num) const
 This function gets the row_num-th column of this matrix. More...
 
void setRow (const unsigned int row_num, const GslVector &row)
 This function copies vector row into the row_num-th column of this matrix. More...
 
void eigen (GslVector &eigenValues, GslMatrix *eigenVectors) const
 This function computes the eigenvalues of a real symmetric matrix. More...
 
void largestEigen (double &eigenValue, GslVector &eigenVector) const
 This function finds largest eigenvalue, namely eigenValue, of this matrix and its corresponding eigenvector, namely eigenVector. More...
 
void smallestEigen (double &eigenValue, GslVector &eigenVector) const
 This function finds smallest eigenvalue, namely eigenValue, of this matrix and its corresponding eigenvector, namely eigenVector. More...
 
Get/Set methods
void cwSet (double value)
 Component-wise set all values to this with value. More...
 
void cwSet (unsigned int rowId, unsigned int colId, const GslMatrix &mat)
 Set the components of which positions are greater than (rowId,colId) with the value of mat(rowId,colId). More...
 
void cwExtract (unsigned int rowId, unsigned int colId, GslMatrix &mat) const
 
void zeroLower (bool includeDiagonal=false)
 This function sets all the entries bellow the main diagonal of this matrix to zero. More...
 
void zeroUpper (bool includeDiagonal=false)
 This function sets all the entries above the main diagonal of this matrix to zero. More...
 
void filterSmallValues (double thresholdValue)
 This function sets to zero (filters) all entries of this matrix which are smaller than thresholdValue. More...
 
void filterLargeValues (double thresholdValue)
 This function sets to zero (filters) all entries of this matrix which are greater than thresholdValue. More...
 
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. More...
 
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. More...
 
void fillWithBlocksDiagonally (unsigned int rowId, unsigned int colId, const std::vector< GslMatrix * > &matrices, bool checkForExactNumRowsMatching, bool checkForExactNumColsMatching)
 This function fills this matrix diagonally with block matrices. More...
 
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. More...
 
void fillWithBlocksHorizontally (unsigned int rowId, unsigned int colId, const std::vector< GslMatrix * > &matrices, bool checkForExactNumRowsMatching, bool checkForExactNumColsMatching)
 This function fills this matrix horizontally with const block. More...
 
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. More...
 
void fillWithBlocksVertically (unsigned int rowId, unsigned int colId, const std::vector< GslMatrix * > &matrices, bool checkForExactNumRowsMatching, bool checkForExactNumColsMatching)
 This function fills this matrix vertically with block matrices. More...
 
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. More...
 
void fillWithTensorProduct (unsigned int rowId, unsigned int colId, const GslMatrix &mat1, const GslVector &vec2, bool checkForExactNumRowsMatching, bool checkForExactNumColsMatching)
 This function calculates the tensor product of matrix mat1 and vector vec2 and stores it in this matrix. More...
 
Miscellaneous methods
gsl_matrix * data ()
 Returns this matrix. More...
 
void mpiSum (const MpiComm &comm, GslMatrix &M_global) const
 
void matlabLinearInterpExtrap (const GslVector &x1Vec, const GslMatrix &y1Mat, const GslVector &x2Vec)
 
I/O methods
void print (std::ostream &os) const
 Print method. Defines the behavior of the ostream << operator inherited from the Object class. More...
 
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. More...
 
void subReadContents (const std::string &fileName, const std::string &fileType, const std::set< unsigned int > &allowedSubEnvIds)
 Read contents of subenvironment from file fileName. More...
 
- Public Member Functions inherited from QUESO::Matrix
 Matrix (const BaseEnvironment &env, const Map &map)
 Shaped constructor. More...
 
virtual ~Matrix ()
 Virtual Destructor. More...
 
const BaseEnvironmentenv () const
 
const Mapmap () const
 
unsigned int numOfProcsForStorage () const
 
void setPrintHorizontally (bool value) const
 Determines whether the matrix should be printed horizontally. More...
 
bool getPrintHorizontally () const
 Checks if matrix will be is printed horizontally. More...
 
void setInDebugMode (bool value) const
 Determines whether QUESO will run through this class in debug mode. More...
 
bool getInDebugMode () const
 Checks if QUESO will run through this class in debug mode. More...
 

Private Member Functions

 GslMatrix ()
 Default Constructor. More...
 
void copy (const GslMatrix &src)
 In this function this matrix receives a copy of matrix src. More...
 
void resetLU ()
 In this function resets the LU decomposition of this matrix, as well as deletes the private member pointers, if existing. More...
 
int internalSvd () const
 This function factorizes the M-by-N matrix A into the singular value decomposition A = U S V^T for M >= N. On output the matrix A is replaced by U. More...
 

Private Attributes

gsl_matrix * m_mat
 GSL matrix, also referred to as this matrix. More...
 
gsl_matrix * m_LU
 GSL matrix for the LU decomposition of m_mat. More...
 
GslMatrixm_inverse
 Inverse matrix of this. More...
 
Mapm_svdColMap
 Mapping for matrices involved in the singular value decomposition (svd) routine. More...
 
GslMatrixm_svdUmat
 m_svdUmat stores the M-by-N orthogonal matrix U after the singular value decomposition of a matrix. More...
 
GslVectorm_svdSvec
 m_svdSvec stores the diagonal of the N-by-N diagonal matrix of singular values S after the singular value decomposition of a matrix. More...
 
GslMatrixm_svdVmat
 m_svdVmat stores the N-by-N orthogonal square matrix V after the singular value decomposition of a matrix. More...
 
GslMatrixm_svdVTmat
 m_svdVmatT stores the transpose of N-by-N orthogonal square matrix V, namely V^T, after the singular value decomposition of a matrix. More...
 
double m_determinant
 The determinant of this matrix. More...
 
double m_lnDeterminant
 The natural logarithm of the determinant of this matrix. More...
 
gsl_permutation * m_permutation
 The permutation matrix of a LU decomposition. More...
 
int m_signum
 m_signum stores the sign of the permutation of the LU decomposition PA = LU. More...
 
bool m_isSingular
 Indicates whether or not this matrix is singular. More...
 

Additional Inherited Members

- Protected Member Functions inherited from QUESO::Matrix
virtual void base_copy (const Matrix &src)
 Copies base data from matrix src to this matrix. More...
 
- Protected Attributes inherited from QUESO::Matrix
const BaseEnvironmentm_env
 QUESO environment variable. More...
 
const Map m_map
 Mapping variable. More...
 
const Mapm_map
 Mapping variable. More...
 
bool m_printHorizontally
 Flag for either or not print this matrix. More...
 
bool m_inDebugMode
 Flag for either or not QUESO is in debug mode. More...
 

Detailed Description

Class for matrix operations using GSL library.

This class creates and provides basic support for matrices of templated type as a specialization of Matrix using GSL matrices, which are defined by an encapsulated gsl_matrix structure.

Definition at line 49 of file GslMatrix.h.

Constructor & Destructor Documentation

QUESO::GslMatrix::GslMatrix ( const BaseEnvironment env,
const Map map,
unsigned int  numCols 
)

Shaped Constructor: creates a shaped matrix with numCols columns.

Definition at line 35 of file GslMatrix.C.

References m_mat.

39  :
40  Matrix (env,map),
41  m_mat (gsl_matrix_calloc(map.NumGlobalElements(),nCols)),
42  m_LU (NULL),
43  m_inverse (NULL),
44  m_svdColMap (NULL),
45  m_svdUmat (NULL),
46  m_svdSvec (NULL),
47  m_svdVmat (NULL),
48  m_svdVTmat (NULL),
49  m_determinant (-INFINITY),
50  m_lnDeterminant(-INFINITY),
51  m_permutation (NULL),
52  m_signum (0),
53  m_isSingular (false)
54 {
55  queso_require_msg(m_mat, "null matrix generated");
56 }
gsl_permutation * m_permutation
The permutation matrix of a LU decomposition.
Definition: GslMatrix.h:455
GslMatrix * m_svdVmat
m_svdVmat stores the N-by-N orthogonal square matrix V after the singular value decomposition of a ma...
Definition: GslMatrix.h:437
const BaseEnvironment & env() const
Definition: Matrix.C:47
int m_signum
m_signum stores the sign of the permutation of the LU decomposition PA = LU.
Definition: GslMatrix.h:461
GslMatrix * m_svdUmat
m_svdUmat stores the M-by-N orthogonal matrix U after the singular value decomposition of a matrix...
Definition: GslMatrix.h:424
bool m_isSingular
Indicates whether or not this matrix is singular.
Definition: GslMatrix.h:464
GslMatrix * m_svdVTmat
m_svdVmatT stores the transpose of N-by-N orthogonal square matrix V, namely V^T, after the singular ...
Definition: GslMatrix.h:443
Map * m_svdColMap
Mapping for matrices involved in the singular value decomposition (svd) routine.
Definition: GslMatrix.h:418
double m_lnDeterminant
The natural logarithm of the determinant of this matrix.
Definition: GslMatrix.h:449
int NumGlobalElements() const
Returns the total number of elements across all processors.
Definition: Map.C:85
gsl_matrix * m_mat
GSL matrix, also referred to as this matrix.
Definition: GslMatrix.h:409
GslVector * m_svdSvec
m_svdSvec stores the diagonal of the N-by-N diagonal matrix of singular values S after the singular v...
Definition: GslMatrix.h:431
double m_determinant
The determinant of this matrix.
Definition: GslMatrix.h:446
gsl_matrix * m_LU
GSL matrix for the LU decomposition of m_mat.
Definition: GslMatrix.h:412
const Map & map() const
Definition: Matrix.C:54
GslMatrix * m_inverse
Inverse matrix of this.
Definition: GslMatrix.h:415
Matrix()
Default constructor.
QUESO::GslMatrix::GslMatrix ( const BaseEnvironment env,
const Map map,
double  diagValue 
)

Shaped Constructor: creates a square matrix with size map.NumGlobalElements() and diagonal values all equal to diagValue.

Definition at line 59 of file GslMatrix.C.

References m_mat.

63  :
64  Matrix (env,map),
65  m_mat (gsl_matrix_calloc(map.NumGlobalElements(),map.NumGlobalElements())),
66  m_LU (NULL),
67  m_inverse (NULL),
68  m_svdColMap (NULL),
69  m_svdUmat (NULL),
70  m_svdSvec (NULL),
71  m_svdVmat (NULL),
72  m_svdVTmat (NULL),
73  m_determinant (-INFINITY),
74  m_lnDeterminant(-INFINITY),
75  m_permutation (NULL),
76  m_signum (0),
77  m_isSingular (false)
78 {
79  queso_require_msg(m_mat, "null matrix generated");
80 
81  for (unsigned int i = 0; i < m_mat->size1; ++i) {
82  (*this)(i,i) = diagValue;
83  }
84 }
gsl_permutation * m_permutation
The permutation matrix of a LU decomposition.
Definition: GslMatrix.h:455
GslMatrix * m_svdVmat
m_svdVmat stores the N-by-N orthogonal square matrix V after the singular value decomposition of a ma...
Definition: GslMatrix.h:437
const BaseEnvironment & env() const
Definition: Matrix.C:47
int m_signum
m_signum stores the sign of the permutation of the LU decomposition PA = LU.
Definition: GslMatrix.h:461
GslMatrix * m_svdUmat
m_svdUmat stores the M-by-N orthogonal matrix U after the singular value decomposition of a matrix...
Definition: GslMatrix.h:424
bool m_isSingular
Indicates whether or not this matrix is singular.
Definition: GslMatrix.h:464
GslMatrix * m_svdVTmat
m_svdVmatT stores the transpose of N-by-N orthogonal square matrix V, namely V^T, after the singular ...
Definition: GslMatrix.h:443
Map * m_svdColMap
Mapping for matrices involved in the singular value decomposition (svd) routine.
Definition: GslMatrix.h:418
double m_lnDeterminant
The natural logarithm of the determinant of this matrix.
Definition: GslMatrix.h:449
int NumGlobalElements() const
Returns the total number of elements across all processors.
Definition: Map.C:85
gsl_matrix * m_mat
GSL matrix, also referred to as this matrix.
Definition: GslMatrix.h:409
GslVector * m_svdSvec
m_svdSvec stores the diagonal of the N-by-N diagonal matrix of singular values S after the singular v...
Definition: GslMatrix.h:431
double m_determinant
The determinant of this matrix.
Definition: GslMatrix.h:446
gsl_matrix * m_LU
GSL matrix for the LU decomposition of m_mat.
Definition: GslMatrix.h:412
const Map & map() const
Definition: Matrix.C:54
GslMatrix * m_inverse
Inverse matrix of this.
Definition: GslMatrix.h:415
Matrix()
Default constructor.
QUESO::GslMatrix::GslMatrix ( const GslVector v,
double  diagValue 
)

Shaped Constructor: creates a square matrix with size v.sizeLocal() and diagonal values all equal to diagValue.

Definition at line 86 of file GslMatrix.C.

References m_mat.

89  :
90  Matrix (v.env(),v.map()),
91  m_mat (gsl_matrix_calloc(v.sizeLocal(),v.sizeLocal())),
92  m_LU (NULL),
93  m_inverse (NULL),
94  m_svdColMap (NULL),
95  m_svdUmat (NULL),
96  m_svdSvec (NULL),
97  m_svdVmat (NULL),
98  m_svdVTmat (NULL),
99  m_determinant (-INFINITY),
100  m_lnDeterminant(-INFINITY),
101  m_permutation (NULL),
102  m_signum (0),
103  m_isSingular (false)
104 {
105  queso_require_msg(m_mat, "null matrix generated");
106 
107  for (unsigned int i = 0; i < m_mat->size1; ++i) {
108  (*this)(i,i) = diagValue;
109  }
110 }
gsl_permutation * m_permutation
The permutation matrix of a LU decomposition.
Definition: GslMatrix.h:455
GslMatrix * m_svdVmat
m_svdVmat stores the N-by-N orthogonal square matrix V after the singular value decomposition of a ma...
Definition: GslMatrix.h:437
int m_signum
m_signum stores the sign of the permutation of the LU decomposition PA = LU.
Definition: GslMatrix.h:461
GslMatrix * m_svdUmat
m_svdUmat stores the M-by-N orthogonal matrix U after the singular value decomposition of a matrix...
Definition: GslMatrix.h:424
bool m_isSingular
Indicates whether or not this matrix is singular.
Definition: GslMatrix.h:464
GslMatrix * m_svdVTmat
m_svdVmatT stores the transpose of N-by-N orthogonal square matrix V, namely V^T, after the singular ...
Definition: GslMatrix.h:443
Map * m_svdColMap
Mapping for matrices involved in the singular value decomposition (svd) routine.
Definition: GslMatrix.h:418
double m_lnDeterminant
The natural logarithm of the determinant of this matrix.
Definition: GslMatrix.h:449
gsl_matrix * m_mat
GSL matrix, also referred to as this matrix.
Definition: GslMatrix.h:409
GslVector * m_svdSvec
m_svdSvec stores the diagonal of the N-by-N diagonal matrix of singular values S after the singular v...
Definition: GslMatrix.h:431
double m_determinant
The determinant of this matrix.
Definition: GslMatrix.h:446
gsl_matrix * m_LU
GSL matrix for the LU decomposition of m_mat.
Definition: GslMatrix.h:412
GslMatrix * m_inverse
Inverse matrix of this.
Definition: GslMatrix.h:415
Matrix()
Default constructor.
QUESO::GslMatrix::GslMatrix ( const GslVector v)

Shaped Constructor: creates a square matrix with size v.sizeLocal().

The diagonal values of this matrix are the elements in vector v.

Definition at line 112 of file GslMatrix.C.

References dim, m_mat, and QUESO::GslVector::sizeLocal().

113  :
114  Matrix (v.env(),v.map()),
115  m_mat (gsl_matrix_calloc(v.sizeLocal(),v.sizeLocal())),
116  m_LU (NULL),
117  m_inverse (NULL),
118  m_svdColMap (NULL),
119  m_svdUmat (NULL),
120  m_svdSvec (NULL),
121  m_svdVmat (NULL),
122  m_svdVTmat (NULL),
123  m_determinant (-INFINITY),
124  m_lnDeterminant(-INFINITY),
125  m_permutation (NULL),
126  m_signum (0),
127  m_isSingular (false)
128 {
129  queso_require_msg(m_mat, "null matrix generated");
130 
131  unsigned int dim = v.sizeLocal();
132  for (unsigned int i = 0; i < dim; ++i) {
133  (*this)(i,i) = v[i];
134  }
135 }
gsl_permutation * m_permutation
The permutation matrix of a LU decomposition.
Definition: GslMatrix.h:455
GslMatrix * m_svdVmat
m_svdVmat stores the N-by-N orthogonal square matrix V after the singular value decomposition of a ma...
Definition: GslMatrix.h:437
int m_signum
m_signum stores the sign of the permutation of the LU decomposition PA = LU.
Definition: GslMatrix.h:461
GslMatrix * m_svdUmat
m_svdUmat stores the M-by-N orthogonal matrix U after the singular value decomposition of a matrix...
Definition: GslMatrix.h:424
bool m_isSingular
Indicates whether or not this matrix is singular.
Definition: GslMatrix.h:464
GslMatrix * m_svdVTmat
m_svdVmatT stores the transpose of N-by-N orthogonal square matrix V, namely V^T, after the singular ...
Definition: GslMatrix.h:443
Map * m_svdColMap
Mapping for matrices involved in the singular value decomposition (svd) routine.
Definition: GslMatrix.h:418
double m_lnDeterminant
The natural logarithm of the determinant of this matrix.
Definition: GslMatrix.h:449
int dim
Definition: ann_test.cpp:472
gsl_matrix * m_mat
GSL matrix, also referred to as this matrix.
Definition: GslMatrix.h:409
GslVector * m_svdSvec
m_svdSvec stores the diagonal of the N-by-N diagonal matrix of singular values S after the singular v...
Definition: GslMatrix.h:431
double m_determinant
The determinant of this matrix.
Definition: GslMatrix.h:446
gsl_matrix * m_LU
GSL matrix for the LU decomposition of m_mat.
Definition: GslMatrix.h:412
GslMatrix * m_inverse
Inverse matrix of this.
Definition: GslMatrix.h:415
Matrix()
Default constructor.
QUESO::GslMatrix::GslMatrix ( const GslMatrix B)

Shaped Constructor: creates a matrix with B.numCols() columns and B.numRowsLocal() rows.

This matrix is a copy of matrix B.

Definition at line 138 of file GslMatrix.C.

References QUESO::Matrix::base_copy(), copy(), and m_mat.

139  :
140  Matrix (B.env(),B.map()),
141  m_mat (gsl_matrix_calloc(B.numRowsLocal(),B.numCols())),
142  m_LU (NULL),
143  m_inverse (NULL),
144  m_svdColMap (NULL),
145  m_svdUmat (NULL),
146  m_svdSvec (NULL),
147  m_svdVmat (NULL),
148  m_svdVTmat (NULL),
149  m_determinant (-INFINITY),
150  m_lnDeterminant(-INFINITY),
151  m_permutation (NULL),
152  m_signum (0),
153  m_isSingular (false)
154 {
155  queso_require_msg(m_mat, "null vector generated");
156  this->Matrix::base_copy(B);
157  this->copy(B);
158 }
gsl_permutation * m_permutation
The permutation matrix of a LU decomposition.
Definition: GslMatrix.h:455
GslMatrix * m_svdVmat
m_svdVmat stores the N-by-N orthogonal square matrix V after the singular value decomposition of a ma...
Definition: GslMatrix.h:437
int m_signum
m_signum stores the sign of the permutation of the LU decomposition PA = LU.
Definition: GslMatrix.h:461
GslMatrix * m_svdUmat
m_svdUmat stores the M-by-N orthogonal matrix U after the singular value decomposition of a matrix...
Definition: GslMatrix.h:424
bool m_isSingular
Indicates whether or not this matrix is singular.
Definition: GslMatrix.h:464
GslMatrix * m_svdVTmat
m_svdVmatT stores the transpose of N-by-N orthogonal square matrix V, namely V^T, after the singular ...
Definition: GslMatrix.h:443
Map * m_svdColMap
Mapping for matrices involved in the singular value decomposition (svd) routine.
Definition: GslMatrix.h:418
double m_lnDeterminant
The natural logarithm of the determinant of this matrix.
Definition: GslMatrix.h:449
void copy(const GslMatrix &src)
In this function this matrix receives a copy of matrix src.
Definition: GslMatrix.C:220
gsl_matrix * m_mat
GSL matrix, also referred to as this matrix.
Definition: GslMatrix.h:409
GslVector * m_svdSvec
m_svdSvec stores the diagonal of the N-by-N diagonal matrix of singular values S after the singular v...
Definition: GslMatrix.h:431
virtual void base_copy(const Matrix &src)
Copies base data from matrix src to this matrix.
Definition: Matrix.C:99
double m_determinant
The determinant of this matrix.
Definition: GslMatrix.h:446
gsl_matrix * m_LU
GSL matrix for the LU decomposition of m_mat.
Definition: GslMatrix.h:412
GslMatrix * m_inverse
Inverse matrix of this.
Definition: GslMatrix.h:415
Matrix()
Default constructor.
QUESO::GslMatrix::~GslMatrix ( )

Destructor.

Definition at line 161 of file GslMatrix.C.

References m_mat, and resetLU().

162 {
163  this->resetLU();
164  if (m_mat) gsl_matrix_free(m_mat);
165 }
void resetLU()
In this function resets the LU decomposition of this matrix, as well as deletes the private member po...
Definition: GslMatrix.C:232
gsl_matrix * m_mat
GSL matrix, also referred to as this matrix.
Definition: GslMatrix.h:409
QUESO::GslMatrix::GslMatrix ( )
private

Default Constructor.

Creates an empty matrix vector of no dimension. It should not be used by user.

Referenced by internalSvd(), and inverse().

Member Function Documentation

int QUESO::GslMatrix::chol ( )
virtual

Computes Cholesky factorization of a real symmetric positive definite matrix this.

In case fails to be symmetric and positive definite, an error will be returned.

Implements QUESO::Matrix.

Definition at line 408 of file GslMatrix.C.

References QUESO::Matrix::m_env, m_mat, QUESO::UQ_MATRIX_IS_NOT_POS_DEFINITE_RC, and QUESO::BaseEnvironment::worldRank().

409 {
410  int iRC;
411  //std::cout << "Calling gsl_linalg_cholesky_decomp()..." << std::endl;
412  gsl_error_handler_t* oldHandler;
413  oldHandler = gsl_set_error_handler_off();
414  iRC = gsl_linalg_cholesky_decomp(m_mat);
415  if (iRC != 0) {
416  std::cerr << "In GslMatrix::chol()"
417  << ": iRC = " << iRC
418  << ", gsl error message = " << gsl_strerror(iRC)
419  << std::endl;
420  std::cerr << "Here is the offending matrix: " << std::endl;
421  std::cerr << *this << std::endl;
422  }
423  gsl_set_error_handler(oldHandler);
424  //std::cout << "Returned from gsl_linalg_cholesky_decomp() with iRC = " << iRC << std::endl;
425  UQ_RC_MACRO(iRC, // Yes, *not* a fatal check on RC
426  m_env.worldRank(),
427  "GslMatrix::chol()",
428  "matrix is not positive definite",
430 
431  return iRC;
432 }
const int UQ_MATRIX_IS_NOT_POS_DEFINITE_RC
Definition: Defines.h:100
int worldRank() const
Returns the same thing as fullRank()
Definition: Environment.C:262
gsl_matrix * m_mat
GSL matrix, also referred to as this matrix.
Definition: GslMatrix.h:409
const BaseEnvironment & m_env
QUESO environment variable.
Definition: Matrix.h:116
void QUESO::GslMatrix::copy ( const GslMatrix src)
private

In this function this matrix receives a copy of matrix src.

Definition at line 220 of file GslMatrix.C.

References m_mat, and resetLU().

Referenced by GslMatrix(), and operator=().

221 {
222  // FIXME - should we be calling Matrix::base_copy here? - RHS
223  this->resetLU();
224  int iRC;
225  iRC = gsl_matrix_memcpy(this->m_mat, src.m_mat);
226  queso_require_msg(!(iRC), "failed");
227 
228  return;
229 }
void resetLU()
In this function resets the LU decomposition of this matrix, as well as deletes the private member po...
Definition: GslMatrix.C:232
gsl_matrix * m_mat
GSL matrix, also referred to as this matrix.
Definition: GslMatrix.h:409
void QUESO::GslMatrix::cwExtract ( unsigned int  rowId,
unsigned int  colId,
GslMatrix mat 
) const

Definition at line 385 of file GslMatrix.C.

References numCols(), and numRowsLocal().

389 {
390  queso_require_less_msg(initialTargetRowId, this->numRowsLocal(), "invalid initialTargetRowId");
391 
392  queso_require_less_equal_msg((initialTargetRowId + mat.numRowsLocal()), this->numRowsLocal(), "invalid vec.numRowsLocal()");
393 
394  queso_require_less_msg(initialTargetColId, this->numCols(), "invalid initialTargetColId");
395 
396  queso_require_less_equal_msg((initialTargetColId + mat.numCols()), this->numCols(), "invalid vec.numCols()");
397 
398  for (unsigned int i = 0; i < mat.numRowsLocal(); ++i) {
399  for (unsigned int j = 0; j < mat.numCols(); ++j) {
400  mat(i,j) = (*this)(initialTargetRowId+i,initialTargetColId+j) ;
401  }
402  }
403 
404  return;
405 }
unsigned int numRowsLocal() const
Returns the local row dimension of this matrix.
Definition: GslMatrix.C:275
unsigned int numCols() const
Returns the column dimension of this matrix.
Definition: GslMatrix.C:287
void QUESO::GslMatrix::cwSet ( double  value)

Component-wise set all values to this with value.

Definition at line 347 of file GslMatrix.C.

References m_mat, numCols(), and numRowsLocal().

348 {
349  unsigned int nRows = this->numRowsLocal();
350  unsigned int nCols = this->numCols();
351 
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;
355  }
356  }
357 
358  return;
359 }
unsigned int numRowsLocal() const
Returns the local row dimension of this matrix.
Definition: GslMatrix.C:275
gsl_matrix * m_mat
GSL matrix, also referred to as this matrix.
Definition: GslMatrix.h:409
unsigned int numCols() const
Returns the column dimension of this matrix.
Definition: GslMatrix.C:287
void QUESO::GslMatrix::cwSet ( unsigned int  rowId,
unsigned int  colId,
const GslMatrix mat 
)

Set the components of which positions are greater than (rowId,colId) with the value of mat(rowId,colId).

Definition at line 362 of file GslMatrix.C.

References numCols(), and numRowsLocal().

366 {
367  queso_require_less_msg(initialTargetRowId, this->numRowsLocal(), "invalid initialTargetRowId");
368 
369  queso_require_less_equal_msg((initialTargetRowId + mat.numRowsLocal()), this->numRowsLocal(), "invalid vec.numRowsLocal()");
370 
371  queso_require_less_msg(initialTargetColId, this->numCols(), "invalid initialTargetColId");
372 
373  queso_require_less_equal_msg((initialTargetColId + mat.numCols()), this->numCols(), "invalid vec.numCols()");
374 
375  for (unsigned int i = 0; i < mat.numRowsLocal(); ++i) {
376  for (unsigned int j = 0; j < mat.numCols(); ++j) {
377  (*this)(initialTargetRowId+i,initialTargetColId+j) = mat(i,j);
378  }
379  }
380 
381  return;
382 }
unsigned int numRowsLocal() const
Returns the local row dimension of this matrix.
Definition: GslMatrix.C:275
unsigned int numCols() const
Returns the column dimension of this matrix.
Definition: GslMatrix.C:287
gsl_matrix * QUESO::GslMatrix::data ( )

Returns this matrix.

Definition at line 1766 of file GslMatrix.C.

References m_mat.

Referenced by internalSvd(), and svdSolve().

1767 {
1768  return m_mat;
1769 }
gsl_matrix * m_mat
GSL matrix, also referred to as this matrix.
Definition: GslMatrix.h:409
double QUESO::GslMatrix::determinant ( ) const

Calculates the determinant of this matrix.

Definition at line 1038 of file GslMatrix.C.

References QUESO::BaseEnvironment::displayVerbosity(), invertMultiply(), m_determinant, QUESO::Matrix::m_env, m_lnDeterminant, m_LU, QUESO::Matrix::m_map, m_signum, and QUESO::BaseEnvironment::subDisplayFile().

1039 {
1040  if (m_determinant == -INFINITY) {
1041  if (m_LU == NULL) {
1042  GslVector tmpB(m_env,m_map);
1043  GslVector tmpX(m_env,m_map);
1044  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 99)) {
1045  *m_env.subDisplayFile() << "In GslMatrix::determinant()"
1046  << ": before 'this->invertMultiply()'"
1047  << std::endl;
1048  }
1049  this->invertMultiply(tmpB,tmpX);
1050  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 99)) {
1051  *m_env.subDisplayFile() << "In GslMatrix::determinant()"
1052  << ": after 'this->invertMultiply()'"
1053  << std::endl;
1054  }
1055  }
1056  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 99)) {
1057  *m_env.subDisplayFile() << "In GslMatrix::determinant()"
1058  << ": before 'gsl_linalg_LU_det()'"
1059  << std::endl;
1060  }
1061  m_determinant = gsl_linalg_LU_det(m_LU,m_signum);
1062  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 99)) {
1063  *m_env.subDisplayFile() << "In GslMatrix::determinant()"
1064  << ": after 'gsl_linalg_LU_det()'"
1065  << std::endl;
1066  }
1067  m_lnDeterminant = gsl_linalg_LU_lndet(m_LU);
1068  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 99)) {
1069  *m_env.subDisplayFile() << "In GslMatrix::determinant()"
1070  << ": after 'gsl_linalg_LU_lndet()'"
1071  << std::endl;
1072  }
1073  }
1074 
1075  return m_determinant;
1076 }
GslVector invertMultiply(const GslVector &b) const
This function calculates the inverse of this matrix and multiplies it with vector b...
Definition: GslMatrix.C:1222
const Map m_map
Mapping variable.
Definition: Matrix.h:121
int m_signum
m_signum stores the sign of the permutation of the LU decomposition PA = LU.
Definition: GslMatrix.h:461
double m_lnDeterminant
The natural logarithm of the determinant of this matrix.
Definition: GslMatrix.h:449
double m_determinant
The determinant of this matrix.
Definition: GslMatrix.h:446
unsigned int displayVerbosity() const
Definition: Environment.C:450
const BaseEnvironment & m_env
QUESO environment variable.
Definition: Matrix.h:116
gsl_matrix * m_LU
GSL matrix for the LU decomposition of m_mat.
Definition: GslMatrix.h:412
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:320
void QUESO::GslMatrix::eigen ( GslVector eigenValues,
GslMatrix eigenVectors 
) const

This function computes the eigenvalues of a real symmetric matrix.

Definition at line 1420 of file GslMatrix.C.

References QUESO::GslVector::data(), m_mat, numRowsLocal(), QUESO::queso_require_equal_to_msg, QUESO::queso_require_not_equal_to_msg, and QUESO::GslVector::sizeLocal().

1421 {
1422  unsigned int n = eigenValues.sizeLocal();
1423 
1424  queso_require_not_equal_to_msg(n, 0, "invalid input vector size");
1425 
1426  if (eigenVectors) {
1427  queso_require_equal_to_msg(eigenValues.sizeLocal(), eigenVectors->numRowsLocal(), "different input vector sizes");
1428  }
1429 
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);
1434  }
1435  else {
1436  gsl_eigen_symmv_workspace* w = gsl_eigen_symmv_alloc((size_t) n);
1437  gsl_eigen_symmv(m_mat,eigenValues.data(),eigenVectors->m_mat,w);
1438  gsl_eigen_symmv_sort(eigenValues.data(),eigenVectors->m_mat,GSL_EIGEN_SORT_VAL_ASC);
1439  gsl_eigen_symmv_free(w);
1440  }
1441 
1442  return;
1443 }
MonteCarloSGOptions::MonteCarloSGOptions(const BaseEnvironment &env, const char *prefix queso_require_not_equal_to_msg)(m_env.optionsInputFileName(), std::string(""), std::string("this constructor is incompatible with the absence of an options input file"))
MonteCarloSGOptions::MonteCarloSGOptions(const BaseEnvironment &env, const char *prefix, const McOptionsValues &alternativeOptionsValues queso_require_equal_to_msg)(m_env.optionsInputFileName(), std::string(""), std::string("this constructor is incompatible with the existence of an options input file"))
gsl_matrix * m_mat
GSL matrix, also referred to as this matrix.
Definition: GslMatrix.h:409
void QUESO::GslMatrix::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.

Definition at line 753 of file GslMatrix.C.

References numCols(), numRowsLocal(), and QUESO::queso_require_equal_to_msg.

759 {
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();
765  }
766  queso_require_greater_equal_msg(this->numRowsLocal(), (initialTargetRowId + sumNumRowsLocals), "too big number of rows");
767  if (checkForExactNumRowsMatching) queso_require_equal_to_msg(this->numRowsLocal(), (initialTargetRowId + sumNumRowsLocals), "inconsistent number of rows");
768  queso_require_greater_equal_msg(this->numCols(), (initialTargetColId + sumNumCols), "too big number of cols");
769  if (checkForExactNumColsMatching) queso_require_equal_to_msg(this->numCols(), (initialTargetColId + sumNumCols), "inconsistent number of cols");
770 
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);
779  }
780  }
781  cumulativeRowId += nRows;
782  cumulativeColId += nCols;
783  }
784 
785  return;
786 }
unsigned int numRowsLocal() const
Returns the local row dimension of this matrix.
Definition: GslMatrix.C:275
MonteCarloSGOptions::MonteCarloSGOptions(const BaseEnvironment &env, const char *prefix, const McOptionsValues &alternativeOptionsValues queso_require_equal_to_msg)(m_env.optionsInputFileName(), std::string(""), std::string("this constructor is incompatible with the existence of an options input file"))
unsigned int numCols() const
Returns the column dimension of this matrix.
Definition: GslMatrix.C:287
void QUESO::GslMatrix::fillWithBlocksDiagonally ( unsigned int  rowId,
unsigned int  colId,
const std::vector< GslMatrix * > &  matrices,
bool  checkForExactNumRowsMatching,
bool  checkForExactNumColsMatching 
)

This function fills this matrix diagonally with block matrices.

Definition at line 789 of file GslMatrix.C.

References numCols(), numRowsLocal(), and QUESO::queso_require_equal_to_msg.

795 {
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();
801  }
802  queso_require_greater_equal_msg(this->numRowsLocal(), (initialTargetRowId + sumNumRowsLocals), "too big number of rows");
803  if (checkForExactNumRowsMatching) queso_require_equal_to_msg(this->numRowsLocal(), (initialTargetRowId + sumNumRowsLocals), "inconsistent number of rows");
804  queso_require_greater_equal_msg(this->numCols(), (initialTargetColId + sumNumCols), "too big number of cols");
805  if (checkForExactNumColsMatching) queso_require_equal_to_msg(this->numCols(), (initialTargetColId + sumNumCols), "inconsistent number of cols");
806 
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);
815  }
816  }
817  cumulativeRowId += nRows;
818  cumulativeColId += nCols;
819  }
820 
821  return;
822 }
unsigned int numRowsLocal() const
Returns the local row dimension of this matrix.
Definition: GslMatrix.C:275
MonteCarloSGOptions::MonteCarloSGOptions(const BaseEnvironment &env, const char *prefix, const McOptionsValues &alternativeOptionsValues queso_require_equal_to_msg)(m_env.optionsInputFileName(), std::string(""), std::string("this constructor is incompatible with the existence of an options input file"))
unsigned int numCols() const
Returns the column dimension of this matrix.
Definition: GslMatrix.C:287
void QUESO::GslMatrix::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.

Definition at line 825 of file GslMatrix.C.

References numCols(), numRowsLocal(), and QUESO::queso_require_equal_to_msg.

831 {
832  unsigned int sumNumCols = 0;
833  for (unsigned int i = 0; i < matrices.size(); ++i) {
834  queso_require_greater_equal_msg(this->numRowsLocal(), (initialTargetRowId + matrices[i]->numRowsLocal()), "too big number of rows");
835  if (checkForExactNumRowsMatching) queso_require_equal_to_msg(this->numRowsLocal(), (initialTargetRowId + matrices[i]->numRowsLocal()), "inconsistent number of rows");
836  sumNumCols += matrices[i]->numCols();
837  }
838  queso_require_greater_equal_msg(this->numCols(), (initialTargetColId + sumNumCols), "too big number of cols");
839  if (checkForExactNumColsMatching) queso_require_equal_to_msg(this->numCols(), (initialTargetColId + sumNumCols), "inconsistent number of cols");
840 
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);
848  }
849  }
850  cumulativeColId += nCols;
851  }
852 
853  return;
854 }
unsigned int numRowsLocal() const
Returns the local row dimension of this matrix.
Definition: GslMatrix.C:275
MonteCarloSGOptions::MonteCarloSGOptions(const BaseEnvironment &env, const char *prefix, const McOptionsValues &alternativeOptionsValues queso_require_equal_to_msg)(m_env.optionsInputFileName(), std::string(""), std::string("this constructor is incompatible with the existence of an options input file"))
unsigned int numCols() const
Returns the column dimension of this matrix.
Definition: GslMatrix.C:287
void QUESO::GslMatrix::fillWithBlocksHorizontally ( unsigned int  rowId,
unsigned int  colId,
const std::vector< GslMatrix * > &  matrices,
bool  checkForExactNumRowsMatching,
bool  checkForExactNumColsMatching 
)

This function fills this matrix horizontally with const block.

Definition at line 857 of file GslMatrix.C.

References numCols(), numRowsLocal(), and QUESO::queso_require_equal_to_msg.

863 {
864  unsigned int sumNumCols = 0;
865  for (unsigned int i = 0; i < matrices.size(); ++i) {
866  queso_require_greater_equal_msg(this->numRowsLocal(), (initialTargetRowId + matrices[i]->numRowsLocal()), "too big number of rows");
867  if (checkForExactNumRowsMatching) queso_require_equal_to_msg(this->numRowsLocal(), (initialTargetRowId + matrices[i]->numRowsLocal()), "inconsistent number of rows");
868  sumNumCols += matrices[i]->numCols();
869  }
870  queso_require_greater_equal_msg(this->numCols(), (initialTargetColId + sumNumCols), "too big number of cols");
871  if (checkForExactNumColsMatching) queso_require_equal_to_msg(this->numCols(), (initialTargetColId + sumNumCols), "inconsistent number of cols");
872 
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);
880  }
881  }
882  cumulativeColId += nCols;
883  }
884 
885  return;
886 }
unsigned int numRowsLocal() const
Returns the local row dimension of this matrix.
Definition: GslMatrix.C:275
MonteCarloSGOptions::MonteCarloSGOptions(const BaseEnvironment &env, const char *prefix, const McOptionsValues &alternativeOptionsValues queso_require_equal_to_msg)(m_env.optionsInputFileName(), std::string(""), std::string("this constructor is incompatible with the existence of an options input file"))
unsigned int numCols() const
Returns the column dimension of this matrix.
Definition: GslMatrix.C:287
void QUESO::GslMatrix::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.

Definition at line 889 of file GslMatrix.C.

References numCols(), numRowsLocal(), and QUESO::queso_require_equal_to_msg.

895 {
896  unsigned int sumNumRows = 0;
897  for (unsigned int i = 0; i < matrices.size(); ++i) {
898  queso_require_greater_equal_msg(this->numCols(), (initialTargetColId + matrices[i]->numCols()), "too big number of cols");
899  if (checkForExactNumColsMatching) queso_require_equal_to_msg(this->numCols(), (initialTargetColId + matrices[i]->numCols()), "inconsistent number of cols");
900  sumNumRows += matrices[i]->numRowsLocal();
901  }
902  queso_require_greater_equal_msg(this->numRowsLocal(), (initialTargetRowId + sumNumRows), "too big number of rows");
903  if (checkForExactNumRowsMatching) queso_require_equal_to_msg(this->numRowsLocal(), (initialTargetRowId + sumNumRows), "inconsistent number of rows");
904 
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);
912  }
913  }
914  cumulativeRowId += nRows;
915  }
916 
917  return;
918 }
unsigned int numRowsLocal() const
Returns the local row dimension of this matrix.
Definition: GslMatrix.C:275
MonteCarloSGOptions::MonteCarloSGOptions(const BaseEnvironment &env, const char *prefix, const McOptionsValues &alternativeOptionsValues queso_require_equal_to_msg)(m_env.optionsInputFileName(), std::string(""), std::string("this constructor is incompatible with the existence of an options input file"))
unsigned int numCols() const
Returns the column dimension of this matrix.
Definition: GslMatrix.C:287
void QUESO::GslMatrix::fillWithBlocksVertically ( unsigned int  rowId,
unsigned int  colId,
const std::vector< GslMatrix * > &  matrices,
bool  checkForExactNumRowsMatching,
bool  checkForExactNumColsMatching 
)

This function fills this matrix vertically with block matrices.

Definition at line 921 of file GslMatrix.C.

References numCols(), numRowsLocal(), and QUESO::queso_require_equal_to_msg.

927 {
928  unsigned int sumNumRows = 0;
929  for (unsigned int i = 0; i < matrices.size(); ++i) {
930  queso_require_greater_equal_msg(this->numCols(), (initialTargetColId + matrices[i]->numCols()), "too big number of cols");
931  if (checkForExactNumColsMatching) queso_require_equal_to_msg(this->numCols(), (initialTargetColId + matrices[i]->numCols()), "inconsistent number of cols");
932  sumNumRows += matrices[i]->numRowsLocal();
933  }
934  queso_require_greater_equal_msg(this->numRowsLocal(), (initialTargetRowId + sumNumRows), "too big number of rows");
935  if (checkForExactNumRowsMatching) queso_require_equal_to_msg(this->numRowsLocal(), (initialTargetRowId + sumNumRows), "inconsistent number of rows");
936 
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);
944  }
945  }
946  cumulativeRowId += nRows;
947  }
948 
949  return;
950 }
unsigned int numRowsLocal() const
Returns the local row dimension of this matrix.
Definition: GslMatrix.C:275
MonteCarloSGOptions::MonteCarloSGOptions(const BaseEnvironment &env, const char *prefix, const McOptionsValues &alternativeOptionsValues queso_require_equal_to_msg)(m_env.optionsInputFileName(), std::string(""), std::string("this constructor is incompatible with the existence of an options input file"))
unsigned int numCols() const
Returns the column dimension of this matrix.
Definition: GslMatrix.C:287
void QUESO::GslMatrix::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.

Definition at line 953 of file GslMatrix.C.

References numCols(), numRowsLocal(), and QUESO::queso_require_equal_to_msg.

960 {
961  queso_require_greater_equal_msg(this->numRowsLocal(), (initialTargetRowId + (mat1.numRowsLocal() * mat2.numRowsLocal())), "too big number of rows");
962  if (checkForExactNumRowsMatching) queso_require_equal_to_msg(this->numRowsLocal(), (initialTargetRowId + (mat1.numRowsLocal() * mat2.numRowsLocal())), "inconsistent number of rows");
963  queso_require_greater_equal_msg(this->numCols(), (initialTargetColId + (mat1.numCols() * mat2.numCols())), "too big number of columns");
964  if (checkForExactNumColsMatching) queso_require_equal_to_msg(this->numCols(), (initialTargetColId + (mat1.numCols() * mat2.numCols())), "inconsistent number of columns");
965 
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);
974  }
975  }
976  }
977  }
978 
979  return;
980 }
unsigned int numRowsLocal() const
Returns the local row dimension of this matrix.
Definition: GslMatrix.C:275
MonteCarloSGOptions::MonteCarloSGOptions(const BaseEnvironment &env, const char *prefix, const McOptionsValues &alternativeOptionsValues queso_require_equal_to_msg)(m_env.optionsInputFileName(), std::string(""), std::string("this constructor is incompatible with the existence of an options input file"))
unsigned int numCols() const
Returns the column dimension of this matrix.
Definition: GslMatrix.C:287
void QUESO::GslMatrix::fillWithTensorProduct ( unsigned int  rowId,
unsigned int  colId,
const GslMatrix mat1,
const GslVector vec2,
bool  checkForExactNumRowsMatching,
bool  checkForExactNumColsMatching 
)

This function calculates the tensor product of matrix mat1 and vector vec2 and stores it in this matrix.

Definition at line 983 of file GslMatrix.C.

References numCols(), numRowsLocal(), QUESO::queso_require_equal_to_msg, and QUESO::GslVector::sizeLocal().

990 {
991  queso_require_greater_equal_msg(this->numRowsLocal(), (initialTargetRowId + (mat1.numRowsLocal() * vec2.sizeLocal())), "too big number of rows");
992  if (checkForExactNumRowsMatching) queso_require_equal_to_msg(this->numRowsLocal(), (initialTargetRowId + (mat1.numRowsLocal() * vec2.sizeLocal())), "inconsistent number of rows");
993  queso_require_greater_equal_msg(this->numCols(), (initialTargetColId + (mat1.numCols() * 1)), "too big number of columns");
994  if (checkForExactNumColsMatching) queso_require_equal_to_msg(this->numCols(), (initialTargetColId + (mat1.numCols() * 1)), "inconsistent number of columns");
995 
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];
1004  }
1005  }
1006  }
1007  }
1008 
1009 
1010  return;
1011 }
unsigned int numRowsLocal() const
Returns the local row dimension of this matrix.
Definition: GslMatrix.C:275
MonteCarloSGOptions::MonteCarloSGOptions(const BaseEnvironment &env, const char *prefix, const McOptionsValues &alternativeOptionsValues queso_require_equal_to_msg)(m_env.optionsInputFileName(), std::string(""), std::string("this constructor is incompatible with the existence of an options input file"))
unsigned int numCols() const
Returns the column dimension of this matrix.
Definition: GslMatrix.C:287
void QUESO::GslMatrix::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.

Definition at line 1014 of file GslMatrix.C.

References numCols(), numRowsLocal(), and QUESO::queso_require_equal_to_msg.

1020 {
1021  unsigned int nRows = mat.numRowsLocal();
1022  unsigned int nCols = mat.numCols();
1023  queso_require_greater_equal_msg(this->numRowsLocal(), (initialTargetRowId + nCols), "too big number of rows");
1024  if (checkForExactNumRowsMatching) queso_require_equal_to_msg(this->numRowsLocal(), (initialTargetRowId + nCols), "inconsistent number of rows");
1025  queso_require_greater_equal_msg(this->numCols(), (initialTargetColId + nRows), "too big number of cols");
1026  if (checkForExactNumColsMatching) queso_require_equal_to_msg(this->numCols(), (initialTargetColId + nRows), "inconsistent number of cols");
1027 
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);
1031  }
1032  }
1033 
1034  return;
1035 }
unsigned int numRowsLocal() const
Returns the local row dimension of this matrix.
Definition: GslMatrix.C:275
MonteCarloSGOptions::MonteCarloSGOptions(const BaseEnvironment &env, const char *prefix, const McOptionsValues &alternativeOptionsValues queso_require_equal_to_msg)(m_env.optionsInputFileName(), std::string(""), std::string("this constructor is incompatible with the existence of an options input file"))
unsigned int numCols() const
Returns the column dimension of this matrix.
Definition: GslMatrix.C:287
void QUESO::GslMatrix::filterLargeValues ( double  thresholdValue)

This function sets to zero (filters) all entries of this matrix which are greater than thresholdValue.

If thresholdValue < 0 then no values will be filtered.

Definition at line 670 of file GslMatrix.C.

References numCols(), and numRowsLocal().

671 {
672  unsigned int nRows = this->numRowsLocal();
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);
677  // If 'thresholdValue' is negative, no values will be filtered
678  if ((aux < 0. ) &&
679  (-thresholdValue > aux)) {
680  (*this)(i,j) = 0.;
681  }
682  if ((aux > 0. ) &&
683  (thresholdValue < aux)) {
684  (*this)(i,j) = 0.;
685  }
686  }
687  }
688 
689  return;
690 }
unsigned int numRowsLocal() const
Returns the local row dimension of this matrix.
Definition: GslMatrix.C:275
unsigned int numCols() const
Returns the column dimension of this matrix.
Definition: GslMatrix.C:287
void QUESO::GslMatrix::filterSmallValues ( double  thresholdValue)

This function sets to zero (filters) all entries of this matrix which are smaller than thresholdValue.

If thresholdValue < 0 then no values will be filtered.

Definition at line 647 of file GslMatrix.C.

References numCols(), and numRowsLocal().

648 {
649  unsigned int nRows = this->numRowsLocal();
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);
654  // If 'thresholdValue' is negative, no values will be filtered
655  if ((aux < 0. ) &&
656  (-thresholdValue < aux)) {
657  (*this)(i,j) = 0.;
658  }
659  if ((aux > 0. ) &&
660  (thresholdValue > aux)) {
661  (*this)(i,j) = 0.;
662  }
663  }
664  }
665 
666  return;
667 }
unsigned int numRowsLocal() const
Returns the local row dimension of this matrix.
Definition: GslMatrix.C:275
unsigned int numCols() const
Returns the column dimension of this matrix.
Definition: GslMatrix.C:287
void QUESO::GslMatrix::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.

Definition at line 1586 of file GslMatrix.C.

References m_mat, numCols(), numRowsLocal(), QUESO::queso_require_equal_to_msg, and QUESO::GslVector::sizeLocal().

Referenced by getColumn(), invertMultiply(), matlabLinearInterpExtrap(), and svdSolve().

1587 {
1588  // Sanity checks
1589  queso_require_less_msg(column_num, this->numCols(), "Specified row number not within range");
1590 
1591  queso_require_equal_to_msg(column.sizeLocal(), this->numRowsLocal(), "column vector not same size as this matrix");
1592 
1593  // Temporary working vector
1594  gsl_vector* gsl_column = gsl_vector_alloc( column.sizeLocal() );
1595 
1596  int error_code = gsl_matrix_get_col( gsl_column, m_mat, column_num );
1597  queso_require_equal_to_msg(error_code, 0, "gsl_matrix_get_col failed");
1598 
1599  // Copy column from gsl matrix into our GslVector object
1600  for( unsigned int i = 0; i < column.sizeLocal(); ++i )
1601  {
1602  column[i] = gsl_vector_get( gsl_column, i );
1603  }
1604 
1605  // Clean up gsl temporaries
1606  gsl_vector_free( gsl_column );
1607 
1608  return;
1609 }
unsigned int numRowsLocal() const
Returns the local row dimension of this matrix.
Definition: GslMatrix.C:275
MonteCarloSGOptions::MonteCarloSGOptions(const BaseEnvironment &env, const char *prefix, const McOptionsValues &alternativeOptionsValues queso_require_equal_to_msg)(m_env.optionsInputFileName(), std::string(""), std::string("this constructor is incompatible with the existence of an options input file"))
gsl_matrix * m_mat
GSL matrix, also referred to as this matrix.
Definition: GslMatrix.h:409
unsigned int numCols() const
Returns the column dimension of this matrix.
Definition: GslMatrix.C:287
GslVector QUESO::GslMatrix::getColumn ( const unsigned int  column_num) const

This function gets the column_num-th column of this matrix.

Definition at line 1651 of file GslMatrix.C.

References getColumn(), QUESO::Matrix::m_env, and QUESO::Matrix::m_map.

1652 {
1653  // We rely on the void getRow's to do sanity checks
1654  // So we don't do them here.
1655 
1656  GslVector column(m_env, m_map);
1657 
1658  this->getColumn( column_num, column );
1659 
1660  return column;
1661 }
const Map m_map
Mapping variable.
Definition: Matrix.h:121
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...
Definition: GslMatrix.C:1586
const BaseEnvironment & m_env
QUESO environment variable.
Definition: Matrix.h:116
void QUESO::GslMatrix::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.

Definition at line 1612 of file GslMatrix.C.

References m_mat, numCols(), numRowsLocal(), QUESO::queso_require_equal_to_msg, and QUESO::GslVector::sizeLocal().

Referenced by getRow().

1613 {
1614  // Sanity checks
1615  queso_require_less_msg(row_num, this->numRowsLocal(), "Specified row number not within range");
1616 
1617  queso_require_equal_to_msg(row.sizeLocal(), this->numCols(), "row vector not same size as this matrix");
1618 
1619  // Temporary working vector
1620  gsl_vector* gsl_row = gsl_vector_alloc( row.sizeLocal() );
1621 
1622  int error_code = gsl_matrix_get_row( gsl_row, m_mat, row_num );
1623  queso_require_equal_to_msg(error_code, 0, "gsl_matrix_get_row failed");
1624 
1625  // Copy row from gsl matrix into our GslVector object
1626  for( unsigned int i = 0; i < row.sizeLocal(); ++i )
1627  {
1628  row[i] = gsl_vector_get( gsl_row, i );
1629  }
1630 
1631  // Clean up gsl temporaries
1632  gsl_vector_free( gsl_row );
1633 
1634  return;
1635 }
unsigned int numRowsLocal() const
Returns the local row dimension of this matrix.
Definition: GslMatrix.C:275
MonteCarloSGOptions::MonteCarloSGOptions(const BaseEnvironment &env, const char *prefix, const McOptionsValues &alternativeOptionsValues queso_require_equal_to_msg)(m_env.optionsInputFileName(), std::string(""), std::string("this constructor is incompatible with the existence of an options input file"))
gsl_matrix * m_mat
GSL matrix, also referred to as this matrix.
Definition: GslMatrix.h:409
unsigned int numCols() const
Returns the column dimension of this matrix.
Definition: GslMatrix.C:287
GslVector QUESO::GslMatrix::getRow ( const unsigned int  row_num) const

This function gets the row_num-th column of this matrix.

Definition at line 1638 of file GslMatrix.C.

References getRow(), QUESO::Matrix::m_env, and QUESO::Matrix::m_map.

1639 {
1640  // We rely on the void getRow's to do sanity checks
1641  // So we don't do them here.
1642 
1643  GslVector row(m_env, m_map);
1644 
1645  this->getRow( row_num, row );
1646 
1647  return row;
1648 }
const Map m_map
Mapping variable.
Definition: Matrix.h:121
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...
Definition: GslMatrix.C:1612
const BaseEnvironment & m_env
QUESO environment variable.
Definition: Matrix.h:116
int QUESO::GslMatrix::internalSvd ( ) const
private

This function factorizes the M-by-N matrix A into the singular value decomposition A = U S V^T for M >= N. On output the matrix A is replaced by U.

Definition at line 532 of file GslMatrix.C.

References QUESO::GslVector::data(), data(), GslMatrix(), QUESO::Matrix::m_env, m_svdColMap, m_svdSvec, m_svdUmat, m_svdVmat, m_svdVTmat, QUESO::Matrix::map(), numCols(), numRowsLocal(), transpose(), QUESO::UQ_MATRIX_SVD_FAILED_RC, and QUESO::BaseEnvironment::worldRank().

Referenced by rank(), svd(), svdMatU(), svdMatV(), and svdSolve().

533 {
534  int iRC = 0;
535 
536  if (m_svdColMap == NULL) {
537  unsigned int nRows = this->numRowsLocal();
538  unsigned int nCols = this->numCols();
539  queso_require_greater_equal_msg(nRows, nCols, "GSL only supports cases where nRows >= nCols");
540 
541  m_svdColMap = new Map(this->numCols(),0,this->map().Comm()); // see 'VectorSpace<.,.>::newMap()' in src/basic/src/GslVectorSpace.C
542  m_svdUmat = new GslMatrix(*this); // Yes, 'this'
543  m_svdSvec = new GslVector(m_env,*m_svdColMap);
544  m_svdVmat = new GslMatrix(*m_svdSvec);
546 
547  //std::cout << "In GslMatrix::internalSvd()"
548  // << ", calling gsl_linalg_SV_decomp_jacobi()..."
549  // << ": nRows = " << nRows
550  // << ", nCols = " << nCols
551  // << std::endl;
552  struct timeval timevalBegin;
553  gettimeofday(&timevalBegin, NULL);
554  gsl_error_handler_t* oldHandler;
555  oldHandler = gsl_set_error_handler_off();
556 #if 1
557  iRC = gsl_linalg_SV_decomp_jacobi(m_svdUmat->data(), m_svdVmat->data(), m_svdSvec->data());
558 #else
559  GslVector vecWork(*m_svdSvec );
560  iRC = gsl_linalg_SV_decomp(m_svdUmat->data(), m_svdVmat->data(), m_svdSvec->data(), vecWork.data());
561 #endif
562  if (iRC != 0) {
563  std::cerr << "In GslMatrix::internalSvd()"
564  << ": iRC = " << iRC
565  << ", gsl error message = " << gsl_strerror(iRC)
566  << std::endl;
567  }
568  gsl_set_error_handler(oldHandler);
569 
570  struct timeval timevalNow;
571  gettimeofday(&timevalNow, NULL);
572  //std::cout << "In GslMatrix::internalSvd()"
573  // << ": returned from gsl_linalg_SV_decomp_jacobi() with iRC = " << iRC
574  // << " after " << timevalNow.tv_sec - timevalBegin.tv_sec
575  // << " seconds"
576  // << std::endl;
577  UQ_RC_MACRO(iRC, // Yes, *not* a fatal check on RC
578  m_env.worldRank(),
579  "GslMatrix::internalSvd()",
580  "matrix svd failed",
583  }
584 
585  return iRC;
586 }
GslMatrix * m_svdVmat
m_svdVmat stores the N-by-N orthogonal square matrix V after the singular value decomposition of a ma...
Definition: GslMatrix.h:437
unsigned int numRowsLocal() const
Returns the local row dimension of this matrix.
Definition: GslMatrix.C:275
GslMatrix * m_svdUmat
m_svdUmat stores the M-by-N orthogonal matrix U after the singular value decomposition of a matrix...
Definition: GslMatrix.h:424
GslMatrix()
Default Constructor.
gsl_matrix * data()
Returns this matrix.
Definition: GslMatrix.C:1766
int worldRank() const
Returns the same thing as fullRank()
Definition: Environment.C:262
GslMatrix * m_svdVTmat
m_svdVmatT stores the transpose of N-by-N orthogonal square matrix V, namely V^T, after the singular ...
Definition: GslMatrix.h:443
Map * m_svdColMap
Mapping for matrices involved in the singular value decomposition (svd) routine.
Definition: GslMatrix.h:418
const int UQ_MATRIX_SVD_FAILED_RC
Definition: Defines.h:103
gsl_vector * data() const
Definition: GslVector.C:970
GslMatrix transpose() const
This function calculated the transpose of this matrix (square).
Definition: GslMatrix.C:693
GslVector * m_svdSvec
m_svdSvec stores the diagonal of the N-by-N diagonal matrix of singular values S after the singular v...
Definition: GslMatrix.h:431
unsigned int numCols() const
Returns the column dimension of this matrix.
Definition: GslMatrix.C:287
const BaseEnvironment & m_env
QUESO environment variable.
Definition: Matrix.h:116
const Map & map() const
Definition: Matrix.C:54
GslMatrix QUESO::GslMatrix::inverse ( ) const

This function calculated the inverse of this matrix (square).

Definition at line 713 of file GslMatrix.C.

References QUESO::BaseEnvironment::checkingLevel(), QUESO::GslVector::cwSet(), QUESO::BaseEnvironment::displayVerbosity(), GslMatrix(), invertMultiply(), lnDeterminant(), QUESO::Matrix::m_env, m_inverse, QUESO::Matrix::m_map, numCols(), numRowsLocal(), QUESO::queso_require_equal_to_msg, and QUESO::BaseEnvironment::subDisplayFile().

714 {
715  unsigned int nRows = this->numRowsLocal();
716  unsigned int nCols = this->numCols();
717 
718  queso_require_equal_to_msg(nRows, nCols, "matrix is not square");
719 
720  if (m_inverse == NULL) {
721  m_inverse = new GslMatrix(m_env,m_map,nCols);
722  GslVector unitVector(m_env,m_map);
723  unitVector.cwSet(0.);
724  GslVector multVector(m_env,m_map);
725  for (unsigned int j = 0; j < nCols; ++j) {
726  if (j > 0) unitVector[j-1] = 0.;
727  unitVector[j] = 1.;
728  this->invertMultiply(unitVector, multVector);
729  for (unsigned int i = 0; i < nRows; ++i) {
730  (*m_inverse)(i,j) = multVector[i];
731  }
732  }
733  }
734  if (m_env.checkingLevel() >= 1) {
735  *m_env.subDisplayFile() << "CHECKING In GslMatrix::inverse()"
736  << ": M.lnDet = " << this->lnDeterminant()
737  << ", M^{-1}.lnDet = " << m_inverse->lnDeterminant()
738  << std::endl;
739  }
740  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 5)) {
741  *m_env.subDisplayFile() << "In GslMatrix::inverse():"
742  << "\n M = " << *this
743  << "\n M^{-1} = " << *m_inverse
744  << "\n M*M^{-1} = " << (*this)*(*m_inverse)
745  << "\n M^{-1}*M = " << (*m_inverse)*(*this)
746  << std::endl;
747  }
748 
749  return *m_inverse;
750 }
GslVector invertMultiply(const GslVector &b) const
This function calculates the inverse of this matrix and multiplies it with vector b...
Definition: GslMatrix.C:1222
const Map m_map
Mapping variable.
Definition: Matrix.h:121
unsigned int numRowsLocal() const
Returns the local row dimension of this matrix.
Definition: GslMatrix.C:275
GslMatrix()
Default Constructor.
unsigned int checkingLevel() const
Access function to private attribute m_checkingLevel.
Definition: Environment.C:464
MonteCarloSGOptions::MonteCarloSGOptions(const BaseEnvironment &env, const char *prefix, const McOptionsValues &alternativeOptionsValues queso_require_equal_to_msg)(m_env.optionsInputFileName(), std::string(""), std::string("this constructor is incompatible with the existence of an options input file"))
double lnDeterminant() const
Calculates the ln(determinant) of this matrix.
Definition: GslMatrix.C:1079
unsigned int displayVerbosity() const
Definition: Environment.C:450
unsigned int numCols() const
Returns the column dimension of this matrix.
Definition: GslMatrix.C:287
const BaseEnvironment & m_env
QUESO environment variable.
Definition: Matrix.h:116
GslMatrix * m_inverse
Inverse matrix of this.
Definition: GslMatrix.h:415
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:320
GslVector QUESO::GslMatrix::invertMultiply ( const GslVector b) const

This function calculates the inverse of this matrix and multiplies it with vector b.

It calls void GslMatrix::invertMultiply(const GslVector& b, GslVector& x) internally.

Definition at line 1222 of file GslMatrix.C.

References QUESO::Matrix::m_env, QUESO::Matrix::m_map, numCols(), QUESO::queso_require_equal_to_msg, and QUESO::GslVector::sizeLocal().

Referenced by determinant(), inverse(), invertMultiply(), and lnDeterminant().

1224 {
1225  queso_require_equal_to_msg(this->numCols(), b.sizeLocal(), "matrix and rhs have incompatible sizes");
1226 
1227  GslVector x(m_env,m_map);
1228  this->invertMultiply(b,x);
1229 
1230  return x;
1231 }
GslVector invertMultiply(const GslVector &b) const
This function calculates the inverse of this matrix and multiplies it with vector b...
Definition: GslMatrix.C:1222
const Map m_map
Mapping variable.
Definition: Matrix.h:121
MonteCarloSGOptions::MonteCarloSGOptions(const BaseEnvironment &env, const char *prefix, const McOptionsValues &alternativeOptionsValues queso_require_equal_to_msg)(m_env.optionsInputFileName(), std::string(""), std::string("this constructor is incompatible with the existence of an options input file"))
unsigned int numCols() const
Returns the column dimension of this matrix.
Definition: GslMatrix.C:287
const BaseEnvironment & m_env
QUESO environment variable.
Definition: Matrix.h:116
void QUESO::GslMatrix::invertMultiply ( const GslVector b,
GslVector x 
) const

This function calculates the inverse of this matrix, multiplies it with vector b and stores the result in vector x.

It checks for a previous LU decomposition of this matrix and does not recompute it if m_MU != NULL .

Definition at line 1234 of file GslMatrix.C.

References QUESO::GslVector::data(), QUESO::BaseEnvironment::displayVerbosity(), QUESO::Matrix::m_env, QUESO::Matrix::m_inDebugMode, m_isSingular, m_LU, m_mat, m_permutation, m_signum, QUESO::GslVector::norm2(), numCols(), numRowsLocal(), QUESO::queso_require_equal_to_msg, QUESO::GslVector::sizeLocal(), and QUESO::BaseEnvironment::subDisplayFile().

1237 {
1238  queso_require_equal_to_msg(this->numCols(), b.sizeLocal(), "matrix and rhs have incompatible sizes");
1239 
1240  queso_require_equal_to_msg(x.sizeLocal(), b.sizeLocal(), "solution and rhs have incompatible sizes");
1241 
1242  int iRC;
1243  if (m_LU == NULL) {
1244  queso_require_msg(!(m_permutation), "m_permutation should be NULL");
1245 
1246  m_LU = gsl_matrix_calloc(this->numRowsLocal(),this->numCols());
1247  queso_require_msg(m_LU, "gsl_matrix_calloc() failed");
1248 
1249  iRC = gsl_matrix_memcpy(m_LU, m_mat);
1250  queso_require_msg(!(iRC), "gsl_matrix_memcpy() failed");
1251 
1252  m_permutation = gsl_permutation_calloc(numCols());
1253  queso_require_msg(m_permutation, "gsl_permutation_calloc() failed");
1254 
1255  if (m_inDebugMode) {
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;
1260  }
1261 
1262  gsl_error_handler_t* oldHandler;
1263  oldHandler = gsl_set_error_handler_off();
1264  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 99)) {
1265  *m_env.subDisplayFile() << "In GslMatrix::invertMultiply()"
1266  << ": before 'gsl_linalg_LU_decomp()'"
1267  << std::endl;
1268  }
1269  iRC = gsl_linalg_LU_decomp(m_LU,m_permutation,&m_signum);
1270  if (iRC != 0) {
1271  std::cerr << "In GslMatrix::invertMultiply()"
1272  << ", after gsl_linalg_LU_decomp()"
1273  << ": iRC = " << iRC
1274  << ", gsl error message = " << gsl_strerror(iRC)
1275  << std::endl;
1276  }
1277  gsl_set_error_handler(oldHandler);
1278  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 99)) {
1279  *m_env.subDisplayFile() << "In GslMatrix::invertMultiply()"
1280  << ": after 'gsl_linalg_LU_decomp()'"
1281  << ", IRC = " << iRC
1282  << std::endl;
1283  }
1284  queso_require_msg(!(iRC), "gsl_linalg_LU_decomp() failed");
1285 
1286  if (m_inDebugMode) {
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;
1291  }
1292  }
1293 
1294  gsl_error_handler_t* oldHandler;
1295  oldHandler = gsl_set_error_handler_off();
1296  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 99)) {
1297  *m_env.subDisplayFile() << "In GslMatrix::invertMultiply()"
1298  << ": before 'gsl_linalg_LU_solve()'"
1299  << std::endl;
1300  }
1301  iRC = gsl_linalg_LU_solve(m_LU,m_permutation,b.data(),x.data());
1302  if (iRC != 0) {
1303  m_isSingular = true;
1304  std::cerr << "In GslMatrix::invertMultiply()"
1305  << ", after gsl_linalg_LU_solve()"
1306  << ": iRC = " << iRC
1307  << ", gsl error message = " << gsl_strerror(iRC)
1308  << std::endl;
1309  }
1310  gsl_set_error_handler(oldHandler);
1311  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 99)) {
1312  *m_env.subDisplayFile() << "In GslMatrix::invertMultiply()"
1313  << ": after 'gsl_linalg_LU_solve()'"
1314  << ", IRC = " << iRC
1315  << std::endl;
1316  }
1317 
1318 
1319  // m_env.worldRank(),
1320  // "GslMatrix::invertMultiply()",
1321  // "gsl_linalg_LU_solve() failed");
1322 
1323  if (m_inDebugMode) {
1324  GslVector tmpVec(b - (*this)*x);
1325  std::cout << "In GslMatrix::invertMultiply()"
1326  << ": ||b - Ax||_2 = " << tmpVec.norm2()
1327  << ": ||b - Ax||_2/||b||_2 = " << tmpVec.norm2()/b.norm2()
1328  << std::endl;
1329  }
1330 
1331  return;
1332 }
gsl_permutation * m_permutation
The permutation matrix of a LU decomposition.
Definition: GslMatrix.h:455
int m_signum
m_signum stores the sign of the permutation of the LU decomposition PA = LU.
Definition: GslMatrix.h:461
unsigned int numRowsLocal() const
Returns the local row dimension of this matrix.
Definition: GslMatrix.C:275
bool m_isSingular
Indicates whether or not this matrix is singular.
Definition: GslMatrix.h:464
MonteCarloSGOptions::MonteCarloSGOptions(const BaseEnvironment &env, const char *prefix, const McOptionsValues &alternativeOptionsValues queso_require_equal_to_msg)(m_env.optionsInputFileName(), std::string(""), std::string("this constructor is incompatible with the existence of an options input file"))
gsl_matrix * m_mat
GSL matrix, also referred to as this matrix.
Definition: GslMatrix.h:409
unsigned int displayVerbosity() const
Definition: Environment.C:450
unsigned int numCols() const
Returns the column dimension of this matrix.
Definition: GslMatrix.C:287
const BaseEnvironment & m_env
QUESO environment variable.
Definition: Matrix.h:116
gsl_matrix * m_LU
GSL matrix for the LU decomposition of m_mat.
Definition: GslMatrix.h:412
bool m_inDebugMode
Flag for either or not QUESO is in debug mode.
Definition: Matrix.h:134
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:320
GslMatrix QUESO::GslMatrix::invertMultiply ( const GslMatrix B) const

This function calculates the inverse of this matrix and multiplies it with matrix B.

It calls void GslMatrix::invertMultiply(const GslMatrix& B, GslMatrix& X) const internally.

Definition at line 1335 of file GslMatrix.C.

References invertMultiply(), QUESO::Matrix::m_env, QUESO::Matrix::m_map, and numCols().

1336 {
1337  GslMatrix X(m_env,m_map,B.numCols());
1338  this->invertMultiply(B,X);
1339 
1340  return X;
1341 }
GslVector invertMultiply(const GslVector &b) const
This function calculates the inverse of this matrix and multiplies it with vector b...
Definition: GslMatrix.C:1222
const Map m_map
Mapping variable.
Definition: Matrix.h:121
GslMatrix()
Default Constructor.
const BaseEnvironment & m_env
QUESO environment variable.
Definition: Matrix.h:116
void QUESO::GslMatrix::invertMultiply ( const GslMatrix B,
GslMatrix X 
) const

This function calculates the inverse of this matrix, multiplies it with matrix B and stores the result in matrix X.

It checks for a previous LU decomposition of this matrix and does not recompute it if m_MU != NULL .

Definition at line 1344 of file GslMatrix.C.

References getColumn(), invertMultiply(), QUESO::Matrix::m_env, QUESO::Matrix::m_map, numCols(), numRowsLocal(), QUESO::queso_require_equal_to_msg, and setColumn().

1345 {
1346  // Sanity Checks
1347  queso_require_equal_to_msg(B.numRowsLocal(), X.numRowsLocal(),
1348  "Matrices B and X are incompatible");
1349  queso_require_equal_to_msg(B.numCols(), X.numCols(),
1350  "Matrices B and X are incompatible");
1351  queso_require_equal_to_msg(this->numRowsLocal(), X.numRowsLocal(),
1352  "This and X matrices are incompatible");
1353 
1354  // Some local variables used within the loop.
1355  GslVector b(m_env, m_map);
1356  GslVector x(m_env, m_map);
1357 
1358  for( unsigned int j = 0; j < B.numCols(); ++j )
1359  {
1360  b = B.getColumn( j );
1361 
1362  //invertMultiply will only do the LU once and store it. So we don't
1363  //need to worry about it doing LU multiple times.
1364  x = this->invertMultiply( b );
1365 
1366  X.setColumn( j, x );
1367  }
1368 
1369  return;
1370 }
GslVector invertMultiply(const GslVector &b) const
This function calculates the inverse of this matrix and multiplies it with vector b...
Definition: GslMatrix.C:1222
const Map m_map
Mapping variable.
Definition: Matrix.h:121
unsigned int numRowsLocal() const
Returns the local row dimension of this matrix.
Definition: GslMatrix.C:275
MonteCarloSGOptions::MonteCarloSGOptions(const BaseEnvironment &env, const char *prefix, const McOptionsValues &alternativeOptionsValues queso_require_equal_to_msg)(m_env.optionsInputFileName(), std::string(""), std::string("this constructor is incompatible with the existence of an options input file"))
const BaseEnvironment & m_env
QUESO environment variable.
Definition: Matrix.h:116
GslVector QUESO::GslMatrix::invertMultiplyForceLU ( const GslVector b) const

This function calculates the inverse of this matrix and multiplies it with vector b.

  It calls void GslMatrix::InvertMultiplyForceLU(const GslVector& b, GslVector& x) const;(const GslVector& b,

GslVector& x) internally.

Definition at line 1373 of file GslMatrix.C.

References QUESO::Matrix::m_env, QUESO::Matrix::m_map, numCols(), QUESO::queso_require_equal_to_msg, and QUESO::GslVector::sizeLocal().

1375 {
1376  queso_require_equal_to_msg(this->numCols(), b.sizeLocal(), "matrix and rhs have incompatible sizes");
1377 
1378  GslVector x(m_env,m_map);
1379  this->invertMultiplyForceLU(b,x);
1380 
1381  return x;
1382 }
GslVector invertMultiplyForceLU(const GslVector &b) const
This function calculates the inverse of this matrix and multiplies it with vector b...
Definition: GslMatrix.C:1373
const Map m_map
Mapping variable.
Definition: Matrix.h:121
MonteCarloSGOptions::MonteCarloSGOptions(const BaseEnvironment &env, const char *prefix, const McOptionsValues &alternativeOptionsValues queso_require_equal_to_msg)(m_env.optionsInputFileName(), std::string(""), std::string("this constructor is incompatible with the existence of an options input file"))
unsigned int numCols() const
Returns the column dimension of this matrix.
Definition: GslMatrix.C:287
const BaseEnvironment & m_env
QUESO environment variable.
Definition: Matrix.h:116
void QUESO::GslMatrix::invertMultiplyForceLU ( const GslVector b,
GslVector x 
) const

This function calculates the inverse of this matrix, multiplies it with vector b and stores the result in vector x.

It recalculates the LU decomposition of this matrix.

Definition at line 1385 of file GslMatrix.C.

References QUESO::GslVector::data(), m_isSingular, m_LU, m_mat, m_permutation, m_signum, numCols(), numRowsLocal(), QUESO::queso_require_equal_to_msg, and QUESO::GslVector::sizeLocal().

1388 {
1389  queso_require_equal_to_msg(this->numCols(), b.sizeLocal(), "matrix and rhs have incompatible sizes");
1390 
1391  queso_require_equal_to_msg(x.sizeLocal(), b.sizeLocal(), "solution and rhs have incompatible sizes");
1392 
1393  int iRC;
1394 
1395  if ( m_LU == NULL ) {
1396  queso_require_msg(!(m_permutation), "m_permutation should be NULL");
1397  m_LU = gsl_matrix_calloc(this->numRowsLocal(),this->numCols());
1398  }
1399  queso_require_msg(m_LU, "gsl_matrix_calloc() failed");
1400 
1401  iRC = gsl_matrix_memcpy(m_LU, m_mat);
1402  queso_require_msg(!(iRC), "gsl_matrix_memcpy() failed");
1403 
1404  if( m_permutation == NULL ) m_permutation = gsl_permutation_calloc(numCols());
1405  queso_require_msg(m_permutation, "gsl_permutation_calloc() failed");
1406 
1407  iRC = gsl_linalg_LU_decomp(m_LU,m_permutation,&m_signum);
1408  queso_require_msg(!(iRC), "gsl_linalg_LU_decomp() failed");
1409 
1410  iRC = gsl_linalg_LU_solve(m_LU,m_permutation,b.data(),x.data());
1411  if (iRC != 0) {
1412  m_isSingular = true;
1413  }
1414  queso_require_msg(!(iRC), "gsl_linalg_LU_solve() failed");
1415 
1416  return;
1417 }
gsl_permutation * m_permutation
The permutation matrix of a LU decomposition.
Definition: GslMatrix.h:455
int m_signum
m_signum stores the sign of the permutation of the LU decomposition PA = LU.
Definition: GslMatrix.h:461
unsigned int numRowsLocal() const
Returns the local row dimension of this matrix.
Definition: GslMatrix.C:275
bool m_isSingular
Indicates whether or not this matrix is singular.
Definition: GslMatrix.h:464
MonteCarloSGOptions::MonteCarloSGOptions(const BaseEnvironment &env, const char *prefix, const McOptionsValues &alternativeOptionsValues queso_require_equal_to_msg)(m_env.optionsInputFileName(), std::string(""), std::string("this constructor is incompatible with the existence of an options input file"))
gsl_matrix * m_mat
GSL matrix, also referred to as this matrix.
Definition: GslMatrix.h:409
unsigned int numCols() const
Returns the column dimension of this matrix.
Definition: GslMatrix.C:287
gsl_matrix * m_LU
GSL matrix for the LU decomposition of m_mat.
Definition: GslMatrix.h:412
void QUESO::GslMatrix::largestEigen ( double &  eigenValue,
GslVector eigenVector 
) const

This function finds largest eigenvalue, namely eigenValue, of this matrix and its corresponding eigenvector, namely eigenVector.

Definition at line 1446 of file GslMatrix.C.

References QUESO::GslVector::abs(), QUESO::Matrix::m_env, QUESO::Matrix::m_map, QUESO::queso_require_not_equal_to_msg, and QUESO::GslVector::sizeLocal().

1447 {
1448 
1449  // Sanity Checks
1450  unsigned int n = eigenVector.sizeLocal();
1451 
1452  queso_require_not_equal_to_msg(n, 0, "invalid input vector size");
1453 
1454  /* The following notation is used:
1455  z = vector used in iteration that ends up being the eigenvector corresponding to the
1456  largest eigenvalue
1457  w = vector used in iteration that we extract the largest eigenvalue from. */
1458 
1459  // Some parameters associated with the algorithm
1460  // TODO: Do we want to add the ability to have these set by the user?
1461  const unsigned int max_num_iterations = 10000;
1462  const double tolerance = 1.0e-13;
1463 
1464  // Create temporary working vectors.
1465  // TODO: We shouldn't have to use these - we ought to be able to work directly
1466  // TODO: with eigenValue and eigenVector. Optimize this once we have regression
1467  // TODO: tests.
1468  GslVector z(m_env, m_map, 1.0 ); // Needs to be initialized to 1.0
1469  GslVector w(m_env, m_map);
1470 
1471  // Some variables we use inside the loop.
1472  int index;
1473  double residual;
1474  double lambda;
1475 
1476  for( unsigned int k = 0; k < max_num_iterations; ++k )
1477  {
1478  w = (*this) * z;
1479 
1480  // For this algorithm, it's crucial to get the maximum in
1481  // absolute value, but then to normalize by the actual value
1482  // and *not* the absolute value.
1483  index = (w.abs()).getMaxValueIndex();
1484 
1485  lambda = w[index];
1486 
1487  z = (1.0/lambda) * w;
1488 
1489  // Here we use the norm of the residual as our convergence check:
1490  // norm( A*x - \lambda*x )
1491  residual = ( (*this)*z - lambda*z ).norm2();
1492 
1493  if( residual < tolerance )
1494  {
1495  eigenValue = lambda;
1496 
1497  // TODO: Do we want to normalize this so eigenVector.norm2() = 1?
1498  eigenVector = z;
1499  return;
1500  }
1501 
1502  }
1503 
1504  // If we reach this point, then we didn't converge. Print error message
1505  // to this effect.
1506  // TODO: We know we failed at this point - need a way to not need a test
1507  // TODO: Just specify the error.
1508  queso_require_less_msg(residual, tolerance, "Maximum num iterations exceeded");
1509 
1510 
1511  return;
1512 }
const Map m_map
Mapping variable.
Definition: Matrix.h:121
MonteCarloSGOptions::MonteCarloSGOptions(const BaseEnvironment &env, const char *prefix queso_require_not_equal_to_msg)(m_env.optionsInputFileName(), std::string(""), std::string("this constructor is incompatible with the absence of an options input file"))
const BaseEnvironment & m_env
QUESO environment variable.
Definition: Matrix.h:116
double QUESO::GslMatrix::lnDeterminant ( ) const

Calculates the ln(determinant) of this matrix.

Definition at line 1079 of file GslMatrix.C.

References QUESO::BaseEnvironment::displayVerbosity(), invertMultiply(), m_determinant, QUESO::Matrix::m_env, m_lnDeterminant, m_LU, QUESO::Matrix::m_map, m_signum, and QUESO::BaseEnvironment::subDisplayFile().

Referenced by inverse().

1080 {
1081  if (m_lnDeterminant == -INFINITY) {
1082  if (m_LU == NULL) {
1083  GslVector tmpB(m_env,m_map);
1084  GslVector tmpX(m_env,m_map);
1085  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 99)) {
1086  *m_env.subDisplayFile() << "In GslMatrix::lnDeterminant()"
1087  << ": before 'this->invertMultiply()'"
1088  << std::endl;
1089  }
1090  this->invertMultiply(tmpB,tmpX);
1091  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 99)) {
1092  *m_env.subDisplayFile() << "In GslMatrix::lnDeterminant()"
1093  << ": after 'this->invertMultiply()'"
1094  << std::endl;
1095  }
1096  }
1097  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 99)) {
1098  *m_env.subDisplayFile() << "In GslMatrix::lnDeterminant()"
1099  << ": before 'gsl_linalg_LU_det()'"
1100  << std::endl;
1101  }
1102  m_determinant = gsl_linalg_LU_det(m_LU,m_signum);
1103  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 99)) {
1104  *m_env.subDisplayFile() << "In GslMatrix::lnDeterminant()"
1105  << ": after 'gsl_linalg_LU_det()'"
1106  << std::endl;
1107  }
1108  m_lnDeterminant = gsl_linalg_LU_lndet(m_LU);
1109  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 99)) {
1110  *m_env.subDisplayFile() << "In GslMatrix::lnDeterminant()"
1111  << ": before 'gsl_linalg_LU_lndet()'"
1112  << std::endl;
1113  }
1114  }
1115 
1116  return m_lnDeterminant;
1117 }
GslVector invertMultiply(const GslVector &b) const
This function calculates the inverse of this matrix and multiplies it with vector b...
Definition: GslMatrix.C:1222
const Map m_map
Mapping variable.
Definition: Matrix.h:121
int m_signum
m_signum stores the sign of the permutation of the LU decomposition PA = LU.
Definition: GslMatrix.h:461
double m_lnDeterminant
The natural logarithm of the determinant of this matrix.
Definition: GslMatrix.h:449
double m_determinant
The determinant of this matrix.
Definition: GslMatrix.h:446
unsigned int displayVerbosity() const
Definition: Environment.C:450
const BaseEnvironment & m_env
QUESO environment variable.
Definition: Matrix.h:116
gsl_matrix * m_LU
GSL matrix for the LU decomposition of m_mat.
Definition: GslMatrix.h:412
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:320
void QUESO::GslMatrix::matlabLinearInterpExtrap ( const GslVector x1Vec,
const GslMatrix y1Mat,
const GslVector x2Vec 
)

Definition at line 1741 of file GslMatrix.C.

References getColumn(), QUESO::GslVector::matlabLinearInterpExtrap(), numCols(), numRowsLocal(), QUESO::queso_require_equal_to_msg, setColumn(), and QUESO::GslVector::sizeLocal().

1745 {
1746  queso_require_greater_msg(x1Vec.sizeLocal(), 1, "invalid 'x1' size");
1747 
1748  queso_require_equal_to_msg(x1Vec.sizeLocal(), y1Mat.numRowsLocal(), "invalid 'x1' and 'y1' sizes");
1749 
1750  queso_require_equal_to_msg(x2Vec.sizeLocal(), this->numRowsLocal(), "invalid 'x2' and 'this' sizes");
1751 
1752  queso_require_equal_to_msg(y1Mat.numCols(), this->numCols(), "invalid 'y1' and 'this' sizes");
1753 
1754  GslVector y1Vec(x1Vec);
1755  GslVector y2Vec(x2Vec);
1756  for (unsigned int colId = 0; colId < y1Mat.numCols(); ++colId) {
1757  y1Mat.getColumn(colId,y1Vec);
1758  y2Vec.matlabLinearInterpExtrap(x1Vec,y1Vec,x2Vec);
1759  this->setColumn(colId,y2Vec);
1760  }
1761 
1762  return;
1763 }
unsigned int numRowsLocal() const
Returns the local row dimension of this matrix.
Definition: GslMatrix.C:275
MonteCarloSGOptions::MonteCarloSGOptions(const BaseEnvironment &env, const char *prefix, const McOptionsValues &alternativeOptionsValues queso_require_equal_to_msg)(m_env.optionsInputFileName(), std::string(""), std::string("this constructor is incompatible with the existence of an options input file"))
unsigned int numCols() const
Returns the column dimension of this matrix.
Definition: GslMatrix.C:287
void setColumn(const unsigned int column_num, const GslVector &column)
This function copies vector column into the column_num-th column of this matrix.
Definition: GslMatrix.C:1681
double QUESO::GslMatrix::max ( ) const

Returns the maximum element value of the matrix.

Definition at line 329 of file GslMatrix.C.

References numCols(), and numRowsLocal().

330 {
331  double value = -INFINITY;
332 
333  unsigned int nRows = this->numRowsLocal();
334  unsigned int nCols = this->numCols();
335  double aux = 0.;
336  for (unsigned int i = 0; i < nRows; i++) {
337  for (unsigned int j = 0; j < nCols; j++) {
338  aux = (*this)(i,j);
339  if (aux > value) value = aux;
340  }
341  }
342 
343  return value;
344 }
unsigned int numRowsLocal() const
Returns the local row dimension of this matrix.
Definition: GslMatrix.C:275
unsigned int numCols() const
Returns the column dimension of this matrix.
Definition: GslMatrix.C:287
void QUESO::GslMatrix::mpiSum ( const MpiComm comm,
GslMatrix M_global 
) const

Definition at line 1698 of file GslMatrix.C.

References QUESO::MpiComm::Allreduce(), QUESO::Matrix::env(), QUESO::BaseEnvironment::fullRank(), numCols(), numRowsLocal(), and QUESO::size.

1699 {
1700  // Sanity Checks
1701  UQ_FATAL_RC_MACRO(((this->numRowsLocal() != M_global.numRowsLocal()) ||
1702  (this->numCols() != M_global.numCols() )),
1703  env().fullRank(),
1704  "GslMatrix::mpiSum()",
1705  "local and global matrices incompatible");
1706 
1707  /* TODO: Probably a better way to handle this unpacking/packing of data */
1708  int size = M_global.numRowsLocal()*M_global.numCols();
1709  std::vector<double> local( size, 0.0 );
1710  std::vector<double> global( size, 0.0 );
1711 
1712  int k;
1713  for( unsigned int i = 0; i < this->numRowsLocal(); i++ )
1714  {
1715  for( unsigned int j = 0; j < this->numCols(); j++ )
1716  {
1717  k = i + j*M_global.numCols();
1718 
1719  local[k] = (*this)(i,j);
1720  }
1721  }
1722 
1723  comm.Allreduce<double>(&local[0], &global[0], size, RawValue_MPI_SUM,
1724  "GslMatrix::mpiSum()",
1725  "failed MPI.Allreduce()");
1726 
1727  for( unsigned int i = 0; i < this->numRowsLocal(); i++ )
1728  {
1729  for( unsigned int j = 0; j < this->numCols(); j++ )
1730  {
1731  k = i + j*M_global.numCols();
1732 
1733  M_global(i,j) = global[k];
1734  }
1735  }
1736 
1737  return;
1738 }
const BaseEnvironment & env() const
Definition: Matrix.C:47
unsigned int numRowsLocal() const
Returns the local row dimension of this matrix.
Definition: GslMatrix.C:275
int fullRank() const
Returns the rank of the MPI process in QUESO&#39;s full communicator.
Definition: Environment.C:268
unsigned int numCols() const
Returns the column dimension of this matrix.
Definition: GslMatrix.C:287
GslVector QUESO::GslMatrix::multiply ( const GslVector x) const

This function multiplies this matrix by vector x and returns the resulting vector.

Definition at line 1155 of file GslMatrix.C.

References QUESO::Matrix::m_env, QUESO::Matrix::m_map, numCols(), QUESO::queso_require_equal_to_msg, and QUESO::GslVector::sizeLocal().

Referenced by multiply(), and QUESO::operator*().

1157 {
1158  queso_require_equal_to_msg(this->numCols(), x.sizeLocal(), "matrix and vector have incompatible sizes");
1159 
1160  GslVector y(m_env,m_map);
1161  this->multiply(x,y);
1162 
1163  return y;
1164 }
const Map m_map
Mapping variable.
Definition: Matrix.h:121
GslVector multiply(const GslVector &x) const
This function multiplies this matrix by vector x and returns the resulting vector.
Definition: GslMatrix.C:1155
MonteCarloSGOptions::MonteCarloSGOptions(const BaseEnvironment &env, const char *prefix, const McOptionsValues &alternativeOptionsValues queso_require_equal_to_msg)(m_env.optionsInputFileName(), std::string(""), std::string("this constructor is incompatible with the existence of an options input file"))
unsigned int numCols() const
Returns the column dimension of this matrix.
Definition: GslMatrix.C:287
const BaseEnvironment & m_env
QUESO environment variable.
Definition: Matrix.h:116
void QUESO::GslMatrix::multiply ( const GslVector x,
GslVector y 
) const

This function multiplies this matrix by vector x and fills.

Definition at line 1167 of file GslMatrix.C.

References numCols(), numRowsLocal(), QUESO::queso_require_equal_to_msg, and QUESO::GslVector::sizeLocal().

1170 {
1171  queso_require_equal_to_msg(this->numCols(), x.sizeLocal(), "matrix and x have incompatible sizes");
1172 
1173  queso_require_equal_to_msg(this->numRowsLocal(), y.sizeLocal(), "matrix and y have incompatible sizes");
1174 
1175  unsigned int sizeX = this->numCols();
1176  unsigned int sizeY = this->numRowsLocal();
1177  for (unsigned int i = 0; i < sizeY; ++i) {
1178  double value = 0.;
1179  for (unsigned int j = 0; j < sizeX; ++j) {
1180  value += (*this)(i,j)*x[j];
1181  }
1182  y[i] = value;
1183  }
1184 
1185  return;
1186 }
unsigned int numRowsLocal() const
Returns the local row dimension of this matrix.
Definition: GslMatrix.C:275
MonteCarloSGOptions::MonteCarloSGOptions(const BaseEnvironment &env, const char *prefix, const McOptionsValues &alternativeOptionsValues queso_require_equal_to_msg)(m_env.optionsInputFileName(), std::string(""), std::string("this constructor is incompatible with the existence of an options input file"))
unsigned int numCols() const
Returns the column dimension of this matrix.
Definition: GslMatrix.C:287
GslMatrix QUESO::GslMatrix::multiply ( const GslMatrix X) const

This function multiplies this matrix by matrix X and returns the resulting matrix.

Definition at line 1189 of file GslMatrix.C.

References QUESO::Matrix::m_env, QUESO::Matrix::m_map, multiply(), and numCols().

1191 {
1192  GslMatrix Y(m_env,m_map,X.numCols());
1193  this->multiply(X,Y);
1194 
1195  return Y;
1196 }
const Map m_map
Mapping variable.
Definition: Matrix.h:121
GslVector multiply(const GslVector &x) const
This function multiplies this matrix by vector x and returns the resulting vector.
Definition: GslMatrix.C:1155
GslMatrix()
Default Constructor.
const BaseEnvironment & m_env
QUESO environment variable.
Definition: Matrix.h:116
void QUESO::GslMatrix::multiply ( const GslMatrix X,
GslMatrix Y 
) const

This function multiplies this matrix by matrix X and fills.

Definition at line 1201 of file GslMatrix.C.

References numCols(), numRowsGlobal(), and QUESO::queso_require_equal_to_msg.

1204 {
1205  queso_require_equal_to_msg(this->numCols(), X.numRowsGlobal(), "matrix and X have incompatible sizes");
1206  queso_require_equal_to_msg(this->numRowsGlobal(), Y.numRowsGlobal(), "matrix and Y have incompatible sizes");
1207  queso_require_equal_to_msg(X.numCols(), Y.numCols(), "X and Y have incompatible sizes");
1208 
1209  const unsigned int m_s = this->numRowsGlobal();
1210  const unsigned int p_s = this->numCols();
1211  const unsigned int n_s = X.numCols();
1212 
1213  for (unsigned int k=0; k<p_s; k++)
1214  for (unsigned int j=0; j<n_s; j++)
1215  if (X(k,j) != 0.)
1216  for (unsigned int i=0; i<m_s; i++)
1217  Y(i,j) += (*this)(i,k) * X(k,j);
1218 }
unsigned int numRowsGlobal() const
Returns the global row dimension of this matrix.
Definition: GslMatrix.C:281
MonteCarloSGOptions::MonteCarloSGOptions(const BaseEnvironment &env, const char *prefix, const McOptionsValues &alternativeOptionsValues queso_require_equal_to_msg)(m_env.optionsInputFileName(), std::string(""), std::string("this constructor is incompatible with the existence of an options input file"))
unsigned int numCols() const
Returns the column dimension of this matrix.
Definition: GslMatrix.C:287
double QUESO::GslMatrix::normFrob ( ) const

Returns the Frobenius norm of this matrix.

Definition at line 293 of file GslMatrix.C.

References numCols(), and numRowsLocal().

294 {
295  double value = 0.;
296 
297  unsigned int nRows = this->numRowsLocal();
298  unsigned int nCols = this->numCols();
299  double aux = 0.;
300  for (unsigned int i = 0; i < nRows; i++) {
301  for (unsigned int j = 0; j < nCols; j++) {
302  aux = (*this)(i,j);
303  value += aux*aux;
304  }
305  }
306 
307  return sqrt(value);
308 }
unsigned int numRowsLocal() const
Returns the local row dimension of this matrix.
Definition: GslMatrix.C:275
unsigned int numCols() const
Returns the column dimension of this matrix.
Definition: GslMatrix.C:287
double QUESO::GslMatrix::normMax ( ) const

Returns the Frobenius norm of this matrix.

Definition at line 311 of file GslMatrix.C.

References numCols(), and numRowsLocal().

312 {
313  double value = 0.;
314 
315  unsigned int nRows = this->numRowsLocal();
316  unsigned int nCols = this->numCols();
317  double aux = 0.;
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;
322  }
323  }
324 
325  return value;
326 }
unsigned int numRowsLocal() const
Returns the local row dimension of this matrix.
Definition: GslMatrix.C:275
unsigned int numCols() const
Returns the column dimension of this matrix.
Definition: GslMatrix.C:287
unsigned int QUESO::GslMatrix::numRowsGlobal ( ) const
virtual

Returns the global row dimension of this matrix.

Implements QUESO::Matrix.

Definition at line 281 of file GslMatrix.C.

References m_mat.

Referenced by multiply(), and transpose().

282 {
283  return m_mat->size1;
284 }
gsl_matrix * m_mat
GSL matrix, also referred to as this matrix.
Definition: GslMatrix.h:409
double& QUESO::GslMatrix::operator() ( unsigned int  i,
unsigned int  j 
)
inline

Element access method (non-const).

Definition at line 104 of file GslMatrix.h.

105  {
106  // Changing the matrix changes its decomposition(s). We'll
107  // invalidate for now; recalculate later if needed.
108  this->resetLU();
109  queso_require_less_msg(i, m_mat->size1, "i is too large");
110  queso_require_less_msg(j, m_mat->size2, "j is too large");
111  return *gsl_matrix_ptr(m_mat,i,j);
112  }
void resetLU()
In this function resets the LU decomposition of this matrix, as well as deletes the private member po...
Definition: GslMatrix.C:232
gsl_matrix * m_mat
GSL matrix, also referred to as this matrix.
Definition: GslMatrix.h:409
const double& QUESO::GslMatrix::operator() ( unsigned int  i,
unsigned int  j 
) const
inline

Element access method (const).

Definition at line 115 of file GslMatrix.h.

References m_mat.

116  {
117  queso_require_less_msg(i, m_mat->size1, "i is too large");
118  queso_require_less_msg(j, m_mat->size2, "j is too large");
119  return *gsl_matrix_ptr(m_mat,i,j);
120  }
gsl_matrix * m_mat
GSL matrix, also referred to as this matrix.
Definition: GslMatrix.h:409
GslMatrix & QUESO::GslMatrix::operator*= ( double  a)

Stores in this the coordinate-wise multiplication of this and a.

Definition at line 177 of file GslMatrix.C.

References m_mat, and resetLU().

178 {
179  this->resetLU();
180  int iRC;
181  iRC = gsl_matrix_scale(m_mat,a);
182  queso_require_msg(!(iRC), "scaling failed");
183  return *this;
184 }
void resetLU()
In this function resets the LU decomposition of this matrix, as well as deletes the private member po...
Definition: GslMatrix.C:232
gsl_matrix * m_mat
GSL matrix, also referred to as this matrix.
Definition: GslMatrix.h:409
GslMatrix & QUESO::GslMatrix::operator+= ( const GslMatrix rhs)

Stores in this the coordinate-wise addition of this and rhs.

Definition at line 196 of file GslMatrix.C.

References m_mat, and resetLU().

197 {
198  this->resetLU();
199  int iRC;
200  iRC = gsl_matrix_add(m_mat,rhs.m_mat);
201  queso_require_msg(!(iRC), "failed");
202 
203  return *this;
204 }
void resetLU()
In this function resets the LU decomposition of this matrix, as well as deletes the private member po...
Definition: GslMatrix.C:232
gsl_matrix * m_mat
GSL matrix, also referred to as this matrix.
Definition: GslMatrix.h:409
GslMatrix & QUESO::GslMatrix::operator-= ( const GslMatrix rhs)

Stores in this the coordinate-wise subtraction of this by rhs.

Definition at line 207 of file GslMatrix.C.

References m_mat, and resetLU().

208 {
209  this->resetLU();
210  int iRC;
211  iRC = gsl_matrix_sub(m_mat,rhs.m_mat);
212  queso_require_msg(!(iRC), "failed");
213 
214  return *this;
215 }
void resetLU()
In this function resets the LU decomposition of this matrix, as well as deletes the private member po...
Definition: GslMatrix.C:232
gsl_matrix * m_mat
GSL matrix, also referred to as this matrix.
Definition: GslMatrix.h:409
GslMatrix & QUESO::GslMatrix::operator/= ( double  a)

Stores in this the coordinate-wise division of this by a.

Definition at line 187 of file GslMatrix.C.

References resetLU().

188 {
189  this->resetLU();
190  *this *= (1./a);
191 
192  return *this;
193 }
void resetLU()
In this function resets the LU decomposition of this matrix, as well as deletes the private member po...
Definition: GslMatrix.C:232
GslMatrix & QUESO::GslMatrix::operator= ( const GslMatrix rhs)

Copies values from matrix rhs to this.

Definition at line 169 of file GslMatrix.C.

References copy(), and resetLU().

170 {
171  this->resetLU();
172  this->copy(obj);
173  return *this;
174 }
void resetLU()
In this function resets the LU decomposition of this matrix, as well as deletes the private member po...
Definition: GslMatrix.C:232
void copy(const GslMatrix &src)
In this function this matrix receives a copy of matrix src.
Definition: GslMatrix.C:220
void QUESO::GslMatrix::print ( std::ostream &  os) const
virtual

Print method. Defines the behavior of the ostream << operator inherited from the Object class.

Implements QUESO::Matrix.

Definition at line 1772 of file GslMatrix.C.

References QUESO::Matrix::m_printHorizontally, numCols(), and numRowsLocal().

Referenced by QUESO::operator<<(), and QUESO::GslBlockMatrix::print().

1773 {
1774  unsigned int nRows = this->numRowsLocal();
1775  unsigned int nCols = this->numCols();
1776 
1777  if (m_printHorizontally) {
1778  for (unsigned int i = 0; i < nRows; ++i) {
1779  for (unsigned int j = 0; j < nCols; ++j) {
1780  os << (*this)(i,j)
1781  << " ";
1782  }
1783  if (i != (nRows-1)) os << "; ";
1784  }
1785  //os << std::endl;
1786  }
1787  else {
1788  for (unsigned int i = 0; i < nRows; ++i) {
1789  for (unsigned int j = 0; j < nCols; ++j) {
1790  os << (*this)(i,j)
1791  << " ";
1792  }
1793  os << std::endl;
1794  }
1795  }
1796 
1797  return;
1798 }
unsigned int numRowsLocal() const
Returns the local row dimension of this matrix.
Definition: GslMatrix.C:275
bool m_printHorizontally
Flag for either or not print this matrix.
Definition: Matrix.h:131
unsigned int numCols() const
Returns the column dimension of this matrix.
Definition: GslMatrix.C:287
unsigned int QUESO::GslMatrix::rank ( double  absoluteZeroThreshold,
double  relativeZeroThreshold 
) const

This function returns the number of singular values of this matrix (rank).

The rank function provides an estimate of the number of linearly independent rows or columns of a full matrix.

Definition at line 1120 of file GslMatrix.C.

References QUESO::BaseEnvironment::displayVerbosity(), internalSvd(), QUESO::Matrix::m_env, m_svdSvec, numCols(), numRowsLocal(), QUESO::GslVector::sizeLocal(), and QUESO::BaseEnvironment::subDisplayFile().

1121 {
1122  int iRC = 0;
1123  iRC = internalSvd();
1124  if (iRC) {}; // just to remove compiler warning
1125 
1126  GslVector relativeVec(*m_svdSvec);
1127  if (relativeVec[0] > 0.) {
1128  relativeVec = (1./relativeVec[0])*relativeVec;
1129  }
1130 
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 )) {
1135  rankValue += 1;
1136  }
1137  }
1138 
1139  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 3)) {
1140  *m_env.subDisplayFile() << "In GslMatrix::rank()"
1141  << ": this->numRowsLocal() = " << this->numRowsLocal()
1142  << ", this->numCols() = " << this->numCols()
1143  << ", absoluteZeroThreshold = " << absoluteZeroThreshold
1144  << ", relativeZeroThreshold = " << relativeZeroThreshold
1145  << ", rankValue = " << rankValue
1146  << ", m_svdSvec = " << *m_svdSvec
1147  << ", relativeVec = " << relativeVec
1148  << std::endl;
1149  }
1150 
1151  return rankValue;
1152 }
unsigned int numRowsLocal() const
Returns the local row dimension of this matrix.
Definition: GslMatrix.C:275
GslVector * m_svdSvec
m_svdSvec stores the diagonal of the N-by-N diagonal matrix of singular values S after the singular v...
Definition: GslMatrix.h:431
unsigned int displayVerbosity() const
Definition: Environment.C:450
unsigned int numCols() const
Returns the column dimension of this matrix.
Definition: GslMatrix.C:287
const BaseEnvironment & m_env
QUESO environment variable.
Definition: Matrix.h:116
int internalSvd() const
This function factorizes the M-by-N matrix A into the singular value decomposition A = U S V^T for M ...
Definition: GslMatrix.C:532
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:320
void QUESO::GslMatrix::resetLU ( )
private

In this function resets the LU decomposition of this matrix, as well as deletes the private member pointers, if existing.

Definition at line 232 of file GslMatrix.C.

References m_determinant, m_inverse, m_isSingular, m_lnDeterminant, m_LU, m_permutation, m_signum, m_svdColMap, m_svdSvec, m_svdUmat, m_svdVmat, and m_svdVTmat.

Referenced by copy(), operator*=(), operator+=(), operator-=(), operator/=(), operator=(), setColumn(), setRow(), zeroLower(), zeroUpper(), and ~GslMatrix().

233 {
234  if (m_LU) {
235  gsl_matrix_free(m_LU);
236  m_LU = NULL;
237  }
238  if (m_inverse) {
239  delete m_inverse;
240  m_inverse = NULL;
241  }
242  if (m_svdColMap) {
243  delete m_svdColMap;
244  m_svdColMap = NULL;
245  }
246  if (m_svdUmat) {
247  delete m_svdUmat;
248  m_svdUmat = NULL;
249  }
250  if (m_svdSvec) {
251  delete m_svdSvec;
252  m_svdSvec = NULL;
253  }
254  if (m_svdVmat) {
255  delete m_svdVmat;
256  m_svdVmat = NULL;
257  }
258  if (m_svdVTmat) {
259  delete m_svdVTmat;
260  m_svdVTmat = NULL;
261  }
262  m_determinant = -INFINITY;
263  m_lnDeterminant = -INFINITY;
264  if (m_permutation) {
265  gsl_permutation_free(m_permutation);
266  m_permutation = NULL;
267  }
268  m_signum = 0;
269  m_isSingular = false;
270 
271  return;
272 }
gsl_permutation * m_permutation
The permutation matrix of a LU decomposition.
Definition: GslMatrix.h:455
GslMatrix * m_svdVmat
m_svdVmat stores the N-by-N orthogonal square matrix V after the singular value decomposition of a ma...
Definition: GslMatrix.h:437
int m_signum
m_signum stores the sign of the permutation of the LU decomposition PA = LU.
Definition: GslMatrix.h:461
GslMatrix * m_svdUmat
m_svdUmat stores the M-by-N orthogonal matrix U after the singular value decomposition of a matrix...
Definition: GslMatrix.h:424
bool m_isSingular
Indicates whether or not this matrix is singular.
Definition: GslMatrix.h:464
GslMatrix * m_svdVTmat
m_svdVmatT stores the transpose of N-by-N orthogonal square matrix V, namely V^T, after the singular ...
Definition: GslMatrix.h:443
Map * m_svdColMap
Mapping for matrices involved in the singular value decomposition (svd) routine.
Definition: GslMatrix.h:418
double m_lnDeterminant
The natural logarithm of the determinant of this matrix.
Definition: GslMatrix.h:449
GslVector * m_svdSvec
m_svdSvec stores the diagonal of the N-by-N diagonal matrix of singular values S after the singular v...
Definition: GslMatrix.h:431
double m_determinant
The determinant of this matrix.
Definition: GslMatrix.h:446
gsl_matrix * m_LU
GSL matrix for the LU decomposition of m_mat.
Definition: GslMatrix.h:412
GslMatrix * m_inverse
Inverse matrix of this.
Definition: GslMatrix.h:415
void QUESO::GslMatrix::setColumn ( const unsigned int  column_num,
const GslVector column 
)

This function copies vector column into the column_num-th column of this matrix.

Definition at line 1681 of file GslMatrix.C.

References QUESO::GslVector::data(), m_mat, numCols(), numRowsLocal(), QUESO::queso_require_equal_to_msg, resetLU(), and QUESO::GslVector::sizeLocal().

Referenced by invertMultiply(), matlabLinearInterpExtrap(), and svdSolve().

1682 {
1683  this->resetLU();
1684  // Sanity checks
1685  queso_require_less_msg(column_num, this->numCols(), "Specified column number not within range");
1686 
1687  queso_require_equal_to_msg(column.sizeLocal(), this->numRowsLocal(), "column vector not same size as this matrix");
1688 
1689  gsl_vector* gsl_column = column.data();
1690 
1691  int error_code = gsl_matrix_set_col( m_mat, column_num, gsl_column );
1692  queso_require_equal_to_msg(error_code, 0, "gsl_matrix_set_col failed");
1693 
1694  return;
1695 }
unsigned int numRowsLocal() const
Returns the local row dimension of this matrix.
Definition: GslMatrix.C:275
void resetLU()
In this function resets the LU decomposition of this matrix, as well as deletes the private member po...
Definition: GslMatrix.C:232
MonteCarloSGOptions::MonteCarloSGOptions(const BaseEnvironment &env, const char *prefix, const McOptionsValues &alternativeOptionsValues queso_require_equal_to_msg)(m_env.optionsInputFileName(), std::string(""), std::string("this constructor is incompatible with the existence of an options input file"))
gsl_matrix * m_mat
GSL matrix, also referred to as this matrix.
Definition: GslMatrix.h:409
unsigned int numCols() const
Returns the column dimension of this matrix.
Definition: GslMatrix.C:287
void QUESO::GslMatrix::setRow ( const unsigned int  row_num,
const GslVector row 
)

This function copies vector row into the row_num-th column of this matrix.

Definition at line 1664 of file GslMatrix.C.

References QUESO::GslVector::data(), m_mat, numCols(), numRowsLocal(), QUESO::queso_require_equal_to_msg, resetLU(), and QUESO::GslVector::sizeLocal().

1665 {
1666  this->resetLU();
1667  // Sanity checks
1668  queso_require_less_msg(row_num, this->numRowsLocal(), "Specified row number not within range");
1669 
1670  queso_require_equal_to_msg(row.sizeLocal(), this->numCols(), "row vector not same size as this matrix");
1671 
1672  gsl_vector* gsl_row = row.data();
1673 
1674  int error_code = gsl_matrix_set_row( m_mat, row_num, gsl_row );
1675  queso_require_equal_to_msg(error_code, 0, "gsl_matrix_set_row failed");
1676 
1677  return;
1678 }
unsigned int numRowsLocal() const
Returns the local row dimension of this matrix.
Definition: GslMatrix.C:275
void resetLU()
In this function resets the LU decomposition of this matrix, as well as deletes the private member po...
Definition: GslMatrix.C:232
MonteCarloSGOptions::MonteCarloSGOptions(const BaseEnvironment &env, const char *prefix, const McOptionsValues &alternativeOptionsValues queso_require_equal_to_msg)(m_env.optionsInputFileName(), std::string(""), std::string("this constructor is incompatible with the existence of an options input file"))
gsl_matrix * m_mat
GSL matrix, also referred to as this matrix.
Definition: GslMatrix.h:409
unsigned int numCols() const
Returns the column dimension of this matrix.
Definition: GslMatrix.C:287
void QUESO::GslMatrix::smallestEigen ( double &  eigenValue,
GslVector eigenVector 
) const

This function finds smallest eigenvalue, namely eigenValue, of this matrix and its corresponding eigenvector, namely eigenVector.

Definition at line 1515 of file GslMatrix.C.

References QUESO::GslVector::abs(), QUESO::Matrix::m_env, QUESO::Matrix::m_map, QUESO::queso_require_not_equal_to_msg, and QUESO::GslVector::sizeLocal().

1516 {
1517  // Sanity Checks
1518  unsigned int n = eigenVector.sizeLocal();
1519 
1520  queso_require_not_equal_to_msg(n, 0, "invalid input vector size");
1521 
1522  /* The following notation is used:
1523  z = vector used in iteration that ends up being the eigenvector corresponding to the
1524  largest eigenvalue
1525  w = vector used in iteration that we extract the largest eigenvalue from. */
1526 
1527  // Some parameters associated with the algorithm
1528  // TODO: Do we want to add the ability to have these set by the user?
1529  const unsigned int max_num_iterations = 1000;
1530  const double tolerance = 1.0e-13;
1531 
1532  // Create temporary working vectors.
1533  // TODO: We shouldn't have to use these - we ought to be able to work directly
1534  // TODO: with eigenValue and eigenVector. Optimize this once we have regression
1535  // TODO: tests.
1536  GslVector z(m_env, m_map, 1.0 ); // Needs to be initialized to 1.0
1537  GslVector w(m_env, m_map);
1538 
1539  // Some variables we use inside the loop.
1540  int index;
1541  double residual;
1542  double one_over_lambda;
1543  double lambda;
1544 
1545  for( unsigned int k = 0; k < max_num_iterations; ++k )
1546  {
1547  w = (*this).invertMultiplyForceLU( z );
1548 
1549  // For this algorithm, it's crucial to get the maximum in
1550  // absolute value, but then to normalize by the actual value
1551  // and *not* the absolute value.
1552  index = (w.abs()).getMaxValueIndex();
1553 
1554  // Remember: Inverse power method solves for max 1/lambda ==> lambda smallest
1555  one_over_lambda = w[index];
1556 
1557  lambda = 1.0/one_over_lambda;
1558 
1559  z = lambda * w;
1560 
1561  // Here we use the norm of the residual as our convergence check:
1562  // norm( A*x - \lambda*x )
1563  residual = ( (*this)*z - lambda*z ).norm2();
1564 
1565  if( residual < tolerance )
1566  {
1567  eigenValue = lambda;
1568 
1569  // TODO: Do we want to normalize this so eigenVector.norm2() = 1?
1570  eigenVector = z;
1571  return;
1572  }
1573 
1574  }
1575 
1576  // If we reach this point, then we didn't converge. Print error message
1577  // to this effect.
1578  // TODO: We know we failed at this point - need a way to not need a test
1579  // TODO: Just specify the error.
1580  queso_require_less_msg(residual, tolerance, "Maximum num iterations exceeded");
1581 
1582  return;
1583 }
const Map m_map
Mapping variable.
Definition: Matrix.h:121
MonteCarloSGOptions::MonteCarloSGOptions(const BaseEnvironment &env, const char *prefix queso_require_not_equal_to_msg)(m_env.optionsInputFileName(), std::string(""), std::string("this constructor is incompatible with the absence of an options input file"))
const BaseEnvironment & m_env
QUESO environment variable.
Definition: Matrix.h:116
void QUESO::GslMatrix::subReadContents ( const std::string &  fileName,
const std::string &  fileType,
const std::set< unsigned int > &  allowedSubEnvIds 
)

Read contents of subenvironment from file fileName.

Definition at line 1841 of file GslMatrix.C.

References QUESO::BaseEnvironment::closeFile(), QUESO::BaseEnvironment::fullRank(), QUESO::FilePtrSetStruct::ifsVar, QUESO::Matrix::m_env, numCols(), QUESO::Matrix::numOfProcsForStorage(), numRowsLocal(), QUESO::BaseEnvironment::openInputFile(), QUESO::queso_require_equal_to_msg, QUESO::BaseEnvironment::subDisplayFile(), and QUESO::BaseEnvironment::subRank().

1845 {
1846  queso_require_greater_equal_msg(m_env.subRank(), 0, "unexpected subRank");
1847 
1848  queso_require_less_equal_msg(this->numOfProcsForStorage(), 1, "implemented just for sequential vectors for now");
1849 
1850  FilePtrSetStruct filePtrSet;
1851  if (m_env.openInputFile(fileName,
1852  fileType, // "m or hdf"
1853  allowedSubEnvIds,
1854  filePtrSet)) {
1855 
1856  // palms
1857  unsigned int nRowsLocal = this->numRowsLocal();
1858 
1859  // In the logic below, the id of a line' begins with value 0 (zero)
1860  unsigned int idOfMyFirstLine = 1;
1861  unsigned int idOfMyLastLine = nRowsLocal;
1862  unsigned int nCols = this->numCols();
1863 
1864  // Read number of matrix rows in the file by taking care of the first line,
1865  // which resembles something like 'variable_name = zeros(n_rows,n_cols);'
1866  std::string tmpString;
1867 
1868  // Read 'variable name' string
1869  *filePtrSet.ifsVar >> tmpString;
1870  //std::cout << "Just read '" << tmpString << "'" << std::endl;
1871 
1872  // Read '=' sign
1873  *filePtrSet.ifsVar >> tmpString;
1874  //std::cout << "Just read '" << tmpString << "'" << std::endl;
1875  queso_require_equal_to_msg(tmpString, std::string("="), std::string("string should be the '=' sign"));
1876 
1877  // Read 'zeros(n_rows,n_cols)' string
1878  *filePtrSet.ifsVar >> tmpString;
1879  //std::cout << "Just read '" << tmpString << "'" << std::endl;
1880  unsigned int posInTmpString = 6;
1881 
1882  // Isolate 'n_rows' in a string
1883  char nRowsString[tmpString.size()-posInTmpString+1];
1884  unsigned int posInRowsString = 0;
1885  do {
1886  queso_require_less_msg(posInTmpString, tmpString.size(), "symbol ',' not found in first line of file");
1887  nRowsString[posInRowsString++] = tmpString[posInTmpString++];
1888  } while (tmpString[posInTmpString] != ',');
1889  nRowsString[posInRowsString] = '\0';
1890 
1891  // Isolate 'n_cols' in a string
1892  posInTmpString++; // Avoid reading ',' char
1893  char nColsString[tmpString.size()-posInTmpString+1];
1894  unsigned int posInColsString = 0;
1895  do {
1896  queso_require_less_msg(posInTmpString, tmpString.size(), "symbol ')' not found in first line of file");
1897  nColsString[posInColsString++] = tmpString[posInTmpString++];
1898  } while (tmpString[posInTmpString] != ')');
1899  nColsString[posInColsString] = '\0';
1900 
1901  // Convert 'n_rows' and 'n_cols' strings to numbers
1902  unsigned int numRowsInFile = (unsigned int) strtod(nRowsString,NULL);
1903  unsigned int numColsInFile = (unsigned int) strtod(nColsString,NULL);
1904  if (m_env.subDisplayFile()) {
1905  *m_env.subDisplayFile() << "In GslMatrix::subReadContents()"
1906  << ": fullRank " << m_env.fullRank()
1907  << ", numRowsInFile = " << numRowsInFile
1908  << ", numColsInFile = " << numColsInFile
1909  << ", nRowsLocal = " << nRowsLocal
1910  << ", nCols = " << nCols
1911  << std::endl;
1912  }
1913 
1914  // Check if [num of rows in file] == [requested matrix row size]
1915  queso_require_equal_to_msg(numRowsInFile, nRowsLocal, "size of vec in file is not big enough");
1916 
1917  // Check if [num of cols in file] == [num cols in current matrix]
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");
1919 
1920  // Code common to any core in a communicator
1921  unsigned int maxCharsPerLine = 64*nCols; // Up to about 60 characters to represent each parameter value
1922 
1923  unsigned int lineId = 0;
1924  while (lineId < idOfMyFirstLine) {
1925  filePtrSet.ifsVar->ignore(maxCharsPerLine,'\n');
1926  lineId++;
1927  };
1928 
1929  if (m_env.subDisplayFile()) {
1930  *m_env.subDisplayFile() << "In GslMatrix::subReadContents()"
1931  << ": beginning to read input actual data"
1932  << std::endl;
1933  }
1934 
1935  // Take care of initial part of the first data line,
1936  // which resembles something like 'variable_name = [value1 value2 ...'
1937  // Read 'variable name' string
1938  *filePtrSet.ifsVar >> tmpString;
1939  //std::cout << "Core 0 just read '" << tmpString << "'" << std::endl;
1940 
1941  // Read '=' sign
1942  *filePtrSet.ifsVar >> tmpString;
1943  //std::cout << "Core 0 just read '" << tmpString << "'" << std::endl;
1944  queso_require_equal_to_msg(tmpString, std::string("="), std::string("in core 0, string should be the '=' sign"));
1945 
1946  // Take into account the ' [' portion
1947  std::streampos tmpPos = filePtrSet.ifsVar->tellg();
1948  filePtrSet.ifsVar->seekg(tmpPos+(std::streampos)2);
1949 
1950  if (m_env.subDisplayFile()) {
1951  *m_env.subDisplayFile() << "In GslMatrix::subReadContents()"
1952  << ": beginning to read lines with numbers only"
1953  << ", lineId = " << lineId
1954  << ", idOfMyFirstLine = " << idOfMyFirstLine
1955  << ", idOfMyLastLine = " << idOfMyLastLine
1956  << std::endl;
1957  }
1958 
1959  double tmpRead;
1960  while (lineId <= idOfMyLastLine) {
1961  for (unsigned int j = 0; j < nCols; ++j) {
1962  *filePtrSet.ifsVar >> tmpRead;
1963  (*this)(lineId-idOfMyFirstLine,j) = tmpRead;
1964  }
1965  lineId++;
1966  };
1967 
1968  m_env.closeFile(filePtrSet,fileType);
1969  }
1970 
1971  return;
1972 }
unsigned int numRowsLocal() const
Returns the local row dimension of this matrix.
Definition: GslMatrix.C:275
unsigned int numOfProcsForStorage() const
Definition: Matrix.C:62
void closeFile(FilePtrSetStruct &filePtrSet, const std::string &fileType) const
Closes the file.
Definition: Environment.C:1084
int subRank() const
Returns the rank of the MPI process in the sub-communicator subComm()
Definition: Environment.C:287
int fullRank() const
Returns the rank of the MPI process in QUESO&#39;s full communicator.
Definition: Environment.C:268
bool openInputFile(const std::string &fileName, const std::string &fileType, const std::set< unsigned int > &allowedSubEnvIds, FilePtrSetStruct &filePtrSet) const
Opens an input file.
Definition: Environment.C:896
MonteCarloSGOptions::MonteCarloSGOptions(const BaseEnvironment &env, const char *prefix, const McOptionsValues &alternativeOptionsValues queso_require_equal_to_msg)(m_env.optionsInputFileName(), std::string(""), std::string("this constructor is incompatible with the existence of an options input file"))
unsigned int numCols() const
Returns the column dimension of this matrix.
Definition: GslMatrix.C:287
const BaseEnvironment & m_env
QUESO environment variable.
Definition: Matrix.h:116
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:320
void QUESO::GslMatrix::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.

Definition at line 1801 of file GslMatrix.C.

References QUESO::BaseEnvironment::closeFile(), QUESO::Matrix::m_env, numCols(), QUESO::Matrix::numOfProcsForStorage(), numRowsLocal(), QUESO::FilePtrSetStruct::ofsVar, QUESO::BaseEnvironment::openOutputFile(), QUESO::BaseEnvironment::subIdString(), and QUESO::BaseEnvironment::subRank().

1806 {
1807  queso_require_greater_equal_msg(m_env.subRank(), 0, "unexpected subRank");
1808 
1809  queso_require_less_equal_msg(this->numOfProcsForStorage(), 1, "implemented just for sequential vectors for now");
1810 
1811  FilePtrSetStruct filePtrSet;
1812  if (m_env.openOutputFile(fileName,
1813  fileType, // "m or hdf"
1814  allowedSubEnvIds,
1815  false,
1816  filePtrSet)) {
1817  unsigned int nRows = this->numRowsLocal();
1818  unsigned int nCols = this->numCols();
1819  *filePtrSet.ofsVar << varNamePrefix << "_sub" << m_env.subIdString() << " = zeros(" << nRows
1820  << "," << nCols
1821  << ");"
1822  << std::endl;
1823  *filePtrSet.ofsVar << varNamePrefix << "_sub" << m_env.subIdString() << " = [";
1824 
1825  for (unsigned int i = 0; i < nRows; ++i) {
1826  for (unsigned int j = 0; j < nCols; ++j) {
1827  *filePtrSet.ofsVar << (*this)(i,j)
1828  << " ";
1829  }
1830  *filePtrSet.ofsVar << "\n";
1831  }
1832  *filePtrSet.ofsVar << "];\n";
1833 
1834  m_env.closeFile(filePtrSet,fileType);
1835  }
1836 
1837  return;
1838 }
unsigned int numRowsLocal() const
Returns the local row dimension of this matrix.
Definition: GslMatrix.C:275
unsigned int numOfProcsForStorage() const
Definition: Matrix.C:62
void closeFile(FilePtrSetStruct &filePtrSet, const std::string &fileType) const
Closes the file.
Definition: Environment.C:1084
int subRank() const
Returns the rank of the MPI process in the sub-communicator subComm()
Definition: Environment.C:287
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.
Definition: Environment.C:348
unsigned int numCols() const
Returns the column dimension of this matrix.
Definition: GslMatrix.C:287
const BaseEnvironment & m_env
QUESO environment variable.
Definition: Matrix.h:116
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...
Definition: Environment.C:521
int QUESO::GslMatrix::svd ( GslMatrix matU,
GslVector vecS,
GslMatrix matVt 
) const

Checks for the dimension of this matrix, matU, VecS and matVt, and calls the protected routine internalSvd to compute the singular values of this.

Definition at line 435 of file GslMatrix.C.

References internalSvd(), m_svdSvec, m_svdUmat, m_svdVTmat, numCols(), numRowsLocal(), QUESO::queso_require_equal_to_msg, and QUESO::GslVector::sizeLocal().

436 {
437  unsigned int nRows = this->numRowsLocal();
438  unsigned int nCols = this->numCols();
439 
440  queso_require_msg(!((matU.numRowsLocal() != nRows) || (matU.numCols() != nCols)), "invalid matU");
441 
442  queso_require_equal_to_msg(vecS.sizeLocal(), nCols, "invalid vecS");
443 
444  queso_require_msg(!((matVt.numRowsLocal() != nCols) || (matVt.numCols() != nCols)), "invalid matVt");
445 
446  int iRC = internalSvd();
447 
448  matU = *m_svdUmat;
449  vecS = *m_svdSvec;
450  matVt = *m_svdVTmat;
451 
452  return iRC;
453 }
unsigned int numRowsLocal() const
Returns the local row dimension of this matrix.
Definition: GslMatrix.C:275
GslMatrix * m_svdUmat
m_svdUmat stores the M-by-N orthogonal matrix U after the singular value decomposition of a matrix...
Definition: GslMatrix.h:424
GslMatrix * m_svdVTmat
m_svdVmatT stores the transpose of N-by-N orthogonal square matrix V, namely V^T, after the singular ...
Definition: GslMatrix.h:443
MonteCarloSGOptions::MonteCarloSGOptions(const BaseEnvironment &env, const char *prefix, const McOptionsValues &alternativeOptionsValues queso_require_equal_to_msg)(m_env.optionsInputFileName(), std::string(""), std::string("this constructor is incompatible with the existence of an options input file"))
GslVector * m_svdSvec
m_svdSvec stores the diagonal of the N-by-N diagonal matrix of singular values S after the singular v...
Definition: GslMatrix.h:431
unsigned int numCols() const
Returns the column dimension of this matrix.
Definition: GslMatrix.C:287
int internalSvd() const
This function factorizes the M-by-N matrix A into the singular value decomposition A = U S V^T for M ...
Definition: GslMatrix.C:532
const GslMatrix & QUESO::GslMatrix::svdMatU ( ) const

This function calls private member GslMatrix::internalSvd() to set a M-by-N orthogonal matrix U of the singular value decomposition (svd) of a general rectangular M-by-N matrix A.

A general rectangular M-by-N matrix A has a singular value decomposition (svd) into the product of an M-by-N orthogonal matrix U, an N-by-N diagonal matrix of singular values S and the transpose of an N-by-N orthogonal square matrix V, A = U S V^T.

Definition at line 512 of file GslMatrix.C.

References internalSvd(), and m_svdUmat.

513 {
514  int iRC = 0;
515  iRC = internalSvd();
516  if (iRC) {}; // just to remove compiler warning
517 
518  return *m_svdUmat;
519 }
GslMatrix * m_svdUmat
m_svdUmat stores the M-by-N orthogonal matrix U after the singular value decomposition of a matrix...
Definition: GslMatrix.h:424
int internalSvd() const
This function factorizes the M-by-N matrix A into the singular value decomposition A = U S V^T for M ...
Definition: GslMatrix.C:532
const GslMatrix & QUESO::GslMatrix::svdMatV ( ) const

This function calls private member GslMatrix::internalSvd() to set a N-by-N orthogonal square matrix V of the singular value decomposition (svd) of a general rectangular M-by-N matrix A.

A general rectangular M-by-N matrix A has a singular value decomposition (svd) into the product of an M-by-N orthogonal matrix U, an N-by-N diagonal matrix of singular values S and the transpose of an N-by-N orthogonal square matrix V, A = U S V^T.

Definition at line 522 of file GslMatrix.C.

References internalSvd(), and m_svdVmat.

523 {
524  int iRC = 0;
525  iRC = internalSvd();
526  if (iRC) {}; // just to remove compiler warning
527 
528  return *m_svdVmat;
529 }
GslMatrix * m_svdVmat
m_svdVmat stores the N-by-N orthogonal square matrix V after the singular value decomposition of a ma...
Definition: GslMatrix.h:437
int internalSvd() const
This function factorizes the M-by-N matrix A into the singular value decomposition A = U S V^T for M ...
Definition: GslMatrix.C:532
int QUESO::GslMatrix::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).

Definition at line 456 of file GslMatrix.C.

References QUESO::GslVector::data(), data(), QUESO::BaseEnvironment::displayVerbosity(), internalSvd(), QUESO::Matrix::m_env, m_svdSvec, m_svdUmat, m_svdVmat, numCols(), numRowsLocal(), QUESO::queso_require_equal_to_msg, QUESO::GslVector::sizeLocal(), and QUESO::BaseEnvironment::subDisplayFile().

Referenced by svdSolve().

457 {
458  unsigned int nRows = this->numRowsLocal();
459  unsigned int nCols = this->numCols();
460 
461  queso_require_equal_to_msg(rhsVec.sizeLocal(), nRows, "invalid rhsVec");
462 
463  queso_require_equal_to_msg(solVec.sizeLocal(), nCols, "invalid solVec");
464 
465  int iRC = internalSvd();
466 
467  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 5)) {
468  *m_env.subDisplayFile() << "In GslMatrix::svdSolve():"
469  << "\n this->numRowsLocal() = " << this->numRowsLocal()
470  << ", this->numCols() = " << this->numCols()
471  << "\n m_svdUmat->numRowsLocal() = " << m_svdUmat->numRowsLocal()
472  << ", m_svdUmat->numCols() = " << m_svdUmat->numCols()
473  << "\n m_svdVmat->numRowsLocal() = " << m_svdVmat->numRowsLocal()
474  << ", m_svdVmat->numCols() = " << m_svdVmat->numCols()
475  << "\n m_svdSvec->sizeLocal() = " << m_svdSvec->sizeLocal()
476  << "\n rhsVec.sizeLocal() = " << rhsVec.sizeLocal()
477  << "\n solVec.sizeLocal() = " << solVec.sizeLocal()
478  << std::endl;
479  }
480 
481  if (iRC == 0) iRC = gsl_linalg_SV_solve(m_svdUmat->data(), m_svdVmat->data(), m_svdSvec->data(), rhsVec.data(), solVec.data());
482 
483  return iRC;
484 }
GslMatrix * m_svdVmat
m_svdVmat stores the N-by-N orthogonal square matrix V after the singular value decomposition of a ma...
Definition: GslMatrix.h:437
unsigned int numRowsLocal() const
Returns the local row dimension of this matrix.
Definition: GslMatrix.C:275
GslMatrix * m_svdUmat
m_svdUmat stores the M-by-N orthogonal matrix U after the singular value decomposition of a matrix...
Definition: GslMatrix.h:424
gsl_matrix * data()
Returns this matrix.
Definition: GslMatrix.C:1766
gsl_vector * data() const
Definition: GslVector.C:970
MonteCarloSGOptions::MonteCarloSGOptions(const BaseEnvironment &env, const char *prefix, const McOptionsValues &alternativeOptionsValues queso_require_equal_to_msg)(m_env.optionsInputFileName(), std::string(""), std::string("this constructor is incompatible with the existence of an options input file"))
GslVector * m_svdSvec
m_svdSvec stores the diagonal of the N-by-N diagonal matrix of singular values S after the singular v...
Definition: GslMatrix.h:431
unsigned int sizeLocal() const
Returns the length of this vector.
Definition: GslVector.C:241
unsigned int displayVerbosity() const
Definition: Environment.C:450
unsigned int numCols() const
Returns the column dimension of this matrix.
Definition: GslMatrix.C:287
const BaseEnvironment & m_env
QUESO environment variable.
Definition: Matrix.h:116
int internalSvd() const
This function factorizes the M-by-N matrix A into the singular value decomposition A = U S V^T for M ...
Definition: GslMatrix.C:532
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:320
int QUESO::GslMatrix::svdSolve ( const GslMatrix rhsMat,
GslMatrix solMat 
) 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=solMat, b=rhsMat).

Definition at line 487 of file GslMatrix.C.

References getColumn(), QUESO::Matrix::m_env, QUESO::Matrix::map(), numCols(), numRowsLocal(), QUESO::queso_require_equal_to_msg, setColumn(), and svdSolve().

488 {
489  unsigned int nRows = this->numRowsLocal();
490  unsigned int nCols = this->numCols();
491 
492  queso_require_equal_to_msg(rhsMat.numRowsLocal(), nRows, "invalid rhsMat");
493 
494  queso_require_equal_to_msg(solMat.numRowsLocal(), nCols, "invalid solMat");
495 
496  queso_require_equal_to_msg(rhsMat.numCols(), solMat.numCols(), "rhsMat and solMat are not compatible");
497 
498  GslVector rhsVec(m_env,rhsMat.map());
499  GslVector solVec(m_env,solMat.map());
500  int iRC = 0;
501  for (unsigned int j = 0; j < rhsMat.numCols(); ++j) {
502  rhsVec = rhsMat.getColumn(j);
503  iRC = this->svdSolve(rhsVec, solVec);
504  if (iRC) break;
505  solMat.setColumn(j,solVec);
506  }
507 
508  return iRC;
509 }
unsigned int numRowsLocal() const
Returns the local row dimension of this matrix.
Definition: GslMatrix.C:275
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).
Definition: GslMatrix.C:456
MonteCarloSGOptions::MonteCarloSGOptions(const BaseEnvironment &env, const char *prefix, const McOptionsValues &alternativeOptionsValues queso_require_equal_to_msg)(m_env.optionsInputFileName(), std::string(""), std::string("this constructor is incompatible with the existence of an options input file"))
unsigned int numCols() const
Returns the column dimension of this matrix.
Definition: GslMatrix.C:287
const BaseEnvironment & m_env
QUESO environment variable.
Definition: Matrix.h:116
GslMatrix QUESO::GslMatrix::transpose ( ) const

This function calculated the transpose of this matrix (square).

Definition at line 693 of file GslMatrix.C.

References QUESO::Map::Comm(), QUESO::Matrix::m_env, QUESO::Matrix::map(), numCols(), and numRowsGlobal().

Referenced by internalSvd().

694 {
695  unsigned int nRows = this->numRowsGlobal();
696  unsigned int nCols = this->numCols();
697 
698  const MpiComm & comm = this->map().Comm();
699  Map serial_map(nCols, 0, comm);
700 
701  GslMatrix mat(m_env,serial_map,nRows);
702 
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);
706  }
707  }
708 
709  return mat;
710 }
unsigned int numRowsGlobal() const
Returns the global row dimension of this matrix.
Definition: GslMatrix.C:281
GslMatrix()
Default Constructor.
unsigned int numCols() const
Returns the column dimension of this matrix.
Definition: GslMatrix.C:287
const MpiComm & Comm() const
Access function for MpiComm communicator.
Definition: Map.C:131
const BaseEnvironment & m_env
QUESO environment variable.
Definition: Matrix.h:116
const Map & map() const
Definition: Matrix.C:54
void QUESO::GslMatrix::zeroLower ( bool  includeDiagonal = false)
virtual

This function sets all the entries bellow the main diagonal of this matrix to zero.

If \c includeDiagonal = false, then only the entries bellow the main diagonal are set to zero;

if includeDiagonal = true, then the elements of the matrix diagonal are also set to zero.

Implements QUESO::Matrix.

Definition at line 591 of file GslMatrix.C.

References numCols(), numRowsLocal(), QUESO::queso_require_equal_to_msg, and resetLU().

592 {
593  unsigned int nRows = this->numRowsLocal();
594  unsigned int nCols = this->numCols();
595 
596  queso_require_equal_to_msg(nRows, nCols, "routine works only for square matrices");
597 
598  this->resetLU();
599 
600  if (includeDiagonal) {
601  for (unsigned int i = 0; i < nRows; i++) {
602  for (unsigned int j = 0; j <= i; j++) {
603  (*this)(i,j) = 0.;
604  }
605  }
606  }
607  else {
608  for (unsigned int i = 0; i < nRows; i++) {
609  for (unsigned int j = 0; j < i; j++) {
610  (*this)(i,j) = 0.;
611  }
612  }
613  }
614 
615  return;
616 }
unsigned int numRowsLocal() const
Returns the local row dimension of this matrix.
Definition: GslMatrix.C:275
void resetLU()
In this function resets the LU decomposition of this matrix, as well as deletes the private member po...
Definition: GslMatrix.C:232
MonteCarloSGOptions::MonteCarloSGOptions(const BaseEnvironment &env, const char *prefix, const McOptionsValues &alternativeOptionsValues queso_require_equal_to_msg)(m_env.optionsInputFileName(), std::string(""), std::string("this constructor is incompatible with the existence of an options input file"))
unsigned int numCols() const
Returns the column dimension of this matrix.
Definition: GslMatrix.C:287
void QUESO::GslMatrix::zeroUpper ( bool  includeDiagonal = false)
virtual

This function sets all the entries above the main diagonal of this matrix to zero.

If \c includeDiagonal = false, then only the entries above the main diagonal are set to zero;

if includeDiagonal = true, then the elements of the matrix diagonal are also set to zero.

Implements QUESO::Matrix.

Definition at line 619 of file GslMatrix.C.

References numCols(), numRowsLocal(), QUESO::queso_require_equal_to_msg, and resetLU().

620 {
621  unsigned int nRows = this->numRowsLocal();
622  unsigned int nCols = this->numCols();
623 
624  queso_require_equal_to_msg(nRows, nCols, "routine works only for square matrices");
625 
626  this->resetLU();
627 
628  if (includeDiagonal) {
629  for (unsigned int i = 0; i < nRows; i++) {
630  for (unsigned int j = i; j < nCols; j++) {
631  (*this)(i,j) = 0.;
632  }
633  }
634  }
635  else {
636  for (unsigned int i = 0; i < nRows; i++) {
637  for (unsigned int j = (i+1); j < nCols; j++) {
638  (*this)(i,j) = 0.;
639  }
640  }
641  }
642 
643  return;
644 }
unsigned int numRowsLocal() const
Returns the local row dimension of this matrix.
Definition: GslMatrix.C:275
void resetLU()
In this function resets the LU decomposition of this matrix, as well as deletes the private member po...
Definition: GslMatrix.C:232
MonteCarloSGOptions::MonteCarloSGOptions(const BaseEnvironment &env, const char *prefix, const McOptionsValues &alternativeOptionsValues queso_require_equal_to_msg)(m_env.optionsInputFileName(), std::string(""), std::string("this constructor is incompatible with the existence of an options input file"))
unsigned int numCols() const
Returns the column dimension of this matrix.
Definition: GslMatrix.C:287

Member Data Documentation

double QUESO::GslMatrix::m_determinant
mutableprivate

The determinant of this matrix.

Definition at line 446 of file GslMatrix.h.

Referenced by determinant(), lnDeterminant(), and resetLU().

GslMatrix* QUESO::GslMatrix::m_inverse
mutableprivate

Inverse matrix of this.

Definition at line 415 of file GslMatrix.h.

Referenced by inverse(), and resetLU().

bool QUESO::GslMatrix::m_isSingular
mutableprivate

Indicates whether or not this matrix is singular.

Definition at line 464 of file GslMatrix.h.

Referenced by invertMultiply(), invertMultiplyForceLU(), and resetLU().

double QUESO::GslMatrix::m_lnDeterminant
mutableprivate

The natural logarithm of the determinant of this matrix.

Definition at line 449 of file GslMatrix.h.

Referenced by determinant(), lnDeterminant(), and resetLU().

gsl_matrix* QUESO::GslMatrix::m_LU
mutableprivate

GSL matrix for the LU decomposition of m_mat.

Definition at line 412 of file GslMatrix.h.

Referenced by determinant(), invertMultiply(), invertMultiplyForceLU(), lnDeterminant(), and resetLU().

gsl_matrix* QUESO::GslMatrix::m_mat
private
gsl_permutation* QUESO::GslMatrix::m_permutation
mutableprivate

The permutation matrix of a LU decomposition.

In the LU decomposition PA = LU, the j-th column of the matrix P is given by the k-th column of the identity matrix, where k = p_j the j-th element of the permutation vector. A is the square matrix of interest and P is the permutation matrix.

Definition at line 455 of file GslMatrix.h.

Referenced by invertMultiply(), invertMultiplyForceLU(), and resetLU().

int QUESO::GslMatrix::m_signum
mutableprivate

m_signum stores the sign of the permutation of the LU decomposition PA = LU.

In the LU decomposition PA = LU, where A is the square matrix of interest and P is the permutation matrix, m_signum has the value (-1)^n, where n is the number of interchanges in the permutation.

Definition at line 461 of file GslMatrix.h.

Referenced by determinant(), invertMultiply(), invertMultiplyForceLU(), lnDeterminant(), and resetLU().

Map* QUESO::GslMatrix::m_svdColMap
mutableprivate

Mapping for matrices involved in the singular value decomposition (svd) routine.

Definition at line 418 of file GslMatrix.h.

Referenced by internalSvd(), and resetLU().

GslVector* QUESO::GslMatrix::m_svdSvec
mutableprivate

m_svdSvec stores the diagonal of the N-by-N diagonal matrix of singular values S after the singular value decomposition of a matrix.

A general rectangular M-by-N matrix A has a singular value decomposition (svd) into the product of an M-by-N orthogonal matrix U, an N-by-N diagonal matrix of singular values S and the transpose of an N-by-N orthogonal square matrix V, A = U S V^T.

Definition at line 431 of file GslMatrix.h.

Referenced by internalSvd(), rank(), resetLU(), svd(), and svdSolve().

GslMatrix* QUESO::GslMatrix::m_svdUmat
mutableprivate

m_svdUmat stores the M-by-N orthogonal matrix U after the singular value decomposition of a matrix.

A general rectangular M-by-N matrix A has a singular value decomposition (svd) into the product of an M-by-N orthogonal matrix U, an N-by-N diagonal matrix of singular values S and the transpose of an N-by-N orthogonal square matrix V, A = U S V^T.

Definition at line 424 of file GslMatrix.h.

Referenced by internalSvd(), resetLU(), svd(), svdMatU(), and svdSolve().

GslMatrix* QUESO::GslMatrix::m_svdVmat
mutableprivate

m_svdVmat stores the N-by-N orthogonal square matrix V after the singular value decomposition of a matrix.

A general rectangular M-by-N matrix A has a singular value decomposition (svd) into the product of an M-by-N orthogonal matrix U, an N-by-N diagonal matrix of singular values S and the transpose of an N-by-N orthogonal square matrix V, A = U S V^T.

Definition at line 437 of file GslMatrix.h.

Referenced by internalSvd(), resetLU(), svdMatV(), and svdSolve().

GslMatrix* QUESO::GslMatrix::m_svdVTmat
mutableprivate

m_svdVmatT stores the transpose of N-by-N orthogonal square matrix V, namely V^T, after the singular value decomposition of a matrix.

A general rectangular M-by-N matrix A has a singular value decomposition (svd) into the product of an M-by-N orthogonal matrix U, an N-by-N diagonal matrix of singular values S and the transpose of an N-by-N orthogonal square matrix V, A = U S V^T.

Definition at line 443 of file GslMatrix.h.

Referenced by internalSvd(), resetLU(), and svd().


The documentation for this class was generated from the following files:

Generated on Tue Jun 5 2018 19:49:08 for queso-0.57.1 by  doxygen 1.8.5