<ol>
-<li> New: class TimerOutput::Scope does automatic scope based enter/
+<li> New: SparseDirectUMFPACK has long had the ability to work with
+BlockSparseMatrix objects, but couldn't deal with BlockVector objects.
+This is now fixed.
+<br>
+(Wolfgang Bangerth, 2013/04/21)
+</li>
+
+<li> New: Class TimerOutput::Scope does automatic scope based enter/
exit_section of a TimerOutput object.
<br>
(Timo Heister, 2013/04/18)
//---------------------------------------------------------------------------
// $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
~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);
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 <i>Ax = b</i>. 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<double> &, const Vector<double> &) const;
+ void vmult (Vector<double> &dst,
+ const Vector<double> &src) const;
/**
- * Not implemented but necessary for compiling.
+ * Same as before, but for block vectors.
*/
- void Tvmult (Vector<double> &, const Vector<double> &) const;
+ void vmult (BlockVector<double> &dst,
+ const BlockVector<double> &src) const;
+
+ /**
+ * Not implemented but necessary for compiling certain other classes.
+ */
+ void Tvmult (Vector<double> &dst,
+ const Vector<double> &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<double> &, const Vector<double> &) const;
+ void vmult_add (Vector<double> &dst,
+ const Vector<double> &src) const;
/**
- * Not implemented but necessary for compiling.
+ * Not implemented but necessary for compiling certain other classes.
+ */
+ void Tvmult_add (Vector<double> &dst,
+ const Vector<double> &src) const;
+
+ /**
+ * @}
+ */
+
+ /**
+ * @name Functions that solve linear systems
+ */
+ /**
+ * @{
*/
- void Tvmult_add (Vector<double> &, const Vector<double> &) const;
/**
* Solve for a certain right hand side vector. This function may be
void solve (Vector<double> &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<double> &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.
void solve (const Matrix &matrix,
Vector<double> &rhs_and_solution);
+ /**
+ * Same as before, but for block vectors.
+ */
+ template <class Matrix>
+ void solve (const Matrix &matrix,
+ BlockVector<double> &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
}
+void
+SparseDirectUMFPACK::solve (BlockVector<double> &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<double> tmp (rhs_and_solution.size());
+ tmp = rhs_and_solution;
+ solve (tmp);
+ rhs_and_solution = tmp;
+}
+
+
template <class Matrix>
void
}
+template <class Matrix>
+void
+SparseDirectUMFPACK::solve (const Matrix &matrix,
+ BlockVector<double> &rhs_and_solution)
+{
+ factorize (matrix);
+ solve (rhs_and_solution);
+}
+
+
#else
}
+
+void
+SparseDirectUMFPACK::solve (BlockVector<double> &) 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 <class Matrix>
void
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 <class Matrix>
+void
+SparseDirectUMFPACK::solve (const Matrix &,
+ BlockVector<double> &)
+{
+ 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
}
+
+void
+SparseDirectUMFPACK::vmult (
+ BlockVector<double> &dst,
+ const BlockVector<double> &src) const
+{
+ dst = src;
+ this->solve(dst);
+}
+
+
void
SparseDirectUMFPACK::Tvmult (
Vector<double> &,
void SparseDirectUMFPACK::solve (const MATRIX &, \
Vector<double> &); \
template \
+ void SparseDirectUMFPACK::solve (const MATRIX &, \
+ BlockVector<double> &); \
+ template \
void SparseDirectUMFPACK::initialize (const MATRIX &, \
const AdditionalData);
--- /dev/null
+//----------------------------------------------------------------------
+// $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 <iomanip>
+#include <fstream>
+#include <cstdlib>
+
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/function.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/lac/block_sparse_matrix.h>
+#include <deal.II/lac/block_sparsity_pattern.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/sparse_direct.h>
+
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/numerics/matrix_tools.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;
+ 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_11/output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-8);
+
+ 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 = 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.
projects / dealii-svn.git / commitdiff