Epetra_SerialDenseSVD.h

Go to the documentation of this file.
00001 
00002 /* Copyright (2001) Sandia Corportation. Under the terms of Contract 
00003  * DE-AC04-94AL85000, there is a non-exclusive license for use of this 
00004  * work by or on behalf of the U.S. Government.  Export of this program
00005  * may require a license from the United States Government. */
00006 
00007 
00008 /* NOTICE:  The United States Government is granted for itself and others
00009  * acting on its behalf a paid-up, nonexclusive, irrevocable worldwide
00010  * license in ths data to reproduce, prepare derivative works, and
00011  * perform publicly and display publicly.  Beginning five (5) years from
00012  * July 25, 2001, the United States Government is granted for itself and
00013  * others acting on its behalf a paid-up, nonexclusive, irrevocable
00014  * worldwide license in this data to reproduce, prepare derivative works,
00015  * distribute copies to the public, perform publicly and display
00016  * publicly, and to permit others to do so.
00017  * 
00018  * NEITHER THE UNITED STATES GOVERNMENT, NOR THE UNITED STATES DEPARTMENT
00019  * OF ENERGY, NOR SANDIA CORPORATION, NOR ANY OF THEIR EMPLOYEES, MAKES
00020  * ANY WARRANTY, EXPRESS OR IMPLIED, OR ASSUMES ANY LEGAL LIABILITY OR
00021  * RESPONSIBILITY FOR THE ACCURACY, COMPLETENESS, OR USEFULNESS OF ANY
00022  * INFORMATION, APPARATUS, PRODUCT, OR PROCESS DISCLOSED, OR REPRESENTS
00023  * THAT ITS USE WOULD NOT INFRINGE PRIVATELY OWNED RIGHTS. */
00024 
00025 #ifndef _EPETRA_SERIALDENSESVD_H_
00026 #define _EPETRA_SERIALDENSESVD_H_
00027 
00028 #include "Epetra_SerialDenseOperator.h"
00029 #include "Epetra_SerialDenseMatrix.h"
00030 #include "Epetra_Object.h" 
00031 #include "Epetra_CompObject.h"
00032 #include "Epetra_BLAS.h"
00033 #include "Epetra_LAPACK.h"
00034 
00035 
00037 
00095 //=========================================================================
00096 class Epetra_SerialDenseSVD : public virtual Epetra_SerialDenseOperator, public Epetra_CompObject, public virtual Epetra_Object, public Epetra_BLAS, public Epetra_LAPACK{
00097   public:
00098   
00100 
00101 
00102   Epetra_SerialDenseSVD();
00103   
00105   virtual ~Epetra_SerialDenseSVD();
00107 
00109 
00110 
00112   int SetMatrix(Epetra_SerialDenseMatrix & A);
00113 
00115 
00118   int SetVectors(Epetra_SerialDenseMatrix & X, Epetra_SerialDenseMatrix & B);
00120 
00122 
00123 
00125 
00127 //  void FactorWithEquilibration(bool Flag) {Equilibrate_ = Flag; return;};
00128 
00130   void SolveWithTranspose(bool Flag) {Transpose_ = Flag; if (Flag) TRANS_ = 'T'; else TRANS_ = 'N'; return;};
00131 
00133 //  void SolveToRefinedSolution(bool Flag) {RefineSolution_ = Flag; return;};
00134 
00135   // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
00136   // Causes all solves to estimate the forward and backward solution error. 
00137   /* Error estimates will be in the arrays FERR and BERR, resp, after the solve step is complete.
00138       These arrays are accessible via the FERR() and BERR() access functions.
00139   */
00140 //  void EstimateSolutionErrors(bool Flag) {EstimateSolutionErrors_ = Flag; return;};
00142 
00144 
00145 
00146   // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
00147   // Computes the SVD factorization of the matrix using the LAPACK routine \e DGESVD.
00148   /* 
00149     \return Integer error code, set to 0 if successful.
00150   */
00151 //  virtual int Factor(void);
00152   virtual int Factor(void);
00153 
00155 
00158   virtual int Solve(void);
00159 
00161 
00164   virtual int Invert( double rthresh = 0.0, double athresh = 0.0 );
00165 
00166   // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
00167   // Computes the scaling vector S(i) = 1/sqrt(A(i,i) of the \e this matrix.
00168   /* 
00169     \return Integer error code, set to 0 if successful. Otherwise returns the LAPACK error code INFO.
00170   */
00171 //  virtual int ComputeEquilibrateScaling(void);
00172 
00173   // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
00174   // Equilibrates the \e this matrix.
00175   /* 
00176     \return Integer error code, set to 0 if successful. Otherwise returns the LAPACK error code INFO.
00177   */
00178 //  virtual int EquilibrateMatrix(void);
00179 
00180   // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
00181   // Equilibrates the current RHS.
00182   /* 
00183     \return Integer error code, set to 0 if successful. Otherwise returns the LAPACK error code INFO.
00184   */
00185 //  int EquilibrateRHS(void);
00186 
00187 
00188   // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
00189   // Apply Iterative Refinement.
00190   /* 
00191     \return Integer error code, set to 0 if successful. Otherwise returns the LAPACK error code INFO.
00192   */
00193 //  virtual int ApplyRefinement(void);
00194 
00195   // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
00196   // Unscales the solution vectors if equilibration was used to solve the system.
00197   /* 
00198     \return Integer error code, set to 0 if successful. Otherwise returns the LAPACK error code INFO.
00199   */
00200 //  int UnequilibrateLHS(void);
00201 
00202   // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
00203   // Returns the reciprocal of the 1-norm condition number of the \e this matrix.
00204   /* 
00205     \param Value Out
00206            On return contains the reciprocal of the 1-norm condition number of the \e this matrix.
00207     
00208     \return Integer error code, set to 0 if successful. Otherwise returns the LAPACK error code INFO.
00209   */
00210 //  virtual int ReciprocalConditionEstimate(double & Value);
00212 
00214 
00215 
00217   bool Transpose() {return(Transpose_);};
00218 
00220   bool Factored() {return(Factored_);};
00221 
00222   // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
00223   // Returns true if factor is equilibrated (factor available via AF() and LDAF()).
00224 //  bool A_Equilibrated() {return(A_Equilibrated_);};
00225 
00226   // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
00227   // Returns true if RHS is equilibrated (RHS available via B() and LDB()).
00228 //  bool B_Equilibrated() {return(B_Equilibrated_);};
00229 
00230   // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
00231   // Returns true if the LAPACK general rules for equilibration suggest you should equilibrate the system.
00232 //  virtual bool ShouldEquilibrate() {ComputeEquilibrateScaling(); return(ShouldEquilibrate_);};
00233 
00234   // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
00235   // Returns true if forward and backward error estimated have been computed (available via FERR() and BERR()).
00236 //  bool SolutionErrorsEstimated() {return(SolutionErrorsEstimated_);};
00237 
00239   bool Inverted() {return(Inverted_);};
00240 
00241   // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
00242   // Returns true if the condition number of the \e this matrix has been computed (value available via ReciprocalConditionEstimate()).
00243 //  bool ReciprocalConditionEstimated() {return(ReciprocalConditionEstimated_);};
00244 
00246   bool Solved() {return(Solved_);};
00247 
00248   // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
00249   // Returns true if the current set of vectors has been refined.
00250 //  bool SolutionRefined() {return(SolutionRefined_);};
00252 
00254 
00255     
00257    Epetra_SerialDenseMatrix * Matrix()  const {return(Matrix_);};
00258        
00259   // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
00260   // Returns pointer to factored matrix (assuming factorization has been performed).
00261 //   Epetra_SerialDenseMatrix * FactoredMatrix()  const {return(Factor_);};
00262 
00264    Epetra_SerialDenseMatrix * InvertedMatrix()  const {return(Inverse_);};
00265 
00267   Epetra_SerialDenseMatrix * LHS()  const {return(LHS_);};
00268     
00270   Epetra_SerialDenseMatrix * RHS()  const {return(RHS_);};
00271     
00273   int M()  const {return(M_);};
00274 
00276   int N()  const {return(N_);};
00277 
00279   double * A()  const {return(A_);};
00280 
00282   int LDA()  const {return(LDA_);};
00283 
00285   double * B()  const {return(B_);};
00286 
00288   int LDB()  const {return(LDB_);};
00289 
00291   int NRHS()  const {return(NRHS_);};
00292 
00294   double * X()  const {return(X_);};
00295 
00297   int LDX()  const {return(LDX_);};
00298 
00299   double * S() const {return(S_);};
00300 
00301   // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
00302   // Returns pointer to the factored matrix (may be the same as A() if factorization done in place).
00303 //  double * AF()  const {return(AF_);};
00304 
00305   // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
00306   // Returns the leading dimension of the factored matrix.
00307 //  int LDAF()  const {return(LDAF_);};
00308 
00310   double * AI()  const {return(AI_);};
00311 
00313   int LDAI()  const {return(LDAI_);};
00314 
00315   // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
00316   // Returns pointer to pivot vector (if factorization has been computed), zero otherwise.
00317 //  int * IPIV()  const {return(IPIV_);};
00318 
00320   double ANORM()  const {return(ANORM_);};
00321 
00322   // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
00323   // Returns the reciprocal of the condition number of the \e this matrix (returns -1 if not yet computed).
00324 //  double RCOND()  const {return(RCOND_);};
00325 
00326   // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
00327   // Ratio of smallest to largest row scale factors for the \e this matrix (returns -1 if not yet computed).
00328   /* If ROWCND() is >= 0.1 and AMAX() is not close to overflow or underflow, then equilibration is not needed.
00329    */
00330 //  double ROWCND()  const {return(ROWCND_);};
00331 
00332   // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
00333   // Ratio of smallest to largest column scale factors for the \e this matrix (returns -1 if not yet computed).
00334   /* If COLCND() is >= 0.1 then equilibration is not needed.
00335    */
00336 //  double COLCND()  const {return(COLCND_);};
00337 
00338   // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
00339   // Returns the absolute value of the largest entry of the \e this matrix (returns -1 if not yet computed).
00340 //  double AMAX()  const {return(AMAX_);};
00341 
00342   // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
00343   // Returns a pointer to the forward error estimates computed by LAPACK.
00344 //  double * FERR()  const {return(FERR_);};
00345 
00346   // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
00347   // Returns a pointer to the backward error estimates computed by LAPACK.
00348 //  double * BERR()  const {return(BERR_);};
00349 
00350   // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
00351   // Returns a pointer to the row scaling vector used for equilibration.
00352 //  double * R()  const {return(R_);};
00353 
00354   // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
00355   // Returns a pointer to the column scale vector used for equilibration.
00356 //  double * C()  const {return(C_);};
00358 
00360 
00361 
00362   virtual void Print(ostream& os) const;
00364 
00366 
00367 
00369 
00378     virtual int SetUseTranspose(bool UseTranspose) { UseTranspose_ = UseTranspose; return (0); }
00379 
00381 
00389     virtual int Apply(const Epetra_SerialDenseMatrix& X, Epetra_SerialDenseMatrix& Y)
00390     { return Y.Multiply( UseTranspose_, false, 1.0, *Matrix(), X, 0.0 ); }
00391 
00393 
00402     virtual int ApplyInverse(const Epetra_SerialDenseMatrix & X, Epetra_SerialDenseMatrix & Y)
00403     { SetVectors(const_cast<Epetra_SerialDenseMatrix&>(X),Y);
00404       SolveWithTranspose(UseTranspose_);
00405       return Solve(); }
00406 
00408     /* Returns the quantity \f$ \| A \|_\infty\f$ such that
00409        \f[\| A \|_\infty = \max_{1\lei\lem} \sum_{j=1}^n |a_{ij}| \f].
00410 
00411        \warning This method must not be called unless HasNormInf() returns true.
00412     */ 
00413     virtual double NormInf() const { return Matrix()->NormInf(); }
00414   
00416     virtual const char * Label() const { return Epetra_Object::Label(); }
00417 
00419     virtual bool UseTranspose() const { return UseTranspose_; }
00420 
00422     virtual bool HasNormInf() const { return true; }
00423 
00425     virtual int RowDim() const { return M(); }
00426 
00428     virtual int ColDim() const { return N(); }
00429 
00431   
00432   void AllocateWORK() {if (WORK_==0) {LWORK_ = 4*N_; WORK_ = new double[LWORK_];} return;};
00433   void AllocateIWORK() {if (IWORK_==0) IWORK_ = new int[N_]; return;};
00434   void InitPointers();
00435   void DeleteArrays();
00436   void ResetMatrix();
00437   void ResetVectors();
00438 
00439 
00440 //  bool Equilibrate_;
00441 //  bool ShouldEquilibrate_;
00442 //  bool A_Equilibrated_;
00443 //  bool B_Equilibrated_;
00444   bool Transpose_;
00445   bool Factored_;
00446 //  bool EstimateSolutionErrors_;
00447 //  bool SolutionErrorsEstimated_;
00448   bool Solved_;
00449   bool Inverted_;
00450 //  bool ReciprocalConditionEstimated_;
00451 //  bool RefineSolution_;
00452 //  bool SolutionRefined_;
00453 
00454   char TRANS_;
00455 
00456   int M_;
00457   int N_;
00458   int Min_MN_;
00459   int NRHS_;
00460   int LDA_;
00461 //  int LDAF_;
00462   int LDAI_;
00463   int LDB_;
00464   int LDX_;
00465   int INFO_;
00466   int LWORK_;
00467 
00468 //  int * IPIV_;
00469   int * IWORK_;
00470 
00471   double ANORM_;
00472 //  double RCOND_;
00473 //  double ROWCND_;
00474 //  double COLCND_;
00475 //  double AMAX_;
00476 
00477   Epetra_SerialDenseMatrix * Matrix_;
00478   Epetra_SerialDenseMatrix * LHS_;
00479   Epetra_SerialDenseMatrix * RHS_;
00480 //  Epetra_SerialDenseMatrix * Factor_;
00481   Epetra_SerialDenseMatrix * Inverse_;
00482   
00483   double * A_;
00484 //  double * FERR_;
00485 //  double * BERR_;
00486 //  double * AF_;
00487   double * AI_;
00488   double * WORK_;
00489 //  double * R_;
00490 //  double * C_;
00491   double * U_;
00492   double * S_;
00493   double * Vt_;
00494 
00495   double * B_;
00496   double * X_;
00497 
00498   bool UseTranspose_;
00499 
00500  private:
00501   // Epetra_SerialDenseSolver copy constructor (put here because we don't want user access)
00502   
00503   Epetra_SerialDenseSVD(const Epetra_SerialDenseSVD& Source);
00504   Epetra_SerialDenseSVD & operator=(const Epetra_SerialDenseSVD& Source);
00505 };
00506 
00507 #endif /* _EPETRA_SERIALDENSESVD_H_ */

Generated on Thu Sep 18 12:37:58 2008 for Epetra Package Browser (Single Doxygen Collection) by doxygen 1.3.9.1