From 36bacdbea22c5888b004e7c503904404499eada3 Mon Sep 17 00:00:00 2001 From: bangerth Date: Sun, 30 Oct 2011 14:32:26 +0000 Subject: [PATCH] Fix SparseMatrix::{T,}mmult in two different regards. git-svn-id: https://svn.dealii.org/trunk@24702 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/doc/news/changes.h | 16 ++- deal.II/include/deal.II/lac/sparse_matrix.h | 28 ++++-- .../deal.II/lac/sparse_matrix.templates.h | 34 +++++-- deal.II/source/lac/sparse_matrix.inst.in | 28 +++--- tests/lac/sparse_matrix_Tmmult_01.cc | 92 ++++++++++++++++++ tests/lac/sparse_matrix_Tmmult_01/cmp/generic | 3 + tests/lac/sparse_matrix_Tmmult_02.cc | 97 +++++++++++++++++++ tests/lac/sparse_matrix_Tmmult_02/cmp/generic | 3 + tests/lac/sparse_matrix_mmult_01.cc | 92 ++++++++++++++++++ tests/lac/sparse_matrix_mmult_01/cmp/generic | 3 + tests/lac/sparse_matrix_mmult_02.cc | 97 +++++++++++++++++++ tests/lac/sparse_matrix_mmult_02/cmp/generic | 3 + 12 files changed, 466 insertions(+), 30 deletions(-) create mode 100644 tests/lac/sparse_matrix_Tmmult_01.cc create mode 100644 tests/lac/sparse_matrix_Tmmult_01/cmp/generic create mode 100644 tests/lac/sparse_matrix_Tmmult_02.cc create mode 100644 tests/lac/sparse_matrix_Tmmult_02/cmp/generic create mode 100644 tests/lac/sparse_matrix_mmult_01.cc create mode 100644 tests/lac/sparse_matrix_mmult_01/cmp/generic create mode 100644 tests/lac/sparse_matrix_mmult_02.cc create mode 100644 tests/lac/sparse_matrix_mmult_02/cmp/generic diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index ff36636034..d66ce3821b 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -48,10 +48,20 @@ enabled due to a missing include file in file

Specific improvements

    +
  1. Fixed: SparseMatrix::mmult and SpareMatrix::Tmmult had a number of +issues that are now fixed: (i) rebuilding the sparsity pattern was allowed +even if several of the matrices involved in these operations shared a +sparsity pattern; (ii) the functions had a vector argument that had a default +value but the default value could not be used because it wasn't used in a +template context deducible by the compiler. +
    +(Wolfgang Bangerth, 2011/10/30) -
  2. New: parallel::distributed::Triangulation::mesh_reconstruction_after_repartitioning -setting which is necessary for save()/load() to be deterministic. Otherwise the matrix assembly -is done in a different order depending on the order of old refinements. +
  3. New: +parallel::distributed::Triangulation::mesh_reconstruction_after_repartitioning +setting which is necessary for save()/load() to be deterministic. Otherwise +the matrix assembly is done in a different order depending on the order of old +refinements.
    (Timo Heister, 2011/10/26) diff --git a/deal.II/include/deal.II/lac/sparse_matrix.h b/deal.II/include/deal.II/lac/sparse_matrix.h index 13ef275822..46274071b9 100644 --- a/deal.II/include/deal.II/lac/sparse_matrix.h +++ b/deal.II/include/deal.II/lac/sparse_matrix.h @@ -1,7 +1,7 @@ //--------------------------------------------------------------------------- // $Id$ // -// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2009, 2010 by the deal.II authors +// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2009, 2010, 2011 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -1504,12 +1504,19 @@ class SparseMatrix : public virtual Subscriptor * and instead uses the sparsity * pattern stored in C. In * that case, make sure that it really - * fits. + * fits. The default is to rebuild the + * sparsity pattern. + * + * @note Rebuilding the sparsity pattern + * requires changing it. This means that + * all other matrices that are associated + * with this sparsity pattern will + * then have invalid entries. */ - template + template void mmult (SparseMatrix &C, const SparseMatrix &B, - const Vector &V = Vector(), + const Vector &V = Vector(), const bool rebuild_sparsity_pattern = true) const; /** @@ -1545,12 +1552,19 @@ class SparseMatrix : public virtual Subscriptor * and instead uses the sparsity * pattern stored in C. In * that case, make sure that it really - * fits. + * fits. The default is to rebuild the + * sparsity pattern. + * + * @note Rebuilding the sparsity pattern + * requires changing it. This means that + * all other matrices that are associated + * with this sparsity pattern will + * then have invalid entries. */ - template + template void Tmmult (SparseMatrix &C, const SparseMatrix &B, - const Vector &V = Vector(), + const Vector &V = Vector(), const bool rebuild_sparsity_pattern = true) const; //@} diff --git a/deal.II/include/deal.II/lac/sparse_matrix.templates.h b/deal.II/include/deal.II/lac/sparse_matrix.templates.h index 1c899618cc..0550a859ff 100644 --- a/deal.II/include/deal.II/lac/sparse_matrix.templates.h +++ b/deal.II/include/deal.II/lac/sparse_matrix.templates.h @@ -897,11 +897,11 @@ SparseMatrix::matrix_scalar_product (const Vector& u, template -template +template void SparseMatrix::mmult (SparseMatrix &C, const SparseMatrix &B, - const Vector &V, + const Vector &V, const bool rebuild_sparsity_C) const { const bool use_vector = V.size() == n() ? true : false; @@ -916,6 +916,17 @@ SparseMatrix::mmult (SparseMatrix &C, // clear previous content of C if (rebuild_sparsity_C == true) { + // we are about to change the sparsity + // pattern of C. this can not work if + // either A or B use the same sparsity + // pattern + Assert (&C.get_sparsity_pattern() != &this->get_sparsity_pattern(), + ExcMessage ("Can't use the same sparsity pattern for " + "different matrices if it is to be rebuilt.")); + Assert (&C.get_sparsity_pattern() != &B.get_sparsity_pattern(), + ExcMessage ("Can't use the same sparsity pattern for " + "different matrices if it is to be rebuilt.")); + // need to change the sparsity pattern of // C, so cast away const-ness. SparsityPattern & sp_C = @@ -1034,11 +1045,11 @@ SparseMatrix::mmult (SparseMatrix &C, template -template +template void SparseMatrix::Tmmult (SparseMatrix &C, const SparseMatrix &B, - const Vector &V, + const Vector &V, const bool rebuild_sparsity_C) const { const bool use_vector = V.size() == m() ? true : false; @@ -1053,8 +1064,19 @@ SparseMatrix::Tmmult (SparseMatrix &C, // clear previous content of C if (rebuild_sparsity_C == true) { - // need to change the sparsity pattern of - // C, so cast away const-ness. + // we are about to change the sparsity + // pattern of C. this can not work if + // either A or B use the same sparsity + // pattern + Assert (&C.get_sparsity_pattern() != &this->get_sparsity_pattern(), + ExcMessage ("Can't use the same sparsity pattern for " + "different matrices if it is to be rebuilt.")); + Assert (&C.get_sparsity_pattern() != &B.get_sparsity_pattern(), + ExcMessage ("Can't use the same sparsity pattern for " + "different matrices if it is to be rebuilt.")); + + // need to change the sparsity pattern of + // C, so cast away const-ness. SparsityPattern & sp_C = *(const_cast(&C.get_sparsity_pattern())); C.clear(); diff --git a/deal.II/source/lac/sparse_matrix.inst.in b/deal.II/source/lac/sparse_matrix.inst.in index f5b3b0eabd..2ce3a691ab 100644 --- a/deal.II/source/lac/sparse_matrix.inst.in +++ b/deal.II/source/lac/sparse_matrix.inst.in @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2006, 2007, 2009 by the deal.II authors +// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2006, 2007, 2009, 2011 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -26,22 +26,22 @@ for (S1, S2 : REAL_SCALARS) template SparseMatrix & SparseMatrix::copy_from (const SparseMatrix &); - template + template void SparseMatrix::copy_from (const FullMatrix &); template void SparseMatrix::add (const S1, const SparseMatrix &); template void SparseMatrix::add (const unsigned int, - const unsigned int, - const unsigned int *, + const unsigned int, + const unsigned int *, const S2 *, - const bool, + const bool, const bool); template void SparseMatrix::set (const unsigned int, - const unsigned int, - const unsigned int *, + const unsigned int, + const unsigned int *, const S2 *, const bool); } @@ -117,7 +117,7 @@ for (S1, S2 : REAL_SCALARS) const S1) const; template void SparseMatrix:: SSOR_step (Vector &, - const Vector &, + const Vector &, const S1) const; } @@ -134,14 +134,14 @@ for (S1, S2, S3 : REAL_SCALARS; Tvmult_add (V1 &, const V2 &) const; } -for (S1, S2, S3, S4: REAL_SCALARS) +for (S1, S2, S3: REAL_SCALARS) { template void SparseMatrix:: - mmult (SparseMatrix &, const SparseMatrix &, const Vector&, + mmult (SparseMatrix &, const SparseMatrix &, const Vector&, const bool) const; template void SparseMatrix:: - Tmmult (SparseMatrix &, const SparseMatrix &, const Vector&, - const bool) const; + Tmmult (SparseMatrix &, const SparseMatrix &, const Vector&, + const bool) const; } @@ -160,7 +160,7 @@ for (S1, S2, S3, S4: REAL_SCALARS) // template SparseMatrix & // SparseMatrix::copy_from (const SparseMatrix &); -// template +// template // void SparseMatrix::copy_from (const FullMatrix &); // template void SparseMatrix::add (const S1, @@ -233,7 +233,7 @@ for (S1, S2, S3, S4: REAL_SCALARS) // const S1) const; // template void SparseMatrix:: // SSOR_step (Vector &, -// const Vector &, +// const Vector &, // const S1) const; // } diff --git a/tests/lac/sparse_matrix_Tmmult_01.cc b/tests/lac/sparse_matrix_Tmmult_01.cc new file mode 100644 index 0000000000..5cfb45991a --- /dev/null +++ b/tests/lac/sparse_matrix_Tmmult_01.cc @@ -0,0 +1,92 @@ +//---------------------------- sparse_matrix_Tmmult_01.cc,v --------------------------- +// $Id$ +// +// Copyright (C) 2011 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. +// +//---------------------------- sparse_matrix_Tmmult_01.cc,v --------------------------- + +// check SparseMatrix::Tmmult. this function has a default argument +// that could previously not be instantiated because it was in a +// non-deduced context but that should not be possible to omit +// +// this test checks SparseMatrix::Tmmult without additional arguments + +#include "../tests.h" +#include +#include +#include +#include +#include + +#include +#include +#include +#include + + +void test (const unsigned int n) +{ + // Create some random full matrices in the + // data structures of a sparse matrix + SparsityPattern sp (n,n); + SparsityPattern C_sp (n,n); + for (unsigned int i=0;i A(sp); + SparseMatrix B(sp); + SparseMatrix C(C_sp); + + for (unsigned int i=0;i x(n), y(n), z(n), tmp(n); + for (unsigned int j=0;j +#include +#include +#include +#include + +#include +#include +#include +#include + + +void test (const unsigned int n) +{ + // Create some random full matrices in the + // data structures of a sparse matrix + SparsityPattern sp (n,n); + SparsityPattern C_sp (n,n); + for (unsigned int i=0;i A(sp); + SparseMatrix B(sp); + SparseMatrix C(C_sp); + + for (unsigned int i=0;i v(n); + for (unsigned int j=0;j x(n), y(n), z(n), tmp(n); + for (unsigned int j=0;j +#include +#include +#include +#include + +#include +#include +#include +#include + + +void test (const unsigned int n) +{ + // Create some random full matrices in the + // data structures of a sparse matrix + SparsityPattern sp (n,n); + SparsityPattern C_sp (n,n); + for (unsigned int i=0;i A(sp); + SparseMatrix B(sp); + SparseMatrix C(C_sp); + + for (unsigned int i=0;i x(n), y(n), z(n), tmp(n); + for (unsigned int j=0;j +#include +#include +#include +#include + +#include +#include +#include +#include + + +void test (const unsigned int n) +{ + // Create some random full matrices in the + // data structures of a sparse matrix + SparsityPattern sp (n,n); + SparsityPattern C_sp (n,n); + for (unsigned int i=0;i A(sp); + SparseMatrix B(sp); + SparseMatrix C(C_sp); + + for (unsigned int i=0;i v(n); + for (unsigned int j=0;j x(n), y(n), z(n), tmp(n); + for (unsigned int j=0;j