From 438b10ec7ca7ed79f31840900656aad45fb238d6 Mon Sep 17 00:00:00 2001 From: wolf Date: Tue, 25 Apr 2006 23:12:39 +0000 Subject: [PATCH] SparseDirectUmfpack now also works for sparse matrices with floats, as well as for block sparse matrices. git-svn-id: https://svn.dealii.org/trunk@12897 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/doc/news/changes.html | 9 +- deal.II/lac/include/lac/sparse_direct.h | 21 ++- deal.II/lac/source/sparse_direct.cc | 93 +++++++--- tests/bits/umfpack_04.cc | 133 ++++++++++++++ tests/bits/umfpack_04/.cvsignore | 1 + tests/bits/umfpack_04/cmp/generic | 25 +++ .../cmp/x86_64-unknown-linux-gnu+gcc3.3 | 25 +++ tests/bits/umfpack_05.cc | 171 ++++++++++++++++++ tests/bits/umfpack_05/.cvsignore | 1 + tests/bits/umfpack_05/cmp/generic | 25 +++ .../cmp/x86_64-unknown-linux-gnu+gcc3.3 | 25 +++ tests/bits/umfpack_06.cc | 171 ++++++++++++++++++ tests/bits/umfpack_06/.cvsignore | 1 + tests/bits/umfpack_06/cmp/generic | 25 +++ .../cmp/x86_64-unknown-linux-gnu+gcc3.3 | 25 +++ 15 files changed, 716 insertions(+), 35 deletions(-) create mode 100644 tests/bits/umfpack_04.cc create mode 100644 tests/bits/umfpack_04/.cvsignore create mode 100644 tests/bits/umfpack_04/cmp/generic create mode 100644 tests/bits/umfpack_04/cmp/x86_64-unknown-linux-gnu+gcc3.3 create mode 100644 tests/bits/umfpack_05.cc create mode 100644 tests/bits/umfpack_05/.cvsignore create mode 100644 tests/bits/umfpack_05/cmp/generic create mode 100644 tests/bits/umfpack_05/cmp/x86_64-unknown-linux-gnu+gcc3.3 create mode 100644 tests/bits/umfpack_06.cc create mode 100644 tests/bits/umfpack_06/.cvsignore create mode 100644 tests/bits/umfpack_06/cmp/generic create mode 100644 tests/bits/umfpack_06/cmp/x86_64-unknown-linux-gnu+gcc3.3 diff --git a/deal.II/doc/news/changes.html b/deal.II/doc/news/changes.html index 3ae3abe42c..f4d0497657 100644 --- a/deal.II/doc/news/changes.html +++ b/deal.II/doc/news/changes.html @@ -442,6 +442,14 @@ inconvenience this causes.
    +
  1. Improved: The SparseDirectUMFPACK solver can now also be + used with sparse matrices with elements of type float, as well + as with block matrices with types float and double. +
    + (WB 2006/04/25) +

    +
  2. New: The function BlockSparsityPattern::row_length adds up the row lengths of the individual blocks of a block matrix for @@ -449,7 +457,6 @@ inconvenience this causes.
    (WB 2006/04/25)

    -
  3. New: There is a new class IdentityMatrix that represents an diff --git a/deal.II/lac/include/lac/sparse_direct.h b/deal.II/lac/include/lac/sparse_direct.h index 3bfa14328b..70b136bdd1 100644 --- a/deal.II/lac/include/lac/sparse_direct.h +++ b/deal.II/lac/include/lac/sparse_direct.h @@ -1027,7 +1027,13 @@ class SparseDirectMA47 : public Subscriptor * from the deal.II ReadMe file. UMFPACK is included courtesy of its author, * Timothy A. Davis. * - * + * + *

    Instantiations

    + * + * There are instantiations of this class for SparseMatrix, + * SparseMatrix, BlockSparseMatrix, and + * BlockSparseMatrix. + * * @ingroup Solvers Preconditioners * @see @ref SoftwareUMFPACK * @@ -1093,14 +1099,16 @@ class SparseDirectUMFPACK : public Subscriptor * this operation, even if subsequent * solves are required. */ - void factorize (const SparseMatrix &matrix); + template + void factorize (const Matrix &matrix); /** * Initialize memory and call * SparseDirectUMFPACK::factorize. */ - void initialize(const SparseMatrix &matrix, - const AdditionalData = AdditionalData()); + template + void initialize(const Matrix &matrix, + const AdditionalData additional_data = AdditionalData()); /** * Preconditioner interface @@ -1168,8 +1176,9 @@ class SparseDirectUMFPACK : public Subscriptor * in place of the right hand * side vector. */ - void solve (const SparseMatrix &matrix, - Vector &rhs_and_solution); + template + void solve (const Matrix &matrix, + Vector &rhs_and_solution); /** * One of the UMFPack routines diff --git a/deal.II/lac/source/sparse_direct.cc b/deal.II/lac/source/sparse_direct.cc index 9534cee1a3..0a3384219f 100644 --- a/deal.II/lac/source/sparse_direct.cc +++ b/deal.II/lac/source/sparse_direct.cc @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 2001, 2002, 2003, 2004, 2005 by the deal.II authors +// Copyright (C) 2001, 2002, 2003, 2004, 2005, 2006 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -16,6 +16,7 @@ #include #include #include +#include #include #include @@ -1679,9 +1680,10 @@ SparseDirectUMFPACK::clear () +template void SparseDirectUMFPACK:: -factorize (const SparseMatrix &matrix) +factorize (const Matrix &matrix) { Assert (matrix.m() == matrix.n(), ExcNotQuadratic()) @@ -1729,16 +1731,37 @@ factorize (const SparseMatrix &matrix) Assert (static_cast(Ap.back()) == Ai.size(), ExcInternalError()); - // then copy over matrix elements + // then copy over matrix + // elements. note that for sparse + // matrices, iterators are sorted + // so that they traverse each row + // from start to end before moving + // on to the next row. however, + // this isn't true for block + // matrices, so we have to do a bit + // of book keeping { - unsigned int index = 0; - for (SparseMatrix::const_iterator p=matrix.begin(); - p!=matrix.end(); ++p, ++index) + // have an array that for each + // row points to the first entry + // not yet written to + std::vector row_pointers = Ap; + + for (typename Matrix::const_iterator p=matrix.begin(); + p!=matrix.end(); ++p) { - Ai[index] = p->column(); - Ax[index] = p->value(); + // write entry into the first + // free one for this row + Ai[row_pointers[p->row()]] = p->column(); + Ax[row_pointers[p->row()]] = p->value(); + + // then move pointer ahead + ++row_pointers[p->row()]; } - Assert (index == Ai.size(), ExcInternalError()); + + // at the end, we should have + // written all rows completely + for (unsigned int i=0; i &matrix) } } } - - int status; @@ -1835,9 +1856,10 @@ SparseDirectUMFPACK::solve (Vector &rhs_and_solution) const +template void -SparseDirectUMFPACK::solve (const SparseMatrix &matrix, - Vector &rhs_and_solution) +SparseDirectUMFPACK::solve (const Matrix &matrix, + Vector &rhs_and_solution) { factorize (matrix); solve (rhs_and_solution); @@ -1860,7 +1882,8 @@ SparseDirectUMFPACK::clear () {} -void SparseDirectUMFPACK::factorize (const SparseMatrix &) +template +void SparseDirectUMFPACK::factorize (const Matrix &) { Assert(false, ExcNeedsUMFPACK()); } @@ -1873,9 +1896,10 @@ SparseDirectUMFPACK::solve (Vector &) const } +template void -SparseDirectUMFPACK::solve (const SparseMatrix &, - Vector &) +SparseDirectUMFPACK::solve (const Matrix &, + Vector &) { Assert(false, ExcNeedsUMFPACK()); } @@ -1884,8 +1908,9 @@ SparseDirectUMFPACK::solve (const SparseMatrix &, #endif +template void -SparseDirectUMFPACK::initialize (const SparseMatrix& M, +SparseDirectUMFPACK::initialize (const Matrix &M, const AdditionalData) { this->factorize(M); @@ -1929,22 +1954,34 @@ SparseDirectUMFPACK::Tvmult_add ( } -// explicit instantiations +// explicit instantiations for SparseMatrixMA27 template -void -SparseDirectMA27::factorize (const SparseMatrix &matrix); +void SparseDirectMA27::factorize (const SparseMatrix &matrix); template -void -SparseDirectMA27::factorize (const SparseMatrix &matrix); +void SparseDirectMA27::factorize (const SparseMatrix &matrix); template -void -SparseDirectMA27::solve (const SparseMatrix &matrix, +void SparseDirectMA27::solve (const SparseMatrix &matrix, Vector &rhs_and_solution); template -void -SparseDirectMA27::solve (const SparseMatrix &matrix, - Vector &rhs_and_solution); - +void SparseDirectMA27::solve (const SparseMatrix &matrix, + Vector &rhs_and_solution); + + +// explicit instantiations for SparseMatrixUMFPACK +#define InstantiateUMFPACK(MATRIX) \ + template \ + void SparseDirectUMFPACK::factorize (const MATRIX &); \ + template \ + void SparseDirectUMFPACK::solve (const MATRIX &, \ + Vector &); \ + template \ + void SparseDirectUMFPACK::initialize (const MATRIX &, \ + const AdditionalData) + +InstantiateUMFPACK(SparseMatrix); +InstantiateUMFPACK(SparseMatrix); +InstantiateUMFPACK(BlockSparseMatrix); +InstantiateUMFPACK(BlockSparseMatrix); diff --git a/tests/bits/umfpack_04.cc b/tests/bits/umfpack_04.cc new file mode 100644 index 0000000000..8d605f4247 --- /dev/null +++ b/tests/bits/umfpack_04.cc @@ -0,0 +1,133 @@ +//---------------------------- umfpack_04.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2002, 2003, 2004, 2005, 2006 by the deal.II authors and Brian Carnes +// +// 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. +// +//---------------------------- umfpack_04.cc --------------------------- + +// test the umfpack sparse direct solver on a mass matrix that is +// slightly modified to make it nonsymmetric. same as umfpack_03, but +// for a SparseMatrix instead of SparseMatrix matrix + +#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; + + SparsityPattern sparsity_pattern; + sparsity_pattern.reinit (dof_handler.n_dofs(), + dof_handler.n_dofs(), + dof_handler.max_couplings_between_dofs()); + DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern); + sparsity_pattern.compress (); + + SparseMatrix B; + B.reinit (sparsity_pattern); + + QGauss qr (2); + MatrixTools::create_mass_matrix (dof_handler, qr, B); + + // scale lower left part of the matrix by + // 1/2 and upper right part by 2 to make + // matrix nonsymmetric + for (SparseMatrix::iterator p=B.begin(); + p!=B.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=B.begin(); + p!=B.end(); ++p) + if (p->column() != p->row()) + Assert (B(p->row(),p->column()) != B(p->column(),p->row()), + ExcInternalError()); + + // 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) + { + Vector solution (dof_handler.n_dofs()); + Vector x (dof_handler.n_dofs()); + Vector b (dof_handler.n_dofs()); + + for (unsigned int j=0; j (); + test<2> (); + test<3> (); +} diff --git a/tests/bits/umfpack_04/.cvsignore b/tests/bits/umfpack_04/.cvsignore new file mode 100644 index 0000000000..d196dfb299 --- /dev/null +++ b/tests/bits/umfpack_04/.cvsignore @@ -0,0 +1 @@ +obj.* exe OK output diff --git a/tests/bits/umfpack_04/cmp/generic b/tests/bits/umfpack_04/cmp/generic new file mode 100644 index 0000000000..b9cb7d0d35 --- /dev/null +++ b/tests/bits/umfpack_04/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 = 1.17531e-10 236911. +DEAL::relative norm distance = 0 +DEAL::absolute norms = 2.35062e-10 473822. +DEAL::3d +DEAL::Number of dofs = 1333 +DEAL::relative norm distance = 0 +DEAL::absolute norms = 1.26673e-10 56165.6 +DEAL::relative norm distance = 0 +DEAL::absolute norms = 2.98388e-10 140414. +DEAL::relative norm distance = 0 +DEAL::absolute norms = 5.96775e-10 280828. diff --git a/tests/bits/umfpack_04/cmp/x86_64-unknown-linux-gnu+gcc3.3 b/tests/bits/umfpack_04/cmp/x86_64-unknown-linux-gnu+gcc3.3 new file mode 100644 index 0000000000..3176fb4ed1 --- /dev/null +++ b/tests/bits/umfpack_04/cmp/x86_64-unknown-linux-gnu+gcc3.3 @@ -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 = 1.23087e-10 236911. +DEAL::relative norm distance = 0 +DEAL::absolute norms = 2.46173e-10 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 = 4.46191e-10 140414. +DEAL::relative norm distance = 0 +DEAL::absolute norms = 8.92381e-10 280828. diff --git a/tests/bits/umfpack_05.cc b/tests/bits/umfpack_05.cc new file mode 100644 index 0000000000..09436b1ff7 --- /dev/null +++ b/tests/bits/umfpack_05.cc @@ -0,0 +1,171 @@ +//---------------------------- umfpack_05.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2002, 2003, 2004, 2005, 2006 by the deal.II authors and Brian Carnes +// +// 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. +// +//---------------------------- umfpack_05.cc --------------------------- + +// test the umfpack sparse direct solver on a mass matrix that is +// slightly modified to make it nonsymmetric. same as umfpack_03, but +// for a BlockSparseMatrix instead of SparseMatrix matrix + +#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 tmp (solution.size()); + tmp = x; + SparseDirectUMFPACK().solve (B, tmp); + x = tmp; + + x -= solution; + deallog << "relative norm distance = " + << x.l2_norm() / solution.l2_norm() + << std::endl; + deallog << "absolute norms = " + << x.l2_norm() << ' ' << solution.l2_norm() + << std::endl; + Assert (x.l2_norm() / solution.l2_norm() < 1e-8, + ExcInternalError()); + } +} + + +int main () +{ + std::ofstream logfile("umfpack_05/output"); + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + test<1> (); + test<2> (); + test<3> (); +} diff --git a/tests/bits/umfpack_05/.cvsignore b/tests/bits/umfpack_05/.cvsignore new file mode 100644 index 0000000000..d196dfb299 --- /dev/null +++ b/tests/bits/umfpack_05/.cvsignore @@ -0,0 +1 @@ +obj.* exe OK output diff --git a/tests/bits/umfpack_05/cmp/generic b/tests/bits/umfpack_05/cmp/generic new file mode 100644 index 0000000000..b9cb7d0d35 --- /dev/null +++ b/tests/bits/umfpack_05/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 = 1.17531e-10 236911. +DEAL::relative norm distance = 0 +DEAL::absolute norms = 2.35062e-10 473822. +DEAL::3d +DEAL::Number of dofs = 1333 +DEAL::relative norm distance = 0 +DEAL::absolute norms = 1.26673e-10 56165.6 +DEAL::relative norm distance = 0 +DEAL::absolute norms = 2.98388e-10 140414. +DEAL::relative norm distance = 0 +DEAL::absolute norms = 5.96775e-10 280828. diff --git a/tests/bits/umfpack_05/cmp/x86_64-unknown-linux-gnu+gcc3.3 b/tests/bits/umfpack_05/cmp/x86_64-unknown-linux-gnu+gcc3.3 new file mode 100644 index 0000000000..92c6c405f8 --- /dev/null +++ b/tests/bits/umfpack_05/cmp/x86_64-unknown-linux-gnu+gcc3.3 @@ -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 = 1.07697e-10 236911. +DEAL::relative norm distance = 0 +DEAL::absolute norms = 2.15395e-10 473822. +DEAL::3d +DEAL::Number of dofs = 1333 +DEAL::relative norm distance = 0 +DEAL::absolute norms = 2.35179e-10 56165.6 +DEAL::relative norm distance = 0 +DEAL::absolute norms = 9.26398e-10 140414. +DEAL::relative norm distance = 0 +DEAL::absolute norms = 1.85280e-09 280828. diff --git a/tests/bits/umfpack_06.cc b/tests/bits/umfpack_06.cc new file mode 100644 index 0000000000..682be140f0 --- /dev/null +++ b/tests/bits/umfpack_06.cc @@ -0,0 +1,171 @@ +//---------------------------- umfpack_05.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2002, 2003, 2004, 2005, 2006 by the deal.II authors and Brian Carnes +// +// 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. +// +//---------------------------- umfpack_05.cc --------------------------- + +// test the umfpack sparse direct solver on a mass matrix that is +// slightly modified to make it nonsymmetric. same as umfpack_03, but +// for a BlockSparseMatrix instead of SparseMatrix matrix + +#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 tmp (solution.size()); + tmp = x; + SparseDirectUMFPACK().solve (B, tmp); + x = tmp; + + x -= solution; + deallog << "relative norm distance = " + << x.l2_norm() / solution.l2_norm() + << std::endl; + deallog << "absolute norms = " + << x.l2_norm() << ' ' << solution.l2_norm() + << std::endl; + Assert (x.l2_norm() / solution.l2_norm() < 1e-8, + ExcInternalError()); + } +} + + +int main () +{ + std::ofstream logfile("umfpack_05/output"); + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + test<1> (); + test<2> (); + test<3> (); +} diff --git a/tests/bits/umfpack_06/.cvsignore b/tests/bits/umfpack_06/.cvsignore new file mode 100644 index 0000000000..d196dfb299 --- /dev/null +++ b/tests/bits/umfpack_06/.cvsignore @@ -0,0 +1 @@ +obj.* exe OK output diff --git a/tests/bits/umfpack_06/cmp/generic b/tests/bits/umfpack_06/cmp/generic new file mode 100644 index 0000000000..b9cb7d0d35 --- /dev/null +++ b/tests/bits/umfpack_06/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 = 1.17531e-10 236911. +DEAL::relative norm distance = 0 +DEAL::absolute norms = 2.35062e-10 473822. +DEAL::3d +DEAL::Number of dofs = 1333 +DEAL::relative norm distance = 0 +DEAL::absolute norms = 1.26673e-10 56165.6 +DEAL::relative norm distance = 0 +DEAL::absolute norms = 2.98388e-10 140414. +DEAL::relative norm distance = 0 +DEAL::absolute norms = 5.96775e-10 280828. diff --git a/tests/bits/umfpack_06/cmp/x86_64-unknown-linux-gnu+gcc3.3 b/tests/bits/umfpack_06/cmp/x86_64-unknown-linux-gnu+gcc3.3 new file mode 100644 index 0000000000..92c6c405f8 --- /dev/null +++ b/tests/bits/umfpack_06/cmp/x86_64-unknown-linux-gnu+gcc3.3 @@ -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 = 1.07697e-10 236911. +DEAL::relative norm distance = 0 +DEAL::absolute norms = 2.15395e-10 473822. +DEAL::3d +DEAL::Number of dofs = 1333 +DEAL::relative norm distance = 0 +DEAL::absolute norms = 2.35179e-10 56165.6 +DEAL::relative norm distance = 0 +DEAL::absolute norms = 9.26398e-10 140414. +DEAL::relative norm distance = 0 +DEAL::absolute norms = 1.85280e-09 280828. -- 2.39.5