Ifpack Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
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 int Invert(void)
Inverts the this matrix.
Epetra_SerialDenseMatrix * RHS_
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.
void GTTRS(const char &TRANS, const OrdinalType &n, const OrdinalType &nrhs, const ScalarType *dl, const ScalarType *d, const ScalarType *du, const ScalarType *du2, const OrdinalType *IPIV, ScalarType *B, const OrdinalType &ldb, OrdinalType *info) const
int SetMatrix(Ifpack_SerialTriDiMatrix &A)
Sets the pointers for coefficient matrix.
void GTTRF(const OrdinalType &n, ScalarType *dl, ScalarType *d, ScalarType *du, ScalarType *du2, OrdinalType *IPIV, OrdinalType *info) const
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()).
virtual double OneNorm() const
Computes the 1-Norm of the this matrix (identical to NormOne() method).
double * A() const
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
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...
#define false
Teuchos::LAPACK< int, double > lapack
int N() const
Returns column dimension of system.
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).
void GETRI(const OrdinalType &n, ScalarType *A, const OrdinalType &lda, const OrdinalType *IPIV, ScalarType *WORK, const OrdinalType &lwork, OrdinalType *info) const
virtual void Print(std::ostream &os) const
Print service methods; defines behavior of ostream &lt;&lt; operator.
void GECON(const char &NORM, const OrdinalType &n, const ScalarType *A, const OrdinalType &lda, const ScalarType &anorm, ScalarType *rcond, ScalarType *WORK, OrdinalType *IWORK, OrdinalType *info) const
Epetra_SerialDenseMatrix * LHS_