IFPACK  Development
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
60 extern "C" {int dfill_diag(int n, NCformat *Astore);}
61 
62 using Teuchos::RefCountPtr;
63 using Teuchos::rcp;
64 
65 #ifdef IFPACK_TEUCHOS_TIME_MONITOR
66 # include "Teuchos_TimeMonitor.hpp"
67 #endif
68 
69 //==============================================================================
70 Ifpack_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 //==============================================================================
101 void 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 //==========================================================================
124 int 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 //==========================================================================
138 template<typename int_type>
139 int 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 
230 int 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 //==========================================================================
248 int 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
315 int 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
346 int 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
358 int 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 //=============================================================================
392 double 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 //=============================================================================
410 std::ostream&
411 Ifpack_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
Epetra_Time::ResetStartTime
void ResetStartTime(void)
Ifpack_SILU::Initialize
int Initialize()
Initialize the preconditioner, does not touch matrix values.
Definition: Ifpack_SILU.cpp:230
Ifpack_SILU::ApplyInverseTime
virtual double ApplyInverseTime() const
Returns the time spent in ApplyInverse().
Definition: Ifpack_SILU.h:255
Ifpack_SILU::ComputeFlops
virtual double ComputeFlops() const
Returns the number of flops in the computation phase.
Definition: Ifpack_SILU.h:265
Ifpack_SILU::ComputeTime
virtual double ComputeTime() const
Returns the time spent in Compute().
Definition: Ifpack_SILU.h:249
Ifpack_SILU::ApplyInverseFlops
virtual double ApplyInverseFlops() const
Returns the number of flops in the application of the preconditioner.
Definition: Ifpack_SILU.h:270
Ifpack_SILU::SetParameters
int SetParameters(Teuchos::ParameterList &parameterlist)
Set parameters using a Teuchos::ParameterList object.
Definition: Ifpack_SILU.cpp:124
Ifpack_SILU::NumInitialize
virtual int NumInitialize() const
Returns the number of calls to Initialize().
Definition: Ifpack_SILU.h:225
Ifpack_SILU::Comm
const Epetra_Comm & Comm() const
Returns the Epetra_BlockMap object associated with the range of this matrix operator.
Definition: Ifpack_SILU.h:213
Ifpack_SILU::Label
const char * Label() const
Returns a character string describing the operator.
Definition: Ifpack_SILU.h:188
Epetra_CrsMatrix
Epetra_MultiVector::Pointers
double ** Pointers() const
Ifpack_SILU::Compute
int Compute()
Compute ILU factors L and U using the specified graph, diagonal perturbation thresholds and relaxatio...
Definition: Ifpack_SILU.cpp:248
Epetra_CrsMatrix::ColMap
const Epetra_Map & ColMap() const
Ifpack_SILU::Print
virtual std::ostream & Print(std::ostream &os) const
Prints on stream basic information about this object.
Definition: Ifpack_SILU.cpp:411
Ifpack_SILU::NumApplyInverse
virtual int NumApplyInverse() const
Returns the number of calls to ApplyInverse().
Definition: Ifpack_SILU.h:237
Epetra_RowMatrix
Epetra_CrsMatrix::IndicesAreContiguous
bool IndicesAreContiguous() const
Epetra_BlockMap::SameAs
bool SameAs(const Epetra_BlockMap &Map) const
Ifpack_SILU::Ifpack_SILU
Ifpack_SILU(Epetra_RowMatrix *A)
Constructor.
Definition: Ifpack_SILU.cpp:70
Ifpack_SILU::ApplyInverse
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.
Definition: Ifpack_SILU.cpp:358
Epetra_MultiVector
Ifpack_SILU::IsInitialized
bool IsInitialized() const
Returns true if the preconditioner has been successfully initialized.
Definition: Ifpack_SILU.h:105
Ifpack_SILU::Condest
double Condest() const
Returns the computed estimated condition number, or -1.0 if not computed.
Definition: Ifpack_SILU.h:177
Epetra_Time::ElapsedTime
double ElapsedTime(void) const
Ifpack_SILU::NumCompute
virtual int NumCompute() const
Returns the number of calls to Compute().
Definition: Ifpack_SILU.h:231
Epetra_CrsMatrix::RowMap
const Epetra_Map & RowMap() const
Epetra_CrsMatrix::ExtractCrsDataPointers
int ExtractCrsDataPointers(int *&IndexOffset, int *&Indices, double *&Values_in) const
Ifpack_SILU::UseTranspose
bool UseTranspose() const
Returns the current UseTranspose setting.
Definition: Ifpack_SILU.h:204
Ifpack_SILU::InitializeTime
virtual double InitializeTime() const
Returns the time spent in Initialize().
Definition: Ifpack_SILU.h:243
Epetra_MultiVector::NumVectors
int NumVectors() const
Ifpack_SILU::IsComputed
bool IsComputed() const
If factor is completed, this query returns true, otherwise it returns false.
Definition: Ifpack_SILU.h:114
Generated by   doxygen 1.8.16