IFPACK  Development
 All Classes Namespaces Files Functions Variables Enumerations Friends Pages
Ifpack_SPARSKIT.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 #ifdef HAVE_IFPACK_SPARSKIT
45 #include "Ifpack_Preconditioner.h"
46 #include "Ifpack_SPARSKIT.h"
47 #include "Ifpack_Condest.h"
48 #include "Ifpack_Utils.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_Util.h"
55 #include "Teuchos_ParameterList.hpp"
56 
57 #define F77_ILUT F77_FUNC(ilut, ILUT)
58 #define F77_ILUD F77_FUNC(ilud, ILUD)
59 #define F77_ILUTP F77_FUNC(ilutp, ILUTP)
60 #define F77_ILUDP F77_FUNC(iludp, ILUDP)
61 #define F77_ILUK F77_FUNC(iluk, ILUK)
62 #define F77_ILU0 F77_FUNC(ilu0, ILU0)
63 #define F77_LUSOL F77_FUNC(lusol, LUSOL)
64 
65 extern "C" {
66  void F77_ILUT(int *, double*, int*, int*, int*, double*,
67  double*, int*, int*, int*, double*, int*, int*);
68  void F77_ILUD(int *, double*, int*, int*, double*, double*,
69  double*, int*, int*, int*, double*, int*, int*);
70  void F77_ILUTP(int *, double*, int*, int*, int*, double*, double*, int*,
71  double*, int*, int*, int*, double*, int*, int*, int*);
72  void F77_ILUDP(int *, double*, int*, int*, double*, double*, double*, int*,
73  double*, int*, int*, int*, double*, int*, int*, int*);
74  void F77_ILUK(int *, double*, int*, int*, int*,
75  double*, int*, int*, int*, int*, double*, int*, int*);
76  void F77_ILU0(int*, double*, int*, int*, double*, int*, int*, int*, int*);
77  void F77_LUSOL(int *, double*, double*, double*, int*, int*);
78 }
79 
80 
81 //==============================================================================
82 Ifpack_SPARSKIT::Ifpack_SPARSKIT(Epetra_RowMatrix* A) :
83  A_(*A),
84  Comm_(A->Comm()),
85  UseTranspose_(false),
86  lfil_(0),
87  droptol_(0.0),
88  tol_(0.0),
89  permtol_(0.1),
90  alph_(0.0),
91  mbloc_(-1),
92  Type_("ILUT"),
93  Condest_(-1.0),
94  IsInitialized_(false),
95  IsComputed_(false),
96  NumInitialize_(0),
97  NumCompute_(0),
98  NumApplyInverse_(0),
99  InitializeTime_(0.0),
100  ComputeTime_(0),
101  ApplyInverseTime_(0),
102  ComputeFlops_(0.0),
103  ApplyInverseFlops_(0.0)
104 {
105  Teuchos::ParameterList List;
106  SetParameters(List);
107 }
108 
109 //==============================================================================
110 Ifpack_SPARSKIT::~Ifpack_SPARSKIT()
111 {
112 }
113 
114 //==========================================================================
115 int Ifpack_SPARSKIT::SetParameters(Teuchos::ParameterList& List)
116 {
117  lfil_ = List.get("fact: sparskit: lfil", lfil_);
118  tol_ = List.get("fact: sparskit: tol", tol_);
119  droptol_ = List.get("fact: sparskit: droptol", droptol_);
120  permtol_ = List.get("fact: sparskit: permtol", permtol_);
121  alph_ = List.get("fact: sparskit: alph", alph_);
122  mbloc_ = List.get("fact: sparskit: mbloc", mbloc_);
123  Type_ = List.get("fact: sparskit: type", Type_);
124 
125  // set label
126  Label_ = "IFPACK SPARSKIT (Type=" + Type_ + ", fill=" +
127  Ifpack_toString(lfil_) + ")";
128 
129  return(0);
130 }
131 
132 //==========================================================================
133 int Ifpack_SPARSKIT::Initialize()
134 {
135  IsInitialized_ = true;
136  IsComputed_ = false;
137  return(0);
138 }
139 
140 //==========================================================================
141 int Ifpack_SPARSKIT::Compute()
142 {
143  if (!IsInitialized())
144  IFPACK_CHK_ERR(Initialize());
145 
146  IsComputed_ = false;
147 
148  // convert the matrix into SPARSKIT format. The matrix is then
149  // free'd after method Compute() returns.
150 
151  // convert the matrix into CSR format. Note that nnz is an over-estimate,
152  // since it contains the nonzeros corresponding to external nodes as well.
153  int n = Matrix().NumMyRows();
154  int nnz = Matrix().NumMyNonzeros();
155 
156  std::vector<double> a(nnz);
157  std::vector<int> ja(nnz);
158  std::vector<int> ia(n + 1);
159 
160  const int MaxNumEntries = Matrix().MaxNumEntries();
161 
162  std::vector<double> Values(MaxNumEntries);
163  std::vector<int> Indices(MaxNumEntries);
164 
165  int count = 0;
166 
167  ia[0] = 1;
168  for (int i = 0 ; i < n ; ++i)
169  {
170  int NumEntries;
171  int NumMyEntries = 0;
172  Matrix().ExtractMyRowCopy(i, MaxNumEntries, NumEntries, &Values[0],
173  &Indices[0]);
174 
175  // NOTE: There might be some issues here with the ILU(0) if the column indices aren't sorted.
176  // The other factorizations are probably OK.
177 
178  for (int j = 0 ; j < NumEntries ; ++j)
179  {
180  if (Indices[j] < n) // skip non-local columns
181  {
182  a[count] = Values[j];
183  ja[count] = Indices[j] + 1; // SPARSKIT is FORTRAN
184  ++count;
185  ++NumMyEntries;
186  }
187  }
188  ia[i + 1] = ia[i] + NumMyEntries;
189  }
190 
191  if (mbloc_ == -1) mbloc_ = n;
192 
193  int iwk;
194 
195  if (Type_ == "ILUT" || Type_ == "ILUTP" || Type_ == "ILUD" ||
196  Type_ == "ILUDP")
197  iwk = nnz + 2 * lfil_ * n;
198  else if (Type_ == "ILUK")
199  iwk = (2 * lfil_ + 1) * nnz + 1;
200  else if (Type_ == "ILU0")
201  iwk = nnz+2;
202 
203  int ierr = 0;
204 
205  alu_.resize(iwk);
206  jlu_.resize(iwk);
207  ju_.resize(n + 1);
208 
209  std::vector<int> jw(n + 1);
210  std::vector<double> w(n + 1);
211 
212  if (Type_ == "ILUT")
213  {
214  jw.resize(2 * n);
215  F77_ILUT(&n, &a[0], &ja[0], &ia[0], &lfil_, &droptol_,
216  &alu_[0], &jlu_[0], &ju_[0], &iwk, &w[0], &jw[0], &ierr);
217  }
218  else if (Type_ == "ILUD")
219  {
220  jw.resize(2 * n);
221  F77_ILUD(&n, &a[0], &ja[0], &ia[0], &alph_, &tol_,
222  &alu_[0], &jlu_[0], &ju_[0], &iwk, &w[0], &jw[0], &ierr);
223  }
224  else if (Type_ == "ILUTP")
225  {
226  jw.resize(2 * n);
227  iperm_.resize(2 * n);
228  F77_ILUTP(&n, &a[0], &ja[0], &ia[0], &lfil_, &droptol_, &permtol_,
229  &mbloc_, &alu_[0], &jlu_[0], &ju_[0], &iwk, &w[0], &jw[0],
230  &iperm_[0], &ierr);
231  for (int i = 0 ; i < n ; ++i)
232  iperm_[i]--;
233  }
234  else if (Type_ == "ILUDP")
235  {
236  jw.resize(2 * n);
237  iperm_.resize(2 * n);
238  F77_ILUDP(&n, &a[0], &ja[0], &ia[0], &alph_, &droptol_, &permtol_,
239  &mbloc_, &alu_[0], &jlu_[0], &ju_[0], &n, &w[0], &jw[0],
240  &iperm_[0], &ierr);
241  for (int i = 0 ; i < n ; ++i)
242  iperm_[i]--;
243  }
244  else if (Type_ == "ILUK")
245  {
246  std::vector<int> levs(iwk);
247  jw.resize(3 * n);
248  F77_ILUK(&n, &a[0], &ja[0], &ia[0], &lfil_,
249  &alu_[0], &jlu_[0], &ju_[0], &levs[0], &iwk, &w[0], &jw[0], &ierr);
250  }
251  else if (Type_ == "ILU0")
252  {
253  // here w is only of size n
254  jw.resize(2 * n);
255  F77_ILU0(&n, &a[0], &ja[0], &ia[0],
256  &alu_[0], &jlu_[0], &ju_[0], &jw[0], &ierr);
257  }
258  IFPACK_CHK_ERR(ierr);
259 
260  IsComputed_ = true;
261  return(0);
262 }
263 
264 //=============================================================================
265 // This function finds Y such that LDU Y = X or U(trans) D L(trans) Y = X for multiple RHS
266 int Ifpack_SPARSKIT::ApplyInverse(const Epetra_MultiVector& X,
267  Epetra_MultiVector& Y) const
268 {
269  if (!IsComputed())
270  IFPACK_CHK_ERR(-3); // compute preconditioner first
271 
272  if (X.NumVectors() != Y.NumVectors())
273  IFPACK_CHK_ERR(-2); // Return error: X and Y not the same size
274 
275  int n = Matrix().NumMyRows();
276 
277  for (int i = 0 ; i < X.NumVectors() ; ++i)
278  F77_LUSOL(&n, (double*)X(i)->Values(), Y(i)->Values(), (double*)&alu_[0],
279  (int*)&jlu_[0], (int*)&ju_[0]);
280 
281  // still need to fix support for permutation
282  if (Type_ == "ILUTP" || Type_ == "ILUDP")
283  {
284  std::vector<double> tmp(n);
285  for (int j = 0 ; j < n ; ++j)
286  tmp[iperm_[j]] = Y[0][j];
287  for (int j = 0 ; j < n ; ++j)
288  Y[0][j] = tmp[j];
289  }
290 
291  ++NumApplyInverse_;
292  return(0);
293 
294 }
295 
296 //=============================================================================
297 double Ifpack_SPARSKIT::Condest(const Ifpack_CondestType CT,
298  const int MaxIters, const double Tol,
299  Epetra_RowMatrix* Matrix)
300 {
301  if (!IsComputed()) // cannot compute right now
302  return(-1.0);
303 
304  if (Condest_ == -1.0)
305  Condest_ = Ifpack_Condest(*this, CT, MaxIters, Tol, Matrix);
306 
307  return(Condest_);
308 }
309 
310 //=============================================================================
311 std::ostream&
312 Ifpack_SPARSKIT::Print(std::ostream& os) const
313 {
314  using std::endl;
315  if (!Comm().MyPID()) {
316  os << endl;
317  os << "================================================================================" << endl;
318  os << "Ifpack_SPARSKIT: " << Label() << endl << endl;
319  os << "Condition number estimate = " << Condest() << endl;
320  os << "Global number of rows = " << A_.NumGlobalRows() << endl;
321  os << endl;
322  os << "Phase # calls Total Time (s) Total MFlops MFlops/s" << endl;
323  os << "----- ------- -------------- ------------ --------" << endl;
324  os << "Initialize() " << std::setw(5) << NumInitialize()
325  << " " << std::setw(15) << InitializeTime()
326  << " 0.0 0.0" << endl;
327  os << "Compute() " << std::setw(5) << NumCompute()
328  << " " << std::setw(15) << ComputeTime()
329  << " " << std::setw(15) << 1.0e-6 * ComputeFlops();
330  if (ComputeTime() != 0.0)
331  os << " " << std::setw(15) << 1.0e-6 * ComputeFlops() / ComputeTime() << endl;
332  else
333  os << " " << std::setw(15) << 0.0 << endl;
334  os << "ApplyInverse() " << std::setw(5) << NumApplyInverse()
335  << " " << std::setw(15) << ApplyInverseTime()
336  << " " << std::setw(15) << 1.0e-6 * ApplyInverseFlops();
337  if (ApplyInverseTime() != 0.0)
338  os << " " << std::setw(15) << 1.0e-6 * ApplyInverseFlops() / ApplyInverseTime() << endl;
339  else
340  os << " " << std::setw(15) << 0.0 << endl;
341  os << "================================================================================" << endl;
342  os << endl;
343  }
344 
345 
346  return(os);
347 }
348 
349 #endif // HAVE_IFPACK_SPARSKIT
std::string Ifpack_toString(const int &x)
Converts an integer to std::string.