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_ */