TeuchosMatrix.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_TEUCHOS_MATRIX_H
#define UQ_TEUCHOS_MATRIX_H

#ifdef QUESO_HAS_TRILINOS
#include <queso/TeuchosVector.h>
#include <Teuchos_SerialDenseVector.hpp>
#include <Teuchos_SerialDenseMatrix.hpp>
#include <Teuchos_LAPACK.hpp>
#endif

#include <queso/Matrix.h>

#ifdef QUESO_HAS_TRILINOS

namespace QUESO {

class TeuchosMatrix : public Matrix
{
public:


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

  TeuchosMatrix(const BaseEnvironment& env,
                       const Map&             map,
                       double                        diagValue);

  TeuchosMatrix(const TeuchosVector&   v,
                           double                    diagValue);


  TeuchosMatrix(const TeuchosVector&   v);

  TeuchosMatrix(const TeuchosMatrix&   B);

 ~TeuchosMatrix();


  TeuchosMatrix& operator= (const TeuchosMatrix& rhs);

  TeuchosMatrix& operator*=(double a);

  TeuchosMatrix& operator/=(double a);

  TeuchosMatrix& operator+=(const TeuchosMatrix& rhs);
  TeuchosMatrix& operator-=(const TeuchosMatrix& 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*           values                 () ;//added by Kemelli on 12/04/12

  int               stride                 () ;//added by Kemelli on 12/04/12

  double            max                       () const;


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

  TeuchosMatrix  transpose                 () const;

  TeuchosMatrix  inverse                   () const;

  double                determinant               () const;

  double                lnDeterminant             () const;


  double            normFrob                  () const;

  double            normMax                   () const;



  int               chol                   () ;

  int               svd                    (TeuchosMatrix& matU, TeuchosVector& vecS, TeuchosMatrix& matVt) const;


  const TeuchosMatrix& svdMatU      () const;


  const TeuchosMatrix& svdMatV      () const;


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


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



  TeuchosVector  multiply                  (const TeuchosVector& x) const;


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


  TeuchosVector  invertMultiply            (const TeuchosVector& b) const;


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


  TeuchosMatrix  invertMultiply            (const TeuchosMatrix& B) const;


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


  TeuchosVector  invertMultiplyForceLU     (const TeuchosVector& b) const;

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

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

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




  void                  cwSet                     (double value);

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


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

  TeuchosVector  getColumn                 (const unsigned int column_num) const;

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

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

  TeuchosVector  getRow                    (const unsigned int row_num) const;

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


  void                  zeroLower              (bool includeDiagonal = false);


  void                  zeroUpper              (bool includeDiagonal = false);


  void               filterSmallValues         (double thresholdValue);


  void               filterLargeValues         (double thresholdValue);

  void               fillWithTranspose         (const TeuchosMatrix& mat);

  void               fillWithBlocksDiagonally  (const std::vector<const TeuchosMatrix* >& matrices);

  void               fillWithBlocksDiagonally  (const std::vector<      TeuchosMatrix* >& matrices);

  void               fillWithBlocksHorizontally(const std::vector<const TeuchosMatrix* >& matrices);

  void               fillWithBlocksHorizontally(const std::vector<      TeuchosMatrix* >& matrices);

  void               fillWithBlocksVertically  (const std::vector<const TeuchosMatrix* >& matrices);

  void               fillWithBlocksVertically  (const std::vector<      TeuchosMatrix* >& matrices);

  void               fillWithTensorProduct     (const TeuchosMatrix& mat1, const TeuchosMatrix& mat2);

  void               fillWithTensorProduct     (const TeuchosMatrix& mat1, const TeuchosVector& vec2);




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

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


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

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

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

private:

  TeuchosMatrix();

  void              copy                      (const TeuchosMatrix& src);

  void              resetLU                   ();

  void              multiply                  (const TeuchosVector& x, TeuchosVector& y) const;


  int               internalSvd               () const;

  Teuchos::SerialDenseMatrix<int,double> m_mat;

  mutable Teuchos::SerialDenseMatrix<int,double> m_LU;

  mutable TeuchosMatrix* m_inverse;

  mutable Map*          m_svdColMap;

  mutable TeuchosMatrix* m_svdUmat;

  mutable TeuchosVector* m_svdSvec;

  mutable TeuchosMatrix* m_svdVmat;

  mutable TeuchosMatrix* m_svdVTmat;

  mutable double                m_determinant;

  mutable double                m_lnDeterminant;

  mutable TeuchosMatrix*  m_permutation;

  mutable int*              v_pivoting;


  mutable int               m_signum;

  mutable bool              m_isSingular;
};

TeuchosMatrix operator*       (double a,                    const TeuchosMatrix& mat);
TeuchosVector operator*       (const TeuchosMatrix& mat, const TeuchosVector& vec);
TeuchosMatrix operator*       (const TeuchosMatrix& m1,  const TeuchosMatrix& m2 );
TeuchosMatrix operator+       (const TeuchosMatrix& m1,  const TeuchosMatrix& m2 );
TeuchosMatrix operator-       (const TeuchosMatrix& m1,  const TeuchosMatrix& m2 );
TeuchosMatrix matrixProduct   (const TeuchosVector& v1,  const TeuchosVector& v2 );
TeuchosMatrix leftDiagScaling (const TeuchosVector& vec, const TeuchosMatrix& mat);
TeuchosMatrix rightDiagScaling(const TeuchosMatrix& mat, const TeuchosVector& vec);
std::ostream&        operator<<      (std::ostream& os,            const TeuchosMatrix& obj);

}  // End namespace QUESO

#endif // ifdef QUESO_HAS_TRILINOS

#endif // UQ_TEUCHOS_MATRIX_H