IFPACK Development
Loading...
Searching...
No Matches
Ifpack_SILU.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_SILU.h"
45#ifdef HAVE_IFPACK_SUPERLU
46
47#include "Ifpack_CondestType.h"
48#include "Epetra_ConfigDefs.h"
49#include "Epetra_Comm.h"
50#include "Epetra_Map.h"
51#include "Epetra_RowMatrix.h"
52#include "Epetra_Vector.h"
53#include "Epetra_MultiVector.h"
54#include "Epetra_CrsGraph.h"
55#include "Epetra_CrsMatrix.h"
56#include "Teuchos_ParameterList.hpp"
57#include "Teuchos_RefCountPtr.hpp"
58
59// SuperLU includes
60extern "C" {int dfill_diag(int n, NCformat *Astore);}
61
62using Teuchos::RefCountPtr;
63using Teuchos::rcp;
64
65#ifdef IFPACK_TEUCHOS_TIME_MONITOR
66# include "Teuchos_TimeMonitor.hpp"
67#endif
68
69//==============================================================================
71 A_(rcp(Matrix_in,false)),
72 Comm_(Matrix_in->Comm()),
73 UseTranspose_(false),
74 NumMyDiagonals_(0),
75 DropTol_(1e-4),
76 FillTol_(1e-2),
77 FillFactor_(10.0),
78 DropRule_(9),
79 Condest_(-1.0),
80 IsInitialized_(false),
81 IsComputed_(false),
82 NumInitialize_(0),
83 NumCompute_(0),
84 NumApplyInverse_(0),
85 InitializeTime_(0.0),
86 ComputeTime_(0.0),
87 ApplyInverseTime_(0.0),
88 Time_(Comm()),
89 etree_(0),
90 perm_r_(0),
91 perm_c_(0)
92{
93 Teuchos::ParameterList List;
94 SetParameters(List);
95 SY_.Store=0;
96}
97
98
99
100//==============================================================================
101void Ifpack_SILU::Destroy()
102{
103 if(IsInitialized_){
104 // SuperLU cleanup
105 StatFree(&stat_);
106
107 Destroy_CompCol_Permuted(&SAc_);
108 Destroy_SuperNode_Matrix(&SL_);
109 Destroy_CompCol_Matrix(&SU_);
110
111 // Make sure NOT to clean up Epetra's memory
112 Destroy_SuperMatrix_Store(&SA_);
113 if(SY_.Store) Destroy_SuperMatrix_Store(&SY_);
114 SY_.Store=0;
115
116 // Cleanup stuff I allocated
117 delete [] etree_;etree_=0;
118 delete [] perm_r_;perm_r_=0;
119 delete [] perm_c_;perm_c_=0;
120 }
121}
122
123//==========================================================================
124int Ifpack_SILU::SetParameters(Teuchos::ParameterList& List)
125{
126 DropTol_=List.get("fact: drop tolerance",DropTol_);
127 FillTol_=List.get("fact: zero pivot threshold",FillTol_);
128 FillFactor_=List.get("fact: maximum fill factor",FillFactor_);
129 DropRule_=List.get("fact: silu drop rule",DropRule_);
130
131 // set label
132 sprintf(Label_, "IFPACK SILU (drop=%d, zpv=%f, ffact=%f, rthr=%f)",
133 DropRule(),FillTol(),FillFactor(),DropTol());
134 return(0);
135}
136
137//==========================================================================
138template<typename int_type>
139int Ifpack_SILU::TInitialize()
140{
141
142#ifdef IFPACK_TEUCHOS_TIME_MONITOR
143 TEUCHOS_FUNC_TIME_MONITOR("Ifpack_SILU::Initialize");
144#endif
145
146 Time_.ResetStartTime();
147
148 // reset this object
149 Destroy();
150
151 IsInitialized_ = false;
152
153 Epetra_CrsMatrix* CrsMatrix;
154 CrsMatrix = dynamic_cast<Epetra_CrsMatrix*>(&*A_);
155
156 if(CrsMatrix && CrsMatrix->RowMap().SameAs(CrsMatrix->ColMap()) && CrsMatrix->IndicesAreContiguous()){
157 // Case #1: Use original matrix
158 Aover_=rcp(CrsMatrix,false);
159 }
160 else if(CrsMatrix && CrsMatrix->IndicesAreContiguous()){
161 // Case #2: Extract using CrsDataPointers w/ contiguous indices
162 int size = A_->MaxNumEntries();
163 int N=A_->NumMyRows();
164 Aover_ = rcp(new Epetra_CrsMatrix(Copy,A_->RowMatrixRowMap(), size));
165 std::vector<int_type> Indices(size);
166 std::vector<double> Values(size);
167
168 int i,j,ct,*rowptr,*colind;
169 double *values;
170 IFPACK_CHK_ERR(CrsMatrix->ExtractCrsDataPointers(rowptr,colind,values));
171
172 // Use the fact that EpetraCrsMatrices always number the off-processor columns *LAST*
173 for(i=0;i<N;i++){
174 for(j=rowptr[i],ct=0;j<rowptr[i+1];j++){
175 if(colind[j]<N){
176 Indices[ct]= (int_type) CrsMatrix->GCID64(colind[j]);
177 Values[ct]=values[j];
178 ct++;
179 }
180 }
181 Aover_->InsertGlobalValues((int_type) CrsMatrix->GRID64(i),ct,&Values[0],&Indices[0]);
182 }
183 IFPACK_CHK_ERR(Aover_->FillComplete(CrsMatrix->RowMap(),CrsMatrix->RowMap()));
184 }
185 else{
186 // Case #3: Extract using copys
187 int size = A_->MaxNumEntries();
188 Aover_ = rcp(new Epetra_CrsMatrix(Copy,A_->RowMatrixRowMap(), size));
189 if (Aover_.get() == 0) IFPACK_CHK_ERR(-5); // memory allocation error
190
191 std::vector<int> Indices1(size);
192 std::vector<int_type> Indices2(size);
193 std::vector<double> Values1(size),Values2(size);
194
195 // extract each row at-a-time, and insert it into
196 // the graph, ignore all off-process entries
197 int N=A_->NumMyRows();
198 for (int i = 0 ; i < N ; ++i) {
199 int NumEntries;
200 int_type GlobalRow = (int_type) A_->RowMatrixRowMap().GID64(i);
201 IFPACK_CHK_ERR(A_->ExtractMyRowCopy(i, size, NumEntries,
202 &Values1[0], &Indices1[0]));
203
204 // convert to global indices, keeping only on-proc entries
205 int ct=0;
206 for (int j=0; j < NumEntries ; ++j) {
207 if(Indices1[j] < N){
208 Indices2[ct] = (int_type) A_->RowMatrixColMap().GID64(Indices1[j]);
209 Values2[ct]=Values1[j];
210 ct++;
211 }
212 }
213 IFPACK_CHK_ERR(Aover_->InsertGlobalValues(GlobalRow,ct,
214 &Values2[0],&Indices2[0]));
215 }
216 IFPACK_CHK_ERR(Aover_->FillComplete(A_->RowMatrixRowMap(),A_->RowMatrixRowMap()));
217 }
218
219 // Finishing touches
220 Aover_->OptimizeStorage();
221 Graph_=rcp(const_cast<Epetra_CrsGraph*>(&Aover_->Graph()),false);
222
223 IsInitialized_ = true;
224 NumInitialize_++;
225 InitializeTime_ += Time_.ElapsedTime();
226
227 return(0);
228}
229
231#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
232 if(A_->RowMatrixRowMap().GlobalIndicesInt()) {
233 return TInitialize<int>();
234 }
235 else
236#endif
237#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
238 if(A_->RowMatrixRowMap().GlobalIndicesLongLong()) {
239 return TInitialize<long long>();
240 }
241 else
242#endif
243 throw "Ifpack_SILU::Initialize: GlobalIndices type unknown for A_";
244}
245
246
247//==========================================================================
249{
250
251#ifdef IFPACK_TEUCHOS_TIME_MONITOR
252 TEUCHOS_FUNC_TIME_MONITOR("Ifpack_SILU::Compute");
253#endif
254
255 if (!IsInitialized())
256 IFPACK_CHK_ERR(Initialize());
257
258 Time_.ResetStartTime();
259 IsComputed_ = false;
260
261 // Initialize the SuperLU statistics & options variables.
262 StatInit(&stat_);
263 ilu_set_default_options(&options_);
264 options_.ILU_DropTol=DropTol_;
265 options_.ILU_FillTol=FillTol_;
266 options_.ILU_DropRule=DropRule_;
267 options_.ILU_FillFactor=FillFactor_;
268
269 // Grab pointers from Aover_
270 int *rowptr,*colind;
271 double *values;
272 IFPACK_CHK_ERR(Aover_->ExtractCrsDataPointers(rowptr,colind,values));
273 int N=Aover_->NumMyRows();
274
275 // Copy the data over to SuperLU land - mark as a transposed CompCol Matrix
276 dCreate_CompCol_Matrix(&SA_,N,N,Aover_->NumMyNonzeros(),
277 values,colind,rowptr,SLU_NC,SLU_D,SLU_GE);
278
279 // Fill any zeros on the diagonal
280 // Commented out for now
281 dfill_diag(N, (NCformat*)SA_.Store);
282
283 // Allocate SLU memory
284 etree_=new int [N];
285 perm_r_=new int[N];
286 perm_c_=new int[N];
287
288 // Get column permutation
289 int permc_spec=options_.ColPerm;
290 if ( permc_spec != MY_PERMC && options_.Fact == DOFACT )
291 get_perm_c(permc_spec,&SA_,perm_c_);
292
293 // Preorder by column permutation
294 sp_preorder(&options_, &SA_, perm_c_, etree_, &SAc_);
295
296 // Call the factorization
297 int panel_size = sp_ienv(1);
298 int relax = sp_ienv(2);
299 int info=0;
300 dgsitrf(&options_,&SAc_,relax,panel_size,etree_,NULL,0,perm_c_,perm_r_,&SL_,&SU_,
301#ifdef HAVE_IFPACK_SUPERLU5_API
302 &lu_,
303#endif
304 &stat_,&info);
305 if(info<0) IFPACK_CHK_ERR(info);
306
307 IsComputed_ = true;
308 NumCompute_++;
309 ComputeTime_ += Time_.ElapsedTime();
310 return 0;
311}
312
313//=============================================================================
314// This function finds Y such that LDU Y = X or U(trans) D L(trans) Y = X for multiple RHS
315int Ifpack_SILU::Solve(bool Trans, const Epetra_MultiVector& X,
316 Epetra_MultiVector& Y) const
317{
318
319#ifdef IFPACK_TEUCHOS_TIME_MONITOR
320 TEUCHOS_FUNC_TIME_MONITOR("Ifpack_SILU::ApplyInverse - Solve");
321#endif
322
323 if (!IsComputed())
324 IFPACK_CHK_ERR(-3);
325 int nrhs=X.NumVectors();
326 int N=A_->NumMyRows();
327
328 // Copy X over to Y
329 Y=X;
330
331 // Move to SuperLU land
332 // NTS: Need to do epetra-style realloc-if-nrhs-changes thing
333 if(SY_.Store) Destroy_SuperMatrix_Store(&SY_);
334 SY_.Store=0;
335 dCreate_Dense_Matrix(&SY_,N,nrhs,Y[0],N,SLU_DN,SLU_D,SLU_GE);
336
337 int info;
338 dgstrs(TRANS,&SL_,&SU_,perm_c_,perm_r_,&SY_,&stat_,&info);
339 if(!info) IFPACK_CHK_ERR(info);
340
341 return(info);
342}
343
344//=============================================================================
345// This function finds X such that LDU Y = X or U(trans) D L(trans) Y = X for multiple RHS
346int Ifpack_SILU::Multiply(bool Trans, const Epetra_MultiVector& X,
347 Epetra_MultiVector& Y) const
348{
349
350 if (!IsComputed())
351 IFPACK_CHK_ERR(-1);
352
353 return(0);
354}
355
356//=============================================================================
357// This function finds X such that LDU Y = X or U(trans) D L(trans) Y = X for multiple RHS
359 Epetra_MultiVector& Y) const
360{
361
362#ifdef IFPACK_TEUCHOS_TIME_MONITOR
363 TEUCHOS_FUNC_TIME_MONITOR("Ifpack_SILU::ApplyInverse");
364#endif
365
366 if (!IsComputed())
367 IFPACK_CHK_ERR(-3);
368
369 if (X.NumVectors() != Y.NumVectors())
370 IFPACK_CHK_ERR(-2);
371
372 Time_.ResetStartTime();
373
374 // AztecOO gives X and Y pointing to the same memory location,
375 // need to create an auxiliary vector, Xcopy
376 Teuchos::RefCountPtr< const Epetra_MultiVector > Xcopy;
377 if (X.Pointers()[0] == Y.Pointers()[0])
378 Xcopy = Teuchos::rcp( new Epetra_MultiVector(X) );
379 else
380 Xcopy = Teuchos::rcp( &X, false );
381
382 IFPACK_CHK_ERR(Solve(Ifpack_SILU::UseTranspose(), *Xcopy, Y));
383
384 ++NumApplyInverse_;
385 ApplyInverseTime_ += Time_.ElapsedTime();
386
387 return(0);
388
389}
390
391//=============================================================================
392double Ifpack_SILU::Condest(const Ifpack_CondestType CT,
393 const int MaxIters, const double Tol,
394 Epetra_RowMatrix* Matrix_in)
395{
396
397#ifdef IFPACK_TEUCHOS_TIME_MONITOR
398 TEUCHOS_FUNC_TIME_MONITOR("Ifpack_SILU::Condest");
399#endif
400
401 if (!IsComputed()) // cannot compute right now
402 return(-1.0);
403
404 Condest_ = Ifpack_Condest(*this, CT, MaxIters, Tol, Matrix_in);
405
406 return(Condest_);
407}
408
409//=============================================================================
410std::ostream&
411Ifpack_SILU::Print(std::ostream& os) const
412{
413 using std::endl;
414
415 if (!Comm().MyPID()) {
416 os << endl;
417 os << "================================================================================" << endl;
418 os << "Ifpack_SILU: " << Label() << endl << endl;
419 os << "Dropping rule = "<< DropRule() << endl;
420 os << "Zero pivot thresh = "<< FillTol() << endl;
421 os << "Max fill factor = "<< FillFactor() << endl;
422 os << "Drop tolerance = "<< DropTol() << endl;
423 os << "Condition number estimate = " << Condest() << endl;
424 os << "Global number of rows = " << A_->NumGlobalRows64() << endl;
425 if (IsComputed_) {
426 // Internal SuperLU info
427 int fnnz=0;
428 if(SL_.Store) fnnz+=((SCformat*)SL_.Store)->nnz;
429 if(SU_.Store) fnnz+=((NCformat*)SU_.Store)->nnz;
430 int lufill=fnnz - A_->NumMyRows();
431 os << "No. of nonzeros in L+U = "<< lufill<<endl;
432 os << "Fill ratio: nnz(F)/nnz(A) = "<<(double)lufill / (double)A_->NumMyNonzeros()<<endl;
433 }
434 os << endl;
435 os << "Phase # calls Total Time (s) Total MFlops MFlops/s" << endl;
436 os << "----- ------- -------------- ------------ --------" << endl;
437 os << "Initialize() " << std::setw(5) << NumInitialize()
438 << " " << std::setw(15) << InitializeTime()
439 << " 0.0 0.0" << endl;
440 os << "Compute() " << std::setw(5) << NumCompute()
441 << " " << std::setw(15) << ComputeTime()
442 << " " << std::setw(15) << 1.0e-6 * ComputeFlops();
443 if (ComputeTime() != 0.0)
444 os << " " << std::setw(15) << 1.0e-6 * ComputeFlops() / ComputeTime() << endl;
445 else
446 os << " " << std::setw(15) << 0.0 << endl;
447 os << "ApplyInverse() " << std::setw(5) << NumApplyInverse()
448 << " " << std::setw(15) << ApplyInverseTime()
449 << " " << std::setw(15) << 1.0e-6 * ApplyInverseFlops();
450 if (ApplyInverseTime() != 0.0)
451 os << " " << std::setw(15) << 1.0e-6 * ApplyInverseFlops() / ApplyInverseTime() << endl;
452 else
453 os << " " << std::setw(15) << 0.0 << endl;
454 os << "================================================================================" << endl;
455 os << endl;
456 }
457
458 return(os);
459}
460
461#endif
bool SameAs(const Epetra_BlockMap &Map) const
const Epetra_Map & RowMap() const
int ExtractCrsDataPointers(int *&IndexOffset, int *&Indices, double *&Values_in) const
const Epetra_Map & ColMap() const
bool IndicesAreContiguous() const
int NumVectors() const
double ** Pointers() const
void ResetStartTime(void)
double ElapsedTime(void) const
double Condest() const
Returns the computed estimated condition number, or -1.0 if not computed.
const Epetra_Comm & Comm() const
Returns the Epetra_BlockMap object associated with the range of this matrix operator.
virtual int NumApplyInverse() const
Returns the number of calls to ApplyInverse().
virtual double ComputeTime() const
Returns the time spent in Compute().
virtual int NumCompute() const
Returns the number of calls to Compute().
int Initialize()
Initialize the preconditioner, does not touch matrix values.
Ifpack_SILU(Epetra_RowMatrix *A)
Constructor.
bool IsInitialized() const
Returns true if the preconditioner has been successfully initialized.
virtual int NumInitialize() const
Returns the number of calls to Initialize().
virtual double ComputeFlops() const
Returns the number of flops in the computation phase.
bool UseTranspose() const
Returns the current UseTranspose setting.
int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Returns the result of a Epetra_Operator inverse applied to an Epetra_MultiVector X in Y.
virtual double InitializeTime() const
Returns the time spent in Initialize().
virtual std::ostream & Print(std::ostream &os) const
Prints on stream basic information about this object.
int SetParameters(Teuchos::ParameterList &parameterlist)
Set parameters using a Teuchos::ParameterList object.
virtual double ApplyInverseTime() const
Returns the time spent in ApplyInverse().
int Compute()
Compute ILU factors L and U using the specified graph, diagonal perturbation thresholds and relaxatio...
bool IsComputed() const
If factor is completed, this query returns true, otherwise it returns false.
virtual double ApplyInverseFlops() const
Returns the number of flops in the application of the preconditioner.
const char * Label() const
Returns a character string describing the operator.