From e1a5e5a8823d3f8f07bbe7cf2d1cec506119ad37 Mon Sep 17 00:00:00 2001 From: Guido Kanschat Date: Sat, 12 Mar 2005 02:51:48 +0000 Subject: [PATCH] new ProductSparseMatrix git-svn-id: https://svn.dealii.org/trunk@10106 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/doc/news/changes.html | 8 +++ deal.II/lac/include/lac/matrix_lib.h | 87 +++++++++++++++++++++++++++ deal.II/lac/source/matrix_lib.cc | 70 +++++++++++++++++++++- tests/lac/Makefile | 21 +++---- tests/lac/matrix_lib.cc | 89 ++++++++++++++++++++++++++++ 5 files changed, 261 insertions(+), 14 deletions(-) create mode 100644 tests/lac/matrix_lib.cc diff --git a/deal.II/doc/news/changes.html b/deal.II/doc/news/changes.html index b351e2cc7e..b0607df23e 100644 --- a/deal.II/doc/news/changes.html +++ b/deal.II/doc/news/changes.html @@ -161,6 +161,14 @@ inconvenience this causes.

lac

    +
  1. + New: The ProductSparseMatrix + implements the product of two rectangular sparse matrices with + the same value_type +
    + (GK, 2005/03/11) +

    +
  2. New: The PreconditionRichardson implements a Richardson preconditioner. diff --git a/deal.II/lac/include/lac/matrix_lib.h b/deal.II/lac/include/lac/matrix_lib.h index f073c1a7b4..e8bfd74245 100644 --- a/deal.II/lac/include/lac/matrix_lib.h +++ b/deal.II/lac/include/lac/matrix_lib.h @@ -9,6 +9,7 @@ template class Vector; template class BlockVector; +template class SparseMatrix; /*! @addtogroup Matrix2 *@{ @@ -96,6 +97,92 @@ class ProductMatrix : public PointerMatrixBase }; +/** + * Poor man's matrix product of two sparse matrices. Stores two + * matrices #m1 and #m2 of arbitrary type SparseMatrix and implements + * matrix-vector multiplications for the product + * M1M2 by performing multiplication with + * both factors consecutively. + * + * The documentation of @ProductMatrix applies with exception that + * these matrices here may be rectangular. + * + * @author Guido Kanschat, 2000, 2001, 2002, 2005 + */ +template +class ProductSparseMatrix : public PointerMatrixBase > +{ + public: + /** + * Define the type of matrices used. + */ + typedef SparseMatrix MatrixType; + + /** + * Define the type of vectors we + * plly this matrix to. + */ + typedef Vector VectorType; + + /** + * Constructor. Additionally to + * the two constituting matrices, a + * memory pool for the auxiliary + * vector must be provided. + */ + ProductSparseMatrix(const MatrixType& m1, + const MatrixType& m2, + VectorMemory& mem); + + /** + * Matrix-vector product w = + * m1 * m2 * v. + */ + virtual void vmult (VectorType& w, + const VectorType& v) const; + + /** + * Tranposed matrix-vector + * product w = m2T + * * m1T * v. + */ + virtual void Tvmult (VectorType& w, + const VectorType& v) const; + + /** + * Adding matrix-vector product + * w += m1 * m2 * v + */ + virtual void vmult_add (VectorType& w, + const VectorType& v) const; + + /** + * Adding, tranposed + * matrix-vector product w += + * m2T * + * m1T * v. + */ + virtual void Tvmult_add (VectorType& w, + const VectorType& v) const; + + private: + /** + * The left matrix of the product. + */ + SmartPointer m1; + + /** + * The right matrix of the product. + */ + SmartPointer m2; + + /** + * Memory for auxiliary vector. + */ + SmartPointer > mem; +}; + + /** * Mean value filter. The @p vmult functions of this matrix filter * out mean values of the vector. If the vector is of type diff --git a/deal.II/lac/source/matrix_lib.cc b/deal.II/lac/source/matrix_lib.cc index 3cf1345ed2..5e16eb1e50 100644 --- a/deal.II/lac/source/matrix_lib.cc +++ b/deal.II/lac/source/matrix_lib.cc @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003 by the deal.II authors +// Copyright (C) 1998 - 2005 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -13,12 +13,80 @@ #include +#include MeanValueFilter::MeanValueFilter(unsigned int component) : component(component) {} + +template +ProductSparseMatrix::ProductSparseMatrix( + const MatrixType& mat1, + const MatrixType& mat2, + VectorMemory& mem) + : + m1(&mat1), m2(&mat2), + mem(&mem) +{ + Assert(mat1.n() == mat2.m(), ExcDimensionMismatch(mat1.n(),mat2.m())); +} + + +template +void +ProductSparseMatrix::vmult (VectorType& dst, const VectorType& src) const +{ + VectorType* v = mem->alloc(); + v->reinit(m1->n()); + m2->vmult (*v, src); + m1->vmult (dst, *v); + mem->free(v); +} + + +template +void +ProductSparseMatrix::vmult_add (VectorType& dst, const VectorType& src) const +{ + VectorType* v = mem->alloc(); + v->reinit(m1->n()); + m2->vmult (*v, src); + m1->vmult_add (dst, *v); + mem->free(v); +} + + +template +void +ProductSparseMatrix::Tvmult (VectorType& dst, const VectorType& src) const +{ + VectorType* v = mem->alloc(); + v->reinit(m1->n()); + m1->Tvmult (*v, src); + m2->Tvmult (dst, *v); + mem->free(v); +} + + +template +void +ProductSparseMatrix::Tvmult_add (VectorType& dst, const VectorType& src) const +{ + VectorType* v = mem->alloc(); + v->reinit(m1->n()); + m1->Tvmult (*v, src); + m2->Tvmult_add (dst, *v); + mem->free(v); +} + + +template class ProductSparseMatrix; +template class ProductSparseMatrix; +template class ProductSparseMatrix; +template class ProductSparseMatrix; + template void MeanValueFilter::filter(Vector&) const; template void MeanValueFilter::filter(Vector&) const; template void MeanValueFilter::filter(BlockVector&) const; diff --git a/tests/lac/Makefile b/tests/lac/Makefile index 4a885957a2..533f165581 100644 --- a/tests/lac/Makefile +++ b/tests/lac/Makefile @@ -1,6 +1,6 @@ ############################################################ # $Id$ -# Copyright (C) 2000, 2001, 2002, 2003, 2004 by the deal.II authors +# Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005 by the deal.II authors ############################################################ ############################################################ @@ -20,18 +20,13 @@ default: run-tests ############################################################ -tests_x = $(sort \ - sparse_matrices \ - block_matrices \ - block_vector \ - block_vector_iterator \ - full_matrix \ - matrix_out \ - solver \ - eigen \ - sparsity_pattern \ - sparse_ilu \ - vector-vector) +tests_x = vector-vector \ + full_matrix sparsity_pattern sparse_matrices \ + block_vector block_vector_iterator block_matrices \ + matrix_lib matrix_out \ + solver eigen \ + sparse_ilu + # from above list of regular expressions, generate the real set of # tests diff --git a/tests/lac/matrix_lib.cc b/tests/lac/matrix_lib.cc new file mode 100644 index 0000000000..5f9e58476d --- /dev/null +++ b/tests/lac/matrix_lib.cc @@ -0,0 +1,89 @@ +//--------------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2005 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//--------------------------------------------------------------------------- + + +#include "../tests.h" +#include +#include +#include +#include +#include + +#include +#include +#include + +template +void check_sparse_product(const SparseMatrix& m1, SparseMatrix& m2) +{ + Vector v(m2.n()); + Vector w(m1.m()); + deallog << "Sizes\t" << v.size() << '\t' << w.size() << std::endl; + + for (unsigned int i=0;i > mem; + + ProductSparseMatrix product(m1, m2, mem); + product.vmult(w,v); + + for (unsigned int i=0;i m1(sparsity1); + SparseMatrix m2(sparsity2); + + for (unsigned int i=0;i<2;++i) + for (unsigned int j=0;j<3;++j) + { + m1.set(i, j, 1.*i-j); + m2.set(j, i, 1.*i-j); + m2.set(j, 2+i, 1.*j-i); + } + check_sparse_product(m1, m2); +} -- 2.39.5