libMesh::PetscLinearSolver< T > Class Template Reference

#include <petsc_linear_solver.h>

Inheritance diagram for libMesh::PetscLinearSolver< T >:

Public Member Functions

 PetscLinearSolver (const libMesh::Parallel::Communicator &comm LIBMESH_CAN_DEFAULT_TO_COMMWORLD)
 
 ~PetscLinearSolver ()
 
void clear ()
 
void init ()
 
void init (PetscMatrix< T > *matrix)
 
virtual void restrict_solve_to (const std::vector< unsigned int > *const dofs, const SubsetSolveMode subset_solve_mode=SUBSET_ZERO)
 
std::pair< unsigned int, Realsolve (SparseMatrix< T > &matrix_in, NumericVector< T > &solution_in, NumericVector< T > &rhs_in, const double tol, const unsigned int m_its)
 
std::pair< unsigned int, Realadjoint_solve (SparseMatrix< T > &matrix_in, NumericVector< T > &solution_in, NumericVector< T > &rhs_in, const double tol, const unsigned int m_its)
 
std::pair< unsigned int, Realsolve (SparseMatrix< T > &matrix, SparseMatrix< T > &preconditioner, NumericVector< T > &solution, NumericVector< T > &rhs, const double tol, const unsigned int m_its)
 
std::pair< unsigned int, Realsolve (const ShellMatrix< T > &shell_matrix, NumericVector< T > &solution_in, NumericVector< T > &rhs_in, const double tol, const unsigned int m_its)
 
virtual std::pair< unsigned
int, Real
solve (const ShellMatrix< T > &shell_matrix, const SparseMatrix< T > &precond_matrix, NumericVector< T > &solution_in, NumericVector< T > &rhs_in, const double tol, const unsigned int m_its)
 
PC pc ()
 
KSP ksp ()
 
void get_residual_history (std::vector< double > &hist)
 
Real get_initial_residual ()
 
virtual void print_converged_reason ()
 
bool initialized () const
 
SolverType solver_type () const
 
void set_solver_type (const SolverType st)
 
PreconditionerType preconditioner_type () const
 
void set_preconditioner_type (const PreconditionerType pct)
 
void attach_preconditioner (Preconditioner< T > *preconditioner)
 
virtual void reuse_preconditioner (bool)
 
bool get_same_preconditioner ()
 
std::pair< unsigned int, Realsolve (SparseMatrix< T > &matrix, SparseMatrix< T > *precond_matrix, NumericVector< T > &, NumericVector< T > &, const double, const unsigned int)
 
std::pair< unsigned int, Realsolve (const ShellMatrix< T > &matrix, const SparseMatrix< T > *precond_matrix, NumericVector< T > &, NumericVector< T > &, const double, const unsigned int)
 
const Parallel::Communicatorcomm () const
 
processor_id_type n_processors () const
 
processor_id_type processor_id () const
 

Static Public Member Functions

static AutoPtr< LinearSolver< T > > build (const libMesh::Parallel::Communicator &comm_in, const SolverPackage solver_package=libMesh::default_solver_package())
 
static std::string get_info ()
 
static void print_info (std::ostream &out=libMesh::out)
 
static unsigned int n_objects ()
 
static void enable_print_counter_info ()
 
static void disable_print_counter_info ()
 

Protected Types

typedef std::map< std::string,
std::pair< unsigned int,
unsigned int > > 
Counts
 

Protected Member Functions

void increment_constructor_count (const std::string &name)
 
void increment_destructor_count (const std::string &name)
 

Protected Attributes

SolverType _solver_type
 
PreconditionerType _preconditioner_type
 
bool _is_initialized
 
Preconditioner< T > * _preconditioner
 
bool same_preconditioner
 
const Parallel::Communicator_communicator
 

Static Protected Attributes

static Counts _counts
 
static Threads::atomic
< unsigned int > 
_n_objects
 
static Threads::spin_mutex _mutex
 
static bool _enable_print_counter = true
 

Private Member Functions

void set_petsc_solver_type ()
 
size_t _restrict_solve_to_is_local_size (void) const
 
void _create_complement_is (const NumericVector< T > &vec_in)
 

Static Private Member Functions

static PetscErrorCode _petsc_shell_matrix_mult (Mat mat, Vec arg, Vec dest)
 
static PetscErrorCode _petsc_shell_matrix_mult_add (Mat mat, Vec arg, Vec add, Vec dest)
 
static PetscErrorCode _petsc_shell_matrix_get_diagonal (Mat mat, Vec dest)
 

Private Attributes

SLES _sles
 
PC _pc
 
KSP _ksp
 
IS _restrict_solve_to_is
 
IS _restrict_solve_to_is_complement
 
SubsetSolveMode _subset_solve_mode
 

Detailed Description

template<typename T>
class libMesh::PetscLinearSolver< T >

This class provides an interface to PETSc iterative solvers that is compatible with the libMesh LinearSolver<>

Author
Benjamin Kirk, 2002-2007

Definition at line 98 of file petsc_linear_solver.h.

Member Typedef Documentation

typedef std::map<std::string, std::pair<unsigned int, unsigned int> > libMesh::ReferenceCounter::Counts
protectedinherited

Data structure to log the information. The log is identified by the class name.

Definition at line 113 of file reference_counter.h.

Constructor & Destructor Documentation

template<typename T >
libMesh::PetscLinearSolver< T >::PetscLinearSolver ( const libMesh::Parallel::Communicator &comm  LIBMESH_CAN_DEFAULT_TO_COMMWORLD)
inline

Constructor. Initializes Petsc data structures

Definition at line 325 of file petsc_linear_solver.h.

References libMesh::LinearSolver< T >::_preconditioner_type, libMeshEnums::BLOCK_JACOBI_PRECOND, libMeshEnums::ILU_PRECOND, and libMesh::ParallelObject::n_processors().

325  :
326  LinearSolver<T>(comm),
327  _restrict_solve_to_is(NULL),
330 {
331  if (this->n_processors() == 1)
333  else
335 }
template<typename T >
libMesh::PetscLinearSolver< T >::~PetscLinearSolver ( )
inline

Destructor.

Definition at line 341 of file petsc_linear_solver.h.

342 {
343  this->clear ();
344 }

Member Function Documentation

template<typename T >
void libMesh::PetscLinearSolver< T >::_create_complement_is ( const NumericVector< T > &  vec_in)
private

Creates _restrict_solve_to_is_complement to contain all indices that are local in vec_in, except those that are contained in _restrict_solve_to_is.

Definition at line 366 of file petsc_linear_solver.h.

References libMesh::ierr, and libMesh::libmesh_assert().

373 {
375 #if PETSC_VERSION_LESS_THAN(3,0,0)
376  // No ISComplement in PETSc 2.3.3
377  libmesh_not_implemented();
378 #else
380  {
381  int ierr = ISComplement(_restrict_solve_to_is,
382  vec_in.first_local_index(),
383  vec_in.last_local_index(),
385  LIBMESH_CHKERRABORT(ierr);
386  }
387 #endif
388 }
template<typename T >
PetscErrorCode libMesh::PetscLinearSolver< T >::_petsc_shell_matrix_get_diagonal ( Mat  mat,
Vec  dest 
)
staticprivate

Internal function if shell matrix mode is used.

Definition at line 1872 of file petsc_linear_solver.C.

References libMesh::CHKERRABORT(), libMesh::ParallelObject::comm(), libMesh::Parallel::Communicator::get(), libMesh::ShellMatrix< T >::get_diagonal(), and libMesh::ierr.

1873 {
1874  /* Get the matrix context. */
1875  PetscErrorCode ierr=0;
1876  void* ctx;
1877  ierr = MatShellGetContext(mat,&ctx);
1878 
1879  /* Get user shell matrix object. */
1880  const ShellMatrix<T>& shell_matrix = *static_cast<const ShellMatrix<T>*>(ctx);
1881  CHKERRABORT(shell_matrix.comm().get(), ierr);
1882 
1883  /* Make \p NumericVector instances around the vector. */
1884  PetscVector<T> dest_global(dest, shell_matrix.comm());
1885 
1886  /* Call the user function. */
1887  shell_matrix.get_diagonal(dest_global);
1888 
1889  return ierr;
1890 }
template<typename T >
PetscErrorCode libMesh::PetscLinearSolver< T >::_petsc_shell_matrix_mult ( Mat  mat,
Vec  arg,
Vec  dest 
)
staticprivate

Internal function if shell matrix mode is used.

Definition at line 1818 of file petsc_linear_solver.C.

References libMesh::CHKERRABORT(), libMesh::ParallelObject::comm(), libMesh::Parallel::Communicator::get(), libMesh::ierr, and libMesh::ShellMatrix< T >::vector_mult().

1819 {
1820  /* Get the matrix context. */
1821  PetscErrorCode ierr=0;
1822  void* ctx;
1823  ierr = MatShellGetContext(mat,&ctx);
1824 
1825  /* Get user shell matrix object. */
1826  const ShellMatrix<T>& shell_matrix = *static_cast<const ShellMatrix<T>*>(ctx);
1827  CHKERRABORT(shell_matrix.comm().get(), ierr);
1828 
1829  /* Make \p NumericVector instances around the vectors. */
1830  PetscVector<T> arg_global(arg, shell_matrix.comm());
1831  PetscVector<T> dest_global(dest, shell_matrix.comm());
1832 
1833  /* Call the user function. */
1834  shell_matrix.vector_mult(dest_global,arg_global);
1835 
1836  return ierr;
1837 }
template<typename T >
PetscErrorCode libMesh::PetscLinearSolver< T >::_petsc_shell_matrix_mult_add ( Mat  mat,
Vec  arg,
Vec  add,
Vec  dest 
)
staticprivate

Internal function if shell matrix mode is used.

Definition at line 1842 of file petsc_linear_solver.C.

References libMesh::CHKERRABORT(), libMesh::ParallelObject::comm(), libMesh::Parallel::Communicator::get(), libMesh::ierr, and libMesh::ShellMatrix< T >::vector_mult_add().

1843 {
1844  /* Get the matrix context. */
1845  PetscErrorCode ierr=0;
1846  void* ctx;
1847  ierr = MatShellGetContext(mat,&ctx);
1848 
1849  /* Get user shell matrix object. */
1850  const ShellMatrix<T>& shell_matrix = *static_cast<const ShellMatrix<T>*>(ctx);
1851  CHKERRABORT(shell_matrix.comm().get(), ierr);
1852 
1853  /* Make \p NumericVector instances around the vectors. */
1854  PetscVector<T> arg_global(arg, shell_matrix.comm());
1855  PetscVector<T> dest_global(dest, shell_matrix.comm());
1856  PetscVector<T> add_global(add, shell_matrix.comm());
1857 
1858  if(add!=arg)
1859  {
1860  arg_global = add_global;
1861  }
1862 
1863  /* Call the user function. */
1864  shell_matrix.vector_mult_add(dest_global,arg_global);
1865 
1866  return ierr;
1867 }
template<typename T >
size_t libMesh::PetscLinearSolver< T >::_restrict_solve_to_is_local_size ( void  ) const
inlineprivate

Internal method that returns the local size of _restrict_solve_to_is.

Definition at line 351 of file petsc_linear_solver.h.

References libMesh::ierr, and libMesh::libmesh_assert().

352 {
354 
355  PetscInt s;
356  int ierr = ISGetLocalSize(_restrict_solve_to_is,&s);
357  LIBMESH_CHKERRABORT(ierr);
358 
359  return static_cast<size_t>(s);
360 }
template<typename T >
std::pair< unsigned int, Real > libMesh::PetscLinearSolver< T >::adjoint_solve ( SparseMatrix< T > &  matrix_in,
NumericVector< T > &  solution_in,
NumericVector< T > &  rhs_in,
const double  tol,
const unsigned int  m_its 
)
virtual

Call the Petsc solver. It calls the method below, using the same matrix for the system and preconditioner matrices.

Reimplemented from libMesh::LinearSolver< T >.

Definition at line 767 of file petsc_linear_solver.C.

References libMesh::PetscMatrix< T >::close(), libMesh::comm, libMesh::ierr, libMesh::TriangleWrapper::init(), libMesh::PetscMatrix< T >::init(), libMesh::NumericVector< T >::local_size(), libMesh::PetscMatrix< T >::mat(), libMesh::START_LOG(), libMesh::STOP_LOG(), libMeshEnums::SUBSET_COPY_RHS, libMeshEnums::SUBSET_DONT_TOUCH, and libMeshEnums::SUBSET_ZERO.

772 {
773  START_LOG("solve()", "PetscLinearSolver");
774 
775  // Make sure the data passed in are really of Petsc types
776  PetscMatrix<T>* matrix = libmesh_cast_ptr<PetscMatrix<T>*>(&matrix_in);
777  // Note that the matrix and precond matrix are the same
778  PetscMatrix<T>* precond = libmesh_cast_ptr<PetscMatrix<T>*>(&matrix_in);
779  PetscVector<T>* solution = libmesh_cast_ptr<PetscVector<T>*>(&solution_in);
780  PetscVector<T>* rhs = libmesh_cast_ptr<PetscVector<T>*>(&rhs_in);
781 
782  this->init (matrix);
783 
785  PetscInt its=0, max_its = static_cast<PetscInt>(m_its);
786  PetscReal final_resid=0.;
787 
788  // Close the matrices and vectors in case this wasn't already done.
789  matrix->close ();
790  precond->close ();
791  solution->close ();
792  rhs->close ();
793 
794 // // If matrix != precond, then this means we have specified a
795 // // special preconditioner, so reset preconditioner type to PCMAT.
796 // if (matrix != precond)
797 // {
798 // this->_preconditioner_type = USER_PRECOND;
799 // this->set_petsc_preconditioner_type ();
800 // }
801 
802 // 2.1.x & earlier style
803 #if PETSC_VERSION_LESS_THAN(2,2,0)
804 
805  if(_restrict_solve_to_is!=NULL)
806  {
807  libmesh_not_implemented();
808  }
809 
810  // Based on http://wolfgang.math.tamu.edu/svn/public/deal.II/branches/MATH676/2008/deal.II/lac/source/petsc_solver.cc, http://tccc.iesl.forth.gr/AMS_EPEAEK/Elements/doc/in_html/petsc/SLES/index.html
811 
812  SLES sles;
813  ierr = SLESCreate (this->comm().get(), &sles);
814  LIBMESH_CHKERRABORT(ierr);
815 
816  ierr = SLESSetOperators (sles, matrix->mat(), precond->mat(), this->same_preconditioner ? SAME_PRECONDITIONER : DIFFERENT_NONZERO_PATTERN);
817  LIBMESH_CHKERRABORT(ierr);
818 
819  KSP ksp;
820  ierr = SLESGetKSP (sles, &ksp);
821  LIBMESH_CHKERRABORT(ierr);
822 
823  ierr = SLESSetUp (sles, rhs->vec(), solution->vec());
824  LIBMESH_CHKERRABORT(ierr);
825 
826  // See http://tccc.iesl.forth.gr/AMS_EPEAEK/Elements/doc/in_html/petsc/KSP/KSPSolveTrans.html#KSPSolveTrans
827  ierr = SLESSolveTrans (ksp, &its);
828  LIBMESH_CHKERRABORT(ierr);
829 
830 // 2.2.0
831 #elif PETSC_VERSION_LESS_THAN(2,2,1)
832 
833  if(_restrict_solve_to_is!=NULL)
834  {
835  libmesh_not_implemented();
836  }
837 
838  // Set operators. The input matrix works as the preconditioning matrix
839  // This was commented earlier but it looks like KSPSetOperators is supported
840  // after PETSc 2.2.0
841  ierr = KSPSetOperators(_ksp, matrix->mat(), precond->mat(),
842  this->same_preconditioner ? SAME_PRECONDITIONER : DIFFERENT_NONZERO_PATTERN);
843  LIBMESH_CHKERRABORT(ierr);
844 
845 
846  // Set the tolerances for the iterative solver. Use the user-supplied
847  // tolerance for the relative residual & leave the others at default values.
848  // Convergence is detected at iteration k if
849  // ||r_k||_2 < max(rtol*||b||_2 , abstol)
850  // where r_k is the residual vector and b is the right-hand side. Note that
851  // it is the *maximum* of the two values, the larger of which will almost
852  // always be rtol*||b||_2.
853  ierr = KSPSetTolerances (_ksp,
854  tol, // rtol = relative decrease in residual (1.e-5)
855  PETSC_DEFAULT, // abstol = absolute convergence tolerance (1.e-50)
856  PETSC_DEFAULT, // dtol = divergence tolerance (1.e+5)
857  max_its);
858  LIBMESH_CHKERRABORT(ierr);
859 
860 
861  // Set the solution vector to use
862  ierr = KSPSetSolution (_ksp, solution->vec());
863  LIBMESH_CHKERRABORT(ierr);
864 
865  // Set the RHS vector to use
866  ierr = KSPSetRhs (_ksp, rhs->vec());
867  LIBMESH_CHKERRABORT(ierr);
868 
869  // Solve the linear system
870  ierr = KSPSolveTranspose (_ksp);
871  LIBMESH_CHKERRABORT(ierr);
872 
873  // Get the number of iterations required for convergence
874  ierr = KSPGetIterationNumber (_ksp, &its);
875  LIBMESH_CHKERRABORT(ierr);
876 
877  // Get the norm of the final residual to return to the user.
878  ierr = KSPGetResidualNorm (_ksp, &final_resid);
879  LIBMESH_CHKERRABORT(ierr);
880 
881 // 2.2.1 & newer style
882 #else
883 
884  Mat submat = NULL;
885  Mat subprecond = NULL;
886  Vec subrhs = NULL;
887  Vec subsolution = NULL;
888  VecScatter scatter = NULL;
889  PetscMatrix<Number>* subprecond_matrix = NULL;
890 
891  // Set operators. Also restrict rhs and solution vector to
892  // subdomain if neccessary.
893  if(_restrict_solve_to_is!=NULL)
894  {
895  size_t is_local_size = this->_restrict_solve_to_is_local_size();
896 
897  ierr = VecCreate(this->comm().get(),&subrhs);
898  LIBMESH_CHKERRABORT(ierr);
899  ierr = VecSetSizes(subrhs,is_local_size,PETSC_DECIDE);
900  LIBMESH_CHKERRABORT(ierr);
901  ierr = VecSetFromOptions(subrhs);
902  LIBMESH_CHKERRABORT(ierr);
903 
904  ierr = VecCreate(this->comm().get(),&subsolution);
905  LIBMESH_CHKERRABORT(ierr);
906  ierr = VecSetSizes(subsolution,is_local_size,PETSC_DECIDE);
907  LIBMESH_CHKERRABORT(ierr);
908  ierr = VecSetFromOptions(subsolution);
909  LIBMESH_CHKERRABORT(ierr);
910 
911  ierr = VecScatterCreate(rhs->vec(),_restrict_solve_to_is, subrhs,NULL, &scatter);
912  LIBMESH_CHKERRABORT(ierr);
913 
914  ierr = VecScatterBegin(scatter,rhs->vec(),subrhs,INSERT_VALUES,SCATTER_FORWARD);
915  LIBMESH_CHKERRABORT(ierr);
916  ierr = VecScatterEnd(scatter,rhs->vec(),subrhs,INSERT_VALUES,SCATTER_FORWARD);
917  LIBMESH_CHKERRABORT(ierr);
918 
919  ierr = VecScatterBegin(scatter,solution->vec(),subsolution,INSERT_VALUES,SCATTER_FORWARD);
920  LIBMESH_CHKERRABORT(ierr);
921  ierr = VecScatterEnd(scatter,solution->vec(),subsolution,INSERT_VALUES,SCATTER_FORWARD);
922  LIBMESH_CHKERRABORT(ierr);
923 
924 #if PETSC_VERSION_LESS_THAN(3,1,0)
925  ierr = MatGetSubMatrix(matrix->mat(),
927  PETSC_DECIDE,MAT_INITIAL_MATRIX,&submat);
928  LIBMESH_CHKERRABORT(ierr);
929  ierr = MatGetSubMatrix(precond->mat(),
931  PETSC_DECIDE,MAT_INITIAL_MATRIX,&subprecond);
932  LIBMESH_CHKERRABORT(ierr);
933 #else
934  ierr = MatGetSubMatrix(matrix->mat(),
936  MAT_INITIAL_MATRIX,&submat);
937  LIBMESH_CHKERRABORT(ierr);
938  ierr = MatGetSubMatrix(precond->mat(),
940  MAT_INITIAL_MATRIX,&subprecond);
941  LIBMESH_CHKERRABORT(ierr);
942 #endif
943 
944  /* Since removing columns of the matrix changes the equation
945  system, we will now change the right hand side to compensate
946  for this. Note that this is not necessary if \p SUBSET_ZERO
947  has been selected. */
949  {
950  _create_complement_is(rhs_in);
951  size_t is_complement_local_size = rhs_in.local_size()-is_local_size;
952 
953  Vec subvec1 = NULL;
954  Mat submat1 = NULL;
955  VecScatter scatter1 = NULL;
956 
957  ierr = VecCreate(this->comm().get(),&subvec1);
958  LIBMESH_CHKERRABORT(ierr);
959  ierr = VecSetSizes(subvec1,is_complement_local_size,PETSC_DECIDE);
960  LIBMESH_CHKERRABORT(ierr);
961  ierr = VecSetFromOptions(subvec1);
962  LIBMESH_CHKERRABORT(ierr);
963 
964  ierr = VecScatterCreate(rhs->vec(),_restrict_solve_to_is_complement, subvec1,NULL, &scatter1);
965  LIBMESH_CHKERRABORT(ierr);
966 
967  ierr = VecScatterBegin(scatter1,_subset_solve_mode==SUBSET_COPY_RHS ? rhs->vec() : solution->vec(),subvec1,INSERT_VALUES,SCATTER_FORWARD);
968  LIBMESH_CHKERRABORT(ierr);
969  ierr = VecScatterEnd(scatter1,_subset_solve_mode==SUBSET_COPY_RHS ? rhs->vec() : solution->vec(),subvec1,INSERT_VALUES,SCATTER_FORWARD);
970  LIBMESH_CHKERRABORT(ierr);
971 
972  ierr = VecScale(subvec1,-1.0);
973  LIBMESH_CHKERRABORT(ierr);
974 
975 #if PETSC_VERSION_LESS_THAN(3,1,0)
976  ierr = MatGetSubMatrix(matrix->mat(),
978  PETSC_DECIDE,MAT_INITIAL_MATRIX,&submat1);
979  LIBMESH_CHKERRABORT(ierr);
980 #else
981  ierr = MatGetSubMatrix(matrix->mat(),
983  MAT_INITIAL_MATRIX,&submat1);
984  LIBMESH_CHKERRABORT(ierr);
985 #endif
986 
987  ierr = MatMultAdd(submat1,subvec1,subrhs,subrhs);
988  LIBMESH_CHKERRABORT(ierr);
989 
990  ierr = LibMeshVecScatterDestroy(&scatter1);
991  LIBMESH_CHKERRABORT(ierr);
992  ierr = LibMeshVecDestroy(&subvec1);
993  LIBMESH_CHKERRABORT(ierr);
994  ierr = LibMeshMatDestroy(&submat1);
995  LIBMESH_CHKERRABORT(ierr);
996  }
997 
998  ierr = KSPSetOperators(_ksp, submat, subprecond,
999  this->same_preconditioner ? SAME_PRECONDITIONER : DIFFERENT_NONZERO_PATTERN);
1000  LIBMESH_CHKERRABORT(ierr);
1001 
1002  if(this->_preconditioner)
1003  {
1004  subprecond_matrix = new PetscMatrix<Number>(subprecond,
1005  this->comm());
1006  this->_preconditioner->set_matrix(*subprecond_matrix);
1007  this->_preconditioner->init();
1008  }
1009  }
1010  else
1011  {
1012  ierr = KSPSetOperators(_ksp, matrix->mat(), precond->mat(),
1013  this->same_preconditioner ? SAME_PRECONDITIONER : DIFFERENT_NONZERO_PATTERN);
1014  LIBMESH_CHKERRABORT(ierr);
1015 
1016  if(this->_preconditioner)
1017  {
1018  this->_preconditioner->set_matrix(matrix_in);
1019  this->_preconditioner->init();
1020  }
1021  }
1022 
1023  // Set the tolerances for the iterative solver. Use the user-supplied
1024  // tolerance for the relative residual & leave the others at default values.
1025  ierr = KSPSetTolerances (_ksp, tol, PETSC_DEFAULT,
1026  PETSC_DEFAULT, max_its);
1027  LIBMESH_CHKERRABORT(ierr);
1028 
1029  // Solve the linear system
1030  if(_restrict_solve_to_is!=NULL)
1031  {
1032  ierr = KSPSolveTranspose (_ksp, subrhs, subsolution);
1033  LIBMESH_CHKERRABORT(ierr);
1034  }
1035  else
1036  {
1037  ierr = KSPSolveTranspose (_ksp, rhs->vec(), solution->vec());
1038  LIBMESH_CHKERRABORT(ierr);
1039  }
1040 
1041  // Get the number of iterations required for convergence
1042  ierr = KSPGetIterationNumber (_ksp, &its);
1043  LIBMESH_CHKERRABORT(ierr);
1044 
1045  // Get the norm of the final residual to return to the user.
1046  ierr = KSPGetResidualNorm (_ksp, &final_resid);
1047  LIBMESH_CHKERRABORT(ierr);
1048 
1049  if(_restrict_solve_to_is!=NULL)
1050  {
1051  switch(_subset_solve_mode)
1052  {
1053  case SUBSET_ZERO:
1054  ierr = VecZeroEntries(solution->vec());
1055  LIBMESH_CHKERRABORT(ierr);
1056  break;
1057 
1058  case SUBSET_COPY_RHS:
1059  ierr = VecCopy(rhs->vec(),solution->vec());
1060  LIBMESH_CHKERRABORT(ierr);
1061  break;
1062 
1063  case SUBSET_DONT_TOUCH:
1064  /* Nothing to do here. */
1065  break;
1066 
1067  }
1068  ierr = VecScatterBegin(scatter,subsolution,solution->vec(),INSERT_VALUES,SCATTER_REVERSE);
1069  LIBMESH_CHKERRABORT(ierr);
1070  ierr = VecScatterEnd(scatter,subsolution,solution->vec(),INSERT_VALUES,SCATTER_REVERSE);
1071  LIBMESH_CHKERRABORT(ierr);
1072 
1073  ierr = LibMeshVecScatterDestroy(&scatter);
1074  LIBMESH_CHKERRABORT(ierr);
1075 
1076  if(this->_preconditioner)
1077  {
1078  /* Before we delete subprecond_matrix, we should give the
1079  _preconditioner a different matrix. */
1080  this->_preconditioner->set_matrix(matrix_in);
1081  this->_preconditioner->init();
1082  delete subprecond_matrix;
1083  subprecond_matrix = NULL;
1084  }
1085 
1086  ierr = LibMeshVecDestroy(&subsolution);
1087  LIBMESH_CHKERRABORT(ierr);
1088  ierr = LibMeshVecDestroy(&subrhs);
1089  LIBMESH_CHKERRABORT(ierr);
1090  ierr = LibMeshMatDestroy(&submat);
1091  LIBMESH_CHKERRABORT(ierr);
1092  ierr = LibMeshMatDestroy(&subprecond);
1093  LIBMESH_CHKERRABORT(ierr);
1094  }
1095 
1096 #endif
1097 
1098  STOP_LOG("solve()", "PetscLinearSolver");
1099  // return the # of its. and the final residual norm.
1100  return std::make_pair(its, final_resid);
1101 }
template<typename T>
void libMesh::LinearSolver< T >::attach_preconditioner ( Preconditioner< T > *  preconditioner)
inherited

Attaches a Preconditioner object to be used

Definition at line 116 of file linear_solver.C.

References libMesh::libMeshPrivateData::_is_initialized, libMesh::err, and libMeshEnums::SHELL_PRECOND.

117 {
118  if(this->_is_initialized)
119  {
120  libMesh::err<<"Preconditioner must be attached before the solver is initialized!"<<std::endl;
121  libmesh_error();
122  }
123 
125  _preconditioner = preconditioner;
126 }
template<typename T >
AutoPtr< LinearSolver< T > > libMesh::LinearSolver< T >::build ( const libMesh::Parallel::Communicator comm_in,
const SolverPackage  solver_package = libMesh::default_solver_package() 
)
staticinherited

Builds a LinearSolver using the linear solver package specified by solver_package

Definition at line 40 of file linear_solver.C.

References libMesh::EIGEN_SOLVERS, libMesh::err, libMesh::LASPACK_SOLVERS, libMeshEnums::PETSC_SOLVERS, and libMesh::TRILINOS_SOLVERS.

Referenced by libMesh::ImplicitSystem::get_linear_solver().

42 {
43  // Build the appropriate solver
44  switch (solver_package)
45  {
46 
47 
48 #ifdef LIBMESH_HAVE_LASPACK
49  case LASPACK_SOLVERS:
50  {
51  AutoPtr<LinearSolver<T> > ap(new LaspackLinearSolver<T>(comm));
52  return ap;
53  }
54 #endif
55 
56 
57 #ifdef LIBMESH_HAVE_PETSC
58  case PETSC_SOLVERS:
59  {
60  AutoPtr<LinearSolver<T> > ap(new PetscLinearSolver<T>(comm));
61  return ap;
62  }
63 #endif
64 
65 
66 #ifdef LIBMESH_HAVE_TRILINOS
67  case TRILINOS_SOLVERS:
68  {
69  AutoPtr<LinearSolver<T> > ap(new AztecLinearSolver<T>(comm));
70  return ap;
71  }
72 #endif
73 
74 
75 #ifdef LIBMESH_HAVE_EIGEN
76  case EIGEN_SOLVERS:
77  {
78  AutoPtr<LinearSolver<T> > ap(new EigenSparseLinearSolver<T>(comm));
79  return ap;
80  }
81 #endif
82 
83  default:
84  libMesh::err << "ERROR: Unrecognized solver package: "
85  << solver_package
86  << std::endl;
87  libmesh_error();
88  }
89 
90  AutoPtr<LinearSolver<T> > ap(NULL);
91  return ap;
92 }
template<typename T >
void libMesh::PetscLinearSolver< T >::clear ( )
virtual

Release all memory and clear data structures.

Reimplemented from libMesh::LinearSolver< T >.

Definition at line 110 of file petsc_linear_solver.C.

References libMesh::libMeshPrivateData::_is_initialized, libMeshEnums::BLOCK_JACOBI_PRECOND, libMeshEnums::GMRES, libMesh::ierr, libMeshEnums::ILU_PRECOND, libMesh::initialized(), and libMesh::n_processors().

111 {
112  if (this->initialized())
113  {
114  /* If we were restricted to some subset, this restriction must
115  be removed and the subset index set destroyed. */
116  if(_restrict_solve_to_is!=NULL)
117  {
118  PetscErrorCode ierr = LibMeshISDestroy(&_restrict_solve_to_is);
119  LIBMESH_CHKERRABORT(ierr);
120  _restrict_solve_to_is = NULL;
121  }
122 
124  {
126  LIBMESH_CHKERRABORT(ierr);
128  }
129 
130  this->_is_initialized = false;
131 
133 
134 #if PETSC_VERSION_LESS_THAN(2,2,0)
135 
136  // 2.1.x & earlier style
137  ierr = SLESDestroy(_sles);
138  LIBMESH_CHKERRABORT(ierr);
139 
140 #else
141 
142  // 2.2.0 & newer style
143  ierr = LibMeshKSPDestroy(&_ksp);
144  LIBMESH_CHKERRABORT(ierr);
145 
146 #endif
147 
148 
149  // Mimic PETSc default solver and preconditioner
150  this->_solver_type = GMRES;
151 
152  if(!this->_preconditioner)
153  {
154  if (this->n_processors() == 1)
156  else
158  }
159  }
160 }
const Parallel::Communicator& libMesh::ParallelObject::comm ( ) const
inlineinherited
Returns
a reference to the Parallel::Communicator object used by this mesh.

Definition at line 86 of file parallel_object.h.

References libMesh::ParallelObject::_communicator.

Referenced by libMesh::__libmesh_petsc_diff_solver_jacobian(), libMesh::__libmesh_petsc_diff_solver_monitor(), libMesh::__libmesh_petsc_diff_solver_residual(), libMesh::__libmesh_petsc_snes_jacobian(), libMesh::__libmesh_petsc_snes_residual(), libMesh::MeshRefinement::_coarsen_elements(), libMesh::ExactSolution::_compute_error(), libMesh::MetisPartitioner::_do_partition(), libMesh::ParmetisPartitioner::_do_repartition(), libMesh::UniformRefinementEstimator::_estimate_error(), libMesh::SlepcEigenSolver< T >::_petsc_shell_matrix_get_diagonal(), libMesh::PetscLinearSolver< T >::_petsc_shell_matrix_get_diagonal(), libMesh::SlepcEigenSolver< T >::_petsc_shell_matrix_mult(), libMesh::PetscLinearSolver< T >::_petsc_shell_matrix_mult(), libMesh::PetscLinearSolver< T >::_petsc_shell_matrix_mult_add(), libMesh::EquationSystems::_read_impl(), libMesh::MeshRefinement::_refine_elements(), libMesh::ParallelMesh::add_elem(), libMesh::ImplicitSystem::add_matrix(), libMesh::ParallelMesh::add_node(), libMesh::System::add_vector(), libMesh::UnstructuredMesh::all_second_order(), libMesh::LaplaceMeshSmoother::allgather_graph(), libMesh::FEMSystem::assemble_qoi(), libMesh::MeshCommunication::assign_global_indices(), libMesh::ParmetisPartitioner::assign_partitioning(), libMesh::DofMap::attach_matrix(), libMesh::MeshTools::bounding_box(), libMesh::System::calculate_norm(), libMesh::MeshRefinement::coarsen_elements(), libMesh::Nemesis_IO_Helper::compute_num_global_elem_blocks(), libMesh::Nemesis_IO_Helper::compute_num_global_nodesets(), libMesh::Nemesis_IO_Helper::compute_num_global_sidesets(), libMesh::Problem_Interface::computeF(), libMesh::Problem_Interface::computeJacobian(), libMesh::Problem_Interface::computePreconditioner(), libMesh::MeshTools::correct_node_proc_ids(), libMesh::MeshCommunication::delete_remote_elements(), libMesh::DofMap::distribute_dofs(), DMlibMeshFunction(), DMlibMeshJacobian(), DMLibMeshSetSystem(), DMVariableBounds_libMesh(), libMesh::MeshRefinement::eliminate_unrefined_patches(), libMesh::WeightedPatchRecoveryErrorEstimator::estimate_error(), libMesh::PatchRecoveryErrorEstimator::estimate_error(), libMesh::JumpErrorEstimator::estimate_error(), libMesh::AdjointRefinementEstimator::estimate_error(), libMesh::MeshRefinement::flag_elements_by_elem_fraction(), libMesh::MeshRefinement::flag_elements_by_error_fraction(), libMesh::MeshRefinement::flag_elements_by_nelem_target(), libMesh::for(), libMesh::CondensedEigenSystem::get_eigenpair(), libMesh::ImplicitSystem::get_linear_solver(), libMesh::LocationMap< T >::init(), libMesh::PetscDiffSolver::init(), libMesh::TimeSolver::init(), libMesh::SystemSubsetBySubdomain::init(), libMesh::EigenSystem::init_data(), libMesh::EigenSystem::init_matrices(), libMesh::ParmetisPartitioner::initialize(), libMesh::MeshTools::libmesh_assert_valid_dof_ids(), libMesh::ParallelMesh::libmesh_assert_valid_parallel_flags(), libMesh::MeshTools::libmesh_assert_valid_procids< Elem >(), libMesh::MeshTools::libmesh_assert_valid_procids< Node >(), libMesh::MeshTools::libmesh_assert_valid_refinement_flags(), libMesh::MeshRefinement::limit_level_mismatch_at_edge(), libMesh::MeshRefinement::limit_level_mismatch_at_node(), libMesh::MeshRefinement::make_coarsening_compatible(), libMesh::MeshCommunication::make_elems_parallel_consistent(), libMesh::MeshRefinement::make_flags_parallel_consistent(), libMesh::MeshCommunication::make_node_ids_parallel_consistent(), libMesh::MeshCommunication::make_node_proc_ids_parallel_consistent(), libMesh::MeshCommunication::make_nodes_parallel_consistent(), libMesh::MeshRefinement::make_refinement_compatible(), libMesh::FEMSystem::mesh_position_set(), libMesh::MeshSerializer::MeshSerializer(), libMesh::ParallelMesh::n_active_elem(), libMesh::MeshTools::n_active_levels(), libMesh::BoundaryInfo::n_boundary_conds(), libMesh::BoundaryInfo::n_edge_conds(), libMesh::CondensedEigenSystem::n_global_non_condensed_dofs(), libMesh::MeshTools::n_levels(), libMesh::BoundaryInfo::n_nodeset_conds(), libMesh::MeshTools::n_p_levels(), libMesh::ParallelMesh::parallel_max_elem_id(), libMesh::ParallelMesh::parallel_max_node_id(), libMesh::ParallelMesh::parallel_n_elem(), libMesh::ParallelMesh::parallel_n_nodes(), libMesh::Partitioner::partition(), libMesh::Partitioner::partition_unpartitioned_elements(), libMesh::System::point_gradient(), libMesh::System::point_hessian(), libMesh::System::point_value(), libMesh::MeshBase::prepare_for_use(), libMesh::System::project_vector(), libMesh::Nemesis_IO::read(), libMesh::XdrIO::read(), libMesh::System::read_header(), libMesh::System::read_legacy_data(), libMesh::System::read_SCALAR_dofs(), libMesh::XdrIO::read_serialized_bc_names(), libMesh::XdrIO::read_serialized_bcs(), libMesh::System::read_serialized_blocked_dof_objects(), libMesh::XdrIO::read_serialized_connectivity(), libMesh::XdrIO::read_serialized_nodes(), libMesh::XdrIO::read_serialized_nodesets(), libMesh::XdrIO::read_serialized_subdomain_names(), libMesh::System::read_serialized_vector(), libMesh::MeshBase::recalculate_n_partitions(), libMesh::MeshRefinement::refine_and_coarsen_elements(), libMesh::MeshRefinement::refine_elements(), libMesh::Partitioner::set_node_processor_ids(), libMesh::DofMap::set_nonlocal_dof_objects(), libMesh::LaplaceMeshSmoother::smooth(), libMesh::MeshBase::subdomain_ids(), libMesh::BoundaryInfo::sync(), libMesh::Parallel::sync_element_data_by_parent_id(), libMesh::MeshRefinement::test_level_one(), libMesh::MeshRefinement::test_unflagged(), libMesh::MeshTools::total_weight(), libMesh::CheckpointIO::write(), libMesh::XdrIO::write(), libMesh::UnstructuredMesh::write(), libMesh::LegacyXdrIO::write_mesh(), libMesh::System::write_SCALAR_dofs(), libMesh::XdrIO::write_serialized_bcs(), libMesh::System::write_serialized_blocked_dof_objects(), libMesh::XdrIO::write_serialized_connectivity(), libMesh::XdrIO::write_serialized_nodes(), libMesh::XdrIO::write_serialized_nodesets(), and libMesh::DivaIO::write_stream().

87  { return _communicator; }
void libMesh::ReferenceCounter::disable_print_counter_info ( )
staticinherited

Definition at line 106 of file reference_counter.C.

References libMesh::ReferenceCounter::_enable_print_counter.

107 {
108  _enable_print_counter = false;
109  return;
110 }
void libMesh::ReferenceCounter::enable_print_counter_info ( )
staticinherited

Methods to enable/disable the reference counter output from print_info()

Definition at line 100 of file reference_counter.C.

References libMesh::ReferenceCounter::_enable_print_counter.

101 {
102  _enable_print_counter = true;
103  return;
104 }
std::string libMesh::ReferenceCounter::get_info ( )
staticinherited

Gets a string containing the reference information.

Definition at line 47 of file reference_counter.C.

References libMesh::ReferenceCounter::_counts, and libMesh::Quality::name().

Referenced by libMesh::ReferenceCounter::print_info().

48 {
49 #if defined(LIBMESH_ENABLE_REFERENCE_COUNTING) && defined(DEBUG)
50 
51  std::ostringstream oss;
52 
53  oss << '\n'
54  << " ---------------------------------------------------------------------------- \n"
55  << "| Reference count information |\n"
56  << " ---------------------------------------------------------------------------- \n";
57 
58  for (Counts::iterator it = _counts.begin();
59  it != _counts.end(); ++it)
60  {
61  const std::string name(it->first);
62  const unsigned int creations = it->second.first;
63  const unsigned int destructions = it->second.second;
64 
65  oss << "| " << name << " reference count information:\n"
66  << "| Creations: " << creations << '\n'
67  << "| Destructions: " << destructions << '\n';
68  }
69 
70  oss << " ---------------------------------------------------------------------------- \n";
71 
72  return oss.str();
73 
74 #else
75 
76  return "";
77 
78 #endif
79 }
template<typename T >
Real libMesh::PetscLinearSolver< T >::get_initial_residual ( )

Returns just the initial residual for the solve just completed with this interface. Use this method instead of the one above if you just want the starting residual and not the entire history.

Definition at line 1722 of file petsc_linear_solver.C.

References libMesh::err, and libMesh::ierr.

1723 {
1724  PetscErrorCode ierr = 0;
1725  PetscInt its = 0;
1726 
1727  // Fill the residual history vector with the residual norms
1728  // Note that GetResidualHistory() does not copy any values, it
1729  // simply sets the pointer p. Note that for some Krylov subspace
1730  // methods, the number of residuals returned in the history
1731  // vector may be different from what you are expecting. For
1732  // example, TFQMR returns two residual values per iteration step.
1733  PetscReal* p;
1734  ierr = KSPGetResidualHistory(_ksp, &p, &its);
1735  LIBMESH_CHKERRABORT(ierr);
1736 
1737  // Check no residual history
1738  if (its == 0)
1739  {
1740  libMesh::err << "No iterations have been performed, returning 0." << std::endl;
1741  return 0.;
1742  }
1743 
1744  // Otherwise, return the value pointed to by p.
1745  return *p;
1746 }
template<typename T >
void libMesh::PetscLinearSolver< T >::get_residual_history ( std::vector< double > &  hist)

Fills the input vector with the sequence of residual norms from the latest iterative solve.

Definition at line 1689 of file petsc_linear_solver.C.

References libMesh::ierr.

1690 {
1691  PetscErrorCode ierr = 0;
1692  PetscInt its = 0;
1693 
1694  // Fill the residual history vector with the residual norms
1695  // Note that GetResidualHistory() does not copy any values, it
1696  // simply sets the pointer p. Note that for some Krylov subspace
1697  // methods, the number of residuals returned in the history
1698  // vector may be different from what you are expecting. For
1699  // example, TFQMR returns two residual values per iteration step.
1700  PetscReal* p;
1701  ierr = KSPGetResidualHistory(_ksp, &p, &its);
1702  LIBMESH_CHKERRABORT(ierr);
1703 
1704  // Check for early return
1705  if (its == 0) return;
1706 
1707  // Create space to store the result
1708  hist.resize(its);
1709 
1710  // Copy history into the vector provided by the user.
1711  for (PetscInt i=0; i<its; ++i)
1712  {
1713  hist[i] = *p;
1714  p++;
1715  }
1716 }
template<typename T >
bool libMesh::LinearSolver< T >::get_same_preconditioner ( )
inlineinherited

Definition at line 290 of file linear_solver.h.

291 {
292  return same_preconditioner;
293 }
void libMesh::ReferenceCounter::increment_constructor_count ( const std::string &  name)
inlineprotectedinherited

Increments the construction counter. Should be called in the constructor of any derived class that will be reference counted.

Definition at line 163 of file reference_counter.h.

References libMesh::ReferenceCounter::_counts, libMesh::Quality::name(), and libMesh::Threads::spin_mtx.

Referenced by libMesh::ReferenceCountedObject< RBParametrized >::ReferenceCountedObject().

164 {
165  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
166  std::pair<unsigned int, unsigned int>& p = _counts[name];
167 
168  p.first++;
169 }
void libMesh::ReferenceCounter::increment_destructor_count ( const std::string &  name)
inlineprotectedinherited

Increments the destruction counter. Should be called in the destructor of any derived class that will be reference counted.

Definition at line 176 of file reference_counter.h.

References libMesh::ReferenceCounter::_counts, libMesh::Quality::name(), and libMesh::Threads::spin_mtx.

Referenced by libMesh::ReferenceCountedObject< RBParametrized >::~ReferenceCountedObject().

177 {
178  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
179  std::pair<unsigned int, unsigned int>& p = _counts[name];
180 
181  p.second++;
182 }
template<typename T >
void libMesh::PetscLinearSolver< T >::init ( )
virtual

Initialize data structures if not done so already.

Implements libMesh::LinearSolver< T >.

Definition at line 165 of file petsc_linear_solver.C.

References libMesh::__libmesh_petsc_preconditioner_apply(), libMesh::__libmesh_petsc_preconditioner_setup(), libMesh::libMeshPrivateData::_is_initialized, libMesh::comm, libMesh::ierr, libMesh::initialized(), and libMesh::PetscPreconditioner< T >::set_petsc_preconditioner_type().

Referenced by libMesh::PetscLinearSolver< T >::ksp(), and libMesh::PetscLinearSolver< T >::pc().

166 {
167  // Initialize the data structures if not done so already.
168  if (!this->initialized())
169  {
170  this->_is_initialized = true;
171 
173 
174 // 2.1.x & earlier style
175 #if PETSC_VERSION_LESS_THAN(2,2,0)
176 
177  // Create the linear solver context
178  ierr = SLESCreate (this->comm().get(), &_sles);
179  LIBMESH_CHKERRABORT(ierr);
180 
181  // Create the Krylov subspace & preconditioner contexts
182  ierr = SLESGetKSP (_sles, &_ksp);
183  LIBMESH_CHKERRABORT(ierr);
184  ierr = SLESGetPC (_sles, &_pc);
185  LIBMESH_CHKERRABORT(ierr);
186 
187  // Set user-specified solver and preconditioner types
188  this->set_petsc_solver_type();
189 
190  // Set the options from user-input
191  // Set runtime options, e.g.,
192  // -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
193  // These options will override those specified above as long as
194  // SLESSetFromOptions() is called _after_ any other customization
195  // routines.
196 
197  ierr = SLESSetFromOptions (_sles);
198  LIBMESH_CHKERRABORT(ierr);
199 
200 // 2.2.0 & newer style
201 #else
202 
203  // Create the linear solver context
204  ierr = KSPCreate (this->comm().get(), &_ksp);
205  LIBMESH_CHKERRABORT(ierr);
206 
207  // Create the preconditioner context
208  ierr = KSPGetPC (_ksp, &_pc);
209  LIBMESH_CHKERRABORT(ierr);
210 
211  // Set user-specified solver and preconditioner types
212  this->set_petsc_solver_type();
213 
214  // Set the options from user-input
215  // Set runtime options, e.g.,
216  // -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
217  // These options will override those specified above as long as
218  // KSPSetFromOptions() is called _after_ any other customization
219  // routines.
220 
221  ierr = KSPSetFromOptions (_ksp);
222  LIBMESH_CHKERRABORT(ierr);
223 
224  // Not sure if this is necessary, or if it is already handled by KSPSetFromOptions?
225  // NOT NECESSARY!!!!
226  //ierr = PCSetFromOptions (_pc);
227  //LIBMESH_CHKERRABORT(ierr);
228 
229 
230 #endif
231 
232  // Have the Krylov subspace method use our good initial guess
233  // rather than 0, unless the user requested a KSPType of
234  // preonly, which complains if asked to use initial guesses.
235 #if PETSC_VERSION_LESS_THAN(3,0,0) || !PETSC_RELEASE_LESS_THAN(3,4,0)
236  // Pre-3.0 and petsc-dev (as of October 2012) use non-const versions
237  KSPType ksp_type;
238 #else
239  const KSPType ksp_type;
240 #endif
241 
242  ierr = KSPGetType (_ksp, &ksp_type);
243  LIBMESH_CHKERRABORT(ierr);
244 
245  if (strcmp(ksp_type, "preonly"))
246  {
247  ierr = KSPSetInitialGuessNonzero (_ksp, PETSC_TRUE);
248  LIBMESH_CHKERRABORT(ierr);
249  }
250 
251  // Notify PETSc of location to store residual history.
252  // This needs to be called before any solves, since
253  // it sets the residual history length to zero. The default
254  // behavior is for PETSc to allocate (internally) an array
255  // of size 1000 to hold the residual norm history.
256  ierr = KSPSetResidualHistory(_ksp,
257  PETSC_NULL, // pointer to the array which holds the history
258  PETSC_DECIDE, // size of the array holding the history
259  PETSC_TRUE); // Whether or not to reset the history for each solve.
260  LIBMESH_CHKERRABORT(ierr);
261 
263 
264  //If there is a preconditioner object we need to set the internal setup and apply routines
265  if(this->_preconditioner)
266  {
267  this->_preconditioner->init();
268  PCShellSetContext(_pc,(void*)this->_preconditioner);
269  PCShellSetSetUp(_pc,__libmesh_petsc_preconditioner_setup);
270  PCShellSetApply(_pc,__libmesh_petsc_preconditioner_apply);
271  }
272  }
273 }
template<typename T >
void libMesh::PetscLinearSolver< T >::init ( PetscMatrix< T > *  matrix)

Initialize data structures if not done so already plus much more

Definition at line 277 of file petsc_linear_solver.C.

References libMesh::__libmesh_petsc_preconditioner_apply(), libMesh::__libmesh_petsc_preconditioner_setup(), libMesh::libMeshPrivateData::_is_initialized, libMesh::comm, libMesh::ierr, libMesh::initialized(), libMesh::PetscMatrix< T >::mat(), and libMesh::PetscPreconditioner< T >::set_petsc_preconditioner_type().

278 {
279  // Initialize the data structures if not done so already.
280  if (!this->initialized())
281  {
282  this->_is_initialized = true;
283 
285 
286 // 2.1.x & earlier style
287 #if PETSC_VERSION_LESS_THAN(2,2,0)
288 
289  // Create the linear solver context
290  ierr = SLESCreate (this->comm().get(), &_sles);
291  LIBMESH_CHKERRABORT(ierr);
292 
293  // Create the Krylov subspace & preconditioner contexts
294  ierr = SLESGetKSP (_sles, &_ksp);
295  LIBMESH_CHKERRABORT(ierr);
296  ierr = SLESGetPC (_sles, &_pc);
297  LIBMESH_CHKERRABORT(ierr);
298 
299  // Set user-specified solver and preconditioner types
300  this->set_petsc_solver_type();
301 
302  // Set the options from user-input
303  // Set runtime options, e.g.,
304  // -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
305  // These options will override those specified above as long as
306  // SLESSetFromOptions() is called _after_ any other customization
307  // routines.
308 
309  ierr = SLESSetFromOptions (_sles);
310  LIBMESH_CHKERRABORT(ierr);
311 
312 // 2.2.0 & newer style
313 #else
314 
315  // Create the linear solver context
316  ierr = KSPCreate (this->comm().get(), &_ksp);
317  LIBMESH_CHKERRABORT(ierr);
318 
319 
320  //ierr = PCCreate (this->comm().get(), &_pc);
321  // LIBMESH_CHKERRABORT(ierr);
322 
323  // Create the preconditioner context
324  ierr = KSPGetPC (_ksp, &_pc);
325  LIBMESH_CHKERRABORT(ierr);
326 
327  // Set operators. The input matrix works as the preconditioning matrix
328  ierr = KSPSetOperators(_ksp, matrix->mat(), matrix->mat(),DIFFERENT_NONZERO_PATTERN);
329  LIBMESH_CHKERRABORT(ierr);
330 
331  // Set user-specified solver and preconditioner types
332  this->set_petsc_solver_type();
333 
334  // Set the options from user-input
335  // Set runtime options, e.g.,
336  // -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
337  // These options will override those specified above as long as
338  // KSPSetFromOptions() is called _after_ any other customization
339  // routines.
340 
341  ierr = KSPSetFromOptions (_ksp);
342  LIBMESH_CHKERRABORT(ierr);
343 
344  // Not sure if this is necessary, or if it is already handled by KSPSetFromOptions?
345  // NOT NECESSARY!!!!
346  //ierr = PCSetFromOptions (_pc);
347  //LIBMESH_CHKERRABORT(ierr);
348 
349 
350 #endif
351 
352  // Have the Krylov subspace method use our good initial guess
353  // rather than 0, unless the user requested a KSPType of
354  // preonly, which complains if asked to use initial guesses.
355 #if PETSC_VERSION_LESS_THAN(3,0,0) || !PETSC_RELEASE_LESS_THAN(3,4,0)
356  KSPType ksp_type;
357 #else
358  const KSPType ksp_type;
359 #endif
360 
361  ierr = KSPGetType (_ksp, &ksp_type);
362  LIBMESH_CHKERRABORT(ierr);
363 
364  if (strcmp(ksp_type, "preonly"))
365  {
366  ierr = KSPSetInitialGuessNonzero (_ksp, PETSC_TRUE);
367  LIBMESH_CHKERRABORT(ierr);
368  }
369 
370  // Notify PETSc of location to store residual history.
371  // This needs to be called before any solves, since
372  // it sets the residual history length to zero. The default
373  // behavior is for PETSc to allocate (internally) an array
374  // of size 1000 to hold the residual norm history.
375  ierr = KSPSetResidualHistory(_ksp,
376  PETSC_NULL, // pointer to the array which holds the history
377  PETSC_DECIDE, // size of the array holding the history
378  PETSC_TRUE); // Whether or not to reset the history for each solve.
379  LIBMESH_CHKERRABORT(ierr);
380 
382  if(this->_preconditioner)
383  {
384  this->_preconditioner->set_matrix(*matrix);
385  this->_preconditioner->init();
386  PCShellSetContext(_pc,(void*)this->_preconditioner);
387  PCShellSetSetUp(_pc,__libmesh_petsc_preconditioner_setup);
388  PCShellSetApply(_pc,__libmesh_petsc_preconditioner_apply);
389  }
390  }
391 }
template<typename T>
bool libMesh::LinearSolver< T >::initialized ( ) const
inlineinherited
Returns
true if the data structures are initialized, false otherwise.

Definition at line 85 of file linear_solver.h.

85 { return _is_initialized; }
template<typename T >
KSP libMesh::PetscLinearSolver< T >::ksp ( )
inline

Returns the raw PETSc ksp context pointer. This is useful if you are for example setting a custom convergence test with KSPSetConvergenceTest().

Definition at line 222 of file petsc_linear_solver.h.

References libMesh::PetscLinearSolver< T >::_ksp, and libMesh::PetscLinearSolver< T >::init().

222 { this->init(); return _ksp; }
static unsigned int libMesh::ReferenceCounter::n_objects ( )
inlinestaticinherited

Prints the number of outstanding (created, but not yet destroyed) objects.

Definition at line 79 of file reference_counter.h.

References libMesh::ReferenceCounter::_n_objects.

80  { return _n_objects; }
processor_id_type libMesh::ParallelObject::n_processors ( ) const
inlineinherited
Returns
the number of processors in the group.

Definition at line 92 of file parallel_object.h.

References libMesh::ParallelObject::_communicator, and libMesh::Parallel::Communicator::size().

Referenced by libMesh::ParmetisPartitioner::_do_repartition(), libMesh::ParallelMesh::add_elem(), libMesh::ParallelMesh::add_node(), libMesh::LaplaceMeshSmoother::allgather_graph(), libMesh::ParmetisPartitioner::assign_partitioning(), libMesh::ParallelMesh::assign_unique_ids(), libMesh::AztecLinearSolver< T >::AztecLinearSolver(), libMesh::ParallelMesh::clear(), libMesh::Nemesis_IO_Helper::compute_border_node_ids(), libMesh::Nemesis_IO_Helper::construct_nemesis_filename(), libMesh::UnstructuredMesh::create_pid_mesh(), libMesh::DofMap::distribute_dofs(), libMesh::DofMap::distribute_local_dofs_node_major(), libMesh::DofMap::distribute_local_dofs_var_major(), libMesh::EnsightIO::EnsightIO(), libMesh::MeshBase::get_info(), libMesh::EquationSystems::init(), libMesh::SystemSubsetBySubdomain::init(), libMesh::ParmetisPartitioner::initialize(), libMesh::Nemesis_IO_Helper::initialize(), libMesh::MeshTools::libmesh_assert_valid_dof_ids(), libMesh::MeshTools::libmesh_assert_valid_procids< Elem >(), libMesh::MeshTools::libmesh_assert_valid_procids< Node >(), libMesh::MeshTools::libmesh_assert_valid_refinement_flags(), libMesh::MeshBase::n_active_elem_on_proc(), libMesh::MeshBase::n_elem_on_proc(), libMesh::MeshBase::n_nodes_on_proc(), libMesh::Partitioner::partition(), libMesh::MeshBase::partition(), libMesh::Partitioner::partition_unpartitioned_elements(), libMesh::PetscLinearSolver< T >::PetscLinearSolver(), libMesh::System::point_gradient(), libMesh::System::point_hessian(), libMesh::System::point_value(), libMesh::MeshTools::processor_bounding_box(), libMesh::System::project_vector(), libMesh::Nemesis_IO::read(), libMesh::CheckpointIO::read(), libMesh::UnstructuredMesh::read(), libMesh::System::read_parallel_data(), libMesh::System::read_SCALAR_dofs(), libMesh::System::read_serialized_blocked_dof_objects(), libMesh::System::read_serialized_vector(), libMesh::Partitioner::repartition(), libMesh::Partitioner::set_node_processor_ids(), libMesh::DofMap::set_nonlocal_dof_objects(), libMesh::BoundaryInfo::sync(), libMesh::ParallelMesh::update_parallel_id_counts(), libMesh::CheckpointIO::write(), libMesh::GMVIO::write_binary(), libMesh::GMVIO::write_discontinuous_gmv(), libMesh::System::write_parallel_data(), libMesh::System::write_SCALAR_dofs(), libMesh::XdrIO::write_serialized_bcs(), libMesh::System::write_serialized_blocked_dof_objects(), libMesh::XdrIO::write_serialized_connectivity(), libMesh::XdrIO::write_serialized_nodes(), and libMesh::XdrIO::write_serialized_nodesets().

93  { return libmesh_cast_int<processor_id_type>(_communicator.size()); }
template<typename T >
PC libMesh::PetscLinearSolver< T >::pc ( )
inline

Returns the raw PETSc preconditioner context pointer. This allows you to specify the PCShellSetApply() and PCShellSetSetUp() functions if you desire. Just don't do anything crazy like calling PCDestroy()!

Definition at line 215 of file petsc_linear_solver.h.

References libMesh::PetscLinearSolver< T >::_pc, and libMesh::PetscLinearSolver< T >::init().

215 { this->init(); return _pc; }
template<typename T >
PreconditionerType libMesh::LinearSolver< T >::preconditioner_type ( ) const
inherited

Returns the type of preconditioner to use.

Definition at line 96 of file linear_solver.C.

97 {
98  if(_preconditioner)
99  return _preconditioner->type();
100 
101  return _preconditioner_type;
102 }
template<typename T >
void libMesh::PetscLinearSolver< T >::print_converged_reason ( )
virtual

Prints a useful message about why the latest linear solve con(di)verged.

Implements libMesh::LinearSolver< T >.

Definition at line 1808 of file petsc_linear_solver.C.

References libMesh::out.

1809 {
1810  KSPConvergedReason reason;
1811  KSPGetConvergedReason(_ksp, &reason);
1812  libMesh::out << "Linear solver convergence/divergence reason: " << KSPConvergedReasons[reason] << std::endl;
1813 }
void libMesh::ReferenceCounter::print_info ( std::ostream &  out = libMesh::out)
staticinherited

Prints the reference information, by default to libMesh::out.

Definition at line 88 of file reference_counter.C.

References libMesh::ReferenceCounter::_enable_print_counter, and libMesh::ReferenceCounter::get_info().

89 {
91 }
processor_id_type libMesh::ParallelObject::processor_id ( ) const
inlineinherited
Returns
the rank of this processor in the group.

Definition at line 98 of file parallel_object.h.

References libMesh::ParallelObject::_communicator, and libMesh::Parallel::Communicator::rank().

Referenced by libMesh::MetisPartitioner::_do_partition(), libMesh::EquationSystems::_read_impl(), libMesh::SerialMesh::active_local_elements_begin(), libMesh::ParallelMesh::active_local_elements_begin(), libMesh::SerialMesh::active_local_elements_end(), libMesh::ParallelMesh::active_local_elements_end(), libMesh::SerialMesh::active_local_subdomain_elements_begin(), libMesh::ParallelMesh::active_local_subdomain_elements_begin(), libMesh::SerialMesh::active_local_subdomain_elements_end(), libMesh::ParallelMesh::active_local_subdomain_elements_end(), libMesh::SerialMesh::active_not_local_elements_begin(), libMesh::ParallelMesh::active_not_local_elements_begin(), libMesh::SerialMesh::active_not_local_elements_end(), libMesh::ParallelMesh::active_not_local_elements_end(), libMesh::ParallelMesh::add_elem(), libMesh::DofMap::add_neighbors_to_send_list(), libMesh::ParallelMesh::add_node(), libMesh::UnstructuredMesh::all_second_order(), libMesh::ParmetisPartitioner::assign_partitioning(), libMesh::EquationSystems::build_discontinuous_solution_vector(), libMesh::Nemesis_IO_Helper::build_element_and_node_maps(), libMesh::ParmetisPartitioner::build_graph(), libMesh::InfElemBuilder::build_inf_elem(), libMesh::DofMap::build_sparsity(), libMesh::ParallelMesh::clear(), libMesh::ExodusII_IO_Helper::close(), libMesh::Nemesis_IO_Helper::compute_border_node_ids(), libMesh::Nemesis_IO_Helper::compute_communication_map_parameters(), libMesh::Nemesis_IO_Helper::compute_internal_and_border_elems_and_internal_nodes(), libMesh::Nemesis_IO_Helper::compute_node_communication_maps(), libMesh::Nemesis_IO_Helper::compute_num_global_elem_blocks(), libMesh::Nemesis_IO_Helper::compute_num_global_nodesets(), libMesh::Nemesis_IO_Helper::compute_num_global_sidesets(), libMesh::Nemesis_IO_Helper::construct_nemesis_filename(), libMesh::ExodusII_IO_Helper::create(), libMesh::DofMap::distribute_dofs(), libMesh::DofMap::distribute_local_dofs_node_major(), libMesh::DofMap::distribute_local_dofs_var_major(), libMesh::DofMap::end_dof(), libMesh::DofMap::end_old_dof(), libMesh::EnsightIO::EnsightIO(), libMesh::UnstructuredMesh::find_neighbors(), libMesh::DofMap::first_dof(), libMesh::DofMap::first_old_dof(), libMesh::Nemesis_IO_Helper::get_cmap_params(), libMesh::Nemesis_IO_Helper::get_eb_info_global(), libMesh::Nemesis_IO_Helper::get_elem_cmap(), libMesh::Nemesis_IO_Helper::get_elem_map(), libMesh::MeshBase::get_info(), libMesh::Nemesis_IO_Helper::get_init_global(), libMesh::Nemesis_IO_Helper::get_init_info(), libMesh::Nemesis_IO_Helper::get_loadbal_param(), libMesh::Nemesis_IO_Helper::get_node_cmap(), libMesh::Nemesis_IO_Helper::get_node_map(), libMesh::Nemesis_IO_Helper::get_ns_param_global(), libMesh::Nemesis_IO_Helper::get_ss_param_global(), libMesh::MeshFunction::gradient(), libMesh::MeshFunction::hessian(), libMesh::SystemSubsetBySubdomain::init(), libMesh::ParmetisPartitioner::initialize(), libMesh::ExodusII_IO_Helper::initialize(), libMesh::ExodusII_IO_Helper::initialize_discontinuous(), libMesh::ExodusII_IO_Helper::initialize_element_variables(), libMesh::ExodusII_IO_Helper::initialize_global_variables(), libMesh::ExodusII_IO_Helper::initialize_nodal_variables(), libMesh::SparsityPattern::Build::join(), libMesh::DofMap::last_dof(), libMesh::MeshTools::libmesh_assert_valid_procids< Elem >(), libMesh::MeshTools::libmesh_assert_valid_procids< Node >(), libMesh::SerialMesh::local_elements_begin(), libMesh::ParallelMesh::local_elements_begin(), libMesh::SerialMesh::local_elements_end(), libMesh::ParallelMesh::local_elements_end(), libMesh::SerialMesh::local_level_elements_begin(), libMesh::ParallelMesh::local_level_elements_begin(), libMesh::SerialMesh::local_level_elements_end(), libMesh::ParallelMesh::local_level_elements_end(), libMesh::SerialMesh::local_nodes_begin(), libMesh::ParallelMesh::local_nodes_begin(), libMesh::SerialMesh::local_nodes_end(), libMesh::ParallelMesh::local_nodes_end(), libMesh::SerialMesh::local_not_level_elements_begin(), libMesh::ParallelMesh::local_not_level_elements_begin(), libMesh::SerialMesh::local_not_level_elements_end(), libMesh::ParallelMesh::local_not_level_elements_end(), libMesh::MeshRefinement::make_coarsening_compatible(), libMesh::MeshBase::n_active_local_elem(), libMesh::BoundaryInfo::n_boundary_conds(), libMesh::BoundaryInfo::n_edge_conds(), libMesh::DofMap::n_local_dofs(), libMesh::System::n_local_dofs(), libMesh::MeshBase::n_local_elem(), libMesh::MeshBase::n_local_nodes(), libMesh::BoundaryInfo::n_nodeset_conds(), libMesh::SerialMesh::not_local_elements_begin(), libMesh::ParallelMesh::not_local_elements_begin(), libMesh::SerialMesh::not_local_elements_end(), libMesh::ParallelMesh::not_local_elements_end(), libMesh::WeightedPatchRecoveryErrorEstimator::EstimateError::operator()(), libMesh::SparsityPattern::Build::operator()(), libMesh::PatchRecoveryErrorEstimator::EstimateError::operator()(), libMesh::MeshFunction::operator()(), libMesh::ParallelMesh::ParallelMesh(), libMesh::System::point_gradient(), libMesh::System::point_hessian(), libMesh::System::point_value(), libMesh::System::project_vector(), libMesh::Nemesis_IO_Helper::put_cmap_params(), libMesh::Nemesis_IO_Helper::put_elem_cmap(), libMesh::Nemesis_IO_Helper::put_elem_map(), libMesh::Nemesis_IO_Helper::put_loadbal_param(), libMesh::Nemesis_IO_Helper::put_node_cmap(), libMesh::Nemesis_IO_Helper::put_node_map(), libMesh::Nemesis_IO::read(), libMesh::CheckpointIO::read(), libMesh::XdrIO::read(), libMesh::UnstructuredMesh::read(), libMesh::CheckpointIO::read_connectivity(), libMesh::ExodusII_IO_Helper::read_elem_num_map(), libMesh::System::read_header(), libMesh::System::read_legacy_data(), libMesh::ExodusII_IO_Helper::read_node_num_map(), libMesh::System::read_parallel_data(), libMesh::System::read_SCALAR_dofs(), libMesh::XdrIO::read_serialized_bc_names(), libMesh::XdrIO::read_serialized_bcs(), libMesh::System::read_serialized_blocked_dof_objects(), libMesh::XdrIO::read_serialized_connectivity(), libMesh::System::read_serialized_data(), libMesh::XdrIO::read_serialized_nodes(), libMesh::XdrIO::read_serialized_nodesets(), libMesh::XdrIO::read_serialized_subdomain_names(), libMesh::System::read_serialized_vector(), libMesh::System::read_serialized_vectors(), libMesh::MeshData::read_xdr(), libMesh::Partitioner::set_node_processor_ids(), libMesh::DofMap::set_nonlocal_dof_objects(), libMesh::LaplaceMeshSmoother::smooth(), libMesh::BoundaryInfo::sync(), libMesh::MeshTools::total_weight(), libMesh::ParallelMesh::update_parallel_id_counts(), libMesh::MeshTools::weight(), libMesh::ExodusII_IO::write(), libMesh::CheckpointIO::write(), libMesh::XdrIO::write(), libMesh::UnstructuredMesh::write(), libMesh::EquationSystems::write(), libMesh::GMVIO::write_discontinuous_gmv(), libMesh::ExodusII_IO::write_element_data(), libMesh::ExodusII_IO_Helper::write_element_values(), libMesh::ExodusII_IO_Helper::write_elements(), libMesh::ExodusII_IO_Helper::write_elements_discontinuous(), libMesh::ExodusII_IO::write_global_data(), libMesh::ExodusII_IO_Helper::write_global_values(), libMesh::System::write_header(), libMesh::ExodusII_IO::write_information_records(), libMesh::ExodusII_IO_Helper::write_information_records(), libMesh::ExodusII_IO_Helper::write_nodal_coordinates(), libMesh::ExodusII_IO_Helper::write_nodal_coordinates_discontinuous(), libMesh::UCDIO::write_nodal_data(), libMesh::ExodusII_IO::write_nodal_data(), libMesh::ExodusII_IO::write_nodal_data_discontinuous(), libMesh::ExodusII_IO_Helper::write_nodal_values(), libMesh::ExodusII_IO_Helper::write_nodesets(), libMesh::Nemesis_IO_Helper::write_nodesets(), libMesh::System::write_parallel_data(), libMesh::System::write_SCALAR_dofs(), libMesh::XdrIO::write_serialized_bc_names(), libMesh::XdrIO::write_serialized_bcs(), libMesh::System::write_serialized_blocked_dof_objects(), libMesh::XdrIO::write_serialized_connectivity(), libMesh::System::write_serialized_data(), libMesh::XdrIO::write_serialized_nodes(), libMesh::XdrIO::write_serialized_nodesets(), libMesh::XdrIO::write_serialized_subdomain_names(), libMesh::System::write_serialized_vector(), libMesh::System::write_serialized_vectors(), libMesh::ExodusII_IO_Helper::write_sidesets(), libMesh::Nemesis_IO_Helper::write_sidesets(), libMesh::ExodusII_IO::write_timestep(), and libMesh::ExodusII_IO_Helper::write_timestep().

99  { return libmesh_cast_int<processor_id_type>(_communicator.rank()); }
template<typename T >
void libMesh::PetscLinearSolver< T >::restrict_solve_to ( const std::vector< unsigned int > *const  dofs,
const SubsetSolveMode  subset_solve_mode = SUBSET_ZERO 
)
virtual

After calling this method, all successive solves will be restricted to the given set of dofs, which must contain local dofs on each processor only and not contain any duplicates. This mode can be disabled by calling this method with dofs being a NULL pointer.

Reimplemented from libMesh::LinearSolver< T >.

Definition at line 397 of file petsc_linear_solver.C.

References libMesh::comm, libMesh::ierr, and PETSC_OWN_POINTER.

399 {
401 
402  /* The preconditioner (in particular if a default preconditioner)
403  will have to be reset. We call this->clear() to do that. This
404  call will also remove and free any previous subset that may have
405  been set before. */
406  this->clear();
407 
408  _subset_solve_mode = subset_solve_mode;
409 
410  if(dofs!=NULL)
411  {
412  PetscInt* petsc_dofs = NULL;
413  ierr = PetscMalloc(dofs->size()*sizeof(PetscInt), &petsc_dofs);
414  LIBMESH_CHKERRABORT(ierr);
415 
416  for(size_t i=0; i<dofs->size(); i++)
417  {
418  petsc_dofs[i] = (*dofs)[i];
419  }
420 
421  ierr = ISCreateLibMesh(this->comm().get(),dofs->size(),petsc_dofs,PETSC_OWN_POINTER,&_restrict_solve_to_is);
422  LIBMESH_CHKERRABORT(ierr);
423  }
424 }
template<typename T >
void libMesh::LinearSolver< T >::reuse_preconditioner ( bool  reuse_flag)
virtualinherited

Definition at line 130 of file linear_solver.C.

Referenced by libMesh::ImplicitSystem::disable_cache().

131  {
132  same_preconditioner = reuse_flag;
133  }
template<typename T >
void libMesh::PetscLinearSolver< T >::set_petsc_solver_type ( )
private

Tells PETSC to use the user-specified solver stored in _solver_type

Definition at line 1752 of file petsc_linear_solver.C.

References libMeshEnums::BICG, libMeshEnums::BICGSTAB, libMeshEnums::CG, libMeshEnums::CGS, libMeshEnums::CHEBYSHEV, libMeshEnums::CR, libMesh::Utility::enum_to_string(), libMesh::err, libMeshEnums::GMRES, libMesh::ierr, libMeshEnums::LSQR, libMeshEnums::MINRES, libMeshEnums::RICHARDSON, libMeshEnums::TCQMR, and libMeshEnums::TFQMR.

1753 {
1754  PetscErrorCode ierr = 0;
1755 
1756  switch (this->_solver_type)
1757  {
1758 
1759  case CG:
1760  ierr = KSPSetType (_ksp, (char*) KSPCG); LIBMESH_CHKERRABORT(ierr); return;
1761 
1762  case CR:
1763  ierr = KSPSetType (_ksp, (char*) KSPCR); LIBMESH_CHKERRABORT(ierr); return;
1764 
1765  case CGS:
1766  ierr = KSPSetType (_ksp, (char*) KSPCGS); LIBMESH_CHKERRABORT(ierr); return;
1767 
1768  case BICG:
1769  ierr = KSPSetType (_ksp, (char*) KSPBICG); LIBMESH_CHKERRABORT(ierr); return;
1770 
1771  case TCQMR:
1772  ierr = KSPSetType (_ksp, (char*) KSPTCQMR); LIBMESH_CHKERRABORT(ierr); return;
1773 
1774  case TFQMR:
1775  ierr = KSPSetType (_ksp, (char*) KSPTFQMR); LIBMESH_CHKERRABORT(ierr); return;
1776 
1777  case LSQR:
1778  ierr = KSPSetType (_ksp, (char*) KSPLSQR); LIBMESH_CHKERRABORT(ierr); return;
1779 
1780  case BICGSTAB:
1781  ierr = KSPSetType (_ksp, (char*) KSPBCGS); LIBMESH_CHKERRABORT(ierr); return;
1782 
1783  case MINRES:
1784  ierr = KSPSetType (_ksp, (char*) KSPMINRES); LIBMESH_CHKERRABORT(ierr); return;
1785 
1786  case GMRES:
1787  ierr = KSPSetType (_ksp, (char*) KSPGMRES); LIBMESH_CHKERRABORT(ierr); return;
1788 
1789  case RICHARDSON:
1790  ierr = KSPSetType (_ksp, (char*) KSPRICHARDSON); LIBMESH_CHKERRABORT(ierr); return;
1791 
1792  case CHEBYSHEV:
1793 #if defined(LIBMESH_HAVE_PETSC) && PETSC_VERSION_LESS_THAN(3,3,0)
1794  ierr = KSPSetType (_ksp, (char*) KSPCHEBYCHEV); LIBMESH_CHKERRABORT(ierr); return;
1795 #else
1796  ierr = KSPSetType (_ksp, (char*) KSPCHEBYSHEV); LIBMESH_CHKERRABORT(ierr); return;
1797 #endif
1798 
1799 
1800  default:
1801  libMesh::err << "ERROR: Unsupported PETSC Solver: "
1802  << Utility::enum_to_string(this->_solver_type) << std::endl
1803  << "Continuing with PETSC defaults" << std::endl;
1804  }
1805 }
template<typename T >
void libMesh::LinearSolver< T >::set_preconditioner_type ( const PreconditionerType  pct)
inherited

Sets the type of preconditioner to use.

Definition at line 106 of file linear_solver.C.

107 {
108  if(_preconditioner)
109  _preconditioner->set_type(pct);
110  else
111  _preconditioner_type = pct;
112 }
template<typename T>
void libMesh::LinearSolver< T >::set_solver_type ( const SolverType  st)
inlineinherited

Sets the type of solver to use.

Definition at line 105 of file linear_solver.h.

106  { _solver_type = st; }
template<typename T >
std::pair<unsigned int, Real> libMesh::PetscLinearSolver< T >::solve ( SparseMatrix< T > &  matrix_in,
NumericVector< T > &  solution_in,
NumericVector< T > &  rhs_in,
const double  tol,
const unsigned int  m_its 
)
inlinevirtual

Call the Petsc solver. It calls the method below, using the same matrix for the system and preconditioner matrices.

Implements libMesh::LinearSolver< T >.

Definition at line 142 of file petsc_linear_solver.h.

147  {
148  return this->solve(matrix_in, matrix_in, solution_in, rhs_in, tol, m_its);
149  }
template<typename T>
std::pair< unsigned int, Real > libMesh::LinearSolver< T >::solve ( SparseMatrix< T > &  matrix,
SparseMatrix< T > *  precond_matrix,
NumericVector< T > &  sol,
NumericVector< T > &  rhs,
const double  tol,
const unsigned int  n_iter 
)
inlineinherited

This function calls the solver "_solver_type" preconditioned with the "_preconditioner_type" preconditioner. The preconditioning matrix is used if it is provided, or the system matrix is used if precond_matrix is null

Definition at line 298 of file linear_solver.h.

304 {
305  if (pc_mat)
306  return this->solve(mat, *pc_mat, sol, rhs, tol, n_iter);
307  else
308  return this->solve(mat, sol, rhs, tol, n_iter);
309 }
template<typename T >
std::pair< unsigned int, Real > libMesh::PetscLinearSolver< T >::solve ( SparseMatrix< T > &  matrix,
SparseMatrix< T > &  preconditioner,
NumericVector< T > &  solution,
NumericVector< T > &  rhs,
const double  tol,
const unsigned int  m_its 
)
virtual

This method allows you to call a linear solver while specifying the matrix to use as the (left) preconditioning matrix. Note that the linear solver will not compute a preconditioner in this case, and will instead premultiply by the matrix you provide.

In PETSc, this is accomplished by calling

PCSetType(_pc, PCMAT);

before invoking KSPSolve(). Note: this functionality is not implemented in the LinearSolver class since there is not a built-in analog to this method for LasPack – You could probably implement it by hand if you wanted.

Implements libMesh::LinearSolver< T >.

Definition at line 430 of file petsc_linear_solver.C.

References libMesh::PetscMatrix< T >::close(), libMesh::comm, libMesh::ierr, libMesh::TriangleWrapper::init(), libMesh::PetscMatrix< T >::init(), libMesh::NumericVector< T >::local_size(), libMesh::PetscMatrix< T >::mat(), libMesh::START_LOG(), libMesh::STOP_LOG(), libMeshEnums::SUBSET_COPY_RHS, libMeshEnums::SUBSET_DONT_TOUCH, and libMeshEnums::SUBSET_ZERO.

436 {
437  START_LOG("solve()", "PetscLinearSolver");
438 
439  // Make sure the data passed in are really of Petsc types
440  PetscMatrix<T>* matrix = libmesh_cast_ptr<PetscMatrix<T>*>(&matrix_in);
441  PetscMatrix<T>* precond = libmesh_cast_ptr<PetscMatrix<T>*>(&precond_in);
442  PetscVector<T>* solution = libmesh_cast_ptr<PetscVector<T>*>(&solution_in);
443  PetscVector<T>* rhs = libmesh_cast_ptr<PetscVector<T>*>(&rhs_in);
444 
445  this->init (matrix);
446 
448  PetscInt its=0, max_its = static_cast<PetscInt>(m_its);
449  PetscReal final_resid=0.;
450 
451  // Close the matrices and vectors in case this wasn't already done.
452  matrix->close ();
453  precond->close ();
454  solution->close ();
455  rhs->close ();
456 
457 // // If matrix != precond, then this means we have specified a
458 // // special preconditioner, so reset preconditioner type to PCMAT.
459 // if (matrix != precond)
460 // {
461 // this->_preconditioner_type = USER_PRECOND;
462 // this->set_petsc_preconditioner_type ();
463 // }
464 
465 // 2.1.x & earlier style
466 #if PETSC_VERSION_LESS_THAN(2,2,0)
467 
468  if(_restrict_solve_to_is!=NULL)
469  {
470  libmesh_not_implemented();
471  }
472 
473  // Set operators. The input matrix works as the preconditioning matrix
474  ierr = SLESSetOperators(_sles, matrix->mat(), precond->mat(),
475  DIFFERENT_NONZERO_PATTERN);
476  LIBMESH_CHKERRABORT(ierr);
477 
478  // Set the tolerances for the iterative solver. Use the user-supplied
479  // tolerance for the relative residual & leave the others at default values.
480  ierr = KSPSetTolerances (_ksp, tol, PETSC_DEFAULT,
481  PETSC_DEFAULT, max_its);
482  LIBMESH_CHKERRABORT(ierr);
483 
484 
485  // Solve the linear system
486  ierr = SLESSolve (_sles, rhs->vec(), solution->vec(), &its);
487  LIBMESH_CHKERRABORT(ierr);
488 
489 
490  // Get the norm of the final residual to return to the user.
491  ierr = KSPGetResidualNorm (_ksp, &final_resid);
492  LIBMESH_CHKERRABORT(ierr);
493 
494 // 2.2.0
495 #elif PETSC_VERSION_LESS_THAN(2,2,1)
496 
497  if(_restrict_solve_to_is!=NULL)
498  {
499  libmesh_not_implemented();
500  }
501 
502  // Set operators. The input matrix works as the preconditioning matrix
503  //ierr = KSPSetOperators(_ksp, matrix->mat(), precond->mat(),
504  // SAME_NONZERO_PATTERN);
505  // LIBMESH_CHKERRABORT(ierr);
506 
507 
508  // Set the tolerances for the iterative solver. Use the user-supplied
509  // tolerance for the relative residual & leave the others at default values.
510  // Convergence is detected at iteration k if
511  // ||r_k||_2 < max(rtol*||b||_2 , abstol)
512  // where r_k is the residual vector and b is the right-hand side. Note that
513  // it is the *maximum* of the two values, the larger of which will almost
514  // always be rtol*||b||_2.
515  ierr = KSPSetTolerances (_ksp,
516  tol, // rtol = relative decrease in residual (1.e-5)
517  PETSC_DEFAULT, // abstol = absolute convergence tolerance (1.e-50)
518  PETSC_DEFAULT, // dtol = divergence tolerance (1.e+5)
519  max_its);
520  LIBMESH_CHKERRABORT(ierr);
521 
522 
523  // Set the solution vector to use
524  ierr = KSPSetSolution (_ksp, solution->vec());
525  LIBMESH_CHKERRABORT(ierr);
526 
527  // Set the RHS vector to use
528  ierr = KSPSetRhs (_ksp, rhs->vec());
529  LIBMESH_CHKERRABORT(ierr);
530 
531  // Solve the linear system
532  ierr = KSPSolve (_ksp);
533  LIBMESH_CHKERRABORT(ierr);
534 
535  // Get the number of iterations required for convergence
536  ierr = KSPGetIterationNumber (_ksp, &its);
537  LIBMESH_CHKERRABORT(ierr);
538 
539  // Get the norm of the final residual to return to the user.
540  ierr = KSPGetResidualNorm (_ksp, &final_resid);
541  LIBMESH_CHKERRABORT(ierr);
542 
543 // 2.2.1 & newer style
544 #else
545 
546  Mat submat = NULL;
547  Mat subprecond = NULL;
548  Vec subrhs = NULL;
549  Vec subsolution = NULL;
550  VecScatter scatter = NULL;
551  PetscMatrix<Number>* subprecond_matrix = NULL;
552 
553  // Set operators. Also restrict rhs and solution vector to
554  // subdomain if neccessary.
555  if(_restrict_solve_to_is!=NULL)
556  {
557  size_t is_local_size = this->_restrict_solve_to_is_local_size();
558 
559  ierr = VecCreate(this->comm().get(),&subrhs);
560  LIBMESH_CHKERRABORT(ierr);
561  ierr = VecSetSizes(subrhs,is_local_size,PETSC_DECIDE);
562  LIBMESH_CHKERRABORT(ierr);
563  ierr = VecSetFromOptions(subrhs);
564  LIBMESH_CHKERRABORT(ierr);
565 
566  ierr = VecCreate(this->comm().get(),&subsolution);
567  LIBMESH_CHKERRABORT(ierr);
568  ierr = VecSetSizes(subsolution,is_local_size,PETSC_DECIDE);
569  LIBMESH_CHKERRABORT(ierr);
570  ierr = VecSetFromOptions(subsolution);
571  LIBMESH_CHKERRABORT(ierr);
572 
573  ierr = VecScatterCreate(rhs->vec(),_restrict_solve_to_is, subrhs,NULL, &scatter);
574  LIBMESH_CHKERRABORT(ierr);
575 
576  ierr = VecScatterBegin(scatter,rhs->vec(),subrhs,INSERT_VALUES,SCATTER_FORWARD);
577  LIBMESH_CHKERRABORT(ierr);
578  ierr = VecScatterEnd(scatter,rhs->vec(),subrhs,INSERT_VALUES,SCATTER_FORWARD);
579  LIBMESH_CHKERRABORT(ierr);
580 
581  ierr = VecScatterBegin(scatter,solution->vec(),subsolution,INSERT_VALUES,SCATTER_FORWARD);
582  LIBMESH_CHKERRABORT(ierr);
583  ierr = VecScatterEnd(scatter,solution->vec(),subsolution,INSERT_VALUES,SCATTER_FORWARD);
584  LIBMESH_CHKERRABORT(ierr);
585 
586 #if PETSC_VERSION_LESS_THAN(3,1,0)
587  ierr = MatGetSubMatrix(matrix->mat(),
589  PETSC_DECIDE,MAT_INITIAL_MATRIX,&submat);
590  LIBMESH_CHKERRABORT(ierr);
591  ierr = MatGetSubMatrix(precond->mat(),
593  PETSC_DECIDE,MAT_INITIAL_MATRIX,&subprecond);
594  LIBMESH_CHKERRABORT(ierr);
595 #else
596  ierr = MatGetSubMatrix(matrix->mat(),
598  MAT_INITIAL_MATRIX,&submat);
599  LIBMESH_CHKERRABORT(ierr);
600  ierr = MatGetSubMatrix(precond->mat(),
602  MAT_INITIAL_MATRIX,&subprecond);
603  LIBMESH_CHKERRABORT(ierr);
604 #endif
605 
606  /* Since removing columns of the matrix changes the equation
607  system, we will now change the right hand side to compensate
608  for this. Note that this is not necessary if \p SUBSET_ZERO
609  has been selected. */
611  {
612  _create_complement_is(rhs_in);
613  size_t is_complement_local_size = rhs_in.local_size()-is_local_size;
614 
615  Vec subvec1 = NULL;
616  Mat submat1 = NULL;
617  VecScatter scatter1 = NULL;
618 
619  ierr = VecCreate(this->comm().get(),&subvec1);
620  LIBMESH_CHKERRABORT(ierr);
621  ierr = VecSetSizes(subvec1,is_complement_local_size,PETSC_DECIDE);
622  LIBMESH_CHKERRABORT(ierr);
623  ierr = VecSetFromOptions(subvec1);
624  LIBMESH_CHKERRABORT(ierr);
625 
626  ierr = VecScatterCreate(rhs->vec(),_restrict_solve_to_is_complement, subvec1,NULL, &scatter1);
627  LIBMESH_CHKERRABORT(ierr);
628 
629  ierr = VecScatterBegin(scatter1,_subset_solve_mode==SUBSET_COPY_RHS ? rhs->vec() : solution->vec(),subvec1,INSERT_VALUES,SCATTER_FORWARD);
630  LIBMESH_CHKERRABORT(ierr);
631  ierr = VecScatterEnd(scatter1,_subset_solve_mode==SUBSET_COPY_RHS ? rhs->vec() : solution->vec(),subvec1,INSERT_VALUES,SCATTER_FORWARD);
632  LIBMESH_CHKERRABORT(ierr);
633 
634  ierr = VecScale(subvec1,-1.0);
635  LIBMESH_CHKERRABORT(ierr);
636 
637 #if PETSC_VERSION_LESS_THAN(3,1,0)
638  ierr = MatGetSubMatrix(matrix->mat(),
640  PETSC_DECIDE,MAT_INITIAL_MATRIX,&submat1);
641  LIBMESH_CHKERRABORT(ierr);
642 #else
643  ierr = MatGetSubMatrix(matrix->mat(),
645  MAT_INITIAL_MATRIX,&submat1);
646  LIBMESH_CHKERRABORT(ierr);
647 #endif
648 
649  ierr = MatMultAdd(submat1,subvec1,subrhs,subrhs);
650  LIBMESH_CHKERRABORT(ierr);
651 
652  ierr = LibMeshVecScatterDestroy(&scatter1);
653  LIBMESH_CHKERRABORT(ierr);
654  ierr = LibMeshVecDestroy(&subvec1);
655  LIBMESH_CHKERRABORT(ierr);
656  ierr = LibMeshMatDestroy(&submat1);
657  LIBMESH_CHKERRABORT(ierr);
658  }
659 
660  ierr = KSPSetOperators(_ksp, submat, subprecond,
661  this->same_preconditioner ? SAME_PRECONDITIONER : DIFFERENT_NONZERO_PATTERN);
662  LIBMESH_CHKERRABORT(ierr);
663 
664  if(this->_preconditioner)
665  {
666  subprecond_matrix = new PetscMatrix<Number>(subprecond,
667  this->comm());
668  this->_preconditioner->set_matrix(*subprecond_matrix);
669  this->_preconditioner->init();
670  }
671  }
672  else
673  {
674  ierr = KSPSetOperators(_ksp, matrix->mat(), precond->mat(),
675  this->same_preconditioner ? SAME_PRECONDITIONER : DIFFERENT_NONZERO_PATTERN);
676  LIBMESH_CHKERRABORT(ierr);
677 
678  if(this->_preconditioner)
679  {
680  this->_preconditioner->set_matrix(matrix_in);
681  this->_preconditioner->init();
682  }
683  }
684 
685  // Set the tolerances for the iterative solver. Use the user-supplied
686  // tolerance for the relative residual & leave the others at default values.
687  ierr = KSPSetTolerances (_ksp, tol, PETSC_DEFAULT,
688  PETSC_DEFAULT, max_its);
689  LIBMESH_CHKERRABORT(ierr);
690 
691  // Solve the linear system
692  if(_restrict_solve_to_is!=NULL)
693  {
694  ierr = KSPSolve (_ksp, subrhs, subsolution);
695  LIBMESH_CHKERRABORT(ierr);
696  }
697  else
698  {
699  ierr = KSPSolve (_ksp, rhs->vec(), solution->vec());
700  LIBMESH_CHKERRABORT(ierr);
701  }
702 
703  // Get the number of iterations required for convergence
704  ierr = KSPGetIterationNumber (_ksp, &its);
705  LIBMESH_CHKERRABORT(ierr);
706 
707  // Get the norm of the final residual to return to the user.
708  ierr = KSPGetResidualNorm (_ksp, &final_resid);
709  LIBMESH_CHKERRABORT(ierr);
710 
711  if(_restrict_solve_to_is!=NULL)
712  {
713  switch(_subset_solve_mode)
714  {
715  case SUBSET_ZERO:
716  ierr = VecZeroEntries(solution->vec());
717  LIBMESH_CHKERRABORT(ierr);
718  break;
719 
720  case SUBSET_COPY_RHS:
721  ierr = VecCopy(rhs->vec(),solution->vec());
722  LIBMESH_CHKERRABORT(ierr);
723  break;
724 
725  case SUBSET_DONT_TOUCH:
726  /* Nothing to do here. */
727  break;
728 
729  }
730  ierr = VecScatterBegin(scatter,subsolution,solution->vec(),INSERT_VALUES,SCATTER_REVERSE);
731  LIBMESH_CHKERRABORT(ierr);
732  ierr = VecScatterEnd(scatter,subsolution,solution->vec(),INSERT_VALUES,SCATTER_REVERSE);
733  LIBMESH_CHKERRABORT(ierr);
734 
735  ierr = LibMeshVecScatterDestroy(&scatter);
736  LIBMESH_CHKERRABORT(ierr);
737 
738  if(this->_preconditioner)
739  {
740  /* Before we delete subprecond_matrix, we should give the
741  _preconditioner a different matrix. */
742  this->_preconditioner->set_matrix(matrix_in);
743  this->_preconditioner->init();
744  delete subprecond_matrix;
745  subprecond_matrix = NULL;
746  }
747 
748  ierr = LibMeshVecDestroy(&subsolution);
749  LIBMESH_CHKERRABORT(ierr);
750  ierr = LibMeshVecDestroy(&subrhs);
751  LIBMESH_CHKERRABORT(ierr);
752  ierr = LibMeshMatDestroy(&submat);
753  LIBMESH_CHKERRABORT(ierr);
754  ierr = LibMeshMatDestroy(&subprecond);
755  LIBMESH_CHKERRABORT(ierr);
756  }
757 
758 #endif
759 
760  STOP_LOG("solve()", "PetscLinearSolver");
761  // return the # of its. and the final residual norm.
762  return std::make_pair(its, final_resid);
763 }
template<typename T >
std::pair< unsigned int, Real > libMesh::PetscLinearSolver< T >::solve ( const ShellMatrix< T > &  shell_matrix,
NumericVector< T > &  solution_in,
NumericVector< T > &  rhs_in,
const double  tol,
const unsigned int  m_its 
)
virtual

This function solves a system whose matrix is a shell matrix.

Implements libMesh::LinearSolver< T >.

Definition at line 1106 of file petsc_linear_solver.C.

References libMesh::PetscVector< T >::close(), libMesh::comm, libMesh::ierr, libMesh::TriangleWrapper::init(), libMesh::libmesh_assert(), libMesh::NumericVector< T >::local_size(), libMesh::out, libMesh::NumericVector< T >::size(), libMesh::START_LOG(), libMesh::STOP_LOG(), libMeshEnums::SUBSET_COPY_RHS, libMeshEnums::SUBSET_DONT_TOUCH, libMeshEnums::SUBSET_ZERO, and libMesh::PetscVector< T >::vec().

1111 {
1112 
1113 #if PETSC_VERSION_LESS_THAN(2,3,1)
1114  // FIXME[JWP]: There will be a bunch of unused variable warnings
1115  // for older PETScs here.
1116  libMesh::out << "This method has been developed with PETSc 2.3.1. "
1117  << "No one has made it backwards compatible with older "
1118  << "versions of PETSc so far; however, it might work "
1119  << "without any change with some older version." << std::endl;
1120  libmesh_error();
1121  return std::make_pair(0,0.0);
1122 
1123 #else
1124 
1125 #if PETSC_VERSION_LESS_THAN(3,1,0)
1126  if(_restrict_solve_to_is!=NULL)
1127  {
1128  libMesh::out << "The current implementation of subset solves with "
1129  << "shell matrices requires PETSc version 3.1 or above. "
1130  << "Older PETSc version do not support automatic "
1131  << "submatrix generation of shell matrices."
1132  << std::endl;
1133  libmesh_error();
1134  }
1135 #endif
1136 
1137  START_LOG("solve()", "PetscLinearSolver");
1138 
1139  // Make sure the data passed in are really of Petsc types
1140  PetscVector<T>* solution = libmesh_cast_ptr<PetscVector<T>*>(&solution_in);
1141  PetscVector<T>* rhs = libmesh_cast_ptr<PetscVector<T>*>(&rhs_in);
1142 
1143  this->init ();
1144 
1145  PetscErrorCode ierr=0;
1146  PetscInt its=0, max_its = static_cast<PetscInt>(m_its);
1147  PetscReal final_resid=0.;
1148 
1149  Mat submat = NULL;
1150  Vec subrhs = NULL;
1151  Vec subsolution = NULL;
1152  VecScatter scatter = NULL;
1153 
1154  // Close the matrices and vectors in case this wasn't already done.
1155  solution->close ();
1156  rhs->close ();
1157 
1158  // Prepare the matrix.
1159  Mat mat;
1160  ierr = MatCreateShell(this->comm().get(),
1161  rhs_in.local_size(),
1162  solution_in.local_size(),
1163  rhs_in.size(),
1164  solution_in.size(),
1165  const_cast<void*>(static_cast<const void*>(&shell_matrix)),
1166  &mat);
1167  /* Note that the const_cast above is only necessary because PETSc
1168  does not accept a const void*. Inside the member function
1169  _petsc_shell_matrix() below, the pointer is casted back to a
1170  const ShellMatrix<T>*. */
1171 
1172  LIBMESH_CHKERRABORT(ierr);
1173  ierr = MatShellSetOperation(mat,MATOP_MULT,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_mult));
1174  LIBMESH_CHKERRABORT(ierr);
1175  ierr = MatShellSetOperation(mat,MATOP_MULT_ADD,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_mult_add));
1176  LIBMESH_CHKERRABORT(ierr);
1177  ierr = MatShellSetOperation(mat,MATOP_GET_DIAGONAL,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_get_diagonal));
1178  LIBMESH_CHKERRABORT(ierr);
1179 
1180  // Restrict rhs and solution vectors and set operators. The input
1181  // matrix works as the preconditioning matrix.
1182  if(_restrict_solve_to_is!=NULL)
1183  {
1184  size_t is_local_size = this->_restrict_solve_to_is_local_size();
1185 
1186  ierr = VecCreate(this->comm().get(),&subrhs);
1187  LIBMESH_CHKERRABORT(ierr);
1188  ierr = VecSetSizes(subrhs,is_local_size,PETSC_DECIDE);
1189  LIBMESH_CHKERRABORT(ierr);
1190  ierr = VecSetFromOptions(subrhs);
1191  LIBMESH_CHKERRABORT(ierr);
1192 
1193  ierr = VecCreate(this->comm().get(),&subsolution);
1194  LIBMESH_CHKERRABORT(ierr);
1195  ierr = VecSetSizes(subsolution,is_local_size,PETSC_DECIDE);
1196  LIBMESH_CHKERRABORT(ierr);
1197  ierr = VecSetFromOptions(subsolution);
1198  LIBMESH_CHKERRABORT(ierr);
1199 
1200  ierr = VecScatterCreate(rhs->vec(),_restrict_solve_to_is, subrhs,NULL, &scatter);
1201  LIBMESH_CHKERRABORT(ierr);
1202 
1203  ierr = VecScatterBegin(scatter,rhs->vec(),subrhs,INSERT_VALUES,SCATTER_FORWARD);
1204  LIBMESH_CHKERRABORT(ierr);
1205  ierr = VecScatterEnd(scatter,rhs->vec(),subrhs,INSERT_VALUES,SCATTER_FORWARD);
1206  LIBMESH_CHKERRABORT(ierr);
1207 
1208  ierr = VecScatterBegin(scatter,solution->vec(),subsolution,INSERT_VALUES,SCATTER_FORWARD);
1209  LIBMESH_CHKERRABORT(ierr);
1210  ierr = VecScatterEnd(scatter,solution->vec(),subsolution,INSERT_VALUES,SCATTER_FORWARD);
1211  LIBMESH_CHKERRABORT(ierr);
1212 
1213 #if PETSC_VERSION_LESS_THAN(3,1,0)
1214  /* This point can't be reached, see above. */
1215  libmesh_assert(false);
1216 #else
1217  ierr = MatGetSubMatrix(mat,
1219  MAT_INITIAL_MATRIX,&submat);
1220  LIBMESH_CHKERRABORT(ierr);
1221 #endif
1222 
1223  /* Since removing columns of the matrix changes the equation
1224  system, we will now change the right hand side to compensate
1225  for this. Note that this is not necessary if \p SUBSET_ZERO
1226  has been selected. */
1228  {
1229  _create_complement_is(rhs_in);
1230  size_t is_complement_local_size = rhs_in.local_size()-is_local_size;
1231 
1232  Vec subvec1 = NULL;
1233  Mat submat1 = NULL;
1234  VecScatter scatter1 = NULL;
1235 
1236  ierr = VecCreate(this->comm().get(),&subvec1);
1237  LIBMESH_CHKERRABORT(ierr);
1238  ierr = VecSetSizes(subvec1,is_complement_local_size,PETSC_DECIDE);
1239  LIBMESH_CHKERRABORT(ierr);
1240  ierr = VecSetFromOptions(subvec1);
1241  LIBMESH_CHKERRABORT(ierr);
1242 
1243  ierr = VecScatterCreate(rhs->vec(),_restrict_solve_to_is_complement, subvec1,NULL, &scatter1);
1244  LIBMESH_CHKERRABORT(ierr);
1245 
1246  ierr = VecScatterBegin(scatter1,_subset_solve_mode==SUBSET_COPY_RHS ? rhs->vec() : solution->vec(),subvec1,INSERT_VALUES,SCATTER_FORWARD);
1247  LIBMESH_CHKERRABORT(ierr);
1248  ierr = VecScatterEnd(scatter1,_subset_solve_mode==SUBSET_COPY_RHS ? rhs->vec() : solution->vec(),subvec1,INSERT_VALUES,SCATTER_FORWARD);
1249  LIBMESH_CHKERRABORT(ierr);
1250 
1251  ierr = VecScale(subvec1,-1.0);
1252  LIBMESH_CHKERRABORT(ierr);
1253 
1254 #if PETSC_VERSION_LESS_THAN(3,1,0)
1255  /* This point can't be reached, see above. */
1256  libmesh_assert(false);
1257 #else
1258  ierr = MatGetSubMatrix(mat,
1260  MAT_INITIAL_MATRIX,&submat1);
1261  LIBMESH_CHKERRABORT(ierr);
1262 #endif
1263 
1264  // The following lines would be correct, but don't work
1265  // correctly in PETSc up to 3.1.0-p5. See discussion in
1266  // petsc-users of Nov 9, 2010.
1267  //
1268  // ierr = MatMultAdd(submat1,subvec1,subrhs,subrhs);
1269  // LIBMESH_CHKERRABORT(ierr);
1270  //
1271  // We workaround by using a temporary vector. Note that the
1272  // fix in PETsc 3.1.0-p6 uses a temporary vector internally,
1273  // so this is no effective performance loss.
1274  Vec subvec2 = NULL;
1275  ierr = VecCreate(this->comm().get(),&subvec2);
1276  LIBMESH_CHKERRABORT(ierr);
1277  ierr = VecSetSizes(subvec2,is_local_size,PETSC_DECIDE);
1278  LIBMESH_CHKERRABORT(ierr);
1279  ierr = VecSetFromOptions(subvec2);
1280  LIBMESH_CHKERRABORT(ierr);
1281  ierr = MatMult(submat1,subvec1,subvec2);
1282  LIBMESH_CHKERRABORT(ierr);
1283  ierr = VecAXPY(subrhs,1.0,subvec2);
1284 
1285  ierr = LibMeshVecScatterDestroy(&scatter1);
1286  LIBMESH_CHKERRABORT(ierr);
1287  ierr = LibMeshVecDestroy(&subvec1);
1288  LIBMESH_CHKERRABORT(ierr);
1289  ierr = LibMeshMatDestroy(&submat1);
1290  LIBMESH_CHKERRABORT(ierr);
1291  }
1292 
1293  ierr = KSPSetOperators(_ksp, submat, submat,
1294  DIFFERENT_NONZERO_PATTERN);
1295  LIBMESH_CHKERRABORT(ierr);
1296  }
1297  else
1298  {
1299  ierr = KSPSetOperators(_ksp, mat, mat,
1300  DIFFERENT_NONZERO_PATTERN);
1301  LIBMESH_CHKERRABORT(ierr);
1302  }
1303 
1304  // Set the tolerances for the iterative solver. Use the user-supplied
1305  // tolerance for the relative residual & leave the others at default values.
1306  ierr = KSPSetTolerances (_ksp, tol, PETSC_DEFAULT,
1307  PETSC_DEFAULT, max_its);
1308  LIBMESH_CHKERRABORT(ierr);
1309 
1310  // Solve the linear system
1311  if(_restrict_solve_to_is!=NULL)
1312  {
1313  ierr = KSPSolve (_ksp, subrhs, subsolution);
1314  LIBMESH_CHKERRABORT(ierr);
1315  }
1316  else
1317  {
1318  ierr = KSPSolve (_ksp, rhs->vec(), solution->vec());
1319  LIBMESH_CHKERRABORT(ierr);
1320  }
1321 
1322  // Get the number of iterations required for convergence
1323  ierr = KSPGetIterationNumber (_ksp, &its);
1324  LIBMESH_CHKERRABORT(ierr);
1325 
1326  // Get the norm of the final residual to return to the user.
1327  ierr = KSPGetResidualNorm (_ksp, &final_resid);
1328  LIBMESH_CHKERRABORT(ierr);
1329 
1330  if(_restrict_solve_to_is!=NULL)
1331  {
1332  switch(_subset_solve_mode)
1333  {
1334  case SUBSET_ZERO:
1335  ierr = VecZeroEntries(solution->vec());
1336  LIBMESH_CHKERRABORT(ierr);
1337  break;
1338 
1339  case SUBSET_COPY_RHS:
1340  ierr = VecCopy(rhs->vec(),solution->vec());
1341  LIBMESH_CHKERRABORT(ierr);
1342  break;
1343 
1344  case SUBSET_DONT_TOUCH:
1345  /* Nothing to do here. */
1346  break;
1347 
1348  }
1349  ierr = VecScatterBegin(scatter,subsolution,solution->vec(),INSERT_VALUES,SCATTER_REVERSE);
1350  LIBMESH_CHKERRABORT(ierr);
1351  ierr = VecScatterEnd(scatter,subsolution,solution->vec(),INSERT_VALUES,SCATTER_REVERSE);
1352  LIBMESH_CHKERRABORT(ierr);
1353 
1354  ierr = LibMeshVecScatterDestroy(&scatter);
1355  LIBMESH_CHKERRABORT(ierr);
1356 
1357  ierr = LibMeshVecDestroy(&subsolution);
1358  LIBMESH_CHKERRABORT(ierr);
1359  ierr = LibMeshVecDestroy(&subrhs);
1360  LIBMESH_CHKERRABORT(ierr);
1361  ierr = LibMeshMatDestroy(&submat);
1362  LIBMESH_CHKERRABORT(ierr);
1363  }
1364 
1365  // Destroy the matrix.
1366  ierr = LibMeshMatDestroy(&mat);
1367  LIBMESH_CHKERRABORT(ierr);
1368 
1369  STOP_LOG("solve()", "PetscLinearSolver");
1370  // return the # of its. and the final residual norm.
1371  return std::make_pair(its, final_resid);
1372 
1373 #endif
1374 
1375 }
template<typename T >
std::pair< unsigned int, Real > libMesh::PetscLinearSolver< T >::solve ( const ShellMatrix< T > &  shell_matrix,
const SparseMatrix< T > &  precond_matrix,
NumericVector< T > &  solution_in,
NumericVector< T > &  rhs_in,
const double  tol,
const unsigned int  m_its 
)
virtual

This function solves a system whose matrix is a shell matrix, but a sparse matrix is used as preconditioning matrix, this allowing other preconditioners than JACOBI.

Implements libMesh::LinearSolver< T >.

Definition at line 1381 of file petsc_linear_solver.C.

References libMesh::PetscMatrix< T >::close(), libMesh::comm, libMesh::ierr, libMesh::TriangleWrapper::init(), libMesh::PetscMatrix< T >::init(), libMesh::libmesh_assert(), libMesh::NumericVector< T >::local_size(), libMesh::out, libMesh::NumericVector< T >::size(), libMesh::START_LOG(), libMesh::STOP_LOG(), libMeshEnums::SUBSET_COPY_RHS, libMeshEnums::SUBSET_DONT_TOUCH, and libMeshEnums::SUBSET_ZERO.

1387 {
1388 
1389 #if PETSC_VERSION_LESS_THAN(2,3,1)
1390  // FIXME[JWP]: There will be a bunch of unused variable warnings
1391  // for older PETScs here.
1392  libMesh::out << "This method has been developed with PETSc 2.3.1. "
1393  << "No one has made it backwards compatible with older "
1394  << "versions of PETSc so far; however, it might work "
1395  << "without any change with some older version." << std::endl;
1396  libmesh_error();
1397  return std::make_pair(0,0.0);
1398 
1399 #else
1400 
1401 #if PETSC_VERSION_LESS_THAN(3,1,0)
1402  if(_restrict_solve_to_is!=NULL)
1403  {
1404  libMesh::out << "The current implementation of subset solves with "
1405  << "shell matrices requires PETSc version 3.1 or above. "
1406  << "Older PETSc version do not support automatic "
1407  << "submatrix generation of shell matrices."
1408  << std::endl;
1409  libmesh_error();
1410  }
1411 #endif
1412 
1413  START_LOG("solve()", "PetscLinearSolver");
1414 
1415  // Make sure the data passed in are really of Petsc types
1416  const PetscMatrix<T>* precond = libmesh_cast_ptr<const PetscMatrix<T>*>(&precond_matrix);
1417  PetscVector<T>* solution = libmesh_cast_ptr<PetscVector<T>*>(&solution_in);
1418  PetscVector<T>* rhs = libmesh_cast_ptr<PetscVector<T>*>(&rhs_in);
1419 
1420  this->init ();
1421 
1422  PetscErrorCode ierr=0;
1423  PetscInt its=0, max_its = static_cast<PetscInt>(m_its);
1424  PetscReal final_resid=0.;
1425 
1426  Mat submat = NULL;
1427  Mat subprecond = NULL;
1428  Vec subrhs = NULL;
1429  Vec subsolution = NULL;
1430  VecScatter scatter = NULL;
1431  PetscMatrix<Number>* subprecond_matrix = NULL;
1432 
1433  // Close the matrices and vectors in case this wasn't already done.
1434  solution->close ();
1435  rhs->close ();
1436 
1437  // Prepare the matrix.
1438  Mat mat;
1439  ierr = MatCreateShell(this->comm().get(),
1440  rhs_in.local_size(),
1441  solution_in.local_size(),
1442  rhs_in.size(),
1443  solution_in.size(),
1444  const_cast<void*>(static_cast<const void*>(&shell_matrix)),
1445  &mat);
1446  /* Note that the const_cast above is only necessary because PETSc
1447  does not accept a const void*. Inside the member function
1448  _petsc_shell_matrix() below, the pointer is casted back to a
1449  const ShellMatrix<T>*. */
1450 
1451  LIBMESH_CHKERRABORT(ierr);
1452  ierr = MatShellSetOperation(mat,MATOP_MULT,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_mult));
1453  LIBMESH_CHKERRABORT(ierr);
1454  ierr = MatShellSetOperation(mat,MATOP_MULT_ADD,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_mult_add));
1455  LIBMESH_CHKERRABORT(ierr);
1456  ierr = MatShellSetOperation(mat,MATOP_GET_DIAGONAL,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_get_diagonal));
1457  LIBMESH_CHKERRABORT(ierr);
1458 
1459  // Restrict rhs and solution vectors and set operators. The input
1460  // matrix works as the preconditioning matrix.
1461  if(_restrict_solve_to_is!=NULL)
1462  {
1463  size_t is_local_size = this->_restrict_solve_to_is_local_size();
1464 
1465  ierr = VecCreate(this->comm().get(),&subrhs);
1466  LIBMESH_CHKERRABORT(ierr);
1467  ierr = VecSetSizes(subrhs,is_local_size,PETSC_DECIDE);
1468  LIBMESH_CHKERRABORT(ierr);
1469  ierr = VecSetFromOptions(subrhs);
1470  LIBMESH_CHKERRABORT(ierr);
1471 
1472  ierr = VecCreate(this->comm().get(),&subsolution);
1473  LIBMESH_CHKERRABORT(ierr);
1474  ierr = VecSetSizes(subsolution,is_local_size,PETSC_DECIDE);
1475  LIBMESH_CHKERRABORT(ierr);
1476  ierr = VecSetFromOptions(subsolution);
1477  LIBMESH_CHKERRABORT(ierr);
1478 
1479  ierr = VecScatterCreate(rhs->vec(),_restrict_solve_to_is, subrhs,NULL, &scatter);
1480  LIBMESH_CHKERRABORT(ierr);
1481 
1482  ierr = VecScatterBegin(scatter,rhs->vec(),subrhs,INSERT_VALUES,SCATTER_FORWARD);
1483  LIBMESH_CHKERRABORT(ierr);
1484  ierr = VecScatterEnd(scatter,rhs->vec(),subrhs,INSERT_VALUES,SCATTER_FORWARD);
1485  LIBMESH_CHKERRABORT(ierr);
1486 
1487  ierr = VecScatterBegin(scatter,solution->vec(),subsolution,INSERT_VALUES,SCATTER_FORWARD);
1488  LIBMESH_CHKERRABORT(ierr);
1489  ierr = VecScatterEnd(scatter,solution->vec(),subsolution,INSERT_VALUES,SCATTER_FORWARD);
1490  LIBMESH_CHKERRABORT(ierr);
1491 
1492 #if PETSC_VERSION_LESS_THAN(3,1,0)
1493  /* This point can't be reached, see above. */
1494  libmesh_assert(false);
1495 #else
1496  ierr = MatGetSubMatrix(mat,
1498  MAT_INITIAL_MATRIX,&submat);
1499  LIBMESH_CHKERRABORT(ierr);
1500  ierr = MatGetSubMatrix(const_cast<PetscMatrix<T>*>(precond)->mat(),
1502  MAT_INITIAL_MATRIX,&subprecond);
1503  LIBMESH_CHKERRABORT(ierr);
1504 #endif
1505 
1506  /* Since removing columns of the matrix changes the equation
1507  system, we will now change the right hand side to compensate
1508  for this. Note that this is not necessary if \p SUBSET_ZERO
1509  has been selected. */
1511  {
1512  _create_complement_is(rhs_in);
1513  size_t is_complement_local_size = rhs_in.local_size()-is_local_size;
1514 
1515  Vec subvec1 = NULL;
1516  Mat submat1 = NULL;
1517  VecScatter scatter1 = NULL;
1518 
1519  ierr = VecCreate(this->comm().get(),&subvec1);
1520  LIBMESH_CHKERRABORT(ierr);
1521  ierr = VecSetSizes(subvec1,is_complement_local_size,PETSC_DECIDE);
1522  LIBMESH_CHKERRABORT(ierr);
1523  ierr = VecSetFromOptions(subvec1);
1524  LIBMESH_CHKERRABORT(ierr);
1525 
1526  ierr = VecScatterCreate(rhs->vec(),_restrict_solve_to_is_complement, subvec1,NULL, &scatter1);
1527  LIBMESH_CHKERRABORT(ierr);
1528 
1529  ierr = VecScatterBegin(scatter1,_subset_solve_mode==SUBSET_COPY_RHS ? rhs->vec() : solution->vec(),subvec1,INSERT_VALUES,SCATTER_FORWARD);
1530  LIBMESH_CHKERRABORT(ierr);
1531  ierr = VecScatterEnd(scatter1,_subset_solve_mode==SUBSET_COPY_RHS ? rhs->vec() : solution->vec(),subvec1,INSERT_VALUES,SCATTER_FORWARD);
1532  LIBMESH_CHKERRABORT(ierr);
1533 
1534  ierr = VecScale(subvec1,-1.0);
1535  LIBMESH_CHKERRABORT(ierr);
1536 
1537 #if PETSC_VERSION_LESS_THAN(3,1,0)
1538  /* This point can't be reached, see above. */
1539  libmesh_assert(false);
1540 #else
1541  ierr = MatGetSubMatrix(mat,
1543  MAT_INITIAL_MATRIX,&submat1);
1544  LIBMESH_CHKERRABORT(ierr);
1545 #endif
1546 
1547  // The following lines would be correct, but don't work
1548  // correctly in PETSc up to 3.1.0-p5. See discussion in
1549  // petsc-users of Nov 9, 2010.
1550  //
1551  // ierr = MatMultAdd(submat1,subvec1,subrhs,subrhs);
1552  // LIBMESH_CHKERRABORT(ierr);
1553  //
1554  // We workaround by using a temporary vector. Note that the
1555  // fix in PETsc 3.1.0-p6 uses a temporary vector internally,
1556  // so this is no effective performance loss.
1557  Vec subvec2 = NULL;
1558  ierr = VecCreate(this->comm().get(),&subvec2);
1559  LIBMESH_CHKERRABORT(ierr);
1560  ierr = VecSetSizes(subvec2,is_local_size,PETSC_DECIDE);
1561  LIBMESH_CHKERRABORT(ierr);
1562  ierr = VecSetFromOptions(subvec2);
1563  LIBMESH_CHKERRABORT(ierr);
1564  ierr = MatMult(submat1,subvec1,subvec2);
1565  LIBMESH_CHKERRABORT(ierr);
1566  ierr = VecAXPY(subrhs,1.0,subvec2);
1567 
1568  ierr = LibMeshVecScatterDestroy(&scatter1);
1569  LIBMESH_CHKERRABORT(ierr);
1570  ierr = LibMeshVecDestroy(&subvec1);
1571  LIBMESH_CHKERRABORT(ierr);
1572  ierr = LibMeshMatDestroy(&submat1);
1573  LIBMESH_CHKERRABORT(ierr);
1574  }
1575 
1576  ierr = KSPSetOperators(_ksp, submat, subprecond,
1577  DIFFERENT_NONZERO_PATTERN);
1578  LIBMESH_CHKERRABORT(ierr);
1579 
1580  if(this->_preconditioner)
1581  {
1582  subprecond_matrix = new PetscMatrix<Number>(subprecond,
1583  this->comm());
1584  this->_preconditioner->set_matrix(*subprecond_matrix);
1585  this->_preconditioner->init();
1586  }
1587  }
1588  else
1589  {
1590  ierr = KSPSetOperators(_ksp, mat, const_cast<PetscMatrix<T>*>(precond)->mat(),
1591  DIFFERENT_NONZERO_PATTERN);
1592  LIBMESH_CHKERRABORT(ierr);
1593 
1594  if(this->_preconditioner)
1595  {
1596  this->_preconditioner->set_matrix(const_cast<SparseMatrix<Number>&>(precond_matrix));
1597  this->_preconditioner->init();
1598  }
1599  }
1600 
1601  // Set the tolerances for the iterative solver. Use the user-supplied
1602  // tolerance for the relative residual & leave the others at default values.
1603  ierr = KSPSetTolerances (_ksp, tol, PETSC_DEFAULT,
1604  PETSC_DEFAULT, max_its);
1605  LIBMESH_CHKERRABORT(ierr);
1606 
1607  // Solve the linear system
1608  if(_restrict_solve_to_is!=NULL)
1609  {
1610  ierr = KSPSolve (_ksp, subrhs, subsolution);
1611  LIBMESH_CHKERRABORT(ierr);
1612  }
1613  else
1614  {
1615  ierr = KSPSolve (_ksp, rhs->vec(), solution->vec());
1616  LIBMESH_CHKERRABORT(ierr);
1617  }
1618 
1619  // Get the number of iterations required for convergence
1620  ierr = KSPGetIterationNumber (_ksp, &its);
1621  LIBMESH_CHKERRABORT(ierr);
1622 
1623  // Get the norm of the final residual to return to the user.
1624  ierr = KSPGetResidualNorm (_ksp, &final_resid);
1625  LIBMESH_CHKERRABORT(ierr);
1626 
1627  if(_restrict_solve_to_is!=NULL)
1628  {
1629  switch(_subset_solve_mode)
1630  {
1631  case SUBSET_ZERO:
1632  ierr = VecZeroEntries(solution->vec());
1633  LIBMESH_CHKERRABORT(ierr);
1634  break;
1635 
1636  case SUBSET_COPY_RHS:
1637  ierr = VecCopy(rhs->vec(),solution->vec());
1638  LIBMESH_CHKERRABORT(ierr);
1639  break;
1640 
1641  case SUBSET_DONT_TOUCH:
1642  /* Nothing to do here. */
1643  break;
1644 
1645  }
1646  ierr = VecScatterBegin(scatter,subsolution,solution->vec(),INSERT_VALUES,SCATTER_REVERSE);
1647  LIBMESH_CHKERRABORT(ierr);
1648  ierr = VecScatterEnd(scatter,subsolution,solution->vec(),INSERT_VALUES,SCATTER_REVERSE);
1649  LIBMESH_CHKERRABORT(ierr);
1650 
1651  ierr = LibMeshVecScatterDestroy(&scatter);
1652  LIBMESH_CHKERRABORT(ierr);
1653 
1654  if(this->_preconditioner)
1655  {
1656  /* Before we delete subprecond_matrix, we should give the
1657  _preconditioner a different matrix. */
1658  this->_preconditioner->set_matrix(const_cast<SparseMatrix<Number>&>(precond_matrix));
1659  this->_preconditioner->init();
1660  delete subprecond_matrix;
1661  subprecond_matrix = NULL;
1662  }
1663 
1664  ierr = LibMeshVecDestroy(&subsolution);
1665  LIBMESH_CHKERRABORT(ierr);
1666  ierr = LibMeshVecDestroy(&subrhs);
1667  LIBMESH_CHKERRABORT(ierr);
1668  ierr = LibMeshMatDestroy(&submat);
1669  LIBMESH_CHKERRABORT(ierr);
1670  ierr = LibMeshMatDestroy(&subprecond);
1671  LIBMESH_CHKERRABORT(ierr);
1672  }
1673 
1674  // Destroy the matrix.
1675  ierr = LibMeshMatDestroy(&mat);
1676  LIBMESH_CHKERRABORT(ierr);
1677 
1678  STOP_LOG("solve()", "PetscLinearSolver");
1679  // return the # of its. and the final residual norm.
1680  return std::make_pair(its, final_resid);
1681 
1682 #endif
1683 
1684 }
template<typename T>
std::pair< unsigned int, Real > libMesh::LinearSolver< T >::solve ( const ShellMatrix< T > &  matrix,
const SparseMatrix< T > *  precond_matrix,
NumericVector< T > &  sol,
NumericVector< T > &  rhs,
const double  tol,
const unsigned int  n_iter 
)
inlineinherited

This function solves a system whose matrix is a shell matrix, but an optional sparse matrix may be used as preconditioning matrix.

Definition at line 315 of file linear_solver.h.

321 {
322  if (pc_mat)
323  return this->solve(mat, *pc_mat, sol, rhs, tol, n_iter);
324  else
325  return this->solve(mat, sol, rhs, tol, n_iter);
326 }
template<typename T>
SolverType libMesh::LinearSolver< T >::solver_type ( ) const
inlineinherited

Returns the type of solver to use.

Definition at line 100 of file linear_solver.h.

100 { return _solver_type; }

Member Data Documentation

ReferenceCounter::Counts libMesh::ReferenceCounter::_counts
staticprotectedinherited
bool libMesh::ReferenceCounter::_enable_print_counter = true
staticprotectedinherited

Flag to control whether reference count information is printed when print_info is called.

Definition at line 137 of file reference_counter.h.

Referenced by libMesh::ReferenceCounter::disable_print_counter_info(), libMesh::ReferenceCounter::enable_print_counter_info(), and libMesh::ReferenceCounter::print_info().

template<typename T>
bool libMesh::LinearSolver< T >::_is_initialized
protectedinherited

Flag indicating if the data structures have been initialized.

Definition at line 246 of file linear_solver.h.

Referenced by libMesh::LinearSolver< Number >::initialized().

template<typename T >
KSP libMesh::PetscLinearSolver< T >::_ksp
private

Krylov subspace context

Definition at line 285 of file petsc_linear_solver.h.

Referenced by libMesh::PetscLinearSolver< T >::ksp().

Threads::spin_mutex libMesh::ReferenceCounter::_mutex
staticprotectedinherited

Mutual exclusion object to enable thread-safe reference counting.

Definition at line 131 of file reference_counter.h.

Threads::atomic< unsigned int > libMesh::ReferenceCounter::_n_objects
staticprotectedinherited

The number of objects. Print the reference count information when the number returns to 0.

Definition at line 126 of file reference_counter.h.

Referenced by libMesh::ReferenceCounter::n_objects(), libMesh::ReferenceCounter::ReferenceCounter(), and libMesh::ReferenceCounter::~ReferenceCounter().

template<typename T >
PC libMesh::PetscLinearSolver< T >::_pc
private

Preconditioner context

Definition at line 280 of file petsc_linear_solver.h.

Referenced by libMesh::PetscLinearSolver< T >::pc().

template<typename T>
Preconditioner<T>* libMesh::LinearSolver< T >::_preconditioner
protectedinherited

Holds the Preconditioner object to be used for the linear solves.

Definition at line 251 of file linear_solver.h.

template<typename T>
PreconditionerType libMesh::LinearSolver< T >::_preconditioner_type
protectedinherited

Enum statitng with type of preconditioner to use.

Definition at line 241 of file linear_solver.h.

Referenced by libMesh::AztecLinearSolver< T >::AztecLinearSolver(), and libMesh::PetscLinearSolver< T >::PetscLinearSolver().

template<typename T >
IS libMesh::PetscLinearSolver< T >::_restrict_solve_to_is
private

PETSc index set containing the dofs on which to solve (NULL means solve on all dofs).

Definition at line 291 of file petsc_linear_solver.h.

template<typename T >
IS libMesh::PetscLinearSolver< T >::_restrict_solve_to_is_complement
private

PETSc index set, complement to _restrict_solve_to_is. This will be created on demand by the method _create_complement_is().

Definition at line 298 of file petsc_linear_solver.h.

template<typename T >
SLES libMesh::PetscLinearSolver< T >::_sles
private

Linear solver context

Definition at line 273 of file petsc_linear_solver.h.

template<typename T>
SolverType libMesh::LinearSolver< T >::_solver_type
protectedinherited

Enum stating which type of iterative solver to use.

Definition at line 236 of file linear_solver.h.

Referenced by libMesh::LinearSolver< Number >::set_solver_type(), and libMesh::LinearSolver< Number >::solver_type().

template<typename T >
SubsetSolveMode libMesh::PetscLinearSolver< T >::_subset_solve_mode
private

If restrict-solve-to-subset mode is active, this member decides what happens with the dofs outside the subset.

Definition at line 317 of file petsc_linear_solver.h.

template<typename T>
bool libMesh::LinearSolver< T >::same_preconditioner
protectedinherited

Boolean flag to indicate whether we want to use an identical preconditioner to the previous solve. This can save substantial work in the cases where the system matrix is the same for successive solves.

Definition at line 259 of file linear_solver.h.


The documentation for this class was generated from the following files:

Site Created By: libMesh Developers
Last modified: February 07 2014 16:58:01 UTC

Hosted By:
SourceForge.net Logo