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,2009,2010,2011,2012,2013 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/Matrix.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_permutation.h>
#include <queso/GslVector.h>

namespace QUESO {

class GslMatrix : public Matrix
{
public:



  GslMatrix();
  
  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);
            
  const double& operator()(unsigned int i, unsigned int j) const;
  
 


  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;


  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:
  void              copy                      (const GslMatrix& src);
  
  void              resetLU                   ();
        
  void              multiply                  (const GslVector& x, GslVector& y) const;
        
  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