From 525c84a2bbbbd187a82c2c3b27727d8016e9876f Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Sun, 21 Apr 2013 16:48:46 +0000 Subject: [PATCH] Make SparseDirectUMFPACK also understand block vectors. git-svn-id: https://svn.dealii.org/trunk@29354 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/doc/news/changes.h | 9 +- deal.II/include/deal.II/lac/sparse_direct.h | 80 +++++++-- deal.II/source/lac/sparse_direct.cc | 55 +++++++ tests/umfpack/umfpack_11.cc | 169 ++++++++++++++++++++ tests/umfpack/umfpack_11/cmp/generic | 25 +++ 5 files changed, 326 insertions(+), 12 deletions(-) create mode 100644 tests/umfpack/umfpack_11.cc create mode 100644 tests/umfpack/umfpack_11/cmp/generic diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 7fcd24d2e4..992cfe9e07 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -111,7 +111,14 @@ this function.
    -
  1. New: class TimerOutput::Scope does automatic scope based enter/ +
  2. New: SparseDirectUMFPACK has long had the ability to work with +BlockSparseMatrix objects, but couldn't deal with BlockVector objects. +This is now fixed. +
    +(Wolfgang Bangerth, 2013/04/21) +
    3. + +
  3. New: Class TimerOutput::Scope does automatic scope based enter/ exit_section of a TimerOutput object.
    (Timo Heister, 2013/04/18) diff --git a/deal.II/include/deal.II/lac/sparse_direct.h b/deal.II/include/deal.II/lac/sparse_direct.h index 195ca99190..da58687d0b 100644 --- a/deal.II/include/deal.II/lac/sparse_direct.h +++ b/deal.II/include/deal.II/lac/sparse_direct.h @@ -1,7 +1,7 @@ //--------------------------------------------------------------------------- // $Id$ // -// Copyright (C) 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012 by the deal.II authors +// Copyright (C) 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -100,8 +100,15 @@ public: ~SparseDirectUMFPACK (); /** - * This function does nothing. It is only here to provide a consistent - * interface. + * @name Setting up a sparse factorization + */ + /** + * @{ + */ + + /** + * This function does nothing. It is only here to provide a interface + * consistent with other sparse direct solvers. */ void initialize (const SparsityPattern &sparsity_pattern); @@ -133,30 +140,65 @@ public: void initialize(const Matrix &matrix, const AdditionalData additional_data = AdditionalData()); + /** + * @} + */ + + /** + * @name Functions that represent the inverse of a matrix + */ + /** + * @{ + */ + /** * Preconditioner interface function. Usually, given the source vector, * this method returns an approximated solution of Ax = b. As this * class provides a wrapper to a direct solver, here it is actually the * exact solution (exact within the range of numerical accuracy of * course). + * + * In other words, this function actually multiplies with the exact + * inverse of the matrix, $A^{-1}$. */ - void vmult (Vector &, const Vector &) const; + void vmult (Vector &dst, + const Vector &src) const; /** - * Not implemented but necessary for compiling. + * Same as before, but for block vectors. */ - void Tvmult (Vector &, const Vector &) const; + void vmult (BlockVector &dst, + const BlockVector &src) const; + + /** + * Not implemented but necessary for compiling certain other classes. + */ + void Tvmult (Vector &dst, + const Vector &src) const; /** * Same as vmult(), but adding to the previous solution. Not implemented - * yet. + * yet but necessary for compiling certain other classes. */ - void vmult_add (Vector &, const Vector &) const; + void vmult_add (Vector &dst, + const Vector &src) const; /** - * Not implemented but necessary for compiling. + * Not implemented but necessary for compiling certain other classes. + */ + void Tvmult_add (Vector &dst, + const Vector &src) const; + + /** + * @} + */ + + /** + * @name Functions that solve linear systems + */ + /** + * @{ */ - void Tvmult_add (Vector &, const Vector &) const; /** * Solve for a certain right hand side vector. This function may be @@ -175,7 +217,12 @@ public: void solve (Vector &rhs_and_solution) const; /** - * Call the two functions factorize and solve in that order, i.e. perform + * Same as before, but for block vectors. + */ + void solve (BlockVector &rhs_and_solution) const; + + /** + * Call the two functions factorize() and solve() in that order, i.e. perform * the whole solution process for the given right hand side vector. * * The solution will be returned in place of the right hand side vector. @@ -184,6 +231,17 @@ public: void solve (const Matrix &matrix, Vector &rhs_and_solution); + /** + * Same as before, but for block vectors. + */ + template + void solve (const Matrix &matrix, + BlockVector &rhs_and_solution); + + /** + * @} + */ + /** * One of the UMFPack routines threw an error. The error code is included * in the output and can be looked up in the UMFPack user manual. The diff --git a/deal.II/source/lac/sparse_direct.cc b/deal.II/source/lac/sparse_direct.cc index 8099df44d3..4df42981df 100644 --- a/deal.II/source/lac/sparse_direct.cc +++ b/deal.II/source/lac/sparse_direct.cc @@ -314,6 +314,19 @@ SparseDirectUMFPACK::solve (Vector &rhs_and_solution) const } +void +SparseDirectUMFPACK::solve (BlockVector &rhs_and_solution) const +{ + // the UMFPACK functions want a contiguous array of elements, so + // there is no way around copying data around. thus, just copy the + // data into a regular vector and back + Vector tmp (rhs_and_solution.size()); + tmp = rhs_and_solution; + solve (tmp); + rhs_and_solution = tmp; +} + + template void @@ -325,6 +338,16 @@ SparseDirectUMFPACK::solve (const Matrix &matrix, } +template +void +SparseDirectUMFPACK::solve (const Matrix &matrix, + BlockVector &rhs_and_solution) +{ + factorize (matrix); + solve (rhs_and_solution); +} + + #else @@ -355,6 +378,14 @@ SparseDirectUMFPACK::solve (Vector &) const } + +void +SparseDirectUMFPACK::solve (BlockVector &) const +{ + AssertThrow(false, ExcMessage("To call this function you need UMFPACK, but configured deal.II without passing the necessary switch to 'cmake'. Please consult the installation instructions in doc/readme.html.")); +} + + template void SparseDirectUMFPACK::solve (const Matrix &, @@ -363,6 +394,16 @@ SparseDirectUMFPACK::solve (const Matrix &, AssertThrow(false, ExcMessage("To call this function you need UMFPACK, but configured deal.II without passing the necessary switch to 'cmake'. Please consult the installation instructions in doc/readme.html.")); } + + +template +void +SparseDirectUMFPACK::solve (const Matrix &, + BlockVector &) +{ + AssertThrow(false, ExcMessage("To call this function you need UMFPACK, but configured deal.II without passing the necessary switch to 'cmake'. Please consult the installation instructions in doc/readme.html.")); +} + #endif @@ -385,6 +426,17 @@ SparseDirectUMFPACK::vmult ( } + +void +SparseDirectUMFPACK::vmult ( + BlockVector &dst, + const BlockVector &src) const +{ + dst = src; + this->solve(dst); +} + + void SparseDirectUMFPACK::Tvmult ( Vector &, @@ -596,6 +648,9 @@ void SparseDirectMUMPS::vmult (Vector &dst, void SparseDirectUMFPACK::solve (const MATRIX &, \ Vector &); \ template \ + void SparseDirectUMFPACK::solve (const MATRIX &, \ + BlockVector &); \ + template \ void SparseDirectUMFPACK::initialize (const MATRIX &, \ const AdditionalData); diff --git a/tests/umfpack/umfpack_11.cc b/tests/umfpack/umfpack_11.cc new file mode 100644 index 0000000000..bbf93cb68a --- /dev/null +++ b/tests/umfpack/umfpack_11.cc @@ -0,0 +1,169 @@ +//---------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2002, 2003, 2004, 2005, 2006, 2009, 2010, 2013 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. +// +//---------------------------------------------------------------------- + +// Like _11 but with block vectors. For reasons lost to history, +// SparseDirectUMFPACK could handle block sparse matrices but not block +// vectors in its solve() and vmult() functions. This was later +// corrected. This test checks this. + +#include "../tests.h" +#include +#include +#include + +#include +#include + +#include +#include + +#include + +#include +#include +#include +#include + +#include +#include + +#include +#include + + +template +void test () +{ + deallog << dim << 'd' << std::endl; + + Triangulation tria; + GridGenerator::hyper_cube (tria,0,1); + tria.refine_global (1); + + // destroy the uniformity of the matrix by + // refining one cell + tria.begin_active()->set_refine_flag (); + tria.execute_coarsening_and_refinement (); + tria.refine_global(8-2*dim); + + FE_Q fe (1); + DoFHandler dof_handler (tria); + dof_handler.distribute_dofs (fe); + + deallog << "Number of dofs = " << dof_handler.n_dofs() << std::endl; + + BlockSparsityPattern sparsity_pattern; + sparsity_pattern.reinit (3, 3); + for (unsigned int i=0; i<3; ++i) + for (unsigned int j=0; j<3; ++j) + sparsity_pattern.block(i,j).reinit (i!=2 ? + dof_handler.n_dofs()/3 : + dof_handler.n_dofs()-2*(dof_handler.n_dofs()/3), + j!=2 ? + dof_handler.n_dofs()/3 : + dof_handler.n_dofs()-2*(dof_handler.n_dofs()/3), + dof_handler.max_couplings_between_dofs(), + false); + sparsity_pattern.collect_sizes(); + DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern); + sparsity_pattern.compress (); + + BlockSparseMatrix B; + B.reinit (sparsity_pattern); + + { + // for some reason, we can't + // create block matrices directly + // in matrixtools... + SparsityPattern xsparsity_pattern; + xsparsity_pattern.reinit (dof_handler.n_dofs(), + dof_handler.n_dofs(), + dof_handler.max_couplings_between_dofs()); + DoFTools::make_sparsity_pattern (dof_handler, xsparsity_pattern); + xsparsity_pattern.compress (); + + SparseMatrix xB; + xB.reinit (xsparsity_pattern); + + QGauss qr (2); + MatrixTools::create_mass_matrix (dof_handler, qr, xB); + + // scale lower left part of the matrix by + // 1/2 and upper right part by 2 to make + // matrix nonsymmetric + for (SparseMatrix::iterator p=xB.begin(); + p!=xB.end(); ++p) + if (p->column() < p->row()) + p->value() = p->value()/2; + else if (p->column() > p->row()) + p->value() = p->value() * 2; + + // check that we've done it right + for (SparseMatrix::iterator p=xB.begin(); + p!=xB.end(); ++p) + if (p->column() != p->row()) + Assert (xB(p->row(),p->column()) != xB(p->column(),p->row()), + ExcInternalError()); + + // now copy stuff over + for (SparseMatrix::const_iterator i = xB.begin(); i != xB.end(); ++i) + B.set (i->row(), i->column(), i->value()); + } + + + // for a number of different solution + // vectors, make up a matching rhs vector + // and check what the UMFPACK solver finds + for (unsigned int i=0; i<3; ++i) + { + BlockVector solution (3); + solution.block(0).reinit(dof_handler.n_dofs()/3); + solution.block(1).reinit(dof_handler.n_dofs()/3); + solution.block(2).reinit(dof_handler.n_dofs()-2*(dof_handler.n_dofs()/3)); + solution.collect_sizes(); + BlockVector x; + x.reinit (solution); + BlockVector b; + b.reinit (solution); + + for (unsigned int j=0; j (); + test<2> (); + test<3> (); +} diff --git a/tests/umfpack/umfpack_11/cmp/generic b/tests/umfpack/umfpack_11/cmp/generic new file mode 100644 index 0000000000..8df3f884c1 --- /dev/null +++ b/tests/umfpack/umfpack_11/cmp/generic @@ -0,0 +1,25 @@ + +DEAL::1d +DEAL::Number of dofs = 193 +DEAL::relative norm distance = 0 +DEAL::absolute norms = 0 3084.00 +DEAL::relative norm distance = 0 +DEAL::absolute norms = 0 7709.99 +DEAL::relative norm distance = 0 +DEAL::absolute norms = 0 15420.0 +DEAL::2d +DEAL::Number of dofs = 1889 +DEAL::relative norm distance = 0 +DEAL::absolute norms = 0 94764.3 +DEAL::relative norm distance = 0 +DEAL::absolute norms = 0 236911. +DEAL::relative norm distance = 0 +DEAL::absolute norms = 0 473822. +DEAL::3d +DEAL::Number of dofs = 1333 +DEAL::relative norm distance = 0 +DEAL::absolute norms = 0 56165.6 +DEAL::relative norm distance = 0 +DEAL::absolute norms = 0 140414. +DEAL::relative norm distance = 0 +DEAL::absolute norms = 0 280828. -- 2.39.5