Epetra Package Browser (Single Doxygen Collection) Development
Epetra_SerialDenseSVD.h
Go to the documentation of this file.
00001 /*
00002 //@HEADER
00003 // ************************************************************************
00004 // 
00005 //               Epetra: Linear Algebra Services Package 
00006 //                 Copyright 2001 Sandia Corporation
00007 // 
00008 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00009 // the U.S. Government retains certain rights in this software.
00010 //
00011 // Redistribution and use in source and binary forms, with or without
00012 // modification, are permitted provided that the following conditions are
00013 // met:
00014 //
00015 // 1. Redistributions of source code must retain the above copyright
00016 // notice, this list of conditions and the following disclaimer.
00017 //
00018 // 2. Redistributions in binary form must reproduce the above copyright
00019 // notice, this list of conditions and the following disclaimer in the
00020 // documentation and/or other materials provided with the distribution.
00021 //
00022 // 3. Neither the name of the Corporation nor the names of the
00023 // contributors may be used to endorse or promote products derived from
00024 // this software without specific prior written permission.
00025 //
00026 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00027 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00028 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00029 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00030 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00031 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00032 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00033 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00034 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00035 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00036 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00037 //
00038 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00039 // 
00040 // ************************************************************************
00041 //@HEADER
00042 */
00043 
00044 #ifndef _EPETRA_SERIALDENSESVD_H_
00045 #define _EPETRA_SERIALDENSESVD_H_
00046 
00047 #include "Epetra_SerialDenseOperator.h"
00048 #include "Epetra_SerialDenseMatrix.h"
00049 #include "Epetra_Object.h" 
00050 #include "Epetra_CompObject.h"
00051 #include "Epetra_BLAS.h"
00052 #include "Epetra_LAPACK.h"
00053 
00054 
00056 
00114 //=========================================================================
00115 class Epetra_SerialDenseSVD : public virtual Epetra_SerialDenseOperator, public Epetra_CompObject, public virtual Epetra_Object, public Epetra_BLAS, public Epetra_LAPACK{
00116   public:
00117   
00119 
00120 
00121   Epetra_SerialDenseSVD();
00122   
00124   virtual ~Epetra_SerialDenseSVD();
00126 
00128 
00129 
00131   int SetMatrix(Epetra_SerialDenseMatrix & A);
00132 
00134 
00137   int SetVectors(Epetra_SerialDenseMatrix & X, Epetra_SerialDenseMatrix & B);
00139 
00141 
00142 
00144 
00146 //  void FactorWithEquilibration(bool Flag) {Equilibrate_ = Flag; return;};
00147 
00149   void SolveWithTranspose(bool Flag) {Transpose_ = Flag; if (Flag) TRANS_ = 'T'; else TRANS_ = 'N'; return;};
00150 
00152 //  void SolveToRefinedSolution(bool Flag) {RefineSolution_ = Flag; return;};
00153 
00154   // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
00155   // Causes all solves to estimate the forward and backward solution error. 
00156   /* Error estimates will be in the arrays FERR and BERR, resp, after the solve step is complete.
00157       These arrays are accessible via the FERR() and BERR() access functions.
00158   */
00159 //  void EstimateSolutionErrors(bool Flag) {EstimateSolutionErrors_ = Flag; return;};
00161 
00163 
00164 
00165   // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
00166   // Computes the SVD factorization of the matrix using the LAPACK routine \e DGESVD.
00167   /* 
00168     \return Integer error code, set to 0 if successful.
00169   */
00170 //  virtual int Factor(void);
00171   virtual int Factor(void);
00172 
00174 
00177   virtual int Solve(void);
00178 
00180 
00183   virtual int Invert( double rthresh = 0.0, double athresh = 0.0 );
00184 
00185   // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
00186   // Computes the scaling vector S(i) = 1/sqrt(A(i,i) of the \e this matrix.
00187   /* 
00188     \return Integer error code, set to 0 if successful. Otherwise returns the LAPACK error code INFO.
00189   */
00190 //  virtual int ComputeEquilibrateScaling(void);
00191 
00192   // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
00193   // Equilibrates the \e this matrix.
00194   /* 
00195     \return Integer error code, set to 0 if successful. Otherwise returns the LAPACK error code INFO.
00196   */
00197 //  virtual int EquilibrateMatrix(void);
00198 
00199   // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
00200   // Equilibrates the current RHS.
00201   /* 
00202     \return Integer error code, set to 0 if successful. Otherwise returns the LAPACK error code INFO.
00203   */
00204 //  int EquilibrateRHS(void);
00205 
00206 
00207   // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
00208   // Apply Iterative Refinement.
00209   /* 
00210     \return Integer error code, set to 0 if successful. Otherwise returns the LAPACK error code INFO.
00211   */
00212 //  virtual int ApplyRefinement(void);
00213 
00214   // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
00215   // Unscales the solution vectors if equilibration was used to solve the system.
00216   /* 
00217     \return Integer error code, set to 0 if successful. Otherwise returns the LAPACK error code INFO.
00218   */
00219 //  int UnequilibrateLHS(void);
00220 
00221   // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
00222   // Returns the reciprocal of the 1-norm condition number of the \e this matrix.
00223   /* 
00224     \param Value Out
00225            On return contains the reciprocal of the 1-norm condition number of the \e this matrix.
00226     
00227     \return Integer error code, set to 0 if successful. Otherwise returns the LAPACK error code INFO.
00228   */
00229 //  virtual int ReciprocalConditionEstimate(double & Value);
00231 
00233 
00234 
00236   bool Transpose() {return(Transpose_);};
00237 
00239   bool Factored() {return(Factored_);};
00240 
00241   // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
00242   // Returns true if factor is equilibrated (factor available via AF() and LDAF()).
00243 //  bool A_Equilibrated() {return(A_Equilibrated_);};
00244 
00245   // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
00246   // Returns true if RHS is equilibrated (RHS available via B() and LDB()).
00247 //  bool B_Equilibrated() {return(B_Equilibrated_);};
00248 
00249   // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
00250   // Returns true if the LAPACK general rules for equilibration suggest you should equilibrate the system.
00251 //  virtual bool ShouldEquilibrate() {ComputeEquilibrateScaling(); return(ShouldEquilibrate_);};
00252 
00253   // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
00254   // Returns true if forward and backward error estimated have been computed (available via FERR() and BERR()).
00255 //  bool SolutionErrorsEstimated() {return(SolutionErrorsEstimated_);};
00256 
00258   bool Inverted() {return(Inverted_);};
00259 
00260   // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
00261   // Returns true if the condition number of the \e this matrix has been computed (value available via ReciprocalConditionEstimate()).
00262 //  bool ReciprocalConditionEstimated() {return(ReciprocalConditionEstimated_);};
00263 
00265   bool Solved() {return(Solved_);};
00266 
00267   // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
00268   // Returns true if the current set of vectors has been refined.
00269 //  bool SolutionRefined() {return(SolutionRefined_);};
00271 
00273 
00274     
00276    Epetra_SerialDenseMatrix * Matrix()  const {return(Matrix_);};
00277        
00278   // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
00279   // Returns pointer to factored matrix (assuming factorization has been performed).
00280 //   Epetra_SerialDenseMatrix * FactoredMatrix()  const {return(Factor_);};
00281 
00283    Epetra_SerialDenseMatrix * InvertedMatrix()  const {return(Inverse_);};
00284 
00286   Epetra_SerialDenseMatrix * LHS()  const {return(LHS_);};
00287     
00289   Epetra_SerialDenseMatrix * RHS()  const {return(RHS_);};
00290     
00292   int M()  const {return(M_);};
00293 
00295   int N()  const {return(N_);};
00296 
00298   double * A()  const {return(A_);};
00299 
00301   int LDA()  const {return(LDA_);};
00302 
00304   double * B()  const {return(B_);};
00305 
00307   int LDB()  const {return(LDB_);};
00308 
00310   int NRHS()  const {return(NRHS_);};
00311 
00313   double * X()  const {return(X_);};
00314 
00316   int LDX()  const {return(LDX_);};
00317 
00318   double * S() const {return(S_);};
00319 
00320   // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
00321   // Returns pointer to the factored matrix (may be the same as A() if factorization done in place).
00322 //  double * AF()  const {return(AF_);};
00323 
00324   // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
00325   // Returns the leading dimension of the factored matrix.
00326 //  int LDAF()  const {return(LDAF_);};
00327 
00329   double * AI()  const {return(AI_);};
00330 
00332   int LDAI()  const {return(LDAI_);};
00333 
00334   // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
00335   // Returns pointer to pivot vector (if factorization has been computed), zero otherwise.
00336 //  int * IPIV()  const {return(IPIV_);};
00337 
00339   double ANORM()  const {return(ANORM_);};
00340 
00341   // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
00342   // Returns the reciprocal of the condition number of the \e this matrix (returns -1 if not yet computed).
00343 //  double RCOND()  const {return(RCOND_);};
00344 
00345   // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
00346   // Ratio of smallest to largest row scale factors for the \e this matrix (returns -1 if not yet computed).
00347   /* If ROWCND() is >= 0.1 and AMAX() is not close to overflow or underflow, then equilibration is not needed.
00348    */
00349 //  double ROWCND()  const {return(ROWCND_);};
00350 
00351   // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
00352   // Ratio of smallest to largest column scale factors for the \e this matrix (returns -1 if not yet computed).
00353   /* If COLCND() is >= 0.1 then equilibration is not needed.
00354    */
00355 //  double COLCND()  const {return(COLCND_);};
00356 
00357   // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
00358   // Returns the absolute value of the largest entry of the \e this matrix (returns -1 if not yet computed).
00359 //  double AMAX()  const {return(AMAX_);};
00360 
00361   // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
00362   // Returns a pointer to the forward error estimates computed by LAPACK.
00363 //  double * FERR()  const {return(FERR_);};
00364 
00365   // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
00366   // Returns a pointer to the backward error estimates computed by LAPACK.
00367 //  double * BERR()  const {return(BERR_);};
00368 
00369   // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
00370   // Returns a pointer to the row scaling vector used for equilibration.
00371 //  double * R()  const {return(R_);};
00372 
00373   // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
00374   // Returns a pointer to the column scale vector used for equilibration.
00375 //  double * C()  const {return(C_);};
00377 
00379 
00380 
00381   virtual void Print(ostream& os) const;
00383 
00385 
00386 
00388 
00397     virtual int SetUseTranspose(bool use_transpose) { UseTranspose_ = use_transpose; return (0); }
00398 
00400 
00408     virtual int Apply(const Epetra_SerialDenseMatrix& Xmat, Epetra_SerialDenseMatrix& Ymat)
00409     { return Ymat.Multiply( UseTranspose_, false, 1.0, *Matrix(), Xmat, 0.0 ); }
00410 
00412 
00421     virtual int ApplyInverse(const Epetra_SerialDenseMatrix & Xmat, Epetra_SerialDenseMatrix & Ymat)
00422     { SetVectors(const_cast<Epetra_SerialDenseMatrix&>(Xmat),Ymat);
00423       SolveWithTranspose(UseTranspose_);
00424       return Solve(); }
00425 
00427     /* Returns the quantity \f$ \| A \|_\infty\f$ such that
00428        \f[\| A \|_\infty = \max_{1\lei\lem} \sum_{j=1}^n |a_{ij}| \f].
00429 
00430        \warning This method must not be called unless HasNormInf() returns true.
00431     */ 
00432     virtual double NormInf() const { return Matrix()->NormInf(); }
00433   
00435     virtual const char * Label() const { return Epetra_Object::Label(); }
00436 
00438     virtual bool UseTranspose() const { return UseTranspose_; }
00439 
00441     virtual bool HasNormInf() const { return true; }
00442 
00444     virtual int RowDim() const { return M(); }
00445 
00447     virtual int ColDim() const { return N(); }
00448 
00450   
00451   void AllocateWORK() {if (WORK_==0) {LWORK_ = 4*N_; WORK_ = new double[LWORK_];} return;};
00452   void AllocateIWORK() {if (IWORK_==0) IWORK_ = new int[N_]; return;};
00453   void InitPointers();
00454   void DeleteArrays();
00455   void ResetMatrix();
00456   void ResetVectors();
00457 
00458 
00459 //  bool Equilibrate_;
00460 //  bool ShouldEquilibrate_;
00461 //  bool A_Equilibrated_;
00462 //  bool B_Equilibrated_;
00463   bool Transpose_;
00464   bool Factored_;
00465 //  bool EstimateSolutionErrors_;
00466 //  bool SolutionErrorsEstimated_;
00467   bool Solved_;
00468   bool Inverted_;
00469 //  bool ReciprocalConditionEstimated_;
00470 //  bool RefineSolution_;
00471 //  bool SolutionRefined_;
00472 
00473   char TRANS_;
00474 
00475   int M_;
00476   int N_;
00477   int Min_MN_;
00478   int NRHS_;
00479   int LDA_;
00480 //  int LDAF_;
00481   int LDAI_;
00482   int LDB_;
00483   int LDX_;
00484   int INFO_;
00485   int LWORK_;
00486 
00487 //  int * IPIV_;
00488   int * IWORK_;
00489 
00490   double ANORM_;
00491 //  double RCOND_;
00492 //  double ROWCND_;
00493 //  double COLCND_;
00494 //  double AMAX_;
00495 
00496   Epetra_SerialDenseMatrix * Matrix_;
00497   Epetra_SerialDenseMatrix * LHS_;
00498   Epetra_SerialDenseMatrix * RHS_;
00499 //  Epetra_SerialDenseMatrix * Factor_;
00500   Epetra_SerialDenseMatrix * Inverse_;
00501   
00502   double * A_;
00503 //  double * FERR_;
00504 //  double * BERR_;
00505 //  double * AF_;
00506   double * AI_;
00507   double * WORK_;
00508 //  double * R_;
00509 //  double * C_;
00510   double * U_;
00511   double * S_;
00512   double * Vt_;
00513 
00514   double * B_;
00515   double * X_;
00516 
00517   bool UseTranspose_;
00518 
00519  private:
00520   // Epetra_SerialDenseSolver copy constructor (put here because we don't want user access)
00521   
00522   Epetra_SerialDenseSVD(const Epetra_SerialDenseSVD& Source);
00523   Epetra_SerialDenseSVD & operator=(const Epetra_SerialDenseSVD& Source);
00524 };
00525 
00526 #endif /* _EPETRA_SERIALDENSESVD_H_ */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines