From 064b4c034a4bb883c2532e9bbc6e7b4d0be872c2 Mon Sep 17 00:00:00 2001 From: turcksin Date: Thu, 5 Dec 2013 22:22:18 +0000 Subject: [PATCH] Start ParalutionWrappers for sparse matrix. git-svn-id: https://svn.dealii.org/branches/branch_paralution@31902 0785d39b-7218-0410-832d-ea1e28bc413d --- .../deal.II/lac/paralution_sparse_matrix.h | 381 ++++++++++++++++++ .../include/deal.II/lac/paralution_vector.h | 2 +- deal.II/include/deal.II/lac/sparse_matrix.h | 2 +- .../deal.II/lac/trilinos_sparse_matrix.h | 2 +- deal.II/source/lac/CMakeLists.txt | 1 + .../source/lac/paralution_sparse_matrix.cc | 92 +++++ 6 files changed, 477 insertions(+), 3 deletions(-) create mode 100644 deal.II/include/deal.II/lac/paralution_sparse_matrix.h create mode 100644 deal.II/source/lac/paralution_sparse_matrix.cc diff --git a/deal.II/include/deal.II/lac/paralution_sparse_matrix.h b/deal.II/include/deal.II/lac/paralution_sparse_matrix.h new file mode 100644 index 0000000000..a856acdaee --- /dev/null +++ b/deal.II/include/deal.II/lac/paralution_sparse_matrix.h @@ -0,0 +1,381 @@ +// --------------------------------------------------------------------- +// $Id: paralution_sparse_matrix.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_matrix_h +#define __deal2__paralution_matrix_h + +#include + +#ifdef DEAL_II_WITH_PARALUTION + +#include +#include +#include +#include + +#include "paralution.hpp" +#include "base/local_matrix.hpp" + +DEAL_II_NAMESPACE_OPEN + +/** + * @addtogroup ParalutionWrappers + *@{ + */ +namespace ParalutionWrappers +{ + /** + * This class implements a wrapper to use the Paralution sparse matrix class + * LocalMatrix. This class is designed for use in either serial + * implementation or as a localized copy on each processors. + * + * The interface of this class is modeled after the existing SparseMatrix + * class in deal.II. It has almost the same member functions, and is often + * exchangeable. However, Paralution LocalMatix only supports float and double. + * + * @ingroup ParalutionWrappers + * @ingroup Matrix1 + * @author Bruno Turcksin, 2013 + */ + template + class SparseMatrix : public Subscriptor + { + public : + /** + * Declare type for container size. + */ + typedef types::global_dof_index size_type; + + /** + * @name 1: Constructors and initialization + */ + //@{ + /** + * Constructor; initializes the matrix to be empty, without any + * structure, i.e. the matrix is not usable at all. This constructor is + * therefore only useful for matrices which are members of a class. All + * other matrices should be created at a point in the data flow where + * all necessary information is available. + * + * You have to initialize the matrix before usage with reinit(const + * SparsityPattern &sparsity_pattern). + */ + SparseMatrix(); + + /** + * Copy constructor. This constructor is only allowed to be called if + * the matrix to be copied is empty. This is for the same reason as for + * the SparsityPattern, see there fore the details. + * + * If you really want to copy a whole matrix, you can do so by using the + * copy_from() function. + */ + //TODO + //SparseMarix(const SparseMatrix &sparse_matrix); + + /** + * Constructo. Takes the given matrix sparsity structure to represent + * the sparsity pattern of this matrix. You can change the sparsity + * pattern later on by calling the reinit(const SparsityPattern&) + * function. + * + * The constructor is marked explicit so as to disallow that someone + * passes a sparsity pattern in place of a spaese matrix to some + * function, where an empty matrix would be generated then. + */ + explicit SparseMatrix(const SparsityPattern &sparsity_pattern); + + /** + * Destructor. Free all memory, but do not release the memory of the + * sparsity structure. + */ + ~SparseMatrix(); + + /** + * Reinitialize the sparse matrix with the given sparsity pattern. The + * latter tells the matrix how many nonzero elements there need to be + * reserved. + * + * The elements of the matrix are set to zero by this function. + */ + void reinit(SparsityPattern const &sparsity_pattern); + + /** + * Release all memory and return to a state just like after having + * called the default constructor. It also forgets the sparsity pattern + * it was previously tied to. + */ + virtual void clear(); + + /** + * This function convert the underlying SparseMatrix to + * Paralution::LocalMatrix. This function frees the SparseMatrix. + */ + void convert_to_paralution_csr(); + //@} + /** + * @name 2: Information on the matrix + */ + //@{ + /** + * Return the dimension of the image space. To remember: the matrix is + * 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; + //@} + /** + * @name 3: Modifying entries + */ + //@{ + /** + * Add value to the element (i,j). Throws an error if + * the entry does not exist or if value is not a finite number. + * Still it is allowed to store zerp values in non-existent fields. + * + * This function only works on the host. + */ + void add(const size_type i, + const size_type j, + const Number value); + + /** + * Add all elements given in a FullMatrix into sparse matrix locations + * given by indices. In other words, this function adds the + * elements in full_matrix to the respective entries in 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 added anyway or these should be + * filtered away and only non-zero data is added. The defaul value is + * true, i.e., zero values won't be added into the matrix. + * + * This function only works on the host. + */ + template + void add(const std::vector &indices, + const FullMatrix &full_matrix, + const bool elide_zero_values = true); + + /** + * Same function as before, but now including the possibility to use + * rectaangular full_matrices and different local-to-global indexing on + * rows and columns, respectively. + * + * This function only works on the host. + */ + template + void add (const std::vector &row_indices, + const std::vector &col_indices, + const FullMatrix &full_matrix, + const bool elide_zero_values = true); + + /** + * 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 added anyway or these should be + * filtered away and only non-zero data is added. The default value is + * true, i.e., zero values won't be added into the matrix. + * + * This function only works on the host. + */ + template + void add (const size_type row, + const std::vector &col_indices, + const std::vector &values, + const bool elide_zero_values = true); + + /** + * Add an array of values given by values in the given global + * matrix row at columns specified by col_indices in the sparse_matrix. + * + * The optional parameter elide_zero_values can be used to + * specify whether zero values should be added anyway or these should be + * filtered away and only non-zero data is added. The default value is + * true>, i.e., zero values won't be added into the matrix. + * + * This function only works on the host. + */ + template + void add (const size_type row, + const size_type n_cols, + const size_type *col_indices, + const Number2 *values, + const bool elide_zero_values = true, + const bool col_indices_are_sorted = false); + + //@} + /** + * @name 4: Access to underlying Paralution data + */ + /** + * Return a constant reference to the underlying Paralution LocalMatrix + * data. + */ + paralution::LocalMatrix const ¶lution_matrix() const; + + /** + * Return a reference to the underlying Paralution LocalMatrix data. + */ + paralution::LocalMatrix& paralution_matrix(); + //@} + + private : + /** + * Underlying Paralution LocalMatrix. + */ + paralution::LocalMatrix local_matrix; + + /** + * Temporary SparseMatrix used to build the Paralution LocalMatrix. + */ + ::dealii::SparseMatrix sparse_matrix; + }; + + + +// ------------------- inline functions -------------- + + template + inline SparseMatrix::SparseMatrix() {} + + + + template + inline SparseMatrix::SparseMatrix(SparsityPattern const &sparsity_pattern) + { + reinit(sparsity_pattern); + } + + + + template + inline SparseMatrix::~SparseMatrix() + { + local_matrix.clear(); + } + + + + template + inline void SparseMatrix::reinit(SparsityPattern const &sparsity_pattern) + { + local_matrix.clear(); + sparse_matrix.reinit(sparsity_pattern); + } + + + + template + inline typename SparseMatrix::size_type SparseMatrix::m() const + { + return local_matrix.get_nrows(); + } + + + + template + inline typename SparseMatrix::size_type SparseMatrix::n() const + { + return local_matrix.get_ncols(); + } + + + + template + inline void SparseMatrix::add(const size_type i, + const size_type j, + const Number value) + { + sparse_matrix.add(i,j,value); + } + + + + template + template + inline void SparseMatrix::add(const std::vector &indices, + const FullMatrix &full_matrix, + const bool elide_zero_values) + { + sparse_matrix.add(indices,full_matrix,elide_zero_values); + } + + + + template + template + inline void SparseMatrix::add(const size_type row, + const std::vector &col_indices, + const std::vector &values, + const bool elide_zero_values) + { + sparse_matrix.add(row,col_indices,values,elide_zero_values); + } + + template + template + inline void SparseMatrix::add(const size_type row, + const size_type n_cols, + const size_type *col_indices, + const Number2 *values, + const bool elide_zero_values, + const bool col_indices_are_sorted) + { + sparse_matrix.add(row,n_cols,col_indices,values,elide_zero_values,col_indices_are_sorted); + } + + + + template + inline paralution::LocalMatrix const &SparseMatrix::paralution_matrix() const + { + return local_matrix; + } + + + + template + inline paralution::LocalMatrix& SparseMatrix::paralution_matrix() + { + return local_matrix; + } +} + + +/*@}*/ + +DEAL_II_NAMESPACE_CLOSE + +#endif // DEAL_II_WITH_PARALUTION + +/*---------------------------- paralution_sparse_matrix.h ---------------------------*/ + +#endif +/*---------------------------- paralution_sparse_matrix.h ---------------------------*/ diff --git a/deal.II/include/deal.II/lac/paralution_vector.h b/deal.II/include/deal.II/lac/paralution_vector.h index 7813163971..d0537b166d 100644 --- a/deal.II/include/deal.II/lac/paralution_vector.h +++ b/deal.II/include/deal.II/lac/paralution_vector.h @@ -40,7 +40,7 @@ namespace ParalutionWrappers * of vector is designed for use in either serial implementation or as a * localized copy on each processor. The implementation of this class is * based on the Paralution vector class LocalVector. The only requirement - * for this class to work is that Trilinos is installed with the same + * for this class to work is that Paralution is installed with the same * compiler as is used for compilation of deal.II. * * The interface of this class is modeled after the existing Vector class in diff --git a/deal.II/include/deal.II/lac/sparse_matrix.h b/deal.II/include/deal.II/lac/sparse_matrix.h index 6f9519ae27..6b4373df6b 100644 --- a/deal.II/include/deal.II/lac/sparse_matrix.h +++ b/deal.II/include/deal.II/lac/sparse_matrix.h @@ -549,7 +549,7 @@ public: }; /** - * @name Constructors and initalization + * @name Constructors and initialization */ //@{ /** diff --git a/deal.II/include/deal.II/lac/trilinos_sparse_matrix.h b/deal.II/include/deal.II/lac/trilinos_sparse_matrix.h index 15c8731d6c..82bf9eaedc 100644 --- a/deal.II/include/deal.II/lac/trilinos_sparse_matrix.h +++ b/deal.II/include/deal.II/lac/trilinos_sparse_matrix.h @@ -496,7 +496,7 @@ namespace TrilinosWrappers * * The interface of this class is modeled after the existing * SparseMatrix class in deal.II. It has almost the same member - * functions, and is often exchangable. However, since Trilinos only + * functions, and is often exchangeable. However, since Trilinos only * supports a single scalar type (double), it is not templated, and only * works with doubles. * diff --git a/deal.II/source/lac/CMakeLists.txt b/deal.II/source/lac/CMakeLists.txt index ebdcde19de..f5dddee3a7 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_sparse_matrix.cc petsc_block_sparse_matrix.cc petsc_full_matrix.cc petsc_matrix_base.cc diff --git a/deal.II/source/lac/paralution_sparse_matrix.cc b/deal.II/source/lac/paralution_sparse_matrix.cc new file mode 100644 index 0000000000..46cadb95d2 --- /dev/null +++ b/deal.II/source/lac/paralution_sparse_matrix.cc @@ -0,0 +1,92 @@ +// --------------------------------------------------------------------- +// $Id: paralution_sparse_matrix.cc 31567 2013-11-06 18:01:36Z turcksin $ +// +// 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 + +DEAL_II_NAMESPACE_OPEN + +namespace ParalutionWrappers +{ + template + void SparseMatrix::convert_to_paralution_csr() + { + local_matrix.AllocateCSR("deal_ii_local_matrix",sparse_matrix.n_nonzero_elements(), + sparse_matrix.m(),sparse_matrix.n()); + + //TODO: replace deprecated functions + std::vector val; + val.reserve(sparse_matrix.n_nonzero_elements()); + typename ::dealii::SparseMatrix::iterator it = sparse_matrix.begin(); + typename ::dealii::SparseMatrix::iterator end_it = sparse_matrix.end(); + for (; it!=end_it; ++it) + val.push_back(it->value()); + + const size_type *colnums_ptr = sparse_matrix.get_sparsity_pattern().get_column_numbers(); + const std::size_t *rowstart_ptr = sparse_matrix.get_sparsity_pattern().get_rowstart_indices(); + std::vector colnums(colnums_ptr,colnums_ptr+sparse_matrix.n_nonzero_elements()); + std::vector rowstart(rowstart_ptr,rowstart_ptr+sparse_matrix.m()+1); + + + // do the copying around of entries so that the diagonal entry is in the + // right place. note that this is easy to detect: since all entries apart + // from the diagonal entry are sorted, we know that the diagonal entry is + // in the wrong place if and only if its column index is larger than the + // column index of the second entry in a row + // + // ignore rows with only one or no entry + for (size_type row=0; row colnums[cursor+1])) + { + std::swap (colnums[cursor], colnums[cursor+1]); + std::swap (val[cursor], val[cursor+1]); + ++cursor; + } + } + + local_matrix.CopyFromCSR(&rowstart[0],&colnums[0],&val[0]); + + // Free the memory used by sparse_matrix. + sparse_matrix.clear(); + } +} + +// Explicit instantiations +namespace ParalutionWrappers +{ + template void SparseMatrix::convert_to_paralution_csr(); + + template void SparseMatrix::convert_to_paralution_csr(); +} + +DEAL_II_NAMESPACE_CLOSE + +#endif // DEAL_II_WITH_PARALUTION -- 2.39.5