Ifpack Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
Ifpack_SerialTriDiSolver.cpp
Go to the documentation of this file.
1//@HEADER
2// ************************************************************************
3//
4// Epetra: Linear Algebra Services Package
5// Copyright 2011 Sandia Corporation
6//
7// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8// the U.S. Government retains certain rights in this software.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38//
39// ************************************************************************
40//@HEADER
41
45#include <iostream>
46
47//=============================================================================
51 Transpose_(false),
52 Factored_(false),
53 EstimateSolutionErrors_(false),
54 SolutionErrorsEstimated_(false),
55 Solved_(false),
56 Inverted_(false),
57 ReciprocalConditionEstimated_(false),
58 RefineSolution_(false),
59 SolutionRefined_(false),
60 TRANS_('N'),
61 N_(0),
62 NRHS_(0),
63 LDA_(0),
64 LDAF_(0),
65 LDB_(0),
66 LDX_(0),
67 INFO_(0),
68 LWORK_(0),
69 IPIV_(0),
70 IWORK_(0),
71 ANORM_(0.0),
72 RCOND_(0.0),
73 ROWCND_(0.0),
74 COLCND_(0.0),
75 AMAX_(0.0),
76 Matrix_(0),
77 LHS_(0),
78 RHS_(0),
79 Factor_(0),
80 A_(0),
81 FERR_(0),
82 BERR_(0),
83 AF_(0),
84 WORK_(0),
85 B_(0),
86 X_(0)
87{
91}
92//=============================================================================
97//=============================================================================
99{
100 IWORK_ = 0;
101 FERR_ = 0;
102 BERR_ = 0;
103 Factor_ =0;
104 Matrix_ =0;
105 AF_ = 0;
106 IPIV_ = 0;
107 WORK_ = 0;
108 INFO_ = 0;
109 LWORK_ = 0;
110}
111//=============================================================================
113{
114 if (IWORK_ != 0) {delete [] IWORK_;IWORK_ = 0;}
115 if (FERR_ != 0) {delete [] FERR_; FERR_ = 0;}
116 if (BERR_ != 0) {delete [] BERR_; BERR_ = 0;}
117 if (Factor_ != Matrix_ && Factor_ != 0) {delete Factor_; Factor_ = 0;}
118 if (Factor_ !=0) Factor_ = 0;
119
120 if (IPIV_ != 0) {delete [] IPIV_;IPIV_ = 0;}
121 if (WORK_ != 0) {delete [] WORK_;WORK_ = 0;}
122
123 if (AF_ !=0) AF_ = 0;
124
125 INFO_ = 0;
126 LWORK_ = 0;
127}
128//=============================================================================
130{
131 DeleteArrays();
132 ResetVectors();
133 Matrix_ = 0;
134 Factor_ = 0;
135 Factored_ = false;
136 Inverted_ = false;
137 N_ = 0;
138 LDA_ = 0;
139 LDAF_ = 0;
140 ANORM_ = -1.0;
141 RCOND_ = -1.0;
142 ROWCND_ = -1.0;
143 COLCND_ = -1.0;
144 AMAX_ = -1.0;
145 A_ = 0;
146
147}
148//=============================================================================
150 ResetMatrix();
151 Matrix_ = &A_in;
152 Factor_ = &A_in;
153 N_ = A_in.N();
154 A_ = A_in.A();
155 LDA_ = A_in.LDA();
156 LDAF_ = LDA_;
157 AF_ = A_in.A();
158 return(0);
159}
160//=============================================================================
162{
163 LHS_ = 0;
164 RHS_ = 0;
165 B_ = 0;
166 X_ = 0;
168 SolutionRefined_ = false;
169 Solved_ = false;
171 NRHS_ = 0;
172 LDB_ = 0;
173 LDX_ = 0;
174}
175//=============================================================================
177{
178 if (B_in.N() != X_in.N()) EPETRA_CHK_ERR(-1);
179 if (B_in.A()==0) EPETRA_CHK_ERR(-2);
180 if (X_in.A()==0) EPETRA_CHK_ERR(-4);
181
182 ResetVectors();
183 LHS_ = &X_in;
184 RHS_ = &B_in;
185 NRHS_ = B_in.N();
186
187 B_ = B_in.A();
188 X_ = X_in.A();
189 return(0);
190}
191//=============================================================================
194 // If the errors are estimated, this implies that the solution must be refined
196 return;
197}
198//=============================================================================
200 if (Factored()) return(0); // Already factored
201 if (Inverted()) EPETRA_CHK_ERR(-100); // Cannot factor inverted matrix
202
203 ANORM_ = Matrix_->OneNorm(); // Compute 1-Norm of A
204
205 // If we want to refine the solution, then the factor must
206 // be stored separatedly from the original matrix
207
209
210 if (A_ == AF_)
211 if (RefineSolution_ ) {
213 F = Factor_;
214 AF_ = Factor_->A();
215 LDAF_ = Factor_->LDA();
216 }
217
218 if (IPIV_==0) IPIV_ = new int[N_]; // Allocated Pivot vector if not already done.
219
220 double * DL_ = F->DL();
221 double * D_ = F->D();
222 double * DU_ = F->DU();
223 double * DU2_ = F->DU2();
224
225 lapack.GTTRF(N_, DL_, D_, DU_, DU2_, IPIV_, &INFO_);
226
227 Factored_ = true;
228 double DN = N_;
229 UpdateFlops( (N_ == 1)? 1. : 4*(DN-1) );
230
232 return(0);
233
234}
235
236//=============================================================================
238 int ierr = 0;
239
240 // We will call one of four routines depending on what services the user wants and
241 // whether or not the matrix has been inverted or factored already.
242 //
243 // If the matrix has been inverted, use DGEMM to compute solution.
244 // Otherwise, if the user want the matrix to be equilibrated or wants a refined solution, we will
245 // call the X interface.
246 // Otherwise, if the matrix is already factored we will call the TRS interface.
247 // Otherwise, if the matrix is unfactored we will call the SV interface.
248
249 if (B_==0) EPETRA_CHK_ERR(-3); // No B
250 if (X_==0) EPETRA_CHK_ERR(-4); // No X
251
252 double DN = N_;
253 double DNRHS = NRHS_;
254 if (Inverted()) {
255
256 EPETRA_CHK_ERR(-101); // don't allow this \cbl
257
258 }
259 else {
260
261 if (!Factored()) Factor(); // Matrix must be factored
262
263 if (B_!=X_) {
264 *LHS_ = *RHS_; // Copy B to X if needed
265 X_ = LHS_->A();
266 }
267
269 if(A_ == AF_)
270 F = Matrix_;
271 else
272 F = Factor_;
273
274 double * DL_ = F->DL();
275 double * D_ = F->D();
276 double * DU_ = F->DU();
277 double * DU2_ = F->DU2();
278
279 lapack.GTTRS(TRANS_,N_,NRHS_,DL_,D_,DU_,DU2_,IPIV_,X_,N_,&INFO_);
280
281 if (INFO_!=0) EPETRA_CHK_ERR(INFO_);
282 UpdateFlops(2.0*4*DN*DNRHS);
283 Solved_ = true;
284
285 }
286 int ierr1=0;
287 if (RefineSolution_ && !Inverted()) ierr1 = ApplyRefinement();
288 if (ierr1!=0) EPETRA_CHK_ERR(ierr1)
289 else
290 EPETRA_CHK_ERR(ierr);
291
292 return(0);
293}
294//=============================================================================
296 {
297 std::cout<<" SerialTriDiSolver::ApplyRefinement this function is not supported"<<std::endl;
298 EPETRA_CHK_ERR(-102);
299 return(0);
300 }
301
302//=============================================================================
304{
305 if (!Factored()) Factor(); // Need matrix factored.
306
307 // Setting LWORK = -1 and calling GETRI will return optimal work space size in
308 AllocateWORK();
309
310 lapack.GETRI ( N_, AF_, LDAF_, IPIV_, WORK_, LWORK_, &INFO_);
311
312 double DN = N_;
313 UpdateFlops((DN*DN*DN));
314 Inverted_ = true;
315 Factored_ = false;
316
318 return(0);
319}
320
321//=============================================================================
323{
324 int ierr = 0;
326 Value = RCOND_;
327 return(0); // Already computed, just return it.
328 }
329
330 if (ANORM_<0.0) ANORM_ = Matrix_->OneNorm();
331 if (!Factored()) ierr = Factor(); // Need matrix factored.
332 if (ierr!=0) EPETRA_CHK_ERR(ierr-2);
333
334 AllocateWORK();
336 // We will assume a one-norm condition number \\ works for TriDi
337 lapack.GECON( '1', N_, AF_, LDAF_, ANORM_, &RCOND_, WORK_, IWORK_, &INFO_);
339 Value = RCOND_;
340 UpdateFlops(2*N_*N_); // Not sure of count
342 return(0);
343}
344//=============================================================================
345void Ifpack_SerialTriDiSolver::Print(std::ostream& os) const {
346
347 if (Matrix_!=0) os << "Solver Matrix" << std::endl << *Matrix_ << std::endl;
348 if (Factor_!=0) os << "Solver Factored Matrix" << std::endl << *Factor_ << std::endl;
349 if (LHS_ !=0) os << "Solver LHS" << std::endl << *LHS_ << std::endl;
350 if (RHS_ !=0) os << "Solver RHS" << std::endl << *RHS_ << std::endl;
351
352}
#define EPETRA_CHK_ERR(a)
void UpdateFlops(int Flops_in) const
double * A() const
Ifpack_SerialTriDiMatrix: A class for constructing and using real double precision general TriDi matr...
virtual double OneNorm() const
Computes the 1-Norm of the this matrix (identical to NormOne() method).
double * A() const
Returns pointer to the this matrix.
double * DL()
Returns pointer to the this matrix.
bool ReciprocalConditionEstimated()
Returns true if the condition number of the this matrix has been computed (value available via Recipr...
bool Inverted()
Returns true if matrix inverse has been computed (inverse available via AF() and LDAF()).
void EstimateSolutionErrors(bool Flag)
Causes all solves to estimate the forward and backward solution error.
virtual int ApplyRefinement(void)
Apply Iterative Refinement.
Epetra_SerialDenseMatrix * LHS_
virtual void Print(std::ostream &os) const
Print service methods; defines behavior of ostream << operator.
bool Factored()
Returns true if matrix is factored (factor available via AF() and LDAF()).
virtual int Invert(void)
Inverts the this matrix.
int SetMatrix(Ifpack_SerialTriDiMatrix &A)
Sets the pointers for coefficient matrix.
virtual int Solve(void)
Computes the solution X to AX = B for the this matrix and the B provided to SetVectors()....
virtual int Factor(void)
Computes the in-place LU factorization of the matrix using the LAPACK routine DGETRF.
Ifpack_SerialTriDiMatrix * Matrix_
Epetra_SerialDenseMatrix * RHS_
Teuchos::LAPACK< int, double > lapack
virtual int ReciprocalConditionEstimate(double &Value)
Unscales the solution vectors if equilibration was used to solve the system.
int SetVectors(Epetra_SerialDenseMatrix &X, Epetra_SerialDenseMatrix &B)
Sets the pointers for left and right hand side vector(s).
Ifpack_SerialTriDiSolver()
Default constructor; matrix should be set using SetMatrix(), LHS and RHS set with SetVectors().
virtual ~Ifpack_SerialTriDiSolver()
Ifpack_SerialTriDiSolver destructor.
Ifpack_SerialTriDiMatrix * Factor_
#define false