/**
* Copy-Constructor. Dimension set to
* that of V, all components are copied
- * from V
+ * from V.
*/
BlockVector (const BlockVector<Number> &V);
const InputIterator end);
/**
- * Destructor. Clears memory
+ * Destructor. Clears memory.
*/
~BlockVector ();
namespace ParalutionWrappers
{
/**
- * Cycle for multigrid preconditioner.
+ * Cycle for multigrid.
*/
enum mg_cycle {V_cycle, W_cycle, K_cycle, F_cycle};
enum mg_interpolation {aggregation, smoothed_aggregation};
/**
- * Precondition for the smoother.
+ * Preconditioner for the smoother of the multigrid.
*/
enum mg_preconditioner {jacobi, sgs, multicolored_sgs, multicolored_sor, ilu,
- ilut, multicolored_ilu
- };
+ ilut, multicolored_ilu};
/**
- * Solver used as smoother and to solve the coarse system.
+ * Solver used as smoother and to solve the coarse system of the multigrid.
*/
enum mg_solver {richardson, cg, cr, bicgstab, gmres, fgmres};
*/
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.
+ * Constructor. @p relaxation_parameter is used only with smoothed
+ * aggregation and is disregarded otherwise. @p over_interpolation is used
+ * only with aggregation and is disregared otherwise.
*/
- bool default_parameters;
+ AdditionalData (const unsigned int verbose = 0,
+ const mg_solver coarse_solver = gmres,
+ const unsigned int n_unknowns_coarse_level = 300,
+ const unsigned int max_iter = 2,
+ const mg_cycle cycle = V_cycle,
+ const mg_interpolation interpolation = smoothed_aggregation,
+ const double coupling_strength = 0.01,
+ const double relaxation_parameter = 2./3.,
+ const double over_interpolation = 1.5);
/**
- * Cycle used by the multigrid preconditioner: V cycle, W cycle, K
- * cycle, or F cycle.
+ * Verbosity: 0 = no output, 1 = print information at the
+ * start and at the end, 3 = print residual at each iteration.
*/
- mg_cycle cycle;
+ unsigned int verbose;
/**
- * Interpolation used by the multigrid preconditioner: smoothed
- * aggregation or aggregation.
+ * Solver of the coarse system.
*/
- mg_interpolation interpolation;
+ mg_solver coarse_solver;
/**
- * Preconditioner for the smoother.
+ * Number of unknowns on the coarsest level.
*/
- mg_preconditioner preconditioner;
+ unsigned int n_unknowns_coarse_level;
/**
- * Smoother.
+ * Maximum number of iterations.
*/
- mg_solver smoother;
+ unsigned int max_iter;
/**
- * Solver of the coarse system.
+ * Cycle used by the multigrid preconditioner: V cycle, W cycle, K
+ * cycle, or F cycle.
*/
- mg_solver coarse_solver;
+ mg_cycle cycle;
/**
- * Verbosity: 0 = no output, 1 = print information at the
- * start and at the end, 3 = print residual at each iteration.
+ * Interpolation used by the multigrid preconditioner: smoothed
+ * aggregation or aggregation.
*/
- unsigned int verbose;
+ mg_interpolation interpolation;
/**
- * Number of unknowns on the coarsest level.
+ * Coupling strength.
*/
- unsigned int n_unknowns_coarse_level;
+ double coupling_strength;
/**
- * Maximum number of iterations.
+ * Relaxation parameter for smoothed interpolation aggregation.
*/
- unsigned int max_iter;
+ double relaxation_parameter;
/**
- * Relaxation parameter for smoothed interpolation aggregation.
+ * Over-interpolation parameter for aggregation.
*/
- double relaxation_parameter;
+ double over_interpolation;
};
/**
private :
/**
- * Free the smoothers, the preconditioners, and the coarse level solver.
+ * Free 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.
*/
namespace ParalutionWrappers
{
/**
- * Base class for solver classes using the Paralution solvers.
+ * Base class for solver classes using the Paralution solvers except AMG which
+ * is not derived from this class.
*
* @ingroup ParalutionWrappers
* @author Bruno Turcksin 2013
bool move_to_accelerator);
/**
- * Reference to the object that controls convergence of the iteratove
+ * Reference to the object that controls convergence of the iterative
* 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.
+ * itself, but we copy the data from this object before starting solution
+ * process, and copy the data back into it afterwards.
*/
SolverControl &solver_control;
};
private:
/**
- * store a copy of the flags for this solver.
+ * Store a copy of the flags for this solver.
*/
const AdditionalData additional_data;
};
private:
/**
- * store a copy of the flags for this solver.
+ * Store a copy of the flags for this solver.
*/
const AdditionalData additional_data;
};
private:
/**
- * store a copy of the flags for this solver.
+ * Store a copy of the flags for this solver.
*/
const AdditionalData additional_data;
};
private:
/**
- * store a copy of the flags for this solver.
+ * Store a copy of the flags for this solver.
*/
const AdditionalData additional_data;
};
private:
/**
- * store a copy of the flags for this solver.
+ * Store a copy of the flags for this solver.
*/
const AdditionalData additional_data;
};
private:
/**
- * store a copy of the flags for this solver.
+ * Store a copy of the flags for this solver.
*/
const AdditionalData additional_data;
};
private:
/**
- * store a copy of the flags for this solver.
+ * 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 SolverAMG
+ {
+ public:
+ /**
+ * Standardized data struct to pipe additional data to the solver.
+ */
+ struct AdditionalData
+ {
+ /**
+ * Constructor. @p relaxation_parameter is used only with smoothed
+ * aggregation and is disregarded otherwise. @p over_interpolation is used
+ * only with aggregation and is disregared otherwise.
+ */
+ AdditionalData (const unsigned int verbose = 0,
+ const mg_solver coarse_solver = gmres,
+ const unsigned int n_unknowns_coarse_level = 300,
+ const mg_solver smoother = richardson,
+ const mg_preconditioner = multicolored_sor,
+ const mg_cycle cycle = V_cycle,
+ const mg_interpolation interpolation = smoothed_aggregation,
+ const double coupling_strength = 0.01,
+ const double relaxation_parameter = 2./3.,
+ const double over_interpolation = 1.5);
+
+ /**
+ * Verbosity: 0 = no output, 1 = print information at the
+ * start and at the end, 3 = print residual at each iteration.
+ */
+ unsigned int verbose;
+
+ /**
+ * Solver of the coarse system.
+ */
+ mg_solver coarse_solver;
+
+ /**
+ * Number of unknowns on the coarsest level.
+ */
+ unsigned int n_unknowns_coarse_level;
+
+ /**
+ * Smoother (pre- and post-smoothers).
+ */
+ mg_solver smoother;
+
+ /**
+ * Preconditioner of the smoother.
+ */
+ mg_preconditioner preconditioner;
+
+ /**
+ * 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;
+
+ /**
+ * Coupling strength.
+ */
+ double coupling_strength;
+
+ /**
+ * Relaxation parameter for smoothed interpolation aggregation.
+ */
+ double relaxation_parameter;
+
+ /**
+ * Over-interpolation parameter for aggregation.
+ */
+ double over_interpolation;
+ };
+
+ /**
+ * Constructor. AdditionalData is a structure that contains additional
+ * flags for tuning this particular solver.
+ */
+ SolverAMG (SolverControl &cn,
+ const AdditionalData &data = AdditionalData());
+
+ /**
+ * Solve the linear system <tt>Ax=b</tt> using the AMG solver of
+ * Paralution. If the flag @p move_to_accelerator is set to true, the
+ * solver, the right-hand side, the matrix, and the solution vector are
+ * moved on the accelerator once the solver is built. The multigrid must be
+ * build on the CPU.
+ */
+ template <typename Number>
+ void solve (SparseMatrix<Number> &A,
+ Vector<Number> &x,
+ Vector<Number> &b,
+ bool move_to_accelerator=false);
+
+ /**
+ * Access to object that controls convergence.
+ */
+ SolverControl &control() const;
+
+ private:
+ /**
+ * Reference to the object that controls convergence of the iterative
+ * solver. In fact, for these Paralution wrappers, Paralution does so
+ * itself, but we copy the data from this object before starting solution
+ * process, and copy the data back into it afterwards.
+ */
+ SolverControl &solver_control;
+
+ /**
+ * Store a copy of the flags for this solver.
*/
const AdditionalData additional_data;
};
* Change the dimension of the vector to @p N. The vector is filled with
* zeros.
*/
- //TODO look to add fast
void reinit(const size_type N);
/**
*/
void move_to_host_async();
+ /**
+ * Print to a stream. @p precision denotes the desired precision with
+ * which values shall be printed, @p scientific whether scientific
+ * notation shall be used. If @p across is @p true then the vector is
+ * printed in a line, while if @p false then the elements are printed on a
+ * separate line each.
+ */
+ void print (std::ostream &out,
+ const unsigned int precision,
+ const bool scientific,
+ const bool across) const;
+
/**
* Synchronize the code when move_to_host_async or move_to_accelerator_async
* is used.
/* -------------------------- 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)
+ const unsigned int max_iter,
+ const mg_cycle cycle,
+ const mg_interpolation interpolation,
+ const double coupling_strength,
+ const double relaxation_parameter,
+ const double over_interpolation)
:
- default_parameters(false),
- cycle(cycle),
- interpolation(interpolation),
- preconditioner(preconditioner),
- smoother(smoother),
- coarse_solver(coarse_solver),
verbose(verbose),
+ coarse_solver(coarse_solver),
n_unknowns_coarse_level(n_unknowns_coarse_level),
max_iter(max_iter),
- relaxation_parameter(relaxation_parameter)
+ cycle(cycle),
+ interpolation(interpolation),
+ coupling_strength(coupling_strength),
+ relaxation_parameter(relaxation_parameter),
+ over_interpolation(over_interpolation)
{}
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>,
{
// 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());
+ 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
+
+ // Set the verbosity of the multigrid.
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);
+ // Free 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);
+ downcasted_ptr->SetOverInterp(additional_data.over_interpolation);
+ 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 coupling strength
+ downcasted_ptr->SetCouplingStrength(additional_data.coupling_strength);
+
+ // 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 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)
- {
+ // Coarse solver
+ switch (additional_data.coarse_solver)
+ {
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::FixedPoint<paralution::LocalMatrix<Number>, paralution::LocalVector<Number>,
+ Number> *cs = new paralution::FixedPoint<paralution::LocalMatrix<Number>,
paralution::LocalVector<Number>, Number>;
- smoothers[i] = sm;
- }
- break;
- }
+ coarse_solver = cs;
+ 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::CG<paralution::LocalMatrix<Number>, paralution::LocalVector<Number>,
+ Number> *cs = new paralution::CG<paralution::LocalMatrix<Number>,
paralution::LocalVector<Number>, Number>;
- smoothers[i] = sm;
- }
- break;
- }
+ coarse_solver = cs;
+ 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::CR<paralution::LocalMatrix<Number>, paralution::LocalVector<Number>,
+ Number> *cs = new paralution::CR<paralution::LocalMatrix<Number>,
paralution::LocalVector<Number>, Number>;
- smoothers[i] = sm;
- }
- break;
- }
+ coarse_solver = cs;
+ 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::BiCGStab<paralution::LocalMatrix<Number>, paralution::LocalVector<Number>,
+ Number> *cs = new paralution::BiCGStab<paralution::LocalMatrix<Number>,
paralution::LocalVector<Number>, Number>;
- smoothers[i] = sm;
- }
- break;
- }
+ coarse_solver = cs;
+ 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::GMRES<paralution::LocalMatrix<Number>, paralution::LocalVector<Number>,
+ Number> *cs = new paralution::GMRES<paralution::LocalMatrix<Number>,
paralution::LocalVector<Number>, Number>;
- smoothers[i] = sm;
- }
- break;
- }
+ coarse_solver = cs;
+ 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::FGMRES<paralution::LocalMatrix<Number>, paralution::LocalVector<Number>,
+ Number> *cs = new paralution::FGMRES<paralution::LocalMatrix<Number>,
paralution::LocalVector<Number>, Number>;
- preconditioners[i] = precond;
- }
- break;
- }
+ coarse_solver = cs;
+ 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."));
- }
- }
+ {
+ AssertThrow(false,ExcMessage("Unknown coarse solver for PreconditionAMG."));
+ }
+ }
+ // Set the coarse solver and the verbosity of the coarse solver.
+ coarse_solver->Verbose(additional_data.verbose);
downcasted_ptr->SetSolver(*coarse_solver);
- }
}
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;
- }
}
}
- SolverControl &SolverBase::control() const
+ SolverControl& SolverBase::control() const
{
return solver_control;
}
this->execute_solve<Number>(solver,A,x,b,preconditioner,move_to_accelerator);
}
+
+
+
+ /* ---------------------- SolverAMG------------------------- */
+
+ SolverAMG::AdditionalData::AdditionalData(const unsigned int verbose,
+ const mg_solver coarse_solver,
+ const unsigned int n_unknowns_coarse_level,
+ const mg_solver smoother,
+ const mg_preconditioner preconditioner,
+ const mg_cycle cycle,
+ const mg_interpolation interpolation,
+ const double coupling_strength,
+ const double relaxation_parameter,
+ const double over_interpolation)
+ :
+ verbose(verbose),
+ coarse_solver(coarse_solver),
+ n_unknowns_coarse_level(n_unknowns_coarse_level),
+ smoother(smoother),
+ preconditioner(preconditioner),
+ cycle(cycle),
+ interpolation(interpolation),
+ coupling_strength(coupling_strength),
+ relaxation_parameter(relaxation_parameter),
+ over_interpolation(over_interpolation)
+ {}
+
+
+
+ SolverAMG::SolverAMG (SolverControl &cn,
+ const AdditionalData &data)
+ :
+ solver_control (cn),
+ additional_data (data)
+ {}
+
+
+
+ template <typename Number>
+ void SolverAMG::solve(SparseMatrix<Number> &A,
+ Vector<Number> &x,
+ Vector<Number> &b,
+ bool move_to_accelerator)
+ {
+ paralution::AMG<paralution::LocalMatrix<Number>,paralution::LocalVector<Number>,Number> solver;
+
+ // Set the system to solve.
+ solver.SetOperator(A.paralution_matrix());
+ solver.SetOperatorFormat(CSR);
+
+ // Set the verbosity of the multigrid.
+ solver.Verbose(additional_data.verbose);
+
+ // Set number of unknowns on coarsest level.
+ solver.SetCoarsestLevel(additional_data.n_unknowns_coarse_level);
+
+ // Set the interpolation type.
+ switch (additional_data.interpolation)
+ {
+ case aggregation :
+ {
+ solver.SetInterpolation(paralution::Aggregation);
+ solver.SetOverInterp(additional_data.over_interpolation);
+ break;
+ }
+ case smoothed_aggregation :
+ {
+ solver.SetInterpolation(paralution::SmoothedAggregation);
+ solver.SetInterpRelax(additional_data.relaxation_parameter);
+ break;
+ }
+ default :
+ {
+ AssertThrow(false,ExcMessage("Unknown interpolation type for PreconditionAMG."));
+ }
+ }
+
+ // Set the coupling strength.
+ solver.SetCouplingStrength(additional_data.coupling_strength);
+
+ // Set the type of cycle
+ switch (additional_data.cycle)
+ {
+ case V_cycle :
+ {
+ solver.SetCycle(paralution::Vcycle);
+ break;
+ }
+ case W_cycle :
+ {
+ solver.SetCycle(paralution::Wcycle);
+ break;
+ }
+ case K_cycle :
+ {
+ solver.SetCycle(paralution::Kcycle);
+ break;
+ }
+ case F_cycle :
+ {
+ solver.SetCycle(paralution::Fcycle);
+ }
+ default :
+ {
+ AssertThrow(false,ExcMessage("Unknown cycle type for PreconditionAMG."));
+ }
+ }
+
+ // Manual smoothers
+ solver.SetManualSmoothers(true);
+
+ // Manual course grid solver
+ solver.SetManualSolver(true);
+
+ // Set grid transfer scaling
+ solver.SetScaling(true);
+
+ // Build the grid hierarchy.
+ solver.BuildHierarchy();
+
+ unsigned int n_levels = solver.GetNumLevels();
+
+ // Smoother for each level
+ paralution::IterativeLinearSolver<paralution::LocalMatrix<Number>,
+ paralution::LocalVector<Number>,Number> **smoothers = NULL;
+ smoothers = new paralution::IterativeLinearSolver<paralution::LocalMatrix<Number>,
+ paralution::LocalVector<Number>,Number>* [n_levels-1];
+
+ // Preconditioner for each smoother
+ paralution::Preconditioner<paralution::LocalMatrix<Number>,
+ paralution::LocalVector<Number>,Number> **preconds = NULL;
+ preconds = new paralution::Preconditioner<paralution::LocalMatrix<Number>,
+ paralution::LocalVector<Number>,Number>* [n_levels-1];
+
+ // Coarse Grid Solver
+ paralution::IterativeLinearSolver<paralution::LocalMatrix<Number>,
+ paralution::LocalVector<Number>,Number> *coarse_solver;
+
+ // Smoother
+ 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> *richardson = NULL;
+ richardson = new paralution::FixedPoint<paralution::LocalMatrix<Number>,
+ paralution::LocalVector<Number>,Number>;
+ smoothers[i] = richardson;
+ smoothers[i]->Verbose(additional_data.verbose);
+ }
+
+ break;
+ }
+ case cg :
+ {
+ for (unsigned int i=0; i<n_levels-1; ++i)
+ {
+ paralution::CG<paralution::LocalMatrix<Number>,paralution::LocalVector<Number>,
+ Number> *cg = NULL;
+ cg = new paralution::CG<paralution::LocalMatrix<Number>,paralution::LocalVector<Number>,
+ Number>;
+ smoothers[i] = cg;
+ smoothers[i]->Verbose(additional_data.verbose);
+ }
+
+ break;
+ }
+ case cr :
+ {
+ for (unsigned int i=0; i<n_levels-1; ++i)
+ {
+ paralution::CR<paralution::LocalMatrix<Number>,paralution::LocalVector<Number>,
+ Number> *cr = NULL;
+ cr = new paralution::CR<paralution::LocalMatrix<Number>,paralution::LocalVector<Number>,
+ Number>;
+ smoothers[i] = cr;
+ smoothers[i]->Verbose(additional_data.verbose);
+ }
+
+ break;
+ }
+ case bicgstab :
+ {
+ for (unsigned int i=0; i<n_levels-1; ++i)
+ {
+ paralution::BiCGStab<paralution::LocalMatrix<Number>,paralution::LocalVector<Number>,
+ Number> *bicgstab = NULL;
+ bicgstab = new paralution::BiCGStab<paralution::LocalMatrix<Number>,
+ paralution::LocalVector<Number>,Number>;
+ smoothers[i] = bicgstab;
+ smoothers[i]->Verbose(additional_data.verbose);
+ }
+
+ break;
+ }
+ case gmres :
+ {
+ for (unsigned int i=0; i<n_levels-1; ++i)
+ {
+ paralution::GMRES<paralution::LocalMatrix<Number>,paralution::LocalVector<Number>,
+ Number> *gmres = NULL;
+ gmres = new paralution::GMRES<paralution::LocalMatrix<Number>,
+ paralution::LocalVector<Number>,Number>;
+ smoothers[i] = gmres;
+ smoothers[i]->Verbose(additional_data.verbose);
+ }
+
+ break;
+ }
+ case fgmres :
+ {
+ for (unsigned int i=0; i<n_levels-1; ++i)
+ {
+ paralution::FGMRES<paralution::LocalMatrix<Number>,paralution::LocalVector<Number>,
+ Number> *fgmres = NULL;
+ fgmres = new paralution::FGMRES<paralution::LocalMatrix<Number>,
+ paralution::LocalVector<Number>,Number>;
+ smoothers[i] = fgmres;
+ smoothers[i]->Verbose(additional_data.verbose);
+ }
+
+ break;
+ }
+ default :
+ {
+ AssertThrow(false,ExcMessage("Unknown smoother for AMG."));
+ }
+ }
+
+ // Preconditioner
+ 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> *jacobi = NULL;
+ jacobi = new paralution::Jacobi<paralution::LocalMatrix<Number>,
+ paralution::LocalVector<Number>,Number>;
+ preconds[i] = jacobi;
+ smoothers[i]->SetPreconditioner(*preconds[i]);
+ }
+
+ break;
+ }
+ case sgs :
+ {
+ for (unsigned int i=0; i<n_levels-1; ++i)
+ {
+ paralution::SGS<paralution::LocalMatrix<Number>,paralution::LocalVector<Number>,
+ Number> *sgs = NULL;
+ sgs = new paralution::SGS<paralution::LocalMatrix<Number>,paralution::LocalVector<Number>,
+ Number>;
+ preconds[i] = sgs;
+ smoothers[i]->SetPreconditioner(*preconds[i]);
+ }
+
+ break;
+ }
+ case multicolored_sgs :
+ {
+ for (unsigned int i=0; i<n_levels-1; ++i)
+ {
+ paralution::MultiColoredSGS<paralution::LocalMatrix<Number>,paralution::LocalVector<Number>,
+ Number> *multicolored_sgs = NULL;
+ multicolored_sgs = new paralution::MultiColoredSGS<paralution::LocalMatrix<Number>,
+ paralution::LocalVector<Number>,Number>;
+ multicolored_sgs->SetPrecondMatrixFormat(ELL);
+ preconds[i] = multicolored_sgs;
+ smoothers[i]->SetPreconditioner(*preconds[i]);
+ }
+
+ break;
+ }
+ case multicolored_sor :
+ {
+ for (unsigned int i=0; i<n_levels-1; ++i)
+ {
+ paralution::MultiColoredGS<paralution::LocalMatrix<Number>,paralution::LocalVector<Number>,
+ Number> *multicolored_sor = NULL;
+ multicolored_sor = new paralution::MultiColoredGS<paralution::LocalMatrix<Number>,
+ paralution::LocalVector<Number>,Number>;
+ multicolored_sor->SetPrecondMatrixFormat(ELL);
+ preconds[i] = multicolored_sor;
+ smoothers[i]->SetPreconditioner(*preconds[i]);
+ }
+
+ break;
+ }
+ case ilu :
+ {
+ for (unsigned int i=0; i<n_levels-1; ++i)
+ {
+ paralution::ILU<paralution::LocalMatrix<Number>,paralution::LocalVector<Number>,
+ Number> *ilu = NULL;
+ ilu = new paralution::ILU<paralution::LocalMatrix<Number>,paralution::LocalVector<Number>,
+ Number>;
+ preconds[i] = ilu;
+ smoothers[i]->SetPreconditioner(*preconds[i]);
+ }
+
+ break;
+ }
+ case ilut :
+ {
+ for (unsigned int i=0; i<n_levels-1; ++i)
+ {
+ paralution::ILUT<paralution::LocalMatrix<Number>,paralution::LocalVector<Number>,
+ Number> *ilut = NULL;
+ ilut = new paralution::ILUT<paralution::LocalMatrix<Number>,
+ paralution::LocalVector<Number>,Number>;
+ preconds[i] = ilut;
+ smoothers[i]->SetPreconditioner(*preconds[i]);
+ }
+
+ break;
+ }
+ case multicolored_ilu :
+ {
+ for (unsigned int i=0; i<n_levels-1; ++i)
+ {
+ paralution::MultiColoredILU<paralution::LocalMatrix<Number>,paralution::LocalVector<Number>,
+ Number> *multicolored_ilu = NULL;
+ multicolored_ilu = new paralution::MultiColoredILU<paralution::LocalMatrix<Number>,
+ paralution::LocalVector<Number>,Number>;
+ multicolored_ilu->SetPrecondMatrixFormat(ELL);
+ preconds[i] = multicolored_ilu;
+ smoothers[i]->SetPreconditioner(*preconds[i]);
+ }
+
+ break;
+ }
+ default :
+ {
+ AssertThrow(false,ExcMessage("Unknown preconditioner for the smoother for AMG."));
+ }
+ }
+
+ // Set the smoothers.
+ solver.SetSmoother(smoothers);
+ solver.SetSmootherPreIter(1);
+ solver.SetSmootherPostIter(2);
+
+ // 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 coarse solver for PreconditionAMG."));
+ }
+ }
+ solver.SetSolver(*coarse_solver);
+
+ // 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)
+ {
+ A.move_to_accelerator();
+ x.move_to_accelerator();
+ b.move_to_accelerator();
+ solver.MoveToAccelerator();
+ }
+ 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()));
+
+ // Free the coarse solver, the preconditioners and the smoothers.
+ if (coarse_solver!=NULL)
+ {
+ delete coarse_solver;
+ coarse_solver = NULL;
+ }
+ if (preconds!=NULL)
+ {
+ for (unsigned int i=0; i<n_levels-1; ++i)
+ delete [] preconds[i];
+ delete [] preconds;
+ preconds = NULL;
+ }
+ if (smoothers!=NULL)
+ {
+ for (unsigned int i=0; i<n_levels-1; ++i)
+ delete [] smoothers[i];
+ delete [] smoothers;
+ smoothers = NULL;
+ }
+ }
+
+
+
+ SolverControl& SolverAMG::control() const
+ {
+ return solver_control;
+ }
}
Vector<double> &b,
PreconditionBase<double> &preconditioner,
bool move_to_accelerator);
+
+ template void SolverAMG::solve<float>(SparseMatrix<float> &A,
+ Vector<float> &x,
+ Vector<float> &b,
+ bool move_to_accelerator);
+
+ template void SolverAMG::solve<double>(SparseMatrix<double> &A,
+ Vector<double> &x,
+ Vector<double> &b,
+ bool move_to_accelerator);
}
DEAL_II_NAMESPACE_CLOSE
local_vector.ScaleAdd2(0.,u.paralution_vector(),a,v.paralution_vector(),b);
local_vector.AddScale(w.paralution_vector(),c);
}
+
+
+
+ template <typename Number>
+ void Vector<Number>::print (std::ostream &out,
+ const unsigned int precision,
+ const bool scientific,
+ const bool across) const
+ {
+ AssertThrow (out, ExcIO());
+
+ // get a representation of the vector and loop over all the elements
+ out.precision (precision);
+ if (scientific)
+ out.setf (std::ios::scientific, std::ios::floatfield);
+ else
+ out.setf (std::ios::fixed, std::ios::floatfield);
+
+ if (across)
+ for (size_type i=0; i<size(); ++i)
+ out << static_cast<double>(local_vector[i]) << ' ';
+ else
+ for (size_type i=0; i<size(); ++i)
+ out << static_cast<double>(local_vector[i]) << std::endl;
+ out << std::endl;
+
+ // restore the representation
+ // of the vector
+ AssertThrow (out, ExcIO());
+ }
}
// Explicit instantiations
projects / dealii-svn.git / commitdiff