Ifpack Package Browser (Single Doxygen Collection)  Development
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 
44 #include "Epetra_SerialDenseVector.h"
45 #include <iostream>
46 
47 //=============================================================================
50  Epetra_BLAS(),
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 {
88  InitPointers();
89  ResetMatrix();
90  ResetVectors();
91 }
92 //=============================================================================
94 {
95  DeleteArrays();
96 }
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;
170  SolutionErrorsEstimated_ = 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();
335  AllocateIWORK();
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 //=============================================================================
345 void 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 }
virtual void Print(std::ostream &os) const
Print service methods; defines behavior of ostream << operator.
virtual int Invert(void)
Inverts the this matrix.
Epetra_SerialDenseMatrix * RHS_
void UpdateFlops(int Flops_in) const
Ifpack_SerialTriDiSolver()
Default constructor; matrix should be set using SetMatrix(), LHS and RHS set with SetVectors()...
virtual int Solve(void)
Computes the solution X to AX = B for the this matrix and the B provided to SetVectors()..
virtual int ApplyRefinement(void)
Apply Iterative Refinement.
virtual int Factor(void)
Computes the in-place LU factorization of the matrix using the LAPACK routine DGETRF.
int SetMatrix(Ifpack_SerialTriDiMatrix &A)
Sets the pointers for coefficient matrix.
int N() const
Returns column dimension of system.
virtual ~Ifpack_SerialTriDiSolver()
Ifpack_SerialTriDiSolver destructor.
bool Inverted()
Returns true if matrix inverse has been computed (inverse available via AF() and LDAF()).
#define EPETRA_CHK_ERR(a)
void EstimateSolutionErrors(bool Flag)
Causes all solves to estimate the forward and backward solution error.
bool Factored()
Returns true if matrix is factored (factor available via AF() and LDAF()).
double * A() const
Returns pointer to the this matrix.
double * DL()
Returns pointer to the this matrix.
Ifpack_SerialTriDiMatrix * Matrix_
virtual int ReciprocalConditionEstimate(double &Value)
Unscales the solution vectors if equilibration was used to solve the system.
double * A() const
bool ReciprocalConditionEstimated()
Returns true if the condition number of the this matrix has been computed (value available via Recipr...
#define false
Teuchos::LAPACK< int, double > lapack
virtual double OneNorm() const
Computes the 1-Norm of the this matrix (identical to NormOne() method).
Ifpack_SerialTriDiMatrix: A class for constructing and using real double precision general TriDi matr...
Ifpack_SerialTriDiMatrix * Factor_
int SetVectors(Epetra_SerialDenseMatrix &X, Epetra_SerialDenseMatrix &B)
Sets the pointers for left and right hand side vector(s).
Epetra_SerialDenseMatrix * LHS_