]> https://gitweb.dealii.org/ - dealii.git/commitdiff
SparseDirectUmfpack now also works for sparse matrices with floats, as well as for...
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Tue, 25 Apr 2006 23:12:39 +0000 (23:12 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Tue, 25 Apr 2006 23:12:39 +0000 (23:12 +0000)
git-svn-id: https://svn.dealii.org/trunk@12897 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/doc/news/changes.html
deal.II/lac/include/lac/sparse_direct.h
deal.II/lac/source/sparse_direct.cc
tests/bits/umfpack_04.cc [new file with mode: 0644]
tests/bits/umfpack_04/cmp/generic [new file with mode: 0644]
tests/bits/umfpack_05.cc [new file with mode: 0644]
tests/bits/umfpack_05/cmp/generic [new file with mode: 0644]
tests/bits/umfpack_06.cc [new file with mode: 0644]
tests/bits/umfpack_06/cmp/generic [new file with mode: 0644]

index 3ae3abe42c09a4043275da7fd4dac07b93db210f..f4d04976575158e01f930ffc02c173f253ca3885 100644 (file)
@@ -442,6 +442,14 @@ inconvenience this causes.
 
 <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
@@ -449,7 +457,6 @@ inconvenience this causes.
        <br> 
        (WB 2006/04/25)
        </p>
-
   
   <li> <p> New: There is a new class <code
        class="class">IdentityMatrix</code> that represents an
index 3bfa14328bd49297e0ffccafbd33097975dbacbc..70b136bdd17002ecfdcd65054a50247b67893921 100644 (file)
@@ -1027,7 +1027,13 @@ class SparseDirectMA47 : public Subscriptor
  * 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
  *
@@ -1093,14 +1099,16 @@ class SparseDirectUMFPACK : public Subscriptor
                                      * 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
@@ -1168,8 +1176,9 @@ class SparseDirectUMFPACK : public Subscriptor
                                      * 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
index 9534cee1a3bc81f8df24f2ef418049cd7d724d72..0a3384219f2c722ed2c9be93751bb25a4a5a9778 100644 (file)
@@ -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 <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>
@@ -1679,9 +1680,10 @@ SparseDirectUMFPACK::clear ()
 
 
 
+template <class Matrix>
 void
 SparseDirectUMFPACK::
-factorize (const SparseMatrix<double> &matrix)
+factorize (const Matrix &matrix)
 {
   Assert (matrix.m() == matrix.n(), ExcNotQuadratic())
   
@@ -1729,16 +1731,37 @@ factorize (const SparseMatrix<double> &matrix)
   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
@@ -1784,8 +1807,6 @@ factorize (const SparseMatrix<double> &matrix)
           }
       }
   }
-  
-  
         
   int status;
 
@@ -1835,9 +1856,10 @@ SparseDirectUMFPACK::solve (Vector<double> &rhs_and_solution) const
 
 
 
+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);
@@ -1860,7 +1882,8 @@ SparseDirectUMFPACK::clear ()
 {}
 
 
-void SparseDirectUMFPACK::factorize (const SparseMatrix<double> &)
+template <class Matrix>
+void SparseDirectUMFPACK::factorize (const Matrix &)
 {
   Assert(false, ExcNeedsUMFPACK());
 }
@@ -1873,9 +1896,10 @@ SparseDirectUMFPACK::solve (Vector<double> &) const
 }
 
 
+template <class Matrix>
 void
-SparseDirectUMFPACK::solve (const SparseMatrix<double> &,
-                            Vector<double>             &)
+SparseDirectUMFPACK::solve (const Matrix   &,
+                            Vector<double> &)
 {
   Assert(false, ExcNeedsUMFPACK());
 }
@@ -1884,8 +1908,9 @@ SparseDirectUMFPACK::solve (const SparseMatrix<double> &,
 #endif
 
 
+template <class Matrix>
 void
-SparseDirectUMFPACK::initialize (const SparseMatrix<double>& 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<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>);
diff --git a/tests/bits/umfpack_04.cc b/tests/bits/umfpack_04.cc
new file mode 100644 (file)
index 0000000..8d605f4
--- /dev/null
@@ -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<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> ();
+}
diff --git a/tests/bits/umfpack_04/cmp/generic b/tests/bits/umfpack_04/cmp/generic
new file mode 100644 (file)
index 0000000..b9cb7d0
--- /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 = 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.cc b/tests/bits/umfpack_05.cc
new file mode 100644 (file)
index 0000000..09436b1
--- /dev/null
@@ -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<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> ();
+}
diff --git a/tests/bits/umfpack_05/cmp/generic b/tests/bits/umfpack_05/cmp/generic
new file mode 100644 (file)
index 0000000..b9cb7d0
--- /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 = 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.cc b/tests/bits/umfpack_06.cc
new file mode 100644 (file)
index 0000000..682be14
--- /dev/null
@@ -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<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> ();
+}
diff --git a/tests/bits/umfpack_06/cmp/generic b/tests/bits/umfpack_06/cmp/generic
new file mode 100644 (file)
index 0000000..b9cb7d0
--- /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 = 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.

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.