From e7d80eb1848542aa1ada000473f18f8fa6c6554e Mon Sep 17 00:00:00 2001
From: Wolfgang Bangerth <bangerth@colostate.edu>
Date: Tue, 9 Jan 2024 15:26:22 -0700
Subject: [PATCH] Also test MUMPS symmetric mode.

---
 tests/petsc/sparse_direct_mumps.cc            | 52 +++++++++++++------
 ...ith_mumps=on.with_64bit_indices=off.output |  8 ++-
 2 files changed, 43 insertions(+), 17 deletions(-)

diff --git a/tests/petsc/sparse_direct_mumps.cc b/tests/petsc/sparse_direct_mumps.cc
index 58910e3e1a..a07ad2d6c9 100644
--- a/tests/petsc/sparse_direct_mumps.cc
+++ b/tests/petsc/sparse_direct_mumps.cc
@@ -40,30 +40,52 @@ main(int argc, char **argv)
 
   Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
   {
-    const unsigned int size = 32;
-    unsigned int       dim  = (size - 1) * (size - 1);
+    const unsigned int n_points = 32;
+    unsigned int       size     = (n_points - 1) * (n_points - 1);
 
-    deallog << "Size " << size << " Unknowns " << dim << std::endl;
+    deallog << "Size " << n_points << " Unknowns " << size << std::endl;
 
     // Make matrix
-    FDMatrix                    testproblem(size, size);
-    PETScWrappers::SparseMatrix A(dim, dim, 5);
+    FDMatrix                    testproblem(n_points, n_points);
+    PETScWrappers::SparseMatrix A(size, size, 5);
     testproblem.five_point(A);
 
-    IndexSet indices(dim);
-    indices.add_range(0, dim);
+    IndexSet indices(size);
+    indices.add_range(0, size);
     PETScWrappers::MPI::Vector f(indices, MPI_COMM_WORLD);
     PETScWrappers::MPI::Vector u(indices, MPI_COMM_WORLD);
-    u = 0.;
-    f = 1.;
     A.compress(VectorOperation::insert);
 
-    SolverControl                    cn;
-    PETScWrappers::SparseDirectMUMPS solver(cn);
-    //    solver.set_symmetric_mode(true);
-    solver.solve(A, u, f);
+    // First test non-symmetric mode:
+    {
+      deallog << "Testing non-symmetric mode" << std::endl;
 
-    PETScWrappers::MPI::Vector tmp(indices, MPI_COMM_WORLD);
-    deallog << "residual = " << A.residual(tmp, u, f) << std::endl;
+      u = 0.;
+      f = 1.;
+
+      SolverControl                    cn;
+      PETScWrappers::SparseDirectMUMPS solver(cn);
+      solver.solve(A, u, f);
+
+      PETScWrappers::MPI::Vector tmp(indices, MPI_COMM_WORLD);
+      deallog << "residual = " << A.residual(tmp, u, f) << std::endl;
+    }
+
+    // Now also test the case where the matrix can be assumed to be
+    // symmetric (which of course it is here):
+    {
+      deallog << "Testing symmetric mode" << std::endl;
+
+      u = 0.;
+      f = 1.;
+
+      SolverControl                    cn;
+      PETScWrappers::SparseDirectMUMPS solver(cn);
+      solver.set_symmetric_mode(true);
+      solver.solve(A, u, f);
+
+      PETScWrappers::MPI::Vector tmp(indices, MPI_COMM_WORLD);
+      deallog << "residual = " << A.residual(tmp, u, f) << std::endl;
+    }
   }
 }
diff --git a/tests/petsc/sparse_direct_mumps.with_petsc_with_mumps=on.with_64bit_indices=off.output b/tests/petsc/sparse_direct_mumps.with_petsc_with_mumps=on.with_64bit_indices=off.output
index 932860edcc..f3462ba980 100644
--- a/tests/petsc/sparse_direct_mumps.with_petsc_with_mumps=on.with_64bit_indices=off.output
+++ b/tests/petsc/sparse_direct_mumps.with_petsc_with_mumps=on.with_64bit_indices=off.output
@@ -1,4 +1,8 @@
 
 DEAL::Size 32 Unknowns 961
-DEAL::Convergence step 1 value 0
-DEAL::residual = 0
+DEAL::Testing non-symmetric mode
+DEAL::Convergence step 1 value 0.000
+DEAL::residual = 7.366e-13
+DEAL::Testing symmetric mode
+DEAL::Convergence step 1 value 0.000
+DEAL::residual = 5.412e-13
-- 
2.39.5