#include <petsc_linear_solver.h>
PetscLinearSolver ()
~PetscLinearSolver ()
void clear ()
void init ()
void init (PetscMatrix< T > *matrix)
std::pair< unsigned int, Real > solve (SparseMatrix< T > &matrix_in, NumericVector< T > &solution_in, NumericVector< T > &rhs_in, const double tol, const unsigned int m_its)
std::pair< unsigned int, Real > solve (SparseMatrix< T > &matrix, SparseMatrix< T > &preconditioner, NumericVector< T > &solution, NumericVector< T > &rhs, const double tol, const unsigned int m_its)
std::pair< unsigned int, Real > solve (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)
static AutoPtr< LinearSolver< T > > build (const SolverPackage solver_package=libMesh::default_solver_package())
static std::string get_info ()
static void print_info ()
static unsigned int n_objects ()
typedef std::map< std::string, std::pair< unsigned int, unsigned int > > Counts
void increment_constructor_count (const std::string &name)
void increment_destructor_count (const std::string &name)
SolverType _solver_type
PreconditionerType _preconditioner_type
bool _is_initialized
Preconditioner< T > * _preconditioner
static Counts _counts
static Threads::atomic< unsigned int > _n_objects
static Threads::spin_mutex _mutex
static PetscErrorCode _petsc_shell_matrix_mult (Mat mat, Vec arg, Vec dest)
static PetscErrorCode _petsc_shell_matrix_get_diagonal (Mat mat, Vec dest)
Author:
Definition at line 93 of file petsc_linear_solver.h.
Definition at line 105 of file reference_counter.h.
Definition at line 257 of file petsc_linear_solver.h.
References libMeshEnums::BLOCK_JACOBI_PRECOND, libMeshEnums::ILU_PRECOND, and libMesh::n_processors().
{
if (libMesh::n_processors() == 1)
this->_preconditioner_type = ILU_PRECOND;
else
this->_preconditioner_type = BLOCK_JACOBI_PRECOND;
}
Definition at line 269 of file petsc_linear_solver.h.
{
this->clear ();
}
Definition at line 886 of file petsc_linear_solver.C.
References libMesh::COMM_WORLD, and ShellMatrix< T >::get_diagonal().
{
/* Get the matrix context. */
int ierr=0;
void* ctx;
ierr = MatShellGetContext(mat,&ctx);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
/* Get user shell matrix object. */
const ShellMatrix<T>& shell_matrix = *static_cast<const ShellMatrix<T>*>(ctx);
/* Make
NumericVector instances around the vector. */
PetscVector<T> dest_global(dest);
/* Call the user function. */
shell_matrix.get_diagonal(dest_global);
return ierr;
}
Definition at line 862 of file petsc_linear_solver.C.
References libMesh::COMM_WORLD, and ShellMatrix< T >::vector_mult().
{
/* Get the matrix context. */
int ierr=0;
void* ctx;
ierr = MatShellGetContext(mat,&ctx);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
/* Get user shell matrix object. */
const ShellMatrix<T>& shell_matrix = *static_cast<const ShellMatrix<T>*>(ctx);
/* Make
NumericVector instances around the vectors. */
PetscVector<T> arg_global(arg);
PetscVector<T> dest_global(dest);
/* Call the user function. */
shell_matrix.vector_mult(dest_global,arg_global);
return ierr;
}
Definition at line 103 of file linear_solver.C.
References libMesh::libMeshPrivateData::_is_initialized, and libMeshEnums::SHELL_PRECOND.
{
if(this->_is_initialized)
{
std::cerr<<'Preconditioner must be attached before the solver is initialized!'<<std::endl;
libmesh_error();
}
_preconditioner_type = SHELL_PRECOND;
_preconditioner = preconditioner;
}
Definition at line 37 of file linear_solver.C.
References LASPACK_SOLVERS, libMeshEnums::PETSC_SOLVERS, and TRILINOS_SOLVERS.
Referenced by LegacyXdrIO::read_mesh().
{
// Build the appropriate solver
switch (solver_package)
{
#ifdef LIBMESH_HAVE_LASPACK
case LASPACK_SOLVERS:
{
AutoPtr<LinearSolver<T> > ap(new LaspackLinearSolver<T>);
return ap;
}
#endif
#ifdef LIBMESH_HAVE_PETSC
case PETSC_SOLVERS:
{
AutoPtr<LinearSolver<T> > ap(new PetscLinearSolver<T>);
return ap;
}
#endif
#ifdef LIBMESH_HAVE_TRILINOS
case TRILINOS_SOLVERS:
{
AutoPtr<LinearSolver<T> > ap(new AztecLinearSolver<T>);
return ap;
}
#endif
default:
std::cerr << 'ERROR: Unrecognized solver package: '
<< solver_package
<< std::endl;
libmesh_error();
}
AutoPtr<LinearSolver<T> > ap(NULL);
return ap;
}
Reimplemented from LinearSolver< T >.
Definition at line 91 of file petsc_linear_solver.C.
References libMesh::libMeshPrivateData::_is_initialized, libMeshEnums::BLOCK_JACOBI_PRECOND, libMesh::COMM_WORLD, libMeshEnums::GMRES, libMeshEnums::ILU_PRECOND, libMesh::initialized(), and libMesh::n_processors().
{
if (this->initialized())
{
this->_is_initialized = false;
int ierr=0;
#if PETSC_VERSION_LESS_THAN(2,2,0)
// 2.1.x & earlier style
ierr = SLESDestroy(_sles);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
#else
// 2.2.0 & newer style
ierr = KSPDestroy(_ksp);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
#endif
// Mimic PETSc default solver and preconditioner
this->_solver_type = GMRES;
if(!this->_preconditioner)
{
if (libMesh::n_processors() == 1)
this->_preconditioner_type = ILU_PRECOND;
else
this->_preconditioner_type = BLOCK_JACOBI_PRECOND;
}
}
}
Definition at line 45 of file reference_counter.C.
References ReferenceCounter::_counts, and Quality::name().
Referenced by ReferenceCounter::print_info().
{
#if defined(LIBMESH_ENABLE_REFERENCE_COUNTING) && defined(DEBUG)
std::ostringstream out;
out << '
<< ' ----------------------------------------------------------------------------
<< '| Reference count information |
<< ' ---------------------------------------------------------------------------- ;
for (Counts::iterator it = _counts.begin();
it != _counts.end(); ++it)
{
const std::string name(it->first);
const unsigned int creations = it->second.first;
const unsigned int destructions = it->second.second;
out << '| ' << name << ' reference count information:
<< '| Creations: ' << creations << '
<< '| Destructions: ' << destructions << ';
}
out << ' ---------------------------------------------------------------------------- ;
return out.str();
#else
return '';
#endif
}
Definition at line 718 of file petsc_linear_solver.C.
References libMesh::COMM_WORLD.
{
int ierr = 0;
int its = 0;
// Fill the residual history vector with the residual norms
// Note that GetResidualHistory() does not copy any values, it
// simply sets the pointer p. Note that for some Krylov subspace
// methods, the number of residuals returned in the history
// vector may be different from what you are expecting. For
// example, TFQMR returns two residual values per iteration step.
PetscReal* p;
ierr = KSPGetResidualHistory(_ksp, &p, &its);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
// Check no residual history
if (its == 0)
{
std::cerr << 'No iterations have been performed, returning 0.' << std::endl;
return 0.;
}
// Otherwise, return the value pointed to by p.
return *p;
}
Definition at line 685 of file petsc_linear_solver.C.
References libMesh::COMM_WORLD.
{
int ierr = 0;
int its = 0;
// Fill the residual history vector with the residual norms
// Note that GetResidualHistory() does not copy any values, it
// simply sets the pointer p. Note that for some Krylov subspace
// methods, the number of residuals returned in the history
// vector may be different from what you are expecting. For
// example, TFQMR returns two residual values per iteration step.
PetscReal* p;
ierr = KSPGetResidualHistory(_ksp, &p, &its);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
// Check for early return
if (its == 0) return;
// Create space to store the result
hist.resize(its);
// Copy history into the vector provided by the user.
for (int i=0; i<its; ++i)
{
hist[i] = *p;
p++;
}
}
Definition at line 149 of file reference_counter.h.
References ReferenceCounter::_counts, Quality::name(), and Threads::spin_mtx.
Referenced by ReferenceCountedObject< Value >::ReferenceCountedObject().
{
Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
std::pair<unsigned int, unsigned int>& p = _counts[name];
p.first++;
}
Definition at line 167 of file reference_counter.h.
References ReferenceCounter::_counts, Quality::name(), and Threads::spin_mtx.
Referenced by ReferenceCountedObject< Value >::~ReferenceCountedObject().
{
Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
std::pair<unsigned int, unsigned int>& p = _counts[name];
p.second++;
}
Implements LinearSolver< T >.
Definition at line 130 of file petsc_linear_solver.C.
References __libmesh_petsc_preconditioner_apply(), __libmesh_petsc_preconditioner_setup(), libMesh::libMeshPrivateData::_is_initialized, libMesh::COMM_WORLD, libMesh::initialized(), and PetscPreconditioner< T >::set_petsc_preconditioner_type().
Referenced by PetscLinearSolver< T >::ksp(), and PetscLinearSolver< T >::pc().
{
// Initialize the data structures if not done so already.
if (!this->initialized())
{
this->_is_initialized = true;
int ierr=0;
// 2.1.x & earlier style
#if PETSC_VERSION_LESS_THAN(2,2,0)
// Create the linear solver context
ierr = SLESCreate (libMesh::COMM_WORLD, &_sles);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
// Create the Krylov subspace & preconditioner contexts
ierr = SLESGetKSP (_sles, &_ksp);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
ierr = SLESGetPC (_sles, &_pc);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
// Have the Krylov subspace method use our good initial guess rather than 0
ierr = KSPSetInitialGuessNonzero (_ksp, PETSC_TRUE);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
// Set user-specified solver and preconditioner types
this->set_petsc_solver_type();
// Set the options from user-input
// Set runtime options, e.g.,
// -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
// These options will override those specified above as long as
// SLESSetFromOptions() is called _after_ any other customization
// routines.
ierr = SLESSetFromOptions (_sles);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
// 2.2.0 & newer style
#else
// Create the linear solver context
ierr = KSPCreate (libMesh::COMM_WORLD, &_ksp);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
// Create the preconditioner context
ierr = KSPGetPC (_ksp, &_pc);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
// Have the Krylov subspace method use our good initial guess rather than 0
ierr = KSPSetInitialGuessNonzero (_ksp, PETSC_TRUE);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
// Set user-specified solver and preconditioner types
this->set_petsc_solver_type();
// Set the options from user-input
// Set runtime options, e.g.,
// -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
// These options will override those specified above as long as
// KSPSetFromOptions() is called _after_ any other customization
// routines.
ierr = KSPSetFromOptions (_ksp);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
// Not sure if this is necessary, or if it is already handled by KSPSetFromOptions?
// NOT NECESSARY!!!!
//ierr = PCSetFromOptions (_pc);
//CHKERRABORT(libMesh::COMM_WORLD,ierr);
#endif
// Notify PETSc of location to store residual history.
// This needs to be called before any solves, since
// it sets the residual history length to zero. The default
// behavior is for PETSc to allocate (internally) an array
// of size 1000 to hold the residual norm history.
ierr = KSPSetResidualHistory(_ksp,
PETSC_NULL, // pointer to the array which holds the history
PETSC_DECIDE, // size of the array holding the history
PETSC_TRUE); // Whether or not to reset the history for each solve.
CHKERRABORT(libMesh::COMM_WORLD,ierr);
PetscPreconditioner<T>::set_petsc_preconditioner_type(this->_preconditioner_type,_pc);
//If there is a preconditioner object we need to set the internal setup and apply routines
if(this->_preconditioner)
{
PCShellSetContext(_pc,(void*)this->_preconditioner);
PCShellSetSetUp(_pc,__libmesh_petsc_preconditioner_setup);
PCShellSetApply(_pc,__libmesh_petsc_preconditioner_apply);
}
}
}
Definition at line 230 of file petsc_linear_solver.C.
References __libmesh_petsc_preconditioner_apply(), __libmesh_petsc_preconditioner_setup(), libMesh::libMeshPrivateData::_is_initialized, libMesh::COMM_WORLD, libMesh::initialized(), PetscMatrix< T >::mat(), and PetscPreconditioner< T >::set_petsc_preconditioner_type().
{
// Initialize the data structures if not done so already.
if (!this->initialized())
{
this->_is_initialized = true;
int ierr=0;
// 2.1.x & earlier style
#if PETSC_VERSION_LESS_THAN(2,2,0)
// Create the linear solver context
ierr = SLESCreate (libMesh::COMM_WORLD, &_sles);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
// Create the Krylov subspace & preconditioner contexts
ierr = SLESGetKSP (_sles, &_ksp);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
ierr = SLESGetPC (_sles, &_pc);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
// Have the Krylov subspace method use our good initial guess rather than 0
ierr = KSPSetInitialGuessNonzero (_ksp, PETSC_TRUE);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
// Set user-specified solver and preconditioner types
this->set_petsc_solver_type();
// Set the options from user-input
// Set runtime options, e.g.,
// -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
// These options will override those specified above as long as
// SLESSetFromOptions() is called _after_ any other customization
// routines.
ierr = SLESSetFromOptions (_sles);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
// 2.2.0 & newer style
#else
// Create the linear solver context
ierr = KSPCreate (libMesh::COMM_WORLD, &_ksp);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
//ierr = PCCreate (libMesh::COMM_WORLD, &_pc);
// CHKERRABORT(libMesh::COMM_WORLD,ierr);
// Create the preconditioner context
ierr = KSPGetPC (_ksp, &_pc);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
// Set operators. The input matrix works as the preconditioning matrix
ierr = KSPSetOperators(_ksp, matrix->mat(), matrix->mat(),SAME_NONZERO_PATTERN);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
// Have the Krylov subspace method use our good initial guess rather than 0
ierr = KSPSetInitialGuessNonzero (_ksp, PETSC_TRUE);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
// Set user-specified solver and preconditioner types
this->set_petsc_solver_type();
// Set the options from user-input
// Set runtime options, e.g.,
// -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
// These options will override those specified above as long as
// KSPSetFromOptions() is called _after_ any other customization
// routines.
ierr = KSPSetFromOptions (_ksp);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
// Not sure if this is necessary, or if it is already handled by KSPSetFromOptions?
// NOT NECESSARY!!!!
//ierr = PCSetFromOptions (_pc);
//CHKERRABORT(libMesh::COMM_WORLD,ierr);
#endif
// Notify PETSc of location to store residual history.
// This needs to be called before any solves, since
// it sets the residual history length to zero. The default
// behavior is for PETSc to allocate (internally) an array
// of size 1000 to hold the residual norm history.
ierr = KSPSetResidualHistory(_ksp,
PETSC_NULL, // pointer to the array which holds the history
PETSC_DECIDE, // size of the array holding the history
PETSC_TRUE); // Whether or not to reset the history for each solve.
CHKERRABORT(libMesh::COMM_WORLD,ierr);
PetscPreconditioner<T>::set_petsc_preconditioner_type(this->_preconditioner_type,_pc);
if(this->_preconditioner)
{
this->_preconditioner->set_matrix(*matrix);
PCShellSetContext(_pc,(void*)this->_preconditioner);
PCShellSetSetUp(_pc,__libmesh_petsc_preconditioner_setup);
PCShellSetApply(_pc,__libmesh_petsc_preconditioner_apply);
}
}
}
Definition at line 78 of file linear_solver.h.
{ return _is_initialized; }
Definition at line 192 of file petsc_linear_solver.h.
References PetscLinearSolver< T >::_ksp, and PetscLinearSolver< T >::init().
{ this->init(); return _ksp; }
Definition at line 76 of file reference_counter.h.
References ReferenceCounter::_n_objects.
Referenced by System::read_serialized_blocked_dof_objects(), and System::write_serialized_blocked_dof_objects().
{ return _n_objects; }
Definition at line 185 of file petsc_linear_solver.h.
References PetscLinearSolver< T >::_pc, and PetscLinearSolver< T >::init().
{ this->init(); return _pc; }
Definition at line 83 of file linear_solver.C.
{
if(_preconditioner)
return _preconditioner->type();
return _preconditioner_type;
}
Implements LinearSolver< T >.
Definition at line 799 of file petsc_linear_solver.C.
{
#if PETSC_VERSION_LESS_THAN(2,3,1)
std::cout << 'This method is currently not supported '
<< '(but may work!) for Petsc 2.3.0 and earlier.' << std::endl;
#else
KSPConvergedReason reason;
KSPGetConvergedReason(_ksp, &reason);
// KSP_CONVERGED_RTOL (residual 2-norm decreased by a factor of rtol, from 2-norm of right hand side)
// KSP_CONVERGED_ATOL (residual 2-norm less than abstol)
// KSP_CONVERGED_ITS (used by the preonly preconditioner that always uses ONE iteration)
// KSP_CONVERGED_STEP_LENGTH
// KSP_DIVERGED_ITS (required more than its to reach convergence)
// KSP_DIVERGED_DTOL (residual norm increased by a factor of divtol)
// KSP_DIVERGED_NAN (residual norm became Not-a-number likely do to 0/0)
// KSP_DIVERGED_BREAKDOWN (generic breakdown in method)
switch (reason)
{
case KSP_CONVERGED_RTOL:
{
std::cout << 'Linear solver converged, relative tolerance reached.' << std::endl;
break;
}
case KSP_CONVERGED_ATOL:
{
std::cout << 'Linear solver converged, absolute tolerance reached.' << std::endl;
break;
}
// Divergence
case KSP_DIVERGED_ITS:
{
std::cout << 'Linear solver diverged, max no. of iterations reached.' << std::endl;
break;
}
case KSP_DIVERGED_DTOL:
{
std::cout << 'Linear solver diverged, residual norm increase by dtol (default 1.e5).' << std::endl;
break;
}
case KSP_DIVERGED_NAN:
{
std::cout << 'Linear solver diverged, residual norm is NaN.' << std::endl;
break;
}
case KSP_DIVERGED_BREAKDOWN:
{
std::cout << 'Linear solver diverged, generic breakdown in the method.' << std::endl;
break;
}
default:
{
std::cout << 'Unknown/unsupported con(di)vergence reason: ' << reason << std::endl;
}
}
#endif
}
Definition at line 83 of file reference_counter.C.
References ReferenceCounter::get_info().
{
#if defined(LIBMESH_ENABLE_REFERENCE_COUNTING) && defined(DEBUG)
std::cout << ReferenceCounter::get_info();
#endif
}
Definition at line 748 of file petsc_linear_solver.C.
References libMeshEnums::BICG, libMeshEnums::BICGSTAB, libMeshEnums::CG, libMeshEnums::CGS, libMeshEnums::CHEBYSHEV, libMesh::COMM_WORLD, libMeshEnums::CR, libMeshEnums::GMRES, libMeshEnums::LSQR, libMeshEnums::MINRES, libMeshEnums::RICHARDSON, libMeshEnums::TCQMR, and libMeshEnums::TFQMR.
{
int ierr = 0;
switch (this->_solver_type)
{
case CG:
ierr = KSPSetType (_ksp, (char*) KSPCG); CHKERRABORT(libMesh::COMM_WORLD,ierr); return;
case CR:
ierr = KSPSetType (_ksp, (char*) KSPCR); CHKERRABORT(libMesh::COMM_WORLD,ierr); return;
case CGS:
ierr = KSPSetType (_ksp, (char*) KSPCGS); CHKERRABORT(libMesh::COMM_WORLD,ierr); return;
case BICG:
ierr = KSPSetType (_ksp, (char*) KSPBICG); CHKERRABORT(libMesh::COMM_WORLD,ierr); return;
case TCQMR:
ierr = KSPSetType (_ksp, (char*) KSPTCQMR); CHKERRABORT(libMesh::COMM_WORLD,ierr); return;
case TFQMR:
ierr = KSPSetType (_ksp, (char*) KSPTFQMR); CHKERRABORT(libMesh::COMM_WORLD,ierr); return;
case LSQR:
ierr = KSPSetType (_ksp, (char*) KSPLSQR); CHKERRABORT(libMesh::COMM_WORLD,ierr); return;
case BICGSTAB:
ierr = KSPSetType (_ksp, (char*) KSPBCGS); CHKERRABORT(libMesh::COMM_WORLD,ierr); return;
case MINRES:
ierr = KSPSetType (_ksp, (char*) KSPMINRES); CHKERRABORT(libMesh::COMM_WORLD,ierr); return;
case GMRES:
ierr = KSPSetType (_ksp, (char*) KSPGMRES); CHKERRABORT(libMesh::COMM_WORLD,ierr); return;
case RICHARDSON:
ierr = KSPSetType (_ksp, (char*) KSPRICHARDSON); CHKERRABORT(libMesh::COMM_WORLD,ierr); return;
case CHEBYSHEV:
ierr = KSPSetType (_ksp, (char*) KSPCHEBYCHEV); CHKERRABORT(libMesh::COMM_WORLD,ierr); return;
default:
std::cerr << 'ERROR: Unsupported PETSC Solver: '
<< this->_solver_type << std::endl
<< 'Continuing with PETSC defaults' << std::endl;
}
}
Definition at line 93 of file linear_solver.C.
{
if(_preconditioner)
_preconditioner->set_type(pct);
else
_preconditioner_type = pct;
}
Definition at line 98 of file linear_solver.h.
{ _solver_type = st; }
Implements LinearSolver< T >.
Definition at line 588 of file petsc_linear_solver.C.
References libMesh::COMM_WORLD, libMesh::init(), NumericVector< T >::local_size(), and NumericVector< T >::size().
{
#if PETSC_VERSION_LESS_THAN(2,3,1)
// FIXME[JWP]: There will be a bunch of unused variable warnings
// for older PETScs here.
std::cout << 'This method has been developed with PETSc 2.3.1. '
<< 'No one has made it backwards compatible with older '
<< 'versions of PETSc so far; however, it might work '
<< 'without any change with some older version.' << std::endl;
libmesh_error();
return std::make_pair(0,0.0);
#else
START_LOG('solve()', 'PetscLinearSolver');
// Make sure the data passed in are really of Petsc types
const PetscMatrix<T>* precond = libmesh_cast_ptr<const PetscMatrix<T>*>(&precond_matrix);
PetscVector<T>* solution = libmesh_cast_ptr<PetscVector<T>*>(&solution_in);
PetscVector<T>* rhs = libmesh_cast_ptr<PetscVector<T>*>(&rhs_in);
this->init ();
int ierr=0;
int its=0, max_its = static_cast<int>(m_its);
PetscReal final_resid=0.;
// Close the matrices and vectors in case this wasn't already done.
solution->close ();
rhs->close ();
// Prepare the matrix.
Mat mat;
ierr = MatCreateShell(libMesh::COMM_WORLD,
rhs_in.local_size(),
solution_in.local_size(),
rhs_in.size(),
solution_in.size(),
const_cast<void*>(static_cast<const void*>(&shell_matrix)),
&mat);
/* Note that the const_cast above is only necessary because PETSc
does not accept a const void*. Inside the member function
_petsc_shell_matrix() below, the pointer is casted back to a
const ShellMatrix<T>*. */
CHKERRABORT(libMesh::COMM_WORLD,ierr);
ierr = MatShellSetOperation(mat,MATOP_MULT,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_mult));
ierr = MatShellSetOperation(mat,MATOP_GET_DIAGONAL,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_get_diagonal));
CHKERRABORT(libMesh::COMM_WORLD,ierr);
// Set operators. The input matrix works as the preconditioning matrix
ierr = KSPSetOperators(_ksp, mat, const_cast<PetscMatrix<T>*>(precond)->mat(),
DIFFERENT_NONZERO_PATTERN);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
if(this->_preconditioner)
this->_preconditioner->set_matrix(const_cast<SparseMatrix<Number>&>(precond_matrix));
// Set the tolerances for the iterative solver. Use the user-supplied
// tolerance for the relative residual & leave the others at default values.
ierr = KSPSetTolerances (_ksp, tol, PETSC_DEFAULT,
PETSC_DEFAULT, max_its);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
// Solve the linear system
ierr = KSPSolve (_ksp, rhs->vec(), solution->vec());
CHKERRABORT(libMesh::COMM_WORLD,ierr);
// Get the number of iterations required for convergence
ierr = KSPGetIterationNumber (_ksp, &its);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
// Get the norm of the final residual to return to the user.
ierr = KSPGetResidualNorm (_ksp, &final_resid);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
// Destroy the matrix.
ierr = MatDestroy(mat);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
STOP_LOG('solve()', 'PetscLinearSolver');
// return the # of its. and the final residual norm.
return std::make_pair(its, final_resid);
#endif
}
Implements LinearSolver< T >.
Definition at line 125 of file petsc_linear_solver.h.
{
return this->solve(matrix_in, matrix_in, solution_in, rhs_in, tol, m_its);
}
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 LinearSolver< T >.
Definition at line 342 of file petsc_linear_solver.C.
References PetscMatrix< T >::close(), libMesh::COMM_WORLD, libMesh::init(), and PetscMatrix< T >::mat().
{
START_LOG('solve()', 'PetscLinearSolver');
// Make sure the data passed in are really of Petsc types
PetscMatrix<T>* matrix = libmesh_cast_ptr<PetscMatrix<T>*>(&matrix_in);
PetscMatrix<T>* precond = libmesh_cast_ptr<PetscMatrix<T>*>(&precond_in);
PetscVector<T>* solution = libmesh_cast_ptr<PetscVector<T>*>(&solution_in);
PetscVector<T>* rhs = libmesh_cast_ptr<PetscVector<T>*>(&rhs_in);
this->init (matrix);
int ierr=0;
int its=0, max_its = static_cast<int>(m_its);
PetscReal final_resid=0.;
// Close the matrices and vectors in case this wasn't already done.
matrix->close ();
precond->close ();
solution->close ();
rhs->close ();
// // If matrix != precond, then this means we have specified a
// // special preconditioner, so reset preconditioner type to PCMAT.
// if (matrix != precond)
// {
// this->_preconditioner_type = USER_PRECOND;
// this->set_petsc_preconditioner_type ();
// }
if(this->_preconditioner)
this->_preconditioner->set_matrix(matrix_in);
// 2.1.x & earlier style
#if PETSC_VERSION_LESS_THAN(2,2,0)
// Set operators. The input matrix works as the preconditioning matrix
ierr = SLESSetOperators(_sles, matrix->mat(), precond->mat(),
SAME_NONZERO_PATTERN);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
// Set the tolerances for the iterative solver. Use the user-supplied
// tolerance for the relative residual & leave the others at default values.
ierr = KSPSetTolerances (_ksp, tol, PETSC_DEFAULT,
PETSC_DEFAULT, max_its);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
// Solve the linear system
ierr = SLESSolve (_sles, rhs->vec(), solution->vec(), &its);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
// Get the norm of the final residual to return to the user.
ierr = KSPGetResidualNorm (_ksp, &final_resid);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
// 2.2.0
#elif PETSC_VERSION_LESS_THAN(2,2,1)
// Set operators. The input matrix works as the preconditioning matrix
//ierr = KSPSetOperators(_ksp, matrix->mat(), precond->mat(),
// SAME_NONZERO_PATTERN);
// CHKERRABORT(libMesh::COMM_WORLD,ierr);
// Set the tolerances for the iterative solver. Use the user-supplied
// tolerance for the relative residual & leave the others at default values.
// Convergence is detected at iteration k if
// ||r_k||_2 < max(rtol*||b||_2 , abstol)
// where r_k is the residual vector and b is the right-hand side. Note that
// it is the *maximum* of the two values, the larger of which will almost
// always be rtol*||b||_2.
ierr = KSPSetTolerances (_ksp,
tol, // rtol = relative decrease in residual (1.e-5)
PETSC_DEFAULT, // abstol = absolute convergence tolerance (1.e-50)
PETSC_DEFAULT, // dtol = divergence tolerance (1.e+5)
max_its);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
// Set the solution vector to use
ierr = KSPSetSolution (_ksp, solution->vec());
CHKERRABORT(libMesh::COMM_WORLD,ierr);
// Set the RHS vector to use
ierr = KSPSetRhs (_ksp, rhs->vec());
CHKERRABORT(libMesh::COMM_WORLD,ierr);
// Solve the linear system
ierr = KSPSolve (_ksp);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
// Get the number of iterations required for convergence
ierr = KSPGetIterationNumber (_ksp, &its);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
// Get the norm of the final residual to return to the user.
ierr = KSPGetResidualNorm (_ksp, &final_resid);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
// 2.2.1 & newer style
#else
// Set operators. The input matrix works as the preconditioning matrix
if(!this->same_preconditioner)
{
ierr = KSPSetOperators(_ksp, matrix->mat(), precond->mat(),
SAME_NONZERO_PATTERN);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
}
else
{
ierr = KSPSetOperators(_ksp, matrix->mat(), precond->mat(),
SAME_PRECONDITIONER);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
}
// Set the tolerances for the iterative solver. Use the user-supplied
// tolerance for the relative residual & leave the others at default values.
ierr = KSPSetTolerances (_ksp, tol, PETSC_DEFAULT,
PETSC_DEFAULT, max_its);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
// Solve the linear system
ierr = KSPSolve (_ksp, rhs->vec(), solution->vec());
CHKERRABORT(libMesh::COMM_WORLD,ierr);
// Get the number of iterations required for convergence
ierr = KSPGetIterationNumber (_ksp, &its);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
// Get the norm of the final residual to return to the user.
ierr = KSPGetResidualNorm (_ksp, &final_resid);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
#endif
STOP_LOG('solve()', 'PetscLinearSolver');
// return the # of its. and the final residual norm.
return std::make_pair(its, final_resid);
}
Implements LinearSolver< T >.
Definition at line 495 of file petsc_linear_solver.C.
References PetscVector< T >::close(), libMesh::COMM_WORLD, libMesh::init(), NumericVector< T >::local_size(), NumericVector< T >::size(), and PetscVector< T >::vec().
{
#if PETSC_VERSION_LESS_THAN(2,3,1)
// FIXME[JWP]: There will be a bunch of unused variable warnings
// for older PETScs here.
std::cout << 'This method has been developed with PETSc 2.3.1. '
<< 'No one has made it backwards compatible with older '
<< 'versions of PETSc so far; however, it might work '
<< 'without any change with some older version.' << std::endl;
libmesh_error();
return std::make_pair(0,0.0);
#else
START_LOG('solve()', 'PetscLinearSolver');
// Make sure the data passed in are really of Petsc types
PetscVector<T>* solution = libmesh_cast_ptr<PetscVector<T>*>(&solution_in);
PetscVector<T>* rhs = libmesh_cast_ptr<PetscVector<T>*>(&rhs_in);
this->init ();
int ierr=0;
int its=0, max_its = static_cast<int>(m_its);
PetscReal final_resid=0.;
// Close the matrices and vectors in case this wasn't already done.
solution->close ();
rhs->close ();
// Prepare the matrix.
Mat mat;
ierr = MatCreateShell(libMesh::COMM_WORLD,
rhs_in.local_size(),
solution_in.local_size(),
rhs_in.size(),
solution_in.size(),
const_cast<void*>(static_cast<const void*>(&shell_matrix)),
&mat);
/* Note that the const_cast above is only necessary because PETSc
does not accept a const void*. Inside the member function
_petsc_shell_matrix() below, the pointer is casted back to a
const ShellMatrix<T>*. */
CHKERRABORT(libMesh::COMM_WORLD,ierr);
ierr = MatShellSetOperation(mat,MATOP_MULT,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_mult));
ierr = MatShellSetOperation(mat,MATOP_GET_DIAGONAL,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_get_diagonal));
CHKERRABORT(libMesh::COMM_WORLD,ierr);
// Set operators. The input matrix works as the preconditioning matrix
ierr = KSPSetOperators(_ksp, mat, mat,
SAME_NONZERO_PATTERN);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
// Set the tolerances for the iterative solver. Use the user-supplied
// tolerance for the relative residual & leave the others at default values.
ierr = KSPSetTolerances (_ksp, tol, PETSC_DEFAULT,
PETSC_DEFAULT, max_its);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
// Solve the linear system
ierr = KSPSolve (_ksp, rhs->vec(), solution->vec());
CHKERRABORT(libMesh::COMM_WORLD,ierr);
// Get the number of iterations required for convergence
ierr = KSPGetIterationNumber (_ksp, &its);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
// Get the norm of the final residual to return to the user.
ierr = KSPGetResidualNorm (_ksp, &final_resid);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
// Destroy the matrix.
ierr = MatDestroy(mat);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
STOP_LOG('solve()', 'PetscLinearSolver');
// return the # of its. and the final residual norm.
return std::make_pair(its, final_resid);
#endif
}
Definition at line 93 of file linear_solver.h.
{ return _solver_type; }
Definition at line 110 of file reference_counter.h.
Referenced by ReferenceCounter::get_info(), ReferenceCounter::increment_constructor_count(), and ReferenceCounter::increment_destructor_count().
Definition at line 200 of file linear_solver.h.
Referenced by LinearSolver< Number >::initialized().
Definition at line 250 of file petsc_linear_solver.h.
Referenced by PetscLinearSolver< T >::ksp().
Definition at line 123 of file reference_counter.h.
Definition at line 118 of file reference_counter.h.
Referenced by ReferenceCounter::n_objects(), ReferenceCounter::ReferenceCounter(), and ReferenceCounter::~ReferenceCounter().
Definition at line 245 of file petsc_linear_solver.h.
Referenced by PetscLinearSolver< T >::pc().
Definition at line 205 of file linear_solver.h.
Definition at line 195 of file linear_solver.h.
Definition at line 238 of file petsc_linear_solver.h.
Definition at line 190 of file linear_solver.h.
Referenced by LinearSolver< Number >::set_solver_type(), and LinearSolver< Number >::solver_type().
Definition at line 182 of file linear_solver.h.
Generated automatically by Doxygen for libMesh from the source code.