Epetra Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
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//=========================================================================
115class 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
120
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
381 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_;
465// bool EstimateSolutionErrors_;
466// bool SolutionErrorsEstimated_;
469// bool ReciprocalConditionEstimated_;
470// bool RefineSolution_;
471// bool SolutionRefined_;
472
473 char TRANS_;
474
475 int M_;
476 int N_;
478 int NRHS_;
479 int LDA_;
480// int LDAF_;
481 int LDAI_;
482 int LDB_;
483 int LDX_;
484 int INFO_;
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_ */
Epetra_BLAS: The Epetra BLAS Wrapper Class.
Definition Epetra_BLAS.h:70
Epetra_CompObject: Functionality and data that is common to all computational classes.
Epetra_LAPACK: The Epetra LAPACK Wrapper Class.
Epetra_Object: The base Epetra class.
virtual void Print(std::ostream &os) const
Print object to an output stream Print method.
virtual const char * Label() const
Epetra_Object Label access funtion.
Epetra_SerialDenseMatrix: A class for constructing and using real double precision general dense matr...
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.
Epetra_SerialDenseOperator: A pure virtual class for using real-valued double-precision operators.
Epetra_SerialDenseSVD: A class for SVDing dense linear problems.
double * A() const
Returns pointer to the this matrix.
virtual const char * Label() const
Returns a character string describing the operator.
Epetra_SerialDenseMatrix * RHS_
Epetra_SerialDenseMatrix * LHS() const
Returns pointer to current LHS.
Epetra_SerialDenseMatrix * Inverse_
Epetra_SerialDenseMatrix * Matrix_
bool Transpose()
Returns true if transpose of this matrix has and will be used.
virtual bool HasNormInf() const
Returns true if the this object can provide an approximate Inf-norm, false otherwise.
int LDX() const
Returns the leading dimension of the solution.
virtual int RowDim() const
Returns the row dimension of operator.
Epetra_SerialDenseMatrix * Matrix() const
Returns pointer to current matrix.
double ANORM() const
Returns the 1-Norm of the this matrix (returns -1 if not yet computed).
Epetra_SerialDenseSVD & operator=(const Epetra_SerialDenseSVD &Source)
double * AI() const
Returns pointer to the inverted matrix (may be the same as A() if factorization done in place).
Epetra_SerialDenseMatrix * RHS() const
Returns pointer to current RHS.
double * B() const
Returns pointer to current RHS.
int NRHS() const
Returns the number of current right hand sides and solution vectors.
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...
int M() const
Returns row dimension of system.
int N() const
Returns column dimension of system.
Epetra_SerialDenseMatrix * LHS_
virtual double NormInf() const
Returns the infinity norm of the global matrix.
Epetra_SerialDenseMatrix * InvertedMatrix() const
Returns pointer to inverted matrix (assuming inverse has been performed).
double * X() const
Returns pointer to current solution.
virtual int ColDim() const
Returns the column dimension of operator.
void SolveWithTranspose(bool Flag)
Causes equilibration to be called just before the matrix factorization as part of the call to Factor.
bool Inverted()
Returns true if matrix inverse has been computed (inverse available via AF() and LDAF()).
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.
Epetra_SerialDenseSVD(const Epetra_SerialDenseSVD &Source)
int LDA() const
Returns the leading dimension of the this matrix.
int LDB() const
Returns the leading dimension of the RHS.
virtual int SetUseTranspose(bool use_transpose)
If set true, transpose of this operator will be applied.
virtual bool UseTranspose() const
Returns the current UseTranspose setting.
int LDAI() const
Returns the leading dimension of the inverted matrix.
bool Factored()
Returns true if matrix is factored (factor available via AF() and LDAF()).
bool Solved()
Returns true if the current set of vectors has been solved.