Epetra Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
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
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();
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();
322 Value = RCOND_;
323 UpdateFlops(2*N_*N_); // Not sure of count
325 return(0);
326}
const double Epetra_Overflow
#define EPETRA_CHK_ERR(a)
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)
void UpdateFlops(int Flops_in) const
Increment Flop count for this object.
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)
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)
void POTRI(const char UPLO, const int N, float *A, const int LDA, int *INFO) const
Epetra_LAPACK inversion for positive definite matrix (SPOTRI)
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)
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 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_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.
Epetra_SerialDenseSolver: A class for solving dense linear problems.
bool ReciprocalConditionEstimated()
Returns true if the condition number of the this matrix has been computed (value available via Recipr...
int EquilibrateRHS(void)
Equilibrates the current RHS.
Epetra_SerialDenseMatrix * Factor_
int SetMatrix(Epetra_SerialDenseMatrix &A)
Sets the pointers for coefficient matrix.
Epetra_SerialDenseMatrix * RHS_
int UnequilibrateLHS(void)
Unscales the solution vectors if equilibration was used to solve the system.
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_SerialDenseMatrix * LHS_
bool Solved()
Returns true if the current set of vectors has been solved.
int Invert(void)
Inverts the this matrix.
int Solve(void)
Computes the solution X to AX = B for the this matrix and the B provided to SetVectors()....
int Factor(void)
Computes the in-place Cholesky factorization of the matrix using the LAPACK routine DPOTRF.
bool ShouldEquilibrate()
Returns true if the LAPACK general rules for equilibration suggest you should equilibrate the system.
Epetra_SerialSymDenseMatrix * SymFactor_
int SetMatrix(Epetra_SerialSymDenseMatrix &A_in)
Sets the pointers for coefficient matrix; special version for symmetric matrices.
Epetra_SerialSpdDenseSolver()
Default constructor; matrix should be set using SetMatrix(), LHS and RHS set with SetVectors().
int EquilibrateMatrix(void)
Equilibrates the this matrix.
virtual ~Epetra_SerialSpdDenseSolver()
Epetra_SerialDenseSolver destructor.
int ComputeEquilibrateScaling(void)
Computes the scaling vector S(i) = 1/sqrt(A(i,i) of the this matrix.
int ReciprocalConditionEstimate(double &Value)
Returns the reciprocal of the 1-norm condition number of the this matrix.
Epetra_SerialSymDenseMatrix * SymMatrix_
int ApplyRefinement(void)
Apply Iterative Refinement.
Epetra_SerialSymDenseMatrix: A class for constructing and using symmetric positive definite dense mat...
void CopyUPLOMat(bool Upper, double *A, int LDA, int NumRows)
char UPLO() const
Returns character value of UPLO used by LAPACK routines.
bool Upper() const
Returns true if upper triangle of this matrix has and will be used.
double OneNorm() const
Computes the 1-Norm of the this matrix (identical to NormOne() method).