]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Make SparseDirectUMFPACK also understand block vectors.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Sun, 21 Apr 2013 16:48:46 +0000 (16:48 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Sun, 21 Apr 2013 16:48:46 +0000 (16:48 +0000)
git-svn-id: https://svn.dealii.org/trunk@29354 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/doc/news/changes.h
deal.II/include/deal.II/lac/sparse_direct.h
deal.II/source/lac/sparse_direct.cc
tests/umfpack/umfpack_11.cc [new file with mode: 0644]
tests/umfpack/umfpack_11/cmp/generic [new file with mode: 0644]

index 7fcd24d2e4e91b11af82d649ef6a74e4529871bb..992cfe9e074114366d59c6e512e33918f4763426 100644 (file)
@@ -111,7 +111,14 @@ this function.
 
 <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)
index 195ca99190f4a6288b43c80acacd575987bd790e..da58687d0b772aa974f0ffa22ca3add7ec8b30fe 100644 (file)
@@ -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 <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
@@ -175,7 +217,12 @@ public:
   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.
@@ -184,6 +231,17 @@ public:
   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
index 8099df44d30280d84788bb0e3e6d4771da8d70d1..4df42981df271c4b41345fad50a704a0f7c1d3cb 100644 (file)
@@ -314,6 +314,19 @@ SparseDirectUMFPACK::solve (Vector<double> &rhs_and_solution) const
 }
 
 
+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
@@ -325,6 +338,16 @@ SparseDirectUMFPACK::solve (const Matrix   &matrix,
 }
 
 
+template <class Matrix>
+void
+SparseDirectUMFPACK::solve (const Matrix   &matrix,
+                            BlockVector<double> &rhs_and_solution)
+{
+  factorize (matrix);
+  solve (rhs_and_solution);
+}
+
+
 #else
 
 
@@ -355,6 +378,14 @@ SparseDirectUMFPACK::solve (Vector<double> &) const
 }
 
 
+
+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 &,
@@ -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 <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
 
 
@@ -385,6 +426,17 @@ SparseDirectUMFPACK::vmult (
 }
 
 
+
+void
+SparseDirectUMFPACK::vmult (
+  BlockVector<double>       &dst,
+  const BlockVector<double> &src) const
+{
+  dst = src;
+  this->solve(dst);
+}
+
+
 void
 SparseDirectUMFPACK::Tvmult (
   Vector<double> &,
@@ -596,6 +648,9 @@ void SparseDirectMUMPS::vmult (Vector<double>       &dst,
   void SparseDirectUMFPACK::solve (const MATRIX   &,    \
                                    Vector<double> &);    \
   template    \
+  void SparseDirectUMFPACK::solve (const MATRIX   &,    \
+                                   BlockVector<double> &);    \
+  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 (file)
index 0000000..bbf93cb
--- /dev/null
@@ -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 <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> ();
+}
diff --git a/tests/umfpack/umfpack_11/cmp/generic b/tests/umfpack/umfpack_11/cmp/generic
new file mode 100644 (file)
index 0000000..8df3f88
--- /dev/null
@@ -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.

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.