/usr/include/libmesh/petsc_linear_solver.h is in libmesh-dev 0.7.1-2ubuntu1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 | // $Id: petsc_linear_solver.h 4278 2011-03-21 15:23:30Z roystgnr $
// The libMesh Finite Element Library.
// Copyright (C) 2002-2008 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
// Lesser General Public License for more details.
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
#ifndef __petsc_linear_solver_h__
#define __petsc_linear_solver_h__
#include "libmesh_config.h"
#ifdef LIBMESH_HAVE_PETSC
// C++ includes
#include <vector>
// Local includes
#include "linear_solver.h"
#include "petsc_macro.h"
/**
* Petsc include files.
*/
EXTERN_C_FOR_PETSC_BEGIN
#if PETSC_VERSION_LESS_THAN(2,2,0)
# include <petscsles.h>
#else
# include <petscksp.h>
#endif
EXTERN_C_FOR_PETSC_END
//--------------------------------------------------------------------
// Functions with C linkage to pass to PETSc. PETSc will call these
// methods as needed for preconditioning
//
// Since they must have C linkage they have no knowledge of a namespace.
// Give them an obscure name to avoid namespace pollution.
extern "C"
{
// Older versions of PETSc do not have the different int typedefs.
// On 64-bit machines, PetscInt may actually be a long long int.
// This change occurred in Petsc-2.2.1.
#if PETSC_VERSION_LESS_THAN(2,2,1)
typedef int PetscErrorCode;
typedef int PetscInt;
#endif
#if PETSC_VERSION_LESS_THAN(3,0,1) && PETSC_VERSION_RELEASE
/**
* This function is called by PETSc to initialize the preconditioner.
* ctx will hold the Preconditioner.
*/
PetscErrorCode __libmesh_petsc_preconditioner_setup (void * ctx);
/**
* This function is called by PETSc to acctually apply the preconditioner.
* ctx will hold the Preconditioner.
*/
PetscErrorCode __libmesh_petsc_preconditioner_apply(void *ctx, Vec x, Vec y);
#else
PetscErrorCode __libmesh_petsc_preconditioner_setup (PC);
PetscErrorCode __libmesh_petsc_preconditioner_apply(PC, Vec x, Vec y);
#endif
} // end extern "C"
namespace libMesh
{
// forward declarations
template <typename T> class PetscMatrix;
/**
* This class provides an interface to PETSc
* iterative solvers that is compatible with the \p libMesh
* \p LinearSolver<>
*
* @author Benjamin Kirk, 2002-2007
*/
template <typename T>
class PetscLinearSolver : public LinearSolver<T>
{
public:
/**
* Constructor. Initializes Petsc data structures
*/
PetscLinearSolver ();
/**
* Destructor.
*/
~PetscLinearSolver ();
/**
* Release all memory and clear data structures.
*/
void clear ();
/**
* Initialize data structures if not done so already.
*/
void init ();
/**
* Initialize data structures if not done so already plus much more
*/
void init (PetscMatrix<T>* matrix);
/**
* 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 \p dofs being a
* \p NULL pointer.
*/
virtual void restrict_solve_to (const std::vector<unsigned int>* const dofs,
const SubsetSolveMode subset_solve_mode=SUBSET_ZERO);
/**
* Call the Petsc solver. It calls the method below, using the
* same matrix for the system and preconditioner matrices.
*/
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)
{
return this->solve(matrix_in, matrix_in, solution_in, rhs_in, tol, m_its);
}
/**
* 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.
*/
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);
/**
* This function solves a system whose matrix is a shell matrix.
*/
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);
/**
* 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.
*/
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);
/**
* 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()!
*/
PC pc() { this->init(); return _pc; }
/**
* Returns the raw PETSc ksp context pointer. This is useful if
* you are for example setting a custom convergence test with
* KSPSetConvergenceTest().
*/
KSP ksp() { this->init(); return _ksp; }
/**
* Fills the input vector with the sequence of residual norms
* from the latest iterative solve.
*/
void get_residual_history(std::vector<double>& hist);
/**
* 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.
*/
Real get_initial_residual();
/**
* Prints a useful message about why the latest linear solve
* con(di)verged.
*/
virtual void print_converged_reason();
private:
/**
* Tells PETSC to use the user-specified solver stored in
* \p _solver_type
*/
void set_petsc_solver_type ();
/**
* Internal function if shell matrix mode is used.
*/
static PetscErrorCode _petsc_shell_matrix_mult(Mat mat, Vec arg, Vec dest);
/**
* Internal function if shell matrix mode is used.
*/
static PetscErrorCode _petsc_shell_matrix_mult_add(Mat mat, Vec arg, Vec add, Vec dest);
/**
* Internal function if shell matrix mode is used.
*/
static PetscErrorCode _petsc_shell_matrix_get_diagonal(Mat mat, Vec dest);
// SLES removed from >= PETSc 2.2.0
#if PETSC_VERSION_LESS_THAN(2,2,0)
/**
* Linear solver context
*/
SLES _sles;
#endif
/**
* Preconditioner context
*/
PC _pc;
/**
* Krylov subspace context
*/
KSP _ksp;
/**
* PETSc index set containing the dofs on which to solve (\p NULL
* means solve on all dofs).
*/
IS _restrict_solve_to_is;
/**
* PETSc index set, complement to \p _restrict_solve_to_is. This
* will be created on demand by the method \p
* _create_complement_is().
*/
IS _restrict_solve_to_is_complement;
/**
* Internal method that returns the local size of \p
* _restrict_solve_to_is.
*/
size_t _restrict_solve_to_is_local_size(void)const;
/**
* Creates \p _restrict_solve_to_is_complement to contain all
* indices that are local in \p vec_in, except those that are
* contained in \p _restrict_solve_to_is.
*/
void _create_complement_is (const NumericVector<T> &vec_in);
/**
* If restrict-solve-to-subset mode is active, this member decides
* what happens with the dofs outside the subset.
*/
SubsetSolveMode _subset_solve_mode;
};
/*----------------------- functions ----------------------------------*/
template <typename T>
inline
PetscLinearSolver<T>::PetscLinearSolver ():
_restrict_solve_to_is(NULL),
_restrict_solve_to_is_complement(NULL),
_subset_solve_mode(SUBSET_ZERO)
{
if (libMesh::n_processors() == 1)
this->_preconditioner_type = ILU_PRECOND;
else
this->_preconditioner_type = BLOCK_JACOBI_PRECOND;
}
template <typename T>
inline
PetscLinearSolver<T>::~PetscLinearSolver ()
{
this->clear ();
}
template <typename T>
inline size_t
PetscLinearSolver<T>::
_restrict_solve_to_is_local_size(void)const
{
libmesh_assert(_restrict_solve_to_is!=NULL);
PetscInt s;
int ierr = ISGetLocalSize(_restrict_solve_to_is,&s);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
return static_cast<size_t>(s);
}
template <typename T>
void
PetscLinearSolver<T>::_create_complement_is (const NumericVector<T> &
#if PETSC_VERSION_LESS_THAN(3,0,0)
// unnamed to avoid compiler "unused parameter" warning
#else
vec_in
#endif
)
{
libmesh_assert(_restrict_solve_to_is!=NULL);
#if PETSC_VERSION_LESS_THAN(3,0,0)
// No ISComplement in PETSc 2.3.3
libmesh_not_implemented();
#else
if(_restrict_solve_to_is_complement==NULL)
{
int ierr = ISComplement(_restrict_solve_to_is,
vec_in.first_local_index(),
vec_in.last_local_index(),
&_restrict_solve_to_is_complement);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
}
#endif
}
} // namespace libMesh
#endif // #ifdef LIBMESH_HAVE_PETSC
#endif // #ifdef __petsc_linear_solver_h__
|