Ifpack Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
Ifpack_SILU.cpp
Go to the documentation of this file.
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//==============================================================================
70Ifpack_SILU::Ifpack_SILU(Epetra_RowMatrix* Matrix_in) :
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
230int Ifpack_SILU::Initialize() {
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//==========================================================================
248int Ifpack_SILU::Compute()
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
358int Ifpack_SILU::ApplyInverse(const Epetra_MultiVector& X,
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
const int N
Ifpack_CondestType
Ifpack_CondestType: enum to define the type of condition number estimate.
double Ifpack_Condest(const Ifpack_Preconditioner &IFP, const Ifpack_CondestType CT, const int MaxIters, const double Tol, Epetra_RowMatrix *Matrix)
#define IFPACK_CHK_ERR(ifpack_err)
int NumVectors() const
double ** Pointers() const
#define false