From: turcksin Date: Thu, 26 Jun 2014 20:49:24 +0000 (+0000) Subject: Refactor AMG preconditioner and add AMG solver. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=6a92637d86ec3a52c5e41a3546b72296d30c33f6;p=dealii-svn.git Refactor AMG preconditioner and add AMG solver. git-svn-id: https://svn.dealii.org/branches/branch_paralution@33091 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/include/deal.II/lac/block_vector.h b/deal.II/include/deal.II/lac/block_vector.h index 77108b5525..64a19b44cc 100644 --- a/deal.II/include/deal.II/lac/block_vector.h +++ b/deal.II/include/deal.II/lac/block_vector.h @@ -112,7 +112,7 @@ public: /** * Copy-Constructor. Dimension set to * that of V, all components are copied - * from V + * from V. */ BlockVector (const BlockVector &V); @@ -195,7 +195,7 @@ public: const InputIterator end); /** - * Destructor. Clears memory + * Destructor. Clears memory. */ ~BlockVector (); diff --git a/deal.II/include/deal.II/lac/paralution_precondition.h b/deal.II/include/deal.II/lac/paralution_precondition.h index f6ffe9e21c..4e94a075cc 100644 --- a/deal.II/include/deal.II/lac/paralution_precondition.h +++ b/deal.II/include/deal.II/lac/paralution_precondition.h @@ -35,7 +35,7 @@ DEAL_II_NAMESPACE_OPEN namespace ParalutionWrappers { /** - * Cycle for multigrid preconditioner. + * Cycle for multigrid. */ enum mg_cycle {V_cycle, W_cycle, K_cycle, F_cycle}; @@ -45,14 +45,13 @@ namespace ParalutionWrappers 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}; @@ -375,79 +374,69 @@ namespace ParalutionWrappers */ 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; }; /** @@ -467,25 +456,10 @@ namespace ParalutionWrappers 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:: - LocalVector,Number> * *smoothers; - /** - * Preconditioners of the smoothers for each level. - */ - paralution::Preconditioner, paralution:: - LocalVector,Number> * *preconditioners; /** * Coarse grid solver. */ diff --git a/deal.II/include/deal.II/lac/paralution_solver.h b/deal.II/include/deal.II/lac/paralution_solver.h index 555d0481c0..9b3a738cdb 100644 --- a/deal.II/include/deal.II/lac/paralution_solver.h +++ b/deal.II/include/deal.II/lac/paralution_solver.h @@ -36,7 +36,8 @@ DEAL_II_NAMESPACE_OPEN 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 @@ -68,10 +69,10 @@ namespace ParalutionWrappers 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; }; @@ -134,7 +135,7 @@ namespace ParalutionWrappers private: /** - * store a copy of the flags for this solver. + * Store a copy of the flags for this solver. */ const AdditionalData additional_data; }; @@ -194,7 +195,7 @@ namespace ParalutionWrappers private: /** - * store a copy of the flags for this solver. + * Store a copy of the flags for this solver. */ const AdditionalData additional_data; }; @@ -249,7 +250,7 @@ namespace ParalutionWrappers private: /** - * store a copy of the flags for this solver. + * Store a copy of the flags for this solver. */ const AdditionalData additional_data; }; @@ -315,7 +316,7 @@ namespace ParalutionWrappers private: /** - * store a copy of the flags for this solver. + * Store a copy of the flags for this solver. */ const AdditionalData additional_data; }; @@ -371,7 +372,7 @@ namespace ParalutionWrappers private: /** - * store a copy of the flags for this solver. + * Store a copy of the flags for this solver. */ const AdditionalData additional_data; }; @@ -435,7 +436,7 @@ namespace ParalutionWrappers private: /** - * store a copy of the flags for this solver. + * Store a copy of the flags for this solver. */ const AdditionalData additional_data; }; @@ -499,7 +500,134 @@ namespace ParalutionWrappers 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 Ax=b 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 + void solve (SparseMatrix &A, + Vector &x, + Vector &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; }; diff --git a/deal.II/include/deal.II/lac/paralution_vector.h b/deal.II/include/deal.II/lac/paralution_vector.h index c7b0b72926..2d857a849c 100644 --- a/deal.II/include/deal.II/lac/paralution_vector.h +++ b/deal.II/include/deal.II/lac/paralution_vector.h @@ -135,7 +135,6 @@ namespace ParalutionWrappers * 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); /** @@ -444,6 +443,18 @@ namespace ParalutionWrappers */ 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. diff --git a/deal.II/source/lac/paralution_precondition.cc b/deal.II/source/lac/paralution_precondition.cc index 3729452c03..2919a9c9c8 100644 --- a/deal.II/source/lac/paralution_precondition.cc +++ b/deal.II/source/lac/paralution_precondition.cc @@ -219,38 +219,26 @@ namespace ParalutionWrappers /* -------------------------- PreconditionAMG --------------------------- */ - template - PreconditionAMG::AdditionalData::AdditionalData(const unsigned int verbose, - const unsigned int max_iter) - : - default_parameters(true), - verbose(verbose), - max_iter(max_iter) - {} - - - template PreconditionAMG::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) {} @@ -258,9 +246,6 @@ namespace ParalutionWrappers template PreconditionAMG::PreconditionAMG(const AdditionalData &additional_data) : - n_levels(1), - smoothers(NULL), - preconditioners(NULL), coarse_solver(NULL) { this->preconditioner.reset(new paralution::AMG, @@ -284,318 +269,137 @@ namespace ParalutionWrappers { // Downcast the preconditioner pointer paralution::AMG,paralution::LocalVector, - Number>* downcasted_ptr = static_cast,paralution::LocalVector,Number>* >(this->preconditioner.get()); + Number>* downcasted_ptr = static_cast,paralution::LocalVector,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::LocalVector, Number>* [n_levels-1]; - switch (additional_data.smoother) - { + // Coarse solver + switch (additional_data.coarse_solver) + { case richardson : - { - for (unsigned int i=0; i, paralution::LocalVector, - Number> *sm = new paralution::FixedPoint, + { + paralution::FixedPoint, paralution::LocalVector, + Number> *cs = new paralution::FixedPoint, paralution::LocalVector, Number>; - smoothers[i] = sm; - } - break; - } + coarse_solver = cs; + break; + } case cg : - { - for (unsigned int i=0; i, paralution::LocalVector, - Number> *sm = new paralution::CG, + { + paralution::CG, paralution::LocalVector, + Number> *cs = new paralution::CG, paralution::LocalVector, Number>; - smoothers[i] = sm; - } - break; - } + coarse_solver = cs; + break; + } case cr : - { - for (unsigned int i=0; i, paralution::LocalVector, - Number> *sm = new paralution::CR, + { + paralution::CR, paralution::LocalVector, + Number> *cs = new paralution::CR, paralution::LocalVector, Number>; - smoothers[i] = sm; - } - break; - } + coarse_solver = cs; + break; + } case bicgstab : - { - for (unsigned int i=0; i, paralution::LocalVector, - Number> *sm = new paralution::BiCGStab, + { + paralution::BiCGStab, paralution::LocalVector, + Number> *cs = new paralution::BiCGStab, paralution::LocalVector, Number>; - smoothers[i] = sm; - } - break; - } + coarse_solver = cs; + break; + } case gmres : - { - for (unsigned int i=0; i, paralution::LocalVector, - Number> *sm = new paralution::GMRES, + { + paralution::GMRES, paralution::LocalVector, + Number> *cs = new paralution::GMRES, paralution::LocalVector, Number>; - smoothers[i] = sm; - } - break; - } + coarse_solver = cs; + break; + } case fgmres : - { - for (unsigned int i=0; i, paralution::LocalVector, - Number> *sm = new paralution::FGMRES, - paralution::LocalVector, Number>; - smoothers[i] = sm; - } - break; - } - default : - { - AssertThrow(false,ExcMessage("Unknown smoother for PreconditionAMG.")); - } - } - - // Preconditioners for each level. - preconditioners = new paralution::Preconditioner, - paralution::LocalVector, Number>*[n_levels-1]; - switch (additional_data.preconditioner) - { - case jacobi : - { - for (unsigned int i=0; i, paralution::LocalVector, - Number> *precond = new paralution::Jacobi, - paralution::LocalVector, Number>; - preconditioners[i] = precond; - } - break; - } - case sgs : - { - for (unsigned int i=0; i, paralution::LocalVector, - Number> *precond = new paralution::SGS, - paralution::LocalVector, Number>; - preconditioners[i] = precond; - } - break; - } - case multicolored_sgs : - { - for (unsigned int i=0; i, - paralution::LocalVector, Number> *precond = - new paralution::MultiColoredSGS, - paralution::LocalVector, Number>; - preconditioners[i] = precond; - } - break; - } - case multicolored_sor : - { - for (unsigned int i=0; i, - paralution::LocalVector, Number> *precond = - new paralution::MultiColoredGS, - paralution::LocalVector, Number>; - preconditioners[i] = precond; - } - break; - } - case ilu : - { - for (unsigned int i=0; i, paralution::LocalVector, - Number> *precond = new paralution::ILU, - paralution::LocalVector, Number>; - preconditioners[i] = precond; - } - break; - } - case ilut : - { - for (unsigned int i=0; i, paralution::LocalVector, - Number> *precond = new paralution::ILUT, - paralution::LocalVector, Number>; - preconditioners[i] = precond; - } - break; - } - case multicolored_ilu : - { - for (unsigned int i=0; i, - paralution::LocalVector, Number> *precond = - new paralution::MultiColoredILU, + { + paralution::FGMRES, paralution::LocalVector, + Number> *cs = new paralution::FGMRES, paralution::LocalVector, 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; iSetPreconditioner(*preconditioners[i]); - smoothers[i]->Verbose(0); - } - downcasted_ptr->SetSmoother(smoothers); - - // Coarse solver - switch (additional_data.coarse_solver) - { - case richardson : - { - paralution::FixedPoint, paralution::LocalVector, - Number> *cs = new paralution::FixedPoint, - paralution::LocalVector, Number>; - coarse_solver = cs; - break; - } - case cg : - { - paralution::CG, paralution::LocalVector, - Number> *cs = new paralution::CG, - paralution::LocalVector, Number>; - coarse_solver = cs; - break; - } - case cr : - { - paralution::CR, paralution::LocalVector, - Number> *cs = new paralution::CR, - paralution::LocalVector, Number>; - coarse_solver = cs; - break; - } - case bicgstab : - { - paralution::BiCGStab, paralution::LocalVector, - Number> *cs = new paralution::BiCGStab, - paralution::LocalVector, Number>; - coarse_solver = cs; - break; - } - case gmres : - { - paralution::GMRES, paralution::LocalVector, - Number> *cs = new paralution::GMRES, - paralution::LocalVector, Number>; - coarse_solver = cs; - break; - } - case fgmres : - { - paralution::FGMRES, paralution::LocalVector, - Number> *cs = new paralution::FGMRES, - paralution::LocalVector, 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); - } } @@ -609,36 +413,6 @@ namespace ParalutionWrappers delete coarse_solver; coarse_solver = NULL; } - - // Free the preconditioners - for (unsigned int i=0; iexecute_solve(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 + void SolverAMG::solve(SparseMatrix &A, + Vector &x, + Vector &b, + bool move_to_accelerator) + { + paralution::AMG,paralution::LocalVector,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::LocalVector,Number> **smoothers = NULL; + smoothers = new paralution::IterativeLinearSolver, + paralution::LocalVector,Number>* [n_levels-1]; + + // Preconditioner for each smoother + paralution::Preconditioner, + paralution::LocalVector,Number> **preconds = NULL; + preconds = new paralution::Preconditioner, + paralution::LocalVector,Number>* [n_levels-1]; + + // Coarse Grid Solver + paralution::IterativeLinearSolver, + paralution::LocalVector,Number> *coarse_solver; + + // Smoother + switch (additional_data.smoother) + { + case richardson : + { + for (unsigned int i=0; i, + paralution::LocalVector,Number> *richardson = NULL; + richardson = new paralution::FixedPoint, + paralution::LocalVector,Number>; + smoothers[i] = richardson; + smoothers[i]->Verbose(additional_data.verbose); + } + + break; + } + case cg : + { + for (unsigned int i=0; i,paralution::LocalVector, + Number> *cg = NULL; + cg = new paralution::CG,paralution::LocalVector, + Number>; + smoothers[i] = cg; + smoothers[i]->Verbose(additional_data.verbose); + } + + break; + } + case cr : + { + for (unsigned int i=0; i,paralution::LocalVector, + Number> *cr = NULL; + cr = new paralution::CR,paralution::LocalVector, + Number>; + smoothers[i] = cr; + smoothers[i]->Verbose(additional_data.verbose); + } + + break; + } + case bicgstab : + { + for (unsigned int i=0; i,paralution::LocalVector, + Number> *bicgstab = NULL; + bicgstab = new paralution::BiCGStab, + paralution::LocalVector,Number>; + smoothers[i] = bicgstab; + smoothers[i]->Verbose(additional_data.verbose); + } + + break; + } + case gmres : + { + for (unsigned int i=0; i,paralution::LocalVector, + Number> *gmres = NULL; + gmres = new paralution::GMRES, + paralution::LocalVector,Number>; + smoothers[i] = gmres; + smoothers[i]->Verbose(additional_data.verbose); + } + + break; + } + case fgmres : + { + for (unsigned int i=0; i,paralution::LocalVector, + Number> *fgmres = NULL; + fgmres = new paralution::FGMRES, + paralution::LocalVector,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,paralution::LocalVector, + Number> *jacobi = NULL; + jacobi = new paralution::Jacobi, + paralution::LocalVector,Number>; + preconds[i] = jacobi; + smoothers[i]->SetPreconditioner(*preconds[i]); + } + + break; + } + case sgs : + { + for (unsigned int i=0; i,paralution::LocalVector, + Number> *sgs = NULL; + sgs = new paralution::SGS,paralution::LocalVector, + Number>; + preconds[i] = sgs; + smoothers[i]->SetPreconditioner(*preconds[i]); + } + + break; + } + case multicolored_sgs : + { + for (unsigned int i=0; i,paralution::LocalVector, + Number> *multicolored_sgs = NULL; + multicolored_sgs = new paralution::MultiColoredSGS, + paralution::LocalVector,Number>; + multicolored_sgs->SetPrecondMatrixFormat(ELL); + preconds[i] = multicolored_sgs; + smoothers[i]->SetPreconditioner(*preconds[i]); + } + + break; + } + case multicolored_sor : + { + for (unsigned int i=0; i,paralution::LocalVector, + Number> *multicolored_sor = NULL; + multicolored_sor = new paralution::MultiColoredGS, + paralution::LocalVector,Number>; + multicolored_sor->SetPrecondMatrixFormat(ELL); + preconds[i] = multicolored_sor; + smoothers[i]->SetPreconditioner(*preconds[i]); + } + + break; + } + case ilu : + { + for (unsigned int i=0; i,paralution::LocalVector, + Number> *ilu = NULL; + ilu = new paralution::ILU,paralution::LocalVector, + Number>; + preconds[i] = ilu; + smoothers[i]->SetPreconditioner(*preconds[i]); + } + + break; + } + case ilut : + { + for (unsigned int i=0; i,paralution::LocalVector, + Number> *ilut = NULL; + ilut = new paralution::ILUT, + paralution::LocalVector,Number>; + preconds[i] = ilut; + smoothers[i]->SetPreconditioner(*preconds[i]); + } + + break; + } + case multicolored_ilu : + { + for (unsigned int i=0; i,paralution::LocalVector, + Number> *multicolored_ilu = NULL; + multicolored_ilu = new paralution::MultiColoredILU, + paralution::LocalVector,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::LocalVector, + Number> *cs = new paralution::FixedPoint, + paralution::LocalVector, Number>; + coarse_solver = cs; + break; + } + case cg : + { + paralution::CG, paralution::LocalVector, + Number> *cs = new paralution::CG, + paralution::LocalVector, Number>; + coarse_solver = cs; + break; + } + case cr : + { + paralution::CR, paralution::LocalVector, + Number> *cs = new paralution::CR, + paralution::LocalVector, Number>; + coarse_solver = cs; + break; + } + case bicgstab : + { + paralution::BiCGStab, paralution::LocalVector, + Number> *cs = new paralution::BiCGStab, + paralution::LocalVector, Number>; + coarse_solver = cs; + break; + } + case gmres : + { + paralution::GMRES, paralution::LocalVector, + Number> *cs = new paralution::GMRES, + paralution::LocalVector, Number>; + coarse_solver = cs; + break; + } + case fgmres : + { + paralution::FGMRES, paralution::LocalVector, + Number> *cs = new paralution::FGMRES, + paralution::LocalVector, 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 &b, PreconditionBase &preconditioner, bool move_to_accelerator); + + template void SolverAMG::solve(SparseMatrix &A, + Vector &x, + Vector &b, + bool move_to_accelerator); + + template void SolverAMG::solve(SparseMatrix &A, + Vector &x, + Vector &b, + bool move_to_accelerator); } DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/source/lac/paralution_vector.cc b/deal.II/source/lac/paralution_vector.cc index 87b7217ab0..bd14accea9 100644 --- a/deal.II/source/lac/paralution_vector.cc +++ b/deal.II/source/lac/paralution_vector.cc @@ -240,6 +240,36 @@ namespace ParalutionWrappers local_vector.ScaleAdd2(0.,u.paralution_vector(),a,v.paralution_vector(),b); local_vector.AddScale(w.paralution_vector(),c); } + + + + template + void Vector::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(local_vector[i]) << ' '; + else + for (size_type i=0; i(local_vector[i]) << std::endl; + out << std::endl; + + // restore the representation + // of the vector + AssertThrow (out, ExcIO()); + } } // Explicit instantiations