Epetra_SerialDenseSVD.h

Go to the documentation of this file.

/*
//@HEADER
// ************************************************************************
//
//               Epetra: Linear Algebra Services Package
//                 Copyright 2011 Sandia Corporation
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// the U.S. Government retains certain rights in this software.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// 3. Neither the name of the Corporation nor the names of the
// contributors may be used to endorse or promote products derived from
// this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
//
// ************************************************************************
//@HEADER
*/
 
#ifndef _EPETRA_SERIALDENSESVD_H_
#define _EPETRA_SERIALDENSESVD_H_
 
#include "Epetra_SerialDenseOperator.h"
#include "Epetra_SerialDenseMatrix.h"
#include "Epetra_Object.h"
#include "Epetra_CompObject.h"
#include "Epetra_BLAS.h"
#include "Epetra_LAPACK.h"
 
 
 
//=========================================================================
class Epetra_SerialDenseSVD : public virtual Epetra_SerialDenseOperator, public Epetra_CompObject, public virtual Epetra_Object, public Epetra_BLAS, public Epetra_LAPACK{
  public:
 
 
  Epetra_SerialDenseSVD();
 
  virtual ~Epetra_SerialDenseSVD();
 
 
 
  int SetMatrix(Epetra_SerialDenseMatrix & A);
 
 
  int SetVectors(Epetra_SerialDenseMatrix & X, Epetra_SerialDenseMatrix & B);
 
 
 
 
//  void FactorWithEquilibration(bool Flag) {Equilibrate_ = Flag; return;};
 
  void SolveWithTranspose(bool Flag) {Transpose_ = Flag; if (Flag) TRANS_ = 'T'; else TRANS_ = 'N'; return;};
 
//  void SolveToRefinedSolution(bool Flag) {RefineSolution_ = Flag; return;};
 
  // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
  // Causes all solves to estimate the forward and backward solution error.
  /* Error estimates will be in the arrays FERR and BERR, resp, after the solve step is complete.
      These arrays are accessible via the FERR() and BERR() access functions.
  */
//  void EstimateSolutionErrors(bool Flag) {EstimateSolutionErrors_ = Flag; return;};
 
 
 
  // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
  // Computes the SVD factorization of the matrix using the LAPACK routine \e DGESVD.
  /*
    \return Integer error code, set to 0 if successful.
  */
//  virtual int Factor(void);
  virtual int Factor(void);
 
 
  virtual int Solve(void);
 
 
  virtual int Invert( double rthresh = 0.0, double athresh = 0.0 );
 
  // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
  // Computes the scaling vector S(i) = 1/sqrt(A(i,i) of the \e this matrix.
  /*
    \return Integer error code, set to 0 if successful. Otherwise returns the LAPACK error code INFO.
  */
//  virtual int ComputeEquilibrateScaling(void);
 
  // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
  // Equilibrates the \e this matrix.
  /*
    \return Integer error code, set to 0 if successful. Otherwise returns the LAPACK error code INFO.
  */
//  virtual int EquilibrateMatrix(void);
 
  // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
  // Equilibrates the current RHS.
  /*
    \return Integer error code, set to 0 if successful. Otherwise returns the LAPACK error code INFO.
  */
//  int EquilibrateRHS(void);
 
 
  // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
  // Apply Iterative Refinement.
  /*
    \return Integer error code, set to 0 if successful. Otherwise returns the LAPACK error code INFO.
  */
//  virtual int ApplyRefinement(void);
 
  // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
  // Unscales the solution vectors if equilibration was used to solve the system.
  /*
    \return Integer error code, set to 0 if successful. Otherwise returns the LAPACK error code INFO.
  */
//  int UnequilibrateLHS(void);
 
  // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
  // Returns the reciprocal of the 1-norm condition number of the \e this matrix.
  /*
    \param Value Out
           On return contains the reciprocal of the 1-norm condition number of the \e this matrix.
 
    \return Integer error code, set to 0 if successful. Otherwise returns the LAPACK error code INFO.
  */
//  virtual int ReciprocalConditionEstimate(double & Value);
 
 
 
  bool Transpose() {return(Transpose_);};
 
  bool Factored() {return(Factored_);};
 
  // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
  // Returns true if factor is equilibrated (factor available via AF() and LDAF()).
//  bool A_Equilibrated() {return(A_Equilibrated_);};
 
  // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
  // Returns true if RHS is equilibrated (RHS available via B() and LDB()).
//  bool B_Equilibrated() {return(B_Equilibrated_);};
 
  // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
  // Returns true if the LAPACK general rules for equilibration suggest you should equilibrate the system.
//  virtual bool ShouldEquilibrate() {ComputeEquilibrateScaling(); return(ShouldEquilibrate_);};
 
  // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
  // Returns true if forward and backward error estimated have been computed (available via FERR() and BERR()).
//  bool SolutionErrorsEstimated() {return(SolutionErrorsEstimated_);};
 
  bool Inverted() {return(Inverted_);};
 
  // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
  // Returns true if the condition number of the \e this matrix has been computed (value available via ReciprocalConditionEstimate()).
//  bool ReciprocalConditionEstimated() {return(ReciprocalConditionEstimated_);};
 
  bool Solved() {return(Solved_);};
 
  // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
  // Returns true if the current set of vectors has been refined.
//  bool SolutionRefined() {return(SolutionRefined_);};
 
 
 
   Epetra_SerialDenseMatrix * Matrix()  const {return(Matrix_);};
 
  // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
  // Returns pointer to factored matrix (assuming factorization has been performed).
//   Epetra_SerialDenseMatrix * FactoredMatrix()  const {return(Factor_);};
 
   Epetra_SerialDenseMatrix * InvertedMatrix()  const {return(Inverse_);};
 
  Epetra_SerialDenseMatrix * LHS()  const {return(LHS_);};
 
  Epetra_SerialDenseMatrix * RHS()  const {return(RHS_);};
 
  int M()  const {return(M_);};
 
  int N()  const {return(N_);};
 
  double * A()  const {return(A_);};
 
  int LDA()  const {return(LDA_);};
 
  double * B()  const {return(B_);};
 
  int LDB()  const {return(LDB_);};
 
  int NRHS()  const {return(NRHS_);};
 
  double * X()  const {return(X_);};
 
  int LDX()  const {return(LDX_);};
 
  double * S() const {return(S_);};
 
  // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
  // Returns pointer to the factored matrix (may be the same as A() if factorization done in place).
//  double * AF()  const {return(AF_);};
 
  // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
  // Returns the leading dimension of the factored matrix.
//  int LDAF()  const {return(LDAF_);};
 
  double * AI()  const {return(AI_);};
 
  int LDAI()  const {return(LDAI_);};
 
  // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
  // Returns pointer to pivot vector (if factorization has been computed), zero otherwise.
//  int * IPIV()  const {return(IPIV_);};
 
  double ANORM()  const {return(ANORM_);};
 
  // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
  // Returns the reciprocal of the condition number of the \e this matrix (returns -1 if not yet computed).
//  double RCOND()  const {return(RCOND_);};
 
  // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
  // Ratio of smallest to largest row scale factors for the \e this matrix (returns -1 if not yet computed).
  /* If ROWCND() is >= 0.1 and AMAX() is not close to overflow or underflow, then equilibration is not needed.
   */
//  double ROWCND()  const {return(ROWCND_);};
 
  // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
  // Ratio of smallest to largest column scale factors for the \e this matrix (returns -1 if not yet computed).
  /* If COLCND() is >= 0.1 then equilibration is not needed.
   */
//  double COLCND()  const {return(COLCND_);};
 
  // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
  // Returns the absolute value of the largest entry of the \e this matrix (returns -1 if not yet computed).
//  double AMAX()  const {return(AMAX_);};
 
  // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
  // Returns a pointer to the forward error estimates computed by LAPACK.
//  double * FERR()  const {return(FERR_);};
 
  // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
  // Returns a pointer to the backward error estimates computed by LAPACK.
//  double * BERR()  const {return(BERR_);};
 
  // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
  // Returns a pointer to the row scaling vector used for equilibration.
//  double * R()  const {return(R_);};
 
  // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
  // Returns a pointer to the column scale vector used for equilibration.
//  double * C()  const {return(C_);};
 
 
  virtual void Print(std::ostream& os) const;
 
 
 
 
    virtual int SetUseTranspose(bool use_transpose) { UseTranspose_ = use_transpose; return (0); }
 
 
    virtual int Apply(const Epetra_SerialDenseMatrix& Xmat, Epetra_SerialDenseMatrix& Ymat)
    { return Ymat.Multiply( UseTranspose_, false, 1.0, *Matrix(), Xmat, 0.0 ); }
 
 
    virtual int ApplyInverse(const Epetra_SerialDenseMatrix & Xmat, Epetra_SerialDenseMatrix & Ymat)
    { SetVectors(const_cast<Epetra_SerialDenseMatrix&>(Xmat),Ymat);
      SolveWithTranspose(UseTranspose_);
      return Solve(); }
 
    /* Returns the quantity \f$ \| A \|_\infty\f$ such that
       \f[\| A \|_\infty = \max_{1\lei\lem} \sum_{j=1}^n |a_{ij}| \f].
 
       \warning This method must not be called unless HasNormInf() returns true.
    */
    virtual double NormInf() const { return Matrix()->NormInf(); }
 
    virtual const char * Label() const { return Epetra_Object::Label(); }
 
    virtual bool UseTranspose() const { return UseTranspose_; }
 
    virtual bool HasNormInf() const { return true; }
 
    virtual int RowDim() const { return M(); }
 
    virtual int ColDim() const { return N(); }
 
 
  void AllocateWORK() {if (WORK_==0) {LWORK_ = 4*N_; WORK_ = new double[LWORK_];} return;};
  void AllocateIWORK() {if (IWORK_==0) IWORK_ = new int[N_]; return;};
  void InitPointers();
  void DeleteArrays();
  void ResetMatrix();
  void ResetVectors();
 
 
//  bool Equilibrate_;
//  bool ShouldEquilibrate_;
//  bool A_Equilibrated_;
//  bool B_Equilibrated_;
  bool Transpose_;
  bool Factored_;
//  bool EstimateSolutionErrors_;
//  bool SolutionErrorsEstimated_;
  bool Solved_;
  bool Inverted_;
//  bool ReciprocalConditionEstimated_;
//  bool RefineSolution_;
//  bool SolutionRefined_;
 
  char TRANS_;
 
  int M_;
  int N_;
  int Min_MN_;
  int NRHS_;
  int LDA_;
//  int LDAF_;
  int LDAI_;
  int LDB_;
  int LDX_;
  int INFO_;
  int LWORK_;
 
//  int * IPIV_;
  int * IWORK_;
 
  double ANORM_;
//  double RCOND_;
//  double ROWCND_;
//  double COLCND_;
//  double AMAX_;
 
  Epetra_SerialDenseMatrix * Matrix_;
  Epetra_SerialDenseMatrix * LHS_;
  Epetra_SerialDenseMatrix * RHS_;
//  Epetra_SerialDenseMatrix * Factor_;
  Epetra_SerialDenseMatrix * Inverse_;
 
  double * A_;
//  double * FERR_;
//  double * BERR_;
//  double * AF_;
  double * AI_;
  double * WORK_;
//  double * R_;
//  double * C_;
  double * U_;
  double * S_;
  double * Vt_;
 
  double * B_;
  double * X_;
 
  bool UseTranspose_;
 
 private:
  // Epetra_SerialDenseSolver copy constructor (put here because we don't want user access)
 
  Epetra_SerialDenseSVD(const Epetra_SerialDenseSVD& Source);
  Epetra_SerialDenseSVD & operator=(const Epetra_SerialDenseSVD& Source);
};
 
#endif /* _EPETRA_SERIALDENSESVD_H_ */