GslMatrix.h

Go to the documentation of this file.

//-----------------------------------------------------------------------bl-
//--------------------------------------------------------------------------
//
// QUESO - a library to support the Quantification of Uncertainty
// for Estimation, Simulation and Optimization
//
// Copyright (C) 2008-2017 The PECOS Development Team
//
// This library is free software; you can redistribute it and/or
// modify it under the terms of the Version 2.1 GNU Lesser General
// Public License as published by the Free Software Foundation.
//
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
// Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc. 51 Franklin Street, Fifth Floor,
// Boston, MA  02110-1301  USA
//
//-----------------------------------------------------------------------el-

#ifndef UQ_GSL_MATRIX_H
#define UQ_GSL_MATRIX_H

#include <queso/Defines.h>
#include <queso/GslVector.h>
#include <queso/Matrix.h>

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_permutation.h>

namespace QUESO {

class GslMatrix : public Matrix
{
public:


  GslMatrix(const BaseEnvironment& env,
                   const Map&             map,
                   unsigned int                  numCols);

  GslMatrix(const BaseEnvironment& env,
                   const Map&             map,
                   double                        diagValue); // MATLAB eye

  GslMatrix(const GslVector&       v,
                   double                        diagValue); // MATLAB eye


  GslMatrix(const GslVector&       v);         // MATLAB diag


  GslMatrix(const GslMatrix&       B);

  ~GslMatrix();



  GslMatrix& operator= (const GslMatrix& rhs);

  GslMatrix& operator*=(double a);

  GslMatrix& operator/=(double a);

  GslMatrix& operator+=(const GslMatrix& rhs);

  GslMatrix& operator-=(const GslMatrix& rhs);



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

  const double& operator()(unsigned int i, unsigned int j) const
  {
    queso_require_less_msg(i, m_mat->size1, "i is too large");
    queso_require_less_msg(j, m_mat->size2, "j is too large");
    return *gsl_matrix_ptr(m_mat,i,j);
  }




  unsigned int      numRowsLocal              () const;

  unsigned int      numRowsGlobal             () const;

  unsigned int      numCols                   () const;

  double            max                       () const;


  unsigned int      rank                      (double absoluteZeroThreshold, double relativeZeroThreshold) const;

  GslMatrix  transpose                 () const;

  GslMatrix  inverse                   () const;


  double            determinant               () const;

  double            lnDeterminant             () const;



  double            normFrob                  () const;

  double            normMax                   () const;








  int               chol                      ();

  int               svd                       (GslMatrix& matU, GslVector& vecS, GslMatrix& matVt) const;


  const GslMatrix& svdMatU                   () const;


  const GslMatrix& svdMatV                   () const;

  int               svdSolve                  (const GslVector& rhsVec, GslVector& solVec) const;

  int               svdSolve                  (const GslMatrix& rhsMat, GslMatrix& solMat) const;



  GslVector  multiply                  (const GslVector& x) const;

  //the preallocated vector \c y.
  void       multiply                  (const GslVector& x, GslVector& y) const;

  GslMatrix  multiply                  (const GslMatrix& X) const;

  //the preallocated zeroed matrix \c Y.
  void       multiply                  (const GslMatrix& X, GslMatrix& Y) const;


  GslVector  invertMultiply            (const GslVector& b) const;


  void              invertMultiply            (const GslVector& b, GslVector& x) const;


  GslMatrix  invertMultiply            (const GslMatrix& B) const;


  void              invertMultiply            (const GslMatrix& B, GslMatrix& X) const;


  GslVector  invertMultiplyForceLU     (const GslVector& b) const;


  void              invertMultiplyForceLU     (const GslVector& b, GslVector& x) const;

  void              getColumn                 (const unsigned int column_num, GslVector& column) const;

  GslVector  getColumn                 (const unsigned int column_num) const;

  void              setColumn                 (const unsigned int column_num, const GslVector& column);

  void              getRow                    (const unsigned int row_num, GslVector& row) const;

  GslVector  getRow                    (const unsigned int row_num) const;

  void              setRow                    (const unsigned int row_num, const GslVector& row);

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

  void              largestEigen              (double& eigenValue, GslVector& eigenVector) const;

  void              smallestEigen             (double& eigenValue, GslVector& eigenVector) const;





  void              cwSet                     (double value);

  void              cwSet                     (unsigned int rowId, unsigned int colId, const GslMatrix& mat);

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


  void              zeroLower                 (bool includeDiagonal = false);


  void              zeroUpper                 (bool includeDiagonal = false);


  void              filterSmallValues         (double thresholdValue);


  void              filterLargeValues         (double thresholdValue);

  void              fillWithTranspose         (unsigned int            rowId,
                                               unsigned int            colId,
                                               const GslMatrix& mat,
                                               bool                    checkForExactNumRowsMatching,
                                               bool                    checkForExactNumColsMatching);

  void              fillWithBlocksDiagonally  (unsigned int                                 rowId,
                                               unsigned int                                 colId,
                                               const std::vector<const GslMatrix* >& matrices,
                                               bool                                         checkForExactNumRowsMatching,
                                               bool                                         checkForExactNumColsMatching);

  void              fillWithBlocksDiagonally  (unsigned int                                 rowId,
                                               unsigned int                                 colId,
                                               const std::vector<      GslMatrix* >& matrices,
                                               bool                                         checkForExactNumRowsMatching,
                                               bool                                         checkForExactNumColsMatching);

  void              fillWithBlocksHorizontally(unsigned int                                 rowId,
                                               unsigned int                                 colId,
                                               const std::vector<const GslMatrix* >& matrices,
                                               bool                                         checkForExactNumRowsMatching,
                                               bool                                         checkForExactNumColsMatching);

  void              fillWithBlocksHorizontally(unsigned int                                 rowId,
                                               unsigned int                                 colId,
                                               const std::vector<      GslMatrix* >& matrices,
                                               bool                                         checkForExactNumRowsMatching,
                                               bool                                         checkForExactNumColsMatching);

  void              fillWithBlocksVertically  (unsigned int                                 rowId,
                                               unsigned int                                 colId,
                                               const std::vector<const GslMatrix* >& matrices,
                                               bool                                         checkForExactNumRowsMatching,
                                               bool                                         checkForExactNumColsMatching);

  void              fillWithBlocksVertically  (unsigned int                                 rowId,
                                               unsigned int                                 colId,
                                               const std::vector<      GslMatrix* >& matrices,
                                               bool                                         checkForExactNumRowsMatching,
                                               bool                                         checkForExactNumColsMatching);

  void              fillWithTensorProduct     (unsigned int            rowId,
                                               unsigned int            colId,
                                               const GslMatrix& mat1,
                                               const GslMatrix& mat2,
                                               bool                    checkForExactNumRowsMatching,
                                               bool                    checkForExactNumColsMatching);

  void              fillWithTensorProduct     (unsigned int            rowId,
                                               unsigned int            colId,
                                               const GslMatrix& mat1,
                                               const GslVector& vec2,
                                               bool                    checkForExactNumRowsMatching,
                                               bool                    checkForExactNumColsMatching);



  gsl_matrix*       data                      ();

  void              mpiSum                    (const MpiComm& comm, GslMatrix& M_global) const;

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





  void              print                     (std::ostream& os) const;

  void              subWriteContents          (const std::string&            varNamePrefix,
                                               const std::string&            fileName,
                                               const std::string&            fileType,
                                               const std::set<unsigned int>& allowedSubEnvIds) const;

  void              subReadContents           (const std::string&            fileName,
                                               const std::string&            fileType,
                                               const std::set<unsigned int>& allowedSubEnvIds);

private:

  GslMatrix();

  void              copy                      (const GslMatrix& src);

  void              resetLU                   ();

  int               internalSvd               () const;

          gsl_matrix*       m_mat;

  mutable gsl_matrix*       m_LU;

  mutable GslMatrix* m_inverse;

  mutable Map*       m_svdColMap;


  mutable GslMatrix* m_svdUmat;


  mutable GslVector* m_svdSvec;


  mutable GslMatrix* m_svdVmat;


  mutable GslMatrix* m_svdVTmat;

  mutable double            m_determinant;

  mutable double            m_lnDeterminant;


  mutable gsl_permutation*  m_permutation;


  mutable int               m_signum;

  mutable bool              m_isSingular;
};

GslMatrix operator*       (double a,                    const GslMatrix& mat);
GslVector operator*       (const GslMatrix& mat, const GslVector& vec);
GslMatrix operator*       (const GslMatrix& m1,  const GslMatrix& m2 );
GslMatrix operator+       (const GslMatrix& m1,  const GslMatrix& m2 );
GslMatrix operator-       (const GslMatrix& m1,  const GslMatrix& m2 );
GslMatrix matrixProduct   (const GslVector& v1,  const GslVector& v2 );

// Row \c i of the returned matrix is equal to row \c i of \c mat multiplied by element \c i of \c v.
GslMatrix leftDiagScaling (const GslVector& vec, const GslMatrix& mat);

// Column \c j of the returned matrix is equal to column \c j of \c mat multiplied by element \c j of \c v.
GslMatrix rightDiagScaling(const GslMatrix& mat, const GslVector& vec);
std::ostream&    operator<<      (std::ostream& os,            const GslMatrix& obj);

}  // End namespace QUESO

#endif // UQ_GSL_MATRIX_H