--- /dev/null
+// ---------------------------------------------------------------------
+// $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 <deal.II/base/config.h>
+
+#ifdef DEAL_II_WITH_PARALUTION
+
+#include <deal.II/lac/exceptions.h>
+#include <deal.II/lac/solver_control.h>
+#include <deal.II/lac/paralution_vector.h>
+#include <deal.II/lac/paralution_sparse_matrix.h>
+
+
+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 <tt>Ax=b</tt> using the CG solver of Paralution.
+ */
+ //TODO: add a flag to move to the accelerator
+ template <typename Number>
+ void solve (const SparseMatrix<Number> &A,
+ Vector<Number> &x,
+ const Vector<Number> &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 <tt>Ax=b</tt> using the BiCGStab solver of
+ * Paralution.
+ */
+ template <typename Number>
+ void solve (const SparseMatrix<Number> &A,
+ Vector<Number> &x,
+ const Vector<Number> &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 <tt>Ax=b</tt> using the GMRES solver of
+ * Paralution.
+ */
+ template <typename Number>
+ void solve (const SparseMatrix<Number> &A,
+ Vector<Number> &x,
+ const Vector<Number> &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 ---------------------------*/
#include <algorithm>
#include "paralution.hpp"
-#include "base/local_matrix.hpp"
DEAL_II_NAMESPACE_OPEN
* 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
* @name 3: Modifying entries
*/
//@{
+ /**
+ * Set the element (<i>i,j</i>) to <tt>value</tt>. Throws an error if the
+ * entry does not exist or if <tt>value</tt> 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 <tt>indices</tt>. In other words, this function writes the
+ * elements in <tt>full_matrix</tt> into the calling matrix, using the
+ * local-to-global indexing specified by <tt>indices</tt> 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 <tt>elide_zero_values</tt> 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 <tt>false</tt>, i.e., even zero values are
+ * treated.
+ */
+ template <typename Number2>
+ void set (const std::vector<size_type> &indices,
+ const FullMatrix<Number2> &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 <typename Number2>
+ void set (const std::vector<size_type> &row_indices,
+ const std::vector<size_type> &col_indices,
+ const FullMatrix<Number2> &full_matrix,
+ const bool elide_zero_values = false);
+
+ /**
+ * Set several elements in the specified row of the matrix with column
+ * indices as given by <tt>col_indices</tt> to the respective value.
+ *
+ * The optional parameter <tt>elide_zero_values</tt> 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 <tt>false</tt>, i.e., even zero values are
+ * treated.
+ */
+ template <typename Number2>
+ void set (const size_type row,
+ const std::vector<size_type> &col_indices,
+ const std::vector<Number2> &values,
+ const bool elide_zero_values = false);
+
+ /**
+ * Set several elements to values given by <tt>values</tt> in a given row in
+ * columns given by col_indices into the sparse matrix.
+ *
+ * The optional parameter <tt>elide_zero_values</tt> can be used to specify
+ * whether zero values should be inserted anyway or they should be filtered
+ * away. The default value is <tt>false</tt>, i.e., even zero values are
+ * inserted/replaced.
+ */
+ template <typename Number2>
+ 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 <tt>value</tt> to the element (<i>i,j</i>). Throws an error if
* the entry does not exist or if <tt>value</tt> is not a finite number.
template <typename Number>
inline SparseMatrix<Number>::~SparseMatrix()
{
- local_matrix.clear();
+ local_matrix.Clear();
}
template <typename Number>
inline void SparseMatrix<Number>::reinit(SparsityPattern const &sparsity_pattern)
{
- local_matrix.clear();
+ local_matrix.Clear();
sparse_matrix.reinit(sparsity_pattern);
}
+ template <typename Number>
+ inline void SparseMatrix<Number>::clear()
+ {
+ local_matrix.clear();
+ sparse_matrix.clear();
+ }
+
+
template <typename Number>
inline typename SparseMatrix<Number>::size_type SparseMatrix<Number>::m() const
+ template <typename Number>
+ inline void SparseMatrix<Number>::set(const size_type i,
+ const size_type j,
+ const Number value)
+ {
+ sparse_matrix.set(i,j,value);
+ }
+
+
+
+ template <typename Number>
+ template <typename Number2>
+ inline void SparseMatrix<Number>::set(const std::vector<size_type> &indices,
+ const FullMatrix<Number2> &full_matrix,
+ const bool elide_zero_values)
+ {
+ sparse_matrix.set(indices,full_matrix,elide_zero_values);
+ }
+
+
+
+ template <typename Number>
+ template <typename Number2>
+ inline void SparseMatrix<Number>::set(const size_type row,
+ const std::vector<size_type> &col_indices,
+ const std::vector<Number2> &values,
+ const bool elide_zero_values)
+ {
+ sparse_matrix.set(row,col_indices,values,elide_zero_values);
+ }
+
+
+
+ template <typename Number>
+ template <typename Number2>
+ inline void SparseMatrix<Number>::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 <typename Number>
inline void SparseMatrix<Number>::add(const size_type i,
const size_type j,
#include <deal.II/base/subscriptor.h>
#include "paralution.hpp"
-#include "base/local_vector.hpp"
DEAL_II_NAMESPACE_OPEN
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
--- /dev/null
+// ---------------------------------------------------------------------
+// $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 <deal.II/lac/paralution_solver.h>
+
+#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 <typename Number>
+ void SolverCG::solve(const SparseMatrix<Number> &A,
+ Vector<Number> &x,
+ const Vector<Number> &b)
+ {
+ paralution::CG<paralution::LocalMatrix<Number>,
+ paralution::LocalVector<Number>,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 <typename Number>
+ void SolverBicgstab::solve(const SparseMatrix<Number> &A,
+ Vector<Number> &x,
+ const Vector<Number> &b)
+ {
+ paralution::BiCGStab<paralution::LocalMatrix<Number>,
+ paralution::LocalVector<Number>,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 <typename Number>
+ void SolverGMRES::solve(const SparseMatrix<Number> &A,
+ Vector<Number> &x,
+ const Vector<Number> &b)
+ {
+ paralution::GMRES<paralution::LocalMatrix<Number>,
+ paralution::LocalVector<Number>,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<float>(const SparseMatrix<float> &A,
+ Vector<float> &x,
+ const Vector<float> &b);
+
+ template void SolverCG::solve<double>(const SparseMatrix<double> &A,
+ Vector<double> &x,
+ const Vector<double> &b);
+
+ template void SolverBicgstab::solve<float>(const SparseMatrix<float> &A,
+ Vector<float> &x,
+ const Vector<float> &b);
+
+ template void SolverBicgstab::solve<double>(const SparseMatrix<double> &A,
+ Vector<double> &x,
+ const Vector<double> &b);
+
+ template void SolverGMRES::solve<float>(const SparseMatrix<float> &A,
+ Vector<float> &x,
+ const Vector<float> &b);
+
+ template void SolverGMRES::solve<double>(const SparseMatrix<double> &A,
+ Vector<double> &x,
+ const Vector<double> &b);
+}
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif // DEAL_II_WITH_PARALUTION
--- /dev/null
+// ---------------------------------------------------------------------
+// $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 <cmath>
+#include <fstream>
+#include <iomanip>
+#include <deal.II/base/logstream.h>
+#include <deal.II/lac/paralution_sparse_matrix.h>
+#include <deal.II/lac/paralution_vector.h>
+#include <deal.II/lac/paralution_solver.h>
+
+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<double> A(structure);
+ testproblem.five_point(A);
+ A.convert_to_paralution_csr();
+
+ ParalutionWrappers::Vector<double> u(dim);
+ ParalutionWrappers::Vector<double> f(dim);
+ for (unsigned int i=0; i<dim; ++i)
+ {
+ u[i] = 0.;
+ f[i] = 1.;
+ }
+ solver.solve(A,u,f);
+ }
+
+ paralution::stop_paralution();
+}
--- /dev/null
+
+DEAL::Size 4 Unknowns 9
+DEAL::Size 12 Unknowns 121
--- /dev/null
+// ---------------------------------------------------------------------
+// $Id: paralution_sparse_matrix_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.
+//
+// ---------------------------------------------------------------------
+
+
+// Test ParalutionWrappers::SparseMatrix
+
+#include "../tests.h"
+#include <deal.II/lac/sparsity_pattern.h>
+#include <deal.II/lac/paralution_sparse_matrix.h>
+#include "paralution.hpp"
+
+template <typename Number>
+void check()
+{
+ ParalutionWrappers::SparseMatrix<Number> 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() <<std::endl;
+ deallog << matrix.n() <<std::endl;
+}
+
+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);
+
+ check<float>();
+ check<double>();
+
+ paralution::stop_paralution();
+}
+
--- /dev/null
+
+DEAL::4
+DEAL::5
+DEAL::4
+DEAL::5
projects / dealii-svn.git / commitdiff