* its own.
*
* The BlockVector is a collection of Vectors (e.g. of either deal.II Vector
-- * objects or PETScWrappers::Vector object). Each of the vectors inside can
-- * have a different size.
++ * objects, PETScWrappers::Vector object, TrilinosWrappers::Vector, or
++ * ParalutionWrappers::Vector). Each of the vectors inside can have a different
++ * size.
*
* The functionality of BlockVector includes everything a Vector can do, plus
* the access to a single Vector inside the BlockVector by
--- /dev/null
- //TODO: preconditioners are missing
+// ---------------------------------------------------------------------
+// $Id: paralution_precondition.h 30040 2013-07-18 17:06:48Z maier $
+//
+// Copyright (C) 2014 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.
+//
+// ---------------------------------------------------------------------
+
+#ifndef __deal2__paralution_precondition_h
+#define __deal2__paralution_precondition_h
+
+
+#include <deal.II/base/config.h>
+
+#ifdef DEAL_II_WITH_PARALUTION
+
++#include <deal.II/base/exceptions.h>
+#include <deal.II/base/std_cxx1x/shared_ptr.h>
+
++#include <vector>
++
+#include <paralution.hpp>
+
+DEAL_II_NAMESPACE_OPEN
+
+
+namespace ParalutionWrappers
+{
++ /**
++ * Cycle for multigrid preconditioner.
++ */
++ enum mg_cycle {V_cycle, W_cycle, K_cycle, F_cycle};
++
++ /**
++ * Interpolation for grid transfer operator.
++ */
++ enum mg_interpolation {aggregation, smoothed_aggregation};
++
++ /**
++ * Precondition for the smoother.
++ */
++ enum mg_preconditioner {jacobi, sgs, multicolored_sgs, multicolored_sor, ilu,
++ ilut, multicolored_ilu
++ };
++
++ /**
++ * Solver used as smoother and to solve the coarse system.
++ */
++ enum mg_solver {richardson, cg, cr, bicgstab, gmres, fgmres};
++
+ /**
+ * The base class for the preconditioner classes using the Paralution
+ * functionality. In Paralution, the preconditioners are not built at the same
+ * time as the solver. Therefore, only the declaration is done in the
+ * preconditiner classes.
+ *
+ * @ingroup ParalutionWrappers
+ * @ingroup Preconditioners
+ * @author Bruno Turcksin, 2014
+ */
- std_cxx1x::shared_ptr<paralution::Preconditioner<paralution::LocalMatrix<Number>,
+ template <typename Number>
+ class PreconditionBase
+ {
+ public :
+ /**
+ * Standardized data struct to pipe additional flags to the
+ * preconditioner.
+ */
+ struct AdditionalData
+ {};
+
+ /**
+ * Constructor. Does not do anything.
+ */
+ PreconditionBase ();
+
+ friend class SolverBase;
+
+ protected :
+ /**
+ * This is a pointer to the preconditioner object that is used when
+ * applying the preconditioner.
+ */
-
- private :
- /**
- * Store a copy of the flags for this
- * particular preconditioner.
- */
- AdditionalData additional_data;
++ std_cxx1x::shared_ptr<paralution::Solver<paralution::LocalMatrix<Number>,
+ paralution::LocalVector<Number>,Number> > preconditioner;
+ };
+
+
+
+ /**
+ * A wrapper for Jacobi preconditioner for Paralution matrices.
+ *
+ * @ingroup ParalutionWrappers
+ * @ingroup Preconditioners
+ * @author Bruno Turcksin, 2014
+ */
+ template <typename Number>
+ class PreconditionJacobi : public PreconditionBase<Number>
+ {
+ public :
+ /**
+ * Standardized data struct to pipe additional flags to the
+ * preconditioner.
+ */
+ struct AdditionalData
+ {};
+
+ /**
+ * Constructor.
+ */
+ PreconditionJacobi ();
-
- private :
- /**
- * Store a copy of the flags for this
- * particular preconditioner.
- */
- AdditionalData additional_data;
-
+ };
+
+
+
+ /**
+ * A wrapper for Symmetric Gauss-Seidel for Paralution matrices.
+ *
+ * @ingroup ParalutionWrappers
+ * @ingroup Preconditioners
+ * @author Bruno Turcksin, 2014
+ */
+ template <typename Number>
+ class PreconditionSGS : public PreconditionBase<Number>
+ {
+ public :
+ /**
+ * Standardized data struct to pipe additional flags to the
+ * preconditioner.
+ */
+ struct AdditionalData
+ {};
+
+ /**
+ * Constructor.
+ */
+ PreconditionSGS ();
-
- private :
- /**
- * Store a copy of the flags for this
- * particular preconditioner.
- */
- AdditionalData additional_data;
-
+ };
+
+
+
+ /**
+ * A wrapper for multi-colored Symmetric Gauss-Seidel for Paralution matrices.
+ *
+ * @ingroup ParalutionWrappers
+ * @ingroup Preconditioners
+ * @author Bruno Turcksin, 2014
+ */
+ template <typename Number>
+ class PreconditionMultiColoredSGS : public PreconditionBase<Number>
+ {
+ public :
+ /**
+ * Standardized data struct to pipe additional flags to the
+ * preconditioner.
+ */
+ struct AdditionalData
+ {};
+
+ /**
+ * Constructor.
+ */
+ PreconditionMultiColoredSGS ();
-
- private :
- /**
- * Store a copy of the flags for this particular preconditioner.
- */
- AdditionalData additional_data;
+ };
+
+
+
+ /**
+ * A wrapper for multi-colored SOR for Paralution matrices.
+ *
+ * @ingroup ParalutionWrappers
+ * @ingroup Preconditioners
+ * @author Bruno Turcksin, 2014
+ */
+ template <typename Number>
+ class PreconditionMultiColoredSOR : public PreconditionBase<Number>
+ {
+ public :
+ /**
+ * Standardized data struct to pipe additional flags to the
+ * preconditioner.
+ */
+ struct AdditionalData
+ {
+ /**
+ * Constructor. By default, set the damping parameter to one.
+ */
+ AdditionalData (const double omega = 1);
+
+ /**
+ * Relaxation parameter.
+ */
+ double omega;
+ };
+
+ /**
+ * Constructor. Take additional flags if there are any.
+ */
+ PreconditionMultiColoredSOR (const AdditionalData &additional_data = AdditionalData());
+
+ /**
+ * This function changes the value of the additional flags.
+ */
+ void initialize (const AdditionalData &additional_data = AdditionalData());
-
- private :
- /**
- * Store a copy of the flags for this particular preconditioner.
- */
- AdditionalData additional_data;
+ };
+
+
+
+ /**
+ * A wrapper for ILU(p) for Paralution matrices.
+ *
+ * @ingroup ParalutionWrappers
+ * @ingroup Preconditioners
+ * @author Bruno Turcksin, 2014
+ */
+ template <typename Number>
+ class PreconditionILU : public PreconditionBase<Number>
+ {
+ public :
+ /**
+ * Standardized data struct to pipe additional flags to the
+ * preconditioner.
+ */
+ struct AdditionalData
+ {
+ /**
+ * Constructor. By default, set the fill-in parameter to zero.
+ */
+ AdditionalData (const unsigned int levels = 0);
+
+ /**
+ * Fill-in parameter.
+ */
+ unsigned int levels;
+ };
+
+ /**
+ * Constructor. Take additional flags if there are any.
+ */
+ PreconditionILU (const AdditionalData &additional_data = AdditionalData());
+
+ /**
+ * This function changes the value of the additional flags.
+ */
+ void initialize (const AdditionalData &additional_data = AdditionalData());
-
- private :
- /**
- * Store a copy of the flags for this particular preconditioner.
- */
- AdditionalData additional_data;
+ };
+
+
+
+ /**
+ * A wrapper for ILUT(t,m) factorization based on threshold (t) and maximum
+ * number of elements per row (m) for Paralution matrices. The preconditioner
+ * can be initialized using only the threshild value or both the threshold
+ * value and the maximum number of elements per row.
+ *
+ * @ingroup ParalutionWrappers
+ * @ingroup Preconditioners
+ * @author Bruno Turcksin, 2014
+ */
+ template <typename Number>
+ class PreconditionILUT : public PreconditionBase<Number>
+ {
+ public :
+ /**
+ * Standardized data struct to pipe additional flags to the
+ * preconditioner.
+ */
+ struct AdditionalData
+ {
+ /**
+ * Constructor. By default, set the threshold to zero.
+ */
+ AdditionalData (const Number threshold = 0.);
+
+ /**
+ * Constructor.
+ */
+ AdditionalData (const Number threshold, const unsigned int max_row);
+
+ /**
+ * Threshold.
+ */
+ Number threshold;
+
+ /**
+ * Maximum number of elements per row.
+ */
+ unsigned int max_row;
+ };
+
+ /**
+ * Constructor. Take additional flags if there are any.
+ */
+ PreconditionILUT (const AdditionalData &additional_data = AdditionalData());
+
+ /**
+ * This function changes the value of the additional flags.
+ */
+ void initialize (const AdditionalData &additional_data = AdditionalData());
- * Store a copy of the flags for this particular preconditioner.
+ };
+
+
+
+ /**
+ * A wrapper for ILU(p,q) for Paralution matrices. ILU(p,q) is based on ILU(p)
+ * with a power(q)-pattern method. Compared to the standard ILU(p), ILU(p,q)
+ * provides a higher degree of parallelism of forward and backward
+ * substitution.
+ *
+ * @ingroup ParalutionWrappers
+ * @ingroup Preconditioners
+ * @author Bruno Turcksin, 2014
+ */
+ template <typename Number>
+ class PreconditionMultiColoredILU : public PreconditionBase<Number>
+ {
+ public :
+ /**
+ * Standardized data struct to pipe additional flags to the
+ * preconditioner.
+ */
+ struct AdditionalData
+ {
+ /**
+ * Constructor. By default, set the fill-in parameter to zero and the
+ * power-pattern to one.
+ */
+ AdditionalData (const unsigned int levels = 0, const unsigned int power = 1);
+
+ /**
+ * Fill-in parameter.
+ */
+ unsigned int levels;
+
+ /**
+ * Power parameter.
+ */
+ unsigned int power;
+ };
+
+ /**
+ * Constructor. Take additional flags if there are any.
+ */
+ PreconditionMultiColoredILU (const AdditionalData &additional_data = AdditionalData());
+
+ /**
+ * This function changes the value of additional flags.
+ */
+ void initialize (const AdditionalData &additional_data = AdditionalData());
++ };
++
++ /**
++ * A simplified wrapper for AMG for Paralution matrices. To use other parameter
++ * combinations use the Paralution directly. AMG can only be built on the host.
++ *
++ * @ingroup ParalutionWrappers
++ * @ingroup Preconditioners
++ * @author Bruno Turcksin, 2014
++ */
++ template <typename Number>
++ class PreconditionAMG : public PreconditionBase<Number>
++ {
++ public :
++ /**
++ * Standardized data struct to pipe additional flags to the
++ * preconditioner.
++ */
++ struct AdditionalData
++ {
++ /**
++ * Constructor. This constructor lets Paralution choose the parameters
++ * of the algebraic multigrid.
++ */
++ AdditionalData (const unsigned int verbose = 0, const unsigned int max_iter = 2);
++
++ /**
++ * Constructor. The relaxation parameter is used only when smoothed
++ * aggregation is selected. Otherwise, this parameter is disregarded.
++ */
++ AdditionalData (const unsigned int verbose,
++ const unsigned int max_iter,
++ const mg_cycle cycle,
++ const mg_interpolation interpolation,
++ const mg_preconditioner preconditioner,
++ const mg_solver smoother,
++ const mg_solver coarse_solver,
++ const unsigned int n_unknowns_coarse_level,
++ const double relaxation_parameter = 1.2);
++
++ /**
++ * The flag is true if the default parameters of Paralution should be
++ * used.
++ */
++ bool default_parameters;
++
++ /**
++ * Cycle used by the multigrid preconditioner: V cycle, W cycle, K
++ * cycle, or F cycle.
++ */
++ mg_cycle cycle;
++
++ /**
++ * Interpolation used by the multigrid preconditioner: smoothed
++ * aggregation or aggregation.
++ */
++ mg_interpolation interpolation;
++
++ /**
++ * Preconditioner for the smoother.
++ */
++ mg_preconditioner preconditioner;
++
++ /**
++ * Smoother.
++ */
++ mg_solver smoother;
++
++ /**
++ * Solver of the coarse system.
++ */
++ mg_solver coarse_solver;
++
++ /**
++ * Verbosity: 0 = no output, 1 = print information at the
++ * start and at the end, 3 = print residual at each iteration.
++ */
++ unsigned int verbose;
++
++ /**
++ * Number of unknowns on the coarsest level.
++ */
++ unsigned int n_unknowns_coarse_level;
++
++ /**
++ * Maximum number of iterations.
++ */
++ unsigned int max_iter;
++
++ /**
++ * Relaxation parameter for smoothed interpolation aggregation.
++ */
++ double relaxation_parameter;
++ };
++
++ /**
++ * Constructor. Take additional flags if there are any.
++ */
++ PreconditionAMG (const AdditionalData &additional_data = AdditionalData());
++
++ /**
++ * Destructor.
++ */
++ ~PreconditionAMG();
++
++ /**
++ * This function changes the value of the additional flags.
++ */
++ void initialize (const AdditionalData &additional_data = AdditionalData());
+
+ private :
+ /**
- AdditionalData additional_data;
++ * Free the smoothers, the preconditioners, and the coarse level solver.
++ */
++ void clear();
++
++ /**
++ * Number of levels of the multigrid.
++ */
++ unsigned int n_levels;
++
++ /**
++ * Smoothers for each level.
++ */
++ paralution::IterativeLinearSolver<paralution::LocalMatrix<Number>, paralution::
++ LocalVector<Number>,Number> * *smoothers;
++ /**
++ * Preconditioners of the smoothers for each level.
++ */
++ paralution::Preconditioner<paralution::LocalMatrix<Number>, paralution::
++ LocalVector<Number>,Number> * *preconditioners;
++ /**
++ * Coarse grid solver.
+ */
++ paralution::IterativeLinearSolver<paralution::LocalMatrix<Number>, paralution::
++ LocalVector<Number>,Number> *coarse_solver;
+ };
+}
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif // DEAL_II_WITH_PARALUTION
+
+/*---------------------------- paralution_precondition.h ---------------------------*/
+
+#endif
+/*---------------------------- paralution_precondition.h ---------------------------*/
--- /dev/null
- const SparseMatrix<Number> &A,
+// ---------------------------------------------------------------------
+// $Id: paralution_solver.h 30040 2013-07-18 17:06:48Z maier $
+//
+// Copyright (C) 2013, 2014 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.
+//
+// ---------------------------------------------------------------------
+
+#ifndef __deal2__paralution_solver_h
+#define __deal2__paralution_solver_h
+
+
+#include <deal.II/base/config.h>
+
+#ifdef DEAL_II_WITH_PARALUTION
+
+#include <deal.II/base/std_cxx1x/shared_ptr.h>
+#include <deal.II/lac/exceptions.h>
+#include <deal.II/lac/solver_control.h>
+#include <deal.II/lac/paralution_precondition.h>
+#include <deal.II/lac/paralution_vector.h>
+#include <deal.II/lac/paralution_sparse_matrix.h>
+
+#include <paralution.hpp>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace ParalutionWrappers
+{
+ /**
+ * Base class for solver classes using the Paralution solvers.
+ *
+ * @ingroup ParalutionWrappers
+ * @author Bruno Turcksin 2013
+ */
+ class SolverBase
+ {
+ public:
+ /**
+ * Constructor. Takes the solver control object and creates the solver.
+ */
+ SolverBase (SolverControl &cn);
+
+ /**
+ * Access to object that controls convergence.
+ */
+ SolverControl &control() const;
+
+ protected:
+ /**
+ * Initialize the solver and solve the system of equations.
+ */
+ template <typename Number>
+ void execute_solve(std_cxx1x::shared_ptr<paralution::IterativeLinearSolver<paralution::
+ LocalMatrix<Number>,paralution::LocalVector<Number>,Number> > solver,
- const Vector<Number> &b,
- const PreconditionBase<Number> &preconditioner,
++ SparseMatrix<Number> &A,
+ Vector<Number> &x,
- * solver is built on the accelerator.
++ Vector<Number> &b,
++ PreconditionBase<Number> &preconditioner,
+ bool move_to_accelerator);
+
+ /**
+ * Reference to the object that controls convergence of the iteratove
+ * solver. In fact, for these Paralution wrappers, Paralution does so
+ * itself, but we copy the data from this object solition process, and
+ * copy the data back into it afterwards.
+ */
+ SolverControl &solver_control;
+ };
+
+
+
+ /**
+ * An implementation of the solver interface using the Paralution Richardson solver
+ * (paralution::fixed_point_iteration).
+ *
+ * @ingroup ParalutionWrappers
+ * @author Bruno Turcksin, 2014
+ */
+ class SolverRichardson : public SolverBase
+ {
+ public:
+ /**
+ * Standardized data struct to pipe additional data to the solver.
+ */
+ struct AdditionalData
+ {
+ /**
+ * Constructor. By default, set the verbosity of the solver to zero and
+ * the relaxation parameter to one.
+ */
+ AdditionalData (const unsigned int verbose = 0,
+ const double relaxation_parameter = 1.);
+
+ /**
+ * Verbosity of the solver: 0 = no output, 1 = print information at the
+ * start and at the end, 3 = print residual at each iteration.
+ */
+ unsigned int verbose;
+
+ /**
+ * Relaxation parameter.
+ */
+ double relaxation_parameter;
+ };
+
+ /**
+ * Constructor. AdditionalData is a structure that contains additional
+ * flags for tuning this particular solver.
+ */
+ SolverRichardson (SolverControl &cn,
+ const AdditionalData &data = AdditionalData());
+
+ /**
+ * Solve the linear system <tt>Ax=b</tt> using the CG solver of
+ * Paralution. If the flag @p move_to_accelerator is set to true, the
- void solve (const SparseMatrix<Number> &A,
- Vector<Number> &x,
- const Vector<Number> &b,
- const PreconditionBase<Number> &preconditioner,
- bool move_to_accelerator=false);
++ * solver, the right-hand side, the matrix, and the solution vector are
++ * moved on the accelerator once the solver is built.
+ */
+ template <typename Number>
- * solver is built on the accelerator.
++ void solve (SparseMatrix<Number> &A,
++ Vector<Number> &x,
++ Vector<Number> &b,
++ PreconditionBase<Number> &preconditioner,
++ bool move_to_accelerator=false);
+
+ private:
+ /**
+ * store a copy of the flags for this solver.
+ */
+ const AdditionalData additional_data;
+ };
+
+
+
+ /**
+ * An implementation of the solver interface using the Paralution CG solver.
+ *
+ * @ingroup ParalutionWrappers
+ * @author Bruno Turcksin, 2013
+ */
+ class SolverCG : public SolverBase
+ {
+ public:
+ /**
+ * Standardized data struct to pipe additional data to the solver.
+ */
+ struct AdditionalData
+ {
+ /**
+ * Constructor. By default, set the verbosity of the solver to zero.
+ */
+ AdditionalData (const unsigned int verbose = 0);
+
+ /**
+ * Verbosity of the solver: 0 = no output, 1 = print information at the
+ * start and at the end, 3 = print residual at each iteration.
+ */
+ unsigned int verbose;
+
+ /**
+ * Relaxation parameter.
+ */
+ double relaxation_parameter;
+ };
+
+ /**
+ * Constructor. AdditionalData is a structure that contains additional
+ * flags for tuning this particular solver.
+ */
+ SolverCG (SolverControl &cn,
+ const AdditionalData &data = AdditionalData());
+
+ /**
+ * Solve the linear system <tt>Ax=b</tt> using the CG solver of
+ * Paralution. If the flag @p move_to_accelerator is set to true, the
- void solve (const SparseMatrix<Number> &A,
- Vector<Number> &x,
- const Vector<Number> &b,
- const PreconditionBase<Number> &preconditioner,
- bool move_to_accelerator=false);
++ * solver, the right-hand side, the matrix, and the solution vector are
++ * moved on the accelerator once the solver is built.
+ */
+ template <typename Number>
- * solver is built on the accelerator.
++ void solve (SparseMatrix<Number> &A,
++ Vector<Number> &x,
++ Vector<Number> &b,
++ PreconditionBase<Number> &preconditioner,
++ bool move_to_accelerator=false);
+
+ private:
+ /**
+ * store a copy of the flags for this solver.
+ */
+ const AdditionalData additional_data;
+ };
+
+
+
+ /**
+ * An implementation of the solver interface using the Paralution CR solver.
+ *
+ * @ingroup ParalutionWrappers
+ * @author Bruno Turcksin, 2014
+ */
+ class SolverCR : public SolverBase
+ {
+ public:
+ /**
+ * Standardized data struct to pipe additional data to the solver.
+ */
+ struct AdditionalData
+ {
+ /**
+ * Constructor. By default, set the verbosity of the solver to zero.
+ */
+ AdditionalData (const unsigned int verbose = 0);
+
+ /**
+ * Verbosity of the solver: 0 = no output, 1 = print information at the
+ * start and at the end, 3 = print residual at each iteration.
+ */
+ unsigned int verbose;
+ };
+
+ /**
+ * Constructor. AdditionalData is a structure that contains additional
+ * flags for tuning this particular solver.
+ */
+ SolverCR (SolverControl &cn,
+ const AdditionalData &data = AdditionalData());
+
+ /**
+ * Solve the linear system <tt>Ax=b</tt> using the CR solver of
+ * Paralution. If the flag @p move_to_accelerator is set to true, the
- void solve (const SparseMatrix<Number> &A,
- Vector<Number> &x,
- const Vector<Number> &b,
- const PreconditionBase<Number> &preconditioner,
- bool move_to_accelerator=false);
++ * solver, the right-hand side, the matrix, and the solution vector are
++ * moved on the accelerator once the solver is built.
+ */
+ template <typename Number>
- * solver is built on the accelerator.
++ void solve (SparseMatrix<Number> &A,
++ Vector<Number> &x,
++ Vector<Number> &b,
++ PreconditionBase<Number> &preconditioner,
++ bool move_to_accelerator=false);
+
+ private:
+ /**
+ * store a copy of the flags for this solver.
+ */
+ const AdditionalData additional_data;
+ };
+
+
+
+ /**
+ * An implementation of the solver interface using the Paralution DPCG solver.
+ * DPCG stands for Deflated Preconditioned Conjugate Gradient. It is a
+ * two-level preconditioned CG algorithm to iteratively solve and
+ * ill-conditioned linear system. Deflation tries to remove the bad
+ * eigenvalues from the spectrum of the preconditioned system.
+ *
+ * @ingroup ParalutionWrappers
+ * @author Bruno Turcksin, 2014
+ */
+ class SolverDPCG : public SolverBase
+ {
+ public:
+ /**
+ * Standardized data struct to pipe additional data to the solver.
+ */
+ struct AdditionalData
+ {
+ /**
+ * Constructor. By default, set the verbosity of the solver to zero and
+ * set the number of deflated vectors to two.
+ */
+ AdditionalData (const unsigned int verbose = 0,
+ const unsigned int n_deflated_vectors = 2);
+
+ /**
+ * Verbosity of the solver: 0 = no output, 1 = print information at the
+ * start and at the end, 3 = print residual at each iteration.
+ */
+ unsigned int verbose;
+
+ /**
+ * Number of deflated vectors.
+ */
+ unsigned int n_deflated_vectors;
+ };
+
+ /**
+ * Constructor. AdditionalData is a structure that contains additional
+ * flags for tuning this particular solver.
+ */
+ SolverDPCG (SolverControl &cn,
+ const AdditionalData &data = AdditionalData());
+
+ /**
+ * Solve the linear system <tt>Ax=b</tt> using the DPCG solver of
+ * Paralution. If the flag @p move_to_accelerator is set to true, the
- void solve (const SparseMatrix<Number> &A,
- Vector<Number> &x,
- const Vector<Number> &b,
- const PreconditionBase<Number> &preconditioner,
- bool move_to_accelerator=false);
++ * solver, the right-hand side, the matrix, and the solution vector are
++ * moved on the accelerator once the solver is built.
+ */
+ template <typename Number>
- * solver is built on the accelerator.
++ void solve (SparseMatrix<Number> &A,
++ Vector<Number> &x,
++ Vector<Number> &b,
++ PreconditionBase<Number> &preconditioner,
++ bool move_to_accelerator=false);
+
+ private:
+ /**
+ * store a copy of the flags for this solver.
+ */
+ const AdditionalData additional_data;
+ };
+
+
+
+ /**
+ * An implementation of the solver interface using the Paralution BiCGStab
+ * solver.
+ *
+ * @ingroup ParalutionWrappers
+ * @author Bruno Turcksin, 2013
+ */
+ class SolverBicgstab : public SolverBase
+ {
+ public:
+ /**
+ * Standardized data struct to pipe additional data to the solver.
+ */
+ struct AdditionalData
+ {
+ /**
+ * Constructor. By default, set the verbosity of the solver to zero.
+ */
+ AdditionalData (const unsigned int verbose = 0);
+
+ /**
+ * Verbosity of the solver: 0 = no output, 1 = print information at the
+ * start and at the end, 3 = print residual at each iteration.
+ */
+ unsigned int verbose;
+ };
+
+ /**
+ * Constructor. AdditionalData is a structure that contains additional
+ * flags for tuning this particular solver.
+ */
+ SolverBicgstab (SolverControl &cn,
+ const AdditionalData &data = AdditionalData());
+
+ /**
+ * Solve the linear system <tt>Ax=b</tt> using the BiCGStab solver of
+ * Paralution. If the flag @p move_to_accelerator is set to true, the
- void solve (const SparseMatrix<Number> &A,
- Vector<Number> &x,
- const Vector<Number> &b,
- const PreconditionBase<Number> &preconditioner,
- bool move_to_accelerator=false);
++ * solver, the right-hand side, the matrix, and the solution vector are
++ * moved on the accelerator once the solver is built.
+ */
+ template <typename Number>
- * solver is built on the accelerator.
++ void solve (SparseMatrix<Number> &A,
++ Vector<Number> &x,
++ Vector<Number> &b,
++ PreconditionBase<Number> &preconditioner,
++ bool move_to_accelerator=false);
+
+ private:
+ /**
+ * store a copy of the flags for this solver.
+ */
+ const AdditionalData additional_data;
+ };
+
+
+
+ /**
+ * An implementation of the solver interface using the Paralution GMRES
+ * solver.
+ *
+ * @ingroup ParalutionWrappers
+ * @author Bruno Turcksin, 2013
+ */
+ class SolverGMRES : public SolverBase
+ {
+ public:
+ /**
+ * Standardized data struct to pipe additional data to the solver.
+ */
+ struct AdditionalData
+ {
+ /**
+ * Constructor. By default, set the verbosity of the solver to zero and
+ * set the size of the Krylov space to 30, i.e. do a restart every 30
+ * iterations.
+ */
+ AdditionalData (const unsigned int verbose = 0,
+ const unsigned int restart_parameter = 30);
+
+ /**
+ * Verbosity of the solver: 0 = no output, 1 = print information at the
+ * start and at the end, 3 = print residual at each iteration.
+ */
+ unsigned int verbose;
+
+ /**
+ * Size of the Krylov space.
+ */
+ unsigned int restart_parameter;
+ };
+
+ /**
+ * Constructor. AdditionalData is a structure that contains additional
+ * flags for tuning this particular solver.
+ */
+ SolverGMRES (SolverControl &cn,
+ const AdditionalData &data = AdditionalData());
+
+ /**
+ * Solve the linear system <tt>Ax=b</tt> using the GMRES solver of
+ * Paralution. If the flag @p move_to_accelerator is set to true, the
- void solve (const SparseMatrix<Number> &A,
- Vector<Number> &x,
- const Vector<Number> &b,
- const PreconditionBase<Number> &preconditioner,
- bool move_to_accelerator=false);
++ * solver, the right-hand side, the matrix, and the solution vector are
++ * moved on the accelerator once the solver is built.
+ */
+ template <typename Number>
- * solver is built on the accelerator.
++ void solve (SparseMatrix<Number> &A,
++ Vector<Number> &x,
++ Vector<Number> &b,
++ PreconditionBase<Number> &preconditioner,
++ bool move_to_accelerator=false);
+
+ private:
+ /**
+ * store a copy of the flags for this solver.
+ */
+ const AdditionalData additional_data;
+ };
+
+
+
+ /**
+ * An implementation of the solver interface using the Paralution GMRES
+ * solver.
+ *
+ * @ingroup ParalutionWrappers
+ * @author Bruno Turcksin, 2014
+ */
+ class SolverFGMRES : public SolverBase
+ {
+ public:
+ /**
+ * Standardized data struct to pipe additional data to the solver.
+ */
+ struct AdditionalData
+ {
+ /**
+ * Constructor. By default, set the verbosity of the solver to zero and
+ * set the size of the Krylov space to 30, i.e. do a restart every 30
+ * iterations.
+ */
+ AdditionalData (const unsigned int verbose = 0,
+ const unsigned int restart_parameter = 30);
+
+ /**
+ * Verbosity of the solver: 0 = no output, 1 = print information at the
+ * start and at the end, 3 = print residual at each iteration.
+ */
+ unsigned int verbose;
+
+ /**
+ * Size of the Krylov space.
+ */
+ unsigned int restart_parameter;
+ };
+
+ /**
+ * Constructor. AdditionalData is a structure that contains additional
+ * flags for tuning this particular solver.
+ */
+ SolverFGMRES (SolverControl &cn,
+ const AdditionalData &data = AdditionalData());
+
+ /**
+ * Solve the linear system <tt>Ax=b</tt> using the GMRES solver of
+ * Paralution. If the flag @p move_to_accelerator is set to true, the
- void solve (const SparseMatrix<Number> &A,
- Vector<Number> &x,
- const Vector<Number> &b,
- const PreconditionBase<Number> &preconditioner,
- bool move_to_accelerator=false);
++ * solver, the right-hand side, the matrix, and the solution vector are
++ * moved on the accelerator once the solver is built.
+ */
+ template <typename Number>
++ void solve (SparseMatrix<Number> &A,
++ Vector<Number> &x,
++ Vector<Number> &b,
++ PreconditionBase<Number> &preconditioner,
++ bool move_to_accelerator=false);
+
+ private:
+ /**
+ * store a copy of the flags for this solver.
+ */
+ const AdditionalData additional_data;
+ };
+}
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif // DEAL_II_WITH_PARALUTION
+
+/*---------------------------- paralution_solver.h ---------------------------*/
+
+#endif
+/*---------------------------- paralution_solver.h ---------------------------*/
--- /dev/null
- :
- additional_data(additional_data)
+// ---------------------------------------------------------------------
+// $Id: paralution_precondition.cc 31349 2013-10-20 19:07:06Z maier $
+//
+// Copyright (C) 2014 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.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/lac/paralution_precondition.h>
+
+#ifdef DEAL_II_WITH_PARALUTION
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace ParalutionWrappers
+{
+ template <typename Number>
+ PreconditionBase<Number>::PreconditionBase()
+ {}
+
+
+
+ /* -------------------------- PreconditionJacobi -------------------------- */
+
+ template <typename Number>
+ PreconditionJacobi<Number>::PreconditionJacobi()
+ {
+ this->preconditioner.reset(new paralution::Jacobi<paralution::LocalMatrix<Number>,
+ paralution::LocalVector<Number>,Number>);
+ }
+
+
+
+ /* -------------------------- PreconditionSGS ----------------------------- */
+
+ template <typename Number>
+ PreconditionSGS<Number>::PreconditionSGS()
+ {
+ this->preconditioner.reset(new paralution::SGS<paralution::LocalMatrix<Number>,
+ paralution::LocalVector<Number>,Number>);
+ }
+
+
+
+ /* -------------------------- PreconditionMultiColoredSGS ----------------- */
+
+ template <typename Number>
+ PreconditionMultiColoredSGS<Number>::PreconditionMultiColoredSGS()
+ {
+ this->preconditioner.reset(new paralution::MultiColoredSGS<paralution::LocalMatrix<Number>,
+ paralution::LocalVector<Number>,Number>);
+ }
+
+
+
+ /* -------------------------- PreconditionMulticoloredSOR ----------------- */
+
+ template <typename Number>
+ PreconditionMultiColoredSOR<Number>::AdditionalData::
+ AdditionalData(const double omega)
+ :
+ omega(omega)
+ {}
+
+
+
+ template <typename Number>
+ PreconditionMultiColoredSOR<Number>::PreconditionMultiColoredSOR(
+ const AdditionalData &additional_data)
- :
- additional_data(additional_data)
+ {
+ this->preconditioner.reset(new paralution::MultiColoredGS<paralution::LocalMatrix<Number>,
+ paralution::LocalVector<Number>,Number>);
+
+ initialize(additional_data);
+ }
+
+
+
+ template <typename Number>
+ void PreconditionMultiColoredSOR<Number>::initialize(const AdditionalData &additional_data)
+ {
+ // Downcast the preconditioner pointer to use SetRelaxation
+ paralution::MultiColoredGS<paralution::LocalMatrix<Number>,paralution::
+ LocalVector<Number>,Number>* downcasted_ptr = static_cast<paralution::
+ MultiColoredGS<paralution::LocalMatrix<Number>,paralution::LocalVector<Number>,
+ Number>* >(this->preconditioner.get());
+
+ downcasted_ptr->SetRelaxation(additional_data.omega);
+ }
+
+
+
+ /* -------------------------- PreconditionILU ----------------------------- */
+
+ template <typename Number>
+ PreconditionILU<Number>::AdditionalData::AdditionalData(const unsigned int levels)
+ :
+ levels(levels)
+ {}
+
+
+
+ template <typename Number>
+ PreconditionILU<Number>::PreconditionILU(const AdditionalData &additional_data)
- :
- additional_data(additional_data)
+ {
+ this->preconditioner.reset(new paralution::ILU<paralution::LocalMatrix<Number>,
+ paralution::LocalVector<Number>,Number>);
+
+ initialize(additional_data);
+ }
+
+
+
+ template <typename Number>
+ void PreconditionILU<Number>::initialize(const AdditionalData &additional_data)
+ {
+ // Downcast the preconditioner pointer to use Set
+ paralution::ILU<paralution::LocalMatrix<Number>,paralution::LocalVector<Number>,
+ Number>* downcasted_ptr = static_cast<paralution::ILU<paralution::
+ LocalMatrix<Number>,paralution::LocalVector<Number>,Number>* >(this->preconditioner.get());
+
+ downcasted_ptr->Set(additional_data.levels);
+ }
+
+
+
+ /* -------------------------- PreconditionILUT ---------------------------- */
+
+ template <typename Number>
+ PreconditionILUT<Number>::AdditionalData::AdditionalData(const Number threshold)
+ :
+ threshold(threshold),
+ max_row(0)
+ {}
+
+
+
+ template <typename Number>
+ PreconditionILUT<Number>::AdditionalData::AdditionalData(const Number threshold,
+ const unsigned int max_row)
+ :
+ threshold(threshold),
+ max_row(max_row)
+ {}
+
+
+
+ template <typename Number>
+ PreconditionILUT<Number>::PreconditionILUT(const AdditionalData &additional_data)
- :
- additional_data(additional_data)
+ {
+ this->preconditioner.reset(new paralution::ILUT<paralution::LocalMatrix<Number>,
+ paralution::LocalVector<Number>,Number>);
+
+ initialize(additional_data);
+ }
+
+
+
+ template <typename Number>
+ void PreconditionILUT<Number>::initialize(const AdditionalData &additional_data)
+ {
+ // Downcast the preconditioner pointer to use Set
+ paralution::ILUT<paralution::LocalMatrix<Number>,paralution::LocalVector<Number>,
+ Number>* downcasted_ptr = static_cast<paralution::ILUT<paralution::
+ LocalMatrix<Number>,paralution::LocalVector<Number>,Number>* >(this->preconditioner.get());
+
+ if (additional_data.max_row==0)
+ downcasted_ptr->Set(additional_data.threshold);
+ else
+ downcasted_ptr->Set(additional_data.threshold,additional_data.max_row);
+ }
+
+
+
+ /* -------------------------- PreconditionMulticoloredILU --------------- */
+
+ template <typename Number>
+ PreconditionMultiColoredILU<Number>::AdditionalData::
+ AdditionalData(const unsigned int levels, const unsigned int power)
+ :
+ levels(levels),
+ power(power)
+ {}
+
+
+
+ template <typename Number>
+ PreconditionMultiColoredILU<Number>::
+ PreconditionMultiColoredILU(const AdditionalData &additional_data)
+ {
+ this->preconditioner.reset(new paralution::MultiColoredILU<paralution::LocalMatrix<Number>,
+ paralution::LocalVector<Number>,Number>);
+
+ initialize(additional_data);
+ }
+
+
+
+ template <typename Number>
+ void PreconditionMultiColoredILU<Number>::initialize(const AdditionalData &additional_data)
+ {
+ // Downcast the preconditioner pointer to use Set
+ paralution::MultiColoredILU<paralution::LocalMatrix<Number>,paralution::LocalVector<Number>,
+ Number>* downcasted_ptr = static_cast<paralution::MultiColoredILU<paralution::
+ LocalMatrix<Number>,paralution::LocalVector<Number>,Number>* >(this->preconditioner.get());
+
+ downcasted_ptr->Set(additional_data.levels,additional_data.power);
+ }
++
++
++
++ /* -------------------------- PreconditionAMG --------------------------- */
++
++ template <typename Number>
++ PreconditionAMG<Number>::AdditionalData::AdditionalData(const unsigned int verbose,
++ const unsigned int max_iter)
++ :
++ default_parameters(true),
++ verbose(verbose),
++ max_iter(max_iter)
++ {}
++
++
++
++ template <typename Number>
++ PreconditionAMG<Number>::AdditionalData::AdditionalData(const unsigned int verbose,
++ const unsigned int max_iter,
++ const mg_cycle cycle,
++ const mg_interpolation interpolation,
++ const mg_preconditioner preconditioner,
++ const mg_solver smoother,
++ const mg_solver coarse_solver,
++ const unsigned int n_unknowns_coarse_level,
++ const double relaxation_parameter)
++ :
++ default_parameters(false),
++ cycle(cycle),
++ interpolation(interpolation),
++ preconditioner(preconditioner),
++ smoother(smoother),
++ coarse_solver(coarse_solver),
++ verbose(verbose),
++ n_unknowns_coarse_level(n_unknowns_coarse_level),
++ max_iter(max_iter),
++ relaxation_parameter(relaxation_parameter)
++ {}
++
++
++
++ template <typename Number>
++ PreconditionAMG<Number>::PreconditionAMG(const AdditionalData &additional_data)
++ :
++ n_levels(1),
++ smoothers(NULL),
++ preconditioners(NULL),
++ coarse_solver(NULL)
++ {
++ this->preconditioner.reset(new paralution::AMG<paralution::LocalMatrix<Number>,
++ paralution::LocalVector<Number>,Number>);
++
++ initialize(additional_data);
++ }
++
++
++
++ template <typename Number>
++ PreconditionAMG<Number>::~PreconditionAMG()
++ {
++ clear();
++ }
++
++
++
++ template <typename Number>
++ void PreconditionAMG<Number>::initialize(const AdditionalData &additional_data)
++ {
++ // Downcast the preconditioner pointer
++ paralution::AMG<paralution::LocalMatrix<Number>,paralution::LocalVector<Number>,
++ Number>* downcasted_ptr = static_cast<paralution::AMG<paralution::
++ LocalMatrix<Number>,paralution::LocalVector<Number>,Number>* >(this->preconditioner.get());
++
++ // Set the maximum number of iterations.
++ downcasted_ptr->InitMaxIter(additional_data.max_iter);
++
++ // Set the verbosity
++ downcasted_ptr->Verbose(additional_data.verbose);
++
++ if (additional_data.default_parameters==false)
++ {
++ // Free the smoothers, the preconditioners, and the coarse level solver,
++ // if necessary.
++ clear();
++
++ // Set the number of unknowns on the coarsest level.
++ downcasted_ptr->SetCoarsestLevel(additional_data.n_unknowns_coarse_level);
++
++ // Set the interpolation type.
++ switch (additional_data.interpolation)
++ {
++ case aggregation :
++ {
++ downcasted_ptr->SetInterpolation(paralution::Aggregation);
++ break;
++ }
++ case smoothed_aggregation :
++ {
++ downcasted_ptr->SetInterpolation(paralution::SmoothedAggregation);
++ downcasted_ptr->SetInterpRelax(additional_data.relaxation_parameter);
++ break;
++ }
++ default :
++ {
++ AssertThrow(false,ExcMessage("Unknown interpolation type for PreconditionAMG."));
++ }
++ }
++
++ // Set the type of cycle
++ switch (additional_data.cycle)
++ {
++ case V_cycle :
++ {
++ downcasted_ptr->SetCycle(paralution::Vcycle);
++ break;
++ }
++ case W_cycle :
++ {
++ downcasted_ptr->SetCycle(paralution::Wcycle);
++ break;
++ }
++ case K_cycle :
++ {
++ downcasted_ptr->SetCycle(paralution::Kcycle);
++ break;
++ }
++ case F_cycle :
++ {
++ downcasted_ptr->SetCycle(paralution::Fcycle);
++ }
++ default :
++ {
++ AssertThrow(false,ExcMessage("Unknown cycle type for PreconditionAMG."));
++ }
++ }
++
++ // Use manual smoothers.
++ downcasted_ptr->SetManualSmoothers(true);
++
++ // Use manual coarse grid solver.
++ downcasted_ptr->SetManualSolver(true);
++
++ // Build the hierarchy of the grids.
++ downcasted_ptr->BuildHierarchy();
++
++ // Get the number of grids/levels.
++ n_levels = downcasted_ptr->GetNumLevels();
++
++ // Smoothers for each level.
++ smoothers = new paralution::IterativeLinearSolver<paralution::LocalMatrix<Number>,
++ paralution::LocalVector<Number>, Number>* [n_levels-1];
++ switch (additional_data.smoother)
++ {
++ case richardson :
++ {
++ for (unsigned int i=0; i<n_levels-1; ++i)
++ {
++ paralution::FixedPoint<paralution::LocalMatrix<Number>, paralution::LocalVector<Number>,
++ Number> *sm = new paralution::FixedPoint<paralution::LocalMatrix<Number>,
++ paralution::LocalVector<Number>, Number>;
++ smoothers[i] = sm;
++ }
++ break;
++ }
++ case cg :
++ {
++ for (unsigned int i=0; i<n_levels-1; ++i)
++ {
++ paralution::CG<paralution::LocalMatrix<Number>, paralution::LocalVector<Number>,
++ Number> *sm = new paralution::CG<paralution::LocalMatrix<Number>,
++ paralution::LocalVector<Number>, Number>;
++ smoothers[i] = sm;
++ }
++ break;
++ }
++ case cr :
++ {
++ for (unsigned int i=0; i<n_levels-1; ++i)
++ {
++ paralution::CR<paralution::LocalMatrix<Number>, paralution::LocalVector<Number>,
++ Number> *sm = new paralution::CR<paralution::LocalMatrix<Number>,
++ paralution::LocalVector<Number>, Number>;
++ smoothers[i] = sm;
++ }
++ break;
++ }
++ case bicgstab :
++ {
++ for (unsigned int i=0; i<n_levels-1; ++i)
++ {
++ paralution::BiCGStab<paralution::LocalMatrix<Number>, paralution::LocalVector<Number>,
++ Number> *sm = new paralution::BiCGStab<paralution::LocalMatrix<Number>,
++ paralution::LocalVector<Number>, Number>;
++ smoothers[i] = sm;
++ }
++ break;
++ }
++ case gmres :
++ {
++ for (unsigned int i=0; i<n_levels-1; ++i)
++ {
++ paralution::GMRES<paralution::LocalMatrix<Number>, paralution::LocalVector<Number>,
++ Number> *sm = new paralution::GMRES<paralution::LocalMatrix<Number>,
++ paralution::LocalVector<Number>, Number>;
++ smoothers[i] = sm;
++ }
++ break;
++ }
++ case fgmres :
++ {
++ for (unsigned int i=0; i<n_levels-1; ++i)
++ {
++ paralution::FGMRES<paralution::LocalMatrix<Number>, paralution::LocalVector<Number>,
++ Number> *sm = new paralution::FGMRES<paralution::LocalMatrix<Number>,
++ paralution::LocalVector<Number>, Number>;
++ smoothers[i] = sm;
++ }
++ break;
++ }
++ default :
++ {
++ AssertThrow(false,ExcMessage("Unknown smoother for PreconditionAMG."));
++ }
++ }
++
++ // Preconditioners for each level.
++ preconditioners = new paralution::Preconditioner<paralution::LocalMatrix<Number>,
++ paralution::LocalVector<Number>, Number>*[n_levels-1];
++ switch (additional_data.preconditioner)
++ {
++ case jacobi :
++ {
++ for (unsigned int i=0; i<n_levels-1; ++i)
++ {
++ paralution::Jacobi<paralution::LocalMatrix<Number>, paralution::LocalVector<Number>,
++ Number> *precond = new paralution::Jacobi<paralution::LocalMatrix<Number>,
++ paralution::LocalVector<Number>, Number>;
++ preconditioners[i] = precond;
++ }
++ break;
++ }
++ case sgs :
++ {
++ for (unsigned int i=0; i<n_levels-1; ++i)
++ {
++ paralution::SGS<paralution::LocalMatrix<Number>, paralution::LocalVector<Number>,
++ Number> *precond = new paralution::SGS<paralution::LocalMatrix<Number>,
++ paralution::LocalVector<Number>, Number>;
++ preconditioners[i] = precond;
++ }
++ break;
++ }
++ case multicolored_sgs :
++ {
++ for (unsigned int i=0; i<n_levels-1; ++i)
++ {
++ paralution::MultiColoredSGS<paralution::LocalMatrix<Number>,
++ paralution::LocalVector<Number>, Number> *precond =
++ new paralution::MultiColoredSGS<paralution::LocalMatrix<Number>,
++ paralution::LocalVector<Number>, Number>;
++ preconditioners[i] = precond;
++ }
++ break;
++ }
++ case multicolored_sor :
++ {
++ for (unsigned int i=0; i<n_levels-1; ++i)
++ {
++ paralution::MultiColoredGS<paralution::LocalMatrix<Number>,
++ paralution::LocalVector<Number>, Number> *precond =
++ new paralution::MultiColoredGS<paralution::LocalMatrix<Number>,
++ paralution::LocalVector<Number>, Number>;
++ preconditioners[i] = precond;
++ }
++ break;
++ }
++ case ilu :
++ {
++ for (unsigned int i=0; i<n_levels-1; ++i)
++ {
++ paralution::ILU<paralution::LocalMatrix<Number>, paralution::LocalVector<Number>,
++ Number> *precond = new paralution::ILU<paralution::LocalMatrix<Number>,
++ paralution::LocalVector<Number>, Number>;
++ preconditioners[i] = precond;
++ }
++ break;
++ }
++ case ilut :
++ {
++ for (unsigned int i=0; i<n_levels-1; ++i)
++ {
++ paralution::ILUT<paralution::LocalMatrix<Number>, paralution::LocalVector<Number>,
++ Number> *precond = new paralution::ILUT<paralution::LocalMatrix<Number>,
++ paralution::LocalVector<Number>, Number>;
++ preconditioners[i] = precond;
++ }
++ break;
++ }
++ case multicolored_ilu :
++ {
++ for (unsigned int i=0; i<n_levels-1; ++i)
++ {
++ paralution::MultiColoredILU<paralution::LocalMatrix<Number>,
++ paralution::LocalVector<Number>, Number> *precond =
++ new paralution::MultiColoredILU<paralution::LocalMatrix<Number>,
++ paralution::LocalVector<Number>, Number>;
++ preconditioners[i] = precond;
++ }
++ break;
++ }
++ default :
++ {
++ AssertThrow(false,ExcMessage("Unknown preconditioner for the smoother for PreconditionAMG."));
++ }
++ }
++
++ for (unsigned int i=0; i<n_levels-1; ++i)
++ {
++ smoothers[i]->SetPreconditioner(*preconditioners[i]);
++ smoothers[i]->Verbose(0);
++ }
++ downcasted_ptr->SetSmoother(smoothers);
++
++ // Coarse solver
++ switch (additional_data.coarse_solver)
++ {
++ case richardson :
++ {
++ paralution::FixedPoint<paralution::LocalMatrix<Number>, paralution::LocalVector<Number>,
++ Number> *cs = new paralution::FixedPoint<paralution::LocalMatrix<Number>,
++ paralution::LocalVector<Number>, Number>;
++ coarse_solver = cs;
++ break;
++ }
++ case cg :
++ {
++ paralution::CG<paralution::LocalMatrix<Number>, paralution::LocalVector<Number>,
++ Number> *cs = new paralution::CG<paralution::LocalMatrix<Number>,
++ paralution::LocalVector<Number>, Number>;
++ coarse_solver = cs;
++ break;
++ }
++ case cr :
++ {
++ paralution::CR<paralution::LocalMatrix<Number>, paralution::LocalVector<Number>,
++ Number> *cs = new paralution::CR<paralution::LocalMatrix<Number>,
++ paralution::LocalVector<Number>, Number>;
++ coarse_solver = cs;
++ break;
++ }
++ case bicgstab :
++ {
++ paralution::BiCGStab<paralution::LocalMatrix<Number>, paralution::LocalVector<Number>,
++ Number> *cs = new paralution::BiCGStab<paralution::LocalMatrix<Number>,
++ paralution::LocalVector<Number>, Number>;
++ coarse_solver = cs;
++ break;
++ }
++ case gmres :
++ {
++ paralution::GMRES<paralution::LocalMatrix<Number>, paralution::LocalVector<Number>,
++ Number> *cs = new paralution::GMRES<paralution::LocalMatrix<Number>,
++ paralution::LocalVector<Number>, Number>;
++ coarse_solver = cs;
++ break;
++ }
++ case fgmres :
++ {
++ paralution::FGMRES<paralution::LocalMatrix<Number>, paralution::LocalVector<Number>,
++ Number> *cs = new paralution::FGMRES<paralution::LocalMatrix<Number>,
++ paralution::LocalVector<Number>, Number>;
++ coarse_solver = cs;
++ break;
++ }
++ default :
++ {
++ AssertThrow(false,ExcMessage("Unknown smoother for PreconditionAMG."));
++ }
++ }
++
++ downcasted_ptr->SetSolver(*coarse_solver);
++ }
++ }
++
++
++
++ template <typename Number>
++ void PreconditionAMG<Number>::clear()
++ {
++ // Free the coarse solver.
++ if (coarse_solver!=NULL)
++ {
++ delete coarse_solver;
++ coarse_solver = NULL;
++ }
++
++ // Free the preconditioners
++ for (unsigned int i=0; i<n_levels-1; ++i)
++ {
++ if (preconditioners[i]!=NULL)
++ {
++ delete preconditioners[i];
++ preconditioners[i] = NULL;
++ }
++ }
++ if (preconditioners!=NULL)
++ {
++ delete [] preconditioners;
++ preconditioners = NULL;
++ }
++
++ // Free the smoothers.
++ for (unsigned int i=0; i<n_levels-1; ++i)
++ {
++ if (smoothers[i]!=NULL)
++ {
++ delete smoothers[i];
++ smoothers[i] = NULL;
++ }
++ }
++ if (smoothers!=NULL)
++ {
++ delete [] smoothers;
++ smoothers = NULL;
++ }
++ }
+}
+
+// Explicit instantiations
+namespace ParalutionWrappers
+{
+ template class PreconditionBase<float>;
+ template class PreconditionBase<double>;
+ template class PreconditionJacobi<float>;
+ template class PreconditionJacobi<double>;
+ template class PreconditionSGS<float>;
+ template class PreconditionSGS<double>;
+ template class PreconditionMultiColoredSGS<float>;
+ template class PreconditionMultiColoredSGS<double>;
+ template class PreconditionMultiColoredSOR<float>;
+ template class PreconditionMultiColoredSOR<double>;
+ template class PreconditionILU<float>;
+ template class PreconditionILU<double>;
+ template class PreconditionILUT<float>;
+ template class PreconditionILUT<double>;
+ template class PreconditionMultiColoredILU<float>;
+ template class PreconditionMultiColoredILU<double>;
++ template class PreconditionAMG<float>;
++ template class PreconditionAMG<double>;
+}
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif // DEAL_II_WITH_PARALUTION
--- /dev/null
- const SparseMatrix<Number> &A,
+// ---------------------------------------------------------------------
+// $Id: paralution_solver.cc 31349 2013-10-20 19:07:06Z maier $
+//
+// Copyright (C) 2013, 2014 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.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/lac/paralution_solver.h>
+
+#ifdef DEAL_II_WITH_PARALUTION
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace ParalutionWrappers
+{
+ SolverBase::SolverBase (SolverControl &cn)
+ :
+ solver_control (cn)
+ {}
+
+
+
+ SolverControl &SolverBase::control() const
+ {
+ return solver_control;
+ }
+
+
+
+ template <typename Number>
+ void SolverBase::execute_solve(std_cxx1x::shared_ptr<paralution::IterativeLinearSolver<paralution::
+ LocalMatrix<Number>,paralution::LocalVector<Number>,Number> > solver,
- const Vector<Number> &b,
- const PreconditionBase<Number> &preconditioner,
++ SparseMatrix<Number> &A,
+ Vector<Number> &x,
- // Set the system to solve
++ Vector<Number> &b,
++ PreconditionBase<Number> &preconditioner,
+ bool move_to_accelerator)
+ {
+ // Set the preconditioner.
+ solver->SetPreconditioner(*(preconditioner.preconditioner));
+
- solver->MoveToAccelerator();
++ // Set the system to solve.
+ solver->SetOperator(A.paralution_matrix());
+
+ // Set absolute tolerance, relative tolerance, divergence tolerance,
+ // maximum number of iterations.
+ solver->Init(solver_control.tolerance(),0.,1.e100,
+ solver_control.max_steps());
+
++ // Build the solver.
++ solver->Build();
++
+ // Move the solver to the accelerator if necessary.
+ if (move_to_accelerator==true)
- solver->Build();
++ {
++ A.move_to_accelerator();
++ x.move_to_accelerator();
++ b.move_to_accelerator();
++ solver->MoveToAccelerator();
++ }
+
- void SolverRichardson::solve(const SparseMatrix<Number> &A,
- Vector<Number> &x,
- const Vector<Number> &b,
- const PreconditionBase<Number> &preconditioner,
- bool move_to_accelerator)
+ solver->Solve(b.paralution_vector(),&(x.paralution_vector()));
++
++ // 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 only exists on the host so the solver
++ // needs to be on the most.
++ solver->MoveToHost();
++
++ solver_control.check (solver->GetIterationCount(), solver->GetCurrentResidual());
++
++ if (solver_control.last_check() != SolverControl::success)
++ AssertThrow(false, SolverControl::NoConvergence (solver_control.last_step(),
++ solver_control.last_value()));
+ }
+
+
+
+ /* --------------------- SolverRichardson----------------- */
+
+ SolverRichardson::AdditionalData::AdditionalData(const unsigned int verbose,
+ const double relaxation_parameter)
+ :
+ verbose(verbose),
+ relaxation_parameter(relaxation_parameter)
+ {}
+
+
+
+ SolverRichardson::SolverRichardson (SolverControl &cn,
+ const AdditionalData &data)
+ :
+ SolverBase (cn),
+ additional_data (data)
+ {}
+
+
+
+ template <typename Number>
- void SolverCG::solve(const SparseMatrix<Number> &A,
- Vector<Number> &x,
- const Vector<Number> &b,
- const PreconditionBase<Number> &preconditioner,
- bool move_to_accelerator)
++ void SolverRichardson::solve(SparseMatrix<Number> &A,
++ Vector<Number> &x,
++ Vector<Number> &b,
++ PreconditionBase<Number> &preconditioner,
++ bool move_to_accelerator)
+ {
+ std_cxx1x::shared_ptr<paralution::FixedPoint<paralution::LocalMatrix<Number>,
+ paralution::LocalVector<Number>,Number> > solver(new paralution::
+ FixedPoint<paralution::LocalMatrix<Number>, paralution::LocalVector<Number>,Number>);
+
+ // Set the verbosity of the solver.
+ solver->Verbose(additional_data.verbose);
+
+ // Set the relaxation parameter.
+ solver->SetRelaxation(additional_data.relaxation_parameter);
+
+ this->execute_solve<Number>(solver,A,x,b,preconditioner,move_to_accelerator);
+ }
+
+
+
+ /* ---------------------- SolverCG ------------------------- */
+
+ SolverCG::AdditionalData::AdditionalData(const unsigned int verbose)
+ :
+ verbose(verbose)
+ {}
+
+
+
+ SolverCG::SolverCG (SolverControl &cn, const AdditionalData &data)
+ :
+ SolverBase (cn),
+ additional_data (data)
+ {}
+
+
+
+ template <typename Number>
- void SolverCR::solve(const SparseMatrix<Number> &A,
- Vector<Number> &x,
- const Vector<Number> &b,
- const PreconditionBase<Number> &preconditioner,
- bool move_to_accelerator)
++ void SolverCG::solve(SparseMatrix<Number> &A,
++ Vector<Number> &x,
++ Vector<Number> &b,
++ PreconditionBase<Number> &preconditioner,
++ bool move_to_accelerator)
+ {
+ std_cxx1x::shared_ptr<paralution::CG<paralution::LocalMatrix<Number>,
+ paralution::LocalVector<Number>,Number> > solver(new paralution::
+ CG<paralution::LocalMatrix<Number>,paralution::LocalVector<Number>,Number>);
+
+ // Set the verbosity of the solver.
+ solver->Verbose(additional_data.verbose);
+
+ this->execute_solve<Number>(solver,A,x,b,preconditioner,move_to_accelerator);
+ }
+
+
+
+ /* ---------------------- SolverCR ------------------------- */
+
+ SolverCR::AdditionalData::AdditionalData(const unsigned int verbose)
+ :
+ verbose(verbose)
+ {}
+
+
+
+ SolverCR::SolverCR (SolverControl &cn, const AdditionalData &data)
+ :
+ SolverBase (cn),
+ additional_data (data)
+ {}
+
+
+
+ template <typename Number>
- const const unsigned int n_deflated_vectors)
++ void SolverCR::solve(SparseMatrix<Number> &A,
++ Vector<Number> &x,
++ Vector<Number> &b,
++ PreconditionBase<Number> &preconditioner,
++ bool move_to_accelerator)
+ {
+ std_cxx1x::shared_ptr<paralution::CR<paralution::LocalMatrix<Number>,
+ paralution::LocalVector<Number>,Number> > solver(new paralution::
+ CR<paralution::LocalMatrix<Number>,paralution::LocalVector<Number>,Number>);
+
+ // Set the verbosity of the solver.
+ solver->Verbose(additional_data.verbose);
+
+ this->execute_solve<Number>(solver,A,x,b,preconditioner,move_to_accelerator);
+ }
+
+
+
+ /* ---------------------- SolverDPCG ----------------------- */
+
+ SolverDPCG::AdditionalData::AdditionalData(const unsigned int verbose,
- void SolverDPCG::solve(const SparseMatrix<Number> &A,
- Vector<Number> &x,
- const Vector<Number> &b,
- const PreconditionBase<Number> &preconditioner,
- bool move_to_accelerator)
++ const unsigned int n_deflated_vectors)
+ :
+ verbose (verbose),
+ n_deflated_vectors(n_deflated_vectors)
+ {}
+
+
+
+ SolverDPCG::SolverDPCG (SolverControl &cn, const AdditionalData &data)
+ :
+ SolverBase (cn),
+ additional_data (data)
+ {}
+
+
+
+ template <typename Number>
- void SolverBicgstab::solve(const SparseMatrix<Number> &A,
- Vector<Number> &x,
- const Vector<Number> &b,
- const PreconditionBase<Number> &preconditioner,
- bool move_to_accelerator)
++ void SolverDPCG::solve(SparseMatrix<Number> &A,
++ Vector<Number> &x,
++ Vector<Number> &b,
++ PreconditionBase<Number> &preconditioner,
++ bool move_to_accelerator)
+ {
+ std_cxx1x::shared_ptr<paralution::DPCG<paralution::LocalMatrix<Number>,
+ paralution::LocalVector<Number>,Number> > solver(new paralution::
+ DPCG<paralution::LocalMatrix<Number>, paralution::LocalVector<Number>,Number>);
+
+ // Set the verbosity of the solver.
+ solver->Verbose(additional_data.verbose);
+
+ // Set the number of deflated vectors.
+ solver->SetNVectors(additional_data.n_deflated_vectors);
+
+ this->execute_solve<Number>(solver,A,x,b,preconditioner,move_to_accelerator);
+ }
+
+
+
+ /* -------------------- SolverBicgstab --------------------- */
+
+ SolverBicgstab::AdditionalData::AdditionalData(const unsigned int verbose)
+ :
+ verbose(verbose)
+ {}
+
+
+
+ SolverBicgstab::SolverBicgstab (SolverControl &cn, const AdditionalData &data)
+ :
+ SolverBase (cn),
+ additional_data (data)
+ {}
+
+
+
+ template <typename Number>
- void SolverGMRES::solve(const SparseMatrix<Number> &A,
- Vector<Number> &x,
- const Vector<Number> &b,
- const PreconditionBase<Number> &preconditioner,
- bool move_to_accelerator)
++ void SolverBicgstab::solve(SparseMatrix<Number> &A,
++ Vector<Number> &x,
++ Vector<Number> &b,
++ PreconditionBase<Number> &preconditioner,
++ bool move_to_accelerator)
+ {
+ std_cxx1x::shared_ptr<paralution::BiCGStab<paralution::LocalMatrix<Number>,
+ paralution::LocalVector<Number>,Number> > solver(new paralution::
+ BiCGStab<paralution::LocalMatrix<Number>,paralution::LocalVector<Number>,Number>);
+
+ // Set the verbosity of the solver.
+ solver->Verbose(additional_data.verbose);
+
+ this->execute_solve<Number>(solver,A,x,b,preconditioner,move_to_accelerator);
+ }
+
+
+
+ /* --------------------- SolverGMRES ----------------------- */
+
+ SolverGMRES::AdditionalData::AdditionalData(const unsigned int verbose,
+ const unsigned int restart_parameter)
+ :
+ verbose(verbose),
+ restart_parameter(restart_parameter)
+ {}
+
+
+
+ SolverGMRES::SolverGMRES (SolverControl &cn,
+ const AdditionalData &data)
+ :
+ SolverBase (cn),
+ additional_data (data)
+ {}
+
+
+
+ template <typename Number>
- void SolverFGMRES::solve(const SparseMatrix<Number> &A,
- Vector<Number> &x,
- const Vector<Number> &b,
- const PreconditionBase<Number> &preconditioner,
- bool move_to_accelerator)
++ void SolverGMRES::solve(SparseMatrix<Number> &A,
++ Vector<Number> &x,
++ Vector<Number> &b,
++ PreconditionBase<Number> &preconditioner,
++ bool move_to_accelerator)
+ {
+ std_cxx1x::shared_ptr<paralution::GMRES<paralution::LocalMatrix<Number>,
+ paralution::LocalVector<Number>,Number> > solver(new paralution::
+ GMRES<paralution::LocalMatrix<Number>, paralution::LocalVector<Number>,Number>);
+
+ // Set the verbosity of the solver.
+ solver->Verbose(additional_data.verbose);
+
+ // Set the restart parameter.
+ solver->SetBasisSize(additional_data.restart_parameter);
+
+ this->execute_solve<Number>(solver,A,x,b,preconditioner,move_to_accelerator);
+ }
+
+
+
+ /* --------------------- SolverFGMRES ---------------------- */
+
+ SolverFGMRES::AdditionalData::AdditionalData(const unsigned int verbose,
+ const unsigned int restart_parameter)
+ :
+ verbose(verbose),
+ restart_parameter(restart_parameter)
+ {}
+
+
+
+ SolverFGMRES::SolverFGMRES (SolverControl &cn,
+ const AdditionalData &data)
+ :
+ SolverBase (cn),
+ additional_data (data)
+ {}
+
+
+
+ template <typename Number>
- template void SolverCG::solve<float>(const SparseMatrix<float> &A,
- Vector<float> &x,
- const Vector<float> &b,
- const PreconditionBase<float> &preconditioner,
- bool move_to_accelerator);
-
- template void SolverCG::solve<double>(const SparseMatrix<double> &A,
- Vector<double> &x,
- const Vector<double> &b,
- const PreconditionBase<double> &preconditioner,
- bool move_to_accelerator);
-
- template void SolverCR::solve<float>(const SparseMatrix<float> &A,
- Vector<float> &x,
- const Vector<float> &b,
- const PreconditionBase<float> &preconditioner,
- bool move_to_accelerator);
-
- template void SolverCR::solve<double>(const SparseMatrix<double> &A,
- Vector<double> &x,
- const Vector<double> &b,
- const PreconditionBase<double> &preconditioner,
- bool move_to_accelerator);
-
- template void SolverDPCG::solve<float>(const SparseMatrix<float> &A,
- Vector<float> &x,
- const Vector<float> &b,
- const PreconditionBase<float> &preconditioner,
- bool move_to_accelerator);
-
- template void SolverDPCG::solve<double>(const SparseMatrix<double> &A,
- Vector<double> &x,
- const Vector<double> &b,
- const PreconditionBase<double> &preconditioner,
- bool move_to_accelerator);
-
- template void SolverBicgstab::solve<float>(const SparseMatrix<float> &A,
- Vector<float> &x,
- const Vector<float> &b,
- const PreconditionBase<float> &preconditioner,
- bool move_to_accelerator);
-
- template void SolverBicgstab::solve<double>(const SparseMatrix<double> &A,
- Vector<double> &x,
- const Vector<double> &b,
- const PreconditionBase<double> &preconditioner,
- bool move_to_accelerator);
-
- template void SolverGMRES::solve<float>(const SparseMatrix<float> &A,
- Vector<float> &x,
- const Vector<float> &b,
- const PreconditionBase<float> &preconditioner,
- bool move_to_accelerator);
-
- template void SolverGMRES::solve<double>(const SparseMatrix<double> &A,
- Vector<double> &x,
- const Vector<double> &b,
- const PreconditionBase<double> &preconditioner,
- bool move_to_accelerator);
-
- template void SolverFGMRES::solve<float>(const SparseMatrix<float> &A,
- Vector<float> &x,
- const Vector<float> &b,
- const PreconditionBase<float> &preconditioner,
- bool move_to_accelerator);
-
- template void SolverFGMRES::solve<double>(const SparseMatrix<double> &A,
- Vector<double> &x,
- const Vector<double> &b,
- const PreconditionBase<double> &preconditioner,
- bool move_to_accelerator);
++ void SolverFGMRES::solve(SparseMatrix<Number> &A,
++ Vector<Number> &x,
++ Vector<Number> &b,
++ PreconditionBase<Number> &preconditioner,
++ bool move_to_accelerator)
+ {
+ std_cxx1x::shared_ptr<paralution::FGMRES<paralution::LocalMatrix<Number>,
+ paralution::LocalVector<Number>,Number> > solver(new paralution::
+ FGMRES<paralution::LocalMatrix<Number>, paralution::LocalVector<Number>,Number>);
+
+ // Set the verbosity of the solver.
+ solver->Verbose(additional_data.verbose);
+
+ // Set the restart parameter.
+ solver->SetBasisSize(additional_data.restart_parameter);
+
+ this->execute_solve<Number>(solver,A,x,b,preconditioner,move_to_accelerator);
+ }
+}
+
+
+
+// Explicit instantiations
+namespace ParalutionWrappers
+{
++ template void SolverCG::solve<float>(SparseMatrix<float> &A,
++ Vector<float> &x,
++ Vector<float> &b,
++ PreconditionBase<float> &preconditioner,
++ bool move_to_accelerator);
++
++ template void SolverCG::solve<double>(SparseMatrix<double> &A,
++ Vector<double> &x,
++ Vector<double> &b,
++ PreconditionBase<double> &preconditioner,
++ bool move_to_accelerator);
++
++ template void SolverCR::solve<float>(SparseMatrix<float> &A,
++ Vector<float> &x,
++ Vector<float> &b,
++ PreconditionBase<float> &preconditioner,
++ bool move_to_accelerator);
++
++ template void SolverCR::solve<double>(SparseMatrix<double> &A,
++ Vector<double> &x,
++ Vector<double> &b,
++ PreconditionBase<double> &preconditioner,
++ bool move_to_accelerator);
++
++ template void SolverDPCG::solve<float>(SparseMatrix<float> &A,
++ Vector<float> &x,
++ Vector<float> &b,
++ PreconditionBase<float> &preconditioner,
++ bool move_to_accelerator);
++
++ template void SolverDPCG::solve<double>(SparseMatrix<double> &A,
++ Vector<double> &x,
++ Vector<double> &b,
++ PreconditionBase<double> &preconditioner,
++ bool move_to_accelerator);
++
++ template void SolverBicgstab::solve<float>(SparseMatrix<float> &A,
++ Vector<float> &x,
++ Vector<float> &b,
++ PreconditionBase<float> &preconditioner,
++ bool move_to_accelerator);
++
++ template void SolverBicgstab::solve<double>(SparseMatrix<double> &A,
++ Vector<double> &x,
++ Vector<double> &b,
++ PreconditionBase<double> &preconditioner,
++ bool move_to_accelerator);
++
++ template void SolverGMRES::solve<float>(SparseMatrix<float> &A,
++ Vector<float> &x,
++ Vector<float> &b,
++ PreconditionBase<float> &preconditioner,
++ bool move_to_accelerator);
++
++ template void SolverGMRES::solve<double>(SparseMatrix<double> &A,
++ Vector<double> &x,
++ Vector<double> &b,
++ PreconditionBase<double> &preconditioner,
++ bool move_to_accelerator);
++
++ template void SolverFGMRES::solve<float>(SparseMatrix<float> &A,
++ Vector<float> &x,
++ Vector<float> &b,
++ PreconditionBase<float> &preconditioner,
++ bool move_to_accelerator);
++
++ template void SolverFGMRES::solve<double>(SparseMatrix<double> &A,
++ Vector<double> &x,
++ Vector<double> &b,
++ PreconditionBase<double> &preconditioner,
++ bool move_to_accelerator);
+}
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif // DEAL_II_WITH_PARALUTION
projects / dealii-svn.git / commitdiff