From 693d0057f9e65b176c4cfbd1bd419dca671dcdb6 Mon Sep 17 00:00:00 2001 From: turcksin Date: Fri, 13 Dec 2013 20:32:36 +0000 Subject: [PATCH] Add some of Paralutions solvers and add tests. git-svn-id: https://svn.dealii.org/branches/branch_paralution@31992 0785d39b-7218-0410-832d-ea1e28bc413d --- .../include/deal.II/lac/paralution_solver.h | 171 ++++++++++++++++++ .../deal.II/lac/paralution_sparse_matrix.h | 134 +++++++++++++- .../include/deal.II/lac/paralution_vector.h | 1 - deal.II/source/lac/CMakeLists.txt | 1 + deal.II/source/lac/paralution_solver.cc | 156 ++++++++++++++++ tests/lac/paralution_solver_01.cc | 70 +++++++ ...ution_solver_01.with_paralution=yes.output | 3 + tests/lac/paralution_sparse_matrix_01.cc | 72 ++++++++ ...parse_matrix_01.with_paralution=yes.output | 5 + 9 files changed, 608 insertions(+), 5 deletions(-) create mode 100644 deal.II/include/deal.II/lac/paralution_solver.h create mode 100644 deal.II/source/lac/paralution_solver.cc create mode 100644 tests/lac/paralution_solver_01.cc create mode 100644 tests/lac/paralution_solver_01.with_paralution=yes.output create mode 100644 tests/lac/paralution_sparse_matrix_01.cc create mode 100644 tests/lac/paralution_sparse_matrix_01.with_paralution=yes.output diff --git a/deal.II/include/deal.II/lac/paralution_solver.h b/deal.II/include/deal.II/lac/paralution_solver.h new file mode 100644 index 0000000000..ff52c026ac --- /dev/null +++ b/deal.II/include/deal.II/lac/paralution_solver.h @@ -0,0 +1,171 @@ +// --------------------------------------------------------------------- +// $Id: paralution_solver.h 30040 2013-07-18 17:06:48Z maier $ +// +// Copyright (C) 2013 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +#ifndef __deal2__paralution_solver_h +#define __deal2__paralution_solver_h + + +#include + +#ifdef DEAL_II_WITH_PARALUTION + +#include +#include +#include +#include + + +DEAL_II_NAMESPACE_OPEN + +namespace ParalutionWrappers +{ + /** + * Base class for solver classes using the Paralution solvers. + * + * @ingroup ParalutionWrappers + * @author Bruno Turcksin 2013 + */ + class SolverBase + { + public: + /** + * Constructor. Takes the solver control object and creates the solver. + */ + SolverBase (SolverControl &cn); + + /** + * Access to object that controls convergence. + */ + SolverControl &control() const; + + protected: + /** + * Reference to the object that controls convergence of the iteratove + * solver. In fact, for these Paralution wrappers, Paralution does so + * itself, but we copy the data from this object solition process, and + * copy the data back into it afterwards. + */ + SolverControl &solver_control; + }; + + + + /** + * An implementation of the solver interface using the Paralution CG solver. + * + * @ingroup ParalutionWrappers + * @author Bruno Turcksin, 2013 + */ + class SolverCG : public SolverBase + { + public: + /** + * Constructor. + */ + SolverCG (SolverControl &cn); + + /** + * Solve the linear system Ax=b using the CG solver of Paralution. + */ + //TODO: add a flag to move to the accelerator + template + void solve (const SparseMatrix &A, + Vector &x, + const Vector &b); + }; + + + + /** + * An implementation of the solver interface using the Paralution BiCGStab + * solver. + */ + class SolverBicgstab : public SolverBase + { + public: + /** + * Constructor. + */ + SolverBicgstab (SolverControl &cn); + + /** + * Solve the linear system Ax=b using the BiCGStab solver of + * Paralution. + */ + template + void solve (const SparseMatrix &A, + Vector &x, + const Vector &b); + }; + + + + /** + * An implementation of the solver interface using the Paralution GMRES + * solver. + */ + class SolverGMRES : public SolverBase + { + public: + /** + * Standardized data struct to pipe additional data to the solver. + */ + struct AdditionalData + { + /** + * Constructor. By default, set the size of the Krylov space to 30, + * i.e. do a restart every 30 iterations. + */ + AdditionalData (const unsigned int restart_parameter = 30); + + /** + * Size of the Krylov space. + */ + unsigned int restart_parameter; + }; + + /** + * Constructor. AdditionalData is a structure that contains additional + * flags for tuning this particulat solver. + */ + SolverGMRES (SolverControl &cn, + const AdditionalData &data = AdditionalData()); + + /** + * Solve the linear system Ax=b using the GMRES solver of + * Paralution. + */ + template + void solve (const SparseMatrix &A, + Vector &x, + const Vector &b); + + private: + /** + * Store a copy of the flags for this solver. + */ + const AdditionalData additional_data; + }; +} + +DEAL_II_NAMESPACE_CLOSE + +#endif // DEAL_II_WITH_PARALUTION + +/*---------------------------- trilinos_solver.h ---------------------------*/ + +#endif +/*---------------------------- trilinos_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 a856acdaee..217d0a90a5 100644 --- a/deal.II/include/deal.II/lac/paralution_sparse_matrix.h +++ b/deal.II/include/deal.II/lac/paralution_sparse_matrix.h @@ -27,7 +27,6 @@ #include #include "paralution.hpp" -#include "base/local_matrix.hpp" DEAL_II_NAMESPACE_OPEN @@ -118,7 +117,7 @@ namespace ParalutionWrappers * called the default constructor. It also forgets the sparsity pattern * it was previously tied to. */ - virtual void clear(); + void clear(); /** * This function convert the underlying SparseMatrix to @@ -149,6 +148,79 @@ namespace ParalutionWrappers * @name 3: Modifying entries */ //@{ + /** + * Set the element (i,j) to value. Throws an error if the + * entry does not exist or if value is not a finite number. Still, + * it is allowed to store zero values in non-existent fields. + */ + void set (const size_type i, + const size_type j, + const Number value); + + /** + * Set all elements given in a FullMatrix into the sparse matrix locations + * given by indices. In other words, this function writes the + * elements in full_matrix into the calling matrix, using the + * local-to-global indexing specified by indices for both the rows + * and the columns of the matrix. This function assumes a quadratic sparse + * matrix and a quadratic full_matrix, the usual situation in FE + * calculations. + * + * The optional parameter elide_zero_values can be used to specify + * whether zero values should be set anyway or they should be filtered away + * (and not change the previous content in the respective element if it + * exists). The default value is false, i.e., even zero values are + * treated. + */ + template + void set (const std::vector &indices, + const FullMatrix &full_matrix, + const bool elide_zero_values = false); + + /** + * Same function as before, but now including the possibility to use + * rectangular full_matrices and different local-to-global indexing on rows + * and columns, respectively. + */ + template + void set (const std::vector &row_indices, + const std::vector &col_indices, + const FullMatrix &full_matrix, + const bool elide_zero_values = false); + + /** + * Set several elements in the specified row of the matrix with column + * indices as given by col_indices to the respective value. + * + * The optional parameter elide_zero_values can be used to specify + * whether zero values should be set anyway or they should be filtered away + * (and not change the previous content in the respective element if it + * exists). The default value is false, i.e., even zero values are + * treated. + */ + template + void set (const size_type row, + const std::vector &col_indices, + const std::vector &values, + const bool elide_zero_values = false); + + /** + * Set several elements to values given by values in a given row in + * columns given by col_indices into the sparse matrix. + * + * The optional parameter elide_zero_values can be used to specify + * whether zero values should be inserted anyway or they should be filtered + * away. The default value is false, i.e., even zero values are + * inserted/replaced. + */ + template + void set (const size_type row, + const size_type n_cols, + const size_type *col_indices, + const Number2 *values, + const bool elide_zero_values = false); + + /** * Add value to the element (i,j). Throws an error if * the entry does not exist or if value is not a finite number. @@ -278,7 +350,7 @@ namespace ParalutionWrappers template inline SparseMatrix::~SparseMatrix() { - local_matrix.clear(); + local_matrix.Clear(); } @@ -286,11 +358,19 @@ namespace ParalutionWrappers template inline void SparseMatrix::reinit(SparsityPattern const &sparsity_pattern) { - local_matrix.clear(); + local_matrix.Clear(); sparse_matrix.reinit(sparsity_pattern); } + template + inline void SparseMatrix::clear() + { + local_matrix.clear(); + sparse_matrix.clear(); + } + + template inline typename SparseMatrix::size_type SparseMatrix::m() const @@ -308,6 +388,52 @@ namespace ParalutionWrappers + template + inline void SparseMatrix::set(const size_type i, + const size_type j, + const Number value) + { + sparse_matrix.set(i,j,value); + } + + + + template + template + inline void SparseMatrix::set(const std::vector &indices, + const FullMatrix &full_matrix, + const bool elide_zero_values) + { + sparse_matrix.set(indices,full_matrix,elide_zero_values); + } + + + + template + template + inline void SparseMatrix::set(const size_type row, + const std::vector &col_indices, + const std::vector &values, + const bool elide_zero_values) + { + sparse_matrix.set(row,col_indices,values,elide_zero_values); + } + + + + template + template + inline void SparseMatrix::set(const size_type row, + const size_type n_cols, + const size_type *col_indices, + const Number2 *values, + const bool elide_zero_values) + { + sparse_matrix.set(row,n_cols,col_indices,values,elide_zero_values); + } + + + template inline void SparseMatrix::add(const size_type i, const size_type j, diff --git a/deal.II/include/deal.II/lac/paralution_vector.h b/deal.II/include/deal.II/lac/paralution_vector.h index 70b3370f9e..00bee9a58e 100644 --- a/deal.II/include/deal.II/lac/paralution_vector.h +++ b/deal.II/include/deal.II/lac/paralution_vector.h @@ -25,7 +25,6 @@ #include #include "paralution.hpp" -#include "base/local_vector.hpp" DEAL_II_NAMESPACE_OPEN diff --git a/deal.II/source/lac/CMakeLists.txt b/deal.II/source/lac/CMakeLists.txt index f5dddee3a7..ddb98aa206 100644 --- a/deal.II/source/lac/CMakeLists.txt +++ b/deal.II/source/lac/CMakeLists.txt @@ -34,6 +34,7 @@ SET(_src matrix_out.cc parallel_vector.cc paralution_sparse_matrix.cc + paralution_solver.cc petsc_block_sparse_matrix.cc petsc_full_matrix.cc petsc_matrix_base.cc diff --git a/deal.II/source/lac/paralution_solver.cc b/deal.II/source/lac/paralution_solver.cc new file mode 100644 index 0000000000..0a7e6564e8 --- /dev/null +++ b/deal.II/source/lac/paralution_solver.cc @@ -0,0 +1,156 @@ +// --------------------------------------------------------------------- +// $Id: paralution_solver.cc 31349 2013-10-20 19:07:06Z maier $ +// +// Copyright (C) 2013 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +#include + +#ifdef DEAL_II_WITH_PARALUTION + +#include "paralution.hpp" + +DEAL_II_NAMESPACE_OPEN + +namespace ParalutionWrappers +{ + SolverBase::SolverBase (SolverControl &cn) + : + solver_control (cn) + {} + + + + SolverControl &SolverBase::control() const + { + return solver_control; + } + + + + /* ---------------------- SolverCG ------------------------ */ + + SolverCG::SolverCG (SolverControl &cn) + : + SolverBase (cn) + {} + + + + template + void SolverCG::solve(const SparseMatrix &A, + Vector &x, + const Vector &b) + { + 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(),1e10,1e10, + solver_control.max_steps()); + solver.Build(); + solver.Solve(b.paralution_vector(),&(x.paralution_vector())); + } + + + + /* -------------------- SolverBicgstab --------------------- */ + + SolverBicgstab::SolverBicgstab (SolverControl &cn) + : + SolverBase (cn) + {} + + + + template + void SolverBicgstab::solve(const SparseMatrix &A, + Vector &x, + const Vector &b) + { + paralution::BiCGStab, + paralution::LocalVector,Number> solver; + // Set absolute tolerance, relative tolerance, divergence tolerance, + // maximum number of iterations. + solver.Init(solver_control.tolerance(),1e10,1e10, + solver_control.max_steps()); + solver.Build(); + solver.Solve(b.paralution_vector(),&(x.paralution_vector())); + } + + + + /* --------------------- SolverGMRES ----------------------- */ + + SolverGMRES::SolverGMRES (SolverControl &cn, + const AdditionalData &data) + : + SolverBase (cn), + additional_data (data) + {} + + + + template + void SolverGMRES::solve(const SparseMatrix &A, + Vector &x, + const Vector &b) + { + paralution::GMRES, + paralution::LocalVector,Number> solver; + // Set absolute tolerance, relative tolerance, divergence tolerance, + // maximum number of iterations. + solver.Init(solver_control.tolerance(),1e10,1e10, + solver_control.max_steps()); + solver.SetBasisSize(additional_data.restart_parameter); + solver.Build(); + solver.Solve(b.paralution_vector(),&(x.paralution_vector())); + } + + +} + + + +// Explicit instantiations +namespace ParalutionWrappers +{ + template void SolverCG::solve(const SparseMatrix &A, + Vector &x, + const Vector &b); + + template void SolverCG::solve(const SparseMatrix &A, + Vector &x, + const Vector &b); + + template void SolverBicgstab::solve(const SparseMatrix &A, + Vector &x, + const Vector &b); + + template void SolverBicgstab::solve(const SparseMatrix &A, + Vector &x, + const Vector &b); + + template void SolverGMRES::solve(const SparseMatrix &A, + Vector &x, + const Vector &b); + + template void SolverGMRES::solve(const SparseMatrix &A, + Vector &x, + const Vector &b); +} + +DEAL_II_NAMESPACE_CLOSE + +#endif // DEAL_II_WITH_PARALUTION diff --git a/tests/lac/paralution_solver_01.cc b/tests/lac/paralution_solver_01.cc new file mode 100644 index 0000000000..fd730e1642 --- /dev/null +++ b/tests/lac/paralution_solver_01.cc @@ -0,0 +1,70 @@ +// --------------------------------------------------------------------- +// $Id: paralution_solver_01.cc 31349 2013-10-20 19:07:06Z maier $ +// +// Copyright (C) 2013 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Create a test for ParalutionWrappers::SolverCG + +#include "../tests.h" +#include "testmatrix.h" +#include +#include +#include +#include +#include +#include +#include + +int main() +{ + paralution::init_paralution(); + + std::ofstream logfile("output"); + deallog << std::fixed; + deallog << std::setprecision(0); + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + SolverControl control(100,1.e-3); + ParalutionWrappers::SolverCG solver(control); + + for (unsigned int size=4; size<=30; size*=3) + { + unsigned int dim = (size-1)*(size-1); + + deallog << "Size "<< size <<" Unknowns "<< dim << std::endl; + + // Make matrix + FDMatrix testproblem(size, size); + SparsityPattern structure(dim, dim, 5); + testproblem.five_point_structure(structure); + structure.compress(); + ParalutionWrappers::SparseMatrix A(structure); + testproblem.five_point(A); + A.convert_to_paralution_csr(); + + ParalutionWrappers::Vector u(dim); + ParalutionWrappers::Vector f(dim); + for (unsigned int i=0; i +#include +#include "paralution.hpp" + +template +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(); + + matrix.reinit(sparsity_pattern); + matrix.add(0,2,3.5); + matrix.add(0,0,1.); + matrix.add(1,0,-2.); + matrix.add(1,3,1.5); + matrix.add(2,4,-2.25); + matrix.add(2,2,-0.5); + matrix.add(3,0,2.); + matrix.add(3,4,0.); + + matrix.convert_to_paralution_csr(); + + deallog << matrix.m() <(); + check(); + + paralution::stop_paralution(); +} + diff --git a/tests/lac/paralution_sparse_matrix_01.with_paralution=yes.output b/tests/lac/paralution_sparse_matrix_01.with_paralution=yes.output new file mode 100644 index 0000000000..6437df48b4 --- /dev/null +++ b/tests/lac/paralution_sparse_matrix_01.with_paralution=yes.output @@ -0,0 +1,5 @@ + +DEAL::4 +DEAL::5 +DEAL::4 +DEAL::5 -- 2.39.5