From c45aedfb1272761698e41c4cfe18979f60fb8d8e Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Thu, 27 Oct 2016 15:20:33 +0200 Subject: [PATCH] Add class for diagonal matrix --- include/deal.II/lac/diagonal_matrix.h | 295 ++++++++++++++++++++++++++ 1 file changed, 295 insertions(+) create mode 100644 include/deal.II/lac/diagonal_matrix.h diff --git a/include/deal.II/lac/diagonal_matrix.h b/include/deal.II/lac/diagonal_matrix.h new file mode 100644 index 0000000000..a9b58f9e9f --- /dev/null +++ b/include/deal.II/lac/diagonal_matrix.h @@ -0,0 +1,295 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 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 dealii__diagonal_matrix_h +#define dealii__diagonal_matrix_h + + +#include + +DEAL_II_NAMESPACE_OPEN + +/** + * This class represents a n x n diagonal matrix based on a vector of + * size n. The matrix-vector products are realized by @p + * VectorType::scale, so the template vector class needs to provide a + * @p scale() method. + * + * @author Martin Kronbichler, 2016 + */ +template > +class DiagonalMatrix : public Subscriptor +{ +public: + typedef typename VectorType::value_type value_type; + typedef typename VectorType::size_type size_type; + + /** + * Constructor. + */ + DiagonalMatrix(); + + /** + * Initialize with a given vector by copying the content of the vector + * @p vec. + */ + void reinit (const VectorType &vec); + + /** + * Returns a reference to the underlying vector for manipulation of the + * entries on the matrix diagonal. + */ + VectorType &get_vector(); + + /** + * Returns a read-only reference to the underlying vector. + */ + const VectorType &get_vector() const; + + /** + * Number of rows of this matrix. This number corresponds to the size of the + * underlying vector. + */ + size_type m () const; + + /** + * Number of columns of this matrix. To remember: this matrix is an n x + * n-matrix. + */ + size_type n () const; + + /** + * Read-only access to a value. This is restricted to the case where + * i==j due to the matrix storage. + * + * If the vector representing the diagonal is distributed with MPI, not all + * of the indices i might actually be accessible. Refer to the method + * get_vector().locally_owned_elements() for the entries that + * actually are accessible. + */ + value_type operator()(const size_type i, + const size_type j) const; + + /** + * Read-write access to a value. This is restricted to the case where + * i==j due to the matrix storage. + * + * If the vector representing the diagonal is distributed with MPI, not all + * of the indices i might actually be accessible. Refer to the method + * get_vector().locally_owned_elements() for the entries that + * actually are accessible. + */ + value_type &operator()(const size_type i, + const size_type j); + + /** + * Add an array of values given by values in the given global + * matrix row at columns specified by col_indices. Due to the storage of + * this matrix, entries are only added to the diagonal of the matrix. All + * other entries are ignored and no exception is thrown. + * + * This function is for a consistent interface with the other matrix classes + * in deal.II and can be used in + * ConstraintMatrix::distribute_local_to_global to get exactly the same + * diagonal as when assembling into a sparse matrix. + */ + 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); + + /** + * Performs a matrix-vector multiplication with the given matrix. + */ + void vmult (VectorType &dst, + const VectorType &src) const; + + /** + * Performs a transpose matrix-vector multiplication with the given + * matrix. Since this represents a diagonal matrix, exactly the same as + * vmult(). + */ + void Tvmult (VectorType &dst, + const VectorType &src) const; + + /** + * Adds the result of a matrix-vector multiplication into the destination + * vector dst. Needs to create a temporary vector, which makes performance + * slower than for @p vmult(). + */ + void vmult_add (VectorType &dst, + const VectorType &src) const; + + /** + * Adds the result of a transpose matrix-vector multiplication into the + * destination vector dst. Needs to create a temporary vector, which makes + * performance slower than for @p Tvmult(). + */ + void Tvmult_add (VectorType &dst, + const VectorType &src) const; + +private: + /** + * The stored vector. + */ + VectorType diagonal; +}; + +/* ---------------------------------- Inline functions ------------------- */ + +#ifndef DOXYGEN + +template +DiagonalMatrix::DiagonalMatrix() + : + Subscriptor() +{} + + + +template +void +DiagonalMatrix::reinit(const VectorType &vec) +{ + diagonal = vec; +} + + + +template +VectorType & +DiagonalMatrix::get_vector() +{ + return diagonal; +} + + + +template +const VectorType & +DiagonalMatrix::get_vector() const +{ + return diagonal; +} + + + +template +typename VectorType::size_type +DiagonalMatrix::m() const +{ + return diagonal.size(); +} + + + +template +typename VectorType::size_type +DiagonalMatrix::n() const +{ + return diagonal.size(); +} + + + +template +typename VectorType::value_type +DiagonalMatrix::operator()(const size_type i, + const size_type j) const +{ + Assert (i==j, ExcIndexRange(j,i,i+1)); + return diagonal(i); +} + + + +template +typename VectorType::value_type & +DiagonalMatrix::operator()(const size_type i, + const size_type j) +{ + Assert (i==j, ExcIndexRange(j,i,i+1)); + return diagonal(i); +} + + + +template +template +void +DiagonalMatrix::add (const size_type row, + const size_type n_cols, + const size_type *col_indices, + const number2 *values, + const bool , + const bool ) +{ + for (size_type i=0; i +void +DiagonalMatrix::vmult(VectorType &dst, + const VectorType &src) const +{ + dst = src; + dst.scale(diagonal); +} + + + +template +void +DiagonalMatrix::Tvmult(VectorType &dst, + const VectorType &src) const +{ + vmult(dst, src); +} + + + +template +void +DiagonalMatrix::vmult_add(VectorType &dst, + const VectorType &src) const +{ + VectorType tmp(src); + tmp.scale(diagonal); + dst += tmp; +} + + + +template +void +DiagonalMatrix::Tvmult_add(VectorType &dst, + const VectorType &src) const +{ + vmult_add(dst, src); +} + + +#endif + +DEAL_II_NAMESPACE_CLOSE + +#endif -- 2.39.5