class VectorBase;
class SparseMatrix;
+ class BlockSparseMatrix;
class SolverBase;
/**
<< ". Check preconditioner and matrix setup.");
friend class SolverBase;
- friend class SolverSaddlePoint;
+ friend class PreconditionStokes;
- protected:
+ //protected:
/**
* This is a pointer to the
* preconditioner object that
--- /dev/null
+//---------------------------------------------------------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2008 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.
+//
+//---------------------------------------------------------------------------
+#ifndef __deal2__trilinos_precondition_block_h
+#define __deal2__trilinos_precondition_block_h
+
+
+#include <base/config.h>
+#include <lac/exceptions.h>
+#include <lac/trilinos_sparse_matrix.h>
+#include <lac/trilinos_block_sparse_matrix.h>
+#include <lac/trilinos_precondition.h>
+
+
+#ifdef DEAL_II_USE_TRILINOS
+
+#include <Teuchos_RefCountPtr.hpp>
+#include <Thyra_LinearOperatorDecl.hpp>
+#include <Epetra_Operator.h>
+
+// some forward declarations
+class Ifpack_Preconditioner;
+
+
+DEAL_II_NAMESPACE_OPEN
+
+/*! @addtogroup TrilinosWrappers
+ *@{
+ */
+
+namespace TrilinosWrappers
+{
+
+/**
+ * This implements an interface class of preconditioners for block
+ * matrices based on Trilinos, which are intended to be used with the
+ * solver classes SolverBlockCG and SolverBlockGMRES. This is based on
+ * the Trilinos Thyra abstract interface of linear operators. This
+ * class does not implement any method, and derived classes will need
+ * to do that job.
+ *
+ * @ingroup TrilinosWrappers
+ * @ingroup Preconditioners
+ * @author Martin Kronbichler, 2008
+ */
+ class PreconditionBlockBase : public Subscriptor
+ {
+ public:
+ /**
+ * Empty constructor.
+ */
+ PreconditionBlockBase();
+
+ /**
+ * Destructor.
+ */
+ ~PreconditionBlockBase();
+
+ protected:
+ /**
+ * Preconditioner object.
+ */
+ Teuchos::RCP<Thyra::ConstLinearOperator<double> > preconditioner;
+
+ friend class SolverBlockBase;
+};
+
+ /**
+ * Internal function that sets
+ * up an inverse matrix.
+ */
+ inline
+ Thyra::ConstLinearOperator<double>
+ inverse_matrix (const SparseMatrix &M,
+ const Epetra_Operator *P,
+ const bool is_symmetric,
+ const unsigned int n_iterations,
+ const double solve_tolerance,
+ const bool output_details);
+
+
+ // Forward declarations.
+
+/**
+ * This class implements a black box preconditioner for saddle points
+ * systems arising from the Stokes or Navier–Stokes equations as
+ * specified by the papers <em>D. Silvester, A. Wathen, Fast iterative
+ * solution of stabilised Stokes systems part II. Using general block
+ * preconditioners, SIAM J. Numer. Anal. 31:1352&ndash1367 (1994)</em>
+ * and <em> D. Kay, D. Loghin, A. Wathen, A preconditioner for the
+ * steady-state Navier–Stokes equations, SIAM
+ * J. Sci. Comput. 24(1):237–256 (2002)</em>, respectively.
+ *
+ * @ingroup TrilinosWrappers
+ * @ingroup Preconditioners
+ * @author Martin Kronbichler, 2008
+ */
+ class PreconditionStokes : public PreconditionBlockBase
+ {
+ public:
+ /**
+ * Standardized data struct to
+ * pipe additional data to the
+ * solver.
+ */
+ struct AdditionalData
+ {
+ /**
+ * Constructor. This sets the
+ * additional data to its
+ * default values. For lazy
+ * initializations, the user
+ * will just set the parameters
+ * in here and let the
+ * <tt>initialize</tt> function
+ * set up good preconditioners,
+ * where some good values have
+ * to be provided here. The
+ * default ones mean that we
+ * assume the velocity-velocity
+ * matrix to be symmetric (only
+ * true for Stokes problems!),
+ * we use a relative inner
+ * tolerance of 1e-9 for the
+ * inversion of matrices, and
+ * take in IC preconditioner
+ * for the inversion of the
+ * pressure mass
+ * matrix. Moreover, there is
+ * no inner CG/GMRES solve
+ * performed on the
+ * velocity-velocity block by
+ * default. Otherwise, we do a
+ * solve with the specified
+ * solver and the number of
+ * iterations and
+ * tolerance. Note that too
+ * many iterations on Av can
+ * slow down the overall
+ * procedure considerably.
+ */
+ AdditionalData (const bool right_preconditioning = false,
+ const bool Av_is_symmetric = true,
+ const std::vector<std::vector<bool> > &Av_null_space = std::vector<std::vector<bool> > (1),
+ const double inner_solve_tolerance = 1e-9,
+ const bool use_ssor_on_mass_matrix = false,
+ const bool velocity_uses_higher_order_elements = false,
+ const bool pressure_uses_higher_order_elements = false,
+ const bool Av_do_solve = false,
+ const unsigned int Av_n_iters = 1,
+ const double Av_tolerance = 1e-3,
+ const bool ouput_details = false);
+
+ /**
+ * This flag specifies whether
+ * the preconditioner should be
+ * build as a left
+ * preconditioner or right
+ * preconditioner. Note that
+ * this setting must be the
+ * same as the one used in the
+ * solver call.
+ */
+ bool right_preconditioning;
+
+ /**
+ * This flag determines whether
+ * the velocity-velocity matrix
+ * is symmetric (and positive
+ * definite) or not. This
+ * choice will influence the
+ * AMG preconditioner that is
+ * built for that system as
+ * well as the inner solve on
+ * Av (if that flag is
+ * enabled).
+ */
+ bool Av_is_symmetric;
+
+ /**
+ * This determines the near
+ * null space for the operator
+ * Av, the vector-valued
+ * velocity-velocity space. In
+ * order to get good
+ * performance, this vector has
+ * to be specified using the
+ * DoFTools::extract_constant_modes()
+ * function.
+ */
+ std::vector<std::vector<bool> > Av_null_space;
+
+ /**
+ * Tolerance level that the
+ * inner solvers on the
+ * pressure mass matrix (and
+ * the pressure Laplace matrix
+ * in case of a Navier-Stokes
+ * problem) should be solve
+ * to. It is essential that
+ * this tolerance is in the
+ * same order of magnitude than
+ * the one chosen for the outer
+ * GMRES solver (or little
+ * less). Otherwise unexpected
+ * variations in the number of
+ * iterations can occur from
+ * solve to solve.
+ *
+ * Default value: <tt>1e-9</tt>
+ * residual relative to initial
+ * residual.
+ */
+ double inner_solve_tolerance;
+
+ /**
+ * Determines whether the inner
+ * solve for the pressure mass
+ * matrix should use an IC
+ * preconditioner (incomplete
+ * Cholesky factorization) or
+ * an SSOR preconditioner. The
+ * former tends to be faster in
+ * most situations, whereas the
+ * latter uses almost no
+ * additional memory besides
+ * the matrix storage.
+ */
+ bool use_ssor_on_mass_matrix;
+
+ /**
+ * Determines whether the
+ * underlying velocity
+ * discretization uses linear
+ * elements or higher order
+ * elements, which has
+ * consequences for the set up
+ * of the internal setup of the
+ * ML AMG preconditioner.
+ */
+ bool velocity_uses_higher_order_elements;
+
+ /**
+ * Determines whether the
+ * underlying pressure
+ * discretization uses linear
+ * elements or higher order
+ * elements, which has
+ * consequences for the set up
+ * of the internal setup of the
+ * ML AMG preconditioner.
+ */
+ bool pressure_uses_higher_order_elements;
+
+ /**
+ * Flag to determine whether we
+ * should do an inner solve on
+ * the velocity-velocity
+ * matrix. By default, this
+ * option is not on, since it
+ * is not necessary for a good
+ * performance of the outer
+ * solver. However, it can be
+ * beneficial to perform some
+ * iterations on this block
+ * only and not on the whole
+ * block, so that the number of
+ * outer iterations is
+ * reduced. See the discussion
+ * in the @step_22 step-22
+ * tutorial program.
+ */
+ unsigned int Av_do_solve;
+
+ /**
+ * The number of iterations in
+ * the inner solve.
+ */
+ unsigned int Av_n_iters;
+
+ /**
+ * The residual tolerance that
+ * is going to be used for the
+ * inner solve on the
+ * velocity-velocity matrix.
+ */
+ double Av_tolerance;
+
+ /**
+ * This defines whether
+ * internal details about the
+ * solve process should be
+ * written to screen. This can
+ * be a lot of information.
+ */
+ bool output_details;
+ };
+
+ /**
+ * Lazy setup of the block
+ * preconditioner for the
+ * (Navier-) Stokes
+ * system. This function takes
+ * the matrices given here and
+ * firs calculates good
+ * preconditioners, i.e.,
+ * algebraic multigrid
+ * preconditioners for the
+ * velocity-velocity matrix
+ * <tt>Av</tt>, that can be
+ * specified to be different
+ * from the (0,0) block in the
+ * system matrix, IC/SSOR on
+ * the pressure mass matrix
+ * <tt>Mp_matrix</tt>, and AMG
+ * on the pressure Laplace
+ * matrix
+ * <tt>Lp_matrix</tt>. Next,
+ * these preconditioners are
+ * used to generate a block
+ * operator.
+ */
+ void initialize (const BlockSparseMatrix &system_matrix,
+ const SparseMatrix &Mp_matrix,
+ const AdditionalData &additional_data,
+ const SparseMatrix &Av_matrix = SparseMatrix(),
+ const SparseMatrix &Lp_matrix = SparseMatrix(),
+ const SparseMatrix &Fp_matrix = SparseMatrix());
+
+ /**
+ * Set up the block
+ * preconditioner for the
+ * (Navier-) Stokes system. In
+ * contrast to the other
+ * <tt>initialize</tt>
+ * function, this one expects
+ * the user to specify
+ * preconditioner objects for
+ * the individual matrices,
+ * that will then be used for
+ * the inner inversion of the
+ * matrices when building the
+ * Schur complement
+ * preconditioner. If no
+ * preconditioner is specified,
+ * the inner solvers will use
+ * an identity preconditioner
+ * object.
+ */
+ void initialize (const BlockSparseMatrix &system_matrix,
+ const SparseMatrix &Mp_matrix,
+ const AdditionalData &additional_data,
+ const PreconditionBase &Mp_preconditioner = PreconditionBase(),
+ const PreconditionBase &Av_preconditioner = PreconditionBase(),
+ const SparseMatrix &Lp_matrix = SparseMatrix(),
+ const PreconditionBase &Lp_preconditioner = PreconditionBase(),
+ const SparseMatrix &Fp_matrix = SparseMatrix());
+
+ /**
+ * This function can be used
+ * for a faster recalculation
+ * of the preconditioner
+ * construction when the system
+ * matrix entries underlying
+ * the preconditioner have
+ * changed but not the sparsity
+ * pattern (this means that the
+ * grid has remained the same
+ * and the matrix structure has
+ * not been
+ * reinitialized). Moreover, it
+ * is assumed that the PDE
+ * parameters have not changed
+ * drastically. This function
+ * is only needed for the lazy
+ * variant of the
+ * preconditioner. In the other
+ * case, the preconditioner can
+ * be modified outside this
+ * function, which will be
+ * recongnized in this
+ * preconditioner as well when
+ * doing a solve, without the
+ * need to restructure anything
+ * here.
+ */
+ void reinit_lazy ();
+
+ /**
+ * Exception.
+ */
+ DeclException1 (ExcNonMatchingMaps,
+ std::string,
+ << "The sparse matrix the preconditioner is based on "
+ << "uses a map that is not compatible to the one in vector"
+ << arg1
+ << ". Check preconditioner and matrix setup.");
+
+ private:
+
+ /**
+ * Pointer to preconditioner
+ * for the mass matrix.
+ */
+ Teuchos::RCP<PreconditionSSOR> Mp_precondition_ssor;
+
+ /**
+ * Pointer to preconditioner
+ * for the mass matrix.
+ */
+ Teuchos::RCP<PreconditionIC> Mp_precondition_ic;
+
+ /**
+ * Pointer to preconditioner
+ * for the velocity-velocity
+ * matrix.
+ */
+ Teuchos::RCP<PreconditionAMG> Av_precondition;
+
+ /**
+ * Pointer to preconditioner
+ * for the pressure Laplace
+ * matrix.
+ */
+ Teuchos::RCP<PreconditionAMG> Lp_precondition;
+ };
+}
+
+/*@}*/
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif // DEAL_II_USE_TRILINOS
+
+/*-------------------------- trilinos_precondition_block.h ---------------------------*/
+
+#endif
+/*-------------------------- trilinos_precondition_block.h ---------------------------*/
* "http://trilinos.sandia.gov/packages/aztecoo/AztecOOUserGuide.pdf">AztecOO
* user guide</a>.
*
+ * This solver class can also be used as a standalone class, where the
+ * respective Krylov method is set via the flag
+ * <tt>solver_name</tt>. This can be done at runtime (e.g., when
+ * parsing the solver from a ParameterList) and is similar to the
+ * deal.II class SolverSelector.
+ *
* @ingroup TrilinosWrappers
* @author Martin Kronbichler, 2008
*/
{
public:
+ /**
+ * Enumeration object that is
+ * set in the constructor of
+ * the derived classes and
+ * tells Trilinos which solver
+ * to use. This option can also
+ * be set in the user program,
+ * so one might use this base
+ * class instead of one of the
+ * specialized derived classes
+ * when the solver should be
+ * set at runtime. Currently
+ * enabled options are:
+ */
+ enum SolverName {cg, cgs, gmres, bicgstab, tfqmr} solver_name;
+
/**
* Standardized data struct to
* pipe additional data to the
/**
* Constructor. Takes the
* solver control object and
- * the MPI communicator over
- * which parallel computations
- * are to happen.
+ * creates the solver.
*/
SolverBase (SolverControl &cn);
+
+ /**
+ * Second constructor. This
+ * constructor takes an enum
+ * object that specifies the
+ * solver name and sets the
+ * appropriate Krylov
+ * method.
+ */
+ SolverBase (const enum SolverName solver_name,
+ SolverControl &cn);
/**
* Destructor.
*/
SolverControl &solver_control;
- /**
- * String object that is set in
- * the constructor of the
- * derived classes and tells
- * Trilinos which solver to
- * use. Currently enabled
- * options are <tt>"cg", "cgs",
- * "gmres", "bicgstab"</tt>,
- * and <tt>"tfqmr"</tt>.
- */
- enum SolverName {cg, cgs, gmres, bicgstab, tfqmr} solver_name;
-
private:
/**
--- /dev/null
+//---------------------------------------------------------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2008 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.
+//
+//---------------------------------------------------------------------------
+#ifndef __deal2__trilinos_block_solver_h
+#define __deal2__trilinos_block_solver_h
+
+#include <base/config.h>
+#include <lac/exceptions.h>
+#include <lac/solver_control.h>
+
+#ifdef DEAL_II_USE_TRILINOS
+
+
+DEAL_II_NAMESPACE_OPEN
+
+
+namespace TrilinosWrappers
+{
+ // forward declarations
+ class BlockSparseMatrix;
+ class BlockVector;
+ namespace MPI
+ {
+ class BlockVector;
+ }
+ class SparseMatrix;
+ class PreconditionBase;
+ class PreconditionBlockBase;
+
+
+/**
+ * Base class for solver classes using the Trilinos solvers on block
+ * matrices based on the Trilinos abstract solver interface Thyra.
+ *
+ * @ingroup TrilinosWrappers
+ * @author Martin Kronbichler, 2008
+ */
+ class SolverBlockBase
+ {
+ public:
+
+ /**
+ * Standardized data struct to
+ * pipe additional data to the
+ * solver.
+ */
+ struct AdditionalData
+ {
+ /**
+ * Restart parameter in case
+ * the derived class is
+ * GMRES.
+ *
+ * TODO: Find a better way for
+ * doing this.
+ */
+ AdditionalData (const unsigned int gmres_restart_parameter = 30,
+ const bool right_preconditioning = false,
+ const bool output_details = false);
+
+ /**
+ * The number of vectors that
+ * should be used in the GMRES
+ * algorithm before the solver
+ * is restarted.
+ */
+ unsigned int gmres_restart_parameter;
+
+ /**
+ * Flag that turns right
+ * preconditioning on.
+ */
+ bool right_preconditioning;
+
+ /**
+ * A flag to determine whether
+ * solver details should be
+ * written to screen.
+ */
+ bool output_details;
+ };
+
+ /**
+ * Constructor. Takes the
+ * solver control object and
+ * the MPI communicator over
+ * which parallel computations
+ * are to happen.
+ */
+ SolverBlockBase (SolverControl &cn);
+
+ /**
+ * Destructor.
+ */
+ virtual ~SolverBlockBase ();
+
+ /**
+ * Solve the linear system
+ * <tt>Ax=b</tt>. Depending on
+ * the information provided by
+ * derived classes and the
+ * object passed as a
+ * preconditioner, one of the
+ * linear solvers and
+ * preconditioners of Trilinos
+ * is chosen. This interface is
+ * designed for Trilinos
+ * running in parallel.
+ */
+ void
+ solve (const BlockSparseMatrix &A,
+ MPI::BlockVector &x,
+ const MPI::BlockVector &b,
+ const PreconditionBlockBase &preconditioner);
+
+ /**
+ * Solve the linear system
+ * <tt>Ax=b</tt>. Depending on
+ * the information provided by
+ * derived classes and the
+ * object passed as a
+ * preconditioner, one of the
+ * linear solvers and
+ * preconditioners of Trilinos
+ * is chosen. This interface is
+ * designed for Trilinos
+ * running in serial.
+ */
+ void
+ solve (const BlockSparseMatrix &A,
+ BlockVector &x,
+ const BlockVector &b,
+ const PreconditionBlockBase &preconditioner);
+
+ /**
+ * Access to object that controls
+ * convergence.
+ */
+ SolverControl & control() const;
+
+ /**
+ * Exception
+ */
+ DeclException1 (ExcTrilinosError,
+ int,
+ << "An error with error number " << arg1
+ << " occured while calling a Trilinos function");
+
+ /**
+ * Exception
+ */
+ DeclException2 (ExcOverlappingMaps,
+ std::string, std::string,
+ << "The Epetra map in the " << arg1 << " " << arg2
+ << " is overlapping between the processors, which"
+ << " is forbidden when doing a solve.");
+
+ /**
+ * Exception.
+ */
+ DeclException1 (ExcNonMatchingMaps,
+ std::string,
+ << "The Epetra map in the vector "
+ << arg1
+ << " does not match the one used in the matrix. "
+ << "Check vector and matrix setup.");
+
+ protected:
+
+ /**
+ * Reference to the object that
+ * controls convergence of the
+ * iterative solver. In fact,
+ * for these Trilinos wrappers,
+ * Trilinos does so itself, but
+ * we copy the data from this
+ * object before starting the
+ * solution process, and copy
+ * the data back into it
+ * afterwards.
+ */
+ SolverControl &solver_control;
+
+ /**
+ * String object that is set in
+ * the constructor of the
+ * derived classes and tells
+ * Trilinos which solver to
+ * use. Currently enabled
+ * options are <tt>"cg",
+ * "gmres"</tt>.
+ */
+ enum SolverBlockName {cg, gmres} solver_name;
+
+ protected:
+
+ /**
+ * Store a copy of the flags for this
+ * particular solver.
+ */
+ AdditionalData additional_data;
+
+ };
+
+
+
+/**
+ * An implementation of the solver interface using the Trilinos CG
+ * solver for block matrices and vectors.
+ *
+ * @ingroup TrilinosWrappers
+ * @author Martin Kronbichler, 2008
+ */
+ class SolverBlockCG : public SolverBlockBase
+ {
+ public:
+
+ /**
+ * Standardized data struct to
+ * pipe additional data to the
+ * solver.
+ */
+ struct AdditionalData
+ {
+ /**
+ * Restart parameter in case
+ * the derived class is
+ * GMRES.
+ *
+ * TODO: Find a better way for
+ * doing this.
+ */
+ AdditionalData (const bool right_preconditioning = false,
+ const bool output_details = false);
+
+ /**
+ * Flag that turns right
+ * preconditioning on.
+ */
+ bool right_preconditioning;
+
+ /**
+ * A flag to determine whether
+ * solver details should be
+ * written to screen.
+ */
+ bool output_details;
+ };
+
+ /**
+ * Constructor. In contrast to
+ * deal.II's own solvers, there is no
+ * need to give a vector memory
+ * object.
+ *
+ * The last argument takes a structure
+ * with additional, solver dependent
+ * flags for tuning.
+ */
+ SolverBlockCG (SolverControl &cn,
+ const AdditionalData &data = AdditionalData());
+
+ protected:
+
+ /**
+ * Store a copy of the flags for this
+ * particular solver.
+ */
+ const AdditionalData data;
+ };
+
+
+
+/**
+ * An implementation of the solver interface using the Trilinos GMRES
+ * solver for block matrices and vectors.
+ *
+ * @author Martin Kronbichler, 2008
+ */
+ class SolverBlockGMRES : public SolverBlockBase
+ {
+ public:
+
+ /**
+ * Standardized data struct to
+ * pipe additional data to the
+ * solver.
+ */
+ struct AdditionalData
+ {
+ /**
+ * Restart parameter in case
+ * the derived class is
+ * GMRES.
+ *
+ * TODO: Find a better way for
+ * doing this.
+ */
+ AdditionalData (const unsigned int gmres_restart_parameter = 30,
+ const bool right_preconditioning = false,
+ const bool output_details = false);
+
+ /**
+ * The number of vectors that
+ * should be used in the GMRES
+ * algorithm before the solver
+ * is restarted.
+ */
+ unsigned int gmres_restart_parameter;
+
+ /**
+ * Flag that turns right
+ * preconditioning on.
+ */
+ bool right_preconditioning;
+
+ /**
+ * A flag to determine whether
+ * solver details should be
+ * written to screen.
+ */
+ bool output_details;
+ };
+
+ /**
+ * Constructor. In contrast to
+ * deal.II's own solvers, there is no
+ * need to give a vector memory
+ * object.
+ *
+ * The last argument takes a structure
+ * with additional, solver dependent
+ * flags for tuning.
+ */
+ SolverBlockGMRES (SolverControl &cn,
+ const AdditionalData &data = AdditionalData());
+
+ protected:
+
+ /**
+ * Store a copy of the flags for this
+ * particular solver.
+ */
+ const AdditionalData data;
+ };
+
+
+}
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif // DEAL_II_USE_TRILINOS
+
+/*---------------------- trilinos_solver_block.h -----------------------*/
+
+#endif
+/*---------------------- trilinos_solver_block.h -----------------------*/
ifpack.release();
ifpack = Teuchos::rcp (Ifpack().Create ("point relaxation",
- &*matrix.matrix,
- additional_data.overlap));
+ &*matrix.matrix,
+ additional_data.overlap));
Assert (&*ifpack != 0, ExcMessage ("Trilinos could not create this "
"preconditioner"));
additional_data.overlap));
Assert (&*ifpack != 0, ExcMessage ("Trilinos could not create this "
- "preconditioner"));
+ "preconditioner"));
int ierr;
additional_data.overlap));
Assert (&*ifpack != 0, ExcMessage ("Trilinos could not create this "
- "preconditioner"));
+ "preconditioner"));
int ierr;
--- /dev/null
+//---------------------------------------------------------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2008 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 <lac/trilinos_precondition_block.h>
+
+#ifdef DEAL_II_USE_TRILINOS
+
+#include <Thyra_EpetraLinearOp.hpp>
+#include <Thyra_EpetraThyraWrappers.hpp>
+
+#include <Thyra_VectorSpaceImpl.hpp>
+#include <Thyra_LinearOperatorDecl.hpp>
+#include <Thyra_DefaultBlockedLinearOpDecl.hpp>
+#include <Thyra_DefaultBlockedLinearOp.hpp>
+#include <Thyra_LinearOpWithSolveFactoryHelpers.hpp>
+#include <Thyra_InverseLinearOperator.hpp>
+#include <Thyra_DefaultInverseLinearOp.hpp>
+#include <Thyra_AztecOOLinearOpWithSolveFactory.hpp>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace TrilinosWrappers
+{
+
+ PreconditionBlockBase::PreconditionBlockBase()
+ {}
+
+
+
+ PreconditionBlockBase::~PreconditionBlockBase()
+ {}
+
+
+
+ /**
+ * Internal function that sets
+ * up an inverse matrix.
+ */
+ inline
+ Thyra::ConstLinearOperator<double>
+ inverse_matrix (const SparseMatrix &input_M,
+ const Epetra_Operator *input_P,
+ const bool is_symmetric,
+ const unsigned int n_iterations,
+ const double solve_tolerance,
+ const bool output_details)
+ {
+
+ // Create discrete
+ // Epetra operator
+ // (Thyra wrapper around the
+ // matrix).
+ Teuchos::RCP<const Thyra::LinearOpBase<double> > tmpM =
+ Teuchos::rcp(new Thyra::EpetraLinearOp(Teuchos::rcp(&*input_M.matrix,
+ false)));
+
+ Teuchos::RCP<const Thyra::LinearOpSourceBase<double> > M
+ = Teuchos::rcp(new Thyra::DefaultLinearOpSource<double>(tmpM), true);
+
+ //Thyra::ConstLinearOperator<double> opM = tmpM;
+
+ // Create a Thyra version of
+ // the preconditioner.
+ Teuchos::RCP<const Thyra::LinearOpBase<double> > tmpP =
+ Teuchos::rcp(new Thyra::EpetraLinearOp(Teuchos::rcp(const_cast<Epetra_Operator*>(input_P),
+ false),
+ Thyra::NOTRANS,
+ Thyra::EPETRA_OP_APPLY_APPLY_INVERSE,
+ Thyra::EPETRA_OP_ADJOINT_SUPPORTED));
+
+ //Thyra::ConstLinearOperator<double> opPM = tmpP;
+
+ Teuchos::RCP<Thyra::PreconditionerBase<double> > PM =
+ Teuchos::rcp (new Thyra::DefaultPreconditioner<double>(), true);
+
+ {
+ Thyra::DefaultPreconditioner<double>
+ *defaultPrec = &Teuchos::dyn_cast<Thyra::DefaultPreconditioner<double> >
+ (*PM);
+
+ (*defaultPrec).initializeRight(tmpP);
+ }
+
+ // Set up the solver.
+ Teuchos::RCP<Teuchos::ParameterList> aztecParams
+ = Teuchos::rcp(new Teuchos::ParameterList("aztecOOFSolverFactory"), true);
+
+ aztecParams->sublist("Forward Solve").set("Max Iterations",
+ (int)n_iterations);
+ aztecParams->sublist("Forward Solve").set("Tolerance",
+ solve_tolerance);
+
+ // aztecOO solver settings
+ if (is_symmetric == true)
+ {
+ aztecParams->sublist("Forward Solve")
+ .sublist("AztecOO Settings").set("Aztec Solver", "CG");
+ }
+ else
+ {
+ aztecParams->sublist("Forward Solve")
+ .sublist("AztecOO Settings").set("Aztec Solver", "GMRES");
+ aztecParams->sublist("Forward Solve")
+ .sublist("AztecOO Settings").set("Size of Krylov Subspace",
+ (int)n_iterations);
+ }
+
+ aztecParams->sublist("Forward Solve")
+ .sublist("AztecOO Settings").set("Aztec Preconditioner", "none");
+ aztecParams->sublist("Forward Solve")
+ .sublist("AztecOO Settings").set("Convergence Test", "r0");
+
+ // This sets if solver details
+ // (the residual in each
+ // iterations) should be
+ // written to screen.
+ if (output_details == true)
+ {
+ aztecParams->sublist("Forward Solve")
+ .sublist("AztecOO Settings").set("Output Frequency", 1);
+ }
+ else
+ {
+ aztecParams->sublist("VerboseObject").set("Verbosity Level", "none");
+ }
+
+ Teuchos::RCP<Thyra::LinearOpWithSolveFactoryBase<double> >
+ aztecLowsFactory = Teuchos::rcp(new Thyra::AztecOOLinearOpWithSolveFactory(), true);
+
+ aztecLowsFactory->setParameterList(aztecParams);
+
+ Teuchos::RCP<Thyra::LinearOpWithSolveBase<double> > rcpAztec
+ = aztecLowsFactory->createOp();
+
+ // Does not work like
+ // this. The preconditioner
+ // will only appear in the
+ // operator, but won't be set
+ // correctly to the right hand
+ // side vector, which makes
+ // the results wrong. So try
+ // something else...
+ //Meros::LinearSolveStrategy<double> azM
+ // = new Meros::AztecSolveStrategy(*(aztecParams.get()));
+ //Thyra::LinearOperator<double> Minverse
+ // = Thyra::inverse(*aztecLowsFactory, opM);
+ //return new Meros::InverseOperator<double>(opM * opPM, azM);
+
+ // Attach the preconditioner
+ // to the solver.
+ aztecLowsFactory->initializePreconditionedOp(M, PM, &*rcpAztec,
+ Thyra::SUPPORT_SOLVE_FORWARD_ONLY);
+
+ // Create an inverse matrix
+ // object and return.
+ Teuchos::RCP<const Thyra::LinearOpBase<double> > Minverse =
+ Teuchos::rcp(new Thyra::DefaultInverseLinearOp<double>(rcpAztec));
+
+ return Minverse;
+ }
+
+
+
+ PreconditionStokes::AdditionalData::
+ AdditionalData (const bool right_preconditioning,
+ const bool Av_is_symmetric,
+ const std::vector<std::vector<bool> > &Av_null_space,
+ const double inner_solve_tolerance,
+ const bool use_ssor_on_mass_matrix,
+ const bool velocity_uses_higher_order_elements,
+ const bool pressure_uses_higher_order_elements,
+ const bool Av_do_solve,
+ const unsigned int Av_n_iters,
+ const double Av_tolerance,
+ const bool output_details)
+ :
+ right_preconditioning (right_preconditioning),
+ Av_is_symmetric (Av_is_symmetric),
+ Av_null_space (Av_null_space),
+ inner_solve_tolerance (inner_solve_tolerance),
+ use_ssor_on_mass_matrix (use_ssor_on_mass_matrix),
+ velocity_uses_higher_order_elements (velocity_uses_higher_order_elements),
+ pressure_uses_higher_order_elements (pressure_uses_higher_order_elements),
+ Av_do_solve (Av_do_solve),
+ Av_n_iters (Av_n_iters),
+ Av_tolerance (Av_tolerance),
+ output_details (output_details)
+ {}
+
+
+
+ void
+ PreconditionStokes::initialize(const BlockSparseMatrix &system_matrix,
+ const SparseMatrix &Mp_matrix,
+ const AdditionalData &additional_data,
+ const SparseMatrix &Av_matrix,
+ const SparseMatrix &Lp_matrix,
+ const SparseMatrix &Fp_matrix)
+ {
+ AssertThrow (system_matrix.n_block_rows() == 2,
+ ExcDimensionMismatch(system_matrix.n_block_rows(),
+ 2));
+ AssertThrow (system_matrix.n_block_cols() == 2,
+ ExcDimensionMismatch(system_matrix.n_block_cols(),
+ 2));
+
+ // Create Mp preconditioner.
+ Teuchos::RCP<PreconditionBase> Mp_precondition;
+ {
+ if (additional_data.use_ssor_on_mass_matrix == true)
+ {
+ Mp_precondition_ssor = Teuchos::rcp (new PreconditionSSOR(), true);
+ Mp_precondition_ssor->initialize (Mp_matrix,
+ PreconditionSSOR::AdditionalData (1.2,0,0));
+ Mp_precondition = Teuchos::rcp(&*Mp_precondition_ssor, false);
+ }
+ else
+ {
+ Mp_precondition_ic = Teuchos::rcp (new PreconditionIC(), true);
+ Mp_precondition_ic->initialize (Mp_matrix);
+ Mp_precondition = Teuchos::rcp(&*Mp_precondition_ic, false);
+ }
+ }
+
+ // Create Av AMG
+ // preconditioner.
+ {
+ Av_precondition = Teuchos::rcp (new PreconditionAMG(), true);
+ PreconditionAMG::AdditionalData av_amg_data;
+ av_amg_data.elliptic = additional_data.Av_is_symmetric;
+ av_amg_data.higher_order_elements = additional_data.velocity_uses_higher_order_elements;
+ av_amg_data.null_space = additional_data.Av_null_space;
+
+ if (Av_matrix.m() == system_matrix.block(0,0).m())
+ Av_precondition->initialize (Av_matrix, av_amg_data);
+ else
+ Av_precondition->initialize (system_matrix.block(0,0), av_amg_data);
+ }
+
+ // Create Lp AMG
+ // preconditioner.
+ if (Lp_matrix.m() == system_matrix.block(1,1).m())
+ {
+ Lp_precondition = Teuchos::rcp (new PreconditionAMG(), true);
+ PreconditionAMG::AdditionalData lp_amg_data;
+ lp_amg_data.elliptic = true;
+ lp_amg_data.higher_order_elements = additional_data.pressure_uses_higher_order_elements;
+
+ Lp_precondition->initialize (Lp_matrix, lp_amg_data);
+ }
+
+ initialize (system_matrix, Mp_matrix, additional_data, *Mp_precondition,
+ *Av_precondition, Lp_matrix, *Lp_precondition, Fp_matrix);
+
+ }
+
+ void
+ PreconditionStokes::initialize (const BlockSparseMatrix &system_matrix,
+ const SparseMatrix &Mp_matrix,
+ const AdditionalData &additional_data,
+ const PreconditionBase &Mp_preconditioner,
+ const PreconditionBase &Av_preconditioner,
+ const SparseMatrix &Lp_matrix,
+ const PreconditionBase &Lp_preconditioner,
+ const SparseMatrix &Fp_matrix)
+ {
+ AssertThrow (system_matrix.n_block_rows() == 2,
+ ExcDimensionMismatch(system_matrix.n_block_rows(),
+ 2));
+ AssertThrow (system_matrix.n_block_cols() == 2,
+ ExcDimensionMismatch(system_matrix.n_block_cols(),
+ 2));
+
+ // Create Av inverse matrix
+ Thyra::ConstLinearOperator<double> Av_inverse;
+ if (additional_data.Av_do_solve)
+ Av_inverse = inverse_matrix(system_matrix.block(0,0),
+ &*Av_preconditioner.preconditioner,
+ additional_data.Av_is_symmetric,
+ additional_data.Av_n_iters,
+ additional_data.Av_tolerance,
+ additional_data.output_details);
+ else
+ {
+ Teuchos::RCP<const Thyra::LinearOpBase<double> >
+ tmpAvp = Thyra::epetraLinearOp(Teuchos::rcp(&*Av_preconditioner.preconditioner,
+ false),
+ Thyra::NOTRANS,
+ Thyra::EPETRA_OP_APPLY_APPLY_INVERSE,
+ Thyra::EPETRA_OP_ADJOINT_SUPPORTED);
+ Av_inverse = tmpAvp;
+ }
+
+
+ // Create Mp inverse matrix.
+ Thyra::ConstLinearOperator<double> Mp_inverse
+ = inverse_matrix(Mp_matrix, &*Mp_preconditioner.preconditioner, true,
+ Mp_matrix.m(),
+ additional_data.inner_solve_tolerance,
+ additional_data.output_details);
+
+ // Create Lp preconditioner
+ // and inverse matrix in case
+ // there is a matrix given.
+ Thyra::ConstLinearOperator<double> Lp_inverse;
+
+ if (Lp_matrix.m() == system_matrix.block(1,1).m())
+ {
+ // Create Lp inverse matrix.
+ Lp_inverse = inverse_matrix(Lp_matrix, &*Lp_preconditioner.preconditioner,
+ true, Lp_matrix.m(),
+ additional_data.inner_solve_tolerance,
+ additional_data.output_details);
+ }
+ else
+ Lp_inverse = Thyra::identity<double> (Mp_inverse.range());
+
+ // Create Fp Thyra operator.
+ Thyra::ConstLinearOperator<double> Fp_op;
+ if (Fp_matrix.m() != system_matrix.block(1,1).m())
+ {
+ Fp_op = Thyra::identity<double> (Mp_inverse.range());
+ }
+ else
+ {
+ Teuchos::RCP<const Thyra::LinearOpBase<double> >
+ tmpFp = Thyra::epetraLinearOp(Teuchos::rcp(&*Fp_matrix.matrix,false));
+ Fp_op = tmpFp;
+ }
+
+ // Build the PCD block
+ // preconditioner factory.
+ // Build identity matrices on
+ // the velocity and pressure
+ // spaces
+ Thyra::ConstLinearOperator<double> Ivel
+ = Thyra::identity<double>(Av_inverse.range());
+ Thyra::ConstLinearOperator<double> Ipress
+ = Thyra::identity<double>(Mp_inverse.domain());
+
+ // Build zero operators. Need
+ // one that is pressure x
+ // velocity and one that is
+ // velocity x pressure
+ Thyra::ConstLinearOperator<double> zero;
+
+ // Build the composed Schur
+ // complement approximation
+ // inverse
+ Thyra::ConstLinearOperator<double> Xinv;
+ if (additional_data.right_preconditioning == true)
+ Xinv = Mp_inverse * Fp_op * Lp_inverse;
+ else
+ Xinv = Lp_inverse * Fp_op * Mp_inverse;
+
+ // Create discrete pressure
+ // gradient and divergence
+ // operators in Thyra format.
+ Teuchos::RCP<const Thyra::LinearOpBase<double> > tmpS01 =
+ Thyra::epetraLinearOp(Teuchos::rcp(&*system_matrix.block(0,1).matrix,
+ false));
+ Thyra::ConstLinearOperator<double> S01 = tmpS01;
+
+ Teuchos::RCP<const Thyra::LinearOpBase<double> > tmpS10 =
+ Thyra::epetraLinearOp(Teuchos::rcp(&*system_matrix.block(1,0).matrix,
+ false));
+ Thyra::ConstLinearOperator<double> S10 = tmpS10;
+
+ // Build the 3 block operators
+ // for the
+ // preconditioner. Here we
+ // have to set different
+ // matrices depending on
+ // whether we do right or left
+ // preconditioning.
+ Thyra::ConstLinearOperator<double> P1
+ = block2x2( Av_inverse, zero, zero, Ipress );
+
+ Thyra::ConstLinearOperator<double> P2;
+ if (additional_data.right_preconditioning == true)
+ P2 = block2x2( Ivel, (-1.0)*S01, zero, Ipress );
+ else
+ P2 = block2x2( Ivel, zero, (-1.0)*S10, Ipress );
+
+ Thyra::ConstLinearOperator<double> P3
+ = block2x2( Ivel, zero, zero, (-1.0)*Xinv );
+
+ if (additional_data.right_preconditioning == true)
+ preconditioner = Teuchos::rcp (new
+ Thyra::ConstLinearOperator<double>(P1 * P2 * P3));
+ else
+ preconditioner = Teuchos::rcp (new
+ Thyra::ConstLinearOperator<double>(P3 * P2 * P1));
+
+ }
+
+
+
+ void
+ PreconditionStokes::reinit_lazy ()
+ {
+ Assert (&*Av_precondition != 0,
+ ExcMessage ("Can only be done when the outer preconditioner owns"
+ " the AMG object!"));
+
+ Av_precondition->reinit();
+ }
+
+} /* end of namespace TrilinosWrappers */
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif // DEAL_II_USE_TRILINOS
SolverBase::SolverBase (SolverControl &cn)
:
- solver_control (cn),
- solver_name (gmres)
+ solver_name (gmres),
+ solver_control (cn)
+ {}
+
+
+
+ SolverBase::SolverBase (const enum SolverBase::SolverName solver_name,
+ SolverControl &cn)
+ :
+ solver_name (solver_name),
+ solver_control (cn)
{}
--- /dev/null
+//---------------------------------------------------------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2008 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 <lac/trilinos_solver_block.h>
+#include <lac/trilinos_block_sparse_matrix.h>
+#include <lac/trilinos_block_vector.h>
+#include <lac/trilinos_precondition_block.h>
+
+#include <base/utilities.h>
+
+#include <cmath>
+
+#ifdef DEAL_II_USE_TRILINOS
+
+#include <Teuchos_RCP.hpp>
+#include <Teuchos_ParameterList.hpp>
+
+#include <Epetra_Vector.h>
+
+#include <Thyra_VectorDecl.hpp>
+#include <Thyra_VectorImpl.hpp>
+#include <Thyra_VectorSpaceImpl.hpp>
+#include <Thyra_MultiVectorBase.hpp>
+#include <Thyra_MultiVectorDefaultBase.hpp>
+#include <Thyra_LinearOperatorDecl.hpp>
+#include <Thyra_LinearOperatorImpl.hpp>
+#include <Thyra_DefaultBlockedLinearOpDecl.hpp>
+#include <Thyra_DefaultBlockedLinearOp.hpp>
+#include <Thyra_EpetraLinearOp.hpp>
+#include <Thyra_EpetraThyraWrappers.hpp>
+#include <Thyra_LinearOpWithSolveFactoryHelpers.hpp>
+#include <Thyra_AztecOOLinearOpWithSolveFactory.hpp>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace TrilinosWrappers
+{
+
+ SolverBlockBase::AdditionalData::
+ AdditionalData (const unsigned int gmres_restart_parameter,
+ const bool right_preconditioning,
+ const bool output_details)
+ :
+ gmres_restart_parameter (gmres_restart_parameter),
+ right_preconditioning (right_preconditioning),
+ output_details (output_details)
+ {}
+
+
+
+ SolverBlockBase::SolverBlockBase (SolverControl &cn)
+ :
+ solver_control (cn),
+ solver_name (gmres)
+ {}
+
+
+
+ SolverBlockBase::~SolverBlockBase ()
+ {}
+
+
+
+ SolverControl &
+ SolverBlockBase::control() const
+ {
+ return solver_control;
+ }
+
+
+
+ void
+ SolverBlockBase::solve (const BlockSparseMatrix &A,
+ BlockVector &x,
+ const BlockVector &b,
+ const PreconditionBlockBase &preconditioner)
+ {
+
+ Assert (x.n_blocks() == b.n_blocks(),
+ ExcDimensionMismatch(x.n_blocks(), b.n_blocks()));
+
+ // Just copy everything over
+ // to a distributed vector...
+ MPI::BlockVector tmpx (x.n_blocks());
+ MPI::BlockVector tmpb (b.n_blocks());
+ for (unsigned int i=0; i<x.n_blocks(); ++i)
+ {
+ tmpx.block(i).reinit (A.block(i,i).matrix->DomainMap(), false);
+ tmpx.block(i) = x.block(i);
+ tmpb.block(i).reinit (A.block(i,i).matrix->RangeMap(), false);
+ tmpb.block(i) = b.block(i);
+ }
+ tmpx.collect_sizes();
+ tmpb.collect_sizes();
+
+ solve (A, tmpx, tmpb, preconditioner);
+ x = tmpx;
+ }
+
+
+
+ void
+ SolverBlockBase::solve (const BlockSparseMatrix &input_A,
+ MPI::BlockVector &input_x,
+ const MPI::BlockVector &input_b,
+ const PreconditionBlockBase &preconditioner)
+ {
+
+ // Get block information about
+ // the preconditioner.
+ Thyra::ConstLinearOperator<double> P = *preconditioner.preconditioner;
+ const unsigned int n_rows = P.numBlockRows();
+ const unsigned int n_cols = P.numBlockCols();
+
+ // Check the dimensions of the
+ // block vectors and matrix
+ // and check if the parallel
+ // distribution of the
+ // matrices and vectors is the
+ // same.
+ AssertThrow (n_rows == n_cols,
+ ExcDimensionMismatch (n_rows, n_cols));
+
+ AssertThrow (input_x.n_blocks() == input_A.n_block_cols(),
+ ExcDimensionMismatch(input_x.n_blocks(),
+ input_A.n_block_cols()));
+ AssertThrow (input_b.n_blocks() == input_A.n_block_rows(),
+ ExcDimensionMismatch(input_b.n_blocks(),
+ input_A.n_block_rows()));
+ AssertThrow (input_A.n_block_rows() == n_rows,
+ ExcDimensionMismatch(input_A.n_block_rows(), n_rows));
+ AssertThrow (input_A.n_block_cols() == n_cols,
+ ExcDimensionMismatch(input_A.n_block_cols(), n_cols));
+
+ for (unsigned int i=0; i<n_rows; ++i)
+ {
+ AssertThrow (input_x.block(i).vector->Map().UniqueGIDs() == true,
+ ExcOverlappingMaps("vector", "x"));
+ AssertThrow (input_b.block(i).vector->Map().UniqueGIDs() == true,
+ ExcOverlappingMaps("vector", "b"));
+ for (unsigned int j=0; j<n_cols; ++j)
+ {
+ {
+ std::ostringstream i_str;
+ i_str << i;
+ std::ostringstream error_component;
+ error_component << "x.block(" << i_str.str() << ")";
+ AssertThrow (input_x.block(j).vector->Map().SameAs(
+ input_A.block(i,j).matrix->DomainMap()) == true,
+ ExcNonMatchingMaps (error_component.str()));
+ }
+ {
+ std::ostringstream i_str;
+ i_str << j;
+ std::ostringstream error_component;
+ error_component << "b.block(" << i_str.str() << ")";
+ AssertThrow (input_b.block(i).vector->Map().SameAs(
+ input_A.block(i,j).matrix->RangeMap()) == true,
+ ExcNonMatchingMaps (error_component.str()));
+ }
+ }
+ }
+
+
+ // Wrap Epetra vectors into
+ // Thyra vectors in order to
+ // be able to work with the
+ // abstract Thyra/Stratimikos
+ // interfaces to the AztecOO
+ // solver.
+ //
+ // As a first step, create a
+ // copy of the vector spaces.
+ std::vector<Teuchos::RCP<const Thyra::VectorSpaceBase<double> > >
+ epetra_vector_spaces;
+ for (unsigned int i=0; i<n_rows; ++i)
+ {
+ Teuchos::RCP<const Thyra::VectorSpaceBase<double> > tmp_space
+ = Thyra::create_VectorSpace(
+ Teuchos::rcp(&input_A.block(i,i).matrix->DomainMap(),
+ false));
+
+ epetra_vector_spaces.push_back(tmp_space);
+ }
+
+ // Convert the block matrix to
+ // a Thyra linear operator
+ // that acts on blocks.
+ Teuchos::RCP<Thyra::PhysicallyBlockedLinearOpBase<double> >
+ tmpA = Teuchos::rcp(new Thyra::DefaultBlockedLinearOp<double>());
+ tmpA->beginBlockFill(n_rows,n_cols);
+ for (unsigned int i=0; i<n_rows; ++i)
+ for (unsigned int j=0; j<n_cols; ++j)
+ {
+ Teuchos::RCP<const Thyra::LinearOpBase<double> >
+ A_ij = Thyra::epetraLinearOp(Teuchos::rcp(&*input_A.block(i,j).matrix,
+ false));
+ tmpA->setBlock(i, j, A_ij);
+ }
+ tmpA->endBlockFill();
+ tmpA->setObjectLabel("Blocked system matrix");
+ Thyra::RCP<const Thyra::LinearOpBase<double> > baseA = tmpA;
+
+ Thyra::ConstLinearOperator<double> A = baseA;
+
+ // Convert block vectors for
+ // right hand side and
+ // solution to Thyra objects.
+ Thyra::Vector<double> rhs = A.range().createMember();
+ Thyra::Vector<double> sol = A.domain().createMember();
+ for (unsigned int i=0; i<n_rows; ++i)
+ {
+ Epetra_Vector *block_ptr = (*input_b.block(i).vector)(0);
+ Teuchos::RCP<Thyra::VectorBase<double> > tmp_rhs_i
+ = Thyra::create_Vector(Teuchos::rcp(block_ptr, false),
+ epetra_vector_spaces[i]);
+ const Thyra::Vector<double> rhs_i = tmp_rhs_i;
+ rhs.setBlock(i, rhs_i);
+
+ block_ptr = (*input_x.block(i).vector)(0);
+ Teuchos::RCP<Thyra::VectorBase<double> > tmp_sol_i
+ = Thyra::create_Vector(Teuchos::rcp(block_ptr, false),
+ epetra_vector_spaces[i]);
+ Thyra::Vector<double> sol_i = tmp_sol_i;
+ sol.setBlock(i, sol_i);
+ }
+
+ // ------- Now build a solver factory for the block problem -------
+
+ // Set up a parameter list for
+ // the AztecOO solver.
+ Teuchos::RCP<Teuchos::ParameterList> aztecBlockParams
+ = rcp(new Teuchos::ParameterList("AztecOO block solver parameters"));
+
+ // Set the number of
+ // iterations and the solver
+ // tolerance from the solver
+ // control object.
+ aztecBlockParams->sublist("Forward Solve").set("Max Iterations",
+ (int)solver_control.max_steps());
+ aztecBlockParams->sublist("Forward Solve").set("Tolerance",
+ solver_control.tolerance());
+ aztecBlockParams->sublist("Forward Solve")
+ .sublist("AztecOO Settings").set("Convergence Test", "no scaling");
+
+ // The settings about the
+ // actual solver that is going
+ // to be used - let this be
+ // determined by the
+ // <tt>solver_name</tt> set in
+ // a derived class.
+ switch (solver_name)
+ {
+ case cg:
+ aztecBlockParams->sublist("Forward Solve")
+ .sublist("AztecOO Settings").set("Aztec Solver", "CG");
+ break;
+ case gmres:
+ aztecBlockParams->sublist("Forward Solve")
+ .sublist("AztecOO Settings").set("Aztec Solver", "GMRES");
+ aztecBlockParams->sublist("Forward Solve")
+ .sublist("AztecOO Settings").set("Size of Krylov Subspace",
+ (int)additional_data.gmres_restart_parameter);
+ break;
+ default:
+ ExcNotImplemented();
+ }
+
+ // We use an externally
+ // provided preconditioner, so
+ // do not set anything here.
+ aztecBlockParams->sublist("Forward Solve")
+ .sublist("AztecOO Settings").set("Aztec Preconditioner", "none");
+
+ if (additional_data.output_details == true)
+ aztecBlockParams->sublist("Forward Solve")
+ .sublist("AztecOO Settings").set("Output Frequency", 1);
+ else
+ aztecBlockParams->sublist("VerboseObject").set("Verbosity Level",
+ "none");
+
+ // Create solve factory from
+ // parameter list.
+ Teuchos::RCP<Thyra::LinearOpWithSolveFactoryBase<double> >
+ aztecBlockLowsFactory = Teuchos::rcp
+ (new Thyra::AztecOOLinearOpWithSolveFactory());
+
+ aztecBlockLowsFactory->setParameterList(aztecBlockParams);
+
+ Teuchos::RCP<Thyra::LinearOpWithSolveBase<double> > rcpBlockAztec
+ = aztecBlockLowsFactory->createOp();
+
+
+ // ---- Set up the preconditioned inverse object and do the solve ----
+
+ // Extract base operators from
+ // block matrix A and block
+ // preconditioner P.
+ Teuchos::RCP<const Thyra::LinearOpSourceBase<double> > A_ptr
+ = Teuchos::rcp(new Thyra::DefaultLinearOpSource<double>(A.constPtr()),
+ true);
+ Teuchos::RCP<Thyra::PreconditionerBase<double> > P_ptr =
+ Teuchos::rcp (new Thyra::DefaultPreconditioner<double>(), true);
+ {
+ Thyra::DefaultPreconditioner<double>
+ *defaultPrec = &Teuchos::dyn_cast<Thyra::DefaultPreconditioner<double> >
+ (*P_ptr);
+ if (additional_data.right_preconditioning == true)
+ (*defaultPrec).initializeRight(P.constPtr());
+ else
+ (*defaultPrec).initializeUnspecified(P.constPtr());
+ }
+
+ // Attach the matrix and the
+ // preconditioner to the
+ // solver object.
+ aztecBlockLowsFactory->initializePreconditionedOp(A_ptr, P_ptr,
+ &*rcpBlockAztec,
+ Thyra::SUPPORT_SOLVE_FORWARD_ONLY);
+
+ // Now do the solve.
+ Thyra::SolveStatus<double>
+ solveStatus = Thyra::solve( *rcpBlockAztec, Thyra::NOTRANS,
+ *rhs.constPtr().get(), &*sol.ptr().get() );
+
+ // Extract the number of
+ // needed iterations from the
+ // solver message.
+ unsigned int n_iterations;
+ {
+ std::string output_text = solveStatus.message;
+ unsigned int position = output_text.find ("using ");
+ const std::pair<int,unsigned int> tmp
+ = Utilities::get_integer_at_position (output_text, position+6);
+ n_iterations = tmp.first;
+ }
+
+ // Do a dummy check of
+ // convergence to tell the
+ // SolverControl object the
+ // number of iterations. We do
+ // not impose the achieved
+ // tolerance here since that
+ // would require additional
+ // work (computing a
+ // residual).
+ solver_control.check (n_iterations, 0);
+
+ }
+
+
+
+
+
+/* ------------------------ SolverCG -------------------------- */
+
+ SolverBlockCG::AdditionalData::
+ AdditionalData (const bool right_preconditioning,
+ const bool output_details)
+ :
+ right_preconditioning (right_preconditioning),
+ output_details (output_details)
+ {}
+
+
+
+ SolverBlockCG::SolverBlockCG (SolverControl &cn,
+ const AdditionalData &data)
+ :
+ SolverBlockBase (cn),
+ data (data)
+ {
+ solver_name = cg;
+ additional_data.right_preconditioning = data.right_preconditioning;
+ additional_data.output_details = data.output_details;
+ }
+
+
+/* ---------------------- SolverGMRES ------------------------- */
+
+ SolverBlockGMRES::AdditionalData::
+ AdditionalData (const unsigned int gmres_restart_parameter,
+ const bool right_preconditioning,
+ const bool output_details)
+ :
+ gmres_restart_parameter (gmres_restart_parameter),
+ right_preconditioning (right_preconditioning),
+ output_details (output_details)
+ {}
+
+
+
+ SolverBlockGMRES::SolverBlockGMRES (SolverControl &cn,
+ const AdditionalData &data)
+ :
+ SolverBlockBase (cn),
+ data (data)
+ {
+ solver_name = gmres;
+ additional_data.gmres_restart_parameter = data.gmres_restart_parameter;
+ additional_data.right_preconditioning = data.right_preconditioning;
+ additional_data.output_details = data.output_details;
+ }
+
+
+}
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif // DEAL_II_USE_PETSC