IFPACK Development
Loading...
Searching...
No Matches
Ifpack_IKLU.cpp
1/*@HEADER
2// ***********************************************************************
3//
4// Ifpack: Object-Oriented Algebraic Preconditioner Package
5// Copyright (2002) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
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*/
42
43#include "Ifpack_ConfigDefs.h"
44#include "Ifpack_Preconditioner.h"
45#include "Ifpack_IKLU.h"
46#include "Ifpack_Condest.h"
47#include "Ifpack_Utils.h"
48#include "Ifpack_HashTable.h"
49#include "Epetra_SerialComm.h"
50#include "Epetra_Comm.h"
51#include "Epetra_Map.h"
52#include "Epetra_RowMatrix.h"
53#include "Epetra_CrsMatrix.h"
54#include "Epetra_Vector.h"
55#include "Epetra_MultiVector.h"
56#include "Epetra_Util.h"
57#include "Teuchos_ParameterList.hpp"
58#include "Teuchos_RefCountPtr.hpp"
59#include <functional>
60#include <algorithm>
61
62using namespace Teuchos;
63
64//==============================================================================
65// FIXME: allocate Comm_ and Time_ the first Initialize() call
67 A_(*A),
68 Comm_(A->Comm()),
69 Condest_(-1.0),
70 Relax_(0.),
71 Athresh_(0.0),
72 Rthresh_(1.0),
73 LevelOfFill_(1.0),
74 DropTolerance_(1e-12),
75 IsInitialized_(false),
76 IsComputed_(false),
77 UseTranspose_(false),
78 NumMyRows_(-1),
79 NumInitialize_(0),
80 NumCompute_(0),
81 NumApplyInverse_(0),
82 InitializeTime_(0.0),
83 ComputeTime_(0.0),
84 ApplyInverseTime_(0.0),
85 ComputeFlops_(0.0),
86 ApplyInverseFlops_(0.0),
87 Time_(Comm()),
88 GlobalNonzeros_(0),
89 csrA_(0),
90 cssS_(0),
91 csrnN_(0)
92{
93 // do nothing here..
94}
95
96//==============================================================================
98{
99 Destroy();
100}
101
102//==============================================================================
103void Ifpack_IKLU::Destroy()
104{
105 IsInitialized_ = false;
106 IsComputed_ = false;
107 if (csrA_)
108 csr_spfree( csrA_ );
109 if (cssS_)
110 csr_sfree( cssS_ );
111 if (csrnN_)
112 csr_nfree( csrnN_ );
113}
114
115//==========================================================================
116int Ifpack_IKLU::SetParameters(Teuchos::ParameterList& List)
117{
118 using std::cerr;
119 using std::endl;
120
121 try
122 {
123 LevelOfFill_ = List.get<double>("fact: ilut level-of-fill", LevelOfFill());
124 if (LevelOfFill_ <= 0.0)
125 IFPACK_CHK_ERR(-2); // must be greater than 0.0
126
127 Athresh_ = List.get<double>("fact: absolute threshold", Athresh_);
128 Rthresh_ = List.get<double>("fact: relative threshold", Rthresh_);
129 Relax_ = List.get<double>("fact: relax value", Relax_);
130 DropTolerance_ = List.get<double>("fact: drop tolerance", DropTolerance_);
131
132 Label_ = "IFPACK IKLU (fill=" + Ifpack_toString(LevelOfFill())
133 + ", relax=" + Ifpack_toString(RelaxValue())
134 + ", athr=" + Ifpack_toString(AbsoluteThreshold())
135 + ", rthr=" + Ifpack_toString(RelativeThreshold())
136 + ", droptol=" + Ifpack_toString(DropTolerance())
137 + ")";
138 return(0);
139 }
140 catch (...)
141 {
142 cerr << "Caught an exception while parsing the parameter list" << endl;
143 cerr << "This typically means that a parameter was set with the" << endl;
144 cerr << "wrong type (for example, int instead of double). " << endl;
145 cerr << "please check the documentation for the type required by each parameer." << endl;
146 IFPACK_CHK_ERR(-1);
147 }
148
149 //return(0); // unreachable
150}
151
152//==========================================================================
154{
155 using std::cerr;
156 using std::cout;
157 using std::endl;
158
159 // delete previously allocated factorization
160 Destroy();
161
162 Time_.ResetStartTime();
163
164 if (A_.Comm().NumProc() != 1) {
165 cout << " There are too many processors !!! " << endl;
166 cerr << "Ifpack_IKLU can be used with Comm().NumProc() == 1" << endl;
167 cerr << "only. This class is a subdomain solver for Ifpack_AdditiveSchwarz," << endl;
168 cerr << "and it is currently not meant to be used otherwise." << endl;
169 exit(EXIT_FAILURE);
170 }
171
172 // check dimensions of input matrix only in serial
173 if (Comm().NumProc() == 1 && Matrix().NumMyRows() != Matrix().NumMyCols())
174 IFPACK_CHK_ERR(-2);
175
176 NumMyRows_ = Matrix().NumMyRows();
177 NumMyNonzeros_ = Matrix().NumMyNonzeros();
178
179 int RowNnz, Length = Matrix().MaxNumEntries();
180 std::vector<int> RowIndices(Length);
181 std::vector<double> RowValues(Length);
182
183 //cout << "Processor " << Comm().MyPID() << " owns " << NumMyRows_ << " rows and has " << NumMyNonzeros_ << " nonzeros " << endl;
184 // get general symbolic structure of the matrix
185 csrA_ = csr_spalloc( NumMyRows_, NumMyRows_, NumMyNonzeros_, 1, 0 );
186
187 // copy the symbolic structure into csrA_
188 int count = 0;
189 csrA_->p[0] = 0;
190 for (int i = 0; i < NumMyRows_; ++i ) {
191
192 IFPACK_CHK_ERR(A_.ExtractMyRowCopy(i,Length,RowNnz,
193 &RowValues[0],&RowIndices[0]));
194 for (int j = 0 ; j < RowNnz ; ++j) {
195 csrA_->j[count++] = RowIndices[j];
196 //cout << "Row = " << i << ", Column = " << RowIndices[j] << ", Value = " << RowValues[j] << endl;
197 }
198 csrA_->p[i+1] = csrA_->p[i] + RowNnz;
199 }
200
201 // Perform symbolic analysis on the current matrix structure
202 int order = 1;
203 cssS_ = csr_sqr( order, csrA_ );
204 for (int i = 0; i < NumMyRows_; ++i ) {
205 cout << "AMD Perm (from inside KLU) [" << i << "] = " << cssS_->q[i] << endl;
206 }
207
208 // nothing else to do here
209 IsInitialized_ = true;
210 ++NumInitialize_;
211 InitializeTime_ += Time_.ElapsedTime();
212
213 return(0);
214}
215
216//==========================================================================
218{
219 public:
220 inline bool operator()(const double& x, const double& y)
221 {
222 return(IFPACK_ABS(x) > IFPACK_ABS(y));
223 }
224};
225
226//==========================================================================
227
229{
230 using std::cout;
231 using std::endl;
232
233 if (!IsInitialized())
234 IFPACK_CHK_ERR(Initialize());
235
236 Time_.ResetStartTime();
237 IsComputed_ = false;
238
239 NumMyRows_ = A_.NumMyRows();
240 int Length = A_.MaxNumEntries();
241
242 bool distributed = (Comm().NumProc() > 1)?true:false;
243#if !defined(EPETRA_NO_32BIT_GLOBAL_INDICES) || !defined(EPETRA_NO_64BIT_GLOBAL_INDICES)
244 if (distributed)
245 {
246 SerialComm_ = rcp(new Epetra_SerialComm);
247 SerialMap_ = rcp(new Epetra_Map(NumMyRows_, 0, *SerialComm_));
248 assert (SerialComm_.get() != 0);
249 assert (SerialMap_.get() != 0);
250 }
251 else
252 SerialMap_ = rcp(const_cast<Epetra_Map*>(&A_.RowMatrixRowMap()), false);
253#endif
254
255 int RowNnz;
256 std::vector<int> RowIndices(Length);
257 std::vector<double> RowValues(Length);
258
259 // copy the values from A_ into csrA_
260 int count = 0;
261 for (int i = 0; i < NumMyRows_; ++i ) {
262
263 IFPACK_CHK_ERR(A_.ExtractMyRowCopy(i,Length,RowNnz,
264 &RowValues[0],&RowIndices[0]));
265 // make sure each row has the same number of nonzeros
266 if (RowNnz != (csrA_->p[i+1]-csrA_->p[i])) {
267 cout << "The number of nonzeros for this row does not math the expected number of nonzeros!!!" << endl;
268 }
269 for (int j = 0 ; j < RowNnz ; ++j) {
270
271 csrA_->x[count++] = RowValues[j];
272 //cout << "Row = " << i << ", Column = " << RowIndices[j] << ", Value = " << RowValues[j] << endl;
273 }
274 }
275
276 // compute the lu factors
277 double tol = 0.1;
278 csrnN_ = csr_lu( &*csrA_, &*cssS_, tol );
279
280 // Create L and U as a view of the information stored in csrnN_->L and csrnN_->U
281 csr* L_tmp = csrnN_->L;
282 csr* U_tmp = csrnN_->U;
283 std::vector<int> numEntriesL( NumMyRows_ ), numEntriesU( NumMyRows_ );
284 for (int i=0; i < NumMyRows_; ++i) {
285 numEntriesL[i] = ( L_tmp->p[i+1] - L_tmp->p[i] );
286 numEntriesU[i] = ( U_tmp->p[i+1] - U_tmp->p[i] );
287 }
288 L_ = rcp(new Epetra_CrsMatrix(View, *SerialMap_, &numEntriesL[0]));
289 U_ = rcp(new Epetra_CrsMatrix(View, *SerialMap_, &numEntriesU[0]));
290
291 // Insert the values into L and U
292#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
293 if(SerialMap_->GlobalIndicesInt()) {
294 for (int i=0; i < NumMyRows_; ++i) {
295 L_->InsertGlobalValues( i, numEntriesL[i], &(L_tmp->x[L_tmp->p[i]]), &(L_tmp->j[L_tmp->p[i]]) );
296 U_->InsertGlobalValues( i, numEntriesU[i], &(U_tmp->x[U_tmp->p[i]]), &(U_tmp->j[U_tmp->p[i]]) );
297 }
298 }
299 else
300#endif
301#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
302 if(SerialMap_->GlobalIndicesLongLong()) {
303
304 const int MaxNumEntries_L_U = std::max(L_->MaxNumEntries(), U_->MaxNumEntries());
305 std::vector<long long> entries(MaxNumEntries_L_U);
306
307 for (int i=0; i < NumMyRows_; ++i) {
308 std::copy(&(L_tmp->j[L_tmp->p[i]]), &(L_tmp->j[L_tmp->p[i]]) + numEntriesL[i], entries.begin());
309 L_->InsertGlobalValues( i, numEntriesL[i], &(L_tmp->x[L_tmp->p[i]]), &(entries[0]) );
310
311 std::copy(&(U_tmp->j[U_tmp->p[i]]), &(U_tmp->j[U_tmp->p[i]]) + numEntriesU[i], entries.begin());
312 U_->InsertGlobalValues( i, numEntriesU[i], &(U_tmp->x[U_tmp->p[i]]), &(entries[0]) );
313 }
314 }
315 else
316#endif
317 throw "Ifpack_IKLU::Compute: GlobalIndices type unknown for SerialMap_";
318
319 IFPACK_CHK_ERR(L_->FillComplete());
320 IFPACK_CHK_ERR(U_->FillComplete());
321
322 long long MyNonzeros = L_->NumGlobalNonzeros64() + U_->NumGlobalNonzeros64();
323 Comm().SumAll(&MyNonzeros, &GlobalNonzeros_, 1);
324
325 IsComputed_ = true;
326
327 ++NumCompute_;
328 ComputeTime_ += Time_.ElapsedTime();
329
330 return(0);
331
332}
333
334//=============================================================================
336 Epetra_MultiVector& Y) const
337{
338 if (!IsComputed())
339 IFPACK_CHK_ERR(-2); // compute preconditioner first
340
341 if (X.NumVectors() != Y.NumVectors())
342 IFPACK_CHK_ERR(-3); // Return error: X and Y not the same size
343
344 Time_.ResetStartTime();
345
346 // NOTE: L_ and U_ are based on SerialMap_, while Xcopy is based
347 // on A.Map()... which are in general different. However, Solve()
348 // does not seem to care... which is fine with me.
349 //
350 // AztecOO gives X and Y pointing to the same memory location,
351 // need to create an auxiliary vector, Xcopy and apply permutation.
352 std::vector<int> invq( NumMyRows_ );
353
354 for (int i=0; i<NumMyRows_; ++i ) {
355 csrnN_->perm[ csrnN_->pinv[i] ] = i;
356 invq[ cssS_->q[i] ] = i;
357 }
358
359 Teuchos::RefCountPtr<Epetra_MultiVector> Xcopy = Teuchos::rcp( new Epetra_MultiVector(X.Map(),X.NumVectors()), false );
360 Teuchos::RefCountPtr<Epetra_MultiVector> Ytemp = Teuchos::rcp( new Epetra_MultiVector(Y.Map(),Y.NumVectors()) );
361
362 for (int i=0; i<NumMyRows_; ++i) {
363 for (int j=0; j<X.NumVectors(); ++j) {
364 Xcopy->ReplaceMyValue( invq[i], j, (*X(j))[i] );
365 }
366 }
367
368 if (!UseTranspose_)
369 {
370 // solves LU Y = X
371 IFPACK_CHK_ERR(L_->Solve(false,false,false,*Xcopy,*Ytemp));
372 IFPACK_CHK_ERR(U_->Solve(true,false,false,*Ytemp,*Ytemp));
373 }
374 else
375 {
376 // solves U(trans) L(trans) Y = X
377 IFPACK_CHK_ERR(U_->Solve(true,true,false,*Xcopy,*Ytemp));
378 IFPACK_CHK_ERR(L_->Solve(false,true,false,*Ytemp,*Ytemp));
379 }
380
381 // Reverse the permutation.
382 for (int i=0; i<NumMyRows_; ++i) {
383 for (int j=0; j<Y.NumVectors(); ++j) {
384 Y.ReplaceMyValue( csrnN_->perm[i], j, (*(*Ytemp)(j))[i] );
385 }
386 }
387
388 ++NumApplyInverse_;
389#ifdef IFPACK_FLOPCOUNTERS
390 ApplyInverseFlops_ += X.NumVectors() * 2 * GlobalNonzeros_;
391#endif
392 ApplyInverseTime_ += Time_.ElapsedTime();
393
394 return(0);
395
396}
397//=============================================================================
398// This function finds X such that LDU Y = X or U(trans) D L(trans) Y = X for multiple RHS
399int Ifpack_IKLU::Apply(const Epetra_MultiVector& /* X */,
400 Epetra_MultiVector& /* Y */) const
401{
402
403 return(-98);
404}
405
406//=============================================================================
407double Ifpack_IKLU::Condest(const Ifpack_CondestType CT,
408 const int MaxIters, const double Tol,
409 Epetra_RowMatrix* Matrix_in)
410{
411 if (!IsComputed()) // cannot compute right now
412 return(-1.0);
413
414 // NOTE: this is computing the *local* condest
415 if (Condest_ == -1.0)
416 Condest_ = Ifpack_Condest(*this, CT, MaxIters, Tol, Matrix_in);
417
418 return(Condest_);
419}
420
421//=============================================================================
422std::ostream&
423Ifpack_IKLU::Print(std::ostream& os) const
424{
425 using std::endl;
426
427 if (!Comm().MyPID()) {
428 os << endl;
429 os << "================================================================================" << endl;
430 os << "Ifpack_IKLU: " << Label() << endl << endl;
431 os << "Level-of-fill = " << LevelOfFill() << endl;
432 os << "Absolute threshold = " << AbsoluteThreshold() << endl;
433 os << "Relative threshold = " << RelativeThreshold() << endl;
434 os << "Relax value = " << RelaxValue() << endl;
435 os << "Condition number estimate = " << Condest() << endl;
436 os << "Global number of rows = " << A_.NumGlobalRows64() << endl;
437 if (IsComputed_) {
438 os << "Number of nonzeros in A = " << A_.NumGlobalNonzeros64() << endl;
439 os << "Number of nonzeros in L + U = " << NumGlobalNonzeros64()
440 << " ( = " << 100.0 * NumGlobalNonzeros64() / A_.NumGlobalNonzeros64()
441 << " % of A)" << endl;
442 os << "nonzeros / rows = "
443 << 1.0 * NumGlobalNonzeros64() / U_->NumGlobalRows64() << endl;
444 }
445 os << endl;
446 os << "Phase # calls Total Time (s) Total MFlops MFlops/s" << endl;
447 os << "----- ------- -------------- ------------ --------" << endl;
448 os << "Initialize() " << std::setw(5) << NumInitialize()
449 << " " << std::setw(15) << InitializeTime()
450 << " 0.0 0.0" << endl;
451 os << "Compute() " << std::setw(5) << NumCompute()
452 << " " << std::setw(15) << ComputeTime()
453 << " " << std::setw(15) << 1.0e-6 * ComputeFlops();
454 if (ComputeTime() != 0.0)
455 os << " " << std::setw(15) << 1.0e-6 * ComputeFlops() / ComputeTime() << endl;
456 else
457 os << " " << std::setw(15) << 0.0 << endl;
458 os << "ApplyInverse() " << std::setw(5) << NumApplyInverse()
459 << " " << std::setw(15) << ApplyInverseTime()
460 << " " << std::setw(15) << 1.0e-6 * ApplyInverseFlops();
461 if (ApplyInverseTime() != 0.0)
462 os << " " << std::setw(15) << 1.0e-6 * ApplyInverseFlops() / ApplyInverseTime() << endl;
463 else
464 os << " " << std::setw(15) << 0.0 << endl;
465 os << "================================================================================" << endl;
466 os << endl;
467 }
468
469 return(os);
470}
View
virtual int NumProc() const=0
virtual int SumAll(double *PartialSums, double *GlobalSums, int Count) const=0
const Epetra_BlockMap & Map() const
int NumVectors() const
int ReplaceMyValue(int MyRow, int VectorIndex, double ScalarValue)
virtual const Epetra_Comm & Comm() const=0
virtual int NumMyRows() const=0
virtual const Epetra_Map & RowMatrixRowMap() const=0
virtual int NumMyNonzeros() const=0
virtual int MaxNumEntries() const=0
virtual int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const=0
void ResetStartTime(void)
double ElapsedTime(void) const
virtual double ApplyInverseTime() const
Returns the time spent in ApplyInverse().
bool IsInitialized() const
Returns true if the preconditioner has been successfully initialized.
virtual int NumCompute() const
Returns the number of calls to Compute().
double DropTolerance() const
Gets the dropping tolerance.
virtual ~Ifpack_IKLU()
Ifpack_IKLU Destructor.
bool IsComputed() const
If factor is completed, this query returns true, otherwise it returns false.
int Compute()
Compute IC factor U using the specified graph, diagonal perturbation thresholds and relaxation parame...
int SetParameters(Teuchos::ParameterList &parameterlis)
Set parameters using a Teuchos::ParameterList object.
const Epetra_RowMatrix & Matrix() const
Returns a reference to the matrix to be preconditioned.
const char * Label() const
Returns the label of this object.
virtual int NumInitialize() const
Returns the number of calls to Initialize().
double AbsoluteThreshold() const
Get absolute threshold value.
Ifpack_IKLU(const Epetra_RowMatrix *A)
Ifpack_IKLU constuctor with variable number of indices per row.
virtual double ComputeFlops() const
Returns the number of flops in the computation phase.
double RelativeThreshold() const
Get relative threshold value.
virtual double InitializeTime() const
Returns the time spent in Initialize().
double Condest() const
Returns the computed estimated condition number, or -1.0 if no computed.
virtual double ComputeTime() const
Returns the time spent in Compute().
int Initialize()
Initialize L and U with values from user matrix A.
double RelaxValue() const
Set relative threshold value.
virtual int NumApplyInverse() const
Returns the number of calls to ApplyInverse().
int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Returns the result of a Ifpack_IKLU forward/back solve on a Epetra_MultiVector X in Y.
const Epetra_Comm & Comm() const
Returns the Epetra_BlockMap object associated with the range of this matrix operator.
virtual std::ostream & Print(std::ostream &os) const
Prints basic information on iostream. This function is used by operator<<.
virtual double ApplyInverseFlops() const
Returns the number of flops in the application of the preconditioner.