<h3>Specific improvements</h3>
<ol>
+ <li> Improved: The Trilinos direct solver, TrilinosWrappers::SolverDirect,
+ now takes a string to select among the available Amesos solvers. Moreover,
+ the solver now also supports deal.II's distributed vectors.
+ <br>
+ (Uwe Köcher, Martin Kronbichler, 2014/07/09)
+ </li>
+
<li> New: There are now three new preconditioner classes
TrilinosWrappers::PreconditionBlockJacobi,
TrilinosWrappers::PreconditionBlockSSOR, and
/**
* Determines which smoother to use for the AMG cycle. Possibilities
* for smoother_type are the following:
- * "Aztec"
- * "IFPACK"
- * "Jacobi"
- * "ML symmetric Gauss-Seidel"
- * "symmetric Gauss-Seidel"
- * "ML Gauss-Seidel"
- * "Gauss-Seidel"
- * "block Gauss-Seidel"
- * "symmetric block Gauss-Seidel"
- * "Chebyshev"
- * "MLS"
- * "Hiptmair"
- * "Amesos-KLU"
- * "Amesos-Superlu"
- * "Amesos-UMFPACK"
- * "Amesos-Superludist"
- * "Amesos-MUMPS"
- * "user-defined"
- * "SuperLU"
- * "IFPACK-Chebyshev"
- * "self"
- * "do-nothing"
- * "IC"
- * "ICT"
- * "ILU"
- * "ILUT"
- * "Block Chebyshev"
- * "IFPACK-Block Chebyshev"
+ * <ol>
+ * <li> "Aztec" </li>
+ * <li> "IFPACK" </li>
+ * <li> "Jacobi" </li>
+ * <li> "ML symmetric Gauss-Seidel" </li>
+ * <li> "symmetric Gauss-Seidel" </li>
+ * <li> "ML Gauss-Seidel" </li>
+ * <li> "Gauss-Seidel" </li>
+ * <li> "block Gauss-Seidel" </li>
+ * <li> "symmetric block Gauss-Seidel" </li>
+ * <li> "Chebyshev" </li>
+ * <li> "MLS" </li>
+ * <li> "Hiptmair" </li>
+ * <li> "Amesos-KLU" </li>
+ * <li> "Amesos-Superlu" </li>
+ * <li> "Amesos-UMFPACK" </li>
+ * <li> "Amesos-Superludist" </li>
+ * <li> "Amesos-MUMPS" </li>
+ * <li> "user-defined" </li>
+ * <li> "SuperLU" </li>
+ * <li> "IFPACK-Chebyshev" </li>
+ * <li> "self" </li>
+ * <li> "do-nothing" </li>
+ * <li> "IC" </li>
+ * <li> "ICT" </li>
+ * <li> "ILU" </li>
+ * <li> "ILUT" </li>
+ * <li> "Block Chebyshev" </li>
+ * <li> "IFPACK-Block Chebyshev" </li>
+ * </ol>
*/
const char* smoother_type;
# include <deal.II/lac/exceptions.h>
# include <deal.II/lac/solver_control.h>
# include <deal.II/lac/vector.h>
+# include <deal.II/lac/parallel_vector.h>
# include <Epetra_LinearProblem.h>
# include <AztecOO.h>
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:
+ * 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;
struct AdditionalData
{
/**
- * Sets the additional data field to
- * the desired output format and puts
- * the 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
- * 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.
+ * TODO: Find a better way for 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
+ * Enables/disables the output of solver details (residual in each
* iterations etc.).
*/
const bool output_solver_details;
/**
- * Restart parameter for GMRES
- * solver.
+ * Restart parameter for GMRES solver.
*/
const unsigned int gmres_restart_parameter;
};
/**
- * Constructor. Takes the
- * solver control object and
- * creates the solver.
+ * Constructor. Takes the solver control object and 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.
+ * 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);
virtual ~SolverBase ();
/**
- * 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.
+ * 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.
*/
void
solve (const SparseMatrix &A,
const PreconditionBase &preconditioner);
/**
- * Solve the linear system
- * <tt>Ax=b</tt> where <tt>A</tt>
- * is an operator. This function
- * can be used for matrix free
- * computation. 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.
+ * Solve the linear system <tt>Ax=b</tt> where <tt>A</tt> is an
+ * operator. This function can be used for matrix free
+ * computation. 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.
*/
void
solve (Epetra_Operator &A,
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.
+ * 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,
const PreconditionBase &preconditioner);
/**
- * Solve the linear system
- * <tt>Ax=b</tt> where <tt>A</tt>
- * is an operator. This function can
- * be used for matric free. 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.
+ * Solve the linear system <tt>Ax=b</tt> where <tt>A</tt> is an
+ * operator. This function can be used for matrix free
+ * computations. 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 (Epetra_Operator &A,
const PreconditionBase &preconditioner);
/**
- * Access to object that controls
- * convergence.
+ * Solve the linear system <tt>Ax=b</tt> for deal.II's parallel
+ * distributed vectors. 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.
+ */
+ void
+ solve (const SparseMatrix &A,
+ dealii::parallel::distributed::Vector<double> &x,
+ const dealii::parallel::distributed::Vector<double> &b,
+ const PreconditionBase &preconditioner);
+
+ /**
+ * Solve the linear system <tt>Ax=b</tt> where <tt>A</tt> is an
+ * operator. This function can be used for matrix free
+ * computation. 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.
+ */
+ void
+ solve (Epetra_Operator &A,
+ dealii::parallel::distributed::Vector<double> &x,
+ const dealii::parallel::distributed::Vector<double> &b,
+ const PreconditionBase &preconditioner);
+
+
+ /**
+ * Access to object that controls convergence.
*/
SolverControl &control() const;
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.
+ * 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;
* The solve function is used to set properly the Epetra_LinearProblem,
* once it is done this function solves the linear problem.
*/
- void execute_solve(const PreconditionBase &preconditioner);
+ void do_solve(const PreconditionBase &preconditioner);
/**
- * A structure that collects
- * the Trilinos sparse matrix,
- * the right hand side vector
- * and the solution vector,
- * which is passed down to the
+ * 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_cxx1x::shared_ptr<Epetra_LinearProblem> linear_problem;
/**
- * A structure that contains
- * the Trilinos solver and
- * preconditioner objects.
+ * A structure that contains the Trilinos solver and preconditioner
+ * objects.
*/
AztecOO solver;
/**
- * Store a copy of the flags for this
- * particular solver.
+ * Store a copy of the flags for this particular solver.
*/
const AdditionalData additional_data;
{
public:
/**
- * Standardized data struct to
- * pipe additional data to the
- * solver.
+ * Standardized data struct to pipe additional data to the solver.
*/
struct AdditionalData
{
/**
- * Sets the additional data field to
- * the desired output format.
+ * 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
+ * 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
- * need to give a vector memory
- * object.
+ * 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
+ * The last argument takes a structure with additional, solver dependent
* flags for tuning.
*/
SolverCG (SolverControl &cn,
protected:
/**
- * Store a copy of the flags for this
- * particular solver.
+ * Store a copy of the flags for this particular solver.
*/
const AdditionalData additional_data;
};
{
public:
/**
- * Standardized data struct to
- * pipe additional data to the
- * solver.
+ * Standardized data struct to pipe additional data to the solver.
*/
struct AdditionalData
{
/**
- * Sets the additional data field to
- * the desired output format.
+ * 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
+ * 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
- * need to give a vector memory
- * object.
+ * 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
+ * The last argument takes a structure with additional, solver dependent
* flags for tuning.
*/
SolverCGS (SolverControl &cn,
protected:
/**
- * Store a copy of the flags for this
- * particular solver.
+ * Store a copy of the flags for this particular solver.
*/
const AdditionalData additional_data;
};
{
public:
/**
- * Standardized data struct to
- * pipe additional data to the
- * solver.
+ * Standardized data struct to pipe additional data to the solver.
*/
struct AdditionalData
{
/**
- * Constructor. By default, set the
- * number of temporary vectors to
- * 30, i.e. do a restart every 30
- * iterations.
+ * Constructor. By default, set the number of temporary vectors to 30,
+ * i.e. do a restart every 30 iterations.
*/
AdditionalData (const bool output_solver_details = false,
const unsigned int restart_parameter = 30);
/**
- * Enables/disables the output of
- * solver details (residual in each
+ * Enables/disables the output of solver details (residual in each
* iterations etc.).
*/
bool output_solver_details;
/**
- * Maximum number of
- * tmp vectors.
+ * Maximum number of tmp vectors.
*/
unsigned int restart_parameter;
};
/**
- * Constructor. In contrast to
- * deal.II's own solvers, there is no
- * need to give a vector memory
- * object.
+ * 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
+ * The last argument takes a structure with additional, solver dependent
* flags for tuning.
*/
SolverGMRES (SolverControl &cn,
protected:
/**
- * Store a copy of the flags for this
- * particular solver.
+ * Store a copy of the flags for this particular solver.
*/
const AdditionalData additional_data;
};
{
public:
/**
- * Standardized data struct to
- * pipe additional data to the
- * solver.
+ * Standardized data struct to pipe additional data to the solver.
*/
struct AdditionalData
{
/**
- * Sets the additional data field to
- * the desired output format.
+ * 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
+ * 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
- * need to give a vector memory
- * object.
+ * 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
+ * The last argument takes a structure with additional, solver dependent
* flags for tuning.
*/
SolverBicgstab (SolverControl &cn,
protected:
/**
- * Store a copy of the flags for this
- * particular solver.
+ * Store a copy of the flags for this particular solver.
*/
const AdditionalData additional_data;
};
{
public:
/**
- * Standardized data struct to
- * pipe additional data to the
- * solver.
+ * Standardized data struct to pipe additional data to the solver.
*/
struct AdditionalData
{
/**
- * Sets the additional data field to
- * the desired output format.
+ * 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
+ * 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
- * need to give a vector memory
- * object.
+ * 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
+ * The last argument takes a structure with additional, solver dependent
* flags for tuning.
*/
SolverTFQMR (SolverControl &cn,
protected:
/**
- * Store a copy of the flags for this
- * particular solver.
+ * Store a copy of the flags for this particular solver.
*/
const AdditionalData additional_data;
};
public:
/**
- * Standardized data struct to
- * pipe additional data to the
- * solver.
+ * Standardized data struct to pipe additional data to the solver.
*/
struct AdditionalData
{
/**
- * Sets the additional data field to
- * the desired output format.
+ * Sets the additional data field to the desired output format.
*/
- AdditionalData (const bool output_solver_details = false);
+ AdditionalData (const bool output_solver_details = false,
+ const std::string &solver_type = "Amesos_Klu");
/**
- * Enables/disables the output of
- * solver details (residual in each
+ * Enables/disables the output of solver details (residual in each
* iterations etc.).
*/
bool output_solver_details;
+
+ /**
+ * Set the solver type (for third party solver support of Trilinos
+ * Amesos package). Possibilities are:
+ * <ol>
+ * <li> "Amesos_Lapack" </li>
+ * <li> "Amesos_Scalapack" </li>
+ * <li> "Amesos_Klu" </li>
+ * <li> "Amesos_Umfpack" </li>
+ * <li> "Amesos_Pardiso" </li>
+ * <li> "Amesos_Taucs" </li>
+ * <li> "Amesos_Superlu" </li>
+ * <li> "Amesos_Superludist" </li>
+ * <li> "Amesos_Dscpack" </li>
+ * <li> "Amesos_Mumps" </li>
+ * </ol>
+ * Note that the availability of these solvers in deal.II depends on
+ * which solvers were set when configuring Trilinos.
+ */
+ std::string solver_type;
};
/**
- * Constructor. Takes the
- * solver control object and
- * creates the solver.
+ * Constructor. Takes the solver control object and creates the solver.
*/
SolverDirect (SolverControl &cn,
const AdditionalData &data = AdditionalData());
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.
+ * Solve the linear system <tt>Ax=b</tt>. Creates a factorization of the
+ * matrix with the package chosen from the additional data structure and
+ * performs the solve. Note that there is no need for a preconditioner
+ * here.
*/
void
solve (const SparseMatrix &A,
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.
+ * Solve the linear system <tt>Ax=b</tt>. This class works with Trilinos
+ * matrices, but takes deal.II serial vectors as argument. Since these
+ * vectors are not distributed, this function does only what you expect in
+ * case the matrix is serial (i.e., locally owned). Otherwise, an
+ * exception will be thrown.
*/
void
solve (const SparseMatrix &A,
dealii::Vector<double> &x,
const dealii::Vector<double> &b);
+ /**
+ * Solve the linear system <tt>Ax=b</tt> for deal.II's own parallel
+ * vectors. Creates a factorization of the matrix with the package chosen
+ * from the additional data structure and performs the solve. Note that
+ * there is no need for a preconditioner here.
+ */
+ void
+ solve (const SparseMatrix &A,
+ dealii::parallel::distributed::Vector<double> &x,
+ const dealii::parallel::distributed::Vector<double> &b);
+
/**
* Access to object that controls
* convergence.
<< " occurred while calling a Trilinos function");
private:
+ /**
+ * Actually performs the operations for solving the linear system,
+ * including the factorization and forward and backward substitution.
+ */
+ void do_solve();
/**
- * 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.
+ * 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
+ * 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_cxx1x::shared_ptr<Epetra_LinearProblem> linear_problem;
/**
- * A structure that contains
- * the Trilinos solver and
- * preconditioner objects.
+ * A structure that contains the Trilinos solver and preconditioner
+ * objects.
*/
std_cxx1x::shared_ptr<Amesos_BaseSolver> solver;
/**
- * Store a copy of the flags for this
- * particular solver.
+ * Store a copy of the flags for this particular solver.
*/
const AdditionalData additional_data;
-
};
}
{
linear_problem.reset();
- // We need an
- // Epetra_LinearProblem object
- // to let the AztecOO solver
- // know about the matrix and
- // vectors.
+ // We need an Epetra_LinearProblem object to let the AztecOO solver know
+ // about the matrix and vectors.
linear_problem.reset
(new Epetra_LinearProblem(const_cast<Epetra_CrsMatrix *>(&A.trilinos_matrix()),
&x.trilinos_vector(),
const_cast<Epetra_MultiVector *>(&b.trilinos_vector())));
- execute_solve(preconditioner);
+ do_solve(preconditioner);
}
{
linear_problem.reset();
- // We need an
- // Epetra_LinearProblem object
- // to let the AztecOO solver
- // know about the matrix and
- // vectors.
+ // We need an Epetra_LinearProblem object to let the AztecOO solver know
+ // about the matrix and vectors.
linear_problem.reset
(new Epetra_LinearProblem(&A,
&x.trilinos_vector(),
const_cast<Epetra_MultiVector *>(&b.trilinos_vector())));
- execute_solve(preconditioner);
+ do_solve(preconditioner);
}
{
linear_problem.reset();
- // In case we call the solver with
- // deal.II vectors, we create views
- // of the vectors in Epetra format.
+ // 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(),
Epetra_Vector ep_x (View, A.domain_partitioner(), x.begin());
Epetra_Vector ep_b (View, A.range_partitioner(), const_cast<double *>(b.begin()));
- // We need an
- // Epetra_LinearProblem object
- // to let the AztecOO solver
- // know about the matrix and
- // vectors.
+ // We need an Epetra_LinearProblem object to let the AztecOO solver know
+ // about the matrix and vectors.
linear_problem.reset (new Epetra_LinearProblem
(const_cast<Epetra_CrsMatrix *>(&A.trilinos_matrix()),
&ep_x, &ep_b));
- execute_solve(preconditioner);
+ do_solve(preconditioner);
}
Epetra_Vector ep_x (View, A.OperatorDomainMap(), x.begin());
Epetra_Vector ep_b (View, A.OperatorRangeMap(), const_cast<double *>(b.begin()));
- // We need an
- // Epetra_LinearProblem object
- // to let the AztecOO solver
- // know about the matrix and
- // vectors.
+ // We need an Epetra_LinearProblem object to let the AztecOO solver know
+ // about the matrix and vectors.
linear_problem.reset (new Epetra_LinearProblem(&A,&ep_x, &ep_b));
- execute_solve(preconditioner);
+ do_solve(preconditioner);
}
void
- SolverBase::execute_solve(const PreconditionBase &preconditioner)
+ SolverBase::solve (const SparseMatrix &A,
+ dealii::parallel::distributed::Vector<double> &x,
+ const dealii::parallel::distributed::Vector<double> &b,
+ const PreconditionBase &preconditioner)
+ {
+ linear_problem.reset();
+
+ // In case we call the solver with deal.II vectors, we create views of the
+ // vectors in Epetra format.
+ AssertDimension (static_cast<TrilinosWrappers::types::int_type>(x.local_size()),
+ A.domain_partitioner().NumMyElements());
+ AssertDimension (static_cast<TrilinosWrappers::types::int_type>(b.local_size()),
+ A.range_partitioner().NumMyElements());
+
+ Epetra_Vector ep_x (View, A.domain_partitioner(), x.begin());
+ Epetra_Vector ep_b (View, A.range_partitioner(), const_cast<double *>(b.begin()));
+
+ // We need an Epetra_LinearProblem object to let the AztecOO solver know
+ // about the matrix and vectors.
+ linear_problem.reset (new Epetra_LinearProblem
+ (const_cast<Epetra_CrsMatrix *>(&A.trilinos_matrix()),
+ &ep_x, &ep_b));
+
+ do_solve(preconditioner);
+ }
+
+
+
+ void
+ SolverBase::solve (Epetra_Operator &A,
+ dealii::parallel::distributed::Vector<double> &x,
+ const dealii::parallel::distributed::Vector<double> &b,
+ const PreconditionBase &preconditioner)
+ {
+ linear_problem.reset();
+
+ AssertDimension (static_cast<TrilinosWrappers::types::int_type>(x.local_size()),
+ A.OperatorDomainMap().NumMyElements());
+ AssertDimension (static_cast<TrilinosWrappers::types::int_type>(b.local_size()),
+ A.OperatorRangeMap().NumMyElements());
+
+ Epetra_Vector ep_x (View, A.OperatorDomainMap(), x.begin());
+ Epetra_Vector ep_b (View, A.OperatorRangeMap(), const_cast<double *>(b.begin()));
+
+ // We need an Epetra_LinearProblem object to let the AztecOO solver know
+ // about the matrix and vectors.
+ linear_problem.reset (new Epetra_LinearProblem(&A,&ep_x, &ep_b));
+
+ do_solve(preconditioner);
+ }
+
+
+
+ void
+ SolverBase::do_solve(const PreconditionBase &preconditioner)
{
int ierr;
- // Next we can allocate the
- // AztecOO solver...
+ // Next we can allocate the AztecOO solver...
solver.SetProblem(*linear_problem);
- // ... and we can specify the
- // solver to be used.
+ // ... and we can specify the solver to be used.
switch (solver_name)
{
case cg:
Assert (false, ExcNotImplemented());
}
- // Introduce the preconditioner,
- // if the identity preconditioner is used,
+ // Introduce the preconditioner, if the identity preconditioner is used,
// the precondioner is set to none, ...
if (preconditioner.preconditioner.use_count()!=0)
{
ierr = solver.Iterate (solver_control.max_steps(),
solver_control.tolerance());
- // report errors in more detail
- // than just by checking whether
- // the return status is zero or
- // greater. the error strings are
- // taken from the implementation
- // of the AztecOO::Iterate
- // function
+ // report errors in more detail than just by checking whether the return
+ // status is zero or greater. the error strings are taken from the
+ // implementation of the AztecOO::Iterate function
switch (ierr)
{
case -1:
"loss of precision"));
case -4:
AssertThrow (false, ExcMessage("AztecOO::Iterate error code -4: "
- "GMRES hessenberg ill-conditioned"));
+ "GMRES Hessenberg ill-conditioned"));
default:
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.
+ // 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)
/* ---------------------- SolverDirect ------------------------ */
SolverDirect::AdditionalData::
- AdditionalData (const bool output_solver_details)
+ AdditionalData (const bool output_solver_details,
+ const std::string &solver_type)
:
- output_solver_details (output_solver_details)
+ output_solver_details (output_solver_details),
+ solver_type(solver_type)
{}
const AdditionalData &data)
:
solver_control (cn),
- additional_data (data.output_solver_details)
+ additional_data (data.output_solver_details,data.solver_type)
{}
void
- SolverDirect::solve (const SparseMatrix &A,
- VectorBase &x,
- const VectorBase &b)
+ SolverDirect::do_solve()
{
- // First set whether we want to print
- // the solver information to screen
- // or not.
+ // Fetch return value of Amesos Solver functions
+ int ierr;
+
+ // 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.reset
- (new Epetra_LinearProblem(const_cast<Epetra_CrsMatrix *>(&A.trilinos_matrix()),
- &x.trilinos_vector(),
- const_cast<Epetra_MultiVector *>(&b.trilinos_vector())));
+ // Next allocate the Amesos solver, this is done in two steps, first we
+ // create a solver Factory and and generate with that the concrete Amesos
+ // solver, if possible.
+ Amesos Factory;
+
+ AssertThrow(
+ Factory.Query(additional_data.solver_type.c_str()),
+ ExcMessage (std::string ("You tried to select the solver type <") +
+ additional_data.solver_type +
+ "> but this solver is not supported by Trilinos either "
+ "because it does not exist, or because Trilinos was not "
+ "configured for its use.")
+ );
- // Next we can allocate the
- // AztecOO solver...
- solver.reset (Amesos().Create("Amesos_Klu", *linear_problem));
+ solver.reset (
+ Factory.Create(additional_data.solver_type.c_str(), *linear_problem)
+ );
verbose_cout << "Starting symbolic factorization" << std::endl;
- solver->SymbolicFactorization();
+ ierr = solver->SymbolicFactorization();
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
verbose_cout << "Starting numeric factorization" << std::endl;
- solver->NumericFactorization();
+ ierr = solver->NumericFactorization();
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
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.
+ ierr = solver->Solve();
+ 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 (0, 0);
if (solver_control.last_check() != SolverControl::success)
}
+ void
+ SolverDirect::solve (const SparseMatrix &A,
+ VectorBase &x,
+ const VectorBase &b)
+ {
+ // We need an Epetra_LinearProblem object to let the Amesos solver know
+ // about the matrix and vectors.
+ linear_problem.reset
+ (new Epetra_LinearProblem(const_cast<Epetra_CrsMatrix *>(&A.trilinos_matrix()),
+ &x.trilinos_vector(),
+ const_cast<Epetra_MultiVector *>(&b.trilinos_vector())));
+
+ do_solve();
+ }
+
+
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.
+ // 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(),
Epetra_Vector ep_x (View, A.domain_partitioner(), x.begin());
Epetra_Vector ep_b (View, A.range_partitioner(), const_cast<double *>(b.begin()));
- // We need an
- // Epetra_LinearProblem object
- // to let the AztecOO solver
- // know about the matrix and
- // vectors.
+ // We need an Epetra_LinearProblem object to let the Amesos solver know
+ // about the matrix and vectors.
linear_problem.reset (new Epetra_LinearProblem
(const_cast<Epetra_CrsMatrix *>(&A.trilinos_matrix()),
&ep_x, &ep_b));
- // Next we can allocate the
- // AztecOO solver...
- solver.reset (Amesos().Create("Amesos_Klu", *linear_problem));
+ do_solve();
+ }
- 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)
- AssertThrow(false, SolverControl::NoConvergence (solver_control.last_step(),
- solver_control.last_value()));
- }
+ void
+ SolverDirect::solve (const SparseMatrix &A,
+ dealii::parallel::distributed::Vector<double> &x,
+ const dealii::parallel::distributed::Vector<double> &b)
+ {
+ AssertDimension (static_cast<TrilinosWrappers::types::int_type>(x.local_size()),
+ A.domain_partitioner().NumMyElements());
+ AssertDimension (static_cast<TrilinosWrappers::types::int_type>(b.local_size()),
+ A.range_partitioner().NumMyElements());
+ Epetra_Vector ep_x (View, A.domain_partitioner(), x.begin());
+ Epetra_Vector ep_b (View, A.range_partitioner(), const_cast<double *>(b.begin()));
+ // We need an Epetra_LinearProblem object to let the Amesos solver know
+ // about the matrix and vectors.
+ linear_problem.reset (new Epetra_LinearProblem
+ (const_cast<Epetra_CrsMatrix *>(&A.trilinos_matrix()),
+ &ep_x, &ep_b));
+ do_solve();
+ }
}
--- /dev/null
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2013 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+// tests Trilinos direct solvers on a 2D Poisson equation for linear elements
+
+#include "../tests.h"
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/lac/compressed_sparsity_pattern.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/trilinos_solver.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/lac/trilinos_precondition.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/base/function.h>
+#include <deal.II/grid/tria.h>
+
+#include <fstream>
+#include <iomanip>
+
+
+template <int dim>
+class Step4
+{
+public:
+ Step4 ();
+ void run ();
+
+private:
+ void make_grid ();
+ void setup_system();
+ void assemble_system ();
+ void solve ();
+
+ Triangulation<dim> triangulation;
+ FE_Q<dim> fe;
+ DoFHandler<dim> dof_handler;
+
+ ConstraintMatrix constraints;
+
+ TrilinosWrappers::SparseMatrix system_matrix;
+
+ Vector<double> solution;
+ Vector<double> system_rhs;
+};
+
+
+template <int dim>
+class RightHandSide : public Function<dim>
+{
+public:
+ RightHandSide () : Function<dim>() {}
+
+ virtual double value (const Point<dim> &p,
+ const unsigned int component = 0) const;
+};
+
+
+
+template <int dim>
+class BoundaryValues : public Function<dim>
+{
+public:
+ BoundaryValues () : Function<dim>() {}
+
+ virtual double value (const Point<dim> &p,
+ const unsigned int component = 0) const;
+};
+
+
+
+
+template <int dim>
+double RightHandSide<dim>::value (const Point<dim> &p,
+ const unsigned int /*component*/) const
+{
+ double return_value = 0;
+ for (unsigned int i=0; i<dim; ++i)
+ return_value += 4*std::pow(p(i), 4);
+
+ return return_value;
+}
+
+
+
+template <int dim>
+double BoundaryValues<dim>::value (const Point<dim> &p,
+ const unsigned int /*component*/) const
+{
+ return p.square();
+}
+
+
+
+template <int dim>
+Step4<dim>::Step4 ()
+ :
+ fe (1),
+ dof_handler (triangulation)
+{}
+
+
+template <int dim>
+void Step4<dim>::make_grid ()
+{
+ GridGenerator::hyper_cube (triangulation, -1, 1);
+ triangulation.refine_global (6);
+}
+
+
+
+template <int dim>
+void Step4<dim>::setup_system ()
+{
+ dof_handler.distribute_dofs (fe);
+
+ constraints.clear();
+ std::map<unsigned int,double> boundary_values;
+ VectorTools::interpolate_boundary_values (dof_handler,
+ 0,
+ BoundaryValues<dim>(),
+ constraints);
+ constraints.close();
+
+ CompressedSparsityPattern c_sparsity(dof_handler.n_dofs());
+ DoFTools::make_sparsity_pattern (dof_handler, c_sparsity, constraints, false);
+ system_matrix.reinit (c_sparsity);
+
+ solution.reinit (dof_handler.n_dofs());
+ system_rhs.reinit (dof_handler.n_dofs());
+}
+
+
+template <int dim>
+void Step4<dim>::assemble_system ()
+{
+ QGauss<dim> quadrature_formula(fe.degree+1);
+
+ const RightHandSide<dim> right_hand_side;
+
+ FEValues<dim> fe_values (fe, quadrature_formula,
+ update_values | update_gradients |
+ update_quadrature_points | update_JxW_values);
+
+ const unsigned int dofs_per_cell = fe.dofs_per_cell;
+ const unsigned int n_q_points = quadrature_formula.size();
+
+ FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
+ Vector<double> cell_rhs (dofs_per_cell);
+
+ std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
+
+ typename DoFHandler<dim>::active_cell_iterator
+ cell = dof_handler.begin_active(),
+ endc = dof_handler.end();
+
+ for (; cell!=endc; ++cell)
+ {
+ fe_values.reinit (cell);
+ cell_matrix = 0;
+ cell_rhs = 0;
+
+ for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ {
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ cell_matrix(i,j) += (fe_values.shape_grad (i, q_point) *
+ fe_values.shape_grad (j, q_point) *
+ fe_values.JxW (q_point));
+
+ cell_rhs(i) += (fe_values.shape_value (i, q_point) *
+ right_hand_side.value (fe_values.quadrature_point (q_point)) *
+ fe_values.JxW (q_point));
+ }
+
+ cell->get_dof_indices (local_dof_indices);
+ constraints.distribute_local_to_global(cell_matrix, cell_rhs,
+ local_dof_indices,
+ system_matrix, system_rhs);
+ }
+ system_matrix.compress(VectorOperation::add);
+}
+
+
+
+template <int dim>
+void Step4<dim>::solve ()
+{
+ // Compute 'reference' solution with CG solver and SSOR preconditioner
+ TrilinosWrappers::PreconditionSSOR preconditioner;
+ solution = 0;
+ SolverControl solver_control (1000, 1e-12);
+ SolverCG<> solver (solver_control);
+ preconditioner.initialize(system_matrix);
+ solver.solve (system_matrix, solution, system_rhs,
+ preconditioner);
+
+ Vector<double> output(solution);
+ {
+ deallog.push("DirectKLU");
+ TrilinosWrappers::SolverDirect::AdditionalData data;
+ data.solver_type = "Amesos_Klu";
+ SolverControl solver_control (1000, 1e-10);
+ TrilinosWrappers::SolverDirect solver(solver_control, data);
+ solver.solve (system_matrix, output, system_rhs);
+ output -= solution;
+ deallog << "Norm of error in direct solve: " << output.l2_norm()
+ << std::endl;
+ deallog.pop();
+ }
+
+ {
+ deallog.push("DirectLAPACK");
+ TrilinosWrappers::SolverDirect::AdditionalData data;
+ data.solver_type = "Amesos_Lapack";
+ SolverControl solver_control (1000, 1e-10);
+ TrilinosWrappers::SolverDirect solver(solver_control, data);
+ solver.solve (system_matrix, output, system_rhs);
+ output -= solution;
+ deallog << "Norm of error in direct solve: " << output.l2_norm()
+ << std::endl;
+ deallog.pop();
+ }
+
+ {
+ deallog.push("DirectDscpack");
+ TrilinosWrappers::SolverDirect::AdditionalData data;
+ data.solver_type = "Amesos_Dscpack";
+ SolverControl solver_control (1000, 1e-10);
+ TrilinosWrappers::SolverDirect solver(solver_control, data);
+ try
+ {
+ solver.solve (system_matrix, output, system_rhs);
+ }
+ catch (dealii::ExceptionBase &exc)
+ {
+ deallog << "Error: " << std::endl;
+ exc.print_info(deallog.get_file_stream());
+ }
+ deallog.pop();
+ }
+
+ {
+ deallog.push("Dummy");
+ TrilinosWrappers::SolverDirect::AdditionalData data;
+ data.solver_type = "DummySolver";
+ SolverControl solver_control (1000, 1e-10);
+ TrilinosWrappers::SolverDirect solver(solver_control, data);
+ try
+ {
+ solver.solve (system_matrix, output, system_rhs);
+ }
+ catch (dealii::ExceptionBase &exc)
+ {
+ deallog << "Error: " << std::endl;
+ exc.print_info(deallog.get_file_stream());
+ }
+ deallog.pop();
+ }
+}
+
+
+
+template <int dim>
+void Step4<dim>::run()
+{
+ make_grid();
+ setup_system();
+ assemble_system();
+ solve();
+}
+
+
+int main (int argc, char **argv)
+{
+ std::ofstream logfile("output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv);
+
+ try
+ {
+ Step4<2> test;
+ test.run();
+ }
+ catch (std::exception &exc)
+ {
+ deallog << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ deallog << "Exception on processing: " << std::endl
+ << exc.what() << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+
+ return 1;
+ }
+ catch (...)
+ {
+ deallog << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ deallog << "Unknown exception!" << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ return 1;
+ };
+}
--- /dev/null
+
+DEAL:cg::Starting value 21.8299
+DEAL:cg::Convergence step 86 value 0
+DEAL:DirectKLU::Starting value 0
+DEAL:DirectKLU::Convergence step 0 value 0
+DEAL:DirectKLU::Norm of error in direct solve: 0
+DEAL:DirectLAPACK::Starting value 0
+DEAL:DirectLAPACK::Convergence step 0 value 0
+DEAL:DirectLAPACK::Norm of error in direct solve: 0
+DEAL:DirectDscpack::Error:
+You tried to select the solver type <Amesos_Dscpack> but this solver is not supported by Trilinos either because it does not exist, or because Trilinos was not configured for its use.
+DEAL:Dummy::Error:
+You tried to select the solver type <DummySolver> but this solver is not supported by Trilinos either because it does not exist, or because Trilinos was not configured for its use.
projects / dealii-svn.git / commitdiff