<ol>
+ <li> <p> Improved: The <code
+ class="class">SparseDirectUMFPACK</code> 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.
+ <br>
+ (WB 2006/04/25)
+ </p>
+
<li> <p> New: The function <code
class="member">BlockSparsityPattern::row_length</code> adds up
the row lengths of the individual blocks of a block matrix for
<br>
(WB 2006/04/25)
</p>
-
<li> <p> New: There is a new class <code
class="class">IdentityMatrix</code> that represents an
* from the deal.II ReadMe file. UMFPACK is included courtesy of its author,
* <a href="http://www.cise.ufl.edu/~davis/">Timothy A. Davis</a>.
*
- *
+ *
+ * <h4>Instantiations</h4>
+ *
+ * There are instantiations of this class for SparseMatrix<double>,
+ * SparseMatrix<float>, BlockSparseMatrix<double>, and
+ * BlockSparseMatrix<float>.
+ *
* @ingroup Solvers Preconditioners
* @see @ref SoftwareUMFPACK
*
* this operation, even if subsequent
* solves are required.
*/
- void factorize (const SparseMatrix<double> &matrix);
+ template <class Matrix>
+ void factorize (const Matrix &matrix);
/**
* Initialize memory and call
* SparseDirectUMFPACK::factorize.
*/
- void initialize(const SparseMatrix<double> &matrix,
- const AdditionalData = AdditionalData());
+ template <class Matrix>
+ void initialize(const Matrix &matrix,
+ const AdditionalData additional_data = AdditionalData());
/**
* Preconditioner interface
* in place of the right hand
* side vector.
*/
- void solve (const SparseMatrix<double> &matrix,
- Vector<double> &rhs_and_solution);
+ template <class Matrix>
+ void solve (const Matrix &matrix,
+ Vector<double> &rhs_and_solution);
/**
* One of the UMFPack routines
// $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
#include <base/memory_consumption.h>
#include <base/thread_management.h>
#include <lac/sparse_matrix.h>
+#include <lac/block_sparse_matrix.h>
#include <lac/vector.h>
#include <iostream>
+template <class Matrix>
void
SparseDirectUMFPACK::
-factorize (const SparseMatrix<double> &matrix)
+factorize (const Matrix &matrix)
{
Assert (matrix.m() == matrix.n(), ExcNotQuadratic())
Assert (static_cast<unsigned int>(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<double>::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<int> 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<Ap.size()-1; ++i)
+ Assert (row_pointers[i] == Ap[i+1], ExcInternalError());
}
// finally do the copying around of entries
}
}
}
-
-
int status;
+template <class Matrix>
void
-SparseDirectUMFPACK::solve (const SparseMatrix<double> &matrix,
- Vector<double> &rhs_and_solution)
+SparseDirectUMFPACK::solve (const Matrix &matrix,
+ Vector<double> &rhs_and_solution)
{
factorize (matrix);
solve (rhs_and_solution);
{}
-void SparseDirectUMFPACK::factorize (const SparseMatrix<double> &)
+template <class Matrix>
+void SparseDirectUMFPACK::factorize (const Matrix &)
{
Assert(false, ExcNeedsUMFPACK());
}
}
+template <class Matrix>
void
-SparseDirectUMFPACK::solve (const SparseMatrix<double> &,
- Vector<double> &)
+SparseDirectUMFPACK::solve (const Matrix &,
+ Vector<double> &)
{
Assert(false, ExcNeedsUMFPACK());
}
#endif
+template <class Matrix>
void
-SparseDirectUMFPACK::initialize (const SparseMatrix<double>& M,
+SparseDirectUMFPACK::initialize (const Matrix &M,
const AdditionalData)
{
this->factorize(M);
}
-// explicit instantiations
+// explicit instantiations for SparseMatrixMA27
template
-void
-SparseDirectMA27::factorize (const SparseMatrix<double> &matrix);
+void SparseDirectMA27::factorize (const SparseMatrix<double> &matrix);
template
-void
-SparseDirectMA27::factorize (const SparseMatrix<float> &matrix);
+void SparseDirectMA27::factorize (const SparseMatrix<float> &matrix);
template
-void
-SparseDirectMA27::solve (const SparseMatrix<double> &matrix,
+void SparseDirectMA27::solve (const SparseMatrix<double> &matrix,
Vector<double> &rhs_and_solution);
template
-void
-SparseDirectMA27::solve (const SparseMatrix<float> &matrix,
- Vector<double> &rhs_and_solution);
-
+void SparseDirectMA27::solve (const SparseMatrix<float> &matrix,
+ Vector<double> &rhs_and_solution);
+
+
+// explicit instantiations for SparseMatrixUMFPACK
+#define InstantiateUMFPACK(MATRIX) \
+ template \
+ void SparseDirectUMFPACK::factorize (const MATRIX &); \
+ template \
+ void SparseDirectUMFPACK::solve (const MATRIX &, \
+ Vector<double> &); \
+ template \
+ void SparseDirectUMFPACK::initialize (const MATRIX &, \
+ const AdditionalData)
+
+InstantiateUMFPACK(SparseMatrix<double>);
+InstantiateUMFPACK(SparseMatrix<float>);
+InstantiateUMFPACK(BlockSparseMatrix<double>);
+InstantiateUMFPACK(BlockSparseMatrix<float>);
--- /dev/null
+//---------------------------- 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<float> instead of SparseMatrix<double> matrix
+
+#include "../tests.h"
+#include <iostream>
+#include <fstream>
+#include <cstdlib>
+
+#include <base/quadrature_lib.h>
+#include <base/function.h>
+
+#include <fe/fe_q.h>
+#include <fe/fe_values.h>
+
+#include <dofs/dof_tools.h>
+
+#include <lac/sparse_matrix.h>
+#include <lac/sparsity_pattern.h>
+#include <lac/vector.h>
+#include <lac/sparse_direct.h>
+
+#include <grid/tria.h>
+#include <grid/grid_generator.h>
+
+#include <numerics/vectors.h>
+#include <numerics/matrices.h>
+
+
+template <int dim>
+void test ()
+{
+ deallog << dim << 'd' << std::endl;
+
+ Triangulation<dim> 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<dim> fe (1);
+ DoFHandler<dim> 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<float> B;
+ B.reinit (sparsity_pattern);
+
+ QGauss<dim> 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<float>::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<float>::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<double> solution (dof_handler.n_dofs());
+ Vector<double> x (dof_handler.n_dofs());
+ Vector<double> b (dof_handler.n_dofs());
+
+ for (unsigned int j=0; j<dof_handler.n_dofs(); ++j)
+ solution(j) = j+j*(i+1)*(i+1);
+
+ B.vmult (b, solution);
+ x = b;
+ SparseDirectUMFPACK().solve (B, x);
+
+ 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_04/output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ test<1> ();
+ test<2> ();
+ test<3> ();
+}
--- /dev/null
+
+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.
--- /dev/null
+//---------------------------- 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<float> instead of SparseMatrix<double> matrix
+
+#include "../tests.h"
+#include <iostream>
+#include <fstream>
+#include <cstdlib>
+
+#include <base/quadrature_lib.h>
+#include <base/function.h>
+
+#include <fe/fe_q.h>
+#include <fe/fe_values.h>
+
+#include <dofs/dof_tools.h>
+
+#include <lac/block_sparse_matrix.h>
+#include <lac/block_sparsity_pattern.h>
+#include <lac/vector.h>
+#include <lac/sparse_direct.h>
+
+#include <grid/tria.h>
+#include <grid/grid_generator.h>
+
+#include <numerics/vectors.h>
+#include <numerics/matrices.h>
+
+
+template <int dim>
+void test ()
+{
+ deallog << dim << 'd' << std::endl;
+
+ Triangulation<dim> 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<dim> fe (1);
+ DoFHandler<dim> 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<float> 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<double> xB;
+ xB.reinit (xsparsity_pattern);
+
+ QGauss<dim> 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<double>::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<double>::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<double>::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<double> 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<double> x;
+ x.reinit (solution);
+ BlockVector<double> b;
+ b.reinit (solution);
+
+ for (unsigned int j=0; j<dof_handler.n_dofs(); ++j)
+ solution(j) = j+j*(i+1)*(i+1);
+
+ B.vmult (b, solution);
+ x = b;
+ Vector<double> 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> ();
+}
--- /dev/null
+
+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.
--- /dev/null
+//---------------------------- 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<double> instead of SparseMatrix<double> matrix
+
+#include "../tests.h"
+#include <iostream>
+#include <fstream>
+#include <cstdlib>
+
+#include <base/quadrature_lib.h>
+#include <base/function.h>
+
+#include <fe/fe_q.h>
+#include <fe/fe_values.h>
+
+#include <dofs/dof_tools.h>
+
+#include <lac/block_sparse_matrix.h>
+#include <lac/block_sparsity_pattern.h>
+#include <lac/vector.h>
+#include <lac/sparse_direct.h>
+
+#include <grid/tria.h>
+#include <grid/grid_generator.h>
+
+#include <numerics/vectors.h>
+#include <numerics/matrices.h>
+
+
+template <int dim>
+void test ()
+{
+ deallog << dim << 'd' << std::endl;
+
+ Triangulation<dim> 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<dim> fe (1);
+ DoFHandler<dim> 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<double> 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<double> xB;
+ xB.reinit (xsparsity_pattern);
+
+ QGauss<dim> 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<double>::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<double>::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<double>::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<double> 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<double> x;
+ x.reinit (solution);
+ BlockVector<double> b;
+ b.reinit (solution);
+
+ for (unsigned int j=0; j<dof_handler.n_dofs(); ++j)
+ solution(j) = j+j*(i+1)*(i+1);
+
+ B.vmult (b, solution);
+ x = b;
+ Vector<double> 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> ();
+}
--- /dev/null
+
+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.