QUESO::GslMatrix Class Reference

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

#include <GslMatrix.h>

Inheritance diagram for QUESO::GslMatrix:

Collaboration diagram for QUESO::GslMatrix:

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
GslMatrix &	operator= (const GslMatrix &rhs)
	Copies values from matrix rhs to this. More...
GslMatrix &	operator*= (double a)
	Stores in this the coordinate-wise multiplication of this and a. More...
GslMatrix &	operator/= (double a)
	Stores in this the coordinate-wise division of this by a. More...
GslMatrix &	operator+= (const GslMatrix &rhs)
	Stores in this the coordinate-wise addition of this and rhs. More...
GslMatrix &	operator-= (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 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. More...
const 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. 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...
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 BaseEnvironment &	env () const
const Map &	map () 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...
void	multiply (const GslVector &x, GslVector &y) const
	This function multiplies this matrix by vector x and stores the resulting vector in y. 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...
GslMatrix *	m_inverse
	Inverse matrix of this. More...
Map *	m_svdColMap
	Mapping for matrices involved in the singular value decomposition (svd) routine. More...
GslMatrix *	m_svdUmat
	m_svdUmat stores the M-by-N orthogonal matrix U after the singular value decomposition of a matrix. More...
GslVector *	m_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...
GslMatrix *	m_svdVmat
	m_svdVmat stores the N-by-N orthogonal square matrix V after the singular value decomposition of a matrix. More...
GslMatrix *	m_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 BaseEnvironment &	m_env
	QUESO environment variable. More...
const Map &	m_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, and queso_require_msg.

  :
  Matrix  (env,map),
  m_mat          (gsl_matrix_calloc(map.NumGlobalElements(),nCols)),
  m_LU           (NULL),
  m_inverse      (NULL),
  m_svdColMap    (NULL),
  m_svdUmat      (NULL),
  m_svdSvec      (NULL),
  m_svdVmat      (NULL),
  m_svdVTmat     (NULL),
  m_determinant  (-INFINITY),
  m_lnDeterminant(-INFINITY),
  m_permutation  (NULL),
  m_signum       (0),
  m_isSingular   (false)
{
  queso_require_msg(m_mat, "null matrix generated");
}

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, and queso_require_msg.

  :
  Matrix  (env,map),
  m_mat          (gsl_matrix_calloc(map.NumGlobalElements(),map.NumGlobalElements())),
  m_LU           (NULL),
  m_inverse      (NULL),
  m_svdColMap    (NULL),
  m_svdUmat      (NULL),
  m_svdSvec      (NULL),
  m_svdVmat      (NULL),
  m_svdVTmat     (NULL),
  m_determinant  (-INFINITY),
  m_lnDeterminant(-INFINITY),
  m_permutation  (NULL),
  m_signum       (0),
  m_isSingular   (false)
{
  queso_require_msg(m_mat, "null matrix generated");

  for (unsigned int i = 0; i < m_mat->size1; ++i) {
    (*this)(i,i) = diagValue;
  }
}

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, and queso_require_msg.

  :
  Matrix  (v.env(),v.map()),
  m_mat          (gsl_matrix_calloc(v.sizeLocal(),v.sizeLocal())),
  m_LU           (NULL),
  m_inverse      (NULL),
  m_svdColMap    (NULL),
  m_svdUmat      (NULL),
  m_svdSvec      (NULL),
  m_svdVmat      (NULL),
  m_svdVTmat     (NULL),
  m_determinant  (-INFINITY),
  m_lnDeterminant(-INFINITY),
  m_permutation  (NULL),
  m_signum       (0),
  m_isSingular   (false)
{
  queso_require_msg(m_mat, "null matrix generated");

  for (unsigned int i = 0; i < m_mat->size1; ++i) {
    (*this)(i,i) = diagValue;
  }
}

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, queso_require_msg, and QUESO::GslVector::sizeLocal().

  :
  Matrix  (v.env(),v.map()),
  m_mat          (gsl_matrix_calloc(v.sizeLocal(),v.sizeLocal())),
  m_LU           (NULL),
  m_inverse      (NULL),
  m_svdColMap    (NULL),
  m_svdUmat      (NULL),
  m_svdSvec      (NULL),
  m_svdVmat      (NULL),
  m_svdVTmat     (NULL),
  m_determinant  (-INFINITY),
  m_lnDeterminant(-INFINITY),
  m_permutation  (NULL),
  m_signum       (0),
  m_isSingular   (false)
{
  queso_require_msg(m_mat, "null matrix generated");

  unsigned int dim = v.sizeLocal();
  for (unsigned int i = 0; i < dim; ++i) {
    (*this)(i,i) = v[i];
  }
}

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(), m_mat, and queso_require_msg.

  :
  Matrix  (B.env(),B.map()),
  m_mat          (gsl_matrix_calloc(B.numRowsLocal(),B.numCols())),
  m_LU           (NULL),
  m_inverse      (NULL),
  m_svdColMap    (NULL),
  m_svdUmat      (NULL),
  m_svdSvec      (NULL),
  m_svdVmat      (NULL),
  m_svdVTmat     (NULL),
  m_determinant  (-INFINITY),
  m_lnDeterminant(-INFINITY),
  m_permutation  (NULL),
  m_signum       (0),
  m_isSingular   (false)
{
  queso_require_msg(m_mat, "null vector generated");
  this->Matrix::base_copy(B);
  this->copy(B);
}

QUESO::GslMatrix::~GslMatrix

(

)

Destructor.

Definition at line 161 of file GslMatrix.C.

References m_mat, and resetLU().

{
  this->resetLU();
  if (m_mat) gsl_matrix_free(m_mat);
}

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, UQ_RC_MACRO, and QUESO::BaseEnvironment::worldRank().

{
  int iRC;
  //std::cout << "Calling gsl_linalg_cholesky_decomp()..." << std::endl;
  gsl_error_handler_t* oldHandler;
  oldHandler = gsl_set_error_handler_off();
  iRC = gsl_linalg_cholesky_decomp(m_mat);
  if (iRC != 0) {
    std::cerr << "In GslMatrix::chol()"
              << ": iRC = " << iRC
              << ", gsl error message = " << gsl_strerror(iRC)
              << std::endl;
    std::cerr << "Here is the offending matrix: " << std::endl;
    std::cerr << *this << std::endl;
  }
  gsl_set_error_handler(oldHandler);
  //std::cout << "Returned from gsl_linalg_cholesky_decomp() with iRC = " << iRC << std::endl;
  UQ_RC_MACRO(iRC, // Yes, *not* a fatal check on RC
              m_env.worldRank(),
              "GslMatrix::chol()",
              "matrix is not positive definite",
              UQ_MATRIX_IS_NOT_POS_DEFINITE_RC);

  return iRC;
}

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, queso_require_msg, and resetLU().

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

{
  // FIXME - should we be calling Matrix::base_copy here? - RHS
  this->resetLU();
  int iRC;
  iRC = gsl_matrix_memcpy(this->m_mat, src.m_mat);
  queso_require_msg(!(iRC), "failed");

  return;
}

void QUESO::GslMatrix::cwExtract	(	unsigned int	rowId,
		unsigned int	colId,
		GslMatrix &	mat
	)		const

Definition at line 385 of file GslMatrix.C.

References numCols(), numRowsLocal(), queso_require_less_equal_msg, and queso_require_less_msg.

{
  queso_require_less_msg(initialTargetRowId, this->numRowsLocal(), "invalid initialTargetRowId");

  queso_require_less_equal_msg((initialTargetRowId + mat.numRowsLocal()), this->numRowsLocal(), "invalid vec.numRowsLocal()");

  queso_require_less_msg(initialTargetColId, this->numCols(), "invalid initialTargetColId");

  queso_require_less_equal_msg((initialTargetColId + mat.numCols()), this->numCols(), "invalid vec.numCols()");

  for (unsigned int i = 0; i < mat.numRowsLocal(); ++i) {
    for (unsigned int j = 0; j < mat.numCols(); ++j) {
      mat(i,j) = (*this)(initialTargetRowId+i,initialTargetColId+j) ;
    }
  }

  return;
}

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().

{
  unsigned int nRows = this->numRowsLocal();
  unsigned int nCols = this->numCols();

  for (unsigned int row = 0; row < nRows; ++row) {
    for (unsigned int col = 0; col < nCols; ++col) {
      *gsl_matrix_ptr(m_mat,row,col) = value;
    }
  }

  return;
}

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(), numRowsLocal(), queso_require_less_equal_msg, and queso_require_less_msg.

{
  queso_require_less_msg(initialTargetRowId, this->numRowsLocal(), "invalid initialTargetRowId");

  queso_require_less_equal_msg((initialTargetRowId + mat.numRowsLocal()), this->numRowsLocal(), "invalid vec.numRowsLocal()");

  queso_require_less_msg(initialTargetColId, this->numCols(), "invalid initialTargetColId");

  queso_require_less_equal_msg((initialTargetColId + mat.numCols()), this->numCols(), "invalid vec.numCols()");

  for (unsigned int i = 0; i < mat.numRowsLocal(); ++i) {
    for (unsigned int j = 0; j < mat.numCols(); ++j) {
      (*this)(initialTargetRowId+i,initialTargetColId+j) = mat(i,j);
    }
  }

  return;
}

gsl_matrix * QUESO::GslMatrix::data

(

)

Returns this matrix.

Definition at line 1731 of file GslMatrix.C.

References m_mat.

Referenced by internalSvd(), and svdSolve().

{
  return m_mat;
}

double QUESO::GslMatrix::determinant

(

)

const

Calculates the determinant of this matrix.

Definition at line 1036 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().

{
  if (m_determinant == -INFINITY) {
    if (m_LU == NULL) {
      GslVector tmpB(m_env,m_map);
      GslVector tmpX(m_env,m_map);
      if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 99)) {
        *m_env.subDisplayFile() << "In GslMatrix::determinant()"
                                << ": before 'this->invertMultiply()'"
                                << std::endl;
      }
      this->invertMultiply(tmpB,tmpX);
      if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 99)) {
        *m_env.subDisplayFile() << "In GslMatrix::determinant()"
                                << ": after 'this->invertMultiply()'"
                                << std::endl;
      }
    }
    if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 99)) {
      *m_env.subDisplayFile() << "In GslMatrix::determinant()"
                              << ": before 'gsl_linalg_LU_det()'"
                              << std::endl;
    }
    m_determinant   = gsl_linalg_LU_det(m_LU,m_signum);
    if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 99)) {
      *m_env.subDisplayFile() << "In GslMatrix::determinant()"
                              << ": after 'gsl_linalg_LU_det()'"
                              << std::endl;
    }
    m_lnDeterminant = gsl_linalg_LU_lndet(m_LU);
    if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 99)) {
      *m_env.subDisplayFile() << "In GslMatrix::determinant()"
                              << ": after 'gsl_linalg_LU_lndet()'"
                              << std::endl;
    }
  }

  return m_determinant;
}

void QUESO::GslMatrix::eigen	(	GslVector &	eigenValues,
		GslMatrix *	eigenVectors
	)		const

This function computes the eigenvalues of a real symmetric matrix.

Definition at line 1385 of file GslMatrix.C.

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

{
  unsigned int n = eigenValues.sizeLocal();

  queso_require_not_equal_to_msg(n, 0, "invalid input vector size");

  if (eigenVectors) {
    queso_require_equal_to_msg(eigenValues.sizeLocal(), eigenVectors->numRowsLocal(), "different input vector sizes");
  }

  if (eigenVectors == NULL) {
    gsl_eigen_symm_workspace* w = gsl_eigen_symm_alloc((size_t) n);
    gsl_eigen_symm(m_mat,eigenValues.data(),w);
    gsl_eigen_symm_free(w);
  }
  else {
    gsl_eigen_symmv_workspace* w = gsl_eigen_symmv_alloc((size_t) n);
    gsl_eigen_symmv(m_mat,eigenValues.data(),eigenVectors->m_mat,w);
    gsl_eigen_symmv_sort(eigenValues.data(),eigenVectors->m_mat,GSL_EIGEN_SORT_VAL_ASC);
    gsl_eigen_symmv_free(w);
  }

  return;
}

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 751 of file GslMatrix.C.

References numCols(), numRowsLocal(), queso_require_equal_to_msg, and queso_require_greater_equal_msg.

{
  unsigned int sumNumRowsLocals = 0;
  unsigned int sumNumCols       = 0;
  for (unsigned int i = 0; i < matrices.size(); ++i) {
    sumNumRowsLocals += matrices[i]->numRowsLocal();
    sumNumCols       += matrices[i]->numCols();
  }
  queso_require_greater_equal_msg(this->numRowsLocal(), (initialTargetRowId + sumNumRowsLocals), "too big number of rows");
  if (checkForExactNumRowsMatching) queso_require_equal_to_msg(this->numRowsLocal(), (initialTargetRowId + sumNumRowsLocals), "inconsistent number of rows");
  queso_require_greater_equal_msg(this->numCols(), (initialTargetColId + sumNumCols), "too big number of cols");
  if (checkForExactNumColsMatching) queso_require_equal_to_msg(this->numCols(), (initialTargetColId + sumNumCols), "inconsistent number of cols");

  unsigned int cumulativeRowId = 0;
  unsigned int cumulativeColId = 0;
  for (unsigned int i = 0; i < matrices.size(); ++i) {
    unsigned int nRows = matrices[i]->numRowsLocal();
    unsigned int nCols = matrices[i]->numCols();
    for (unsigned int rowId = 0; rowId < nRows; ++rowId) {
      for (unsigned int colId = 0; colId < nCols; ++colId) {
        (*this)(initialTargetRowId + cumulativeRowId + rowId, initialTargetColId + cumulativeColId + colId) = (*(matrices[i]))(rowId,colId);
      }
    }
    cumulativeRowId += nRows;
    cumulativeColId += nCols;
  }

  return;
}

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 787 of file GslMatrix.C.

References numCols(), numRowsLocal(), queso_require_equal_to_msg, and queso_require_greater_equal_msg.

{
  unsigned int sumNumRowsLocals = 0;
  unsigned int sumNumCols       = 0;
  for (unsigned int i = 0; i < matrices.size(); ++i) {
    sumNumRowsLocals += matrices[i]->numRowsLocal();
    sumNumCols       += matrices[i]->numCols();
  }
  queso_require_greater_equal_msg(this->numRowsLocal(), (initialTargetRowId + sumNumRowsLocals), "too big number of rows");
  if (checkForExactNumRowsMatching) queso_require_equal_to_msg(this->numRowsLocal(), (initialTargetRowId + sumNumRowsLocals), "inconsistent number of rows");
  queso_require_greater_equal_msg(this->numCols(), (initialTargetColId + sumNumCols), "too big number of cols");
  if (checkForExactNumColsMatching) queso_require_equal_to_msg(this->numCols(), (initialTargetColId + sumNumCols), "inconsistent number of cols");

  unsigned int cumulativeRowId = 0;
  unsigned int cumulativeColId = 0;
  for (unsigned int i = 0; i < matrices.size(); ++i) {
    unsigned int nRows = matrices[i]->numRowsLocal();
    unsigned int nCols = matrices[i]->numCols();
    for (unsigned int rowId = 0; rowId < nRows; ++rowId) {
      for (unsigned int colId = 0; colId < nCols; ++colId) {
        (*this)(initialTargetRowId + cumulativeRowId + rowId, initialTargetColId + cumulativeColId + colId) = (*(matrices[i]))(rowId,colId);
      }
    }
    cumulativeRowId += nRows;
    cumulativeColId += nCols;
  }

  return;
}

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 823 of file GslMatrix.C.

References numCols(), numRowsLocal(), queso_require_equal_to_msg, and queso_require_greater_equal_msg.

{
  unsigned int sumNumCols = 0;
  for (unsigned int i = 0; i < matrices.size(); ++i) {
    queso_require_greater_equal_msg(this->numRowsLocal(), (initialTargetRowId + matrices[i]->numRowsLocal()), "too big number of rows");
    if (checkForExactNumRowsMatching) queso_require_equal_to_msg(this->numRowsLocal(), (initialTargetRowId + matrices[i]->numRowsLocal()), "inconsistent number of rows");
    sumNumCols += matrices[i]->numCols();
  }
  queso_require_greater_equal_msg(this->numCols(), (initialTargetColId + sumNumCols), "too big number of cols");
  if (checkForExactNumColsMatching) queso_require_equal_to_msg(this->numCols(), (initialTargetColId + sumNumCols), "inconsistent number of cols");

  unsigned int cumulativeColId = 0;
  for (unsigned int i = 0; i < matrices.size(); ++i) {
    unsigned int nRows = matrices[i]->numRowsLocal();
    unsigned int nCols = matrices[i]->numCols();
    for (unsigned int rowId = 0; rowId < nRows; ++rowId) {
      for (unsigned int colId = 0; colId < nCols; ++colId) {
        (*this)(initialTargetRowId + rowId, initialTargetColId + cumulativeColId + colId) = (*(matrices[i]))(rowId,colId);
      }
    }
    cumulativeColId += nCols;
  }

  return;
}

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 855 of file GslMatrix.C.

References numCols(), numRowsLocal(), queso_require_equal_to_msg, and queso_require_greater_equal_msg.

{
  unsigned int sumNumCols = 0;
  for (unsigned int i = 0; i < matrices.size(); ++i) {
    queso_require_greater_equal_msg(this->numRowsLocal(), (initialTargetRowId + matrices[i]->numRowsLocal()), "too big number of rows");
    if (checkForExactNumRowsMatching) queso_require_equal_to_msg(this->numRowsLocal(), (initialTargetRowId + matrices[i]->numRowsLocal()), "inconsistent number of rows");
    sumNumCols += matrices[i]->numCols();
  }
  queso_require_greater_equal_msg(this->numCols(), (initialTargetColId + sumNumCols), "too big number of cols");
  if (checkForExactNumColsMatching) queso_require_equal_to_msg(this->numCols(), (initialTargetColId + sumNumCols), "inconsistent number of cols");

  unsigned int cumulativeColId = 0;
  for (unsigned int i = 0; i < matrices.size(); ++i) {
    unsigned int nRows = matrices[i]->numRowsLocal();
    unsigned int nCols = matrices[i]->numCols();
    for (unsigned int rowId = 0; rowId < nRows; ++rowId) {
      for (unsigned int colId = 0; colId < nCols; ++colId) {
        (*this)(initialTargetRowId + rowId, initialTargetColId + cumulativeColId + colId) = (*(matrices[i]))(rowId,colId);
      }
    }
    cumulativeColId += nCols;
  }

  return;
}

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 887 of file GslMatrix.C.

References numCols(), numRowsLocal(), queso_require_equal_to_msg, and queso_require_greater_equal_msg.

{
  unsigned int sumNumRows = 0;
  for (unsigned int i = 0; i < matrices.size(); ++i) {
    queso_require_greater_equal_msg(this->numCols(), (initialTargetColId + matrices[i]->numCols()), "too big number of cols");
    if (checkForExactNumColsMatching) queso_require_equal_to_msg(this->numCols(), (initialTargetColId + matrices[i]->numCols()), "inconsistent number of cols");
    sumNumRows += matrices[i]->numRowsLocal();
  }
  queso_require_greater_equal_msg(this->numRowsLocal(), (initialTargetRowId + sumNumRows), "too big number of rows");
  if (checkForExactNumRowsMatching) queso_require_equal_to_msg(this->numRowsLocal(), (initialTargetRowId + sumNumRows), "inconsistent number of rows");

  unsigned int cumulativeRowId = 0;
  for (unsigned int i = 0; i < matrices.size(); ++i) {
    unsigned int nRows = matrices[i]->numRowsLocal();
    unsigned int nCols = matrices[i]->numCols();
    for (unsigned int rowId = 0; rowId < nRows; ++rowId) {
      for (unsigned int colId = 0; colId < nCols; ++colId) {
        (*this)(initialTargetRowId + cumulativeRowId + rowId, initialTargetColId + colId) = (*(matrices[i]))(rowId,colId);
      }
    }
    cumulativeRowId += nRows;
  }

  return;
}

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 919 of file GslMatrix.C.

References numCols(), numRowsLocal(), queso_require_equal_to_msg, and queso_require_greater_equal_msg.

{
  unsigned int sumNumRows = 0;
  for (unsigned int i = 0; i < matrices.size(); ++i) {
    queso_require_greater_equal_msg(this->numCols(), (initialTargetColId + matrices[i]->numCols()), "too big number of cols");
    if (checkForExactNumColsMatching) queso_require_equal_to_msg(this->numCols(), (initialTargetColId + matrices[i]->numCols()), "inconsistent number of cols");
    sumNumRows += matrices[i]->numRowsLocal();
  }
  queso_require_greater_equal_msg(this->numRowsLocal(), (initialTargetRowId + sumNumRows), "too big number of rows");
  if (checkForExactNumRowsMatching) queso_require_equal_to_msg(this->numRowsLocal(), (initialTargetRowId + sumNumRows), "inconsistent number of rows");

  unsigned int cumulativeRowId = 0;
  for (unsigned int i = 0; i < matrices.size(); ++i) {
    unsigned int nRows = matrices[i]->numRowsLocal();
    unsigned int nCols = matrices[i]->numCols();
    for (unsigned int rowId = 0; rowId < nRows; ++rowId) {
      for (unsigned int colId = 0; colId < nCols; ++colId) {
        (*this)(initialTargetRowId + cumulativeRowId + rowId, initialTargetColId + colId) = (*(matrices[i]))(rowId,colId);
      }
    }
    cumulativeRowId += nRows;
  }

  return;
}

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 951 of file GslMatrix.C.

References numCols(), numRowsLocal(), queso_require_equal_to_msg, and queso_require_greater_equal_msg.

{
  queso_require_greater_equal_msg(this->numRowsLocal(), (initialTargetRowId + (mat1.numRowsLocal() * mat2.numRowsLocal())), "too big number of rows");
  if (checkForExactNumRowsMatching) queso_require_equal_to_msg(this->numRowsLocal(), (initialTargetRowId + (mat1.numRowsLocal() * mat2.numRowsLocal())), "inconsistent number of rows");
  queso_require_greater_equal_msg(this->numCols(), (initialTargetColId + (mat1.numCols() * mat2.numCols())), "too big number of columns");
  if (checkForExactNumColsMatching) queso_require_equal_to_msg(this->numCols(), (initialTargetColId + (mat1.numCols() * mat2.numCols())), "inconsistent number of columns");

  for (unsigned int rowId1 = 0; rowId1 < mat1.numRowsLocal(); ++rowId1) {
    for (unsigned int colId1 = 0; colId1 < mat1.numCols(); ++colId1) {
      double multiplicativeFactor = mat1(rowId1,colId1);
      unsigned int targetRowId = rowId1 * mat2.numRowsLocal();
      unsigned int targetColId = colId1 * mat2.numCols();
      for (unsigned int rowId2 = 0; rowId2 < mat2.numRowsLocal(); ++rowId2) {
        for (unsigned int colId2 = 0; colId2 < mat2.numCols(); ++colId2) {
          (*this)(initialTargetRowId + targetRowId + rowId2, initialTargetColId + targetColId + colId2) = multiplicativeFactor * mat2(rowId2,colId2);
        }
      }
    }
  }

  return;
}

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 981 of file GslMatrix.C.

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

{
  queso_require_greater_equal_msg(this->numRowsLocal(), (initialTargetRowId + (mat1.numRowsLocal() * vec2.sizeLocal())), "too big number of rows");
  if (checkForExactNumRowsMatching) queso_require_equal_to_msg(this->numRowsLocal(), (initialTargetRowId + (mat1.numRowsLocal() * vec2.sizeLocal())), "inconsistent number of rows");
  queso_require_greater_equal_msg(this->numCols(), (initialTargetColId + (mat1.numCols() * 1)), "too big number of columns");
  if (checkForExactNumColsMatching) queso_require_equal_to_msg(this->numCols(), (initialTargetColId + (mat1.numCols() * 1)), "inconsistent number of columns");

  for (unsigned int rowId1 = 0; rowId1 < mat1.numRowsLocal(); ++rowId1) {
    for (unsigned int colId1 = 0; colId1 < mat1.numCols(); ++colId1) {
      double multiplicativeFactor = mat1(rowId1,colId1);
      unsigned int targetRowId = rowId1 * vec2.sizeLocal();
      unsigned int targetColId = colId1 * 1;
      for (unsigned int rowId2 = 0; rowId2 < vec2.sizeLocal(); ++rowId2) {
        for (unsigned int colId2 = 0; colId2 < 1; ++colId2) {
          (*this)(initialTargetRowId + targetRowId + rowId2, initialTargetColId + targetColId + colId2) = multiplicativeFactor * vec2[rowId2];
        }
      }
    }
  }


  return;
}

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 1012 of file GslMatrix.C.

References numCols(), numRowsLocal(), queso_require_equal_to_msg, and queso_require_greater_equal_msg.

{
  unsigned int nRows = mat.numRowsLocal();
  unsigned int nCols = mat.numCols();
  queso_require_greater_equal_msg(this->numRowsLocal(), (initialTargetRowId + nCols), "too big number of rows");
  if (checkForExactNumRowsMatching) queso_require_equal_to_msg(this->numRowsLocal(), (initialTargetRowId + nCols), "inconsistent number of rows");
  queso_require_greater_equal_msg(this->numCols(), (initialTargetColId + nRows), "too big number of cols");
  if (checkForExactNumColsMatching) queso_require_equal_to_msg(this->numCols(), (initialTargetColId + nRows), "inconsistent number of cols");

  for (unsigned int row = 0; row < nRows; ++row) {
    for (unsigned int col = 0; col < nCols; ++col) {
      (*this)(initialTargetRowId + col, initialTargetColId + row) = mat(row,col);
    }
  }

  return;
}

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().

{
  unsigned int nRows = this->numRowsLocal();
  unsigned int nCols = this->numCols();
  for (unsigned int i = 0; i < nRows; ++i) {
    for (unsigned int j = 0; j < nCols; ++j) {
      double aux = (*this)(i,j);
      // If 'thresholdValue' is negative, no values will be filtered
      if ((aux             < 0. ) &&
          (-thresholdValue > aux)) {
        (*this)(i,j) = 0.;
      }
      if ((aux            > 0. ) &&
          (thresholdValue < aux)) {
        (*this)(i,j) = 0.;
      }
    }
  }

  return;
}

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().

{
  unsigned int nRows = this->numRowsLocal();
  unsigned int nCols = this->numCols();
  for (unsigned int i = 0; i < nRows; ++i) {
    for (unsigned int j = 0; j < nCols; ++j) {
      double aux = (*this)(i,j);
      // If 'thresholdValue' is negative, no values will be filtered
      if ((aux             < 0. ) &&
          (-thresholdValue < aux)) {
        (*this)(i,j) = 0.;
      }
      if ((aux            > 0. ) &&
          (thresholdValue > aux)) {
        (*this)(i,j) = 0.;
      }
    }
  }

  return;
}

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 1551 of file GslMatrix.C.

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

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

{
  // Sanity checks
  queso_require_less_msg(column_num, this->numCols(), "Specified row number not within range");

  queso_require_equal_to_msg(column.sizeLocal(), this->numRowsLocal(), "column vector not same size as this matrix");

  // Temporary working vector
  gsl_vector* gsl_column = gsl_vector_alloc( column.sizeLocal() );

  int error_code = gsl_matrix_get_col( gsl_column, m_mat, column_num );
  queso_require_equal_to_msg(error_code, 0, "gsl_matrix_get_col failed");

  // Copy column from gsl matrix into our GslVector object
  for( unsigned int i = 0; i < column.sizeLocal(); ++i )
    {
      column[i] = gsl_vector_get( gsl_column, i );
    }

  // Clean up gsl temporaries
  gsl_vector_free( gsl_column );

  return;
}

GslVector QUESO::GslMatrix::getColumn

(

const unsigned int

column_num

)

const

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

Definition at line 1616 of file GslMatrix.C.

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

{
  // We rely on the void getRow's to do sanity checks
  // So we don't do them here.

  GslVector column(m_env, m_map);

  this->getColumn( column_num, column );

  return column;
}

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 1577 of file GslMatrix.C.

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

Referenced by getRow().

{
  // Sanity checks
  queso_require_less_msg(row_num, this->numRowsLocal(), "Specified row number not within range");

  queso_require_equal_to_msg(row.sizeLocal(), this->numCols(), "row vector not same size as this matrix");

  // Temporary working vector
  gsl_vector* gsl_row = gsl_vector_alloc( row.sizeLocal() );

  int error_code = gsl_matrix_get_row( gsl_row, m_mat, row_num );
  queso_require_equal_to_msg(error_code, 0, "gsl_matrix_get_row failed");

  // Copy row from gsl matrix into our GslVector object
  for( unsigned int i = 0; i < row.sizeLocal(); ++i )
    {
      row[i] = gsl_vector_get( gsl_row, i );
    }

  // Clean up gsl temporaries
  gsl_vector_free( gsl_row );

  return;
}

GslVector QUESO::GslMatrix::getRow

(

const unsigned int

row_num

)

const

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

Definition at line 1603 of file GslMatrix.C.

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

{
  // We rely on the void getRow's to do sanity checks
  // So we don't do them here.

  GslVector row(m_env, m_map);

  this->getRow( row_num, row );

  return row;
}

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(), queso_require_greater_equal_msg, transpose(), QUESO::UQ_MATRIX_SVD_FAILED_RC, UQ_RC_MACRO, and QUESO::BaseEnvironment::worldRank().

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

{
  int iRC = 0;

  if (m_svdColMap == NULL) {
    unsigned int nRows = this->numRowsLocal();
    unsigned int nCols = this->numCols();
    queso_require_greater_equal_msg(nRows, nCols, "GSL only supports cases where nRows >= nCols");

    m_svdColMap = new Map(this->numCols(),0,this->map().Comm()); // see 'VectorSpace<.,.>::newMap()' in src/basic/src/GslVectorSpace.C
    m_svdUmat   = new GslMatrix(*this); // Yes, 'this'
    m_svdSvec   = new GslVector(m_env,*m_svdColMap);
    m_svdVmat   = new GslMatrix(*m_svdSvec);
    m_svdVTmat  = new GslMatrix(*m_svdSvec);

    //std::cout << "In GslMatrix::internalSvd()"
    //          << ", calling gsl_linalg_SV_decomp_jacobi()..."
    //          << ": nRows = " << nRows
    //          << ", nCols = " << nCols
    //          << std::endl;
    struct timeval timevalBegin;
    gettimeofday(&timevalBegin, NULL);
    gsl_error_handler_t* oldHandler;
    oldHandler = gsl_set_error_handler_off();
#if 1
    iRC = gsl_linalg_SV_decomp_jacobi(m_svdUmat->data(), m_svdVmat->data(), m_svdSvec->data());
#else
    GslVector vecWork(*m_svdSvec );
    iRC = gsl_linalg_SV_decomp(m_svdUmat->data(), m_svdVmat->data(), m_svdSvec->data(), vecWork.data());
#endif
    if (iRC != 0) {
      std::cerr << "In GslMatrix::internalSvd()"
                << ": iRC = " << iRC
                << ", gsl error message = " << gsl_strerror(iRC)
                << std::endl;
    }
    gsl_set_error_handler(oldHandler);

    struct timeval timevalNow;
    gettimeofday(&timevalNow, NULL);
    //std::cout << "In GslMatrix::internalSvd()"
    //          << ": returned from gsl_linalg_SV_decomp_jacobi() with iRC = " << iRC
    //          << " after " << timevalNow.tv_sec - timevalBegin.tv_sec
    //          << " seconds"
    //          << std::endl;
    UQ_RC_MACRO(iRC, // Yes, *not* a fatal check on RC
                m_env.worldRank(),
                "GslMatrix::internalSvd()",
                "matrix svd failed",
                UQ_MATRIX_SVD_FAILED_RC);
    *m_svdVTmat = m_svdVmat->transpose();
  }

  return iRC;
}

GslMatrix QUESO::GslMatrix::inverse

(

)

const

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

Definition at line 711 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_require_equal_to_msg, and QUESO::BaseEnvironment::subDisplayFile().

{
  unsigned int nRows = this->numRowsLocal();
  unsigned int nCols = this->numCols();

  queso_require_equal_to_msg(nRows, nCols, "matrix is not square");

  if (m_inverse == NULL) {
    m_inverse = new GslMatrix(m_env,m_map,nCols);
    GslVector unitVector(m_env,m_map);
    unitVector.cwSet(0.);
    GslVector multVector(m_env,m_map);
    for (unsigned int j = 0; j < nCols; ++j) {
      if (j > 0) unitVector[j-1] = 0.;
      unitVector[j] = 1.;
      this->invertMultiply(unitVector, multVector);
      for (unsigned int i = 0; i < nRows; ++i) {
        (*m_inverse)(i,j) = multVector[i];
      }
    }
  }
  if (m_env.checkingLevel() >= 1) {
    *m_env.subDisplayFile() << "CHECKING In GslMatrix::inverse()"
                            << ": M.lnDet = "      << this->lnDeterminant()
                            << ", M^{-1}.lnDet = " << m_inverse->lnDeterminant()
                            << std::endl;
  }
  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 5)) {
    *m_env.subDisplayFile() << "In GslMatrix::inverse():"
                            << "\n M = "        << *this
                            << "\n M^{-1} = "   << *m_inverse
                            << "\n M*M^{-1} = " << (*this)*(*m_inverse)
                            << "\n M^{-1}*M = " << (*m_inverse)*(*this)
                            << std::endl;
  }

  return *m_inverse;
}

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 1187 of file GslMatrix.C.

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

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

{
  queso_require_equal_to_msg(this->numCols(), b.sizeLocal(), "matrix and rhs have incompatible sizes");

  GslVector x(m_env,m_map);
  this->invertMultiply(b,x);

  return x;
}

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 1199 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_require_equal_to_msg, queso_require_msg, QUESO::GslVector::sizeLocal(), and QUESO::BaseEnvironment::subDisplayFile().

{
  queso_require_equal_to_msg(this->numCols(), b.sizeLocal(), "matrix and rhs have incompatible sizes");

  queso_require_equal_to_msg(x.sizeLocal(), b.sizeLocal(), "solution and rhs have incompatible sizes");

  int iRC;
  if (m_LU == NULL) {
    queso_require_msg(!(m_permutation), "m_permutation should be NULL");

    m_LU = gsl_matrix_calloc(this->numRowsLocal(),this->numCols());
    queso_require_msg(m_LU, "gsl_matrix_calloc() failed");

    iRC = gsl_matrix_memcpy(m_LU, m_mat);
    queso_require_msg(!(iRC), "gsl_matrix_memcpy() failed");

    m_permutation = gsl_permutation_calloc(numCols());
    queso_require_msg(m_permutation, "gsl_permutation_calloc() failed");

    if (m_inDebugMode) {
      std::cout << "In GslMatrix::invertMultiply()"
                << ": before LU decomposition, m_LU = ";
      gsl_matrix_fprintf(stdout, m_LU, "%f");
      std::cout << std::endl;
    }

    gsl_error_handler_t* oldHandler;
    oldHandler = gsl_set_error_handler_off();
    if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 99)) {
      *m_env.subDisplayFile() << "In GslMatrix::invertMultiply()"
                              << ": before 'gsl_linalg_LU_decomp()'"
                              << std::endl;
    }
    iRC = gsl_linalg_LU_decomp(m_LU,m_permutation,&m_signum);
    if (iRC != 0) {
      std::cerr << "In GslMatrix::invertMultiply()"
                << ", after gsl_linalg_LU_decomp()"
                << ": iRC = " << iRC
                << ", gsl error message = " << gsl_strerror(iRC)
                << std::endl;
    }
    gsl_set_error_handler(oldHandler);
    if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 99)) {
      *m_env.subDisplayFile() << "In GslMatrix::invertMultiply()"
                              << ": after 'gsl_linalg_LU_decomp()'"
                              << ", IRC = " << iRC
                              << std::endl;
    }
    queso_require_msg(!(iRC), "gsl_linalg_LU_decomp() failed");

    if (m_inDebugMode) {
      std::cout << "In GslMatrix::invertMultiply()"
                << ": after LU decomposition, m_LU = ";
      gsl_matrix_fprintf(stdout, m_LU, "%f");
      std::cout << std::endl;
    }
  }

  gsl_error_handler_t* oldHandler;
  oldHandler = gsl_set_error_handler_off();
  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 99)) {
    *m_env.subDisplayFile() << "In GslMatrix::invertMultiply()"
                            << ": before 'gsl_linalg_LU_solve()'"
                            << std::endl;
  }
  iRC = gsl_linalg_LU_solve(m_LU,m_permutation,b.data(),x.data());
  if (iRC != 0) {
    m_isSingular = true;
    std::cerr << "In GslMatrix::invertMultiply()"
              << ", after gsl_linalg_LU_solve()"
              << ": iRC = " << iRC
              << ", gsl error message = " << gsl_strerror(iRC)
              << std::endl;
  }
  gsl_set_error_handler(oldHandler);
  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 99)) {
    *m_env.subDisplayFile() << "In GslMatrix::invertMultiply()"
                            << ": after 'gsl_linalg_LU_solve()'"
                            << ", IRC = " << iRC
                            << std::endl;
  }


  //                  m_env.worldRank(),
  //                  "GslMatrix::invertMultiply()",
  //                  "gsl_linalg_LU_solve() failed");

  if (m_inDebugMode) {
    GslVector tmpVec(b - (*this)*x);
    std::cout << "In GslMatrix::invertMultiply()"
              << ": ||b - Ax||_2 = "         << tmpVec.norm2()
              << ": ||b - Ax||_2/||b||_2 = " << tmpVec.norm2()/b.norm2()
              << std::endl;
  }

  return;
}

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 1300 of file GslMatrix.C.

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

{
  GslMatrix X(m_env,m_map,B.numCols());
  this->invertMultiply(B,X);

  return X;
}

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 1309 of file GslMatrix.C.

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

{
  // Sanity Checks
  queso_require_equal_to_msg(B.numRowsLocal(), X.numRowsLocal(),
                             "Matrices B and X are incompatible");
  queso_require_equal_to_msg(B.numCols(),      X.numCols(),
                             "Matrices B and X are incompatible");
  queso_require_equal_to_msg(this->numRowsLocal(), X.numRowsLocal(),
                             "This and X matrices are incompatible");

  // Some local variables used within the loop.
  GslVector b(m_env, m_map);
  GslVector x(m_env, m_map);

  for( unsigned int j = 0; j < B.numCols(); ++j )
    {
      b = B.getColumn( j );

      //invertMultiply will only do the LU once and store it. So we don't
      //need to worry about it doing LU multiple times.
      x = this->invertMultiply( b );

      X.setColumn( j, x );
    }

  return;
}

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 1338 of file GslMatrix.C.

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

{
  queso_require_equal_to_msg(this->numCols(), b.sizeLocal(), "matrix and rhs have incompatible sizes");

  GslVector x(m_env,m_map);
  this->invertMultiplyForceLU(b,x);

  return x;
}

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 1350 of file GslMatrix.C.

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

{
  queso_require_equal_to_msg(this->numCols(), b.sizeLocal(), "matrix and rhs have incompatible sizes");

  queso_require_equal_to_msg(x.sizeLocal(), b.sizeLocal(), "solution and rhs have incompatible sizes");

  int iRC;

  if ( m_LU == NULL ) {
    queso_require_msg(!(m_permutation), "m_permutation should be NULL");
    m_LU = gsl_matrix_calloc(this->numRowsLocal(),this->numCols());
  }
  queso_require_msg(m_LU, "gsl_matrix_calloc() failed");

  iRC = gsl_matrix_memcpy(m_LU, m_mat);
  queso_require_msg(!(iRC), "gsl_matrix_memcpy() failed");

  if( m_permutation == NULL ) m_permutation = gsl_permutation_calloc(numCols());
  queso_require_msg(m_permutation, "gsl_permutation_calloc() failed");

  iRC = gsl_linalg_LU_decomp(m_LU,m_permutation,&m_signum);
  queso_require_msg(!(iRC), "gsl_linalg_LU_decomp() failed");

  iRC = gsl_linalg_LU_solve(m_LU,m_permutation,b.data(),x.data());
  if (iRC != 0) {
    m_isSingular = true;
  }
  queso_require_msg(!(iRC), "gsl_linalg_LU_solve() failed");

  return;
}

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 1411 of file GslMatrix.C.

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

{

  // Sanity Checks
  unsigned int n = eigenVector.sizeLocal();

  queso_require_not_equal_to_msg(n, 0, "invalid input vector size");

  /* The following notation is used:
     z = vector used in iteration that ends up being the eigenvector corresponding to the
         largest eigenvalue
     w = vector used in iteration that we extract the largest eigenvalue from.  */

  // Some parameters associated with the algorithm
  // TODO: Do we want to add the ability to have these set by the user?
  const unsigned int max_num_iterations = 10000;
  const double tolerance = 1.0e-13;

  // Create temporary working vectors.
  // TODO: We shouldn't have to use these - we ought to be able to work directly
  // TODO: with eigenValue and eigenVector. Optimize this once we have regression
  // TODO: tests.
  GslVector z(m_env, m_map, 1.0 ); // Needs to be initialized to 1.0
  GslVector w(m_env, m_map);

  // Some variables we use inside the loop.
  int index;
  double residual;
  double lambda;

  for( unsigned int k = 0; k < max_num_iterations; ++k )
    {
      w = (*this) * z;

      // For this algorithm, it's crucial to get the maximum in
      // absolute value, but then to normalize by the actual value
      // and *not* the absolute value.
      index = (w.abs()).getMaxValueIndex();

      lambda = w[index];

      z = (1.0/lambda) * w;

      // Here we use the norm of the residual as our convergence check:
      // norm( A*x - \lambda*x )
      residual = ( (*this)*z - lambda*z ).norm2();

      if( residual < tolerance )
        {
          eigenValue = lambda;

          // TODO: Do we want to normalize this so eigenVector.norm2() = 1?
          eigenVector = z;
          return;
        }

    }

  // If we reach this point, then we didn't converge. Print error message
  // to this effect.
  // TODO: We know we failed at this point - need a way to not need a test
  // TODO: Just specify the error.
  queso_require_less_msg(residual, tolerance, "Maximum num iterations exceeded");


  return;
}

double QUESO::GslMatrix::lnDeterminant

(

)

const

Calculates the ln(determinant) of this matrix.

Definition at line 1077 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().

{
  if (m_lnDeterminant == -INFINITY) {
    if (m_LU == NULL) {
      GslVector tmpB(m_env,m_map);
      GslVector tmpX(m_env,m_map);
      if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 99)) {
        *m_env.subDisplayFile() << "In GslMatrix::lnDeterminant()"
                                << ": before 'this->invertMultiply()'"
                                << std::endl;
      }
      this->invertMultiply(tmpB,tmpX);
      if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 99)) {
        *m_env.subDisplayFile() << "In GslMatrix::lnDeterminant()"
                                << ": after 'this->invertMultiply()'"
                                << std::endl;
      }
    }
    if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 99)) {
      *m_env.subDisplayFile() << "In GslMatrix::lnDeterminant()"
                              << ": before 'gsl_linalg_LU_det()'"
                              << std::endl;
    }
    m_determinant   = gsl_linalg_LU_det(m_LU,m_signum);
    if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 99)) {
      *m_env.subDisplayFile() << "In GslMatrix::lnDeterminant()"
                              << ": after 'gsl_linalg_LU_det()'"
                              << std::endl;
    }
    m_lnDeterminant = gsl_linalg_LU_lndet(m_LU);
    if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 99)) {
      *m_env.subDisplayFile() << "In GslMatrix::lnDeterminant()"
                              << ": before 'gsl_linalg_LU_lndet()'"
                              << std::endl;
    }
  }

  return m_lnDeterminant;
}

void QUESO::GslMatrix::matlabLinearInterpExtrap	(	const GslVector &	x1Vec,
		const GslMatrix &	y1Mat,
		const GslVector &	x2Vec
	)

Definition at line 1706 of file GslMatrix.C.

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

{
  queso_require_greater_msg(x1Vec.sizeLocal(), 1, "invalid 'x1' size");

  queso_require_equal_to_msg(x1Vec.sizeLocal(), y1Mat.numRowsLocal(), "invalid 'x1' and 'y1' sizes");

  queso_require_equal_to_msg(x2Vec.sizeLocal(), this->numRowsLocal(), "invalid 'x2' and 'this' sizes");

  queso_require_equal_to_msg(y1Mat.numCols(), this->numCols(), "invalid 'y1' and 'this' sizes");

  GslVector y1Vec(x1Vec);
  GslVector y2Vec(x2Vec);
  for (unsigned int colId = 0; colId < y1Mat.numCols(); ++colId) {
    y1Mat.getColumn(colId,y1Vec);
    y2Vec.matlabLinearInterpExtrap(x1Vec,y1Vec,x2Vec);
    this->setColumn(colId,y2Vec);
  }

  return;
}

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().

{
  double value = -INFINITY;

  unsigned int nRows = this->numRowsLocal();
  unsigned int nCols = this->numCols();
  double aux = 0.;
  for (unsigned int i = 0; i < nRows; i++) {
    for (unsigned int j = 0; j < nCols; j++) {
      aux = (*this)(i,j);
      if (aux > value) value = aux;
    }
  }

  return value;
}

void QUESO::GslMatrix::mpiSum	(	const MpiComm &	comm,
		GslMatrix &	M_global
	)		const

Definition at line 1663 of file GslMatrix.C.

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

{
  // Sanity Checks
  UQ_FATAL_RC_MACRO(((this->numRowsLocal() != M_global.numRowsLocal()) ||
                     (this->numCols()      != M_global.numCols()     )),
                    env().fullRank(),
                    "GslMatrix::mpiSum()",
                    "local and global matrices incompatible");

  /* TODO: Probably a better way to handle this unpacking/packing of data */
  int size = M_global.numRowsLocal()*M_global.numCols();
  std::vector<double> local( size, 0.0 );
  std::vector<double> global( size, 0.0 );

  int k;
  for( unsigned int i = 0; i < this->numRowsLocal(); i++ )
    {
      for( unsigned int j = 0; j < this->numCols(); j++ )
        {
          k = i + j*M_global.numCols();

          local[k] = (*this)(i,j);
        }
    }

  comm.Allreduce<double>(&local[0], &global[0], size, RawValue_MPI_SUM,
                 "GslMatrix::mpiSum()",
                 "failed MPI.Allreduce()");

  for( unsigned int i = 0; i < this->numRowsLocal(); i++ )
    {
      for( unsigned int j = 0; j < this->numCols(); j++ )
        {
          k = i + j*M_global.numCols();

          M_global(i,j) = global[k];
        }
    }

  return;
}

GslVector QUESO::GslMatrix::multiply

(

const GslVector &

x

)

const

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

Definition at line 1153 of file GslMatrix.C.

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

Referenced by QUESO::operator*().

{
  queso_require_equal_to_msg(this->numCols(), x.sizeLocal(), "matrix and vector have incompatible sizes");

  GslVector y(m_env,m_map);
  this->multiply(x,y);

  return y;
}

void QUESO::GslMatrix::multiply ( const GslVector &  x, GslVector &  y  ) const

private

This function multiplies this matrix by vector x and stores the resulting vector in y.

Definition at line 1165 of file GslMatrix.C.

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

{
  queso_require_equal_to_msg(this->numCols(), x.sizeLocal(), "matrix and x have incompatible sizes");

  queso_require_equal_to_msg(this->numRowsLocal(), y.sizeLocal(), "matrix and y have incompatible sizes");

  unsigned int sizeX = this->numCols();
  unsigned int sizeY = this->numRowsLocal();
  for (unsigned int i = 0; i < sizeY; ++i) {
    double value = 0.;
    for (unsigned int j = 0; j < sizeX; ++j) {
      value += (*this)(i,j)*x[j];
    }
    y[i] = value;
  }

  return;
}

double QUESO::GslMatrix::normFrob

(

)

const

Returns the Frobenius norm of this matrix.

Definition at line 293 of file GslMatrix.C.

References numCols(), and numRowsLocal().

{
  double value = 0.;

  unsigned int nRows = this->numRowsLocal();
  unsigned int nCols = this->numCols();
  double aux = 0.;
  for (unsigned int i = 0; i < nRows; i++) {
    for (unsigned int j = 0; j < nCols; j++) {
      aux = (*this)(i,j);
      value += aux*aux;
    }
  }

  return sqrt(value);
}

double QUESO::GslMatrix::normMax

(

)

const

Returns the Frobenius norm of this matrix.

Definition at line 311 of file GslMatrix.C.

References numCols(), and numRowsLocal().

{
  double value = 0.;

  unsigned int nRows = this->numRowsLocal();
  unsigned int nCols = this->numCols();
  double aux = 0.;
  for (unsigned int i = 0; i < nRows; i++) {
    for (unsigned int j = 0; j < nCols; j++) {
      aux = fabs((*this)(i,j));
      if (aux > value) value = aux;
    }
  }

  return value;
}

unsigned int QUESO::GslMatrix::numCols ( ) const

virtual

Returns the column dimension of this matrix.

Implements QUESO::Matrix.

Definition at line 287 of file GslMatrix.C.

References m_mat.

Referenced by cwExtract(), cwSet(), fillWithBlocksDiagonally(), fillWithBlocksHorizontally(), fillWithBlocksVertically(), fillWithTensorProduct(), fillWithTranspose(), filterLargeValues(), filterSmallValues(), getColumn(), getRow(), internalSvd(), inverse(), invertMultiply(), invertMultiplyForceLU(), QUESO::leftDiagScaling(), matlabLinearInterpExtrap(), max(), mpiSum(), multiply(), normFrob(), normMax(), QUESO::operator*(), print(), rank(), QUESO::rightDiagScaling(), setColumn(), setRow(), subReadContents(), subWriteContents(), svd(), svdSolve(), transpose(), zeroLower(), and zeroUpper().

{
  return m_mat->size2;
}

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.

{
  return m_mat->size1;
}

unsigned int QUESO::GslMatrix::numRowsLocal ( ) const

virtual

Returns the local row dimension of this matrix.

Implements QUESO::Matrix.

Definition at line 275 of file GslMatrix.C.

References m_mat.

Referenced by cwExtract(), cwSet(), eigen(), fillWithBlocksDiagonally(), fillWithBlocksHorizontally(), fillWithBlocksVertically(), fillWithTensorProduct(), fillWithTranspose(), filterLargeValues(), filterSmallValues(), QUESO::GaussianLikelihoodBlockDiagonalCovariance< V, M >::GaussianLikelihoodBlockDiagonalCovariance(), QUESO::GaussianLikelihoodBlockDiagonalCovarianceRandomCoefficients< V, M >::GaussianLikelihoodBlockDiagonalCovarianceRandomCoefficients(), getColumn(), getRow(), internalSvd(), inverse(), invertMultiply(), invertMultiplyForceLU(), QUESO::leftDiagScaling(), matlabLinearInterpExtrap(), max(), mpiSum(), multiply(), normFrob(), normMax(), QUESO::operator*(), print(), rank(), QUESO::rightDiagScaling(), setColumn(), setRow(), subReadContents(), subWriteContents(), svd(), svdSolve(), transpose(), zeroLower(), and zeroUpper().

{
  return m_mat->size1;
}

double& QUESO::GslMatrix::operator() ( unsigned int  i, unsigned int  j  )

inline

Element access method (non-const).

Definition at line 104 of file GslMatrix.h.

  {
    // Changing the matrix changes its decomposition(s).  We'll
    // invalidate for now; recalculate later if needed.
    this->resetLU();
    queso_require_less_msg(i, m_mat->size1, "i is too large");
    queso_require_less_msg(j, m_mat->size2, "j is too large");
    return *gsl_matrix_ptr(m_mat,i,j);
  }

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, and queso_require_less_msg.

  {
    queso_require_less_msg(i, m_mat->size1, "i is too large");
    queso_require_less_msg(j, m_mat->size2, "j is too large");
    return *gsl_matrix_ptr(m_mat,i,j);
  }

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, queso_require_msg, and resetLU().

{
  this->resetLU();
  int iRC;
  iRC = gsl_matrix_scale(m_mat,a);
  queso_require_msg(!(iRC), "scaling failed");
  return *this;
}

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, queso_require_msg, and resetLU().

{
  this->resetLU();
  int iRC;
  iRC = gsl_matrix_add(m_mat,rhs.m_mat);
  queso_require_msg(!(iRC), "failed");

  return *this;
}

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, queso_require_msg, and resetLU().

{
  this->resetLU();
  int iRC;
  iRC = gsl_matrix_sub(m_mat,rhs.m_mat);
  queso_require_msg(!(iRC), "failed");

  return *this;
}

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().

{
  this->resetLU();
  *this *= (1./a);

  return *this;
}

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().

{
  this->resetLU();
  this->copy(obj);
  return *this;
}

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 1737 of file GslMatrix.C.

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

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

{
  unsigned int nRows = this->numRowsLocal();
  unsigned int nCols = this->numCols();

  if (m_printHorizontally) {
    for (unsigned int i = 0; i < nRows; ++i) {
      for (unsigned int j = 0; j < nCols; ++j) {
        os << (*this)(i,j)
           << " ";
      }
      if (i != (nRows-1)) os << "; ";
    }
    //os << std::endl;
  }
  else {
    for (unsigned int i = 0; i < nRows; ++i) {
      for (unsigned int j = 0; j < nCols; ++j) {
        os << (*this)(i,j)
           << " ";
      }
      os << std::endl;
    }
  }

  return;
}

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 1118 of file GslMatrix.C.

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

{
  int iRC = 0;
  iRC = internalSvd();
  if (iRC) {}; // just to remove compiler warning

  GslVector relativeVec(*m_svdSvec);
  if (relativeVec[0] > 0.) {
    relativeVec = (1./relativeVec[0])*relativeVec;
  }

  unsigned int rankValue = 0;
  for (unsigned int i = 0; i < relativeVec.sizeLocal(); ++i) {
    if (( (*m_svdSvec)[i] >= absoluteZeroThreshold ) &&
        ( relativeVec [i] >= relativeZeroThreshold )) {
       rankValue += 1;
    }
  }

  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 3)) {
    *m_env.subDisplayFile() << "In GslMatrix::rank()"
                            << ": this->numRowsLocal() = "  << this->numRowsLocal()
                            << ", this->numCols() = "       << this->numCols()
                            << ", absoluteZeroThreshold = " << absoluteZeroThreshold
                            << ", relativeZeroThreshold = " << relativeZeroThreshold
                            << ", rankValue = "             << rankValue
                            << ", m_svdSvec = "             << *m_svdSvec
                            << ", relativeVec = "           << relativeVec
                            << std::endl;
  }

  return rankValue;
}

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().

{
  if (m_LU) {
    gsl_matrix_free(m_LU);
    m_LU = NULL;
  }
  if (m_inverse) {
    delete m_inverse;
    m_inverse = NULL;
  }
  if (m_svdColMap) {
    delete m_svdColMap;
    m_svdColMap = NULL;
  }
  if (m_svdUmat) {
    delete m_svdUmat;
    m_svdUmat = NULL;
  }
  if (m_svdSvec) {
    delete m_svdSvec;
    m_svdSvec = NULL;
  }
  if (m_svdVmat) {
    delete m_svdVmat;
    m_svdVmat = NULL;
  }
  if (m_svdVTmat) {
    delete m_svdVTmat;
    m_svdVTmat = NULL;
  }
  m_determinant   = -INFINITY;
  m_lnDeterminant = -INFINITY;
  if (m_permutation) {
    gsl_permutation_free(m_permutation);
    m_permutation = NULL;
  }
  m_signum = 0;
  m_isSingular = false;

  return;
}

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 1646 of file GslMatrix.C.

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

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

{
  this->resetLU();
  // Sanity checks
  queso_require_less_msg(column_num, this->numCols(), "Specified column number not within range");

  queso_require_equal_to_msg(column.sizeLocal(), this->numRowsLocal(), "column vector not same size as this matrix");

  gsl_vector* gsl_column = column.data();

  int error_code = gsl_matrix_set_col( m_mat, column_num, gsl_column );
  queso_require_equal_to_msg(error_code, 0, "gsl_matrix_set_col failed");

  return;
}

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 1629 of file GslMatrix.C.

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

{
  this->resetLU();
  // Sanity checks
  queso_require_less_msg(row_num, this->numRowsLocal(), "Specified row number not within range");

  queso_require_equal_to_msg(row.sizeLocal(), this->numCols(), "row vector not same size as this matrix");

  gsl_vector* gsl_row = row.data();

  int error_code = gsl_matrix_set_row( m_mat, row_num, gsl_row );
  queso_require_equal_to_msg(error_code, 0, "gsl_matrix_set_row failed");

  return;
}

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 1480 of file GslMatrix.C.

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

{
  // Sanity Checks
  unsigned int n = eigenVector.sizeLocal();

  queso_require_not_equal_to_msg(n, 0, "invalid input vector size");

  /* The following notation is used:
     z = vector used in iteration that ends up being the eigenvector corresponding to the
         largest eigenvalue
     w = vector used in iteration that we extract the largest eigenvalue from.  */

  // Some parameters associated with the algorithm
  // TODO: Do we want to add the ability to have these set by the user?
  const unsigned int max_num_iterations = 1000;
  const double tolerance = 1.0e-13;

  // Create temporary working vectors.
  // TODO: We shouldn't have to use these - we ought to be able to work directly
  // TODO: with eigenValue and eigenVector. Optimize this once we have regression
  // TODO: tests.
  GslVector z(m_env, m_map, 1.0 ); // Needs to be initialized to 1.0
  GslVector w(m_env, m_map);

  // Some variables we use inside the loop.
  int index;
  double residual;
  double one_over_lambda;
  double lambda;

  for( unsigned int k = 0; k < max_num_iterations; ++k )
    {
      w = (*this).invertMultiplyForceLU( z );

      // For this algorithm, it's crucial to get the maximum in
      // absolute value, but then to normalize by the actual value
      // and *not* the absolute value.
      index = (w.abs()).getMaxValueIndex();

      // Remember: Inverse power method solves for max 1/lambda ==> lambda smallest
      one_over_lambda = w[index];

      lambda = 1.0/one_over_lambda;

      z = lambda * w;

      // Here we use the norm of the residual as our convergence check:
      // norm( A*x - \lambda*x )
      residual = ( (*this)*z - lambda*z ).norm2();

      if( residual < tolerance )
        {
          eigenValue = lambda;

          // TODO: Do we want to normalize this so eigenVector.norm2() = 1?
          eigenVector = z;
          return;
        }

    }

  // If we reach this point, then we didn't converge. Print error message
  // to this effect.
  // TODO: We know we failed at this point - need a way to not need a test
  // TODO: Just specify the error.
  queso_require_less_msg(residual, tolerance, "Maximum num iterations exceeded");

  return;
}

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 1806 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_require_equal_to_msg, queso_require_greater_equal_msg, queso_require_less_equal_msg, queso_require_less_msg, QUESO::BaseEnvironment::subDisplayFile(), and QUESO::BaseEnvironment::subRank().

{
  queso_require_greater_equal_msg(m_env.subRank(), 0, "unexpected subRank");

  queso_require_less_equal_msg(this->numOfProcsForStorage(), 1, "implemented just for sequential vectors for now");

  FilePtrSetStruct filePtrSet;
  if (m_env.openInputFile(fileName,
                          fileType, // "m or hdf"
                          allowedSubEnvIds,
                          filePtrSet)) {

    // palms
    unsigned int nRowsLocal = this->numRowsLocal();

    // In the logic below, the id of a line' begins with value 0 (zero)
    unsigned int idOfMyFirstLine = 1;
    unsigned int idOfMyLastLine = nRowsLocal;
    unsigned int nCols = this->numCols();

    // Read number of matrix rows in the file by taking care of the first line,
    // which resembles something like 'variable_name = zeros(n_rows,n_cols);'
    std::string tmpString;

    // Read 'variable name' string
    *filePtrSet.ifsVar >> tmpString;
    //std::cout << "Just read '" << tmpString << "'" << std::endl;

    // Read '=' sign
    *filePtrSet.ifsVar >> tmpString;
    //std::cout << "Just read '" << tmpString << "'" << std::endl;
    queso_require_equal_to_msg(tmpString, "=", "string should be the '=' sign");

    // Read 'zeros(n_rows,n_cols)' string
    *filePtrSet.ifsVar >> tmpString;
    //std::cout << "Just read '" << tmpString << "'" << std::endl;
    unsigned int posInTmpString = 6;

    // Isolate 'n_rows' in a string
    char nRowsString[tmpString.size()-posInTmpString+1];
    unsigned int posInRowsString = 0;
    do {
      queso_require_less_msg(posInTmpString, tmpString.size(), "symbol ',' not found in first line of file");
      nRowsString[posInRowsString++] = tmpString[posInTmpString++];
    } while (tmpString[posInTmpString] != ',');
    nRowsString[posInRowsString] = '\0';

    // Isolate 'n_cols' in a string
    posInTmpString++; // Avoid reading ',' char
    char nColsString[tmpString.size()-posInTmpString+1];
    unsigned int posInColsString = 0;
    do {
      queso_require_less_msg(posInTmpString, tmpString.size(), "symbol ')' not found in first line of file");
      nColsString[posInColsString++] = tmpString[posInTmpString++];
    } while (tmpString[posInTmpString] != ')');
    nColsString[posInColsString] = '\0';

    // Convert 'n_rows' and 'n_cols' strings to numbers
    unsigned int numRowsInFile = (unsigned int) strtod(nRowsString,NULL);
    unsigned int numColsInFile = (unsigned int) strtod(nColsString,NULL);
    if (m_env.subDisplayFile()) {
      *m_env.subDisplayFile() << "In GslMatrix::subReadContents()"
                              << ": fullRank "        << m_env.fullRank()
                              << ", numRowsInFile = " << numRowsInFile
                              << ", numColsInFile = " << numColsInFile
                              << ", nRowsLocal = "    << nRowsLocal
                              << ", nCols = "         << nCols
                              << std::endl;
    }

    // Check if [num of rows in file] == [requested matrix row size]
    queso_require_equal_to_msg(numRowsInFile, nRowsLocal, "size of vec in file is not big enough");

    // Check if [num of cols in file] == [num cols in current matrix]
    queso_require_equal_to_msg(numColsInFile, nCols, "number of parameters of vec in file is different than number of parameters in this vec object");

    // Code common to any core in a communicator
    unsigned int maxCharsPerLine = 64*nCols; // Up to about 60 characters to represent each parameter value

    unsigned int lineId = 0;
    while (lineId < idOfMyFirstLine) {
      filePtrSet.ifsVar->ignore(maxCharsPerLine,'\n');
      lineId++;
    };

    if (m_env.subDisplayFile()) {
      *m_env.subDisplayFile() << "In GslMatrix::subReadContents()"
                              << ": beginning to read input actual data"
                              << std::endl;
    }

    // Take care of initial part of the first data line,
    // which resembles something like 'variable_name = [value1 value2 ...'
    // Read 'variable name' string
    *filePtrSet.ifsVar >> tmpString;
    //std::cout << "Core 0 just read '" << tmpString << "'" << std::endl;

    // Read '=' sign
    *filePtrSet.ifsVar >> tmpString;
    //std::cout << "Core 0 just read '" << tmpString << "'" << std::endl;
    queso_require_equal_to_msg(tmpString, "=", "in core 0, string should be the '=' sign");

    // Take into account the ' [' portion
    std::streampos tmpPos = filePtrSet.ifsVar->tellg();
    filePtrSet.ifsVar->seekg(tmpPos+(std::streampos)2);

    if (m_env.subDisplayFile()) {
      *m_env.subDisplayFile() << "In GslMatrix::subReadContents()"
                              << ": beginning to read lines with numbers only"
                              << ", lineId = "          << lineId
                              << ", idOfMyFirstLine = " << idOfMyFirstLine
                              << ", idOfMyLastLine = "  << idOfMyLastLine
                              << std::endl;
    }

    double tmpRead;
    while (lineId <= idOfMyLastLine) {
      for (unsigned int j = 0; j < nCols; ++j) {
        *filePtrSet.ifsVar >> tmpRead;
        (*this)(lineId-idOfMyFirstLine,j) = tmpRead;
      }
      lineId++;
    };

    m_env.closeFile(filePtrSet,fileType);
  }

  return;
}

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 1766 of file GslMatrix.C.

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

{
  queso_require_greater_equal_msg(m_env.subRank(), 0, "unexpected subRank");

  queso_require_less_equal_msg(this->numOfProcsForStorage(), 1, "implemented just for sequential vectors for now");

  FilePtrSetStruct filePtrSet;
  if (m_env.openOutputFile(fileName,
                           fileType, // "m or hdf"
                           allowedSubEnvIds,
                           false,
                           filePtrSet)) {
    unsigned int nRows = this->numRowsLocal();
    unsigned int nCols = this->numCols();
    *filePtrSet.ofsVar << varNamePrefix << "_sub" << m_env.subIdString() << " = zeros(" << nRows
                       << ","                                                           << nCols
                       << ");"
                       << std::endl;
    *filePtrSet.ofsVar << varNamePrefix << "_sub" << m_env.subIdString() << " = [";

    for (unsigned int i = 0; i < nRows; ++i) {
      for (unsigned int j = 0; j < nCols; ++j) {
        *filePtrSet.ofsVar << (*this)(i,j)
                << " ";
      }
      *filePtrSet.ofsVar << "\n";
    }
    *filePtrSet.ofsVar << "];\n";

    m_env.closeFile(filePtrSet,fileType);
  }

  return;
}

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_require_equal_to_msg, queso_require_msg, and QUESO::GslVector::sizeLocal().

{
  unsigned int nRows = this->numRowsLocal();
  unsigned int nCols = this->numCols();

  queso_require_msg(!((matU.numRowsLocal() != nRows) || (matU.numCols() != nCols)), "invalid matU");

  queso_require_equal_to_msg(vecS.sizeLocal(), nCols, "invalid vecS");

  queso_require_msg(!((matVt.numRowsLocal() != nCols) || (matVt.numCols() != nCols)), "invalid matVt");

  int iRC = internalSvd();

  matU  = *m_svdUmat;
  vecS  = *m_svdSvec;
  matVt = *m_svdVTmat;

  return iRC;
}

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.

{
  int iRC = 0;
  iRC = internalSvd();
  if (iRC) {}; // just to remove compiler warning

  return *m_svdUmat;
}

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.

{
  int iRC = 0;
  iRC = internalSvd();
  if (iRC) {}; // just to remove compiler warning

  return *m_svdVmat;
}

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_require_equal_to_msg, QUESO::GslVector::sizeLocal(), and QUESO::BaseEnvironment::subDisplayFile().

Referenced by svdSolve().

{
  unsigned int nRows = this->numRowsLocal();
  unsigned int nCols = this->numCols();

  queso_require_equal_to_msg(rhsVec.sizeLocal(), nRows, "invalid rhsVec");

  queso_require_equal_to_msg(solVec.sizeLocal(), nCols, "invalid solVec");

  int iRC = internalSvd();

  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 5)) {
    *m_env.subDisplayFile() << "In GslMatrix::svdSolve():"
                            << "\n this->numRowsLocal()      = " << this->numRowsLocal()
                            << ", this->numCols()      = "       << this->numCols()
                            << "\n m_svdUmat->numRowsLocal() = " << m_svdUmat->numRowsLocal()
                            << ", m_svdUmat->numCols() = "       << m_svdUmat->numCols()
                            << "\n m_svdVmat->numRowsLocal() = " << m_svdVmat->numRowsLocal()
                            << ", m_svdVmat->numCols() = "       << m_svdVmat->numCols()
                            << "\n m_svdSvec->sizeLocal()    = " << m_svdSvec->sizeLocal()
                            << "\n rhsVec.sizeLocal()        = " << rhsVec.sizeLocal()
                            << "\n solVec.sizeLocal()        = " << solVec.sizeLocal()
                            << std::endl;
  }

  if (iRC == 0) iRC = gsl_linalg_SV_solve(m_svdUmat->data(), m_svdVmat->data(), m_svdSvec->data(), rhsVec.data(), solVec.data());

  return iRC;
}

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_require_equal_to_msg, setColumn(), and svdSolve().

{
  unsigned int nRows = this->numRowsLocal();
  unsigned int nCols = this->numCols();

  queso_require_equal_to_msg(rhsMat.numRowsLocal(), nRows, "invalid rhsMat");

  queso_require_equal_to_msg(solMat.numRowsLocal(), nCols, "invalid solMat");

  queso_require_equal_to_msg(rhsMat.numCols(), solMat.numCols(), "rhsMat and solMat are not compatible");

  GslVector rhsVec(m_env,rhsMat.map());
  GslVector solVec(m_env,solMat.map());
  int iRC = 0;
  for (unsigned int j = 0; j < rhsMat.numCols(); ++j) {
    rhsVec = rhsMat.getColumn(j);
    iRC = this->svdSolve(rhsVec, solVec);
    if (iRC) break;
    solMat.setColumn(j,solVec);
  }

  return iRC;
}

GslMatrix QUESO::GslMatrix::transpose

(

)

const

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

Definition at line 693 of file GslMatrix.C.

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

Referenced by internalSvd().

{
  unsigned int nRows = this->numRowsLocal();
  unsigned int nCols = this->numCols();

  queso_require_equal_to_msg(nRows, nCols, "routine works only for square matrices");

  GslMatrix mat(m_env,m_map,nCols);
  for (unsigned int row = 0; row < nRows; ++row) {
    for (unsigned int col = 0; col < nCols; ++col) {
      mat(row,col) = (*this)(col,row);
    }
  }

  return mat;
}

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_require_equal_to_msg, and resetLU().

{
  unsigned int nRows = this->numRowsLocal();
  unsigned int nCols = this->numCols();

  queso_require_equal_to_msg(nRows, nCols, "routine works only for square matrices");

  this->resetLU();

  if (includeDiagonal) {
    for (unsigned int i = 0; i < nRows; i++) {
      for (unsigned int j = 0; j <= i; j++) {
        (*this)(i,j) = 0.;
      }
    }
  }
  else {
    for (unsigned int i = 0; i < nRows; i++) {
      for (unsigned int j = 0; j < i; j++) {
        (*this)(i,j) = 0.;
      }
    }
  }

  return;
}

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_require_equal_to_msg, and resetLU().

{
  unsigned int nRows = this->numRowsLocal();
  unsigned int nCols = this->numCols();

  queso_require_equal_to_msg(nRows, nCols, "routine works only for square matrices");

  this->resetLU();

  if (includeDiagonal) {
    for (unsigned int i = 0; i < nRows; i++) {
      for (unsigned int j = i; j < nCols; j++) {
        (*this)(i,j) = 0.;
      }
    }
  }
  else {
    for (unsigned int i = 0; i < nRows; i++) {
      for (unsigned int j = (i+1); j < nCols; j++) {
        (*this)(i,j) = 0.;
      }
    }
  }

  return;
}

Member Data Documentation

double QUESO::GslMatrix::m_determinant

mutableprivate

The determinant of this matrix.

Definition at line 438 of file GslMatrix.h.

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

GslMatrix* QUESO::GslMatrix::m_inverse

mutableprivate

Inverse matrix of this.

Definition at line 407 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 456 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 441 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 404 of file GslMatrix.h.

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

gsl_matrix* QUESO::GslMatrix::m_mat

private

GSL matrix, also referred to as this matrix.

Definition at line 401 of file GslMatrix.h.

Referenced by chol(), copy(), cwSet(), data(), eigen(), getColumn(), getRow(), GslMatrix(), invertMultiply(), invertMultiplyForceLU(), numCols(), numRowsGlobal(), numRowsLocal(), operator()(), operator*=(), operator+=(), operator-=(), setColumn(), setRow(), and ~GslMatrix().

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 447 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 453 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 410 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 423 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 416 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 429 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 435 of file GslMatrix.h.

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

---

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

• src/core/inc/GslMatrix.h
• src/core/src/GslMatrix.C