From f11759ad4875f828f4edaced8751acf44b74c8e8 Mon Sep 17 00:00:00 2001
From: Denis Davydov <davydden@gmail.com>
Date: Thu, 13 Oct 2016 10:39:48 +0200
Subject: [PATCH] add a quick test for umfpack

---
 tests/quick_tests/CMakeLists.txt |   5 ++
 tests/quick_tests/umfpack.cc     | 121 +++++++++++++++++++++++++++++++
 2 files changed, 126 insertions(+)
 create mode 100644 tests/quick_tests/umfpack.cc

diff --git a/tests/quick_tests/CMakeLists.txt b/tests/quick_tests/CMakeLists.txt
index c6ecfc44dc..6338fa0df9 100644
--- a/tests/quick_tests/CMakeLists.txt
+++ b/tests/quick_tests/CMakeLists.txt
@@ -126,6 +126,11 @@ IF (DEAL_II_WITH_LAPACK)
   make_quicktest("lapack" ${_mybuild} "")
 ENDIF()
 
+# Test Umfpack
+IF (DEAL_II_WITH_UMFPACK)
+  make_quicktest("umfpack" ${_mybuild} "")
+ENDIF()
+
 
 
 # A custom test target:
diff --git a/tests/quick_tests/umfpack.cc b/tests/quick_tests/umfpack.cc
new file mode 100644
index 0000000000..3ba1112943
--- /dev/null
+++ b/tests/quick_tests/umfpack.cc
@@ -0,0 +1,121 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// copy-paste from umfpack/umfpack_01 for dim==2
+// test the umfpack sparse direct solver on a mass matrix.
+// test of the transpose as well
+
+#include "../tests.h"
+#include <iostream>
+#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/sparse_matrix.h>
+#include <deal.II/lac/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 (bool transpose = false)
+{
+  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<double> B;
+  B.reinit (sparsity_pattern);
+
+  QGauss<dim> qr (2);
+  MatrixTools::create_mass_matrix (dof_handler, qr, B);
+
+  // 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);
+
+      if (transpose)
+        B.Tvmult (b, solution);
+      else
+        B.vmult (b, solution);
+
+      x = b;
+      SparseDirectUMFPACK().solve (B, x, transpose);
+
+      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("output");
+  deallog.attach(logfile);
+  deallog.threshold_double(1.e-9);
+
+  test<2> ();
+  test<2> (/*transpose =*/ true);
+}
-- 
2.39.5