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