*/
void reinit (const unsigned int num_blocks);
- /**
- * This reinit function is
- * meant to be used for
- * parallel calculations where
- * some non-local data has to
- * be used. The typical
- * situation where one needs
- * this function is the call of
- * the
- * FEValues<dim>::get_function_values
- * function (or of some
- * derivatives) in
- * parallel. Since it is
- * usually faster to retrieve
- * the data in advance, this
- * function can be called
- * before the assembly forks
- * out to the different
- * processors. What this
- * function does is the
- * following: It takes the
- * information in the columns
- * of the given matrix and
- * looks which data couples
- * between the different
- * processors. That data is
- * then queried from the input
- * vector. Note that you should
- * not write to the resulting
- * vector any more, since the
- * some data can be stored
- * several times on different
- * processors, leading to
- * unpredictable results. In
- * particular, such a vector
- * cannot be used for
- * matrix-vector products as
- * for example done during the
- * solution of linear systems.
- */
+ /**
+ * This reinit function is meant to
+ * be used for parallel
+ * calculations where some
+ * non-local data has to be
+ * used. The typical situation
+ * where one needs this function is
+ * the call of the
+ * FEValues<dim>::get_function_values
+ * function (or of some
+ * derivatives) in parallel. Since
+ * it is usually faster to retrieve
+ * the data in advance, this
+ * function can be called before
+ * the assembly forks out to the
+ * different processors. What this
+ * function does is the following:
+ * It takes the information in the
+ * columns of the given matrix and
+ * looks which data couples between
+ * the different processors. That
+ * data is then queried from the
+ * input vector. Note that you
+ * should not write to the
+ * resulting vector any more, since
+ * the some data can be stored
+ * several times on different
+ * processors, leading to
+ * unpredictable results. In
+ * particular, such a vector cannot
+ * be used for matrix-vector
+ * products as for example done
+ * during the solution of linear
+ * systems.
+ */
void import_nonlocal_data_for_fe (const TrilinosWrappers::BlockSparseMatrix &m,
const BlockVector &v);
+
+
+
+/**
+ * A wrapper class for a thresholded incomplete LU factorization (ILU-T)
+ * preconditioner for Trilinos matrices. This preconditioner works both in
+ * serial and in parallel, depending on the matrix it is based on. In
+ * general, an incomplete factorization does not take all fill-in elements
+ * that would appear in a full factorization (that is the basis for a direct
+ * solve). For the ILU-T precondtioner, the parameter <tt>ilut_drop</tt>
+ * lets the user specify which elements should be dropped (i.e., should not
+ * be part of the incomplete decomposition). Trilinos calculates first the
+ * complete factorization for one row, and then skips those elements that
+ * are lower than the threshold. This is the main difference to the
+ * non-thresholded ILU preconditioner, where the parameter
+ * <tt>ilut_fill</tt> governs the incomplete factorization structure. This
+ * parameter is available here as well, but provides only some extra
+ * information here.
+ *
+ * The AdditionalData data structure allows to set preconditioner
+ * options. Besides the fill-in arguments, these options are some options
+ * for perturbations (see the documentation of the AdditionalData structure
+ * for details), and a parameter <tt>overlap</tt> that determines if and how
+ * much overlap there should be between the matrix partitions on the various
+ * MPI processes. The default settings are 0 for the additional fill-in, 0
+ * for the absolute augmentation tolerance, 1 for the relative augmentation
+ * tolerance, 0 for the overlap.
+ *
+ * Note that a parallel applicatoin of the ILU-T preconditioner is
+ * actually a block-Jacobi preconditioner with block size equal to the
+ * local matrix size. Spoken more technically, this parallel operation
+ * is an <a
+ * href="http://en.wikipedia.org/wiki/Additive_Schwarz_method">additive
+ * Schwarz method</a> with an ILU <em>approximate solve</em> as inner
+ * solver, based on the (outer) parallel partitioning.
+ *
+ * @ingroup TrilinosWrappers
+ * @ingroup Preconditioners
+ * @author Martin Kronbichler, 2009
+ */
+ class PreconditionILUT : public PreconditionBase
+ {
+ public:
+ /**
+ * Standardized data struct to pipe
+ * additional parameters to the
+ * preconditioner. The Trilinos ILU-T
+ * decomposition allows for some
+ * fill-in, so it actually is a
+ * threshold incomplete LU
+ * factorization. The amount of
+ * fill-in, and hence, the amount of
+ * memory used by this
+ * preconditioner, is controlled by
+ * the parameters <tt>ilut_drop</tt>
+ * and <tt>ilut_fill</tt>, which
+ * specifies a threshold about which
+ * values should form the incomplete
+ * factorization and the level of
+ * additional fill-in. When forming
+ * the preconditioner, for certain
+ * problems bad conditioning (or just
+ * bad luck) can cause the
+ * preconditioner to be very poorly
+ * conditioned. Hence it can help to
+ * add diagonal perturbations to the
+ * original matrix and form the
+ * preconditioner for this slightly
+ * better matrix. <tt>ilut_atol</tt>
+ * is an absolute perturbation that
+ * is added to the diagonal before
+ * forming the prec, and
+ * <tt>ilu_rtol</tt> is a scaling
+ * factor $rtol \geq 1$. The last
+ * parameter specifies the overlap of
+ * the partitions when the
+ * preconditioner runs in parallel.
+ */
+ struct AdditionalData
+ {
+ /**
+ * Constructor. By default, no
+ * element will be dropped, the level
+ * of extra fill-ins is set to be
+ * zero (just use the matrix
+ * structure, do not generate any
+ * additional fill-in except the one
+ * that results from non-dropping
+ * large elements), the tolerance
+ * level are 0 and 1, respectively,
+ * and the overlap in case of a
+ * parallel execution is zero. This
+ * overlap in a block-application of
+ * the ILU in the parallel case makes
+ * the preconditioner a so-called
+ * additive Schwarz preconditioner.
+ */
+ AdditionalData (const double ilut_drop = 0.,
+ const unsigned int ilut_fill = 0,
+ const double ilut_atol = 0.,
+ const double ilut_rtol = 1.,
+ const unsigned int overlap = 0);
+
+ /**
+ * This specifies the relative size
+ * of elements which should be
+ * dropped when forming an incomplete
+ * LU decomposition with threshold.
+ */
+ double ilut_drop;
+
+ /**
+ * This specifies the amount of
+ * additional fill-in elements
+ * besides the sparse matrix
+ * structure. When
+ * <tt>ilu_fill</tt> is large,
+ * this means that many
+ * fill-ins will be added, so
+ * that the ILU preconditioner
+ * comes closer to a (direct)
+ * sparse LU
+ * decomposition. Note,
+ * however, that this will
+ * drastically increase the
+ * memory requirement,
+ * especially when the
+ * preconditioner is used in
+ * 3D.
+ */
+ unsigned int ilut_fill;
+
+ /**
+ * This specifies the amount of
+ * an absolute perturbation
+ * that will be added to the
+ * diagonal of the matrix,
+ * which sometimes can help to
+ * get better preconditioners.
+ */
+ double ilut_atol;
+
+ /**
+ * This specifies the factor by
+ * which the diagonal of the
+ * matrix will be scaled, which
+ * sometimes can help to get
+ * better preconditioners.
+ */
+ double ilut_rtol;
+
+ /**
+ * This determines how large
+ * the overlap of the local
+ * matrix portions on each
+ * processor in a parallel
+ * application should be.
+ */
+ unsigned int overlap;
+ };
+
+ /**
+ * Initialize function. Takes
+ * the matrix which is used to
+ * form the preconditioner, and
+ * additional flags if there
+ * are any.
+ */
+ void initialize (const SparseMatrix &matrix,
+ const AdditionalData &additional_data = AdditionalData());
+
+ private:
+ /**
+ * This is a pointer to the
+ * Ifpack data contained in
+ * this preconditioner.
+ */
+ Teuchos::RCP<Ifpack_Preconditioner> ifpack;
+ };
+
+
+
/**
* A wrapper class for a sparse direct LU decomposition on parallel
* blocks for Trilinos matrices. When run in serial, this corresponds
#include <lac/solver_control.h>
#include <boost/shared_ptr.hpp>
+#include <lac/vector.h>
#ifdef DEAL_II_USE_TRILINOS
#include <Epetra_LinearProblem.h>
#include <AztecOO.h>
#include <Epetra_Operator.h>
+#include <Amesos.h>
DEAL_II_NAMESPACE_OPEN
-
namespace TrilinosWrappers
{
// forward declarations
* deal.II class SolverSelector.
*
* @ingroup TrilinosWrappers
- * @author Martin Kronbichler, 2008
+ * @author Martin Kronbichler, 2008, 2009
*/
class SolverBase
{
* pipe additional data to the
* solver.
*/
+
struct AdditionalData
{
/**
- * Restart parameter in case
- * the derived class is
- * GMRES.
+ * Sets the additional data field to
+ * the desired output format and puts
+ * the restart parameter in case the
+ * derived class is GMRES.
*
* TODO: Find a better way for
- * doing this.
+ * setting the GMRES restart
+ * parameter since it is quite
+ * inelegant to set a specific option
+ * of one solver in the base class
+ * for all solvers.
+ */
+ AdditionalData (const bool output_solver_details = false,
+ const unsigned int gmres_restart_parameter = 30);
+
+ /**
+ * Enables/disables the output of
+ * solver details (residual in each
+ * iterations etc.).
*/
- AdditionalData (const unsigned int gmres_restart_parameter = 30);
+ const bool output_solver_details;
+ /**
+ * Restart parameter for GMRES
+ * solver.
+ */
const unsigned int gmres_restart_parameter;
};
const VectorBase &b,
const PreconditionBase &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 class
+ * works with matrices according to
+ * the TrilinosWrappers format, but
+ * can take deal.II vectors as
+ * argument. Since deal.II are serial
+ * vectors (not distributed), this
+ * function does only what you expect
+ * in case the matrix is locally
+ * owned. Otherwise, an exception
+ * will be thrown.
+ */
+ void
+ solve (const SparseMatrix &A,
+ dealii::Vector<double> &x,
+ const dealii::Vector<double> &b,
+ const PreconditionBase &preconditioner);
+
/**
* Access to object that controls
* convergence.
* pipe additional data to the
* solver.
*/
+
struct AdditionalData
- {};
-
+ {
+ /**
+ * Sets the additional data field to
+ * the desired output format.
+ */
+ AdditionalData (const bool output_solver_details = false);
+
+ /**
+ * Enables/disables the output of
+ * solver details (residual in each
+ * iterations etc.).
+ */
+ bool output_solver_details;
+ };
+
/**
* Constructor. In contrast to
* deal.II's own solvers, there is no
* solver.
*/
struct AdditionalData
- {};
+ {
+ /**
+ * Sets the additional data field to
+ * the desired output format.
+ */
+ AdditionalData (const bool output_solver_details = false);
+
+ /**
+ * Enables/disables the output of
+ * solver details (residual in each
+ * iterations etc.).
+ */
+ bool output_solver_details;
+ };
/**
* Constructor. In contrast to
* 30, i.e. do a restart every 30
* iterations.
*/
- AdditionalData (const unsigned int restart_parameter = 30);
-
- /**
- * Maximum number of
- * tmp vectors.
- */
+ AdditionalData (const bool output_solver_details = false,
+ const unsigned int restart_parameter = 30);
+
+ /**
+ * Enables/disables the output of
+ * solver details (residual in each
+ * iterations etc.).
+ */
+ bool output_solver_details;
+
+ /**
+ * Maximum number of
+ * tmp vectors.
+ */
unsigned int restart_parameter;
};
* solver.
*/
struct AdditionalData
- {};
+ {
+ /**
+ * Sets the additional data field to
+ * the desired output format.
+ */
+ AdditionalData (const bool output_solver_details = false);
+
+ /**
+ * Enables/disables the output of
+ * solver details (residual in each
+ * iterations etc.).
+ */
+ bool output_solver_details;
+ };
/**
* Constructor. In contrast to
* solver.
*/
struct AdditionalData
- {};
+ {
+ /**
+ * Sets the additional data field to
+ * the desired output format.
+ */
+ AdditionalData (const bool output_solver_details = false);
+
+ /**
+ * Enables/disables the output of
+ * solver details (residual in each
+ * iterations etc.).
+ */
+ bool output_solver_details;
+ };
/**
* Constructor. In contrast to
const AdditionalData additional_data;
};
+
+
+/**
+ * An implementation of the Trilinos KLU direct solver (using the Amesos
+ * package).
+ *
+ * @ingroup TrilinosWrappers
+ * @author Martin Kronbichler, 2009
+ */
+ class SolverDirect
+ {
+ public:
+
+ /**
+ * Standardized data struct to
+ * pipe additional data to the
+ * solver.
+ */
+
+ struct AdditionalData
+ {
+ /**
+ * Sets the additional data field to
+ * the desired output format.
+ */
+ AdditionalData (const bool output_solver_details = false);
+
+ /**
+ * Enables/disables the output of
+ * solver details (residual in each
+ * iterations etc.).
+ */
+ bool output_solver_details;
+ };
+
+ /**
+ * Constructor. Takes the
+ * solver control object and
+ * creates the solver.
+ */
+ SolverDirect (SolverControl &cn,
+ const AdditionalData &data = AdditionalData());
+
+ /**
+ * Destructor.
+ */
+ virtual ~SolverDirect ();
+
+ /**
+ * Solve the linear system
+ * <tt>Ax=b</tt>. Creates a KLU
+ * factorization of the matrix and
+ * performs the solve. Note that
+ * there is no need for a
+ * preconditioner here.
+ */
+ void
+ solve (const SparseMatrix &A,
+ VectorBase &x,
+ const VectorBase &b);
+
+ /**
+ * 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 class
+ * works with matrices according to
+ * the TrilinosWrappers format, but
+ * can take deal.II vectors as
+ * argument. Since deal.II are serial
+ * vectors (not distributed), this
+ * function does only what you expect
+ * in case the matrix is locally
+ * owned. Otherwise, an exception
+ * will be thrown.
+ */
+ void
+ solve (const SparseMatrix &A,
+ dealii::Vector<double> &x,
+ const dealii::Vector<double> &b);
+
+ /**
+ * 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");
+
+ private:
+
+ /**
+ * 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;
+
+ /**
+ * A structure that collects
+ * the Trilinos sparse matrix,
+ * the right hand side vector
+ * and the solution vector,
+ * which is passed down to the
+ * Trilinos solver.
+ */
+ std::auto_ptr<Epetra_LinearProblem> linear_problem;
+
+ /**
+ * A structure that contains
+ * the Trilinos solver and
+ * preconditioner objects.
+ */
+ std::auto_ptr<Amesos_BaseSolver> solver;
+
+ /**
+ * Store a copy of the flags for this
+ * particular solver.
+ */
+ const AdditionalData additional_data;
+
+ };
+
}
DEAL_II_NAMESPACE_CLOSE
Assert (values.m() == values.n(), ExcNotQuadratic());
for (unsigned int i=0; i<indices.size(); ++i)
- add (indices[i], indices.size(), &indices[0], &values(i,0),
+ add (indices[i], indices.size(), &indices[0], &values(i,0),
elide_zero_values);
}
ExcDimensionMismatch(col_indices.size(), values.n()));
for (unsigned int i=0; i<row_indices.size(); ++i)
- add (row_indices[i], col_indices.size(), &col_indices[0], &values(i,0),
- elide_zero_values);
+ add (row_indices[i], col_indices.size(), &col_indices[0],
+ &values(i,0), elide_zero_values);
}
{
// forward declaration
class VectorBase;
+ class SparseMatrix;
/**
void clear ();
/**
- * Reinit functionality, sets
- * the dimension and the
- * possibly parallel
- * partitioning (Epetra_Map) of
- * the calling vector the value
- * of the input vector.
+ * Reinit functionality, sets the
+ * dimension and possibly the
+ * parallel partitioning (Epetra_Map)
+ * of the calling vector to the
+ * settings of the input vector.
*/
void reinit (const VectorBase &v,
const bool fast = false);
this->components.resize(this->n_blocks());
for (unsigned int i=0;i<this->n_blocks();++i)
- block(i).clear();
+ components[i].clear();
collect_sizes();
}
+/* -------------------------- PreconditionILUT -------------------------- */
+
+ PreconditionILUT::AdditionalData::
+ AdditionalData (const double ilut_drop,
+ const unsigned int ilut_fill,
+ const double ilut_atol,
+ const double ilut_rtol,
+ const unsigned int overlap)
+ :
+ ilut_drop (ilut_drop),
+ ilut_fill (ilut_fill),
+ ilut_atol (ilut_atol),
+ ilut_rtol (ilut_rtol),
+ overlap (overlap)
+ {}
+
+
+
+ void
+ PreconditionILUT::initialize (const SparseMatrix &matrix,
+ const AdditionalData &additional_data)
+ {
+ preconditioner.release();
+ ifpack.release();
+
+ ifpack = Teuchos::rcp (Ifpack().Create ("ILUT", &*matrix.matrix,
+ additional_data.overlap));
+
+ Assert (&*ifpack != 0, ExcMessage ("Trilinos could not create this "
+ "preconditioner"));
+
+ int ierr;
+
+ Teuchos::ParameterList parameter_list;
+ parameter_list.set ("fact: drop value",additional_data.ilut_drop);
+ parameter_list.set ("fact: level-of-fill",(int)additional_data.ilut_fill);
+ parameter_list.set ("fact: absolute threshold",additional_data.ilut_atol);
+ parameter_list.set ("fact: relative threshold",additional_data.ilut_rtol);
+ parameter_list.set ("schwarz: combine mode", "Add");
+
+ ierr = ifpack->SetParameters(parameter_list);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+ ierr = ifpack->Initialize();
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+ ierr = ifpack->Compute();
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+ preconditioner = Teuchos::rcp (&*ifpack, false);
+ }
+
+
+
/* ---------------------- PreconditionBlockDirect --------------------- */
PreconditionBlockwiseDirect::AdditionalData::
//---------------------------------------------------------------------------
+#include <base/conditional_ostream.h>
#include <lac/trilinos_solver.h>
#include <lac/trilinos_sparse_matrix.h>
#include <lac/trilinos_vector_base.h>
#ifdef DEAL_II_USE_TRILINOS
-#include <Epetra_Operator.h>
-
DEAL_II_NAMESPACE_OPEN
namespace TrilinosWrappers
{
- SolverBase::AdditionalData::AdditionalData (const unsigned int gmres_restart_parameter)
+ SolverBase::AdditionalData::AdditionalData (const bool output_solver_details,
+ const unsigned int gmres_restart_parameter)
:
+ output_solver_details (output_solver_details),
gmres_restart_parameter (gmres_restart_parameter)
{}
-
-
+
+
SolverBase::SolverBase (SolverControl &cn)
:
solver_name (gmres),
solver_control (cn)
{}
-
+
SolverBase::SolverBase (const enum SolverBase::SolverName solver_name,
SolverControl &cn)
solver_control (cn)
{}
-
+
SolverBase::~SolverBase ()
{}
-
+
SolverControl &
SolverBase::control() const
AssertThrow (ierr == 0, ExcTrilinosError(ierr));
// ... set some options, ...
- solver.SetAztecOption (AZ_output, AZ_none);
+ solver.SetAztecOption (AZ_output, additional_data.output_solver_details ?
+ AZ_all : AZ_none);
+ solver.SetAztecOption (AZ_conv, AZ_noscaled);
+
+ // ... and then solve!
+ ierr = solver.Iterate (solver_control.max_steps(),
+ solver_control.tolerance());
+ AssertThrow (ierr >= 0, ExcTrilinosError(ierr));
+
+ // Finally, let the deal.II
+ // SolverControl object know
+ // what has happened. If the
+ // solve succeeded, the status
+ // of the solver control will
+ // turn into
+ // SolverControl::success.
+ solver_control.check (solver.NumIters(), solver.TrueResidual());
+
+ if (solver_control.last_check() != SolverControl::success)
+ throw SolverControl::NoConvergence (solver_control.last_step(),
+ solver_control.last_value());
+ }
+
+
+
+ void
+ SolverBase::solve (const SparseMatrix &A,
+ dealii::Vector<double> &x,
+ const dealii::Vector<double> &b,
+ const PreconditionBase &preconditioner)
+ {
+ int ierr;
+
+ linear_problem.reset();
+
+ // In case we call the solver with
+ // deal.II vectors, we create views
+ // of the vectors in Epetra format.
+ Assert (x.size() == A.n(),
+ ExcDimensionMismatch(x.size(), A.n()));
+ Assert (b.size() == A.m(),
+ ExcDimensionMismatch(b.size(), A.m()));
+ Assert (A.local_range ().second == A.m(),
+ ExcMessage ("Can only work in serial when using deal.II vectors."));
+ Assert (A.matrix->Filled(),
+ ExcMessage ("Matrix is not compressed. Call compress() method."));
+
+ Epetra_Vector ep_x (View, A.matrix->DomainMap(), x.begin());
+ Epetra_Vector ep_b (View, A.matrix->RangeMap(), const_cast<double*>(b.begin()));
+
+ // We need an
+ // Epetra_LinearProblem object
+ // to let the AztecOO solver
+ // know about the matrix and
+ // vectors.
+ linear_problem = std::auto_ptr<Epetra_LinearProblem>
+ (new Epetra_LinearProblem(&*(A.matrix), &ep_x, &ep_b));
+
+ // Next we can allocate the
+ // AztecOO solver...
+ solver.SetProblem(*linear_problem);
+
+ // ... and we can specify the
+ // solver to be used.
+ switch (solver_name)
+ {
+ case cg:
+ solver.SetAztecOption(AZ_solver, AZ_cg);
+ break;
+ case cgs:
+ solver.SetAztecOption(AZ_solver, AZ_cgs);
+ break;
+ case gmres:
+ solver.SetAztecOption(AZ_solver, AZ_gmres);
+ solver.SetAztecOption(AZ_kspace, additional_data.gmres_restart_parameter);
+ break;
+ case bicgstab:
+ solver.SetAztecOption(AZ_solver, AZ_bicgstab);
+ break;
+ case tfqmr:
+ solver.SetAztecOption(AZ_solver, AZ_tfqmr);
+ break;
+ default:
+ Assert (false, ExcNotImplemented());
+ }
+
+ // Introduce the
+ // preconditioner, ...
+ ierr = solver.SetPrecOperator (const_cast<Epetra_Operator*>
+ (&*preconditioner.preconditioner));
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+ // ... set some options, ...
+ solver.SetAztecOption (AZ_output, additional_data.output_solver_details ?
+ AZ_all : AZ_none);
solver.SetAztecOption (AZ_conv, AZ_noscaled);
// ... and then solve!
/* ---------------------- SolverCG ------------------------ */
+ SolverCG::AdditionalData::
+ AdditionalData (const bool output_solver_details)
+ :
+ output_solver_details (output_solver_details)
+ {}
+
+
+
SolverCG::SolverCG (SolverControl &cn,
const AdditionalData &data)
:
SolverBase (cn),
- additional_data (data)
+ additional_data (data.output_solver_details)
{
solver_name = cg;
}
/* ---------------------- SolverGMRES ------------------------ */
SolverGMRES::AdditionalData::
- AdditionalData (const unsigned int restart_parameter)
+ AdditionalData (const bool output_solver_details,
+ const unsigned int restart_parameter)
:
+ output_solver_details (output_solver_details),
restart_parameter (restart_parameter)
{}
const AdditionalData &data)
:
SolverBase (cn),
- additional_data (data.restart_parameter)
+ additional_data (data.output_solver_details,
+ data.restart_parameter)
{
solver_name = gmres;
}
/* ---------------------- SolverBicgstab ------------------------ */
+ SolverBicgstab::AdditionalData::
+ AdditionalData (const bool output_solver_details)
+ :
+ output_solver_details (output_solver_details)
+ {}
+
+
+
+
SolverBicgstab::SolverBicgstab (SolverControl &cn,
const AdditionalData &data)
:
SolverBase (cn),
- additional_data (data)
+ additional_data (data.output_solver_details)
{
solver_name = bicgstab;
}
/* ---------------------- SolverCGS ------------------------ */
+ SolverCGS::AdditionalData::
+ AdditionalData (const bool output_solver_details)
+ :
+ output_solver_details (output_solver_details)
+ {}
+
+
+
+
SolverCGS::SolverCGS (SolverControl &cn,
const AdditionalData &data)
:
SolverBase (cn),
- additional_data (data)
+ additional_data (data.output_solver_details)
{
solver_name = cgs;
}
/* ---------------------- SolverTFQMR ------------------------ */
+ SolverTFQMR::AdditionalData::
+ AdditionalData (const bool output_solver_details)
+ :
+ output_solver_details (output_solver_details)
+ {}
+
+
+
SolverTFQMR::SolverTFQMR (SolverControl &cn,
const AdditionalData &data)
:
SolverBase (cn),
- additional_data (data)
+ additional_data (data.output_solver_details)
{
solver_name = tfqmr;
}
+
+
+/* ---------------------- SolverDirect ------------------------ */
+
+ SolverDirect::AdditionalData::
+ AdditionalData (const bool output_solver_details)
+ :
+ output_solver_details (output_solver_details)
+ {}
+
+
+
+
+ SolverDirect::SolverDirect (SolverControl &cn,
+ const AdditionalData &data)
+ :
+ solver_control (cn),
+ additional_data (data.output_solver_details)
+ {}
+
+
+
+ SolverDirect::~SolverDirect ()
+ {}
+
+
+
+ SolverControl &
+ SolverDirect::control() const
+ {
+ return solver_control;
+ }
+
+
+
+ void
+ SolverDirect::solve (const SparseMatrix &A,
+ VectorBase &x,
+ const VectorBase &b)
+ {
+ // First set whether we want to print
+ // the solver information to screen
+ // or not.
+ ConditionalOStream verbose_cout (std::cout,
+ additional_data.output_solver_details);
+
+ linear_problem.reset();
+ solver.reset();
+
+ // We need an
+ // Epetra_LinearProblem object
+ // to let the AztecOO solver
+ // know about the matrix and
+ // vectors.
+ linear_problem = std::auto_ptr<Epetra_LinearProblem> (
+ new Epetra_LinearProblem(&*(A.matrix), &*x.vector,
+ &*b.vector));
+
+ // Next we can allocate the
+ // AztecOO solver...
+ solver = std::auto_ptr<Amesos_BaseSolver> (Amesos().Create("Amesos_Klu",
+ *linear_problem));
+
+ verbose_cout << "Starting symbolic factorization" << std::endl;
+ solver->SymbolicFactorization();
+
+ verbose_cout << "Starting numeric factorization" << std::endl;
+ solver->NumericFactorization();
+
+ verbose_cout << "Starting solve" << std::endl;
+ solver->Solve();
+
+ // Finally, let the deal.II
+ // SolverControl object know
+ // what has happened. If the
+ // solve succeeded, the status
+ // of the solver control will
+ // turn into
+ // SolverControl::success.
+ solver_control.check (0, 0);
+
+ if (solver_control.last_check() != SolverControl::success)
+ throw SolverControl::NoConvergence (solver_control.last_step(),
+ solver_control.last_value());
+ }
+
+
+
+ void
+ SolverDirect::solve (const SparseMatrix &A,
+ dealii::Vector<double> &x,
+ const dealii::Vector<double> &b)
+ {
+ // First set whether we want to print
+ // the solver information to screen
+ // or not.
+ ConditionalOStream verbose_cout (std::cout,
+ additional_data.output_solver_details);
+
+ linear_problem.reset();
+ solver.reset();
+
+
+ // In case we call the solver with
+ // deal.II vectors, we create views
+ // of the vectors in Epetra format.
+ Assert (x.size() == A.n(),
+ ExcDimensionMismatch(x.size(), A.n()));
+ Assert (b.size() == A.m(),
+ ExcDimensionMismatch(b.size(), A.m()));
+ Assert (A.local_range ().second == A.m(),
+ ExcMessage ("Can only work in serial when using deal.II vectors."));
+ Epetra_Vector ep_x (View, A.matrix->DomainMap(), x.begin());
+ Epetra_Vector ep_b (View, A.matrix->RangeMap(), const_cast<double*>(b.begin()));
+
+ // We need an
+ // Epetra_LinearProblem object
+ // to let the AztecOO solver
+ // know about the matrix and
+ // vectors.
+ linear_problem = std::auto_ptr<Epetra_LinearProblem>
+ (new Epetra_LinearProblem(&*(A.matrix), &ep_x, &ep_b));
+
+ // Next we can allocate the
+ // AztecOO solver...
+ solver = std::auto_ptr<Amesos_BaseSolver> (Amesos().Create("Amesos_Klu",
+ *linear_problem));
+
+ verbose_cout << "Starting symbolic factorization" << std::endl;
+ solver->SymbolicFactorization();
+
+ verbose_cout << "Starting numeric factorization" << std::endl;
+ solver->NumericFactorization();
+
+ verbose_cout << "Starting solve" << std::endl;
+ solver->Solve();
+
+ // Finally, let the deal.II
+ // SolverControl object know
+ // what has happened. If the
+ // solve succeeded, the status
+ // of the solver control will
+ // turn into
+ // SolverControl::success.
+ solver_control.check (0, 0);
+
+ if (solver_control.last_check() != SolverControl::success)
+ throw SolverControl::NoConvergence (solver_control.last_step(),
+ solver_control.last_value());
+ }
+
+
+
+
}
DEAL_II_NAMESPACE_CLOSE
projects / dealii-svn.git / commitdiff