Epetra Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
Epetra_SerialDenseSolver.cpp
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
45
46//=============================================================================
51 Equilibrate_(false),
52 ShouldEquilibrate_(false),
53 A_Equilibrated_(false),
54 B_Equilibrated_(false),
55 Transpose_(false),
56 Factored_(false),
57 EstimateSolutionErrors_(false),
58 SolutionErrorsEstimated_(false),
59 Solved_(false),
60 Inverted_(false),
61 ReciprocalConditionEstimated_(false),
62 RefineSolution_(false),
63 SolutionRefined_(false),
64 TRANS_('N'),
65 M_(0),
66 N_(0),
67 Min_MN_(0),
68 NRHS_(0),
69 LDA_(0),
70 LDAF_(0),
71 LDB_(0),
72 LDX_(0),
73 INFO_(0),
74 LWORK_(0),
75 IPIV_(0),
76 IWORK_(0),
77 ANORM_(0.0),
78 RCOND_(0.0),
79 ROWCND_(0.0),
80 COLCND_(0.0),
81 AMAX_(0.0),
82 Matrix_(0),
83 LHS_(0),
84 RHS_(0),
85 Factor_(0),
86 A_(0),
87 FERR_(0),
88 BERR_(0),
89 AF_(0),
90 WORK_(0),
91 R_(0),
92 C_(0),
93 B_(0),
94 X_(0)
95{
99}
100//=============================================================================
105//=============================================================================
107{
108 IWORK_ = 0;
109 FERR_ = 0;
110 BERR_ = 0;
111 Factor_ =0;
112 Matrix_ =0;
113 AF_ = 0;
114 IPIV_ = 0;
115 WORK_ = 0;
116 R_ = 0;
117 C_ = 0;
118 INFO_ = 0;
119 LWORK_ = 0;
120}
121//=============================================================================
123{
124 if (IWORK_ != 0) {delete [] IWORK_; IWORK_ = 0;}
125 if (FERR_ != 0) {delete [] FERR_; FERR_ = 0;}
126 if (BERR_ != 0) {delete [] BERR_; BERR_ = 0;}
127 if (Factor_ != Matrix_ && Factor_ != 0) {delete Factor_; Factor_ = 0;}
128 if (Factor_ !=0) Factor_ = 0;
129 if (AF_ !=0) AF_ = 0;
130 if (IPIV_ != 0) {delete [] IPIV_;IPIV_ = 0;}
131 if (WORK_ != 0) {delete [] WORK_;WORK_ = 0;}
132 if (R_ != 0 && R_ != C_) {delete [] R_;R_ = 0;}
133 if (R_ != 0) R_ = 0;
134 if (C_ != 0) {delete [] C_;C_ = 0;}
135 INFO_ = 0;
136 LWORK_ = 0;
137}
138//=============================================================================
140{
141 DeleteArrays();
142 ResetVectors();
143 Matrix_ = 0;
144 Factor_ = 0;
145 A_Equilibrated_ = false;
146 Factored_ = false;
147 Inverted_ = false;
148 M_ = 0;
149 N_ = 0;
150 Min_MN_ = 0;
151 LDA_ = 0;
152 LDAF_ = 0;
153 ANORM_ = -1.0;
154 RCOND_ = -1.0;
155 ROWCND_ = -1.0;
156 COLCND_ = -1.0;
157 AMAX_ = -1.0;
158 A_ = 0;
159
160}
161//=============================================================================
163 ResetMatrix();
164 Matrix_ = &A_in;
165 Factor_ = &A_in;
166 M_ = A_in.M();
167 N_ = A_in.N();
169 LDA_ = A_in.LDA();
170 LDAF_ = LDA_;
171 A_ = A_in.A();
172 AF_ = A_in.A();
173 return(0);
174}
175//=============================================================================
177{
178 LHS_ = 0;
179 RHS_ = 0;
180 B_ = 0;
181 X_ = 0;
183 SolutionRefined_ = false;
184 Solved_ = false;
186 B_Equilibrated_ = false;
187 NRHS_ = 0;
188 LDB_ = 0;
189 LDX_ = 0;
190}
191//=============================================================================
193{
194 if (B_in.M()!=X_in.M() || B_in.N() != X_in.N()) EPETRA_CHK_ERR(-1);
195 if (B_in.A()==0) EPETRA_CHK_ERR(-2);
196 if (B_in.LDA()<1) EPETRA_CHK_ERR(-3);
197 if (X_in.A()==0) EPETRA_CHK_ERR(-4);
198 if (X_in.LDA()<1) EPETRA_CHK_ERR(-5);
199
200 ResetVectors();
201 LHS_ = &X_in;
202 RHS_ = &B_in;
203 NRHS_ = B_in.N();
204
205 B_ = B_in.A();
206 LDB_ = B_in.LDA();
207 X_ = X_in.A();
208 LDX_ = X_in.LDA();
209 return(0);
210}
211//=============================================================================
214 // If the errors are estimated, this implies that the solution must be refined
216 return;
217}
218//=============================================================================
220 if (Factored()) return(0); // Already factored
221 if (Inverted()) EPETRA_CHK_ERR(-100); // Cannot factor inverted matrix
222 int ierr = 0;
223
224 ANORM_ = Matrix_->OneNorm(); // Compute 1-Norm of A
225
226
227 // If we want to refine the solution, then the factor must
228 // be stored separatedly from the original matrix
229
230 if (A_ == AF_)
231 if (RefineSolution_ ) {
233 AF_ = Factor_->A();
234 LDAF_ = Factor_->LDA();
235 }
236
237 if (Equilibrate_) ierr = EquilibrateMatrix();
238
239 if (ierr!=0) EPETRA_CHK_ERR(ierr-2);
240
241 if (IPIV_==0) IPIV_ = new int[Min_MN_]; // Allocated Pivot vector if not already done.
242
243 GETRF (M_, N_, AF_, LDAF_, IPIV_, &INFO_);
244
245 Factored_ = true;
246 double DN = N_;
247 UpdateFlops(2.0*(DN*DN*DN)/3.0);
248
250 return(0);
251
252}
253
254//=============================================================================
256 int ierr = 0;
257
258 // We will call one of four routines depending on what services the user wants and
259 // whether or not the matrix has been inverted or factored already.
260 //
261 // If the matrix has been inverted, use DGEMM to compute solution.
262 // Otherwise, if the user want the matrix to be equilibrated or wants a refined solution, we will
263 // call the X interface.
264 // Otherwise, if the matrix is already factored we will call the TRS interface.
265 // Otherwise, if the matrix is unfactored we will call the SV interface.
266
267
268
269 if (Equilibrate_) {
270 ierr = EquilibrateRHS();
271 B_Equilibrated_ = true;
272 }
273 EPETRA_CHK_ERR(ierr);
274 if (A_Equilibrated_ && !B_Equilibrated_) EPETRA_CHK_ERR(-1); // Matrix and vectors must be similarly scaled
276 if (B_==0) EPETRA_CHK_ERR(-3); // No B
277 if (X_==0) EPETRA_CHK_ERR(-4); // No X
278
279 if (ShouldEquilibrate() && !A_Equilibrated_) ierr = 1; // Warn that the system should be equilibrated.
280
281 double DN = N_;
282 double DNRHS = NRHS_;
283 if (Inverted()) {
284
285 if (B_==X_) EPETRA_CHK_ERR(-100); // B and X must be different for this case
286
287 GEMM(TRANS_, 'N', N_, NRHS_, N_, 1.0, AF_, LDAF_,
288 B_, LDB_, 0.0, X_, LDX_);
289 if (INFO_!=0) EPETRA_CHK_ERR(INFO_);
290 UpdateFlops(2.0*DN*DN*DNRHS);
291 Solved_ = true;
292 }
293 else {
294
295 if (!Factored()) Factor(); // Matrix must be factored
296
297 if (B_!=X_) {
298 *LHS_ = *RHS_; // Copy B to X if needed
299 X_ = LHS_->A(); LDX_ = LHS_->LDA();
300 }
302 if (INFO_!=0) EPETRA_CHK_ERR(INFO_);
303 UpdateFlops(2.0*DN*DN*DNRHS);
304 Solved_ = true;
305
306 }
307 int ierr1=0;
308 if (RefineSolution_ && !Inverted()) ierr1 = ApplyRefinement();
309 if (ierr1!=0) EPETRA_CHK_ERR(ierr1)
310 else
311 EPETRA_CHK_ERR(ierr);
312
313 if (Equilibrate_) ierr1 = UnequilibrateLHS();
314 EPETRA_CHK_ERR(ierr1);
315 return(0);
316}
317//=============================================================================
319{
320 double DN = N_;
321 double DNRHS = NRHS_;
322 if (!Solved()) EPETRA_CHK_ERR(-100); // Must have an existing solution
323 if (A_==AF_) EPETRA_CHK_ERR(-101); // Cannot apply refine if no original copy of A.
324
325 if (FERR_ != 0) delete [] FERR_; // Always start with a fresh copy of FERR_, since NRHS_ may change
326 FERR_ = new double[NRHS_];
327 if (BERR_ != 0) delete [] BERR_; // Always start with a fresh copy of BERR_, since NRHS_ may change
328 BERR_ = new double[NRHS_];
329 AllocateWORK();
331
333 B_, LDB_, X_, LDX_, FERR_, BERR_,
334 WORK_, IWORK_, &INFO_);
335
336
339 SolutionRefined_ = true;
340
341 UpdateFlops(2.0*DN*DN*DNRHS); // Not sure of count
342
344 return(0);
345
346}
347
348//=============================================================================
350 if (R_!=0) return(0); // Already computed
351
352 double DM = M_;
353 double DN = N_;
354 R_ = new double[M_];
355 C_ = new double[N_];
356
357 GEEQU (M_, N_, AF_, LDAF_, R_, C_, &ROWCND_, &COLCND_, &AMAX_, &INFO_);
358 if (INFO_ != 0) EPETRA_CHK_ERR(INFO_);
359
360 if (COLCND_<0.1 || ROWCND_<0.1 || AMAX_ < Epetra_Underflow || AMAX_ > Epetra_Overflow) ShouldEquilibrate_ = true;
361
362 UpdateFlops(4.0*DM*DN);
363
364 return(0);
365}
366
367//=============================================================================
369{
370 int i, j;
371 int ierr = 0;
372
373 double DN = N_;
374 double DM = M_;
375
376 if (A_Equilibrated_) return(0); // Already done
377 if (R_==0) ierr = ComputeEquilibrateScaling(); // Compute R and C if needed
378 if (ierr!=0) EPETRA_CHK_ERR(ierr);
379 if (A_==AF_) {
380 double * ptr;
381 for (j=0; j<N_; j++) {
382 ptr = A_ + j*LDA_;
383 double s1 = C_[j];
384 for (i=0; i<M_; i++) {
385 *ptr = *ptr*s1*R_[i];
386 ptr++;
387 }
388 }
389 UpdateFlops(2.0*DM*DN);
390 }
391 else {
392 double * ptr;
393 double * ptr1;
394 for (j=0; j<N_; j++) {
395 ptr = A_ + j*LDA_;
396 ptr1 = AF_ + j*LDAF_;
397 double s1 = C_[j];
398 for (i=0; i<M_; i++) {
399 *ptr = *ptr*s1*R_[i];
400 ptr++;
401 *ptr1 = *ptr1*s1*R_[i];
402 ptr1++;
403 }
404 }
405 UpdateFlops(4.0*DM*DN);
406 }
407
408 A_Equilibrated_ = true;
409
410 return(0);
411}
412
413//=============================================================================
415{
416 int i, j;
417 int ierr = 0;
418
419 if (B_Equilibrated_) return(0); // Already done
420 if (R_==0) ierr = ComputeEquilibrateScaling(); // Compute R and C if needed
421 if (ierr!=0) EPETRA_CHK_ERR(ierr);
422
423 double * R_tmp = R_;
424 if (Transpose_) R_tmp = C_;
425
426 double * ptr;
427 for (j=0; j<NRHS_; j++) {
428 ptr = B_ + j*LDB_;
429 for (i=0; i<M_; i++) {
430 *ptr = *ptr*R_tmp[i];
431 ptr++;
432 }
433 }
434
435
436 B_Equilibrated_ = true;
437 UpdateFlops((double) N_*(double) NRHS_);
438
439 return(0);
440}
441
442//=============================================================================
444{
445 int i, j;
446
447 if (!B_Equilibrated_) return(0); // Nothing to do
448
449 double * C_tmp = C_;
450 if (Transpose_) C_tmp = R_;
451
452 double * ptr;
453 for (j=0; j<NRHS_; j++) {
454 ptr = X_ + j*LDX_;
455 for (i=0; i<N_; i++) {
456 *ptr = *ptr*C_tmp[i];
457 ptr++;
458 }
459 }
460
461
462 UpdateFlops((double) N_ *(double) NRHS_);
463
464 return(0);
465}
466
467//=============================================================================
469{
470 if (!Factored()) Factor(); // Need matrix factored.
471
472 /* This section work with LAPACK Version 3.0 only
473 // Setting LWORK = -1 and calling GETRI will return optimal work space size in WORK_TMP
474 int LWORK_TMP = -1;
475 double WORK_TMP;
476 GETRI ( N_, AF_, LDAF_, IPIV_, &WORK_TMP, &LWORK_TMP, &INFO_);
477 LWORK_TMP = WORK_TMP; // Convert to integer
478 if (LWORK_TMP>LWORK_) {
479 if (WORK_!=0) delete [] WORK_;
480 LWORK_ = LWORK_TMP;
481 WORK_ = new double[LWORK_];
482 }
483 */
484 // This section will work with any version of LAPACK
485 AllocateWORK();
486
487 GETRI ( N_, AF_, LDAF_, IPIV_, WORK_, &LWORK_, &INFO_);
488
489 double DN = N_;
490 UpdateFlops((DN*DN*DN));
491 Inverted_ = true;
492 Factored_ = false;
493
495 return(0);
496}
497
498//=============================================================================
500{
501 int ierr = 0;
503 Value = RCOND_;
504 return(0); // Already computed, just return it.
505 }
506
507 if (ANORM_<0.0) ANORM_ = Matrix_->OneNorm();
508 if (!Factored()) ierr = Factor(); // Need matrix factored.
509 if (ierr!=0) EPETRA_CHK_ERR(ierr-2);
510
511 AllocateWORK();
513 // We will assume a one-norm condition number
514 GECON( '1', N_, AF_, LDAF_, ANORM_, &RCOND_, WORK_, IWORK_, &INFO_);
516 Value = RCOND_;
517 UpdateFlops(2*N_*N_); // Not sure of count
519 return(0);
520}
521//=============================================================================
522void Epetra_SerialDenseSolver::Print(std::ostream& os) const {
523
524 if (Matrix_!=0) os << "Solver Matrix" << std::endl << *Matrix_ << std::endl;
525 if (Factor_!=0) os << "Solver Factored Matrix" << std::endl << *Factor_ << std::endl;
526 if (LHS_ !=0) os << "Solver LHS" << std::endl << *LHS_ << std::endl;
527 if (RHS_ !=0) os << "Solver RHS" << std::endl << *RHS_ << std::endl;
528
529}
const double Epetra_Overflow
#define EPETRA_MIN(x, y)
#define EPETRA_CHK_ERR(a)
Epetra_BLAS: The Epetra BLAS Wrapper Class.
Definition Epetra_BLAS.h:70
void GEMM(const char TRANSA, const char TRANSB, const int M, const int N, const int K, const float ALPHA, const float *A, const int LDA, const float *B, const int LDB, const float BETA, float *C, const int LDC) const
Epetra_BLAS matrix-matrix multiply function (SGEMM)
Epetra_CompObject: Functionality and data that is common to all computational classes.
void UpdateFlops(int Flops_in) const
Increment Flop count for this object.
Epetra_LAPACK: The Epetra LAPACK Wrapper Class.
void GEEQU(const int M, const int N, const float *A, const int LDA, float *R, float *C, float *ROWCND, float *COLCND, float *AMAX, int *INFO) const
Epetra_LAPACK equilibration for general matrix (SGEEQU)
void GECON(const char NORM, const int N, const float *A, const int LDA, const float ANORM, float *RCOND, float *WORK, int *IWORK, int *INFO) const
Epetra_LAPACK condition number estimator for general matrix (SGECON)
void GERFS(const char TRANS, const int N, const int NRHS, const float *A, const int LDA, const float *AF, const int LDAF, const int *IPIV, const float *B, const int LDB, float *X, const int LDX, float *FERR, float *BERR, float *WORK, int *IWORK, int *INFO) const
Epetra_LAPACK Refine solution (GERFS)
void GETRF(const int M, const int N, float *A, const int LDA, int *IPIV, int *INFO) const
Epetra_LAPACK factorization for general matrix (SGETRF)
void GETRI(const int N, float *A, const int LDA, int *IPIV, float *WORK, const int *LWORK, int *INFO) const
Epetra_LAPACK inversion for general matrix (SGETRI)
void GETRS(const char TRANS, const int N, const int NRHS, const float *A, const int LDA, const int *IPIV, float *X, const int LDX, int *INFO) const
Epetra_LAPACK solve (after factorization) for general matrix (SGETRS)
Epetra_SerialDenseMatrix: A class for constructing and using real double precision general dense matr...
double * A() const
Returns pointer to the this matrix.
int LDA() const
Returns the leading dimension of the this matrix.
virtual double OneNorm() const
Computes the 1-Norm of the this matrix (identical to NormOne() method).
int M() const
Returns row dimension of system.
int N() const
Returns column dimension of system.
bool ReciprocalConditionEstimated()
Returns true if the condition number of the this matrix has been computed (value available via Recipr...
virtual int Solve(void)
Computes the solution X to AX = B for the this matrix and the B provided to SetVectors()....
int EquilibrateRHS(void)
Equilibrates the current RHS.
Epetra_SerialDenseMatrix * Factor_
virtual void Print(std::ostream &os) const
Print service methods; defines behavior of ostream << operator.
int SetMatrix(Epetra_SerialDenseMatrix &A)
Sets the pointers for coefficient matrix.
Epetra_SerialDenseMatrix * Matrix_
virtual int ReciprocalConditionEstimate(double &Value)
Returns the reciprocal of the 1-norm condition number of the this matrix.
virtual int ApplyRefinement(void)
Apply Iterative Refinement.
void EstimateSolutionErrors(bool Flag)
Causes all solves to estimate the forward and backward solution error.
virtual ~Epetra_SerialDenseSolver()
Epetra_SerialDenseSolver destructor.
Epetra_SerialDenseMatrix * RHS_
int UnequilibrateLHS(void)
Unscales the solution vectors if equilibration was used to solve the system.
virtual int ComputeEquilibrateScaling(void)
Computes the scaling vector S(i) = 1/sqrt(A(i,i)) of the this matrix.
virtual bool ShouldEquilibrate()
Returns true if the LAPACK general rules for equilibration suggest you should equilibrate the system.
virtual int Factor(void)
Computes the in-place LU factorization of the matrix using the LAPACK routine DGETRF.
bool Factored()
Returns true if matrix is factored (factor available via AF() and LDAF()).
bool Inverted()
Returns true if matrix inverse has been computed (inverse available via AF() and LDAF()).
Epetra_SerialDenseSolver()
Default constructor; matrix should be set using SetMatrix(), LHS and RHS set with SetVectors().
Epetra_SerialDenseMatrix * LHS_
virtual int EquilibrateMatrix(void)
Equilibrates the this matrix.
virtual int Invert(void)
Inverts the this matrix.
bool Solved()
Returns true if the current set of vectors has been solved.
int SetVectors(Epetra_SerialDenseMatrix &X, Epetra_SerialDenseMatrix &B)
Sets the pointers for left and right hand side vector(s).