petsc_linear_solver.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2014 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 #include "libmesh/libmesh_common.h"
19 
20 #ifdef LIBMESH_HAVE_PETSC
21 
22 
23 // C++ includes
24 #include <string.h>
25 
26 // Local Includes
29 #include "libmesh/shell_matrix.h"
30 #include "libmesh/petsc_matrix.h"
32 #include "libmesh/petsc_vector.h"
33 #include "libmesh/string_to_enum.h"
34 
35 namespace libMesh
36 {
37 
38 extern "C"
39 {
40 #if PETSC_VERSION_LESS_THAN(2,2,1)
41  typedef int PetscErrorCode;
42  typedef int PetscInt;
43 #endif
44 
45 
46 #if PETSC_RELEASE_LESS_THAN(3,0,1)
47  PetscErrorCode __libmesh_petsc_preconditioner_setup (void * ctx)
48  {
49  Preconditioner<Number> * preconditioner = static_cast<Preconditioner<Number>*>(ctx);
50 
51  if(!preconditioner->initialized())
52  {
53  err<<"Preconditioner not initialized! Make sure you call init() before solve!"<<std::endl;
54  libmesh_error();
55  }
56 
57  preconditioner->setup();
58 
59  return 0;
60  }
61 
62 
63  PetscErrorCode __libmesh_petsc_preconditioner_apply(void *ctx, Vec x, Vec y)
64  {
65  Preconditioner<Number> * preconditioner = static_cast<Preconditioner<Number>*>(ctx);
66 
67  PetscVector<Number> x_vec(x, preconditioner->comm());
68  PetscVector<Number> y_vec(y, preconditioner->comm());
69 
70  preconditioner->apply(x_vec,y_vec);
71 
72  return 0;
73  }
74 #else
75  PetscErrorCode __libmesh_petsc_preconditioner_setup (PC pc)
76  {
77  void *ctx;
78  PetscErrorCode ierr = PCShellGetContext(pc,&ctx);CHKERRQ(ierr);
79  Preconditioner<Number> * preconditioner = static_cast<Preconditioner<Number>*>(ctx);
80 
81  if(!preconditioner->initialized())
82  {
83  err<<"Preconditioner not initialized! Make sure you call init() before solve!"<<std::endl;
84  libmesh_error();
85  }
86 
87  preconditioner->setup();
88 
89  return 0;
90  }
91 
92  PetscErrorCode __libmesh_petsc_preconditioner_apply(PC pc, Vec x, Vec y)
93  {
94  void *ctx;
95  PetscErrorCode ierr = PCShellGetContext(pc,&ctx);CHKERRQ(ierr);
96  Preconditioner<Number> * preconditioner = static_cast<Preconditioner<Number>*>(ctx);
97 
98  PetscVector<Number> x_vec(x, preconditioner->comm());
99  PetscVector<Number> y_vec(y, preconditioner->comm());
100 
101  preconditioner->apply(x_vec,y_vec);
102 
103  return 0;
104  }
105 #endif
106 } // end extern "C"
107 
108 /*----------------------- functions ----------------------------------*/
109 template <typename T>
111 {
112  if (this->initialized())
113  {
114  /* If we were restricted to some subset, this restriction must
115  be removed and the subset index set destroyed. */
116  if(_restrict_solve_to_is!=NULL)
117  {
118  PetscErrorCode ierr = LibMeshISDestroy(&_restrict_solve_to_is);
119  LIBMESH_CHKERRABORT(ierr);
120  _restrict_solve_to_is = NULL;
121  }
122 
123  if(_restrict_solve_to_is_complement!=NULL)
124  {
125  PetscErrorCode ierr = LibMeshISDestroy(&_restrict_solve_to_is_complement);
126  LIBMESH_CHKERRABORT(ierr);
127  _restrict_solve_to_is_complement = NULL;
128  }
129 
130  this->_is_initialized = false;
131 
133 
134 #if PETSC_VERSION_LESS_THAN(2,2,0)
135 
136  // 2.1.x & earlier style
137  ierr = SLESDestroy(_sles);
138  LIBMESH_CHKERRABORT(ierr);
139 
140 #else
141 
142  // 2.2.0 & newer style
143  ierr = LibMeshKSPDestroy(&_ksp);
144  LIBMESH_CHKERRABORT(ierr);
145 
146 #endif
147 
148 
149  // Mimic PETSc default solver and preconditioner
150  this->_solver_type = GMRES;
151 
152  if(!this->_preconditioner)
153  {
154  if (this->n_processors() == 1)
155  this->_preconditioner_type = ILU_PRECOND;
156  else
157  this->_preconditioner_type = BLOCK_JACOBI_PRECOND;
158  }
159  }
160 }
161 
162 
163 
164 template <typename T>
166 {
167  // Initialize the data structures if not done so already.
168  if (!this->initialized())
169  {
170  this->_is_initialized = true;
171 
173 
174 // 2.1.x & earlier style
175 #if PETSC_VERSION_LESS_THAN(2,2,0)
176 
177  // Create the linear solver context
178  ierr = SLESCreate (this->comm().get(), &_sles);
179  LIBMESH_CHKERRABORT(ierr);
180 
181  // Create the Krylov subspace & preconditioner contexts
182  ierr = SLESGetKSP (_sles, &_ksp);
183  LIBMESH_CHKERRABORT(ierr);
184  ierr = SLESGetPC (_sles, &_pc);
185  LIBMESH_CHKERRABORT(ierr);
186 
187  // Set user-specified solver and preconditioner types
188  this->set_petsc_solver_type();
189 
190  // Set the options from user-input
191  // Set runtime options, e.g.,
192  // -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
193  // These options will override those specified above as long as
194  // SLESSetFromOptions() is called _after_ any other customization
195  // routines.
196 
197  ierr = SLESSetFromOptions (_sles);
198  LIBMESH_CHKERRABORT(ierr);
199 
200 // 2.2.0 & newer style
201 #else
202 
203  // Create the linear solver context
204  ierr = KSPCreate (this->comm().get(), &_ksp);
205  LIBMESH_CHKERRABORT(ierr);
206 
207  // Create the preconditioner context
208  ierr = KSPGetPC (_ksp, &_pc);
209  LIBMESH_CHKERRABORT(ierr);
210 
211  // Set user-specified solver and preconditioner types
212  this->set_petsc_solver_type();
213 
214  // Set the options from user-input
215  // Set runtime options, e.g.,
216  // -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
217  // These options will override those specified above as long as
218  // KSPSetFromOptions() is called _after_ any other customization
219  // routines.
220 
221  ierr = KSPSetFromOptions (_ksp);
222  LIBMESH_CHKERRABORT(ierr);
223 
224  // Not sure if this is necessary, or if it is already handled by KSPSetFromOptions?
225  // NOT NECESSARY!!!!
226  //ierr = PCSetFromOptions (_pc);
227  //LIBMESH_CHKERRABORT(ierr);
228 
229 
230 #endif
231 
232  // Have the Krylov subspace method use our good initial guess
233  // rather than 0, unless the user requested a KSPType of
234  // preonly, which complains if asked to use initial guesses.
235 #if PETSC_VERSION_LESS_THAN(3,0,0) || !PETSC_RELEASE_LESS_THAN(3,4,0)
236  // Pre-3.0 and petsc-dev (as of October 2012) use non-const versions
237  KSPType ksp_type;
238 #else
239  const KSPType ksp_type;
240 #endif
241 
242  ierr = KSPGetType (_ksp, &ksp_type);
243  LIBMESH_CHKERRABORT(ierr);
244 
245  if (strcmp(ksp_type, "preonly"))
246  {
247  ierr = KSPSetInitialGuessNonzero (_ksp, PETSC_TRUE);
248  LIBMESH_CHKERRABORT(ierr);
249  }
250 
251  // Notify PETSc of location to store residual history.
252  // This needs to be called before any solves, since
253  // it sets the residual history length to zero. The default
254  // behavior is for PETSc to allocate (internally) an array
255  // of size 1000 to hold the residual norm history.
256  ierr = KSPSetResidualHistory(_ksp,
257  PETSC_NULL, // pointer to the array which holds the history
258  PETSC_DECIDE, // size of the array holding the history
259  PETSC_TRUE); // Whether or not to reset the history for each solve.
260  LIBMESH_CHKERRABORT(ierr);
261 
262  PetscPreconditioner<T>::set_petsc_preconditioner_type(this->_preconditioner_type,_pc);
263 
264  //If there is a preconditioner object we need to set the internal setup and apply routines
265  if(this->_preconditioner)
266  {
267  this->_preconditioner->init();
268  PCShellSetContext(_pc,(void*)this->_preconditioner);
269  PCShellSetSetUp(_pc,__libmesh_petsc_preconditioner_setup);
270  PCShellSetApply(_pc,__libmesh_petsc_preconditioner_apply);
271  }
272  }
273 }
274 
275 
276 template <typename T>
278 {
279  // Initialize the data structures if not done so already.
280  if (!this->initialized())
281  {
282  this->_is_initialized = true;
283 
285 
286 // 2.1.x & earlier style
287 #if PETSC_VERSION_LESS_THAN(2,2,0)
288 
289  // Create the linear solver context
290  ierr = SLESCreate (this->comm().get(), &_sles);
291  LIBMESH_CHKERRABORT(ierr);
292 
293  // Create the Krylov subspace & preconditioner contexts
294  ierr = SLESGetKSP (_sles, &_ksp);
295  LIBMESH_CHKERRABORT(ierr);
296  ierr = SLESGetPC (_sles, &_pc);
297  LIBMESH_CHKERRABORT(ierr);
298 
299  // Set user-specified solver and preconditioner types
300  this->set_petsc_solver_type();
301 
302  // Set the options from user-input
303  // Set runtime options, e.g.,
304  // -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
305  // These options will override those specified above as long as
306  // SLESSetFromOptions() is called _after_ any other customization
307  // routines.
308 
309  ierr = SLESSetFromOptions (_sles);
310  LIBMESH_CHKERRABORT(ierr);
311 
312 // 2.2.0 & newer style
313 #else
314 
315  // Create the linear solver context
316  ierr = KSPCreate (this->comm().get(), &_ksp);
317  LIBMESH_CHKERRABORT(ierr);
318 
319 
320  //ierr = PCCreate (this->comm().get(), &_pc);
321  // LIBMESH_CHKERRABORT(ierr);
322 
323  // Create the preconditioner context
324  ierr = KSPGetPC (_ksp, &_pc);
325  LIBMESH_CHKERRABORT(ierr);
326 
327  // Set operators. The input matrix works as the preconditioning matrix
328  ierr = KSPSetOperators(_ksp, matrix->mat(), matrix->mat(),DIFFERENT_NONZERO_PATTERN);
329  LIBMESH_CHKERRABORT(ierr);
330 
331  // Set user-specified solver and preconditioner types
332  this->set_petsc_solver_type();
333 
334  // Set the options from user-input
335  // Set runtime options, e.g.,
336  // -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
337  // These options will override those specified above as long as
338  // KSPSetFromOptions() is called _after_ any other customization
339  // routines.
340 
341  ierr = KSPSetFromOptions (_ksp);
342  LIBMESH_CHKERRABORT(ierr);
343 
344  // Not sure if this is necessary, or if it is already handled by KSPSetFromOptions?
345  // NOT NECESSARY!!!!
346  //ierr = PCSetFromOptions (_pc);
347  //LIBMESH_CHKERRABORT(ierr);
348 
349 
350 #endif
351 
352  // Have the Krylov subspace method use our good initial guess
353  // rather than 0, unless the user requested a KSPType of
354  // preonly, which complains if asked to use initial guesses.
355 #if PETSC_VERSION_LESS_THAN(3,0,0) || !PETSC_RELEASE_LESS_THAN(3,4,0)
356  KSPType ksp_type;
357 #else
358  const KSPType ksp_type;
359 #endif
360 
361  ierr = KSPGetType (_ksp, &ksp_type);
362  LIBMESH_CHKERRABORT(ierr);
363 
364  if (strcmp(ksp_type, "preonly"))
365  {
366  ierr = KSPSetInitialGuessNonzero (_ksp, PETSC_TRUE);
367  LIBMESH_CHKERRABORT(ierr);
368  }
369 
370  // Notify PETSc of location to store residual history.
371  // This needs to be called before any solves, since
372  // it sets the residual history length to zero. The default
373  // behavior is for PETSc to allocate (internally) an array
374  // of size 1000 to hold the residual norm history.
375  ierr = KSPSetResidualHistory(_ksp,
376  PETSC_NULL, // pointer to the array which holds the history
377  PETSC_DECIDE, // size of the array holding the history
378  PETSC_TRUE); // Whether or not to reset the history for each solve.
379  LIBMESH_CHKERRABORT(ierr);
380 
381  PetscPreconditioner<T>::set_petsc_preconditioner_type(this->_preconditioner_type,_pc);
382  if(this->_preconditioner)
383  {
384  this->_preconditioner->set_matrix(*matrix);
385  this->_preconditioner->init();
386  PCShellSetContext(_pc,(void*)this->_preconditioner);
387  PCShellSetSetUp(_pc,__libmesh_petsc_preconditioner_setup);
388  PCShellSetApply(_pc,__libmesh_petsc_preconditioner_apply);
389  }
390  }
391 }
392 
393 
394 
395 template <typename T>
396 void
397 PetscLinearSolver<T>::restrict_solve_to (const std::vector<unsigned int>* const dofs,
398  const SubsetSolveMode subset_solve_mode)
399 {
401 
402  /* The preconditioner (in particular if a default preconditioner)
403  will have to be reset. We call this->clear() to do that. This
404  call will also remove and free any previous subset that may have
405  been set before. */
406  this->clear();
407 
408  _subset_solve_mode = subset_solve_mode;
409 
410  if(dofs!=NULL)
411  {
412  PetscInt* petsc_dofs = NULL;
413  ierr = PetscMalloc(dofs->size()*sizeof(PetscInt), &petsc_dofs);
414  LIBMESH_CHKERRABORT(ierr);
415 
416  for(size_t i=0; i<dofs->size(); i++)
417  {
418  petsc_dofs[i] = (*dofs)[i];
419  }
420 
421  ierr = ISCreateLibMesh(this->comm().get(),dofs->size(),petsc_dofs,PETSC_OWN_POINTER,&_restrict_solve_to_is);
422  LIBMESH_CHKERRABORT(ierr);
423  }
424 }
425 
426 
427 
428 template <typename T>
429 std::pair<unsigned int, Real>
431  SparseMatrix<T>& precond_in,
432  NumericVector<T>& solution_in,
433  NumericVector<T>& rhs_in,
434  const double tol,
435  const unsigned int m_its)
436 {
437  START_LOG("solve()", "PetscLinearSolver");
438 
439  // Make sure the data passed in are really of Petsc types
440  PetscMatrix<T>* matrix = libmesh_cast_ptr<PetscMatrix<T>*>(&matrix_in);
441  PetscMatrix<T>* precond = libmesh_cast_ptr<PetscMatrix<T>*>(&precond_in);
442  PetscVector<T>* solution = libmesh_cast_ptr<PetscVector<T>*>(&solution_in);
443  PetscVector<T>* rhs = libmesh_cast_ptr<PetscVector<T>*>(&rhs_in);
444 
445  this->init (matrix);
446 
448  PetscInt its=0, max_its = static_cast<PetscInt>(m_its);
449  PetscReal final_resid=0.;
450 
451  // Close the matrices and vectors in case this wasn't already done.
452  matrix->close ();
453  precond->close ();
454  solution->close ();
455  rhs->close ();
456 
457 // // If matrix != precond, then this means we have specified a
458 // // special preconditioner, so reset preconditioner type to PCMAT.
459 // if (matrix != precond)
460 // {
461 // this->_preconditioner_type = USER_PRECOND;
462 // this->set_petsc_preconditioner_type ();
463 // }
464 
465 // 2.1.x & earlier style
466 #if PETSC_VERSION_LESS_THAN(2,2,0)
467 
468  if(_restrict_solve_to_is!=NULL)
469  {
470  libmesh_not_implemented();
471  }
472 
473  // Set operators. The input matrix works as the preconditioning matrix
474  ierr = SLESSetOperators(_sles, matrix->mat(), precond->mat(),
475  DIFFERENT_NONZERO_PATTERN);
476  LIBMESH_CHKERRABORT(ierr);
477 
478  // Set the tolerances for the iterative solver. Use the user-supplied
479  // tolerance for the relative residual & leave the others at default values.
480  ierr = KSPSetTolerances (_ksp, tol, PETSC_DEFAULT,
481  PETSC_DEFAULT, max_its);
482  LIBMESH_CHKERRABORT(ierr);
483 
484 
485  // Solve the linear system
486  ierr = SLESSolve (_sles, rhs->vec(), solution->vec(), &its);
487  LIBMESH_CHKERRABORT(ierr);
488 
489 
490  // Get the norm of the final residual to return to the user.
491  ierr = KSPGetResidualNorm (_ksp, &final_resid);
492  LIBMESH_CHKERRABORT(ierr);
493 
494 // 2.2.0
495 #elif PETSC_VERSION_LESS_THAN(2,2,1)
496 
497  if(_restrict_solve_to_is!=NULL)
498  {
499  libmesh_not_implemented();
500  }
501 
502  // Set operators. The input matrix works as the preconditioning matrix
503  //ierr = KSPSetOperators(_ksp, matrix->mat(), precond->mat(),
504  // SAME_NONZERO_PATTERN);
505  // LIBMESH_CHKERRABORT(ierr);
506 
507 
508  // Set the tolerances for the iterative solver. Use the user-supplied
509  // tolerance for the relative residual & leave the others at default values.
510  // Convergence is detected at iteration k if
511  // ||r_k||_2 < max(rtol*||b||_2 , abstol)
512  // where r_k is the residual vector and b is the right-hand side. Note that
513  // it is the *maximum* of the two values, the larger of which will almost
514  // always be rtol*||b||_2.
515  ierr = KSPSetTolerances (_ksp,
516  tol, // rtol = relative decrease in residual (1.e-5)
517  PETSC_DEFAULT, // abstol = absolute convergence tolerance (1.e-50)
518  PETSC_DEFAULT, // dtol = divergence tolerance (1.e+5)
519  max_its);
520  LIBMESH_CHKERRABORT(ierr);
521 
522 
523  // Set the solution vector to use
524  ierr = KSPSetSolution (_ksp, solution->vec());
525  LIBMESH_CHKERRABORT(ierr);
526 
527  // Set the RHS vector to use
528  ierr = KSPSetRhs (_ksp, rhs->vec());
529  LIBMESH_CHKERRABORT(ierr);
530 
531  // Solve the linear system
532  ierr = KSPSolve (_ksp);
533  LIBMESH_CHKERRABORT(ierr);
534 
535  // Get the number of iterations required for convergence
536  ierr = KSPGetIterationNumber (_ksp, &its);
537  LIBMESH_CHKERRABORT(ierr);
538 
539  // Get the norm of the final residual to return to the user.
540  ierr = KSPGetResidualNorm (_ksp, &final_resid);
541  LIBMESH_CHKERRABORT(ierr);
542 
543 // 2.2.1 & newer style
544 #else
545 
546  Mat submat = NULL;
547  Mat subprecond = NULL;
548  Vec subrhs = NULL;
549  Vec subsolution = NULL;
550  VecScatter scatter = NULL;
551  PetscMatrix<Number>* subprecond_matrix = NULL;
552 
553  // Set operators. Also restrict rhs and solution vector to
554  // subdomain if neccessary.
555  if(_restrict_solve_to_is!=NULL)
556  {
557  size_t is_local_size = this->_restrict_solve_to_is_local_size();
558 
559  ierr = VecCreate(this->comm().get(),&subrhs);
560  LIBMESH_CHKERRABORT(ierr);
561  ierr = VecSetSizes(subrhs,is_local_size,PETSC_DECIDE);
562  LIBMESH_CHKERRABORT(ierr);
563  ierr = VecSetFromOptions(subrhs);
564  LIBMESH_CHKERRABORT(ierr);
565 
566  ierr = VecCreate(this->comm().get(),&subsolution);
567  LIBMESH_CHKERRABORT(ierr);
568  ierr = VecSetSizes(subsolution,is_local_size,PETSC_DECIDE);
569  LIBMESH_CHKERRABORT(ierr);
570  ierr = VecSetFromOptions(subsolution);
571  LIBMESH_CHKERRABORT(ierr);
572 
573  ierr = VecScatterCreate(rhs->vec(),_restrict_solve_to_is, subrhs,NULL, &scatter);
574  LIBMESH_CHKERRABORT(ierr);
575 
576  ierr = VecScatterBegin(scatter,rhs->vec(),subrhs,INSERT_VALUES,SCATTER_FORWARD);
577  LIBMESH_CHKERRABORT(ierr);
578  ierr = VecScatterEnd(scatter,rhs->vec(),subrhs,INSERT_VALUES,SCATTER_FORWARD);
579  LIBMESH_CHKERRABORT(ierr);
580 
581  ierr = VecScatterBegin(scatter,solution->vec(),subsolution,INSERT_VALUES,SCATTER_FORWARD);
582  LIBMESH_CHKERRABORT(ierr);
583  ierr = VecScatterEnd(scatter,solution->vec(),subsolution,INSERT_VALUES,SCATTER_FORWARD);
584  LIBMESH_CHKERRABORT(ierr);
585 
586 #if PETSC_VERSION_LESS_THAN(3,1,0)
587  ierr = MatGetSubMatrix(matrix->mat(),
588  _restrict_solve_to_is,_restrict_solve_to_is,
589  PETSC_DECIDE,MAT_INITIAL_MATRIX,&submat);
590  LIBMESH_CHKERRABORT(ierr);
591  ierr = MatGetSubMatrix(precond->mat(),
592  _restrict_solve_to_is,_restrict_solve_to_is,
593  PETSC_DECIDE,MAT_INITIAL_MATRIX,&subprecond);
594  LIBMESH_CHKERRABORT(ierr);
595 #else
596  ierr = MatGetSubMatrix(matrix->mat(),
597  _restrict_solve_to_is,_restrict_solve_to_is,
598  MAT_INITIAL_MATRIX,&submat);
599  LIBMESH_CHKERRABORT(ierr);
600  ierr = MatGetSubMatrix(precond->mat(),
601  _restrict_solve_to_is,_restrict_solve_to_is,
602  MAT_INITIAL_MATRIX,&subprecond);
603  LIBMESH_CHKERRABORT(ierr);
604 #endif
605 
606  /* Since removing columns of the matrix changes the equation
607  system, we will now change the right hand side to compensate
608  for this. Note that this is not necessary if \p SUBSET_ZERO
609  has been selected. */
610  if(_subset_solve_mode!=SUBSET_ZERO)
611  {
612  _create_complement_is(rhs_in);
613  size_t is_complement_local_size = rhs_in.local_size()-is_local_size;
614 
615  Vec subvec1 = NULL;
616  Mat submat1 = NULL;
617  VecScatter scatter1 = NULL;
618 
619  ierr = VecCreate(this->comm().get(),&subvec1);
620  LIBMESH_CHKERRABORT(ierr);
621  ierr = VecSetSizes(subvec1,is_complement_local_size,PETSC_DECIDE);
622  LIBMESH_CHKERRABORT(ierr);
623  ierr = VecSetFromOptions(subvec1);
624  LIBMESH_CHKERRABORT(ierr);
625 
626  ierr = VecScatterCreate(rhs->vec(),_restrict_solve_to_is_complement, subvec1,NULL, &scatter1);
627  LIBMESH_CHKERRABORT(ierr);
628 
629  ierr = VecScatterBegin(scatter1,_subset_solve_mode==SUBSET_COPY_RHS ? rhs->vec() : solution->vec(),subvec1,INSERT_VALUES,SCATTER_FORWARD);
630  LIBMESH_CHKERRABORT(ierr);
631  ierr = VecScatterEnd(scatter1,_subset_solve_mode==SUBSET_COPY_RHS ? rhs->vec() : solution->vec(),subvec1,INSERT_VALUES,SCATTER_FORWARD);
632  LIBMESH_CHKERRABORT(ierr);
633 
634  ierr = VecScale(subvec1,-1.0);
635  LIBMESH_CHKERRABORT(ierr);
636 
637 #if PETSC_VERSION_LESS_THAN(3,1,0)
638  ierr = MatGetSubMatrix(matrix->mat(),
639  _restrict_solve_to_is,_restrict_solve_to_is_complement,
640  PETSC_DECIDE,MAT_INITIAL_MATRIX,&submat1);
641  LIBMESH_CHKERRABORT(ierr);
642 #else
643  ierr = MatGetSubMatrix(matrix->mat(),
644  _restrict_solve_to_is,_restrict_solve_to_is_complement,
645  MAT_INITIAL_MATRIX,&submat1);
646  LIBMESH_CHKERRABORT(ierr);
647 #endif
648 
649  ierr = MatMultAdd(submat1,subvec1,subrhs,subrhs);
650  LIBMESH_CHKERRABORT(ierr);
651 
652  ierr = LibMeshVecScatterDestroy(&scatter1);
653  LIBMESH_CHKERRABORT(ierr);
654  ierr = LibMeshVecDestroy(&subvec1);
655  LIBMESH_CHKERRABORT(ierr);
656  ierr = LibMeshMatDestroy(&submat1);
657  LIBMESH_CHKERRABORT(ierr);
658  }
659 
660  ierr = KSPSetOperators(_ksp, submat, subprecond,
661  this->same_preconditioner ? SAME_PRECONDITIONER : DIFFERENT_NONZERO_PATTERN);
662  LIBMESH_CHKERRABORT(ierr);
663 
664  if(this->_preconditioner)
665  {
666  subprecond_matrix = new PetscMatrix<Number>(subprecond,
667  this->comm());
668  this->_preconditioner->set_matrix(*subprecond_matrix);
669  this->_preconditioner->init();
670  }
671  }
672  else
673  {
674  ierr = KSPSetOperators(_ksp, matrix->mat(), precond->mat(),
675  this->same_preconditioner ? SAME_PRECONDITIONER : DIFFERENT_NONZERO_PATTERN);
676  LIBMESH_CHKERRABORT(ierr);
677 
678  if(this->_preconditioner)
679  {
680  this->_preconditioner->set_matrix(matrix_in);
681  this->_preconditioner->init();
682  }
683  }
684 
685  // Set the tolerances for the iterative solver. Use the user-supplied
686  // tolerance for the relative residual & leave the others at default values.
687  ierr = KSPSetTolerances (_ksp, tol, PETSC_DEFAULT,
688  PETSC_DEFAULT, max_its);
689  LIBMESH_CHKERRABORT(ierr);
690 
691  // Solve the linear system
692  if(_restrict_solve_to_is!=NULL)
693  {
694  ierr = KSPSolve (_ksp, subrhs, subsolution);
695  LIBMESH_CHKERRABORT(ierr);
696  }
697  else
698  {
699  ierr = KSPSolve (_ksp, rhs->vec(), solution->vec());
700  LIBMESH_CHKERRABORT(ierr);
701  }
702 
703  // Get the number of iterations required for convergence
704  ierr = KSPGetIterationNumber (_ksp, &its);
705  LIBMESH_CHKERRABORT(ierr);
706 
707  // Get the norm of the final residual to return to the user.
708  ierr = KSPGetResidualNorm (_ksp, &final_resid);
709  LIBMESH_CHKERRABORT(ierr);
710 
711  if(_restrict_solve_to_is!=NULL)
712  {
713  switch(_subset_solve_mode)
714  {
715  case SUBSET_ZERO:
716  ierr = VecZeroEntries(solution->vec());
717  LIBMESH_CHKERRABORT(ierr);
718  break;
719 
720  case SUBSET_COPY_RHS:
721  ierr = VecCopy(rhs->vec(),solution->vec());
722  LIBMESH_CHKERRABORT(ierr);
723  break;
724 
725  case SUBSET_DONT_TOUCH:
726  /* Nothing to do here. */
727  break;
728 
729  }
730  ierr = VecScatterBegin(scatter,subsolution,solution->vec(),INSERT_VALUES,SCATTER_REVERSE);
731  LIBMESH_CHKERRABORT(ierr);
732  ierr = VecScatterEnd(scatter,subsolution,solution->vec(),INSERT_VALUES,SCATTER_REVERSE);
733  LIBMESH_CHKERRABORT(ierr);
734 
735  ierr = LibMeshVecScatterDestroy(&scatter);
736  LIBMESH_CHKERRABORT(ierr);
737 
738  if(this->_preconditioner)
739  {
740  /* Before we delete subprecond_matrix, we should give the
741  _preconditioner a different matrix. */
742  this->_preconditioner->set_matrix(matrix_in);
743  this->_preconditioner->init();
744  delete subprecond_matrix;
745  subprecond_matrix = NULL;
746  }
747 
748  ierr = LibMeshVecDestroy(&subsolution);
749  LIBMESH_CHKERRABORT(ierr);
750  ierr = LibMeshVecDestroy(&subrhs);
751  LIBMESH_CHKERRABORT(ierr);
752  ierr = LibMeshMatDestroy(&submat);
753  LIBMESH_CHKERRABORT(ierr);
754  ierr = LibMeshMatDestroy(&subprecond);
755  LIBMESH_CHKERRABORT(ierr);
756  }
757 
758 #endif
759 
760  STOP_LOG("solve()", "PetscLinearSolver");
761  // return the # of its. and the final residual norm.
762  return std::make_pair(its, final_resid);
763 }
764 
765 template <typename T>
766 std::pair<unsigned int, Real>
768  NumericVector<T>& solution_in,
769  NumericVector<T>& rhs_in,
770  const double tol,
771  const unsigned int m_its)
772 {
773  START_LOG("solve()", "PetscLinearSolver");
774 
775  // Make sure the data passed in are really of Petsc types
776  PetscMatrix<T>* matrix = libmesh_cast_ptr<PetscMatrix<T>*>(&matrix_in);
777  // Note that the matrix and precond matrix are the same
778  PetscMatrix<T>* precond = libmesh_cast_ptr<PetscMatrix<T>*>(&matrix_in);
779  PetscVector<T>* solution = libmesh_cast_ptr<PetscVector<T>*>(&solution_in);
780  PetscVector<T>* rhs = libmesh_cast_ptr<PetscVector<T>*>(&rhs_in);
781 
782  this->init (matrix);
783 
785  PetscInt its=0, max_its = static_cast<PetscInt>(m_its);
786  PetscReal final_resid=0.;
787 
788  // Close the matrices and vectors in case this wasn't already done.
789  matrix->close ();
790  precond->close ();
791  solution->close ();
792  rhs->close ();
793 
794 // // If matrix != precond, then this means we have specified a
795 // // special preconditioner, so reset preconditioner type to PCMAT.
796 // if (matrix != precond)
797 // {
798 // this->_preconditioner_type = USER_PRECOND;
799 // this->set_petsc_preconditioner_type ();
800 // }
801 
802 // 2.1.x & earlier style
803 #if PETSC_VERSION_LESS_THAN(2,2,0)
804 
805  if(_restrict_solve_to_is!=NULL)
806  {
807  libmesh_not_implemented();
808  }
809 
810  // Based on http://wolfgang.math.tamu.edu/svn/public/deal.II/branches/MATH676/2008/deal.II/lac/source/petsc_solver.cc, http://tccc.iesl.forth.gr/AMS_EPEAEK/Elements/doc/in_html/petsc/SLES/index.html
811 
812  SLES sles;
813  ierr = SLESCreate (this->comm().get(), &sles);
814  LIBMESH_CHKERRABORT(ierr);
815 
816  ierr = SLESSetOperators (sles, matrix->mat(), precond->mat(), this->same_preconditioner ? SAME_PRECONDITIONER : DIFFERENT_NONZERO_PATTERN);
817  LIBMESH_CHKERRABORT(ierr);
818 
819  KSP ksp;
820  ierr = SLESGetKSP (sles, &ksp);
821  LIBMESH_CHKERRABORT(ierr);
822 
823  ierr = SLESSetUp (sles, rhs->vec(), solution->vec());
824  LIBMESH_CHKERRABORT(ierr);
825 
826  // See http://tccc.iesl.forth.gr/AMS_EPEAEK/Elements/doc/in_html/petsc/KSP/KSPSolveTrans.html#KSPSolveTrans
827  ierr = SLESSolveTrans (ksp, &its);
828  LIBMESH_CHKERRABORT(ierr);
829 
830 // 2.2.0
831 #elif PETSC_VERSION_LESS_THAN(2,2,1)
832 
833  if(_restrict_solve_to_is!=NULL)
834  {
835  libmesh_not_implemented();
836  }
837 
838  // Set operators. The input matrix works as the preconditioning matrix
839  // This was commented earlier but it looks like KSPSetOperators is supported
840  // after PETSc 2.2.0
841  ierr = KSPSetOperators(_ksp, matrix->mat(), precond->mat(),
842  this->same_preconditioner ? SAME_PRECONDITIONER : DIFFERENT_NONZERO_PATTERN);
843  LIBMESH_CHKERRABORT(ierr);
844 
845 
846  // Set the tolerances for the iterative solver. Use the user-supplied
847  // tolerance for the relative residual & leave the others at default values.
848  // Convergence is detected at iteration k if
849  // ||r_k||_2 < max(rtol*||b||_2 , abstol)
850  // where r_k is the residual vector and b is the right-hand side. Note that
851  // it is the *maximum* of the two values, the larger of which will almost
852  // always be rtol*||b||_2.
853  ierr = KSPSetTolerances (_ksp,
854  tol, // rtol = relative decrease in residual (1.e-5)
855  PETSC_DEFAULT, // abstol = absolute convergence tolerance (1.e-50)
856  PETSC_DEFAULT, // dtol = divergence tolerance (1.e+5)
857  max_its);
858  LIBMESH_CHKERRABORT(ierr);
859 
860 
861  // Set the solution vector to use
862  ierr = KSPSetSolution (_ksp, solution->vec());
863  LIBMESH_CHKERRABORT(ierr);
864 
865  // Set the RHS vector to use
866  ierr = KSPSetRhs (_ksp, rhs->vec());
867  LIBMESH_CHKERRABORT(ierr);
868 
869  // Solve the linear system
870  ierr = KSPSolveTranspose (_ksp);
871  LIBMESH_CHKERRABORT(ierr);
872 
873  // Get the number of iterations required for convergence
874  ierr = KSPGetIterationNumber (_ksp, &its);
875  LIBMESH_CHKERRABORT(ierr);
876 
877  // Get the norm of the final residual to return to the user.
878  ierr = KSPGetResidualNorm (_ksp, &final_resid);
879  LIBMESH_CHKERRABORT(ierr);
880 
881 // 2.2.1 & newer style
882 #else
883 
884  Mat submat = NULL;
885  Mat subprecond = NULL;
886  Vec subrhs = NULL;
887  Vec subsolution = NULL;
888  VecScatter scatter = NULL;
889  PetscMatrix<Number>* subprecond_matrix = NULL;
890 
891  // Set operators. Also restrict rhs and solution vector to
892  // subdomain if neccessary.
893  if(_restrict_solve_to_is!=NULL)
894  {
895  size_t is_local_size = this->_restrict_solve_to_is_local_size();
896 
897  ierr = VecCreate(this->comm().get(),&subrhs);
898  LIBMESH_CHKERRABORT(ierr);
899  ierr = VecSetSizes(subrhs,is_local_size,PETSC_DECIDE);
900  LIBMESH_CHKERRABORT(ierr);
901  ierr = VecSetFromOptions(subrhs);
902  LIBMESH_CHKERRABORT(ierr);
903 
904  ierr = VecCreate(this->comm().get(),&subsolution);
905  LIBMESH_CHKERRABORT(ierr);
906  ierr = VecSetSizes(subsolution,is_local_size,PETSC_DECIDE);
907  LIBMESH_CHKERRABORT(ierr);
908  ierr = VecSetFromOptions(subsolution);
909  LIBMESH_CHKERRABORT(ierr);
910 
911  ierr = VecScatterCreate(rhs->vec(),_restrict_solve_to_is, subrhs,NULL, &scatter);
912  LIBMESH_CHKERRABORT(ierr);
913 
914  ierr = VecScatterBegin(scatter,rhs->vec(),subrhs,INSERT_VALUES,SCATTER_FORWARD);
915  LIBMESH_CHKERRABORT(ierr);
916  ierr = VecScatterEnd(scatter,rhs->vec(),subrhs,INSERT_VALUES,SCATTER_FORWARD);
917  LIBMESH_CHKERRABORT(ierr);
918 
919  ierr = VecScatterBegin(scatter,solution->vec(),subsolution,INSERT_VALUES,SCATTER_FORWARD);
920  LIBMESH_CHKERRABORT(ierr);
921  ierr = VecScatterEnd(scatter,solution->vec(),subsolution,INSERT_VALUES,SCATTER_FORWARD);
922  LIBMESH_CHKERRABORT(ierr);
923 
924 #if PETSC_VERSION_LESS_THAN(3,1,0)
925  ierr = MatGetSubMatrix(matrix->mat(),
926  _restrict_solve_to_is,_restrict_solve_to_is,
927  PETSC_DECIDE,MAT_INITIAL_MATRIX,&submat);
928  LIBMESH_CHKERRABORT(ierr);
929  ierr = MatGetSubMatrix(precond->mat(),
930  _restrict_solve_to_is,_restrict_solve_to_is,
931  PETSC_DECIDE,MAT_INITIAL_MATRIX,&subprecond);
932  LIBMESH_CHKERRABORT(ierr);
933 #else
934  ierr = MatGetSubMatrix(matrix->mat(),
935  _restrict_solve_to_is,_restrict_solve_to_is,
936  MAT_INITIAL_MATRIX,&submat);
937  LIBMESH_CHKERRABORT(ierr);
938  ierr = MatGetSubMatrix(precond->mat(),
939  _restrict_solve_to_is,_restrict_solve_to_is,
940  MAT_INITIAL_MATRIX,&subprecond);
941  LIBMESH_CHKERRABORT(ierr);
942 #endif
943 
944  /* Since removing columns of the matrix changes the equation
945  system, we will now change the right hand side to compensate
946  for this. Note that this is not necessary if \p SUBSET_ZERO
947  has been selected. */
948  if(_subset_solve_mode!=SUBSET_ZERO)
949  {
950  _create_complement_is(rhs_in);
951  size_t is_complement_local_size = rhs_in.local_size()-is_local_size;
952 
953  Vec subvec1 = NULL;
954  Mat submat1 = NULL;
955  VecScatter scatter1 = NULL;
956 
957  ierr = VecCreate(this->comm().get(),&subvec1);
958  LIBMESH_CHKERRABORT(ierr);
959  ierr = VecSetSizes(subvec1,is_complement_local_size,PETSC_DECIDE);
960  LIBMESH_CHKERRABORT(ierr);
961  ierr = VecSetFromOptions(subvec1);
962  LIBMESH_CHKERRABORT(ierr);
963 
964  ierr = VecScatterCreate(rhs->vec(),_restrict_solve_to_is_complement, subvec1,NULL, &scatter1);
965  LIBMESH_CHKERRABORT(ierr);
966 
967  ierr = VecScatterBegin(scatter1,_subset_solve_mode==SUBSET_COPY_RHS ? rhs->vec() : solution->vec(),subvec1,INSERT_VALUES,SCATTER_FORWARD);
968  LIBMESH_CHKERRABORT(ierr);
969  ierr = VecScatterEnd(scatter1,_subset_solve_mode==SUBSET_COPY_RHS ? rhs->vec() : solution->vec(),subvec1,INSERT_VALUES,SCATTER_FORWARD);
970  LIBMESH_CHKERRABORT(ierr);
971 
972  ierr = VecScale(subvec1,-1.0);
973  LIBMESH_CHKERRABORT(ierr);
974 
975 #if PETSC_VERSION_LESS_THAN(3,1,0)
976  ierr = MatGetSubMatrix(matrix->mat(),
977  _restrict_solve_to_is,_restrict_solve_to_is_complement,
978  PETSC_DECIDE,MAT_INITIAL_MATRIX,&submat1);
979  LIBMESH_CHKERRABORT(ierr);
980 #else
981  ierr = MatGetSubMatrix(matrix->mat(),
982  _restrict_solve_to_is,_restrict_solve_to_is_complement,
983  MAT_INITIAL_MATRIX,&submat1);
984  LIBMESH_CHKERRABORT(ierr);
985 #endif
986 
987  ierr = MatMultAdd(submat1,subvec1,subrhs,subrhs);
988  LIBMESH_CHKERRABORT(ierr);
989 
990  ierr = LibMeshVecScatterDestroy(&scatter1);
991  LIBMESH_CHKERRABORT(ierr);
992  ierr = LibMeshVecDestroy(&subvec1);
993  LIBMESH_CHKERRABORT(ierr);
994  ierr = LibMeshMatDestroy(&submat1);
995  LIBMESH_CHKERRABORT(ierr);
996  }
997 
998  ierr = KSPSetOperators(_ksp, submat, subprecond,
999  this->same_preconditioner ? SAME_PRECONDITIONER : DIFFERENT_NONZERO_PATTERN);
1000  LIBMESH_CHKERRABORT(ierr);
1001 
1002  if(this->_preconditioner)
1003  {
1004  subprecond_matrix = new PetscMatrix<Number>(subprecond,
1005  this->comm());
1006  this->_preconditioner->set_matrix(*subprecond_matrix);
1007  this->_preconditioner->init();
1008  }
1009  }
1010  else
1011  {
1012  ierr = KSPSetOperators(_ksp, matrix->mat(), precond->mat(),
1013  this->same_preconditioner ? SAME_PRECONDITIONER : DIFFERENT_NONZERO_PATTERN);
1014  LIBMESH_CHKERRABORT(ierr);
1015 
1016  if(this->_preconditioner)
1017  {
1018  this->_preconditioner->set_matrix(matrix_in);
1019  this->_preconditioner->init();
1020  }
1021  }
1022 
1023  // Set the tolerances for the iterative solver. Use the user-supplied
1024  // tolerance for the relative residual & leave the others at default values.
1025  ierr = KSPSetTolerances (_ksp, tol, PETSC_DEFAULT,
1026  PETSC_DEFAULT, max_its);
1027  LIBMESH_CHKERRABORT(ierr);
1028 
1029  // Solve the linear system
1030  if(_restrict_solve_to_is!=NULL)
1031  {
1032  ierr = KSPSolveTranspose (_ksp, subrhs, subsolution);
1033  LIBMESH_CHKERRABORT(ierr);
1034  }
1035  else
1036  {
1037  ierr = KSPSolveTranspose (_ksp, rhs->vec(), solution->vec());
1038  LIBMESH_CHKERRABORT(ierr);
1039  }
1040 
1041  // Get the number of iterations required for convergence
1042  ierr = KSPGetIterationNumber (_ksp, &its);
1043  LIBMESH_CHKERRABORT(ierr);
1044 
1045  // Get the norm of the final residual to return to the user.
1046  ierr = KSPGetResidualNorm (_ksp, &final_resid);
1047  LIBMESH_CHKERRABORT(ierr);
1048 
1049  if(_restrict_solve_to_is!=NULL)
1050  {
1051  switch(_subset_solve_mode)
1052  {
1053  case SUBSET_ZERO:
1054  ierr = VecZeroEntries(solution->vec());
1055  LIBMESH_CHKERRABORT(ierr);
1056  break;
1057 
1058  case SUBSET_COPY_RHS:
1059  ierr = VecCopy(rhs->vec(),solution->vec());
1060  LIBMESH_CHKERRABORT(ierr);
1061  break;
1062 
1063  case SUBSET_DONT_TOUCH:
1064  /* Nothing to do here. */
1065  break;
1066 
1067  }
1068  ierr = VecScatterBegin(scatter,subsolution,solution->vec(),INSERT_VALUES,SCATTER_REVERSE);
1069  LIBMESH_CHKERRABORT(ierr);
1070  ierr = VecScatterEnd(scatter,subsolution,solution->vec(),INSERT_VALUES,SCATTER_REVERSE);
1071  LIBMESH_CHKERRABORT(ierr);
1072 
1073  ierr = LibMeshVecScatterDestroy(&scatter);
1074  LIBMESH_CHKERRABORT(ierr);
1075 
1076  if(this->_preconditioner)
1077  {
1078  /* Before we delete subprecond_matrix, we should give the
1079  _preconditioner a different matrix. */
1080  this->_preconditioner->set_matrix(matrix_in);
1081  this->_preconditioner->init();
1082  delete subprecond_matrix;
1083  subprecond_matrix = NULL;
1084  }
1085 
1086  ierr = LibMeshVecDestroy(&subsolution);
1087  LIBMESH_CHKERRABORT(ierr);
1088  ierr = LibMeshVecDestroy(&subrhs);
1089  LIBMESH_CHKERRABORT(ierr);
1090  ierr = LibMeshMatDestroy(&submat);
1091  LIBMESH_CHKERRABORT(ierr);
1092  ierr = LibMeshMatDestroy(&subprecond);
1093  LIBMESH_CHKERRABORT(ierr);
1094  }
1095 
1096 #endif
1097 
1098  STOP_LOG("solve()", "PetscLinearSolver");
1099  // return the # of its. and the final residual norm.
1100  return std::make_pair(its, final_resid);
1101 }
1102 
1103 
1104 template <typename T>
1105 std::pair<unsigned int, Real>
1107  NumericVector<T>& solution_in,
1108  NumericVector<T>& rhs_in,
1109  const double tol,
1110  const unsigned int m_its)
1111 {
1112 
1113 #if PETSC_VERSION_LESS_THAN(2,3,1)
1114  // FIXME[JWP]: There will be a bunch of unused variable warnings
1115  // for older PETScs here.
1116  libMesh::out << "This method has been developed with PETSc 2.3.1. "
1117  << "No one has made it backwards compatible with older "
1118  << "versions of PETSc so far; however, it might work "
1119  << "without any change with some older version." << std::endl;
1120  libmesh_error();
1121  return std::make_pair(0,0.0);
1122 
1123 #else
1124 
1125 #if PETSC_VERSION_LESS_THAN(3,1,0)
1126  if(_restrict_solve_to_is!=NULL)
1127  {
1128  libMesh::out << "The current implementation of subset solves with "
1129  << "shell matrices requires PETSc version 3.1 or above. "
1130  << "Older PETSc version do not support automatic "
1131  << "submatrix generation of shell matrices."
1132  << std::endl;
1133  libmesh_error();
1134  }
1135 #endif
1136 
1137  START_LOG("solve()", "PetscLinearSolver");
1138 
1139  // Make sure the data passed in are really of Petsc types
1140  PetscVector<T>* solution = libmesh_cast_ptr<PetscVector<T>*>(&solution_in);
1141  PetscVector<T>* rhs = libmesh_cast_ptr<PetscVector<T>*>(&rhs_in);
1142 
1143  this->init ();
1144 
1145  PetscErrorCode ierr=0;
1146  PetscInt its=0, max_its = static_cast<PetscInt>(m_its);
1147  PetscReal final_resid=0.;
1148 
1149  Mat submat = NULL;
1150  Vec subrhs = NULL;
1151  Vec subsolution = NULL;
1152  VecScatter scatter = NULL;
1153 
1154  // Close the matrices and vectors in case this wasn't already done.
1155  solution->close ();
1156  rhs->close ();
1157 
1158  // Prepare the matrix.
1159  Mat mat;
1160  ierr = MatCreateShell(this->comm().get(),
1161  rhs_in.local_size(),
1162  solution_in.local_size(),
1163  rhs_in.size(),
1164  solution_in.size(),
1165  const_cast<void*>(static_cast<const void*>(&shell_matrix)),
1166  &mat);
1167  /* Note that the const_cast above is only necessary because PETSc
1168  does not accept a const void*. Inside the member function
1169  _petsc_shell_matrix() below, the pointer is casted back to a
1170  const ShellMatrix<T>*. */
1171 
1172  LIBMESH_CHKERRABORT(ierr);
1173  ierr = MatShellSetOperation(mat,MATOP_MULT,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_mult));
1174  LIBMESH_CHKERRABORT(ierr);
1175  ierr = MatShellSetOperation(mat,MATOP_MULT_ADD,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_mult_add));
1176  LIBMESH_CHKERRABORT(ierr);
1177  ierr = MatShellSetOperation(mat,MATOP_GET_DIAGONAL,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_get_diagonal));
1178  LIBMESH_CHKERRABORT(ierr);
1179 
1180  // Restrict rhs and solution vectors and set operators. The input
1181  // matrix works as the preconditioning matrix.
1182  if(_restrict_solve_to_is!=NULL)
1183  {
1184  size_t is_local_size = this->_restrict_solve_to_is_local_size();
1185 
1186  ierr = VecCreate(this->comm().get(),&subrhs);
1187  LIBMESH_CHKERRABORT(ierr);
1188  ierr = VecSetSizes(subrhs,is_local_size,PETSC_DECIDE);
1189  LIBMESH_CHKERRABORT(ierr);
1190  ierr = VecSetFromOptions(subrhs);
1191  LIBMESH_CHKERRABORT(ierr);
1192 
1193  ierr = VecCreate(this->comm().get(),&subsolution);
1194  LIBMESH_CHKERRABORT(ierr);
1195  ierr = VecSetSizes(subsolution,is_local_size,PETSC_DECIDE);
1196  LIBMESH_CHKERRABORT(ierr);
1197  ierr = VecSetFromOptions(subsolution);
1198  LIBMESH_CHKERRABORT(ierr);
1199 
1200  ierr = VecScatterCreate(rhs->vec(),_restrict_solve_to_is, subrhs,NULL, &scatter);
1201  LIBMESH_CHKERRABORT(ierr);
1202 
1203  ierr = VecScatterBegin(scatter,rhs->vec(),subrhs,INSERT_VALUES,SCATTER_FORWARD);
1204  LIBMESH_CHKERRABORT(ierr);
1205  ierr = VecScatterEnd(scatter,rhs->vec(),subrhs,INSERT_VALUES,SCATTER_FORWARD);
1206  LIBMESH_CHKERRABORT(ierr);
1207 
1208  ierr = VecScatterBegin(scatter,solution->vec(),subsolution,INSERT_VALUES,SCATTER_FORWARD);
1209  LIBMESH_CHKERRABORT(ierr);
1210  ierr = VecScatterEnd(scatter,solution->vec(),subsolution,INSERT_VALUES,SCATTER_FORWARD);
1211  LIBMESH_CHKERRABORT(ierr);
1212 
1213 #if PETSC_VERSION_LESS_THAN(3,1,0)
1214  /* This point can't be reached, see above. */
1215  libmesh_assert(false);
1216 #else
1217  ierr = MatGetSubMatrix(mat,
1218  _restrict_solve_to_is,_restrict_solve_to_is,
1219  MAT_INITIAL_MATRIX,&submat);
1220  LIBMESH_CHKERRABORT(ierr);
1221 #endif
1222 
1223  /* Since removing columns of the matrix changes the equation
1224  system, we will now change the right hand side to compensate
1225  for this. Note that this is not necessary if \p SUBSET_ZERO
1226  has been selected. */
1227  if(_subset_solve_mode!=SUBSET_ZERO)
1228  {
1229  _create_complement_is(rhs_in);
1230  size_t is_complement_local_size = rhs_in.local_size()-is_local_size;
1231 
1232  Vec subvec1 = NULL;
1233  Mat submat1 = NULL;
1234  VecScatter scatter1 = NULL;
1235 
1236  ierr = VecCreate(this->comm().get(),&subvec1);
1237  LIBMESH_CHKERRABORT(ierr);
1238  ierr = VecSetSizes(subvec1,is_complement_local_size,PETSC_DECIDE);
1239  LIBMESH_CHKERRABORT(ierr);
1240  ierr = VecSetFromOptions(subvec1);
1241  LIBMESH_CHKERRABORT(ierr);
1242 
1243  ierr = VecScatterCreate(rhs->vec(),_restrict_solve_to_is_complement, subvec1,NULL, &scatter1);
1244  LIBMESH_CHKERRABORT(ierr);
1245 
1246  ierr = VecScatterBegin(scatter1,_subset_solve_mode==SUBSET_COPY_RHS ? rhs->vec() : solution->vec(),subvec1,INSERT_VALUES,SCATTER_FORWARD);
1247  LIBMESH_CHKERRABORT(ierr);
1248  ierr = VecScatterEnd(scatter1,_subset_solve_mode==SUBSET_COPY_RHS ? rhs->vec() : solution->vec(),subvec1,INSERT_VALUES,SCATTER_FORWARD);
1249  LIBMESH_CHKERRABORT(ierr);
1250 
1251  ierr = VecScale(subvec1,-1.0);
1252  LIBMESH_CHKERRABORT(ierr);
1253 
1254 #if PETSC_VERSION_LESS_THAN(3,1,0)
1255  /* This point can't be reached, see above. */
1256  libmesh_assert(false);
1257 #else
1258  ierr = MatGetSubMatrix(mat,
1259  _restrict_solve_to_is,_restrict_solve_to_is_complement,
1260  MAT_INITIAL_MATRIX,&submat1);
1261  LIBMESH_CHKERRABORT(ierr);
1262 #endif
1263 
1264  // The following lines would be correct, but don't work
1265  // correctly in PETSc up to 3.1.0-p5. See discussion in
1266  // petsc-users of Nov 9, 2010.
1267  //
1268  // ierr = MatMultAdd(submat1,subvec1,subrhs,subrhs);
1269  // LIBMESH_CHKERRABORT(ierr);
1270  //
1271  // We workaround by using a temporary vector. Note that the
1272  // fix in PETsc 3.1.0-p6 uses a temporary vector internally,
1273  // so this is no effective performance loss.
1274  Vec subvec2 = NULL;
1275  ierr = VecCreate(this->comm().get(),&subvec2);
1276  LIBMESH_CHKERRABORT(ierr);
1277  ierr = VecSetSizes(subvec2,is_local_size,PETSC_DECIDE);
1278  LIBMESH_CHKERRABORT(ierr);
1279  ierr = VecSetFromOptions(subvec2);
1280  LIBMESH_CHKERRABORT(ierr);
1281  ierr = MatMult(submat1,subvec1,subvec2);
1282  LIBMESH_CHKERRABORT(ierr);
1283  ierr = VecAXPY(subrhs,1.0,subvec2);
1284 
1285  ierr = LibMeshVecScatterDestroy(&scatter1);
1286  LIBMESH_CHKERRABORT(ierr);
1287  ierr = LibMeshVecDestroy(&subvec1);
1288  LIBMESH_CHKERRABORT(ierr);
1289  ierr = LibMeshMatDestroy(&submat1);
1290  LIBMESH_CHKERRABORT(ierr);
1291  }
1292 
1293  ierr = KSPSetOperators(_ksp, submat, submat,
1294  DIFFERENT_NONZERO_PATTERN);
1295  LIBMESH_CHKERRABORT(ierr);
1296  }
1297  else
1298  {
1299  ierr = KSPSetOperators(_ksp, mat, mat,
1300  DIFFERENT_NONZERO_PATTERN);
1301  LIBMESH_CHKERRABORT(ierr);
1302  }
1303 
1304  // Set the tolerances for the iterative solver. Use the user-supplied
1305  // tolerance for the relative residual & leave the others at default values.
1306  ierr = KSPSetTolerances (_ksp, tol, PETSC_DEFAULT,
1307  PETSC_DEFAULT, max_its);
1308  LIBMESH_CHKERRABORT(ierr);
1309 
1310  // Solve the linear system
1311  if(_restrict_solve_to_is!=NULL)
1312  {
1313  ierr = KSPSolve (_ksp, subrhs, subsolution);
1314  LIBMESH_CHKERRABORT(ierr);
1315  }
1316  else
1317  {
1318  ierr = KSPSolve (_ksp, rhs->vec(), solution->vec());
1319  LIBMESH_CHKERRABORT(ierr);
1320  }
1321 
1322  // Get the number of iterations required for convergence
1323  ierr = KSPGetIterationNumber (_ksp, &its);
1324  LIBMESH_CHKERRABORT(ierr);
1325 
1326  // Get the norm of the final residual to return to the user.
1327  ierr = KSPGetResidualNorm (_ksp, &final_resid);
1328  LIBMESH_CHKERRABORT(ierr);
1329 
1330  if(_restrict_solve_to_is!=NULL)
1331  {
1332  switch(_subset_solve_mode)
1333  {
1334  case SUBSET_ZERO:
1335  ierr = VecZeroEntries(solution->vec());
1336  LIBMESH_CHKERRABORT(ierr);
1337  break;
1338 
1339  case SUBSET_COPY_RHS:
1340  ierr = VecCopy(rhs->vec(),solution->vec());
1341  LIBMESH_CHKERRABORT(ierr);
1342  break;
1343 
1344  case SUBSET_DONT_TOUCH:
1345  /* Nothing to do here. */
1346  break;
1347 
1348  }
1349  ierr = VecScatterBegin(scatter,subsolution,solution->vec(),INSERT_VALUES,SCATTER_REVERSE);
1350  LIBMESH_CHKERRABORT(ierr);
1351  ierr = VecScatterEnd(scatter,subsolution,solution->vec(),INSERT_VALUES,SCATTER_REVERSE);
1352  LIBMESH_CHKERRABORT(ierr);
1353 
1354  ierr = LibMeshVecScatterDestroy(&scatter);
1355  LIBMESH_CHKERRABORT(ierr);
1356 
1357  ierr = LibMeshVecDestroy(&subsolution);
1358  LIBMESH_CHKERRABORT(ierr);
1359  ierr = LibMeshVecDestroy(&subrhs);
1360  LIBMESH_CHKERRABORT(ierr);
1361  ierr = LibMeshMatDestroy(&submat);
1362  LIBMESH_CHKERRABORT(ierr);
1363  }
1364 
1365  // Destroy the matrix.
1366  ierr = LibMeshMatDestroy(&mat);
1367  LIBMESH_CHKERRABORT(ierr);
1368 
1369  STOP_LOG("solve()", "PetscLinearSolver");
1370  // return the # of its. and the final residual norm.
1371  return std::make_pair(its, final_resid);
1372 
1373 #endif
1374 
1375 }
1376 
1377 
1378 
1379 template <typename T>
1380 std::pair<unsigned int, Real>
1382  const SparseMatrix<T>& precond_matrix,
1383  NumericVector<T> &solution_in,
1384  NumericVector<T> &rhs_in,
1385  const double tol,
1386  const unsigned int m_its)
1387 {
1388 
1389 #if PETSC_VERSION_LESS_THAN(2,3,1)
1390  // FIXME[JWP]: There will be a bunch of unused variable warnings
1391  // for older PETScs here.
1392  libMesh::out << "This method has been developed with PETSc 2.3.1. "
1393  << "No one has made it backwards compatible with older "
1394  << "versions of PETSc so far; however, it might work "
1395  << "without any change with some older version." << std::endl;
1396  libmesh_error();
1397  return std::make_pair(0,0.0);
1398 
1399 #else
1400 
1401 #if PETSC_VERSION_LESS_THAN(3,1,0)
1402  if(_restrict_solve_to_is!=NULL)
1403  {
1404  libMesh::out << "The current implementation of subset solves with "
1405  << "shell matrices requires PETSc version 3.1 or above. "
1406  << "Older PETSc version do not support automatic "
1407  << "submatrix generation of shell matrices."
1408  << std::endl;
1409  libmesh_error();
1410  }
1411 #endif
1412 
1413  START_LOG("solve()", "PetscLinearSolver");
1414 
1415  // Make sure the data passed in are really of Petsc types
1416  const PetscMatrix<T>* precond = libmesh_cast_ptr<const PetscMatrix<T>*>(&precond_matrix);
1417  PetscVector<T>* solution = libmesh_cast_ptr<PetscVector<T>*>(&solution_in);
1418  PetscVector<T>* rhs = libmesh_cast_ptr<PetscVector<T>*>(&rhs_in);
1419 
1420  this->init ();
1421 
1422  PetscErrorCode ierr=0;
1423  PetscInt its=0, max_its = static_cast<PetscInt>(m_its);
1424  PetscReal final_resid=0.;
1425 
1426  Mat submat = NULL;
1427  Mat subprecond = NULL;
1428  Vec subrhs = NULL;
1429  Vec subsolution = NULL;
1430  VecScatter scatter = NULL;
1431  PetscMatrix<Number>* subprecond_matrix = NULL;
1432 
1433  // Close the matrices and vectors in case this wasn't already done.
1434  solution->close ();
1435  rhs->close ();
1436 
1437  // Prepare the matrix.
1438  Mat mat;
1439  ierr = MatCreateShell(this->comm().get(),
1440  rhs_in.local_size(),
1441  solution_in.local_size(),
1442  rhs_in.size(),
1443  solution_in.size(),
1444  const_cast<void*>(static_cast<const void*>(&shell_matrix)),
1445  &mat);
1446  /* Note that the const_cast above is only necessary because PETSc
1447  does not accept a const void*. Inside the member function
1448  _petsc_shell_matrix() below, the pointer is casted back to a
1449  const ShellMatrix<T>*. */
1450 
1451  LIBMESH_CHKERRABORT(ierr);
1452  ierr = MatShellSetOperation(mat,MATOP_MULT,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_mult));
1453  LIBMESH_CHKERRABORT(ierr);
1454  ierr = MatShellSetOperation(mat,MATOP_MULT_ADD,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_mult_add));
1455  LIBMESH_CHKERRABORT(ierr);
1456  ierr = MatShellSetOperation(mat,MATOP_GET_DIAGONAL,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_get_diagonal));
1457  LIBMESH_CHKERRABORT(ierr);
1458 
1459  // Restrict rhs and solution vectors and set operators. The input
1460  // matrix works as the preconditioning matrix.
1461  if(_restrict_solve_to_is!=NULL)
1462  {
1463  size_t is_local_size = this->_restrict_solve_to_is_local_size();
1464 
1465  ierr = VecCreate(this->comm().get(),&subrhs);
1466  LIBMESH_CHKERRABORT(ierr);
1467  ierr = VecSetSizes(subrhs,is_local_size,PETSC_DECIDE);
1468  LIBMESH_CHKERRABORT(ierr);
1469  ierr = VecSetFromOptions(subrhs);
1470  LIBMESH_CHKERRABORT(ierr);
1471 
1472  ierr = VecCreate(this->comm().get(),&subsolution);
1473  LIBMESH_CHKERRABORT(ierr);
1474  ierr = VecSetSizes(subsolution,is_local_size,PETSC_DECIDE);
1475  LIBMESH_CHKERRABORT(ierr);
1476  ierr = VecSetFromOptions(subsolution);
1477  LIBMESH_CHKERRABORT(ierr);
1478 
1479  ierr = VecScatterCreate(rhs->vec(),_restrict_solve_to_is, subrhs,NULL, &scatter);
1480  LIBMESH_CHKERRABORT(ierr);
1481 
1482  ierr = VecScatterBegin(scatter,rhs->vec(),subrhs,INSERT_VALUES,SCATTER_FORWARD);
1483  LIBMESH_CHKERRABORT(ierr);
1484  ierr = VecScatterEnd(scatter,rhs->vec(),subrhs,INSERT_VALUES,SCATTER_FORWARD);
1485  LIBMESH_CHKERRABORT(ierr);
1486 
1487  ierr = VecScatterBegin(scatter,solution->vec(),subsolution,INSERT_VALUES,SCATTER_FORWARD);
1488  LIBMESH_CHKERRABORT(ierr);
1489  ierr = VecScatterEnd(scatter,solution->vec(),subsolution,INSERT_VALUES,SCATTER_FORWARD);
1490  LIBMESH_CHKERRABORT(ierr);
1491 
1492 #if PETSC_VERSION_LESS_THAN(3,1,0)
1493  /* This point can't be reached, see above. */
1494  libmesh_assert(false);
1495 #else
1496  ierr = MatGetSubMatrix(mat,
1497  _restrict_solve_to_is,_restrict_solve_to_is,
1498  MAT_INITIAL_MATRIX,&submat);
1499  LIBMESH_CHKERRABORT(ierr);
1500  ierr = MatGetSubMatrix(const_cast<PetscMatrix<T>*>(precond)->mat(),
1501  _restrict_solve_to_is,_restrict_solve_to_is,
1502  MAT_INITIAL_MATRIX,&subprecond);
1503  LIBMESH_CHKERRABORT(ierr);
1504 #endif
1505 
1506  /* Since removing columns of the matrix changes the equation
1507  system, we will now change the right hand side to compensate
1508  for this. Note that this is not necessary if \p SUBSET_ZERO
1509  has been selected. */
1510  if(_subset_solve_mode!=SUBSET_ZERO)
1511  {
1512  _create_complement_is(rhs_in);
1513  size_t is_complement_local_size = rhs_in.local_size()-is_local_size;
1514 
1515  Vec subvec1 = NULL;
1516  Mat submat1 = NULL;
1517  VecScatter scatter1 = NULL;
1518 
1519  ierr = VecCreate(this->comm().get(),&subvec1);
1520  LIBMESH_CHKERRABORT(ierr);
1521  ierr = VecSetSizes(subvec1,is_complement_local_size,PETSC_DECIDE);
1522  LIBMESH_CHKERRABORT(ierr);
1523  ierr = VecSetFromOptions(subvec1);
1524  LIBMESH_CHKERRABORT(ierr);
1525 
1526  ierr = VecScatterCreate(rhs->vec(),_restrict_solve_to_is_complement, subvec1,NULL, &scatter1);
1527  LIBMESH_CHKERRABORT(ierr);
1528 
1529  ierr = VecScatterBegin(scatter1,_subset_solve_mode==SUBSET_COPY_RHS ? rhs->vec() : solution->vec(),subvec1,INSERT_VALUES,SCATTER_FORWARD);
1530  LIBMESH_CHKERRABORT(ierr);
1531  ierr = VecScatterEnd(scatter1,_subset_solve_mode==SUBSET_COPY_RHS ? rhs->vec() : solution->vec(),subvec1,INSERT_VALUES,SCATTER_FORWARD);
1532  LIBMESH_CHKERRABORT(ierr);
1533 
1534  ierr = VecScale(subvec1,-1.0);
1535  LIBMESH_CHKERRABORT(ierr);
1536 
1537 #if PETSC_VERSION_LESS_THAN(3,1,0)
1538  /* This point can't be reached, see above. */
1539  libmesh_assert(false);
1540 #else
1541  ierr = MatGetSubMatrix(mat,
1542  _restrict_solve_to_is,_restrict_solve_to_is_complement,
1543  MAT_INITIAL_MATRIX,&submat1);
1544  LIBMESH_CHKERRABORT(ierr);
1545 #endif
1546 
1547  // The following lines would be correct, but don't work
1548  // correctly in PETSc up to 3.1.0-p5. See discussion in
1549  // petsc-users of Nov 9, 2010.
1550  //
1551  // ierr = MatMultAdd(submat1,subvec1,subrhs,subrhs);
1552  // LIBMESH_CHKERRABORT(ierr);
1553  //
1554  // We workaround by using a temporary vector. Note that the
1555  // fix in PETsc 3.1.0-p6 uses a temporary vector internally,
1556  // so this is no effective performance loss.
1557  Vec subvec2 = NULL;
1558  ierr = VecCreate(this->comm().get(),&subvec2);
1559  LIBMESH_CHKERRABORT(ierr);
1560  ierr = VecSetSizes(subvec2,is_local_size,PETSC_DECIDE);
1561  LIBMESH_CHKERRABORT(ierr);
1562  ierr = VecSetFromOptions(subvec2);
1563  LIBMESH_CHKERRABORT(ierr);
1564  ierr = MatMult(submat1,subvec1,subvec2);
1565  LIBMESH_CHKERRABORT(ierr);
1566  ierr = VecAXPY(subrhs,1.0,subvec2);
1567 
1568  ierr = LibMeshVecScatterDestroy(&scatter1);
1569  LIBMESH_CHKERRABORT(ierr);
1570  ierr = LibMeshVecDestroy(&subvec1);
1571  LIBMESH_CHKERRABORT(ierr);
1572  ierr = LibMeshMatDestroy(&submat1);
1573  LIBMESH_CHKERRABORT(ierr);
1574  }
1575 
1576  ierr = KSPSetOperators(_ksp, submat, subprecond,
1577  DIFFERENT_NONZERO_PATTERN);
1578  LIBMESH_CHKERRABORT(ierr);
1579 
1580  if(this->_preconditioner)
1581  {
1582  subprecond_matrix = new PetscMatrix<Number>(subprecond,
1583  this->comm());
1584  this->_preconditioner->set_matrix(*subprecond_matrix);
1585  this->_preconditioner->init();
1586  }
1587  }
1588  else
1589  {
1590  ierr = KSPSetOperators(_ksp, mat, const_cast<PetscMatrix<T>*>(precond)->mat(),
1591  DIFFERENT_NONZERO_PATTERN);
1592  LIBMESH_CHKERRABORT(ierr);
1593 
1594  if(this->_preconditioner)
1595  {
1596  this->_preconditioner->set_matrix(const_cast<SparseMatrix<Number>&>(precond_matrix));
1597  this->_preconditioner->init();
1598  }
1599  }
1600 
1601  // Set the tolerances for the iterative solver. Use the user-supplied
1602  // tolerance for the relative residual & leave the others at default values.
1603  ierr = KSPSetTolerances (_ksp, tol, PETSC_DEFAULT,
1604  PETSC_DEFAULT, max_its);
1605  LIBMESH_CHKERRABORT(ierr);
1606 
1607  // Solve the linear system
1608  if(_restrict_solve_to_is!=NULL)
1609  {
1610  ierr = KSPSolve (_ksp, subrhs, subsolution);
1611  LIBMESH_CHKERRABORT(ierr);
1612  }
1613  else
1614  {
1615  ierr = KSPSolve (_ksp, rhs->vec(), solution->vec());
1616  LIBMESH_CHKERRABORT(ierr);
1617  }
1618 
1619  // Get the number of iterations required for convergence
1620  ierr = KSPGetIterationNumber (_ksp, &its);
1621  LIBMESH_CHKERRABORT(ierr);
1622 
1623  // Get the norm of the final residual to return to the user.
1624  ierr = KSPGetResidualNorm (_ksp, &final_resid);
1625  LIBMESH_CHKERRABORT(ierr);
1626 
1627  if(_restrict_solve_to_is!=NULL)
1628  {
1629  switch(_subset_solve_mode)
1630  {
1631  case SUBSET_ZERO:
1632  ierr = VecZeroEntries(solution->vec());
1633  LIBMESH_CHKERRABORT(ierr);
1634  break;
1635 
1636  case SUBSET_COPY_RHS:
1637  ierr = VecCopy(rhs->vec(),solution->vec());
1638  LIBMESH_CHKERRABORT(ierr);
1639  break;
1640 
1641  case SUBSET_DONT_TOUCH:
1642  /* Nothing to do here. */
1643  break;
1644 
1645  }
1646  ierr = VecScatterBegin(scatter,subsolution,solution->vec(),INSERT_VALUES,SCATTER_REVERSE);
1647  LIBMESH_CHKERRABORT(ierr);
1648  ierr = VecScatterEnd(scatter,subsolution,solution->vec(),INSERT_VALUES,SCATTER_REVERSE);
1649  LIBMESH_CHKERRABORT(ierr);
1650 
1651  ierr = LibMeshVecScatterDestroy(&scatter);
1652  LIBMESH_CHKERRABORT(ierr);
1653 
1654  if(this->_preconditioner)
1655  {
1656  /* Before we delete subprecond_matrix, we should give the
1657  _preconditioner a different matrix. */
1658  this->_preconditioner->set_matrix(const_cast<SparseMatrix<Number>&>(precond_matrix));
1659  this->_preconditioner->init();
1660  delete subprecond_matrix;
1661  subprecond_matrix = NULL;
1662  }
1663 
1664  ierr = LibMeshVecDestroy(&subsolution);
1665  LIBMESH_CHKERRABORT(ierr);
1666  ierr = LibMeshVecDestroy(&subrhs);
1667  LIBMESH_CHKERRABORT(ierr);
1668  ierr = LibMeshMatDestroy(&submat);
1669  LIBMESH_CHKERRABORT(ierr);
1670  ierr = LibMeshMatDestroy(&subprecond);
1671  LIBMESH_CHKERRABORT(ierr);
1672  }
1673 
1674  // Destroy the matrix.
1675  ierr = LibMeshMatDestroy(&mat);
1676  LIBMESH_CHKERRABORT(ierr);
1677 
1678  STOP_LOG("solve()", "PetscLinearSolver");
1679  // return the # of its. and the final residual norm.
1680  return std::make_pair(its, final_resid);
1681 
1682 #endif
1683 
1684 }
1685 
1686 
1687 
1688 template <typename T>
1689 void PetscLinearSolver<T>::get_residual_history(std::vector<double>& hist)
1690 {
1691  PetscErrorCode ierr = 0;
1692  PetscInt its = 0;
1693 
1694  // Fill the residual history vector with the residual norms
1695  // Note that GetResidualHistory() does not copy any values, it
1696  // simply sets the pointer p. Note that for some Krylov subspace
1697  // methods, the number of residuals returned in the history
1698  // vector may be different from what you are expecting. For
1699  // example, TFQMR returns two residual values per iteration step.
1700  PetscReal* p;
1701  ierr = KSPGetResidualHistory(_ksp, &p, &its);
1702  LIBMESH_CHKERRABORT(ierr);
1703 
1704  // Check for early return
1705  if (its == 0) return;
1706 
1707  // Create space to store the result
1708  hist.resize(its);
1709 
1710  // Copy history into the vector provided by the user.
1711  for (PetscInt i=0; i<its; ++i)
1712  {
1713  hist[i] = *p;
1714  p++;
1715  }
1716 }
1717 
1718 
1719 
1720 
1721 template <typename T>
1723 {
1724  PetscErrorCode ierr = 0;
1725  PetscInt its = 0;
1726 
1727  // Fill the residual history vector with the residual norms
1728  // Note that GetResidualHistory() does not copy any values, it
1729  // simply sets the pointer p. Note that for some Krylov subspace
1730  // methods, the number of residuals returned in the history
1731  // vector may be different from what you are expecting. For
1732  // example, TFQMR returns two residual values per iteration step.
1733  PetscReal* p;
1734  ierr = KSPGetResidualHistory(_ksp, &p, &its);
1735  LIBMESH_CHKERRABORT(ierr);
1736 
1737  // Check no residual history
1738  if (its == 0)
1739  {
1740  libMesh::err << "No iterations have been performed, returning 0." << std::endl;
1741  return 0.;
1742  }
1743 
1744  // Otherwise, return the value pointed to by p.
1745  return *p;
1746 }
1747 
1748 
1749 
1750 
1751 template <typename T>
1753 {
1754  PetscErrorCode ierr = 0;
1755 
1756  switch (this->_solver_type)
1757  {
1758 
1759  case CG:
1760  ierr = KSPSetType (_ksp, (char*) KSPCG); LIBMESH_CHKERRABORT(ierr); return;
1761 
1762  case CR:
1763  ierr = KSPSetType (_ksp, (char*) KSPCR); LIBMESH_CHKERRABORT(ierr); return;
1764 
1765  case CGS:
1766  ierr = KSPSetType (_ksp, (char*) KSPCGS); LIBMESH_CHKERRABORT(ierr); return;
1767 
1768  case BICG:
1769  ierr = KSPSetType (_ksp, (char*) KSPBICG); LIBMESH_CHKERRABORT(ierr); return;
1770 
1771  case TCQMR:
1772  ierr = KSPSetType (_ksp, (char*) KSPTCQMR); LIBMESH_CHKERRABORT(ierr); return;
1773 
1774  case TFQMR:
1775  ierr = KSPSetType (_ksp, (char*) KSPTFQMR); LIBMESH_CHKERRABORT(ierr); return;
1776 
1777  case LSQR:
1778  ierr = KSPSetType (_ksp, (char*) KSPLSQR); LIBMESH_CHKERRABORT(ierr); return;
1779 
1780  case BICGSTAB:
1781  ierr = KSPSetType (_ksp, (char*) KSPBCGS); LIBMESH_CHKERRABORT(ierr); return;
1782 
1783  case MINRES:
1784  ierr = KSPSetType (_ksp, (char*) KSPMINRES); LIBMESH_CHKERRABORT(ierr); return;
1785 
1786  case GMRES:
1787  ierr = KSPSetType (_ksp, (char*) KSPGMRES); LIBMESH_CHKERRABORT(ierr); return;
1788 
1789  case RICHARDSON:
1790  ierr = KSPSetType (_ksp, (char*) KSPRICHARDSON); LIBMESH_CHKERRABORT(ierr); return;
1791 
1792  case CHEBYSHEV:
1793 #if defined(LIBMESH_HAVE_PETSC) && PETSC_VERSION_LESS_THAN(3,3,0)
1794  ierr = KSPSetType (_ksp, (char*) KSPCHEBYCHEV); LIBMESH_CHKERRABORT(ierr); return;
1795 #else
1796  ierr = KSPSetType (_ksp, (char*) KSPCHEBYSHEV); LIBMESH_CHKERRABORT(ierr); return;
1797 #endif
1798 
1799 
1800  default:
1801  libMesh::err << "ERROR: Unsupported PETSC Solver: "
1802  << Utility::enum_to_string(this->_solver_type) << std::endl
1803  << "Continuing with PETSC defaults" << std::endl;
1804  }
1805 }
1806 
1807 template <typename T>
1809 {
1810  KSPConvergedReason reason;
1811  KSPGetConvergedReason(_ksp, &reason);
1812  libMesh::out << "Linear solver convergence/divergence reason: " << KSPConvergedReasons[reason] << std::endl;
1813 }
1814 
1815 
1816 
1817 template <typename T>
1819 {
1820  /* Get the matrix context. */
1821  PetscErrorCode ierr=0;
1822  void* ctx;
1823  ierr = MatShellGetContext(mat,&ctx);
1824 
1825  /* Get user shell matrix object. */
1826  const ShellMatrix<T>& shell_matrix = *static_cast<const ShellMatrix<T>*>(ctx);
1827  CHKERRABORT(shell_matrix.comm().get(), ierr);
1828 
1829  /* Make \p NumericVector instances around the vectors. */
1830  PetscVector<T> arg_global(arg, shell_matrix.comm());
1831  PetscVector<T> dest_global(dest, shell_matrix.comm());
1832 
1833  /* Call the user function. */
1834  shell_matrix.vector_mult(dest_global,arg_global);
1835 
1836  return ierr;
1837 }
1838 
1839 
1840 
1841 template <typename T>
1843 {
1844  /* Get the matrix context. */
1845  PetscErrorCode ierr=0;
1846  void* ctx;
1847  ierr = MatShellGetContext(mat,&ctx);
1848 
1849  /* Get user shell matrix object. */
1850  const ShellMatrix<T>& shell_matrix = *static_cast<const ShellMatrix<T>*>(ctx);
1851  CHKERRABORT(shell_matrix.comm().get(), ierr);
1852 
1853  /* Make \p NumericVector instances around the vectors. */
1854  PetscVector<T> arg_global(arg, shell_matrix.comm());
1855  PetscVector<T> dest_global(dest, shell_matrix.comm());
1856  PetscVector<T> add_global(add, shell_matrix.comm());
1857 
1858  if(add!=arg)
1859  {
1860  arg_global = add_global;
1861  }
1862 
1863  /* Call the user function. */
1864  shell_matrix.vector_mult_add(dest_global,arg_global);
1865 
1866  return ierr;
1867 }
1868 
1869 
1870 
1871 template <typename T>
1873 {
1874  /* Get the matrix context. */
1875  PetscErrorCode ierr=0;
1876  void* ctx;
1877  ierr = MatShellGetContext(mat,&ctx);
1878 
1879  /* Get user shell matrix object. */
1880  const ShellMatrix<T>& shell_matrix = *static_cast<const ShellMatrix<T>*>(ctx);
1881  CHKERRABORT(shell_matrix.comm().get(), ierr);
1882 
1883  /* Make \p NumericVector instances around the vector. */
1884  PetscVector<T> dest_global(dest, shell_matrix.comm());
1885 
1886  /* Call the user function. */
1887  shell_matrix.get_diagonal(dest_global);
1888 
1889  return ierr;
1890 }
1891 
1892 
1893 
1894 //------------------------------------------------------------------
1895 // Explicit instantiations
1896 template class PetscLinearSolver<Number>;
1897 
1898 } // namespace libMesh
1899 
1900 
1901 
1902 #endif // #ifdef LIBMESH_HAVE_PETSC

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:06 UTC

Hosted By:
SourceForge.net Logo