IFPACK Development
Loading...
Searching...
No Matches
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
65extern "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//==============================================================================
82Ifpack_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//==============================================================================
110Ifpack_SPARSKIT::~Ifpack_SPARSKIT()
111{
112}
113
114//==========================================================================
115int 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//==========================================================================
133int Ifpack_SPARSKIT::Initialize()
134{
135 IsInitialized_ = true;
136 IsComputed_ = false;
137 return(0);
138}
139
140//==========================================================================
141int 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
266int 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//=============================================================================
297double 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//=============================================================================
311std::ostream&
312Ifpack_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
int NumVectors() const