Epetra Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Epetra_SerialSpdDenseSolver.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 
46 
47 //=============================================================================
50  SCOND_(-1.0),
51  SymMatrix_(0),
52  SymFactor_(0)
53 {
54 }
55 //=============================================================================
57 {
58  if (SymFactor_ != SymMatrix_ && SymFactor_ != 0) {
59  delete SymFactor_; SymFactor_ = 0; Factor_ = 0;
60  }
61 }
62 //=============================================================================
64 
65  SymMatrix_=&A_in;
66  SymFactor_=&A_in;
67  SCOND_ = -1.0;
68  // Also call SerialDensematrix set method
70 }
71 //=============================================================================
73  if (Factored()) return(0); // Already factored
74  if (Inverted()) EPETRA_CHK_ERR(-100); // Cannot factor inverted matrix
75  int ierr = 0;
76 
78 
79 
80  // If we want to refine the solution, then the factor must
81  // be stored separatedly from the original matrix
82 
83  if (A_ == AF_)
84  if (RefineSolution_ ) {
87  AF_ = SymFactor_->A();
88  LDAF_ = SymFactor_->LDA();
89  }
90  if (Equilibrate_) ierr = EquilibrateMatrix();
91 
92  if (ierr!=0) EPETRA_CHK_ERR(ierr-2);
93 
94  POTRF (SymMatrix_->UPLO(), N_, AF_, LDAF_, &INFO_);
95  Factored_ = true;
96  double DN = N_;
97  UpdateFlops((DN*DN*DN)/3.0);
98 
100  return(0);
101 
102 }
103 
104 //=============================================================================
106  int ierr = 0;
107 
108  // We will call one of four routines depending on what services the user wants and
109  // whether or not the matrix has been inverted or factored already.
110  //
111  // If the matrix has been inverted, use DGEMM to compute solution.
112  // Otherwise, if the user want the matrix to be equilibrated or wants a refined solution, we will
113  // call the X interface.
114  // Otherwise, if the matrix is already factored we will call the TRS interface.
115  // Otherwise, if the matrix is unfactored we will call the SV interface.
116 
117  if (Equilibrate_) {
119  B_Equilibrated_ = true;
120  }
121  EPETRA_CHK_ERR(ierr);
122  if (A_Equilibrated_ && !B_Equilibrated_) EPETRA_CHK_ERR(-1); // Matrix and vectors must be similarly scaled
124  if (B_==0) EPETRA_CHK_ERR(-3); // No B
125  if (X_==0) EPETRA_CHK_ERR(-4); // No B
126 
127  if (ShouldEquilibrate() && !A_Equilibrated_) ierr = 1; // Warn that the system should be equilibrated.
128 
129  double DN = N_;
130  double DNRHS = NRHS_;
131  if (Inverted()) {
132 
133  if (B_==X_) EPETRA_CHK_ERR(-100); // B and X must be different for this case
134  GEMM('N', 'N', N_, NRHS_, N_, 1.0, AF_, LDAF_,
135  B_, LDB_, 0.0, X_, LDX_);
136  if (INFO_!=0) EPETRA_CHK_ERR(INFO_);
137  UpdateFlops(2.0*DN*DN*DNRHS);
138  Solved_ = true;
139  }
140  else {
141 
142  if (!Factored()) Factor(); // Matrix must be factored
143  if (B_!=X_) {
144  *LHS_ = *RHS_; // Copy B to X if needed
145  X_ = LHS_->A(); LDX_ = LHS_->LDA();
146  }
147 
149  if (INFO_!=0) EPETRA_CHK_ERR(INFO_);
150  UpdateFlops(2.0*DN*DN*DNRHS);
151  Solved_ = true;
152 
153  }
154  int ierr1=0;
155  if (RefineSolution_) ierr1 = ApplyRefinement();
156  if (ierr1!=0) {
157  EPETRA_CHK_ERR(ierr1);
158  }
159  else {
160  EPETRA_CHK_ERR(ierr);
161  }
163  EPETRA_CHK_ERR(ierr1);
164  return(0);
165 }
166 //=============================================================================
168 {
169  double DN = N_;
170  double DNRHS = NRHS_;
171  if (!Solved()) EPETRA_CHK_ERR(-100); // Must have an existing solution
172  if (A_==AF_) EPETRA_CHK_ERR(-101); // Cannot apply refine if no original copy of A.
173 
174  if (FERR_ != 0) delete [] FERR_; // Always start with a fresh copy of FERR_, since NRHS_ may change
175  FERR_ = new double[NRHS_];
176  if (BERR_ != 0) delete [] BERR_; // Always start with a fresh copy of FERR_, since NRHS_ may change
177  BERR_ = new double[NRHS_];
178  AllocateWORK();
179  AllocateIWORK();
180 
182  B_, LDB_, X_, LDX_, FERR_, BERR_,
183  WORK_, IWORK_, &INFO_);
184 
185 
188  SolutionRefined_ = true;
189 
190  UpdateFlops(2.0*DN*DN*DNRHS); // Not sure of count
191 
193  return(0);
194 
195 }
196 
197 //=============================================================================
199  if (R_!=0) return(0); // Already computed
200 
201  double DN = N_;
202  R_ = new double[N_];
203  C_ = R_;
204 
205  POEQU (N_, AF_, LDAF_, R_, &SCOND_, &AMAX_, &INFO_);
206  if (INFO_ != 0) EPETRA_CHK_ERR(INFO_);
207 
208  if (SCOND_<0.1 || AMAX_ < Epetra_Underflow || AMAX_ > Epetra_Overflow) ShouldEquilibrate_ = true;
209 
210  C_ = R_; // Set column scaling pointer so we can use EquilibrateRHS and UnequilibrateLHS from base class
211  UpdateFlops(2.0*DN*DN);
212 
213  return(0);
214 }
215 
216 //=============================================================================
218  int i, j;
219  int ierr = 0;
220 
221  if (A_Equilibrated_) return(0); // Already done
222  if (R_==0) ierr = ComputeEquilibrateScaling(); // Compute S if needed
223  if (ierr!=0) EPETRA_CHK_ERR(ierr);
224  if (SymMatrix_->Upper()) {
225  if (A_==AF_) {
226  double * ptr;
227  for (j=0; j<N_; j++) {
228  ptr = A_ + j*LDA_;
229  double s1 = R_[j];
230  for (i=0; i<=j; i++) {
231  *ptr = *ptr*s1*R_[i];
232  ptr++;
233  }
234  }
235  }
236  else {
237  double * ptr;
238  double * ptr1;
239  for (j=0; j<N_; j++) {
240  ptr = A_ + j*LDA_;
241  ptr1 = AF_ + j*LDAF_;
242  double s1 = R_[j];
243  for (i=0; i<=j; i++) {
244  *ptr = *ptr*s1*R_[i];
245  ptr++;
246  *ptr1 = *ptr1*s1*R_[i];
247  ptr1++;
248  }
249  }
250  }
251  }
252  else {
253  if (A_==AF_) {
254  double * ptr;
255  for (j=0; j<N_; j++) {
256  ptr = A_ + j + j*LDA_;
257  double s1 = R_[j];
258  for (i=j; i<N_; i++) {
259  *ptr = *ptr*s1*R_[i];
260  ptr++;
261  }
262  }
263  }
264  else {
265  double * ptr;
266  double * ptr1;
267  for (j=0; j<N_; j++) {
268  ptr = A_ + j + j*LDA_;
269  ptr1 = AF_ + j + j*LDAF_;
270  double s1 = R_[j];
271  for (i=j; i<N_; i++) {
272  *ptr = *ptr*s1*R_[i];
273  ptr++;
274  *ptr1 = *ptr1*s1*R_[i];
275  ptr1++;
276  }
277  }
278  }
279  }
280  A_Equilibrated_ = true;
281  double NumFlops = (double) ((N_+1)*N_/2);
282  if (A_==AF_) NumFlops += NumFlops;
283  UpdateFlops(NumFlops);
284 
285  return(0);
286 }
287 
288 //=============================================================================
290 {
291  if (!Factored()) Factor(); // Need matrix factored.
292  POTRI ( SymMatrix_->UPLO(), N_, AF_, LDAF_, &INFO_);
293  // Copy lower/upper triangle to upper/lower triangle: make full inverse
295  double DN = N_;
296  UpdateFlops((DN*DN*DN));
297  Inverted_ = true;
298  Factored_ = false;
299 
301  return(0);
302 }
303 
304 //=============================================================================
306 {
307  int ierr = 0;
309  Value = RCOND_;
310  return(0); // Already computed, just return it.
311  }
312 
313  if (ANORM_<0.0) ANORM_ = SymMatrix_->OneNorm();
314  if (ierr!=0) EPETRA_CHK_ERR(ierr-1);
315  if (!Factored()) ierr = Factor(); // Need matrix factored.
316  if (ierr!=0) EPETRA_CHK_ERR(ierr-2);
317 
318  AllocateWORK();
319  AllocateIWORK();
322  Value = RCOND_;
323  UpdateFlops(2*N_*N_); // Not sure of count
325  return(0);
326 }
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_SerialSymDenseMatrix * SymFactor_
int Invert(void)
Inverts the this matrix.
double OneNorm() const
Computes the 1-Norm of the this matrix (identical to NormOne() method).
Epetra_SerialSpdDenseSolver()
Default constructor; matrix should be set using SetMatrix(), LHS and RHS set with SetVectors()...
virtual ~Epetra_SerialSpdDenseSolver()
Epetra_SerialDenseSolver destructor.
#define EPETRA_CHK_ERR(a)
int SetMatrix(Epetra_SerialSymDenseMatrix &A_in)
Sets the pointers for coefficient matrix; special version for symmetric matrices. ...
Epetra_SerialDenseMatrix: A class for constructing and using real double precision general dense matr...
Epetra_SerialSymDenseMatrix: A class for constructing and using symmetric positive definite dense mat...
int UnequilibrateLHS(void)
Unscales the solution vectors if equilibration was used to solve the system.
void POCON(const char UPLO, 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 positive definite matrix (SPOCON)
Epetra_SerialDenseMatrix * RHS_
void CopyUPLOMat(bool Upper, double *A, int LDA, int NumRows)
double * A() const
Returns pointer to the this matrix.
int ApplyRefinement(void)
Apply Iterative Refinement.
int Solve(void)
Computes the solution X to AX = B for the this matrix and the B provided to SetVectors()..
bool ShouldEquilibrate()
Returns true if the LAPACK general rules for equilibration suggest you should equilibrate the system...
int ReciprocalConditionEstimate(double &Value)
Returns the reciprocal of the 1-norm condition number of the this matrix.
void POEQU(const int N, const float *A, const int LDA, float *S, float *SCOND, float *AMAX, int *INFO) const
Epetra_LAPACK equilibration for positive definite matrix (SPOEQU)
Epetra_SerialSymDenseMatrix * SymMatrix_
char UPLO() const
Returns character value of UPLO used by LAPACK routines.
bool Inverted()
Returns true if matrix inverse has been computed (inverse available via AF() and LDAF()).
bool ReciprocalConditionEstimated()
Returns true if the condition number of the this matrix has been computed (value available via Recipr...
int ComputeEquilibrateScaling(void)
Computes the scaling vector S(i) = 1/sqrt(A(i,i) of the this matrix.
bool Upper() const
Returns true if upper triangle of this matrix has and will be used.
int SetMatrix(Epetra_SerialDenseMatrix &A)
Sets the pointers for coefficient matrix.
int LDA() const
Returns the leading dimension of the this matrix.
void POTRF(const char UPLO, const int N, float *A, const int LDA, int *INFO) const
Epetra_LAPACK factorization for positive definite matrix (SPOTRF)
void PORFS(const char UPLO, const int N, const int NRHS, const float *A, const int LDA, const float *AF, const int LDAF, const float *B, const int LDB, float *X, const int LDX, float *FERR, float *BERR, float *WORK, int *IWORK, int *INFO) const
Epetra_LAPACK solve driver for positive definite matrix (SPOSVX)
Epetra_SerialDenseMatrix * Factor_
void POTRS(const char UPLO, const int N, const int NRHS, const float *A, const int LDA, float *X, const int LDX, int *INFO) const
Epetra_LAPACK solve (after factorization) for positive definite matrix (SPOTRS)
const double Epetra_Overflow
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.
int EquilibrateRHS(void)
Equilibrates the current RHS.
Epetra_SerialDenseMatrix * LHS_
int Factor(void)
Computes the in-place Cholesky factorization of the matrix using the LAPACK routine DPOTRF...
void POTRI(const char UPLO, const int N, float *A, const int LDA, int *INFO) const
Epetra_LAPACK inversion for positive definite matrix (SPOTRI)
Epetra_SerialDenseSolver: A class for solving dense linear problems.
void UpdateFlops(int Flops_in) const
Increment Flop count for this object.
int EquilibrateMatrix(void)
Equilibrates the this matrix.