Epetra Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Epetra_SerialDenseSVD.h
Go to the documentation of this file.
1 /*
2 //@HEADER
3 // ************************************************************************
4 //
5 // Epetra: Linear Algebra Services Package
6 // Copyright 2011 Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
39 //
40 // ************************************************************************
41 //@HEADER
42 */
43 
44 #ifndef _EPETRA_SERIALDENSESVD_H_
45 #define _EPETRA_SERIALDENSESVD_H_
46 
49 #include "Epetra_Object.h"
50 #include "Epetra_CompObject.h"
51 #include "Epetra_BLAS.h"
52 #include "Epetra_LAPACK.h"
53 
54 
56 
114 //=========================================================================
115 class EPETRA_LIB_DLL_EXPORT Epetra_SerialDenseSVD : public virtual Epetra_SerialDenseOperator, public Epetra_CompObject, public virtual Epetra_Object, public Epetra_BLAS, public Epetra_LAPACK{
116  public:
117 
119 
122 
124  virtual ~Epetra_SerialDenseSVD();
126 
128 
129 
131  int SetMatrix(Epetra_SerialDenseMatrix & A);
132 
134 
137  int SetVectors(Epetra_SerialDenseMatrix & X, Epetra_SerialDenseMatrix & B);
139 
141 
142 
144 
146 // void FactorWithEquilibration(bool Flag) {Equilibrate_ = Flag; return;};
147 
149  void SolveWithTranspose(bool Flag) {Transpose_ = Flag; if (Flag) TRANS_ = 'T'; else TRANS_ = 'N'; return;};
150 
152 // void SolveToRefinedSolution(bool Flag) {RefineSolution_ = Flag; return;};
153 
154  // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
155  // Causes all solves to estimate the forward and backward solution error.
156  /* Error estimates will be in the arrays FERR and BERR, resp, after the solve step is complete.
157  These arrays are accessible via the FERR() and BERR() access functions.
158  */
159 // void EstimateSolutionErrors(bool Flag) {EstimateSolutionErrors_ = Flag; return;};
161 
163 
164 
165  // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
166  // Computes the SVD factorization of the matrix using the LAPACK routine \e DGESVD.
167  /*
168  \return Integer error code, set to 0 if successful.
169  */
170 // virtual int Factor(void);
171  virtual int Factor(void);
172 
174 
177  virtual int Solve(void);
178 
180 
183  virtual int Invert( double rthresh = 0.0, double athresh = 0.0 );
184 
185  // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
186  // Computes the scaling vector S(i) = 1/sqrt(A(i,i) of the \e this matrix.
187  /*
188  \return Integer error code, set to 0 if successful. Otherwise returns the LAPACK error code INFO.
189  */
190 // virtual int ComputeEquilibrateScaling(void);
191 
192  // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
193  // Equilibrates the \e this matrix.
194  /*
195  \return Integer error code, set to 0 if successful. Otherwise returns the LAPACK error code INFO.
196  */
197 // virtual int EquilibrateMatrix(void);
198 
199  // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
200  // Equilibrates the current RHS.
201  /*
202  \return Integer error code, set to 0 if successful. Otherwise returns the LAPACK error code INFO.
203  */
204 // int EquilibrateRHS(void);
205 
206 
207  // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
208  // Apply Iterative Refinement.
209  /*
210  \return Integer error code, set to 0 if successful. Otherwise returns the LAPACK error code INFO.
211  */
212 // virtual int ApplyRefinement(void);
213 
214  // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
215  // Unscales the solution vectors if equilibration was used to solve the system.
216  /*
217  \return Integer error code, set to 0 if successful. Otherwise returns the LAPACK error code INFO.
218  */
219 // int UnequilibrateLHS(void);
220 
221  // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
222  // Returns the reciprocal of the 1-norm condition number of the \e this matrix.
223  /*
224  \param Value Out
225  On return contains the reciprocal of the 1-norm condition number of the \e this matrix.
226 
227  \return Integer error code, set to 0 if successful. Otherwise returns the LAPACK error code INFO.
228  */
229 // virtual int ReciprocalConditionEstimate(double & Value);
231 
233 
234 
236  bool Transpose() {return(Transpose_);};
237 
239  bool Factored() {return(Factored_);};
240 
241  // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
242  // Returns true if factor is equilibrated (factor available via AF() and LDAF()).
243 // bool A_Equilibrated() {return(A_Equilibrated_);};
244 
245  // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
246  // Returns true if RHS is equilibrated (RHS available via B() and LDB()).
247 // bool B_Equilibrated() {return(B_Equilibrated_);};
248 
249  // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
250  // Returns true if the LAPACK general rules for equilibration suggest you should equilibrate the system.
251 // virtual bool ShouldEquilibrate() {ComputeEquilibrateScaling(); return(ShouldEquilibrate_);};
252 
253  // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
254  // Returns true if forward and backward error estimated have been computed (available via FERR() and BERR()).
255 // bool SolutionErrorsEstimated() {return(SolutionErrorsEstimated_);};
256 
258  bool Inverted() {return(Inverted_);};
259 
260  // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
261  // Returns true if the condition number of the \e this matrix has been computed (value available via ReciprocalConditionEstimate()).
262 // bool ReciprocalConditionEstimated() {return(ReciprocalConditionEstimated_);};
263 
265  bool Solved() {return(Solved_);};
266 
267  // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
268  // Returns true if the current set of vectors has been refined.
269 // bool SolutionRefined() {return(SolutionRefined_);};
271 
273 
274 
276  Epetra_SerialDenseMatrix * Matrix() const {return(Matrix_);};
277 
278  // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
279  // Returns pointer to factored matrix (assuming factorization has been performed).
280 // Epetra_SerialDenseMatrix * FactoredMatrix() const {return(Factor_);};
281 
283  Epetra_SerialDenseMatrix * InvertedMatrix() const {return(Inverse_);};
284 
286  Epetra_SerialDenseMatrix * LHS() const {return(LHS_);};
287 
289  Epetra_SerialDenseMatrix * RHS() const {return(RHS_);};
290 
292  int M() const {return(M_);};
293 
295  int N() const {return(N_);};
296 
298  double * A() const {return(A_);};
299 
301  int LDA() const {return(LDA_);};
302 
304  double * B() const {return(B_);};
305 
307  int LDB() const {return(LDB_);};
308 
310  int NRHS() const {return(NRHS_);};
311 
313  double * X() const {return(X_);};
314 
316  int LDX() const {return(LDX_);};
317 
318  double * S() const {return(S_);};
319 
320  // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
321  // Returns pointer to the factored matrix (may be the same as A() if factorization done in place).
322 // double * AF() const {return(AF_);};
323 
324  // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
325  // Returns the leading dimension of the factored matrix.
326 // int LDAF() const {return(LDAF_);};
327 
329  double * AI() const {return(AI_);};
330 
332  int LDAI() const {return(LDAI_);};
333 
334  // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
335  // Returns pointer to pivot vector (if factorization has been computed), zero otherwise.
336 // int * IPIV() const {return(IPIV_);};
337 
339  double ANORM() const {return(ANORM_);};
340 
341  // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
342  // Returns the reciprocal of the condition number of the \e this matrix (returns -1 if not yet computed).
343 // double RCOND() const {return(RCOND_);};
344 
345  // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
346  // Ratio of smallest to largest row scale factors for the \e this matrix (returns -1 if not yet computed).
347  /* If ROWCND() is >= 0.1 and AMAX() is not close to overflow or underflow, then equilibration is not needed.
348  */
349 // double ROWCND() const {return(ROWCND_);};
350 
351  // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
352  // Ratio of smallest to largest column scale factors for the \e this matrix (returns -1 if not yet computed).
353  /* If COLCND() is >= 0.1 then equilibration is not needed.
354  */
355 // double COLCND() const {return(COLCND_);};
356 
357  // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
358  // Returns the absolute value of the largest entry of the \e this matrix (returns -1 if not yet computed).
359 // double AMAX() const {return(AMAX_);};
360 
361  // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
362  // Returns a pointer to the forward error estimates computed by LAPACK.
363 // double * FERR() const {return(FERR_);};
364 
365  // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
366  // Returns a pointer to the backward error estimates computed by LAPACK.
367 // double * BERR() const {return(BERR_);};
368 
369  // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
370  // Returns a pointer to the row scaling vector used for equilibration.
371 // double * R() const {return(R_);};
372 
373  // NOTE: doxygen-style documentation needs to be re-enabled if this function is re-enabled
374  // Returns a pointer to the column scale vector used for equilibration.
375 // double * C() const {return(C_);};
377 
379 
380  virtual void Print(std::ostream& os) const;
383 
385 
386 
388 
397  virtual int SetUseTranspose(bool use_transpose) { UseTranspose_ = use_transpose; return (0); }
398 
400 
409  { return Ymat.Multiply( UseTranspose_, false, 1.0, *Matrix(), Xmat, 0.0 ); }
410 
412 
422  { SetVectors(const_cast<Epetra_SerialDenseMatrix&>(Xmat),Ymat);
423  SolveWithTranspose(UseTranspose_);
424  return Solve(); }
425 
427  /* Returns the quantity \f$ \| A \|_\infty\f$ such that
428  \f[\| A \|_\infty = \max_{1\lei\lem} \sum_{j=1}^n |a_{ij}| \f].
429 
430  \warning This method must not be called unless HasNormInf() returns true.
431  */
432  virtual double NormInf() const { return Matrix()->NormInf(); }
433 
435  virtual const char * Label() const { return Epetra_Object::Label(); }
436 
438  virtual bool UseTranspose() const { return UseTranspose_; }
439 
441  virtual bool HasNormInf() const { return true; }
442 
444  virtual int RowDim() const { return M(); }
445 
447  virtual int ColDim() const { return N(); }
448 
450 
451  void AllocateWORK() {if (WORK_==0) {LWORK_ = 4*N_; WORK_ = new double[LWORK_];} return;};
452  void AllocateIWORK() {if (IWORK_==0) IWORK_ = new int[N_]; return;};
453  void InitPointers();
454  void DeleteArrays();
455  void ResetMatrix();
456  void ResetVectors();
457 
458 
459 // bool Equilibrate_;
460 // bool ShouldEquilibrate_;
461 // bool A_Equilibrated_;
462 // bool B_Equilibrated_;
464  bool Factored_;
465 // bool EstimateSolutionErrors_;
466 // bool SolutionErrorsEstimated_;
467  bool Solved_;
468  bool Inverted_;
469 // bool ReciprocalConditionEstimated_;
470 // bool RefineSolution_;
471 // bool SolutionRefined_;
472 
473  char TRANS_;
474 
475  int M_;
476  int N_;
477  int Min_MN_;
478  int NRHS_;
479  int LDA_;
480 // int LDAF_;
481  int LDAI_;
482  int LDB_;
483  int LDX_;
484  int INFO_;
485  int LWORK_;
486 
487 // int * IPIV_;
488  int * IWORK_;
489 
490  double ANORM_;
491 // double RCOND_;
492 // double ROWCND_;
493 // double COLCND_;
494 // double AMAX_;
495 
499 // Epetra_SerialDenseMatrix * Factor_;
501 
502  double * A_;
503 // double * FERR_;
504 // double * BERR_;
505 // double * AF_;
506  double * AI_;
507  double * WORK_;
508 // double * R_;
509 // double * C_;
510  double * U_;
511  double * S_;
512  double * Vt_;
513 
514  double * B_;
515  double * X_;
516 
518 
519  private:
520  // Epetra_SerialDenseSolver copy constructor (put here because we don't want user access)
521 
524 };
525 
526 #endif /* _EPETRA_SERIALDENSESVD_H_ */
double * B() const
Returns pointer to current RHS.
virtual void Print(std::ostream &os) const
Print object to an output stream Print method.
virtual bool UseTranspose() const
Returns the current UseTranspose setting.
double * A() const
Returns pointer to the this matrix.
Epetra_SerialDenseMatrix * Inverse_
Epetra_SerialDenseMatrix * LHS() const
Returns pointer to current LHS.
Epetra_SerialDenseSVD: A class for SVDing dense linear problems.
double * X() const
Returns pointer to current solution.
Epetra_CompObject & operator=(const Epetra_CompObject &src)
Epetra_SerialDenseMatrix: A class for constructing and using real double precision general dense matr...
double * AI() const
Returns pointer to the inverted matrix (may be the same as A() if factorization done in place)...
int NRHS() const
Returns the number of current right hand sides and solution vectors.
bool Solved()
Returns true if the current set of vectors has been solved.
Epetra_SerialDenseMatrix * Matrix_
virtual int RowDim() const
Returns the row dimension of operator.
int N() const
Returns column dimension of system.
double ANORM() const
Returns the 1-Norm of the this matrix (returns -1 if not yet computed).
int LDA() const
Returns the leading dimension of the this matrix.
Epetra_SerialDenseMatrix * RHS() const
Returns pointer to current RHS.
Epetra_SerialDenseMatrix * InvertedMatrix() const
Returns pointer to inverted matrix (assuming inverse has been performed).
Epetra_BLAS: The Epetra BLAS Wrapper Class.
Definition: Epetra_BLAS.h:70
virtual const char * Label() const
Epetra_Object Label access funtion.
Epetra_Object: The base Epetra class.
Definition: Epetra_Object.h:57
Epetra_SerialDenseMatrix * RHS_
Epetra_CompObject: Functionality and data that is common to all computational classes.
virtual const char * Label() const
Returns a character string describing the operator.
Epetra_LAPACK: The Epetra LAPACK Wrapper Class.
Definition: Epetra_LAPACK.h:67
Epetra_SerialDenseOperator: A pure virtual class for using real-valued double-precision operators...
virtual bool HasNormInf() const
Returns true if the this object can provide an approximate Inf-norm, false otherwise.
int Multiply(char TransA, char TransB, double ScalarAB, const Epetra_SerialDenseMatrix &A, const Epetra_SerialDenseMatrix &B, double ScalarThis)
Matrix-Matrix multiplication, this = ScalarThis*this + ScalarAB*A*B.
int LDX() const
Returns the leading dimension of the solution.
bool Inverted()
Returns true if matrix inverse has been computed (inverse available via AF() and LDAF()).
virtual int ApplyInverse(const Epetra_SerialDenseMatrix &Xmat, Epetra_SerialDenseMatrix &Ymat)
Returns the result of a Epetra_SerialDenseOperator inverse applied to an Epetra_SerialDenseMatrix X i...
Epetra_SerialDenseMatrix * Matrix() const
Returns pointer to current matrix.
virtual int Apply(const Epetra_SerialDenseMatrix &Xmat, Epetra_SerialDenseMatrix &Ymat)
Returns the result of a Epetra_SerialDenseOperator applied to a Epetra_SerialDenseMatrix X in Y...
int Solve(int, TYPE *, TYPE *, TYPE *)
void SolveWithTranspose(bool Flag)
Causes equilibration to be called just before the matrix factorization as part of the call to Factor...
virtual int ColDim() const
Returns the column dimension of operator.
bool Transpose()
Returns true if transpose of this matrix has and will be used.
virtual double NormInf() const
Returns the infinity norm of the global matrix.
int LDAI() const
Returns the leading dimension of the inverted matrix.
Epetra_SerialDenseMatrix * LHS_
int M() const
Returns row dimension of system.
virtual int SetUseTranspose(bool use_transpose)
If set true, transpose of this operator will be applied.
bool Factored()
Returns true if matrix is factored (factor available via AF() and LDAF()).
int LDB() const
Returns the leading dimension of the RHS.