]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add tests.
authorWolfgang Bangerth <bangerth@colostate.edu>
Mon, 30 Dec 2019 19:08:22 +0000 (12:08 -0700)
committerWolfgang Bangerth <bangerth@colostate.edu>
Mon, 6 Jan 2020 19:29:50 +0000 (12:29 -0700)
tests/umfpack/umfpack_03.cc
tests/umfpack/umfpack_03_complex.cc [new file with mode: 0644]
tests/umfpack/umfpack_03_complex.output [new file with mode: 0644]
tests/umfpack/umfpack_03_complex_complex.cc [new file with mode: 0644]
tests/umfpack/umfpack_03_complex_complex.output [new file with mode: 0644]
tests/umfpack/umfpack_09_complex.cc [new file with mode: 0644]
tests/umfpack/umfpack_09_complex.output [new file with mode: 0644]
tests/umfpack/umfpack_09_complex_complex_complex_nonblock.cc [new file with mode: 0644]
tests/umfpack/umfpack_09_complex_complex_complex_nonblock.output [new file with mode: 0644]

index f6bf236c96bd2e11df1ad0aef6dfe6f31f955496..c43739fa11b9560f824ab33f7feab91eed67edf7 100644 (file)
@@ -43,7 +43,7 @@
 
 template <int dim>
 void
-test(bool transpose = false)
+test(const bool transpose = false)
 {
   deallog << dim << 'd' << std::endl;
 
diff --git a/tests/umfpack/umfpack_03_complex.cc b/tests/umfpack/umfpack_03_complex.cc
new file mode 100644 (file)
index 0000000..6b3103f
--- /dev/null
@@ -0,0 +1,149 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2002 - 2018 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.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+// Like the _03 test, but with complex matrix and vector. That said,
+// the matrix just so happens to have zero imaginary components, even
+// though we store its entries as complex numbers.
+
+#include <deal.II/base/function.h>
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/lac/sparse_direct.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/sparsity_pattern.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <iostream>
+
+#include "../tests.h"
+
+
+template <int dim>
+void
+test(const 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<std::complex<double>> B;
+  B.reinit(sparsity_pattern);
+
+  // Create a real sparse matrix and copy it onto B
+  {
+    SparseMatrix<double> BB;
+    BB.reinit(sparsity_pattern);
+
+    QGauss<dim> qr(2);
+    MatrixTools::create_mass_matrix(dof_handler, qr, BB);
+
+    // 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 = BB.begin(); p != BB.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 = BB.begin(); p != BB.end(); ++p)
+      if (p->column() != p->row())
+        AssertThrow(B(p->row(), p->column()) != BB(p->column(), p->row()),
+                    ExcInternalError());
+
+    for (SparseMatrix<double>::iterator p = BB.begin(); p != BB.end(); ++p)
+      B.set(p->row(), p->column(), {p->value(), 0});
+  }
+
+
+  // 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<std::complex<double>> solution(dof_handler.n_dofs());
+      Vector<std::complex<double>> x(dof_handler.n_dofs());
+      Vector<std::complex<double>> b(dof_handler.n_dofs());
+
+      for (unsigned int j = 0; j < dof_handler.n_dofs(); ++j)
+        solution(j) =
+          1. * (j + j * (i + 1) * (i + 1)) *
+          std::complex<double>(1. / std::sqrt(2.), 1. / std::sqrt(2.));
+
+      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()
+{
+  initlog();
+
+  test<1>();
+  test<2>();
+  test<3>();
+
+  test<1>(/*transpose =*/true);
+  test<2>(/*transpose =*/true);
+  test<3>(/*transpose =*/true);
+}
diff --git a/tests/umfpack/umfpack_03_complex.output b/tests/umfpack/umfpack_03_complex.output
new file mode 100644 (file)
index 0000000..5cbb5f8
--- /dev/null
@@ -0,0 +1,49 @@
+
+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.
+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.
diff --git a/tests/umfpack/umfpack_03_complex_complex.cc b/tests/umfpack/umfpack_03_complex_complex.cc
new file mode 100644 (file)
index 0000000..2c75e27
--- /dev/null
@@ -0,0 +1,148 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2002 - 2018 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.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+// Like the _03_complex, but this time with nonzero real and imaginary
+// parts.
+
+#include <deal.II/base/function.h>
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/lac/sparse_direct.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/sparsity_pattern.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <iostream>
+
+#include "../tests.h"
+
+
+template <int dim>
+void
+test(const 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<std::complex<double>> B;
+  B.reinit(sparsity_pattern);
+
+  // Create a real sparse matrix and copy it onto B
+  {
+    SparseMatrix<double> BB;
+    BB.reinit(sparsity_pattern);
+
+    QGauss<dim> qr(2);
+    MatrixTools::create_mass_matrix(dof_handler, qr, BB);
+
+    // 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 = BB.begin(); p != BB.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 = BB.begin(); p != BB.end(); ++p)
+      if (p->column() != p->row())
+        AssertThrow(B(p->row(), p->column()) != BB(p->column(), p->row()),
+                    ExcInternalError());
+
+    for (SparseMatrix<double>::iterator p = BB.begin(); p != BB.end(); ++p)
+      B.set(p->row(), p->column(), {p->value(), 2 * p->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)
+    {
+      Vector<std::complex<double>> solution(dof_handler.n_dofs());
+      Vector<std::complex<double>> x(dof_handler.n_dofs());
+      Vector<std::complex<double>> b(dof_handler.n_dofs());
+
+      for (unsigned int j = 0; j < dof_handler.n_dofs(); ++j)
+        solution(j) =
+          1. * (j + j * (i + 1) * (i + 1)) *
+          std::complex<double>(1. / std::sqrt(5.), 2. / std::sqrt(5.));
+
+      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()
+{
+  initlog();
+
+  test<1>();
+  test<2>();
+  test<3>();
+
+  test<1>(/*transpose =*/true);
+  test<2>(/*transpose =*/true);
+  test<3>(/*transpose =*/true);
+}
diff --git a/tests/umfpack/umfpack_03_complex_complex.output b/tests/umfpack/umfpack_03_complex_complex.output
new file mode 100644 (file)
index 0000000..5cbb5f8
--- /dev/null
@@ -0,0 +1,49 @@
+
+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.
+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.
diff --git a/tests/umfpack/umfpack_09_complex.cc b/tests/umfpack/umfpack_09_complex.cc
new file mode 100644 (file)
index 0000000..d7ed802
--- /dev/null
@@ -0,0 +1,180 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2002 - 2018 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.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+// This test is like the _09 test in that it uses a real-valued
+// matrix, but it then solves a linear system with a complex-valued
+// right hand side and corresponding complex-valued solution vector.
+
+#include <deal.II/base/function.h>
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/lac/block_sparse_matrix.h>
+#include <deal.II/lac/block_sparsity_pattern.h>
+#include <deal.II/lac/sparse_direct.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include "../tests.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());
+  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())
+        AssertThrow(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<std::complex<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();
+
+      // Pick a solution vector. Normalize it in such a way that we
+      // get the exact same norm as for the regular _09 test. Since
+      // the actual choice of rhs shouldn't matter for this test, we
+      // are free to normalize as we see fit.
+      for (unsigned int j = 0; j < dof_handler.n_dofs(); ++j)
+        {
+          solution(j) =
+            std::complex<double>(1 + j + j * (i + 1) * (i + 1), // real part
+                                 1 + i +
+                                   i * (j + 1) * (j + 1)); // imaginary part
+          solution(j) *=
+            1. * (j + j * (i + 1) * (i + 1)) / std::fabs(solution(j));
+        }
+
+      // Then choose as rhs for the linear system the vector
+      //   b = B*solution
+      // so that later the solution of the linear system Bx=b
+      // is again
+      //   x = solution
+      BlockVector<std::complex<double>> x;
+      x.reinit(solution);
+      BlockVector<std::complex<double>> b;
+      b.reinit(solution);
+
+      B.vmult(b, solution);
+      x = b;
+      SparseDirectUMFPACK().solve(B, x);
+
+      // Check that we really got what we expected
+      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()
+{
+  initlog();
+
+  test<1>();
+  test<2>();
+  test<3>();
+}
diff --git a/tests/umfpack/umfpack_09_complex.output b/tests/umfpack/umfpack_09_complex.output
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.
diff --git a/tests/umfpack/umfpack_09_complex_complex_complex_nonblock.cc b/tests/umfpack/umfpack_09_complex_complex_complex_nonblock.cc
new file mode 100644 (file)
index 0000000..0f9f48e
--- /dev/null
@@ -0,0 +1,157 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2002 - 2018 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.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+// This test is like the _09_complex_complex test, but indeed stores a
+// matrix with complex entries.
+
+#include <deal.II/base/function.h>
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/lac/block_sparse_matrix.h>
+#include <deal.II/lac/block_sparsity_pattern.h>
+#include <deal.II/lac/sparse_direct.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include "../tests.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(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<std::complex<double>> B;
+  B.reinit(sparsity_pattern);
+
+  {
+    SparseMatrix<double> BB;
+    BB.reinit(sparsity_pattern);
+
+    QGauss<dim> qr(2);
+    MatrixTools::create_mass_matrix(dof_handler, qr, BB);
+
+    // // 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 = BB.begin(); p != BB.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 = BB.begin(); p != BB.end(); ++p)
+    //   if (p->column() != p->row())
+    //     AssertThrow(BB(p->row(), p->column()) != BB(p->column(), p->row()),
+    //                 ExcInternalError());
+
+    // Now copy stuff over and make sure that what we have is indeed
+    // complex-valued
+    for (auto i = BB.begin(); i != BB.end(); ++i)
+      B.set(i->row(), i->column(), {i->value(), 2. * 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)
+    {
+      Vector<std::complex<double>> solution(dof_handler.n_dofs());
+
+      // Pick a solution vector. Normalize it in such a way that we
+      // get the exact same norm as for the regular _09 test. Since
+      // the actual choice of rhs shouldn't matter for this test, we
+      // are free to normalize as we see fit.
+      for (unsigned int j = 0; j < dof_handler.n_dofs(); ++j)
+        solution(j) =
+          1. * (j + j * (i + 1) * (i + 1)) *
+          std::complex<double>(1. / std::sqrt(5.), 2. / std::sqrt(5.));
+
+      // Then choose as rhs for the linear system the vector
+      //   b = B*solution
+      // so that later the solution of the linear system Bx=b
+      // is again
+      //   x = solution
+      Vector<std::complex<double>> x;
+      x.reinit(solution);
+      Vector<std::complex<double>> b;
+      b.reinit(solution);
+
+      B.Tvmult(b, solution);
+      x = b;
+      SparseDirectUMFPACK().solve(B, x, true);
+
+      // Now check that the residual is zero or close to
+      Vector<std::complex<double>> res;
+      res.reinit(solution);
+      B.residual(res, x, b);
+      Assert(res.l2_norm() <= 1e-8 * b.l2_norm(), ExcInternalError());
+
+      // Check that we really got what we expected
+      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()
+{
+  initlog();
+
+  test<1>();
+  test<2>();
+  test<3>();
+}
diff --git a/tests/umfpack/umfpack_09_complex_complex_complex_nonblock.output b/tests/umfpack/umfpack_09_complex_complex_complex_nonblock.output
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.