← Back to branch summary

~njansson/dolfin/hpc

« back to all changes in this revision

Viewing changes to src/kernel/la/PETScKrylovSolver.cpp

  • Committer: Johannes Ring
  • Date: 2008-03-05 22:43:06 UTC
  • Revision ID: johannr@simula.no-20080305224306-2npsdyhfdpl2esji
The BIG commit!

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/kernel/la/PETScKrylovSolver.cpp
1
 
// Copyright (C) 2005 Johan Jansson.
2
 
// Licensed under the GNU LGPL Version 2.1.
3
 
//
4
 
// Modified by Anders Logg 2005-2006.
5
 
// Modified by Garth N. Wells 2005-2006.
6
 
//
7
 
// First added:  2005-12-02
8
 
// Last changed: 2007-07-31
9
 
 
10
 
#ifdef HAVE_PETSC_H
11
 
 
12
 
#include <private/pcimpl.h>
13
 
 
14
 
#include <dolfin/dolfin_log.h>
15
 
#include <dolfin/PETScKrylovSolver.h>
16
 
 
17
 
#include <dolfin/PETScMatrix.h>
18
 
#include <dolfin/PETScVector.h>
19
 
#include <dolfin/PETScKrylovMatrix.h>
20
 
 
21
 
using namespace dolfin;
22
 
 
23
 
// Monitor function
24
 
namespace dolfin
25
 
{
26
 
  int monitor(KSP ksp, int iteration, real rnorm, void *mctx)
27
 
  {
28
 
    message("Iteration %d: residual = %g", iteration, rnorm);
29
 
    return 0;
30
 
  }
31
 
}
32
 
 
33
 
//-----------------------------------------------------------------------------
34
 
PETScKrylovSolver::PETScKrylovSolver(KrylovMethod method, Preconditioner pc)
35
 
  : PETScLinearSolver(),
36
 
    method(method), pc_petsc(pc), pc_dolfin(0),
37
 
    ksp(0), M(0), N(0), parameters_read(false), pc_set(false)
38
 
{
39
 
  // Do nothing
40
 
}
41
 
//-----------------------------------------------------------------------------
42
 
PETScKrylovSolver::PETScKrylovSolver(KrylovMethod method,
43
 
                                     PETScPreconditioner& preconditioner)
44
 
  : PETScLinearSolver(),
45
 
    method(method), pc_petsc(default_pc), pc_dolfin(&preconditioner),
46
 
    ksp(0), M(0), N(0), parameters_read(false), pc_set(false)
47
 
{
48
 
  // Do nothing
49
 
}
50
 
//-----------------------------------------------------------------------------
51
 
PETScKrylovSolver::~PETScKrylovSolver()
52
 
{
53
 
  // Destroy solver environment.
54
 
  if ( ksp ) KSPDestroy(ksp);
55
 
}
56
 
//-----------------------------------------------------------------------------
57
 
dolfin::uint PETScKrylovSolver::solve(const PETScMatrix& A, PETScVector& x, const PETScVector& b)
58
 
{
59
 
  // Check dimensions
60
 
  uint M = A.size(0);
61
 
  uint N = A.size(1);
62
 
  if ( N != b.size() )
63
 
    error("Non-matching dimensions for linear system.");
64
 
 
65
 
  // Write a message
66
 
  if ( get("Krylov report") )
67
 
    message("Solving linear system of size %d x %d (Krylov solver).", M, N);
68
 
 
69
 
  // Reinitialize KSP solver if necessary
70
 
  init(M, N);
71
 
 
72
 
  // Reinitialize solution vector if necessary
73
 
  x.init(M);
74
 
 
75
 
  // Read parameters if not done
76
 
  if ( !parameters_read )
77
 
    readParameters();
78
 
 
79
 
  // Solve linear system
80
 
  KSPSetOperators(ksp, A.mat(), A.mat(), SAME_NONZERO_PATTERN);
81
 
 
82
 
  // FIXME: Preconditioner being set here to avoid PETSc bug with Hypre.
83
 
  //        See explanation inside PETScKrylovSolver:init().
84
 
  if( !pc_set )
85
 
  { 
86
 
    setPETScPreconditioner();
87
 
    pc_set = true;   
88
 
  }
89
 
 
90
 
  KSPSolve(ksp, b.vec(), x.vec());
91
 
 
92
 
  // Check if the solution converged
93
 
  KSPConvergedReason reason;
94
 
  KSPGetConvergedReason(ksp, &reason);
95
 
  if ( reason < 0 )
96
 
    error("Krylov solver did not converge.");
97
 
 
98
 
  // Get the number of iterations
99
 
  int num_iterations = 0;
100
 
  KSPGetIterationNumber(ksp, &num_iterations);
101
 
 
102
 
  // Report results
103
 
  writeReport(num_iterations);
104
 
 
105
 
  return num_iterations;
106
 
}
107
 
//-----------------------------------------------------------------------------
108
 
dolfin::uint PETScKrylovSolver::solve(const PETScKrylovMatrix& A, PETScVector& x, const PETScVector& b)
109
 
{
110
 
  // Check dimensions
111
 
  uint M = A.size(0);
112
 
  uint N = A.size(1);
113
 
  if ( N != b.size() )
114
 
    error("Non-matching dimensions for linear system.");
115
 
  
116
 
  // Write a message
117
 
  if ( get("Krylov report") )
118
 
    message("Solving virtual linear system of size %d x %d (Krylov solver).", M, N);
119
 
 
120
 
  // Reinitialize KSP solver if necessary
121
 
  init(M, N);
122
 
 
123
 
  // Reinitialize solution vector if necessary
124
 
  x.init(M);
125
 
 
126
 
  // Read parameters if not done
127
 
  if ( !parameters_read )
128
 
    readParameters();
129
 
 
130
 
  // Don't use preconditioner that can't handle virtual (shell) matrix
131
 
  if ( !pc_dolfin )
132
 
  {
133
 
    PC pc;
134
 
    KSPGetPC(ksp, &pc);
135
 
    PCSetType(pc, PCNONE);
136
 
  }
137
 
 
138
 
  // Solve linear system
139
 
  KSPSetOperators(ksp, A.mat(), A.mat(), DIFFERENT_NONZERO_PATTERN);
140
 
  KSPSolve(ksp, b.vec(), x.vec());  
141
 
 
142
 
  // Check if the solution converged
143
 
  KSPConvergedReason reason;
144
 
  KSPGetConvergedReason(ksp, &reason);
145
 
  if ( reason < 0 )
146
 
    error("Krylov solver did not converge.");
147
 
 
148
 
  // Get the number of iterations
149
 
  int num_iterations = 0;
150
 
  KSPGetIterationNumber(ksp, &num_iterations);
151
 
  
152
 
  // Report results
153
 
  writeReport(num_iterations);
154
 
 
155
 
  return num_iterations;
156
 
}
157
 
//-----------------------------------------------------------------------------
158
 
void PETScKrylovSolver::disp() const
159
 
{
160
 
  KSPView(ksp, PETSC_VIEWER_STDOUT_WORLD);
161
 
}
162
 
//-----------------------------------------------------------------------------
163
 
void PETScKrylovSolver::init(uint M, uint N)
164
 
{
165
 
  // Check if we need to reinitialize
166
 
  if ( ksp != 0 && M == this->M && N == this->N )
167
 
    return;
168
 
 
169
 
  // Save size of system
170
 
  this->M = M;
171
 
  this->N = N;
172
 
 
173
 
  // Destroy old solver environment if necessary
174
 
  if ( ksp != 0 )
175
 
    KSPDestroy(ksp);
176
 
 
177
 
  // Set up solver environment
178
 
  KSPCreate(PETSC_COMM_SELF, &ksp);
179
 
  KSPSetFromOptions(ksp);  
180
 
  //KSPSetInitialGuessNonzero(ksp, PETSC_TRUE);
181
 
 
182
 
  // Set solver
183
 
  setSolver();
184
 
 
185
 
  // FIXME: The preconditioner is being set in solve() due to a PETSc bug
186
 
  //        when using Hypre preconditioner. The problem can be avoided by
187
 
  //        setting the preconditioner after KSPSetOperators(). This will be
188
 
  //        fixed in PETSc, the the preconditioner can be set here again.
189
 
  // Set preconditioner
190
 
//  setPETScPreconditioner();
191
 
}
192
 
//-----------------------------------------------------------------------------
193
 
void PETScKrylovSolver::readParameters()
194
 
{
195
 
  // Don't do anything if not initialized
196
 
  if ( !ksp )
197
 
    return;
198
 
 
199
 
  // Set monitor
200
 
  if ( get("Krylov monitor convergence") )
201
 
    KSPMonitorSet(ksp, monitor, 0, 0);
202
 
 
203
 
  // Set tolerances
204
 
  KSPSetTolerances(ksp,
205
 
                   get("Krylov relative tolerance"),
206
 
                   get("Krylov absolute tolerance"),
207
 
                   get("Krylov divergence limit"),
208
 
                   get("Krylov maximum iterations"));
209
 
 
210
 
  // Set nonzero shift for preconditioner
211
 
  if ( !pc_dolfin )
212
 
  {
213
 
    PC pc;
214
 
    KSPGetPC(ksp, &pc);
215
 
    PCFactorSetShiftNonzero(pc, get("Krylov shift nonzero"));
216
 
  }
217
 
 
218
 
  // Remember that we have read parameters
219
 
  parameters_read = true;
220
 
}
221
 
//-----------------------------------------------------------------------------
222
 
void PETScKrylovSolver::setSolver()
223
 
{
224
 
  // Don't do anything for default method
225
 
  if ( method == default_method )
226
 
    return;
227
 
  
228
 
  // Set PETSc Krylov solver
229
 
  KSPType ksp_type = getType(method);
230
 
  KSPSetType(ksp, ksp_type);
231
 
}
232
 
//-----------------------------------------------------------------------------
233
 
void PETScKrylovSolver::setPETScPreconditioner()
234
 
{
235
 
  // Treat special case DOLFIN user-defined preconditioner
236
 
  if ( pc_dolfin )
237
 
  {
238
 
    PETScPreconditioner::setup(ksp, *pc_dolfin);
239
 
    return;
240
 
  }
241
 
 
242
 
  // Treat special case default preconditioner (do nothing)
243
 
  if ( pc_petsc == default_pc )
244
 
    return;
245
 
 
246
 
  // Get PETSc PC pointer
247
 
  PC pc;
248
 
  KSPGetPC(ksp, &pc);
249
 
 
250
 
  // Treat special case Hypre AMG preconditioner
251
 
  if ( pc_petsc == amg )
252
 
  {  
253
 
#if PETSC_HAVE_HYPRE
254
 
    PCSetType(pc, PCHYPRE);
255
 
    PCHYPRESetType(pc, "boomeramg");
256
 
#else
257
 
    warning("PETSc has not been compiled with the HYPRE library for   "
258
 
                   "algerbraic multigrid. Default PETSc solver will be used. "
259
 
                   "For performance, installation of HYPRE is recommended.   "
260
 
                   "See the DOLFIN user manual for more information.");
261
 
#endif
262
 
    return;
263
 
  }
264
 
 
265
 
  // Set preconditioner
266
 
  PCSetType(pc, PETScPreconditioner::getType(pc_petsc));
267
 
}
268
 
//-----------------------------------------------------------------------------
269
 
void PETScKrylovSolver::writeReport(int num_iterations)
270
 
{
271
 
  // Check if we should write the report
272
 
  bool report = get("Krylov report");
273
 
  if ( !report )
274
 
    return;
275
 
    
276
 
  // Get name of solver
277
 
  KSPType ksp_type;
278
 
  KSPGetType(ksp, &ksp_type);
279
 
 
280
 
  // Get name of preconditioner
281
 
  PC pc;
282
 
  KSPGetPC(ksp, &pc);
283
 
  PCType pc_type;
284
 
  PCGetType(pc, &pc_type);
285
 
 
286
 
  // Report number of iterations and solver type
287
 
  message("Krylov solver (%s, %s) converged in %d iterations.",
288
 
              ksp_type, pc_type, num_iterations);
289
 
}
290
 
//-----------------------------------------------------------------------------
291
 
KSPType PETScKrylovSolver::getType(KrylovMethod method) const
292
 
{
293
 
  switch (method)
294
 
  {
295
 
  case bicgstab:
296
 
    return KSPBCGS;
297
 
  case cg:
298
 
    return KSPCG;
299
 
  case default_method:
300
 
    return "default";
301
 
  case gmres:
302
 
    return KSPGMRES;
303
 
  default:
304
 
    warning("Requested Krylov method unknown. Using GMRES.");
305
 
    return KSPGMRES;
306
 
  }
307
 
}
308
 
//-----------------------------------------------------------------------------
309
 
 
310
 
#endif