From d6023dcf2a88651796593790f58f4b7457f7a736 Mon Sep 17 00:00:00 2001 From: turcksin Date: Mon, 2 Jun 2014 15:10:58 +0000 Subject: [PATCH] Add wrappers for Paralution preconditioners. git-svn-id: https://svn.dealii.org/branches/branch_paralution@33006 0785d39b-7218-0410-832d-ea1e28bc413d --- .../deal.II/lac/paralution_precondition.h | 392 ++++++++++++++++++ .../include/deal.II/lac/paralution_solver.h | 48 ++- .../deal.II/lac/paralution_sparse_matrix.h | 6 +- .../include/deal.II/lac/paralution_vector.h | 2 +- deal.II/source/lac/CMakeLists.txt | 1 + deal.II/source/lac/paralution_precondition.cc | 233 +++++++++++ deal.II/source/lac/paralution_solver.cc | 173 ++++---- .../source/lac/paralution_sparse_matrix.cc | 1 + tests/lac/paralution_solver_01.cc | 6 +- tests/lac/paralution_sparse_matrix_01.cc | 22 +- 10 files changed, 775 insertions(+), 109 deletions(-) create mode 100644 deal.II/include/deal.II/lac/paralution_precondition.h create mode 100644 deal.II/source/lac/paralution_precondition.cc diff --git a/deal.II/include/deal.II/lac/paralution_precondition.h b/deal.II/include/deal.II/lac/paralution_precondition.h new file mode 100644 index 0000000000..5f07ed9128 --- /dev/null +++ b/deal.II/include/deal.II/lac/paralution_precondition.h @@ -0,0 +1,392 @@ +// --------------------------------------------------------------------- +// $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 + +#ifdef DEAL_II_WITH_PARALUTION + +#include + +#include + +DEAL_II_NAMESPACE_OPEN + + +namespace ParalutionWrappers +{ + /** + * 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 + */ + //TODO: preconditioners are missing + template + 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. + */ + std_cxx1x::shared_ptr, + paralution::LocalVector,Number> > preconditioner; + }; + + + + /** + * A wrapper for Jacobi preconditioner for Paralution matrices. + * + * @ingroup ParalutionWrappers + * @ingroup Preconditioners + * @author Bruno Turcksin, 2014 + */ + template + class PreconditionJacobi : public PreconditionBase + { + 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 + class PreconditionSGS : public PreconditionBase + { + 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 + class PreconditionMultiColoredSGS : public PreconditionBase + { + 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 + class PreconditionMultiColoredSOR : public PreconditionBase + { + 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 + class PreconditionILU : public PreconditionBase + { + 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 + class PreconditionILUT : public PreconditionBase + { + 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()); + + private : + /** + * Store a copy of the flags for this particular preconditioner. + */ + AdditionalData additional_data; + }; + + + + /** + * 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 + class PreconditionMultiColoredILU : public PreconditionBase + { + 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()); + + private : + /** + * Store a copy of the flags for this particular preconditioner. + */ + AdditionalData additional_data; + }; +} + + +DEAL_II_NAMESPACE_CLOSE + +#endif // DEAL_II_WITH_PARALUTION + +/*---------------------------- paralution_precondition.h ---------------------------*/ + +#endif +/*---------------------------- paralution_precondition.h ---------------------------*/ diff --git a/deal.II/include/deal.II/lac/paralution_solver.h b/deal.II/include/deal.II/lac/paralution_solver.h index 1e559b2b6b..03ff4a06b2 100644 --- a/deal.II/include/deal.II/lac/paralution_solver.h +++ b/deal.II/include/deal.II/lac/paralution_solver.h @@ -1,7 +1,7 @@ // --------------------------------------------------------------------- // $Id: paralution_solver.h 30040 2013-07-18 17:06:48Z maier $ // -// Copyright (C) 2013 by the deal.II authors +// Copyright (C) 2013, 2014 by the deal.II authors // // This file is part of the deal.II library. // @@ -22,11 +22,14 @@ #ifdef DEAL_II_WITH_PARALUTION +#include #include #include +#include #include #include +#include DEAL_II_NAMESPACE_OPEN @@ -52,6 +55,18 @@ namespace ParalutionWrappers SolverControl &control() const; protected: + /** + * Initialize the solver and solve the system of equations. + */ + template + void execute_solve(std_cxx1x::shared_ptr,paralution::LocalVector,Number> > solver, + const SparseMatrix &A, + Vector &x, + const Vector &b, + const PreconditionBase &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 @@ -83,10 +98,11 @@ namespace ParalutionWrappers * solver is built on the accelerator. */ template - void solve (const SparseMatrix &A, - Vector &x, - const Vector &b, - bool move_to_accelerator=false); + void solve (const SparseMatrix &A, + Vector &x, + const Vector &b, + const PreconditionBase &preconditioner, + bool move_to_accelerator=false); }; @@ -109,10 +125,11 @@ namespace ParalutionWrappers * solver is built on the accelerator. */ template - void solve (const SparseMatrix &A, - Vector &x, - const Vector &b, - bool move_to_accelerator=false); + void solve (const SparseMatrix &A, + Vector &x, + const Vector &b, + const PreconditionBase &preconditioner, + bool move_to_accelerator=false); }; @@ -154,10 +171,11 @@ namespace ParalutionWrappers * solver is built on the accelerator. */ template - void solve (const SparseMatrix &A, - Vector &x, - const Vector &b, - bool move_to_accelerator=false); + void solve (const SparseMatrix &A, + Vector &x, + const Vector &b, + const PreconditionBase &preconditioner, + bool move_to_accelerator=false); private: /** @@ -171,7 +189,7 @@ DEAL_II_NAMESPACE_CLOSE #endif // DEAL_II_WITH_PARALUTION -/*---------------------------- trilinos_solver.h ---------------------------*/ +/*---------------------------- paralution_solver.h ---------------------------*/ #endif -/*---------------------------- trilinos_solver.h ---------------------------*/ +/*---------------------------- paralution_solver.h ---------------------------*/ diff --git a/deal.II/include/deal.II/lac/paralution_sparse_matrix.h b/deal.II/include/deal.II/lac/paralution_sparse_matrix.h index 6eca8f0ae2..20861a3028 100644 --- a/deal.II/include/deal.II/lac/paralution_sparse_matrix.h +++ b/deal.II/include/deal.II/lac/paralution_sparse_matrix.h @@ -26,7 +26,7 @@ #include #include -#include "paralution.hpp" +#include DEAL_II_NAMESPACE_OPEN @@ -46,7 +46,7 @@ namespace ParalutionWrappers * exchangeable. However, Paralution LocalMatix only supports float and double. * * @ingroup ParalutionWrappers - * @ingroup Matrix1 + * @ingroup Matrix * @author Bruno Turcksin, 2013 */ template @@ -135,14 +135,12 @@ namespace ParalutionWrappers * of dimension $m \times n$. This function works only after * convert_to_paralution_csr has been called. */ - //TODO make the function work all the time size_type m() const; /** * Return the dimension of the range space. To remember: the matrix is * of dimension $m \times n$. */ - //TODO make the function work all the time size_type n() const; //@} /** diff --git a/deal.II/include/deal.II/lac/paralution_vector.h b/deal.II/include/deal.II/lac/paralution_vector.h index 366c1ee724..c8136a3866 100644 --- a/deal.II/include/deal.II/lac/paralution_vector.h +++ b/deal.II/include/deal.II/lac/paralution_vector.h @@ -27,7 +27,7 @@ #include #include -#include "paralution.hpp" +#include DEAL_II_NAMESPACE_OPEN diff --git a/deal.II/source/lac/CMakeLists.txt b/deal.II/source/lac/CMakeLists.txt index ddb98aa206..259004d6c2 100644 --- a/deal.II/source/lac/CMakeLists.txt +++ b/deal.II/source/lac/CMakeLists.txt @@ -33,6 +33,7 @@ SET(_src matrix_lib.cc matrix_out.cc parallel_vector.cc + paralution_precondition.cc paralution_sparse_matrix.cc paralution_solver.cc petsc_block_sparse_matrix.cc diff --git a/deal.II/source/lac/paralution_precondition.cc b/deal.II/source/lac/paralution_precondition.cc new file mode 100644 index 0000000000..9502c32452 --- /dev/null +++ b/deal.II/source/lac/paralution_precondition.cc @@ -0,0 +1,233 @@ +// --------------------------------------------------------------------- +// $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 + +#ifdef DEAL_II_WITH_PARALUTION + +DEAL_II_NAMESPACE_OPEN + +namespace ParalutionWrappers +{ + template + PreconditionBase::PreconditionBase() + {} + + + + /* -------------------------- PreconditionJacobi -------------------------- */ + + template + PreconditionJacobi::PreconditionJacobi() + { + this->preconditioner.reset(new paralution::Jacobi, + paralution::LocalVector,Number>); + } + + + + /* -------------------------- PreconditionSGS ----------------------------- */ + + template + PreconditionSGS::PreconditionSGS() + { + this->preconditioner.reset(new paralution::SGS, + paralution::LocalVector,Number>); + } + + + + /* -------------------------- PreconditionMultiColoredSGS ----------------- */ + + template + PreconditionMultiColoredSGS::PreconditionMultiColoredSGS() + { + this->preconditioner.reset(new paralution::MultiColoredSGS, + paralution::LocalVector,Number>); + } + + + + /* -------------------------- PreconditionMulticoloredSOR ----------------- */ + + template + PreconditionMultiColoredSOR::AdditionalData:: + AdditionalData(const double omega) + : + omega(omega) + {} + + + template + PreconditionMultiColoredSOR::PreconditionMultiColoredSOR( + const AdditionalData &additional_data) + : + additional_data(additional_data) + { + this->preconditioner.reset(new paralution::MultiColoredGS, + paralution::LocalVector,Number>); + + initialize(additional_data); + } + + template + void PreconditionMultiColoredSOR::initialize(const AdditionalData &additional_data) + { + // Downcast the preconditioner pointer to use SetRelaxation + paralution::MultiColoredGS,paralution:: + LocalVector,Number>* downcasted_ptr = static_cast,paralution::LocalVector, + Number>* >(this->preconditioner.get()); + + downcasted_ptr->SetRelaxation(additional_data.omega); + } + + + + /* -------------------------- PreconditionILU ----------------------------- */ + + template + PreconditionILU::AdditionalData::AdditionalData(const unsigned int levels) + : + levels(levels) + {} + + template + PreconditionILU::PreconditionILU(const AdditionalData &additional_data) + : + additional_data(additional_data) + { + this->preconditioner.reset(new paralution::ILU, + paralution::LocalVector,Number>); + + initialize(additional_data); + } + + template + void PreconditionILU::initialize(const AdditionalData &additional_data) + { + // Downcast the preconditioner pointer to use Set + paralution::ILU,paralution::LocalVector, + Number>* downcasted_ptr = static_cast,paralution::LocalVector,Number>* >(this->preconditioner.get()); + + downcasted_ptr->Set(additional_data.levels); + } + + + + /* -------------------------- PreconditionILUT ---------------------------- */ + + template + PreconditionILUT::AdditionalData::AdditionalData(const Number threshold) + : + threshold(threshold), + max_row(0) + {} + + template + PreconditionILUT::AdditionalData::AdditionalData(const Number threshold, + const unsigned int max_row) + : + threshold(threshold), + max_row(max_row) + {} + + template + PreconditionILUT::PreconditionILUT(const AdditionalData &additional_data) + : + additional_data(additional_data) + { + this->preconditioner.reset(new paralution::ILUT, + paralution::LocalVector,Number>); + + initialize(additional_data); + } + + template + void PreconditionILUT::initialize(const AdditionalData &additional_data) + { + // Downcast the preconditioner pointer to use Set + paralution::ILUT,paralution::LocalVector, + Number>* downcasted_ptr = static_cast,paralution::LocalVector,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 + PreconditionMultiColoredILU::AdditionalData:: + AdditionalData(const unsigned int levels, const unsigned int power) + : + levels(levels), + power(power) + {} + + template + PreconditionMultiColoredILU:: + PreconditionMultiColoredILU(const AdditionalData &additional_data) + : + additional_data(additional_data) + { + this->preconditioner.reset(new paralution::MultiColoredILU, + paralution::LocalVector,Number>); + + initialize(additional_data); + } + + template + void PreconditionMultiColoredILU::initialize(const AdditionalData &additional_data) + { + // Downcast the preconditioner pointer to use Set + paralution::MultiColoredILU,paralution::LocalVector, + Number>* downcasted_ptr = static_cast,paralution::LocalVector,Number>* >(this->preconditioner.get()); + + downcasted_ptr->Set(additional_data.levels,additional_data.power); + } +} + +// Explicit instantiations +namespace ParalutionWrappers +{ + template class PreconditionBase; + template class PreconditionBase; + template class PreconditionJacobi; + template class PreconditionJacobi; + template class PreconditionSGS; + template class PreconditionSGS; + template class PreconditionMultiColoredSGS; + template class PreconditionMultiColoredSGS; + template class PreconditionMultiColoredSOR; + template class PreconditionMultiColoredSOR; + template class PreconditionILU; + template class PreconditionILU; + template class PreconditionILUT; + template class PreconditionILUT; + template class PreconditionMultiColoredILU; + template class PreconditionMultiColoredILU; +} + +DEAL_II_NAMESPACE_CLOSE + +#endif // DEAL_II_WITH_PARALUTION diff --git a/deal.II/source/lac/paralution_solver.cc b/deal.II/source/lac/paralution_solver.cc index 3bf9e3668b..59c999830b 100644 --- a/deal.II/source/lac/paralution_solver.cc +++ b/deal.II/source/lac/paralution_solver.cc @@ -1,7 +1,7 @@ // --------------------------------------------------------------------- // $Id: paralution_solver.cc 31349 2013-10-20 19:07:06Z maier $ // -// Copyright (C) 2013 by the deal.II authors +// Copyright (C) 2013, 2014 by the deal.II authors // // This file is part of the deal.II library. // @@ -18,8 +18,6 @@ #ifdef DEAL_II_WITH_PARALUTION -#include "paralution.hpp" - DEAL_II_NAMESPACE_OPEN namespace ParalutionWrappers @@ -36,6 +34,34 @@ namespace ParalutionWrappers return solver_control; } + template + void SolverBase::execute_solve(std_cxx1x::shared_ptr,paralution::LocalVector,Number> > solver, + const SparseMatrix &A, + Vector &x, + const Vector &b, + const PreconditionBase &preconditioner, + bool move_to_accelerator) + { + // Set the preconditioner. + solver->SetPreconditioner(*(preconditioner.preconditioner)); + + // 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()); + + // Move the solver to the accelerator if necessary. + if (move_to_accelerator==true) + solver->MoveToAccelerator(); + + solver->Build(); + solver->Solve(b.paralution_vector(),&(x.paralution_vector())); + } + /* ---------------------- SolverCG ------------------------ */ @@ -48,22 +74,17 @@ namespace ParalutionWrappers template - void SolverCG::solve(const SparseMatrix &A, - Vector &x, - const Vector &b, - bool move_to_accelerator) + void SolverCG::solve(const SparseMatrix &A, + Vector &x, + const Vector &b, + const PreconditionBase &preconditioner, + bool move_to_accelerator) { - paralution::CG, - paralution::LocalVector,Number> solver; - 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()); - if (move_to_accelerator==true) - solver.MoveToAccelerator(); - solver.Build(); - solver.Solve(b.paralution_vector(),&(x.paralution_vector())); + std_cxx1x::shared_ptr, + paralution::LocalVector,Number> > solver(new paralution:: + CG,paralution::LocalVector,Number>); + + this->execute_solve(solver,A,x,b,preconditioner,move_to_accelerator); } @@ -78,21 +99,17 @@ namespace ParalutionWrappers template - void SolverBicgstab::solve(const SparseMatrix &A, - Vector &x, - const Vector &b, - bool move_to_accelerator) + void SolverBicgstab::solve(const SparseMatrix &A, + Vector &x, + const Vector &b, + const PreconditionBase &preconditioner, + bool move_to_accelerator) { - paralution::BiCGStab, - paralution::LocalVector,Number> solver; - // Set absolute tolerance, relative tolerance, divergence tolerance, - // maximum number of iterations. - solver.Init(solver_control.tolerance(),0.,1.e100, - solver_control.max_steps()); - if (move_to_accelerator==true) - solver.MoveToAccelerator(); - solver.Build(); - solver.Solve(b.paralution_vector(),&(x.paralution_vector())); + std_cxx1x::shared_ptr, + paralution::LocalVector,Number> > solver(new paralution:: + BiCGStab,paralution::LocalVector,Number>); + + this->execute_solve(solver,A,x,b,preconditioner,move_to_accelerator); } @@ -109,25 +126,21 @@ namespace ParalutionWrappers template - void SolverGMRES::solve(const SparseMatrix &A, - Vector &x, - const Vector &b, - bool move_to_accelerator) + void SolverGMRES::solve(const SparseMatrix &A, + Vector &x, + const Vector &b, + const PreconditionBase &preconditioner, + bool move_to_accelerator) { - paralution::GMRES, - paralution::LocalVector,Number> solver; - // Set absolute tolerance, relative tolerance, divergence tolerance, - // maximum number of iterations. - solver.Init(solver_control.tolerance(),0.,1.e100, - solver_control.max_steps()); - solver.SetBasisSize(additional_data.restart_parameter); - if (move_to_accelerator==true) - solver.MoveToAccelerator(); - solver.Build(); - solver.Solve(b.paralution_vector(),&(x.paralution_vector())); - } + std_cxx1x::shared_ptr, + paralution::LocalVector,Number> > solver(new paralution:: + GMRES, paralution::LocalVector,Number>); + // Set the restart parameter. + solver->SetBasisSize(additional_data.restart_parameter); + this->execute_solve(solver,A,x,b,preconditioner,move_to_accelerator); + } } @@ -135,35 +148,41 @@ namespace ParalutionWrappers // Explicit instantiations namespace ParalutionWrappers { - template void SolverCG::solve(const SparseMatrix &A, - Vector &x, - const Vector &b, - bool move_to_accelerator); - - template void SolverCG::solve(const SparseMatrix &A, - Vector &x, - const Vector &b, - bool move_to_accelerator); - - template void SolverBicgstab::solve(const SparseMatrix &A, - Vector &x, - const Vector &b, - bool move_to_accelerator); - - template void SolverBicgstab::solve(const SparseMatrix &A, - Vector &x, - const Vector &b, - bool move_to_accelerator); - - template void SolverGMRES::solve(const SparseMatrix &A, - Vector &x, - const Vector &b, - bool move_to_accelerator); - - template void SolverGMRES::solve(const SparseMatrix &A, - Vector &x, - const Vector &b, - bool move_to_accelerator); + template void SolverCG::solve(const SparseMatrix &A, + Vector &x, + const Vector &b, + const PreconditionBase &preconditioner, + bool move_to_accelerator); + + template void SolverCG::solve(const SparseMatrix &A, + Vector &x, + const Vector &b, + const PreconditionBase &preconditioner, + bool move_to_accelerator); + + template void SolverBicgstab::solve(const SparseMatrix &A, + Vector &x, + const Vector &b, + const PreconditionBase &preconditioner, + bool move_to_accelerator); + + template void SolverBicgstab::solve(const SparseMatrix &A, + Vector &x, + const Vector &b, + const PreconditionBase &preconditioner, + bool move_to_accelerator); + + template void SolverGMRES::solve(const SparseMatrix &A, + Vector &x, + const Vector &b, + const PreconditionBase &preconditioner, + bool move_to_accelerator); + + template void SolverGMRES::solve(const SparseMatrix &A, + Vector &x, + const Vector &b, + const PreconditionBase &preconditioner, + bool move_to_accelerator); } DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/source/lac/paralution_sparse_matrix.cc b/deal.II/source/lac/paralution_sparse_matrix.cc index 46cadb95d2..38d924f3d2 100644 --- a/deal.II/source/lac/paralution_sparse_matrix.cc +++ b/deal.II/source/lac/paralution_sparse_matrix.cc @@ -76,6 +76,7 @@ namespace ParalutionWrappers // Free the memory used by sparse_matrix. sparse_matrix.clear(); + is_local_matrix = true; } } diff --git a/tests/lac/paralution_solver_01.cc b/tests/lac/paralution_solver_01.cc index fd730e1642..720f95e074 100644 --- a/tests/lac/paralution_solver_01.cc +++ b/tests/lac/paralution_solver_01.cc @@ -26,6 +26,7 @@ #include #include #include +#include int main() { @@ -56,6 +57,9 @@ int main() testproblem.five_point(A); A.convert_to_paralution_csr(); + //Preconditioner + ParalutionWrappers::PreconditionJacobi p; + ParalutionWrappers::Vector u(dim); ParalutionWrappers::Vector f(dim); for (unsigned int i=0; i void check() { - ParalutionWrappers::SparseMatrix matrix; - SparsityPattern pattern(4,5,2); - pattern.add(0,2); - pattern.add(0,0); - pattern.add(1,0); - pattern.add(1,3); - pattern.add(2,4); - pattern.add(2,2); - pattern.add(3,0); - pattern.add(3,4); - pattern.compress(); + SparsityPattern sparsity_pattern(4,5,2); + sparsity_pattern.add(0,2); + sparsity_pattern.add(0,0); + sparsity_pattern.add(1,0); + sparsity_pattern.add(1,3); + sparsity_pattern.add(2,4); + sparsity_pattern.add(2,2); + sparsity_pattern.add(3,0); + sparsity_pattern.add(3,4); + sparsity_pattern.compress(); + ParalutionWrappers::SparseMatrix matrix; matrix.reinit(sparsity_pattern); matrix.add(0,2,3.5); matrix.add(0,0,1.); -- 2.39.5