QUESO::TeuchosMatrix Class Reference

Class for matrix operations using Teuchos (Trilinos). More...

#include <TeuchosMatrix.h>

Inheritance diagram for QUESO::TeuchosMatrix:

Public Member Functions
Constructor/Destructor methods
	TeuchosMatrix (const BaseEnvironment &env, const Map &map, unsigned int numCols)
	Shaped Constructor: creates a shaped matrix with numCols columns. More...
	TeuchosMatrix (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...
	TeuchosMatrix (const TeuchosVector &v, double diagValue)
	Shaped Constructor: creates a square matrix with size v.sizeLocal() and diagonal values all equal to diagValue. More...
	TeuchosMatrix (const TeuchosVector &v)
	Shaped Constructor: creates a square matrix with size v.sizeLocal(). More...
	TeuchosMatrix (const TeuchosMatrix &B)
	Shaped Constructor: this matrix is a copy of matrix B. More...
	~TeuchosMatrix ()
	Destructor. More...
Set methods
TeuchosMatrix &	operator= (const TeuchosMatrix &rhs)
	Copies values from matrix rhs to this. More...
TeuchosMatrix &	operator*= (double a)
	Stores in this the coordinate-wise multiplication of this and a. More...
TeuchosMatrix &	operator/= (double a)
	Stores in this the coordinate-wise division of this by a. More...
TeuchosMatrix &	operator+= (const TeuchosMatrix &rhs)
	Stores in this the coordinate-wise addition of this and rhs. More...
TeuchosMatrix &	operator-= (const TeuchosMatrix &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 *	values ()
	Returns a pointer to the first element of this matrix. More...
int	stride ()
	Returns the stride between the columns of this matrix in memory. 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...
TeuchosMatrix	transpose () const
	This function calculates the transpose of this matrix (square). More...
TeuchosMatrix	inverse () const
	This function calculates 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 this, a real symmetric positive definite matrix. More...
int	svd (TeuchosMatrix &matU, TeuchosVector &vecS, TeuchosMatrix &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 TeuchosMatrix &	svdMatU () const
	This function calls private member TeuchosMatrix::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 TeuchosMatrix &	svdMatV () const
	This function calls private member TeuchosMatrix::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 TeuchosVector &rhsVec, TeuchosVector &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 TeuchosMatrix::svd (x=solVec, b=rhsVec). More...
int	svdSolve (const TeuchosMatrix &rhsMat, TeuchosMatrix &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 TeuchosMatrix::svd (x=solMat, b=rhsMat). More...
TeuchosVector	multiply (const TeuchosVector &x) const
	This function multiplies this matrix by vector x and returns a vector. More...
void	invertMultiply (const TeuchosVector &b, TeuchosVector &x) const
	This function calculates the inverse of this matrix, multiplies it with vector b and stores the result in vector x. More...
TeuchosVector	invertMultiply (const TeuchosVector &b) const
	This function calculates the inverse of this matrix and multiplies it with vector b. More...
void	invertMultiply (const TeuchosMatrix &B, TeuchosMatrix &X) const
	This function calculates the inverse of this matrix, multiplies it with matrix B and stores the result in matrix X. More...
TeuchosMatrix	invertMultiply (const TeuchosMatrix &B) const
	This function calculates the inverse of this matrix and multiplies it with matrix B. More...
void	invertMultiplyForceLU (const TeuchosVector &b, TeuchosVector &x) const
	This function calculates the inverse of this matrix, multiplies it with vector b and stores the result in vector x. More...
TeuchosVector	invertMultiplyForceLU (const TeuchosVector &b) const
	This function calculates the inverse of this matrix and multiplies it with vector b. More...
void	eigen (TeuchosVector &eigenValues, TeuchosMatrix *eigenVectors) const
	This function computes the eigenvalues of a real symmetric matrix. More...
void	largestEigen (double &eigenValue, TeuchosVector &eigenVector) const
	This function finds largest eigenvalue, namely eigenValue, of this matrix and its corresponding eigenvector, namely eigenVector. More...
void	smallestEigen (double &eigenValue, TeuchosVector &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 TeuchosMatrix &mat)
	Set the components of which positions are greater than (rowId,colId) with the value of mat(rowId,colId). More...
void	getColumn (const unsigned int column_num, TeuchosVector &column) const
	This function gets the column_num-th column of this matrix and stores it into vector column. More...
TeuchosVector	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 TeuchosVector &column)
	This function copies vector column into the column_num-th column of this matrix. More...
void	getRow (const unsigned int row_num, TeuchosVector &row) const
	This function gets the row_num-th column of this matrix and stores it into vector row. More...
TeuchosVector	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 TeuchosVector &row)
	This function copies vector row into the row_num-th column of this matrix. More...
void	zeroLower (bool includeDiagonal=false)
	This function sets all the entries above the main diagonal (inclusive or not) of this matrix to zero. More...
void	zeroUpper (bool includeDiagonal=false)
	This function sets all the entries bellow the main diagonal (inclusive or not) 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 (const TeuchosMatrix &mat)
	This function stores the transpose of this matrix into this matrix. More...
void	fillWithBlocksDiagonally (const std::vector< const TeuchosMatrix * > &matrices)
	This function fills this matrix diagonally with const block matrices. More...
void	fillWithBlocksDiagonally (const std::vector< TeuchosMatrix * > &matrices)
	This function fills this matrix diagonally with block matrices. More...
void	fillWithBlocksHorizontally (const std::vector< const TeuchosMatrix * > &matrices)
	This function fills this matrix horizontally with const block matrices. More...
void	fillWithBlocksHorizontally (const std::vector< TeuchosMatrix * > &matrices)
	This function fills this matrix horizontally with block matrices. More...
void	fillWithBlocksVertically (const std::vector< const TeuchosMatrix * > &matrices)
	This function fills this matrix vertically with const block matrices. More...
void	fillWithBlocksVertically (const std::vector< TeuchosMatrix * > &matrices)
	This function fills this matrix vertically with block matrices. More...
void	fillWithTensorProduct (const TeuchosMatrix &mat1, const TeuchosMatrix &mat2)
	This function calculates the tensor product of matrices mat1 and mat2 and stores it in this matrix. More...
void	fillWithTensorProduct (const TeuchosMatrix &mat1, const TeuchosVector &vec2)
	This function calculates the tensor product of matrix mat1 and vector vec2 and stores it in this matrix. More...
Miscellaneous methods
void	mpiSum (const MpiComm &comm, TeuchosMatrix &M_global) const
void	matlabLinearInterpExtrap (const TeuchosVector &x1Vec, const TeuchosMatrix &y1Mat, const TeuchosVector &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	subReadContents (const std::string &fileName, const std::string &fileType, const std::set< unsigned int > &allowedSubEnvIds)
	Read contents of subenvironment from file fileName. 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...
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
	TeuchosMatrix ()
	Default Constructor. More...
void	copy (const TeuchosMatrix &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 TeuchosVector &x, TeuchosVector &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
Teuchos::SerialDenseMatrix < int, double >	m_mat
	Teuchos matrix, also referred to as this matrix. More...
Teuchos::SerialDenseMatrix < int, double >	m_LU
	Teuchos matrix for the LU decomposition of m_mat. More...
TeuchosMatrix *	m_inverse
	Stores the inverse of this matrix. More...
Map *	m_svdColMap
	Mapping for matrices involved in the singular value decomposition (svd) routine. More...
TeuchosMatrix *	m_svdUmat
	m_svdUmat stores the M-by-N orthogonal matrix U after the singular value decomposition of a matrix. More...
TeuchosVector *	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...
TeuchosMatrix *	m_svdVmat
	m_svdVmat stores the N-by-N orthogonal square matrix V after the singular value decomposition of a matrix. More...
TeuchosMatrix *	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...
TeuchosMatrix *	m_permutation
int *	v_pivoting
	The pivoting vector of a LU decomposition. More...
int	m_signum
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...
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 Teuchos (Trilinos).

This class creates and provides basic support for matrices of templated type as a specialization of Matrix using Teuchos matrices (from Trilinos), which are defined by an encapsulated Teuchos::SerialDenseMatrix structure.

Definition at line 53 of file TeuchosMatrix.h.

Constructor & Destructor Documentation

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

Shaped Constructor: creates a shaped matrix with numCols columns.

Definition at line 43 of file TeuchosMatrix.C.

References m_LU, m_mat, and QUESO::Map::NumGlobalElements().

  :
  Matrix  (env,map),
  m_inverse      (NULL),
  m_svdColMap    (NULL),
  m_svdUmat      (NULL),
  m_svdSvec      (NULL),
  m_svdVmat      (NULL),
  m_svdVTmat     (NULL),
  m_determinant  (-INFINITY),
  m_lnDeterminant(-INFINITY),
  v_pivoting     (NULL),
  m_signum       (0),
  m_isSingular   (false)
{
  m_mat.shape(map.NumGlobalElements(),nCols);
  m_LU.shape(0,0);
}

QUESO::TeuchosMatrix::TeuchosMatrix	(	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 67 of file TeuchosMatrix.C.

References m_LU, m_mat, and QUESO::Map::NumGlobalElements().

  :
  Matrix  (env,map),
  m_inverse      (NULL),
  m_svdColMap    (NULL),
  m_svdUmat      (NULL),
  m_svdSvec      (NULL),
  m_svdVmat      (NULL),
  m_svdVTmat     (NULL),
  m_determinant  (-INFINITY),
  m_lnDeterminant(-INFINITY),
  v_pivoting     (NULL),
  m_signum       (0),
  m_isSingular   (false)
{
  m_mat.shape    (map.NumGlobalElements(),map.NumGlobalElements());
  m_LU.shape(0,0);

  for (unsigned int i = 0; i < (unsigned int) m_mat.numRows(); ++i) {
    m_mat(i,i) = diagValue;
  }
}

QUESO::TeuchosMatrix::TeuchosMatrix	(	const TeuchosVector &	v,
		double	diagValue
	)

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

Definition at line 95 of file TeuchosMatrix.C.

References m_LU, m_mat, and QUESO::TeuchosVector::sizeLocal().

  :
  Matrix  (v.env(),v.map()),
  m_inverse      (NULL),
  m_svdColMap    (NULL),
  m_svdUmat      (NULL),
  m_svdSvec      (NULL),
  m_svdVmat      (NULL),
  m_svdVTmat     (NULL),
  m_determinant  (-INFINITY),
  m_lnDeterminant(-INFINITY),
  v_pivoting     (NULL),
  m_signum       (0),
  m_isSingular   (false)
{
 m_mat.shape    (v.sizeLocal(),v.sizeLocal());
 m_LU.shape(0,0);

 for (unsigned int i = 0; i < (unsigned int) m_mat.numRows(); ++i) {
    m_mat(i,i) = diagValue;
  }
}

QUESO::TeuchosMatrix::TeuchosMatrix

(

const TeuchosVector &

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 123 of file TeuchosMatrix.C.

References dim, m_LU, m_mat, and QUESO::TeuchosVector::sizeLocal().

  :
  Matrix  (v.env(),v.map()),
  m_inverse      (NULL),
  m_svdColMap    (NULL),
  m_svdUmat      (NULL),
  m_svdSvec      (NULL),
  m_svdVmat      (NULL),
  m_svdVTmat     (NULL),
  m_determinant  (-INFINITY),
  m_lnDeterminant(-INFINITY),
  v_pivoting     (NULL),
  m_signum       (0),
  m_isSingular   (false)
{
  m_mat.shape    (v.sizeLocal(),v.sizeLocal());
  m_LU.shape(0,0);

  unsigned int dim =  v.sizeLocal();

  for (unsigned int i = 0; i < dim; ++i) {
    m_mat(i,i) = v[i];
  }
}

QUESO::TeuchosMatrix::TeuchosMatrix

(

const TeuchosMatrix &

B

)

Shaped Constructor: this matrix is a copy of matrix B.

Definition at line 150 of file TeuchosMatrix.C.

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

  :
  Matrix  (B.env(),B.map()),
  m_inverse      (NULL),
  m_svdColMap    (NULL),
  m_svdUmat      (NULL),
  m_svdSvec      (NULL),
  m_svdVmat      (NULL),
  m_svdVTmat     (NULL),
  m_determinant  (-INFINITY),
  m_lnDeterminant(-INFINITY),
  v_pivoting     (NULL),
  m_signum       (0),
  m_isSingular   (false)
{
  m_mat.shape    (B.numRowsLocal(),B.numCols());
  m_LU.shape(0,0);

  this->Matrix::base_copy(B);
  this->copy(B);
}

QUESO::TeuchosMatrix::~TeuchosMatrix

(

)

Destructor.

Definition at line 174 of file TeuchosMatrix.C.

References resetLU().

{
  this->resetLU();
}

QUESO::TeuchosMatrix::TeuchosMatrix ( )

private

Default Constructor.

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

Referenced by internalSvd(), inverse(), and svdSolve().

Member Function Documentation

int QUESO::TeuchosMatrix::chol ( )

virtual

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

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

Implements QUESO::Matrix.

Definition at line 568 of file TeuchosMatrix.C.

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

{
  int return_success =0 ;
/*  If UPLO = 'L', the leading N-by-N lower triangular part of A contains the lower
 *          triangular part of the matrix A (lets call it L), and the strictly upper
 *          triangular part of A is not referenced. Thus, the upper triangular part
 *          of A +++must be manually+++ overwritten with L^T. */

  Teuchos::LAPACK<int, double> lapack;
  int info;//= 0:  successful exit
           //< 0:  if INFO = -i, the i-th argument had an illegal value
           //> 0:  if INFO = i, the leading minor of order i is not
           //      positive definite, and the factorization could not be
           //      completed.
  char UPLO = 'U';
          //= 'U':  Upper triangle of A is stored;
          //= 'L':  Lower triangle of A is stored.

  lapack.POTRF (UPLO, m_mat.numRows(), m_mat.values(), m_mat.stride(), &info);

  // Overwriting the upper triangular part of the input matrix A with  L^T
  //(the diagonal terms are identical for both L and L^T)

  for (int i=0;i<m_mat.numRows();i++){
    for (int j=i+1;j<m_mat.numCols();j++)
       m_mat(i,j) = m_mat(j,i) ;
  }

  if (info != 0) {
    std::cerr << "In TeuchosMtrix::chol()"
              << ": INFO = " << info
              << ",\n INFO < 0:  if INFO = -i, the i-th argument had an illegal value."
              << "\n INFO > 0:  if INFO = i, the leading minor of order i is not "
              << " positive definite, and the factorization could not be completed."
              << std::endl;
    return_success =1 ;
  }

  UQ_RC_MACRO(info, // Yes, *not* a fatal check on RC
              m_env.worldRank(),
              "TeuchosMatrix::chol()",
              "matrix is not positive definite",
              UQ_MATRIX_IS_NOT_POS_DEFINITE_RC);

  return return_success;
};

void QUESO::TeuchosMatrix::copy ( const TeuchosMatrix &  src )

private

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

Definition at line 1931 of file TeuchosMatrix.C.

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

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

{
  // FIXME - should we be calling Matrix::base_copy here? - RHS
  this->resetLU();
  unsigned int i,j, nrows=src.numRowsLocal(), ncols=src.numCols();

  for(i=0; i< nrows ; i++)
    for (j = 0; j < ncols; j++)
      m_mat(i,j) = src(i,j);

  return;
}

void QUESO::TeuchosMatrix::cwSet

(

double

value

)

Component-wise set all values to this with value.

Definition at line 1189 of file TeuchosMatrix.C.

References m_mat.

{
  m_mat.putScalar(value);
  return;
}

void QUESO::TeuchosMatrix::cwSet	(	unsigned int	rowId,
		unsigned int	colId,
		const TeuchosMatrix &	mat
	)

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

Definition at line 1198 of file TeuchosMatrix.C.

References numCols(), and numRowsLocal().

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

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

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

  queso_require_less_equal_msg((colId + 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)(rowId+i,colId+j) = mat(i,j);
    }
  }

  return;
}

double QUESO::TeuchosMatrix::determinant

(

)

const

Calculates the determinant of this matrix.

Definition at line 442 of file TeuchosMatrix.C.

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

{
  if (m_determinant == -INFINITY)
  {
    if(m_LU.numRows() ==0 && m_LU.numCols() ==0)  //dummy
    {
      TeuchosVector tmpB(m_env,m_map);
      TeuchosVector tmpX(m_env,m_map);
      if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 99)) {
        *m_env.subDisplayFile() << "In TeuchosMatrix::lnDeterminant()"
                                << ": before 'this->invertMultiply()'"
                                << std::endl;
      }
      this->invertMultiply(tmpB,tmpX);
      if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 99)) {
        *m_env.subDisplayFile() << "In TeuchosMatrix::lnDeterminant()"
                                << ": after 'this->invertMultiply()'"
                                << std::endl;
      }
    }
    if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 99)) {
      *m_env.subDisplayFile() << "In TeuchosMatrix::lnDeterminant()"
                              << ": before computing det"
                              << std::endl;
    }

    double det   = 1.0;
    double lnDet = 0.0;
    for (int i=0;i<m_LU.numCols();i++) {
      det   *= m_LU(i,i);
      lnDet += std::log(m_LU(i,i));
    }

    m_determinant   = det;
    m_lnDeterminant = lnDet;

    if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 99)) {
      *m_env.subDisplayFile() << "In TeuchosMatrix::lnDeterminant()"
                              << ": after computing det"
                              << std::endl;
    }
  }

  return m_determinant;
}

void QUESO::TeuchosMatrix::eigen	(	TeuchosVector &	eigenValues,
		TeuchosMatrix *	eigenVectors
	)		const

This function computes the eigenvalues of a real symmetric matrix.

Definition at line 1039 of file TeuchosMatrix.C.

References m_mat, numRowsLocal(), QUESO::queso_require_equal_to_msg, QUESO::queso_require_not_equal_to_msg, and QUESO::TeuchosVector::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");
  }

  // At the end of execution, Lapack funcion destroys the input matrix
  // So, lets not use m_mat, but instead, a copy of it
  Teuchos::SerialDenseMatrix<int,double> copy_m_mat;
  copy_m_mat = m_mat;

  Teuchos::LAPACK<int, double> lapack;
  int lwork = 3*n -1;
  int info;
  char UPLO = 'L';
  double *W, *WORK;

  W = new double[n];
  WORK = new double[lwork];

  if (eigenVectors == NULL) {
    char JOBZ = 'N';//eigenvalues only

    lapack.SYEV(JOBZ, UPLO, copy_m_mat.numRows(), copy_m_mat.values(), copy_m_mat.stride(), &W[0], &WORK[0], lwork, &info);

    for (unsigned int i=0; i< n; i++)
      eigenValues[i] = W[i];

    queso_require_equal_to_msg(info, 0, "invalid input vector size");
  }
  else {
        char JOBZ = 'V';//eigenvalues and eigenvectors

    lapack.SYEV(JOBZ, UPLO, copy_m_mat.numRows(), copy_m_mat.values(), copy_m_mat.stride(), &W[0], &WORK[0], lwork, &info);

        for (unsigned int i=0; i< n; i++)
        eigenValues[i] = W[i];

        eigenVectors->m_mat = copy_m_mat;
  }

  if (info != 0) {
    std::cerr << "In TeuchosMtrix::eigen()"
              << ": INFO = " << info
              << ",\n INFO < 0:  if INFO = -i, the i-th argument had an illegal value."
              << "\n INFO > 0:  if INFO = i, the algorithm failed to converge; i off-diagonal "
                      << " elements of an intermediate tridiagonal form did not converge to zero."
              << std::endl;
  }

  return;
}

void QUESO::TeuchosMatrix::fillWithBlocksDiagonally

(

const std::vector< const TeuchosMatrix * > &

matrices

)

This function fills this matrix diagonally with const block matrices.

Definition at line 1439 of file TeuchosMatrix.C.

References numCols(), numRowsLocal(), and QUESO::queso_require_equal_to_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_equal_to_msg(this->numRowsLocal(), sumNumRowsLocals, "inconsistent local number of rows");
  queso_require_equal_to_msg(this->numCols(), 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)(cumulativeRowId + rowId, cumulativeColId + colId) = (*(matrices[i]))(rowId,colId);
      }
    }
    cumulativeRowId += nRows;
    cumulativeColId += nCols;
  }

  return;
}

void QUESO::TeuchosMatrix::fillWithBlocksDiagonally

(

const std::vector< TeuchosMatrix * > &

matrices

)

This function fills this matrix diagonally with block matrices.

Definition at line 1469 of file TeuchosMatrix.C.

References numCols(), numRowsLocal(), and QUESO::queso_require_equal_to_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_equal_to_msg(this->numRowsLocal(), sumNumRowsLocals, "inconsistent local number of rows");
  queso_require_equal_to_msg(this->numCols(), 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)(cumulativeRowId + rowId, cumulativeColId + colId) = (*(matrices[i]))(rowId,colId);
      }
    }
    cumulativeRowId += nRows;
    cumulativeColId += nCols;
  }

  return;
}

void QUESO::TeuchosMatrix::fillWithBlocksHorizontally

(

const std::vector< const TeuchosMatrix * > &

matrices

)

This function fills this matrix horizontally with const block matrices.

Definition at line 1499 of file TeuchosMatrix.C.

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

{
  unsigned int sumNumCols = 0;
  for (unsigned int i = 0; i < matrices.size(); ++i) {
    queso_require_equal_to_msg(this->numRowsLocal(), matrices[i]->numRowsLocal(), "inconsistent local number of rows");
    sumNumCols += matrices[i]->numCols();
  }
  queso_require_equal_to_msg(this->numCols(), 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)(rowId, cumulativeColId + colId) = (*(matrices[i]))(rowId,colId);
      }
    }
    cumulativeColId += nCols;
  }

  return;
}

void QUESO::TeuchosMatrix::fillWithBlocksHorizontally

(

const std::vector< TeuchosMatrix * > &

matrices

)

This function fills this matrix horizontally with block matrices.

Definition at line 1526 of file TeuchosMatrix.C.

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

{
  unsigned int sumNumCols = 0;
  for (unsigned int i = 0; i < matrices.size(); ++i) {
    queso_require_equal_to_msg(this->numRowsLocal(), matrices[i]->numRowsLocal(), "inconsistent local number of rows");
    sumNumCols += matrices[i]->numCols();
  }
  queso_require_equal_to_msg(this->numCols(), 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)(rowId, cumulativeColId + colId) = (*(matrices[i]))(rowId,colId);
      }
    }
    cumulativeColId += nCols;
  }

  return;
}

void QUESO::TeuchosMatrix::fillWithBlocksVertically

(

const std::vector< const TeuchosMatrix * > &

matrices

)

This function fills this matrix vertically with const block matrices.

Definition at line 1553 of file TeuchosMatrix.C.

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

{
  unsigned int sumNumRows = 0;
  for (unsigned int i = 0; i < matrices.size(); ++i) {
    queso_require_equal_to_msg(this->numCols(), matrices[i]->numCols(), "inconsistent local number of cols");
    sumNumRows += matrices[i]->numRowsLocal();
  }
  queso_require_equal_to_msg(this->numRowsLocal(), 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)(cumulativeRowId + rowId, colId) = (*(matrices[i]))(rowId,colId);
      }
    }
    cumulativeRowId += nRows;
  }
  return;
}

void QUESO::TeuchosMatrix::fillWithBlocksVertically

(

const std::vector< TeuchosMatrix * > &

matrices

)

This function fills this matrix vertically with block matrices.

Definition at line 1579 of file TeuchosMatrix.C.

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

{
  unsigned int sumNumRows = 0;
  for (unsigned int i = 0; i < matrices.size(); ++i) {
    queso_require_equal_to_msg(this->numCols(), matrices[i]->numCols(), "inconsistent local number of cols");
    sumNumRows += matrices[i]->numRowsLocal();
  }
  queso_require_equal_to_msg(this->numRowsLocal(), 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)(cumulativeRowId + rowId, colId) = (*(matrices[i]))(rowId,colId);
      }
    }
    cumulativeRowId += nRows;
  }
  return;
}

void QUESO::TeuchosMatrix::fillWithTensorProduct	(	const TeuchosMatrix &	mat1,
		const TeuchosMatrix &	mat2
	)

This function calculates the tensor product of matrices mat1 and mat2 and stores it in this matrix.

Definition at line 1605 of file TeuchosMatrix.C.

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

{
  queso_require_equal_to_msg(this->numRowsLocal(), (mat1.numRowsLocal() * mat2.numRowsLocal()), "inconsistent local number of rows");
  queso_require_equal_to_msg(this->numCols(), (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)(targetRowId + rowId2, targetColId + colId2) = multiplicativeFactor * mat2(rowId2,colId2);
        }
      }
    }
  }
  return;
}

void QUESO::TeuchosMatrix::fillWithTensorProduct	(	const TeuchosMatrix &	mat1,
		const TeuchosVector &	vec2
	)

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

Definition at line 1628 of file TeuchosMatrix.C.

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

{
  queso_require_equal_to_msg(this->numRowsLocal(), (mat1.numRowsLocal() * vec2.sizeLocal()), "inconsistent local number of rows");
  queso_require_equal_to_msg(this->numCols(), (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)(targetRowId + rowId2, targetColId + colId2) = multiplicativeFactor * vec2[rowId2];
        }
      }
    }
  }
  return;
}

void QUESO::TeuchosMatrix::fillWithTranspose

(

const TeuchosMatrix &

mat

)

This function stores the transpose of this matrix into this matrix.

Definition at line 1420 of file TeuchosMatrix.C.

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

{
  unsigned int nRows = mat.numRowsLocal();
  unsigned int nCols = mat.numCols();
  queso_require_equal_to_msg(this->numRowsLocal(), nCols, "inconsistent local number of rows");
  queso_require_equal_to_msg(this->numCols(), nRows, "inconsistent number of cols");

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

  return;
}

void QUESO::TeuchosMatrix::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 1397 of file TeuchosMatrix.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::TeuchosMatrix::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 1375 of file TeuchosMatrix.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::TeuchosMatrix::getColumn	(	const unsigned int	column_num,
		TeuchosVector &	column
	)		const

This function gets the column_num-th column of this matrix and stores it into vector column.

Definition at line 1220 of file TeuchosMatrix.C.

References m_mat, numCols(), numRowsLocal(), QUESO::queso_require_equal_to_msg, and QUESO::TeuchosVector::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 pointer
  const double* temp_ptr ;

  // get the column_num- of matrix m_mat
  temp_ptr = m_mat[column_num];

  // Copy column from Teuchos matrix into our TeuchosVector object
  for (unsigned int i=0; i< column.sizeLocal();i++)
    column[i] = temp_ptr[i];

  return;
}

TeuchosVector QUESO::TeuchosMatrix::getColumn

(

const unsigned int

column_num

)

const

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

Definition at line 1242 of file TeuchosMatrix.C.

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

{
  // We rely on the void getColumn's to do sanity checks
  // So we don't do them here.
  TeuchosVector column(m_env, m_map);
  this->getColumn( column_num, column );
  return column;
}

void QUESO::TeuchosMatrix::getRow	(	const unsigned int	row_num,
		TeuchosVector &	row
	)		const

This function gets the row_num-th column of this matrix and stores it into vector row.

Definition at line 1272 of file TeuchosMatrix.C.

References m_mat, numRowsLocal(), QUESO::queso_require_equal_to_msg, and QUESO::TeuchosVector::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->numRowsLocal(), "row vector not same size as this matrix");

  // Copy row from Teuchos matrix into our TeuchosVector object
  for (unsigned int i=0; i< row.sizeLocal();i++)
    row[i] = m_mat(row_num,i);

  return;
}

TeuchosVector QUESO::TeuchosMatrix::getRow

(

const unsigned int

row_num

)

const

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

Definition at line 1289 of file TeuchosMatrix.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.
  TeuchosVector row(m_env, m_map);
  this->getRow( row_num, row );
  return row;
}

int QUESO::TeuchosMatrix::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.

This function uses Teuchos GESVD computes the singular value decomposition (SVD) of a real M-by-N matrix A, optionally computing the left and/or right singular vectors. The SVD is written A = U * SIGMA * transpose(V), where SIGMA is an M-by-N matrix which is zero except for its min(m,n) diagonal elements, U is an M-by-M orthogonal matrix, and V is an N-by-N orthogonal matrix. The diagonal elements of SIGMA are the singular values of A; they are real and non-negative, and are returned in descending order. The first min(m,n) columns of U and V are the left and right singular vectors of A. Note that the routine returns V**T, not V.

Definition at line 2015 of file TeuchosMatrix.C.

References QUESO::Matrix::m_env, m_mat, m_svdColMap, m_svdSvec, m_svdUmat, m_svdVmat, m_svdVTmat, QUESO::Matrix::map(), numCols(), numRowsLocal(), stride(), TeuchosMatrix(), values(), and QUESO::TeuchosVector::values().

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

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

    m_svdColMap = new Map(this->numCols(),0,this->map().Comm()); // see 'VectorSpace<.,.>::newMap()'
  //in src/basic/src/TeuchosVectorSpace.C //old comment already existent in GslMatrix
    m_svdUmat   = new TeuchosMatrix(*this); // Yes, 'this'
    m_svdSvec   = new TeuchosVector(m_env,*m_svdColMap);
    m_svdVmat   = new TeuchosMatrix(*m_svdSvec);
    m_svdVTmat  = new TeuchosMatrix(*m_svdSvec);

    int minRowsCols, maxRowsCols;

    if (nRows>=nCols) { minRowsCols = nCols; maxRowsCols = nRows; } else { minRowsCols = nRows; maxRowsCols = nCols; }

    char jobu, jobvt;
    int  lwork, info;
    double  *work, *rwork;

    jobu = 'S';
    jobvt = 'S';

    lwork = 15*maxRowsCols; // Set up the work array, larger than needed.
    work = new double[lwork];

    int aux1= 5*minRowsCols+7, aux2= 2*maxRowsCols+2*minRowsCols+1;
    int aux_dim;

    if (aux1>=aux2) { aux_dim = minRowsCols*aux1; } else {aux_dim = minRowsCols*aux2; }

    rwork = new double[aux_dim];

    Teuchos::LAPACK<int, double> lapack;

    lapack.GESVD(jobu,jobvt,m_mat.numRows(),m_mat.numCols(),m_mat.values(),m_mat.stride(),
                m_svdSvec->values(),m_svdUmat->values(),m_svdUmat->stride(),m_svdVTmat->values(),
                m_svdVTmat->stride(),&work[0],lwork,&rwork[0],&info);
  }
  return 0;
}

TeuchosMatrix QUESO::TeuchosMatrix::inverse

(

)

const

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

Definition at line 392 of file TeuchosMatrix.C.

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

{
  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 TeuchosMatrix(m_env,m_map,nCols);
    TeuchosVector unitVector(m_env,m_map);
    unitVector.cwSet(0.);
    TeuchosVector 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 TeuchosMatrix::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 TeuchosMatrix::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;
  }

/* cout << "Inverse ---------------------\n";
   for (int i=0; i<m_inverse->numRowsLocal(); i++) {
     for (int j=0; j<m_inverse->numCols(); j++)
       cout<< m_inverse->m_mat(i,j) << " ";
     cout << endl;
   }
   cout << endl;
*/
  return *m_inverse;
}

void QUESO::TeuchosMatrix::invertMultiply	(	const TeuchosVector &	b,
		TeuchosVector &	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 private attribute m_LU != NULL .

Definition at line 791 of file TeuchosMatrix.C.

References QUESO::BaseEnvironment::displayVerbosity(), QUESO::Matrix::m_env, QUESO::Matrix::m_inDebugMode, m_isSingular, m_LU, m_mat, QUESO::TeuchosVector::norm2(), numCols(), QUESO::queso_require_equal_to_msg, QUESO::queso_require_not_equal_to_msg, QUESO::TeuchosVector::sizeLocal(), QUESO::BaseEnvironment::subDisplayFile(), and v_pivoting.

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

{
  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");

  if (m_LU.numCols() == 0 && m_LU.numRows() == 0)
  {
  queso_require_msg(!(v_pivoting), "v_pivoting should be NULL");

  //allocate m_LU and v_pivoting
  m_LU = m_mat;
  v_pivoting =(int *) malloc(sizeof(int)*m_LU.numCols() );

  queso_require_not_equal_to_msg(m_LU.numCols(), 0 && m_LU.numRows() == 0, "malloc() failed");

  queso_require_msg(v_pivoting, "malloc() failed");

  if (m_inDebugMode) {
    std::cout << "In TeuchosMatrix::invertMultiply()"
              << ": before LU decomposition, m_LU = ";
              for (int i=0;i<3;i++){
                for (int j=0;j<3;j++)
                  cout << m_LU(i,j) <<"\t" ;
                cout << endl;
              }
    }

  // Perform an LU factorization of matrix m_LU. Checked 12/06/12
  Teuchos::LAPACK<int, double> lapack;
  int info;

  lapack.GETRF( m_LU.numRows(), m_LU.numCols(), m_LU.values(), m_LU.stride(), v_pivoting, &info );

  if (info != 0) {
    std::cerr   << "In TeuchosMatrix::invertMultiply()"
                << ", after lapack.GETRF"
                << ": INFO = " << info
                << ",\nINFO < 0:  if INFO = -i, the i-th argument had an illegal value.\n"
                << "INFO > 0:  if INFO = i, U(i,i) is exactly zero. The factorization \n"
                << "has been completed, but the factor U is exactly singular, and division \n"
                << "by zero will occur if it is used to solve a system of equations."
                << std::endl;
  }
  queso_require_msg(!(info), "GETRF() failed");

  if (info >  0)
     m_isSingular = true;

  if (m_inDebugMode) {
    std::cout << "In TeuchosMatrix::invertMultiply()"
              << ": after LU decomposition, m_LU = ";
              for (int i=0;i<3;i++){
                for (int j=0;j<3;j++)
                  cout << m_LU(i,j) <<"\t" ;
                cout << endl;
              }
      std::cout << std::endl;
    }
  }
  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 99)) {
    *m_env.subDisplayFile() << "In TeuchosMatrix::invertMultiply()"
                            << ": before 'lapack.GETRS()'"
                            << std::endl;
  }

  // Solve the linear system.
  Teuchos::LAPACK<int, double> lapack;
  int NRHS = 1; // NRHS: number of right hand sides, i.e., the number
                                // of columns of the matrix B. In this case, vector b.
  char TRANS = 'N';  // 'N':  A * x= B  (No transpose). Specifies the
                                     // form of the system of equations.
  int info02;

  //GETRS expects the matrix to be already factored in LU and uses the
  //same ipiv vector, which are the pivot indices of the LU factorization

  x=b;
  lapack.GETRS(TRANS, m_LU.numRows(), NRHS, m_LU.values(), m_LU.stride(), v_pivoting, &x[0],x.sizeLocal(), &info02 );

  if (info02 != 0) {
      std::cerr << "In TeuchosMatrix::invertMultiply()"
                << ", after lapack.GETRS - solve LU system"
                << ": INFO = " << info02
                << ",\nINFO < 0:  if INFO = -i, the i-th argument had an illegal value.\n"
                << std::endl;
  }
  queso_require_msg(!(info02), "GETRS() failed");
 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 99)) {
    *m_env.subDisplayFile() << "In TeuchosMatrix::invertMultiply()"
                            << ", after lapack.GETRS() - solve LU system."
                            << std::endl;
 }
 if (m_inDebugMode) {
    TeuchosVector tmpVec(b - (*this)*x);
    std::cout << "In TeuchosMatrix::invertMultiply()"
              << ": ||b - Ax||_2 = "         << tmpVec.norm2()
              << ": ||b - Ax||_2/||b||_2 = " << tmpVec.norm2()/b.norm2()
              << std::endl;
 }
 return;
}

TeuchosVector QUESO::TeuchosMatrix::invertMultiply

(

const TeuchosVector &

b

)

const

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

It calls void TeuchosMatrix::invertMultiply(const TeuchosVector& b, TeuchosVector& x) internally.

Definition at line 778 of file TeuchosMatrix.C.

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

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

  this->invertMultiply(b,x);

  return x;
}

void QUESO::TeuchosMatrix::invertMultiply	(	const TeuchosMatrix &	B,
		TeuchosMatrix &	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 private attribute m_LU != NULL .

Definition at line 906 of file TeuchosMatrix.C.

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

{
  // Sanity Checks
  queso_require_msg(!((B.numRowsLocal() != X.numRowsLocal()) || (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.
  TeuchosVector b(m_env, m_map);
  TeuchosVector x(m_env, m_map);

  for( unsigned int j = 0; j < B.numCols(); ++j ) {
    b = B.getColumn(j);
    //invertMultiply will only do the LU factorization once and store it.
    x = this->invertMultiply(b);
    X.setColumn(j,x);
  }
  return;
}

TeuchosMatrix QUESO::TeuchosMatrix::invertMultiply

(

const TeuchosMatrix &

B

)

const

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

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

Definition at line 896 of file TeuchosMatrix.C.

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

{
  TeuchosMatrix X(m_env,m_map,B.numCols());
  this->invertMultiply(B,X);
  return X;
}

void QUESO::TeuchosMatrix::invertMultiplyForceLU	(	const TeuchosVector &	b,
		TeuchosVector &	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 941 of file TeuchosMatrix.C.

References m_isSingular, m_LU, m_mat, numCols(), QUESO::queso_require_equal_to_msg, QUESO::queso_require_not_equal_to_msg, QUESO::TeuchosVector::sizeLocal(), and v_pivoting.

Referenced by invertMultiplyForceLU().

{
  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");

  if (m_LU.numCols() == 0 && m_LU.numRows() == 0) {
    queso_require_msg(!(v_pivoting), "v_pivoting should be NULL");
  }

  //allocate m_LU , yes outside the if above
  m_LU = m_mat;

  queso_require_not_equal_to_msg(m_LU.numCols(), 0 && m_LU.numRows() == 0, "Teuchos atttribuition m_LU = m_mat failed");

  //allocate v_pivoting
  if ( v_pivoting == NULL )
    v_pivoting =(int *) malloc(sizeof(int)*m_LU.numCols() );

  queso_require_msg(v_pivoting, "malloc() for v_pivoting failed");

  // Perform an LU factorization of matrix m_LU. Checked 12/06/12
  Teuchos::LAPACK<int, double> lapack;
  int info;

  lapack.GETRF( m_LU.numRows(), m_LU.numCols(), m_LU.values(), m_LU.stride(), v_pivoting, &info );

  if (info != 0) {
    std::cerr   << "In TeuchosMatrix::invertMultiply()"
                << ", after lapack.GETRF"
                << ": INFO = " << info
                << ",\nINFO < 0:  if INFO = -i, the i-th argument had an illegal value.\n"
                << "INFO > 0:  if INFO = i, U(i,i) is exactly zero. The factorization \n"
                << "has been completed, but the factor U is exactly singular, and division \n"
                << "by zero will occur if it is used to solve a system of equations."
                << std::endl;
  }

  queso_require_msg(!(info), "GETRF() failed");

  if (info >  0)
     m_isSingular = true;

  // Solve the linear system.
  int NRHS = 1; // NRHS: number of right hand sides, i.e., the number
                                // of columns of the matrix B. In this case, vector b.
  char TRANS = 'N';  // 'N':  A * x= B  (No transpose). Specifies the
                                     // form of the system of equations.
  int info02;

  //GETRS expects the matrix to be already factored in LU and uses the
  //same ipiv vector, which are the pivot indices of the LU factorization

  x=b;
  lapack.GETRS(TRANS, m_LU.numRows(), NRHS, m_LU.values(), m_LU.stride(), v_pivoting, &x[0],x.sizeLocal(), &info02 );

  if (info02 != 0) {
      std::cerr << "In TeuchosMatrix::invertMultiply()"
                << ", after lapack.GETRS - solve LU system"
                << ": INFO = " << info02
                << ",\nINFO < 0:  if INFO = -i, the i-th argument had an illegal value.\n"
                << std::endl;
    }

    queso_require_msg(!(info02), "GETRS() failed");

 return;
}

TeuchosVector QUESO::TeuchosMatrix::invertMultiplyForceLU

(

const TeuchosVector &

b

)

const

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

  It calls void TeuchosMatrix::InvertMultiplyForceLU(const TeuchosVector& b, TeuchosVector& x) const;(const TeuchosVector& b,

TeuchosVector& x) internally.

Definition at line 928 of file TeuchosMatrix.C.

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

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

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

  return x;
}

void QUESO::TeuchosMatrix::largestEigen	(	double &	eigenValue,
		TeuchosVector &	eigenVector
	)		const

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

Definition at line 1097 of file TeuchosMatrix.C.

References m_mat, QUESO::queso_require_not_equal_to_msg, and QUESO::TeuchosVector::sizeLocal().

{
  unsigned int n = eigenVector.sizeLocal();

  queso_require_not_equal_to_msg(n, 0, "invalid input vector size");
  Teuchos::LAPACK<int, double> lapack;

  //SYEV destroys the input matrix, so lets operate in a copy
  Teuchos::SerialDenseMatrix<int,double> copy_m_mat;
  copy_m_mat = m_mat;

  int lwork = 3*n -1;
  int info;
  char UPLO = 'L';
  char JOBZ = 'V'; //eigenvalues and eigenvectors
  double *W, *WORK;

  W = new double[n];
  WORK = new double[lwork];

  lapack.SYEV(JOBZ, UPLO, copy_m_mat.numRows(), copy_m_mat.values(), copy_m_mat.stride(), &W[0], &WORK[0], lwork, &info);

  if (info != 0) {
    std::cerr << "In TeuchosMtrix::largestEigen()"
              << ": INFO = " << info
              << ",\n INFO < 0:  if INFO = -i, the i-th argument had an illegal value."
              << "\n INFO > 0:  if INFO = i, the algorithm failed to converge; i off-diagonal "
                      << " elements of an intermediate tridiagonal form did not converge to zero."
              << std::endl;
  }
  // If INFO = 0, W contains the eigenvalues in ascending order.
  // Thus the largest eigenvalue is in W[n-1].
  eigenValue = W[n-1];

  // Eigenvector associated to the largest eigenvalue.
  // Stored in the n-th column of matrix copy_m_mat.
  for (int i=0; i< copy_m_mat.numRows(); i++)
    eigenVector[i] = copy_m_mat(i,n-1);

  return;
}

double QUESO::TeuchosMatrix::lnDeterminant

(

)

const

Calculates the ln(determinant) of this matrix.

Definition at line 491 of file TeuchosMatrix.C.

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

Referenced by inverse().

{
  if (m_lnDeterminant == -INFINITY)
  {
    if(m_LU.numRows() ==0 && m_LU.numCols() ==0)  //dummy
    {
      TeuchosVector tmpB(m_env,m_map);
      TeuchosVector tmpX(m_env,m_map);
      if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 99)) {
        *m_env.subDisplayFile() << "In TeuchosMatrix::lnDeterminant()"
                                << ": before 'this->invertMultiply()'"
                                << std::endl;
      }
      this->invertMultiply(tmpB,tmpX);
      if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 99)) {
        *m_env.subDisplayFile() << "In TeuchosMatrix::lnDeterminant()"
                                << ": after 'this->invertMultiply()'"
                                << std::endl;
      }
    }
    if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 99)) {
      *m_env.subDisplayFile() << "In TeuchosMatrix::lnDeterminant()"
                              << ": before computing lnDet"
                              << std::endl;
    }

    double det   = 1.0;
    double lnDet = 0.0;
    for (int i=0;i<m_LU.numCols();i++) {
      det   *= m_LU(i,i);
      lnDet += std::log(m_LU(i,i));
    }

    m_determinant   = det;
    m_lnDeterminant = lnDet;

    if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 99)) {
      *m_env.subDisplayFile() << "In TeuchosMatrix::lnDeterminant()"
                              << ": after computing lnDet"
                              << std::endl;
    }
  }

  return m_lnDeterminant;
}

void QUESO::TeuchosMatrix::matlabLinearInterpExtrap	(	const TeuchosVector &	x1Vec,
		const TeuchosMatrix &	y1Mat,
		const TeuchosVector &	x2Vec
	)

Definition at line 1693 of file TeuchosMatrix.C.

References getColumn(), QUESO::TeuchosVector::matlabLinearInterpExtrap(), numCols(), numRowsLocal(), QUESO::queso_require_equal_to_msg, setColumn(), and QUESO::TeuchosVector::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");

  TeuchosVector y1Vec(x1Vec);
  TeuchosVector 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::TeuchosMatrix::max

(

)

const

Returns the maximum element value of the matrix.

Definition at line 315 of file TeuchosMatrix.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::TeuchosMatrix::mpiSum	(	const MpiComm &	comm,
		TeuchosMatrix &	M_global
	)		const

Definition at line 1653 of file TeuchosMatrix.C.

References QUESO::MpiComm::Allreduce(), numCols(), numRowsLocal(), QUESO::queso_require_equal_to_msg, and QUESO::size.

{
  // Sanity Checks
  queso_require_equal_to_msg(this->numRowsLocal(),
                             M_global.numRowsLocal(),
                             "local and global matrices incompatible");
  queso_require_equal_to_msg(this->numCols(),
                             M_global.numCols(),
                             "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,
                 "TeuchosMatrix::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;
}

TeuchosVector QUESO::TeuchosMatrix::multiply

(

const TeuchosVector &

x

)

const

This function multiplies this matrix by vector x and returns a vector.

Definition at line 765 of file TeuchosMatrix.C.

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

Referenced by QUESO::operator*().

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

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

  return y;
}

void QUESO::TeuchosMatrix::multiply ( const TeuchosVector &  x, TeuchosVector &  y  ) const

private

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

Definition at line 1993 of file TeuchosMatrix.C.

References numCols(), numRowsLocal(), QUESO::queso_require_equal_to_msg, and QUESO::TeuchosVector::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::TeuchosMatrix::normFrob

(

)

const

Returns the Frobenius norm of this matrix.

Definition at line 541 of file TeuchosMatrix.C.

References m_mat.

{
  return m_mat.normFrobenius ();
}

double QUESO::TeuchosMatrix::normMax

(

)

const

Returns the Frobenius norm of this matrix.

Definition at line 548 of file TeuchosMatrix.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::TeuchosMatrix::numCols ( ) const

virtual

Returns the column dimension of this matrix.

Implements QUESO::Matrix.

Definition at line 290 of file TeuchosMatrix.C.

References m_mat.

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

{
  return (unsigned int) m_mat.numCols();
};

unsigned int QUESO::TeuchosMatrix::numRowsGlobal ( ) const

virtual

Returns the global row dimension of this matrix.

Implements QUESO::Matrix.

Definition at line 282 of file TeuchosMatrix.C.

References m_mat.

{
  return (unsigned int) m_mat.numRows();
}

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

virtual

Returns the local row dimension of this matrix.

Implements QUESO::Matrix.

Definition at line 274 of file TeuchosMatrix.C.

References m_mat.

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

{
  return (unsigned int) m_mat.numRows();
}

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

Element access method (non-const).

Definition at line 244 of file TeuchosMatrix.C.

References m_mat, and resetLU().

{
  this->resetLU();
  if ((i >= (unsigned int) m_mat.numRows()) ||
      (j >= (unsigned int) m_mat.numCols())) {
    std::cerr << "In TeuchosMatrix::operator(i,j) (non-const)"
              << ": i = " << i
              << ", j = " << j
              << ", m_mat.numRows() = " << (unsigned int) m_mat.numRows()
              << ", m_mat.numCols() = " << (unsigned int) m_mat.numCols()
              << std::endl;
    queso_require_less_msg(i, (unsigned int) m_mat.numRows(), "i is too large");
    queso_require_less_msg(j, (unsigned int) m_mat.numCols(), "j is too large");
  }
  return m_mat(i,j);
}

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

Element access method (const).

Definition at line 262 of file TeuchosMatrix.C.

References m_mat.

{
  queso_require_less_msg(i, (unsigned int) m_mat.numRows(), "i is too large");
  queso_require_less_msg(j, (unsigned int) m_mat.numCols(), "j is too large");
  return m_mat(i,j);
}

TeuchosMatrix & QUESO::TeuchosMatrix::operator*=

(

double

a

)

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

Definition at line 192 of file TeuchosMatrix.C.

References m_mat, and resetLU().

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

TeuchosMatrix & QUESO::TeuchosMatrix::operator+=

(

const TeuchosMatrix &

rhs

)

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

Definition at line 212 of file TeuchosMatrix.C.

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

{
  this->resetLU();

  unsigned int i,j, nrows=rhs.numRowsLocal(), ncols=rhs.numCols();

  for(i=0; i< nrows ; i++)
    for (j = 0; j < ncols; j++)
      (*this)(i,j) += rhs(i,j);

  return *this;
}

TeuchosMatrix & QUESO::TeuchosMatrix::operator-=

(

const TeuchosMatrix &

rhs

)

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

Definition at line 227 of file TeuchosMatrix.C.

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

{
  this->resetLU();

  unsigned int i,j, nrows=rhs.numRowsLocal(), ncols=rhs.numCols();

  for(i=0; i< nrows ; i++)
    for (j = 0; j < ncols; j++)
      (*this)(i,j) -= rhs(i,j);

  return *this;
}

TeuchosMatrix & QUESO::TeuchosMatrix::operator/=

(

double

a

)

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

Definition at line 203 of file TeuchosMatrix.C.

References m_mat, and resetLU().

{
  this->resetLU();
  m_mat.scale(1./a);
  return *this;
}

TeuchosMatrix & QUESO::TeuchosMatrix::operator=

(

const TeuchosMatrix &

rhs

)

Copies values from matrix rhs to this.

Definition at line 183 of file TeuchosMatrix.C.

References copy(), and resetLU().

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

void QUESO::TeuchosMatrix::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 1721 of file TeuchosMatrix.C.

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

Referenced by QUESO::operator<<().

{
  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::TeuchosMatrix::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 335 of file TeuchosMatrix.C.

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

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

  TeuchosVector 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 TeuchosMatrix::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::TeuchosMatrix::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 1946 of file TeuchosMatrix.C.

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

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

{
  if (m_LU.numCols() >0 || m_LU.numRows() > 0) {
    m_LU.reshape(0,0); //Kemelli, 12/06/12, dummy
  }
  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 (v_pivoting) {  //Kemelli added 12/09/12
    free(v_pivoting);
    v_pivoting = NULL;

  }
  m_signum = 0;
  m_isSingular = false;

  return;
}

void QUESO::TeuchosMatrix::setColumn	(	const unsigned int	column_num,
		const TeuchosVector &	column
	)

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

Definition at line 1255 of file TeuchosMatrix.C.

References m_mat, numCols(), numRowsLocal(), QUESO::queso_require_equal_to_msg, resetLU(), and QUESO::TeuchosVector::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");

  for (unsigned int i =0; i < column.sizeLocal(); i++)
    m_mat(i,column_num) = column[i];

  return;
}

void QUESO::TeuchosMatrix::setRow	(	const unsigned int	row_num,
		const TeuchosVector &	row
	)

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

Definition at line 1301 of file TeuchosMatrix.C.

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

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

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

  // Copy our TeuchosVector object to our Teuchos Matrix
  for (unsigned int i=0; i< row.sizeLocal();i++)
     m_mat(row_num,i) = row[i] ;

  return;
}

void QUESO::TeuchosMatrix::smallestEigen	(	double &	eigenValue,
		TeuchosVector &	eigenVector
	)		const

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

Definition at line 1141 of file TeuchosMatrix.C.

References m_mat, QUESO::queso_require_not_equal_to_msg, and QUESO::TeuchosVector::sizeLocal().

{
  unsigned int n = eigenVector.sizeLocal();

  queso_require_not_equal_to_msg(n, 0, "invalid input vector size");
  Teuchos::LAPACK<int, double> lapack;

  //SYEV destroys the input matrix, so lets operate in a copy
  Teuchos::SerialDenseMatrix<int,double> copy_m_mat;
  copy_m_mat = m_mat;

  int lwork = 3*n -1;
  int info;
  char UPLO = 'L';
  char JOBZ = 'V';//eigenvalues and eigenvectors
  double *W, *WORK;

  W = new double[n];
  WORK = new double[lwork];

  lapack.SYEV(JOBZ, UPLO, copy_m_mat.numRows(), copy_m_mat.values(), copy_m_mat.stride(), &W[0], &WORK[0], lwork, &info);

  if (info != 0) {
    std::cerr << "In TeuchosMtrix::smallestEigen()"
              << ": INFO = " << info
              << ",\n INFO < 0:  if INFO = -i, the i-th argument had an illegal value."
              << "\n INFO > 0:  if INFO = i, the algorithm failed to converge; i off-diagonal "
                      << " elements of an intermediate tridiagonal form did not converge to zero."
              << std::endl;
  }

  // If INFO = 0, W contains the eigenvalues in ascending order.
  // Thus the smallest eigenvalue is in W[0].
  eigenValue = W[0];

  // Eigenvector associated to the smallest eigenvalue.
  // Stored in the n-th column of matrix copy_m_mat.
  for (int i=0; i< copy_m_mat.numRows(); i++)
    eigenVector[i] = copy_m_mat(i,0);

  return;
}

int QUESO::TeuchosMatrix::stride

(

)

Returns the stride between the columns of this matrix in memory.

Definition at line 306 of file TeuchosMatrix.C.

References m_mat.

Referenced by internalSvd(), and svdSolve().

{
  return  m_mat.stride();
};

void QUESO::TeuchosMatrix::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 1751 of file TeuchosMatrix.C.

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

{
  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;

    // Read '=' sign
    *filePtrSet.ifsVar >> tmpString;

    queso_require_equal_to_msg(tmpString, "=", "string should be the '=' sign");

    // Read 'zeros(n_rows,n_cols)' string
    *filePtrSet.ifsVar >> tmpString;

    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 TeuchosMatrix::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 TeuchosMatrix::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 TeuchosMatrix::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::TeuchosMatrix::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 1885 of file TeuchosMatrix.C.

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

{
  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::TeuchosMatrix::svd	(	TeuchosMatrix &	matU,
		TeuchosVector &	vecS,
		TeuchosMatrix &	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 617 of file TeuchosMatrix.C.

References internalSvd(), m_svdSvec, m_svdUmat, m_svdVTmat, numCols(), numRowsLocal(), QUESO::queso_require_equal_to_msg, and QUESO::TeuchosVector::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 TeuchosMatrix & QUESO::TeuchosMatrix::svdMatU

(

)

const

This function calls private member TeuchosMatrix::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 639 of file TeuchosMatrix.C.

References internalSvd(), and m_svdUmat.

Referenced by svdSolve().

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

const TeuchosMatrix & QUESO::TeuchosMatrix::svdMatV

(

)

const

This function calls private member TeuchosMatrix::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 649 of file TeuchosMatrix.C.

References internalSvd(), and m_svdVmat.

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

int QUESO::TeuchosMatrix::svdSolve	(	const TeuchosVector &	rhsVec,
		TeuchosVector &	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 TeuchosMatrix::svd (x=solVec, b=rhsVec).

An orthogonal matrix A has a norm-preserving property, i.e. for any vector v,

\[ ||Av|| = ||v|| \]

Then:

\[ min( ||Ax - b||^2) = min(||Ax - b||) = min(||U D V^T x - b||) = min(||D V x - U b||)\]

Substituting

\[ y = V^T x \]

and

\[ b' = U^T b \]

gives us

\[ Dy = b'\]

with D being a diagonal matrix. Or,

\[ y = inv(D) U^T b \]

and we only have to solve the linear system:

\[ V^T x = y\]

Definition at line 667 of file TeuchosMatrix.C.

References QUESO::BaseEnvironment::displayVerbosity(), internalSvd(), QUESO::Matrix::m_env, m_svdSvec, m_svdUmat, m_svdVmat, m_svdVTmat, numCols(), numRowsLocal(), QUESO::queso_require_equal_to_msg, QUESO::TeuchosVector::sizeLocal(), stride(), QUESO::BaseEnvironment::subDisplayFile(), svdMatU(), TeuchosMatrix(), transpose(), values(), and QUESO::TeuchosVector::values().

Referenced by svdSolve().

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

  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 TeuchosMatrix::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)
 {
   TeuchosMatrix  invD = TeuchosMatrix(solVec);
   TeuchosMatrix  auxMatrix = TeuchosMatrix(solVec);

   for (i=0; i<nRows; i++){
     invD(i,i) = 1./(m_svdSvec->values()[i]);
   }

  // GESV: For a system Ax=b, on entry, b contains the right-hand side b;
  // on exit it contains the solution x.
  // Thus, intead of doing y = inv(D)*UT*rhsVec = auxMatrix*rhsVec, with
  // auxMatrix = inv(D)*UT; and then assigning solVec=y (requirement for using GESV)
  // lets do solVec= auxMatrix*rhsVec, and save one step in the calculation

    auxMatrix = invD * svdMatU().transpose();
    solVec = auxMatrix*rhsVec;

  // solve the linear system VT * solVec = y
  // GESV changes the values of the matrix, so lets make a copy of it and use it.

    TeuchosMatrix* aux_m_svdVTmat  = new TeuchosMatrix(*m_svdVTmat);

    int ipiv[0], info;
    Teuchos::LAPACK<int, double> lapack;
    lapack.GESV(nCols, 1,  aux_m_svdVTmat->values(),  aux_m_svdVTmat->stride(), ipiv, &solVec[0], solVec.sizeLocal(), &info );

  /* GESV output INFO: = 0:  successful exit
   *          < 0:  if INFO = -i, the i-th argument had an illegal value
   *          > 0:  if INFO = i, U(i,i) is exactly zero.  The factorization
   *                has been completed, but the factor U is exactly
   *                singular, so the solution could not be computed.   */

    iRC = info;
    delete aux_m_svdVTmat;
  }
  return iRC;
}

int QUESO::TeuchosMatrix::svdSolve	(	const TeuchosMatrix &	rhsMat,
		TeuchosMatrix &	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 TeuchosMatrix::svd (x=solMat, b=rhsMat).

Note that solMat is a matrix. Thus, to solve the system of equations ' * solMath = rhsMat', this method calls svdSolve(const TeuchosVector& rhsVec, TeuchosVector& solVec) passing one column of solMat with its respective column of rhsMat, .

Definition at line 735 of file TeuchosMatrix.C.

References getColumn(), QUESO::Matrix::m_env, QUESO::Matrix::map(), numCols(), numRowsLocal(), QUESO::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");

  TeuchosVector rhsVec(m_env,rhsMat.map());
  TeuchosVector 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;
}

TeuchosMatrix QUESO::TeuchosMatrix::transpose

(

)

const

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

Definition at line 372 of file TeuchosMatrix.C.

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

Referenced by svdSolve().

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

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

  TeuchosMatrix 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;
}

double * QUESO::TeuchosMatrix::values

(

)

Returns a pointer to the first element of this matrix.

Definition at line 298 of file TeuchosMatrix.C.

References m_mat.

Referenced by internalSvd(), and svdSolve().

{
  return  m_mat.values();
};

void QUESO::TeuchosMatrix::zeroLower ( bool  includeDiagonal = false )

virtual

This function sets all the entries above the main diagonal (inclusive or not) 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 1318 of file TeuchosMatrix.C.

References numCols(), numRowsLocal(), QUESO::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::TeuchosMatrix::zeroUpper ( bool  includeDiagonal = false )

virtual

This function sets all the entries bellow the main diagonal (inclusive or not) 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 1347 of file TeuchosMatrix.C.

References numCols(), numRowsLocal(), QUESO::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::TeuchosMatrix::m_determinant

mutableprivate

The determinant of this matrix.

Definition at line 376 of file TeuchosMatrix.h.

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

TeuchosMatrix* QUESO::TeuchosMatrix::m_inverse

mutableprivate

Stores the inverse of this matrix.

Definition at line 358 of file TeuchosMatrix.h.

Referenced by inverse(), and resetLU().

bool QUESO::TeuchosMatrix::m_isSingular

mutableprivate

Indicates whether or not this matrix is singular.

Definition at line 390 of file TeuchosMatrix.h.

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

double QUESO::TeuchosMatrix::m_lnDeterminant

mutableprivate

The natural logarithm of the determinant of this matrix.

Definition at line 379 of file TeuchosMatrix.h.

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

Teuchos::SerialDenseMatrix<int,double> QUESO::TeuchosMatrix::m_LU

mutableprivate

Teuchos matrix for the LU decomposition of m_mat.

Definition at line 355 of file TeuchosMatrix.h.

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

Teuchos::SerialDenseMatrix<int,double> QUESO::TeuchosMatrix::m_mat

private

Teuchos matrix, also referred to as this matrix.

Definition at line 352 of file TeuchosMatrix.h.

Referenced by chol(), copy(), cwSet(), eigen(), getColumn(), getRow(), internalSvd(), invertMultiply(), invertMultiplyForceLU(), largestEigen(), normFrob(), numCols(), numRowsGlobal(), numRowsLocal(), operator()(), operator*=(), operator/=(), setColumn(), setRow(), smallestEigen(), stride(), TeuchosMatrix(), and values().

TeuchosMatrix* QUESO::TeuchosMatrix::m_permutation

mutableprivate

Definition at line 381 of file TeuchosMatrix.h.

int QUESO::TeuchosMatrix::m_signum

mutableprivate

Definition at line 387 of file TeuchosMatrix.h.

Referenced by resetLU().

Map* QUESO::TeuchosMatrix::m_svdColMap

mutableprivate

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

Definition at line 361 of file TeuchosMatrix.h.

Referenced by internalSvd(), and resetLU().

TeuchosVector* QUESO::TeuchosMatrix::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.

Definition at line 367 of file TeuchosMatrix.h.

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

TeuchosMatrix* QUESO::TeuchosMatrix::m_svdUmat

mutableprivate

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

Definition at line 364 of file TeuchosMatrix.h.

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

TeuchosMatrix* QUESO::TeuchosMatrix::m_svdVmat

mutableprivate

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

Definition at line 370 of file TeuchosMatrix.h.

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

TeuchosMatrix* QUESO::TeuchosMatrix::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.

Definition at line 373 of file TeuchosMatrix.h.

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

int* QUESO::TeuchosMatrix::v_pivoting

mutableprivate

The pivoting vector of a LU decomposition.

Definition at line 384 of file TeuchosMatrix.h.

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

---

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

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