// $Id: petsc_matrix_free.cc 26043 2012-09-24 19:25:57Z steigemann $
// Version: $Name$
//
// Copyright (C) 2012 by the deal.II authors
//
// This file is subject to QPL and may not be distributed
// without copyright and license information. Please refer
// to the file deal.II/doc/license.html for the text and
// further information on this license.
//
//---------------------------------------------------------------------------
#include <deal.II/lac/petsc_matrix_free.h>
#ifdef DEAL_II_USE_PETSC
DEAL_II_NAMESPACE_OPEN
namespace PETScWrappers
{
MatrixFree::MatrixFree ()
: communicator (PETSC_COMM_SELF)
{
const int m=0;
do_reinit (m, m, m, m);
}
MatrixFree::MatrixFree (const MPI_Comm &communicator,
const unsigned int m,
const unsigned int n,
const unsigned int local_rows,
const unsigned int local_columns)
: communicator (communicator)
{
do_reinit (m, n, local_rows, local_columns);
}
MatrixFree::MatrixFree (const MPI_Comm &communicator,
const unsigned int m,
const unsigned int n,
const std::vector<unsigned int> &local_rows_per_process,
const std::vector<unsigned int> &local_columns_per_process,
const unsigned int this_process)
: communicator (communicator)
{
Assert (local_rows_per_process.size() == local_columns_per_process.size(),
ExcDimensionMismatch (local_rows_per_process.size(),
local_columns_per_process.size()));
Assert (this_process < local_rows_per_process.size(),
ExcInternalError());
do_reinit (m, n,
local_rows_per_process[this_process],
local_columns_per_process[this_process]);
}
MatrixFree::MatrixFree (const unsigned int m,
const unsigned int n,
const unsigned int local_rows,
const unsigned int local_columns)
: communicator (MPI_COMM_WORLD)
{
do_reinit (m, n, local_rows, local_columns);
}
MatrixFree::MatrixFree (const unsigned int m,
const unsigned int n,
const std::vector<unsigned int> &local_rows_per_process,
const std::vector<unsigned int> &local_columns_per_process,
const unsigned int this_process)
: communicator (MPI_COMM_WORLD)
{
Assert (local_rows_per_process.size() == local_columns_per_process.size(),
ExcDimensionMismatch (local_rows_per_process.size(),
local_columns_per_process.size()));
Assert (this_process < local_rows_per_process.size(),
ExcInternalError());
do_reinit (m, n,
local_rows_per_process[this_process],
local_columns_per_process[this_process]);
}
void MatrixFree::reinit (const MPI_Comm &communicator,
const unsigned int m,
const unsigned int n,
const unsigned int local_rows,
const unsigned int local_columns)
{
this->communicator = communicator;
// destroy the matrix and
// generate a new one
#if DEAL_II_PETSC_VERSION_LT(3,2,0)
int ierr = MatDestroy (matrix);
#else
int ierr = MatDestroy (&matrix);
#endif
AssertThrow (ierr == 0, ExcPETScError(ierr));
do_reinit (m, n, local_rows, local_columns);
}
void MatrixFree::reinit (const MPI_Comm &communicator,
const unsigned int m,
const unsigned int n,
const std::vector<unsigned int> &local_rows_per_process,
const std::vector<unsigned int> &local_columns_per_process,
const unsigned int this_process)
{
Assert (local_rows_per_process.size() == local_columns_per_process.size(),
ExcDimensionMismatch (local_rows_per_process.size(),
local_columns_per_process.size()));
Assert (this_process < local_rows_per_process.size(),
ExcInternalError());
this->communicator = communicator;
#if DEAL_II_PETSC_VERSION_LT(3,2,0)
int ierr = MatDestroy (matrix);
#else
int ierr = MatDestroy (&matrix);
#endif
AssertThrow (ierr == 0, ExcPETScError(ierr));
do_reinit (m, n,
local_rows_per_process[this_process],
local_columns_per_process[this_process]);
}
void MatrixFree::reinit (const unsigned int m,
const unsigned int n,
const unsigned int local_rows,
const unsigned int local_columns)
{
reinit (MPI_COMM_WORLD, m, n, local_rows, local_columns);
}
void MatrixFree::reinit (const unsigned int m,
const unsigned int n,
const std::vector<unsigned int> &local_rows_per_process,
const std::vector<unsigned int> &local_columns_per_process,
const unsigned int this_process)
{
reinit (MPI_COMM_WORLD, m, n, local_rows_per_process, local_columns_per_process, this_process);
}
void MatrixFree::clear ()
{
#if DEAL_II_PETSC_VERSION_LT(3,2,0)
int ierr = MatDestroy (matrix);
#else
int ierr = MatDestroy (&matrix);
#endif
AssertThrow (ierr == 0, ExcPETScError(ierr));
const int m=0;
do_reinit (m, m, m, m);
}
void MatrixFree::vmult (Vec &dst, const Vec &src) const
{
//TODO: Translate the given PETSc Vec* vector into a deal.II
// vector so we can call the vmult function with the usual
// interface; then convert back. This could be much more
// efficient, if the PETScWrappers::*::Vector classes
// had a way to simply generate such a vector object from
// a given PETSc Vec* object without allocating new memory
// and without taking ownership of the Vec*
VectorBase *x = 0;
VectorBase *y = 0;
// because we do not know,
// if dst and src are sequential
// or distributed vectors,
// we ask for the vector-type
// and reinit x and y with
// dealii::PETScWrappers::*::Vector:
const char *vec_type;
int ierr = VecGetType (src, &vec_type);
PetscInt local_size;
ierr = VecGetLocalSize (src, &local_size);
AssertThrow (ierr == 0, ExcPETScError(ierr));
if (strcmp(vec_type,"mpi") == 0)
{
PetscInt size;
ierr = VecGetSize (src, &size);
AssertThrow (ierr == 0, ExcPETScError(ierr));
x = new PETScWrappers::MPI::Vector (this->get_mpi_communicator (), size, local_size);
y = new PETScWrappers::MPI::Vector (this->get_mpi_communicator (), size, local_size);
}
else if (strcmp(vec_type,"seq") == 0)
{
x = new PETScWrappers::Vector (local_size);
y = new PETScWrappers::Vector (local_size);
}
else
AssertThrow (false, ExcMessage("PETScWrappers::MPI::MatrixFree::do_matrix_vector_action: "
"This only works for Petsc Vec Type = VECMPI | VECSEQ"));
// copy src to x
x->equ(1., PETScWrappers::VectorBase(src));
// and call vmult(x,y) which must
// be reimplemented in derived classes
vmult (*y, *x);
y->compress();
// copy the result back to dst
ierr = VecCopy (&(*(*y)), dst);
AssertThrow (ierr == 0, ExcPETScError(ierr));
delete (x);
delete (y);
}
int MatrixFree::matrix_free_mult (Mat A, Vec src, Vec dst)
{
// create a pointer to this MatrixFree
// object and link the given matrix A
// to the matrix-vector multiplication
// of this MatrixFree object,
MatrixFree *this_object;
int ierr = MatShellGetContext (A, &this_object);
AssertThrow (ierr == 0, ExcPETScError(ierr));
// call vmult of this object:
this_object->vmult (dst, src);
return (0);
}
void MatrixFree::do_reinit (const unsigned int m,
const unsigned int n,
const unsigned int local_rows,
const unsigned int local_columns)
{
Assert (local_rows <= m, ExcDimensionMismatch (local_rows, m));
Assert (local_columns <= n, ExcDimensionMismatch (local_columns, n));
int ierr;
// create a PETSc MatShell matrix-type
// object of dimension m x n and local size
// local_rows x local_columns
ierr = MatCreateShell(communicator, local_rows, local_columns, m, n, (void*)this, &matrix);
AssertThrow (ierr == 0, ExcPETScError(ierr));
// register the MatrixFree::matrix_free_mult function
// as the matrix multiplication used by this matrix
ierr = MatShellSetOperation (matrix, MATOP_MULT,
(void(*)(void))&dealii::PETScWrappers::MatrixFree::matrix_free_mult);
AssertThrow (ierr == 0, ExcPETScError(ierr));
ierr = MatSetFromOptions (matrix);
AssertThrow (ierr == 0, ExcPETScError(ierr));
}
}
DEAL_II_NAMESPACE_CLOSE
#endif // DEAL_II_USE_PETSC
git-svn-id: https://svn.dealii.org/trunk@26688
0785d39b-7218-0410-832d-
ea1e28bc413d
/**
* Milne-rule. Closed Newton-Cotes formula, exact for polynomials of degree 5.
- * See Stoer: Einführung in die Numerische Mathematik I, p. 102
+ * See Stoer: Einf�hrung in die Numerische Mathematik I, p. 102
*/
template <int dim>
class QMilne : public Quadrature<dim>
/**
* Weddle-rule. Closed Newton-Cotes formula, exact for polynomials of degree 7.
- * See Stoer: Einführung in die Numerische Mathematik I, p. 102
+ * See Stoer: Einf�hrung in die Numerische Mathematik I, p. 102
*/
template <int dim>
class QWeddle : public Quadrature<dim>
+/**
+ * Gauss Quadrature Formula with $1/R^{3/2}$ weighting function. This formula
+ * can be used to to integrate $1/R^{3/2} \ f(x)$ on the reference
+ * element $[0,1]^2$, where $f$ is a smooth function without
+ * singularities, and $R$ is the distance from the point $x$ to the vertex
+ * $\xi$, given at construction time by specifying its index. Notice that
+ * this distance is evaluated in the reference element.
+ *
+ * This quadrature formula is a specialization of QGaussOneOverR.
+ * We apply a second transformation $R = t^2$ to cancel a singularity
+ * of order $R^{1/2}$:
+ * \f[
+ * \int_0^1 \int_0^1 \frac{1}{R^{1/2}} f(x,y) \biggr|_{x < y} dxdy =
+ * \int_0^1 \int_0^{r(\pi/4 v)} \frac{1}{R^{1/2}} f(R \cos(\pi/4 v), R \sin(\pi/4 v)) \frac{\pi}{4} R dr dv =
+ * 2 \int_0^1 \int_0^{\sqrt{r(\pi/4 v)}} f(t^2 \cos(\pi/4 v), t^2 \sin(\pi/4 v)) \frac{\pi}{4} t^2 dt dv =
+ * 2 \int_0^1 \int_0^1 f(t^2, t^2 \tan(\pi/4 v)) \frac{\pi}{4} \frac{t^2}{(\cos(\pi/4 v))^{3/2}} dt dv
+ * \f]
+ *
+ * Upon construction it is possible to specify wether we want the
+ * singularity removed, or not. In other words, this quadrature can be
+ * used to integrate $g(x) = 1/R^{3/2}\ f(x)$, or simply $f(x)$, with the $1/R^{3/2}$
+ * factor already included in the quadrature weights.
+ */
+template<int dim>
+class QGaussOneOverRThreeHalfs : public Quadrature<dim>
+{
+ public:
+ /**
+ * The constructor takes three arguments: the order of the Gauss
+ * formula, the index of the vertex where the singularity is
+ * located, and whether we include the weighting singular function
+ * inside the quadrature, or we leave it in the user function to
+ * be integrated. Notice that this constructor only works for the
+ * vertices of the quadrilateral.
+ *
+ * Traditionally, quadrature formulas include their weighting
+ * function, and the last argument is set to false by
+ * default. There are cases, however, where this is undesirable
+ * (for example when you only know that your singularity has the
+ * same order of 1/R, but cannot be written exactly in this
+ * way).
+ *
+ * In other words, you can use this function in either of
+ * the following way, obtaining the same result:
+ *
+ * @code
+ * QGaussOneOverRThreeHalfs singular_quad(order, vertex_id, false);
+ * // This will produce the integral of f(x)/R^{3/2}
+ * for(unsigned int i=0; i<singular_quad.size(); ++i)
+ * integral += f(singular_quad.point(i))*singular_quad.weight(i);
+ *
+ * // And the same here
+ * QGaussOneOverRThreeHalfs singular_quad_noR(order, vertex_id, true);
+ *
+ * // This also will produce the integral of f(x)/R^{3/2}, but 1/R^{3/2} has to
+ * // be specified.
+ * for(unsigned int i=0; i<singular_quad.size(); ++i) {
+ * double R = (singular_quad_noR.point(i)-cell->vertex(vertex_id)).norm();
+ * integral += f(singular_quad_noR.point(i))*singular_quad_noR.weight(i)/R^{3/2};
+ * }
+ * @endcode
+ */
+ QGaussOneOverRThreeHalfs(const unsigned int n,
+ const unsigned int vertex_index,
+ const bool factor_out_singular_weight=false);
+};
+
+
+
/*@}*/
/* -------------- declaration of explicit specializations ------------- */
template <> QGaussLog<1>::QGaussLog (const unsigned int n, const bool revert);
template <> QGaussLogR<1>::QGaussLogR (const unsigned int n, const Point<1> x0, const double alpha, const bool flag);
template <> QGaussOneOverR<2>::QGaussOneOverR (const unsigned int n, const unsigned int index, const bool flag);
+template <> QGaussOneOverRThreeHalfs<2>::QGaussOneOverRThreeHalfs (const unsigned int n, const unsigned int index, const bool flag);
*/
void write_tex (std::ostream &file, const bool with_header=true) const;
+ /**
+ * Write table as a tex file
+ * and use the booktabs package.
+ * Values are written in
+ * math mode: $value$
+ * If with_header is set to false
+ * (it is true by default), then
+ * no "\documentclass{...}",
+ * "\begin{document}" and
+ * "\end{document}" are used. In
+ * this way the file can be
+ * included into an existing tex
+ * file using a command like
+ * "\input{table_file}".
+ */
+ void write_results (std::ostream &file, const bool with_header=true) const;
+
/**
* Read or write the data of this
* object to or from a stream for
}
}
+
+
template <typename T>
void TableHandler::add_value (const std::string &key,
const T value)
bool output_details;
};
-
-
/**
* Empty Constructor. You need to call
* initialize() before using this
*/
AdditionalData additional_data;
};
+
+
+
+/**
+ * A class that implements the interface to use the ParaSails sparse
+ * approximate inverse preconditioner from the HYPRE suite. Note that
+ * PETSc has to be configured with HYPRE (e.g. with --download-hypre=1).
+ *
+ * @ingroup PETScWrappers
+ * @author Martin Steigemann, 2011
+ */
+ class PreconditionParaSails : public PreconditionerBase
+ {
+ public:
+ /**
+ * Standardized data struct to
+ * pipe additional flags to the
+ * preconditioner.
+ */
+ struct AdditionalData
+ {
+ /**
+ * Constructor.
+ */
+ AdditionalData (
+ const unsigned int symmetric = 0,
+ const unsigned int n_levels = 1,
+ const double threshold = 0.1,
+ const double filter = 0.05,
+ const double load_bal = 0.,
+ const bool output_details = false
+ );
+
+ /**
+ * This parameter has the following meanings,
+ * to indicate the symmetry and definiteness
+ * of the problem, and to specify the type
+ * of the preconditioner to construct:
+ * <ul>
+ * <li> @p 0: nonsymmetric and/or indefinite problem, and nonsymmetric preconditioner
+ * <li> @p 1: SPD problem, and SPD (factored) preconditioner
+ * <li> @p 2: nonsymmetric, definite problem, and SPD (factored) preconditioner
+ * </ul>
+ *
+ */
+ unsigned int symmetric;
+
+ unsigned int n_levels;
+
+ double threshold;
+
+ double filter;
+
+ double load_bal;
+
+ /**
+ * Setting this flag to true
+ * produces output from HYPRE,
+ * when the preconditioner
+ * is constructed.
+ */
+ bool output_details;
+ };
+
+
+
+ /**
+ * Empty Constructor. You need to call
+ * initialize() before using this
+ * object.
+ */
+ PreconditionParaSails ();
+
+ /**
+ * Constructor. Take the matrix which
+ * is used to form the preconditioner,
+ * and additional flags if there are
+ * any.
+ */
+ PreconditionParaSails (const MatrixBase &matrix,
+ const AdditionalData &additional_data = AdditionalData());
+
+ /**
+ * Initializes the preconditioner
+ * object and calculate all data that
+ * is necessary for applying it in a
+ * solver. This function is
+ * automatically called when calling
+ * the constructor with the same
+ * arguments and is only used if you
+ * create the preconditioner without
+ * arguments.
+ */
+ void initialize (const MatrixBase &matrix,
+ const AdditionalData &additional_data = AdditionalData());
+
+ protected:
+ /**
+ * Store a copy of the flags for this
+ * particular preconditioner.
+ */
+ AdditionalData additional_data;
+ };
+
+
+
+/**
+ * A class that implements the interface to use the scalable implementation
+ * of the Parallel ILU algorithm in the Euclid library from the HYPRE suite.
+ * Note that PETSc has to be configured with HYPRE (e.g. with --download-hypre=1).
+ *
+ * @ingroup PETScWrappers
+ * @author Martin Steigemann, 2011
+ */
+ class PreconditionEuclid : public PreconditionerBase
+ {
+ public:
+ /**
+ * Standardized data struct to
+ * pipe additional flags to the
+ * preconditioner.
+ */
+ struct AdditionalData
+ {
+ /**
+ * Constructor.
+ */
+ AdditionalData (
+ const unsigned int level = 1,
+ const bool use_block_jacobi = false,
+ const bool output_details = false
+ );
+
+ /**
+ * Factorization level for ILU(k). Default is 1.
+ * For 2D convection-diffusion and similar problems,
+ * fastest solution time is typically obtained with
+ * levels 4 through 8. For 3D problems, fastest solution
+ * time is typically obtained with level 1.
+ */
+ unsigned int level;
+
+ /**
+ * Use Block Jacobi ILU preconditioning
+ * instead of PILU. Default is false.
+ * If subdomains contain relatively few nodes
+ * (less than 1000), or the problem is not
+ * well partitioned, Block Jacobi ILU
+ * may give faster solution time then PILU.
+ */
+ unsigned int use_block_jacobi;
+
+ /**
+ * Setting this flag to true
+ * produces debug output from
+ * HYPRE, when the preconditioner
+ * is constructed.
+ */
+ bool output_details;
+ };
+
+
+
+ /**
+ * Empty Constructor. You need to call
+ * initialize() before using this
+ * object.
+ */
+ PreconditionEuclid ();
+
+ /**
+ * Constructor. Take the matrix which
+ * is used to form the preconditioner,
+ * and additional flags if there are
+ * any.
+ */
+ PreconditionEuclid (const MatrixBase &matrix,
+ const AdditionalData &additional_data = AdditionalData());
+
+ /**
+ * Initializes the preconditioner
+ * object and calculate all data that
+ * is necessary for applying it in a
+ * solver. This function is
+ * automatically called when calling
+ * the constructor with the same
+ * arguments and is only used if you
+ * create the preconditioner without
+ * arguments.
+ */
+ void initialize (const MatrixBase &matrix,
+ const AdditionalData &additional_data = AdditionalData());
+
+ protected:
+ /**
+ * Store a copy of the flags for this
+ * particular preconditioner.
+ */
+ AdditionalData additional_data;
+ };
+
+
+
+/**
+ * A class that implements a non-preconditioned Krylov method.
+ *
+ * @ingroup PETScWrappers
+ * @author Martin Steigemann, 2011
+ */
+ class PreconditionNone : public PreconditionerBase
+ {
+ public:
+ /**
+ * Standardized data struct to
+ * pipe additional flags to the
+ * preconditioner.
+ */
+ struct AdditionalData
+ {};
+
+ /**
+ * Empty Constructor. You need to call
+ * initialize() before using this
+ * object.
+ */
+ PreconditionNone ();
+
+ /**
+ * Constructor. Take the matrix which
+ * is used to form the preconditioner,
+ * and additional flags if there are
+ * any.
+ */
+ PreconditionNone (const MatrixBase &matrix,
+ const AdditionalData &additional_data = AdditionalData());
+
+ /**
+ * Initializes the preconditioner
+ * object and calculate all data that
+ * is necessary for applying it in a
+ * solver. This function is
+ * automatically called when calling
+ * the constructor with the same
+ * arguments and is only used if you
+ * create the preconditioner without
+ * arguments.
+ */
+ void initialize (const MatrixBase &matrix,
+ const AdditionalData &additional_data = AdditionalData());
+
+ protected:
+ /**
+ * Store a copy of the flags for this
+ * particular preconditioner.
+ */
+ AdditionalData additional_data;
+ };
}
const PreconditionerBase &preconditioner);
+
+ void
+ solve (const MatrixBase &A,
+ VectorBase &x,
+ const VectorBase &b,
+ const PreconditionerBase &preconditioner,
+ const std::vector<VectorBase> &nullspace);
+
/**
* Resets the contained preconditioner
* and solver object. See class
<< "An error with error number " << arg1
<< " occurred while calling a PETSc function");
+ DeclException1 (ExcPETScSolverError,
+ char*,
+ << "PETSc solver failed: " << arg1);
+
protected:
/**
}
+
+template<>
+QGaussOneOverRThreeHalfs<2>::QGaussOneOverRThreeHalfs(const unsigned int n,
+ const unsigned int vertex_index,
+ const bool factor_out_singularity) :
+ Quadrature<2>(2*n*n)
+{
+ // This version of the constructor
+ // works only for the 4
+ // vertices. If you need a more
+ // general one, you should use the
+ // one with the Point<2> in the
+ // constructor.
+ Assert(vertex_index <4, ExcIndexRange(vertex_index, 0, 4));
+
+ // Start with the gauss quadrature
+ // formula on the (u,v) reference
+ // element.
+ QGauss<2> gauss(n);
+
+ Assert(gauss.size() == n*n, ExcInternalError());
+ Assert(vertex_index < 4, ExcIndexRange(vertex_index, 0, 4));
+
+ // We create only the first one. All other pieces are rotation of
+ // this one.
+ // In this case the transformation is
+ //
+ // (x,y) = (u*u, u*u tan(pi/4 v))
+ //
+ // with Jacobian
+ //
+ // J = 2 pi/4 R*R sqrt(cos(pi/4 v))
+ //
+ // And we get rid of R to take into account the singularity,
+ // unless specified differently in the constructor.
+ std::vector<Point<2> > &ps = this->quadrature_points;
+ std::vector<double> &ws = this->weights;
+ double pi4 = numbers::PI/4;
+
+ for(unsigned int q=0; q<gauss.size(); ++q) {
+ const Point<2> &gp = gauss.point(q);
+ ps[q][0] = gp[0]*gp[0];
+ ps[q][1] = gp[0]*gp[0]*std::tan(pi4 *gp[1]);
+ ws[q] = 2.*gauss.weight(q)*pi4/std::sqrt(std::cos(pi4 *gp[1]));
+ if(factor_out_singularity) {
+ const double abs_value = (ps[q]-GeometryInfo<2>::unit_cell_vertex(0)).norm();
+ ws[q] *= abs_value;
+ ws[q] *= std::sqrt(abs_value);
+ }
+ // The other half of the quadrilateral is symmetric with
+ // respect to xy plane.
+ ws[gauss.size()+q] = ws[q];
+ ps[gauss.size()+q][0] = ps[q][1];
+ ps[gauss.size()+q][1] = ps[q][0];
+ }
+
+ // Now we distribute these vertices in the correct manner
+ double theta = 0;
+ switch(vertex_index) {
+ case 0:
+ theta = 0;
+ break;
+ case 1:
+ theta = numbers::PI/2;
+ break;
+ case 2:
+ theta = -numbers::PI/2;
+ break;
+ case 3:
+ theta = numbers::PI;
+ break;
+ }
+
+ double R00 = std::cos(theta), R01 = -std::sin(theta);
+ double R10 = std::sin(theta), R11 = std::cos(theta);
+
+ if(vertex_index != 0)
+ for(unsigned int q=0; q<size(); ++q) {
+ double x = ps[q][0]-.5, y = ps[q][1]-.5;
+
+ ps[q][0] = R00*x + R01*y + .5;
+ ps[q][1] = R10*x + R11*y + .5;
+ }
+}
+
+
+
// construct the quadrature formulae in higher dimensions by
// tensor product of lower dimensions
}
+void TableHandler::write_results (std::ostream &out, const bool with_header) const
+{
+ bool math_mode = true;
+
+ AssertThrow (out, ExcIO());
+ if (with_header)
+ out << "\\documentclass[10pt]{report}" << std::endl
+ << "\\usepackage{float}" << std::endl
+ << "\\usepackage{booktabs}" << std::endl
+ << "\\usepackage{amsmath}" << std::endl << std::endl
+ << "\\begin{document}" << std::endl;
+
+ out << "\\begin{table}" << std::endl
+ << "\\begin{center}" << std::endl
+ << "\\begin{tabular}{@{}";
+
+ std::vector<std::string> sel_columns;
+ get_selected_columns(sel_columns);
+
+ // write the column formats
+ for (unsigned int j=0; j<column_order.size(); ++j)
+ {
+ std::string key=column_order[j];
+ // avoid `supercolumns[key]'
+ const std::map<std::string, std::vector<std::string> >::const_iterator
+ super_iter=supercolumns.find(key);
+
+ if (super_iter!=supercolumns.end())
+ {
+ const unsigned int n_subcolumns=super_iter->second.size();
+ for (unsigned int k=0; k<n_subcolumns; ++k)
+ {
+ // avoid `columns[supercolumns[key]]'
+ const std::map<std::string, Column>::const_iterator
+ col_iter=columns.find(super_iter->second[k]);
+ Assert(col_iter!=columns.end(), ExcInternalError());
+
+ out << col_iter->second.tex_format;
+ }
+ }
+ else
+ {
+ // avoid `columns[key]';
+ const std::map<std::string, Column>::const_iterator
+ col_iter=columns.find(key);
+ Assert(col_iter!=columns.end(), ExcInternalError());
+ out << col_iter->second.tex_format;
+ }
+ }
+ out << "@{}} \\toprule" << std::endl;
+
+ // write the caption line of the table
+ for (unsigned int j=0; j<column_order.size(); ++j)
+ {
+ std::string key=column_order[j];
+ const std::map<std::string, std::vector<std::string> >::const_iterator
+ super_iter=supercolumns.find(key);
+
+ if (super_iter!=supercolumns.end())
+ {
+ const unsigned int n_subcolumns=super_iter->second.size();
+ // avoid use of `tex_supercaptions[key]'
+ std::map<std::string,std::string>::const_iterator
+ tex_super_cap_iter=tex_supercaptions.find(key);
+ out << std::endl << "\\multicolumn{" << n_subcolumns << "}{c}{"
+ << tex_super_cap_iter->second << "}";
+ }
+ else
+ {
+ // col_iter->second=columns[col];
+ const std::map<std::string, Column>::const_iterator
+ col_iter=columns.find(key);
+ Assert(col_iter!=columns.end(), ExcInternalError());
+ out << col_iter->second.tex_caption;
+ }
+ if (j<column_order.size()-1)
+ out << " & ";
+ }
+ out << "\\\\ \\midrule" << std::endl;
+
+ // write the n rows
+ const unsigned int nrows=n_rows();
+ for (unsigned int i=0; i<nrows; ++i)
+ {
+ const unsigned int n_cols=sel_columns.size();
+
+ for (unsigned int j=0; j<n_cols; ++j)
+ {
+ std::string key=sel_columns[j];
+ // avoid `column[key]'
+ const std::map<std::string, Column>::const_iterator
+ col_iter=columns.find(key);
+ Assert(col_iter!=columns.end(), ExcInternalError());
+
+ const Column &column=col_iter->second;
+
+ out << std::setprecision(column.precision);
+
+ if (col_iter->second.scientific)
+ out.setf(std::ios::scientific, std::ios::floatfield);
+ else
+ out.setf(std::ios::fixed, std::ios::floatfield);
+
+ if (math_mode)
+ out << "$";
+
+ out << column.entries[i].value;
+
+ if (math_mode)
+ out << "$";
+
+ if (j<n_cols-1)
+ out << " & ";
+ }
+
+ out << (i == nrows-1 ? "\\\\ \\bottomrule" : "\\\\") << std::endl;
+ }
+
+ out << "\\end{tabular}" << std::endl
+ << "\\end{center}" << std::endl;
+ if(tex_table_caption!="")
+ out << "\\caption{" << tex_table_caption << "}" << std::endl;
+ if(tex_table_label!="")
+ out << "\\caption{" << tex_table_label << "}" << std::endl;
+ out << "\\end{table}" << std::endl;
+ if (with_header)
+ out << "\\end{document}" << std::endl;
+}
+
+
unsigned int TableHandler::n_rows() const
{
if (columns.size() == 0)
template <int dim, int spacedim>
void
Triangulation<dim,spacedim>::
- copy_triangulation (const dealii::Triangulation<dim, spacedim> &)
+ copy_triangulation (const dealii::Triangulation<dim, spacedim> &old_tria)
{
- Assert (false, ExcNotImplemented());
+ clear();
+
+ try
+ {
+ dealii::Triangulation<dim,spacedim>::
+ copy_triangulation (old_tria);
+ }
+ catch (const typename dealii::Triangulation<dim,spacedim>::DistortedCellList &)
+ {
+ // the underlying
+ // triangulation should not
+ // be checking for
+ // distorted cells
+ AssertThrow (false, ExcInternalError());
+ }
+
+ // note that now we have some content in
+ // the p4est objects and call the
+ // functions that do the actual work
+ // (which are dimension dependent, so
+ // separate)
+ triangulation_has_content = true;
+
+ Assert (old_tria.n_levels() == 1,
+ ExcMessage ("Parallel distributed triangulations can only be copied, "
+ "if they are not refined!"));
+
+ if (dynamic_cast<const dealii::parallel::distributed::Triangulation<dim,spacedim> *>(&old_tria) != 0)
+ {
+ Assert (!(dynamic_cast<const parallel::distributed::Triangulation<dim,spacedim>&>
+ (old_tria).refinement_in_progress),
+ ExcMessage ("Parallel distributed triangulations can only "
+ "be copied, if no refinement is in progress!"));
+
+ coarse_cell_to_p4est_tree_permutation =
+ dynamic_cast<const parallel::distributed::Triangulation<dim,spacedim>&>
+ (old_tria).coarse_cell_to_p4est_tree_permutation;
+
+ p4est_tree_to_coarse_cell_permutation =
+ dynamic_cast<const parallel::distributed::Triangulation<dim,spacedim>&>
+ (old_tria).p4est_tree_to_coarse_cell_permutation;
+
+ attached_data_size =
+ dynamic_cast<const parallel::distributed::Triangulation<dim,spacedim>&>
+ (old_tria).attached_data_size;
+
+ n_attached_datas =
+ dynamic_cast<const parallel::distributed::Triangulation<dim,spacedim>&>
+ (old_tria).n_attached_datas;
+ }
+ else
+ {
+ setup_coarse_cell_to_p4est_tree_permutation ();
+ };
+
+ copy_new_triangulation_to_p4est (dealii::internal::int2type<dim>());
+
+ try
+ {
+ copy_local_forest_to_triangulation ();
+ }
+ catch (const typename Triangulation<dim>::DistortedCellList &)
+ {
+ // the underlying
+ // triangulation should not
+ // be checking for
+ // distorted cells
+ AssertThrow (false, ExcInternalError());
+ }
+
+ update_number_cache ();
}
}
+/* ----------------- PreconditionParaSails -------------------- */
+
+ PreconditionParaSails::AdditionalData::
+ AdditionalData(const unsigned int symmetric,
+ const unsigned int n_levels,
+ const double threshold,
+ const double filter,
+ const double load_bal,
+ const bool output_details)
+ :
+ symmetric(symmetric),
+ n_levels(n_levels),
+ threshold(threshold),
+ filter(filter),
+ load_bal(load_bal),
+ output_details(output_details)
+ {}
+
+
+ PreconditionParaSails::PreconditionParaSails ()
+ {}
+
+
+ PreconditionParaSails::PreconditionParaSails (const MatrixBase &matrix,
+ const AdditionalData &additional_data)
+ {
+ initialize(matrix, additional_data);
+ }
+
+
+ void
+ PreconditionParaSails::initialize (const MatrixBase &matrix_,
+ const AdditionalData &additional_data_)
+ {
+ matrix = static_cast<Mat>(matrix_);
+ additional_data = additional_data_;
+
+#ifdef PETSC_HAVE_HYPRE
+ create_pc();
+
+ int ierr;
+ ierr = PCSetType (pc, const_cast<char *>(PCHYPRE));
+ AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+ ierr = PCHYPRESetType(pc, "parasails");
+ AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+ if (additional_data.output_details)
+ PetscOptionsSetValue("-pc_hypre_parasails_logging","1");
+
+ Assert ((additional_data.symmetric == 0 ||
+ additional_data.symmetric == 1 ||
+ additional_data.symmetric == 2),
+ ExcMessage("ParaSails parameter symmetric can only be equal to 0, 1, 2!"));
+
+ std::stringstream ssStream;
+
+ switch (additional_data.symmetric)
+ {
+ case 0:
+ {
+ ssStream << "nonsymmetric";
+ break;
+ }
+
+ case 1:
+ {
+ ssStream << "SPD";
+ break;
+ }
+
+ case 2:
+ {
+ ssStream << "nonsymmetric,SPD";
+ break;
+ }
+
+ default:
+ Assert (false,
+ ExcMessage("ParaSails parameter symmetric can only be equal to 0, 1, 2!"));
+ };
+
+ PetscOptionsSetValue("-pc_hypre_parasails_sym",ssStream.str().c_str());
+
+ PetscOptionsSetValue("-pc_hypre_parasails_nlevels",
+ Utilities::int_to_string(
+ additional_data.n_levels
+ ).c_str());
+
+ ssStream.str(""); // empty the stringstream
+ ssStream << additional_data.threshold;
+ PetscOptionsSetValue("-pc_hypre_parasails_thresh", ssStream.str().c_str());
+
+ ssStream.str(""); // empty the stringstream
+ ssStream << additional_data.filter;
+ PetscOptionsSetValue("-pc_hypre_parasails_filter", ssStream.str().c_str());
+
+ ssStream.str(""); // empty the stringstream
+ ssStream << additional_data.load_bal;
+ PetscOptionsSetValue("-pc_hypre_parasails_loadbal", ssStream.str().c_str());
+
+ ierr = PCSetFromOptions (pc);
+ AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+ ierr = PCSetUp (pc);
+ AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+#else // PETSC_HAVE_HYPRE
+ (void)pc;
+ Assert (false,
+ ExcMessage ("Your PETSc installation does not include a copy of "
+ "the hypre package necessary for this preconditioner."));
+#endif
+ }
+
+
+/* ----------------- PreconditionEuclid ----------------------- */
+
+ PreconditionEuclid::AdditionalData::
+ AdditionalData(const unsigned int level,
+ const bool use_block_jacobi,
+ const bool output_details)
+ :
+ level(level),
+ use_block_jacobi(use_block_jacobi),
+ output_details(output_details)
+ {}
+
+
+ PreconditionEuclid::PreconditionEuclid ()
+ {}
+
+
+ PreconditionEuclid::PreconditionEuclid (const MatrixBase &matrix,
+ const AdditionalData &additional_data)
+ {
+ initialize(matrix, additional_data);
+ }
+
+
+ void
+ PreconditionEuclid::initialize (const MatrixBase &matrix_,
+ const AdditionalData &additional_data_)
+ {
+ matrix = static_cast<Mat>(matrix_);
+ additional_data = additional_data_;
+
+#ifdef PETSC_HAVE_HYPRE
+ create_pc();
+
+ int ierr;
+ ierr = PCSetType (pc, const_cast<char *>(PCHYPRE));
+ AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+ ierr = PCHYPRESetType(pc, "euclid");
+ AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+ PetscOptionsSetValue("-pc_hypre_euclid_levels",
+ Utilities::int_to_string(
+ additional_data.level
+ ).c_str());
+
+ if (additional_data.use_block_jacobi)
+ PetscOptionsSetValue("-pc_hypre_euclid_bj","1");
+
+ if (additional_data.output_details)
+ PetscOptionsSetValue("-pc_hypre_euclid_print_statistics","1");
+
+ ierr = PCSetFromOptions (pc);
+ AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+ ierr = PCSetUp (pc);
+ AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+#else // PETSC_HAVE_HYPRE
+ (void)pc;
+ Assert (false,
+ ExcMessage ("Your PETSc installation does not include a copy of "
+ "the hypre package necessary for this preconditioner."));
+#endif
+ }
+
+
+/* ----------------- PreconditionNone ------------------------- */
+
+ PreconditionNone::PreconditionNone ()
+ {}
+
+
+ PreconditionNone::PreconditionNone (const MatrixBase &matrix,
+ const AdditionalData &additional_data)
+ {
+ initialize(matrix, additional_data);
+ }
+
+
+ void
+ PreconditionNone::initialize (const MatrixBase &matrix_,
+ const AdditionalData &additional_data_)
+ {
+ matrix = static_cast<Mat>(matrix_);
+ additional_data = additional_data_;
+
+ create_pc();
+
+ int ierr;
+ ierr = PCSetType (pc, const_cast<char *>(PCNONE));
+ AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+ ierr = PCSetFromOptions (pc);
+ AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+ ierr = PCSetUp (pc);
+ AssertThrow (ierr == 0, ExcPETScError(ierr));
+ }
+
+
/* ----------------- PreconditionLU -------------------- */
PreconditionLU::AdditionalData::
}
+ void
+ SolverBase::solve (const MatrixBase &A,
+ VectorBase &x,
+ const VectorBase &b,
+ const PreconditionerBase &preconditioner,
+ const std::vector<VectorBase> &nullspace)
+ {
+ int ierr;
+ // first create a solver object if this
+ // is necessary
+ if (solver_data.get() == 0)
+ {
+ solver_data.reset (new SolverData());
+
+ ierr = KSPCreate (mpi_communicator, &solver_data->ksp);
+ AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+ // set the matrices involved. the
+ // last argument is irrelevant here,
+ // since we use the solver only once
+ // anyway
+ ierr = KSPSetOperators (solver_data->ksp, A, preconditioner,
+ SAME_PRECONDITIONER);
+ AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+ // let derived classes set the solver
+ // type, and the preconditioning
+ // object set the type of
+ // preconditioner
+ set_solver_type (solver_data->ksp);
+
+ ierr = KSPSetPC (solver_data->ksp, preconditioner.get_pc());
+ AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+ PC prec;
+ ierr = KSPGetPC (solver_data->ksp, &prec);
+ ierr = PCFactorSetShiftType (prec, MAT_SHIFT_POSITIVE_DEFINITE);
+
+ // then a convergence monitor
+ // function. that function simply
+ // checks with the solver_control
+ // object we have in this object for
+ // convergence
+#if DEAL_II_PETSC_VERSION_LT(3,0,0)
+ KSPSetConvergenceTest (solver_data->ksp, &convergence_test,
+ reinterpret_cast<void *>(&solver_control));
+#else
+ KSPSetConvergenceTest (solver_data->ksp, &convergence_test,
+ reinterpret_cast<void *>(&solver_control),
+ PETSC_NULL);
+#endif
+
+ }
+
+ // then do the real work: set up solver
+ // internal data and solve the
+ // system.
+ ierr = KSPSetUp (solver_data->ksp);
+ AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+ unsigned int dim_nullsp = nullspace.size();
+
+ if (dim_nullsp > 0)
+ {
+ Vec * nullsp_basis;
+
+ PetscMalloc (dim_nullsp*sizeof(Vec), &nullsp_basis);
+
+ for (unsigned int i=0; i<dim_nullsp; ++i)
+ nullsp_basis[i] = &(*nullspace[i]);
+
+ MatNullSpace nullsp;
+
+ ierr = MatNullSpaceCreate (mpi_communicator, PETSC_FALSE, dim_nullsp, nullsp_basis, &nullsp);
+ AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+ PetscBool is_nullspace;
+
+ ierr = MatNullSpaceTest (nullsp, A, &is_nullspace);
+ AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+ ierr = KSPSetNullSpace (solver_data->ksp, nullsp);
+ AssertThrow (ierr == 0, ExcPETScError(ierr));
+ };
+
+ ierr = KSPSolve (solver_data->ksp, b, x);
+ AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+ KSPConvergedReason reason;
+ ierr = KSPGetConvergedReason (solver_data->ksp, &reason);
+ AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+ if (reason < 0)
+ {
+ std::ostringstream error_code;
+
+ if (reason == KSP_DIVERGED_NULL)
+ error_code << "KSP_DIVERGED_NULL";
+ else if (reason == KSP_DIVERGED_ITS)
+ error_code << "KSP_DIVERGED_ITS";
+ else if (reason == KSP_DIVERGED_ITS)
+ error_code << "KSP_DIVERGED_DTOL";
+ else if (reason == KSP_DIVERGED_BREAKDOWN)
+ error_code << "KSP_DIVERGED_BREAKDOWN";
+ else if (reason == KSP_DIVERGED_BREAKDOWN_BICG)
+ error_code << "KSP_DIVERGED_BREAKDOWN_BICG";
+ else if (reason == KSP_DIVERGED_NONSYMMETRIC)
+ error_code << "KSP_DIVERGED_NONSYMMETRIC";
+ else if (reason == KSP_DIVERGED_INDEFINITE_PC)
+ error_code << "KSP_DIVERGED_INDEFINITE_PC";
+ else if (reason == KSP_DIVERGED_NAN)
+ error_code << "KSP_DIVERGED_NAN";
+ else if (reason == KSP_DIVERGED_INDEFINITE_MAT)
+ error_code << "KSP_DIVERGED_INDEFINITE_MAT";
+ else
+ error_code << "Unknown Error";
+
+ AssertThrow (false, ExcPETScSolverError(error_code.str().c_str()));
+ };
+
+ // do not destroy solver object
+// solver_data.reset ();
+
+ // in case of failure: throw
+ // exception
+ if (solver_control.last_check() != SolverControl::success)
+ throw SolverControl::NoConvergence (solver_control.last_step(),
+ solver_control.last_value());
+ // otherwise exit as normal
+ }
+
+
void
SolverBase::set_prefix(const std::string &prefix)
{
// honor the initial guess in the
// solution vector. do so here as well:
KSPSetInitialGuessNonzero (ksp, PETSC_TRUE);
+
+ KSPSetTolerances(ksp, PETSC_DEFAULT, this->solver_control.tolerance(),
+ PETSC_DEFAULT, this->solver_control.max_steps()+1);
}
// honor the initial guess in the
// solution vector. do so here as well:
KSPSetInitialGuessNonzero (ksp, PETSC_TRUE);
+
+ KSPSetTolerances(ksp, PETSC_DEFAULT, this->solver_control.tolerance(),
+ PETSC_DEFAULT, this->solver_control.max_steps()+1);
}
projects / dealii-svn.git / commitdiff