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