From c35fe568d7052abf79e17f4154820df5f54741a3 Mon Sep 17 00:00:00 2001
From: Charlie Vanaret <vanaret@zib.de>
Date: Tue, 19 Nov 2024 14:16:56 +0100
Subject: [PATCH] Trying out other MUMPS parameters

---
 .../SymmetricIndefiniteLinearSystem.hpp           |  6 +++---
 uno/solvers/MUMPS/MUMPSSolver.cpp                 | 15 +++++++++------
 uno/solvers/MUMPS/MUMPSSolver.hpp                 |  2 +-
 3 files changed, 13 insertions(+), 10 deletions(-)

diff --git a/uno/linear_algebra/SymmetricIndefiniteLinearSystem.hpp b/uno/linear_algebra/SymmetricIndefiniteLinearSystem.hpp
index 86f9b741..ff93243c 100644
--- a/uno/linear_algebra/SymmetricIndefiniteLinearSystem.hpp
+++ b/uno/linear_algebra/SymmetricIndefiniteLinearSystem.hpp
@@ -91,7 +91,7 @@ namespace uno {
    template <typename ElementType>
    void SymmetricIndefiniteLinearSystem<ElementType>::factorize_matrix(DirectSymmetricIndefiniteLinearSolver<size_t, ElementType>& linear_solver,
          const WarmstartInformation& warmstart_information) {
-      if (warmstart_information.problem_structure_changed) {
+      if (true) {
          DEBUG << "Performing symbolic factorization of the indefinite system\n";
          linear_solver.do_symbolic_factorization(this->matrix);
       }
@@ -111,7 +111,7 @@ namespace uno {
       DEBUG << "Number of attempts: " << number_attempts << "\n\n";
 
       auto [number_pos_eigenvalues, number_neg_eigenvalues, number_zero_eigenvalues] = linear_solver.get_inertia();
-      DEBUG << "Expected inertia (" << size_primal_block << ", " << size_dual_block << ", 0), ";
+      DEBUG << "Expected inertia  (" << size_primal_block << ", " << size_dual_block << ", 0)\n";
       DEBUG << "Estimated inertia (" << number_pos_eigenvalues << ", " << number_neg_eigenvalues << ", " << number_zero_eigenvalues << ")\n";
 
       if (number_pos_eigenvalues == size_primal_block && number_neg_eigenvalues == size_dual_block && number_zero_eigenvalues == 0) {
@@ -148,7 +148,7 @@ namespace uno {
          DEBUG << "Number of attempts: " << number_attempts << "\n";
 
          std::tie(number_pos_eigenvalues, number_neg_eigenvalues, number_zero_eigenvalues) = linear_solver.get_inertia();
-         DEBUG << "Expected inertia (" << size_primal_block << ", " << size_dual_block << ", 0), ";
+         DEBUG << "Expected inertia  (" << size_primal_block << ", " << size_dual_block << ", 0)\n";
          DEBUG << "Estimated inertia (" << number_pos_eigenvalues << ", " << number_neg_eigenvalues << ", " << number_zero_eigenvalues << ")\n";
 
          if (number_pos_eigenvalues == size_primal_block && number_neg_eigenvalues == size_dual_block && number_zero_eigenvalues == 0) {
diff --git a/uno/solvers/MUMPS/MUMPSSolver.cpp b/uno/solvers/MUMPS/MUMPSSolver.cpp
index 8a754c3d..c19e4648 100644
--- a/uno/solvers/MUMPS/MUMPSSolver.cpp
+++ b/uno/solvers/MUMPS/MUMPSSolver.cpp
@@ -25,9 +25,12 @@ namespace uno {
       dmumps_c(&this->mumps_structure);
       // control parameters
       this->mumps_structure.icntl[0] = -1;
-      this->mumps_structure.icntl[1] = -1;
-      this->mumps_structure.icntl[2] = -1;
+      this->mumps_structure.icntl[1] = 0;
+      this->mumps_structure.icntl[2] = 0;
       this->mumps_structure.icntl[3] = 0;
+      this->mumps_structure.icntl[6] = 7; // pivot order
+      this->mumps_structure.icntl[7] = 6; // scaling
+      this->mumps_structure.icntl[9] = 0; // no iterative refinement
       this->mumps_structure.icntl[12] = 1;
       this->mumps_structure.icntl[23] = 1; // ICNTL(24) controls the detection of “null pivot rows”
    }
@@ -44,11 +47,11 @@ namespace uno {
    }
 
    void MUMPSSolver::do_symbolic_factorization(const SymmetricMatrix<size_t, double>& matrix) {
-      this->save_matrix_to_local_format(matrix);
+      this->mumps_structure.job = MUMPSSolver::JOB_ANALYSIS;
+      this->save_sparsity_to_local_format(matrix);
       this->mumps_structure.n = static_cast<int>(matrix.dimension());
       this->mumps_structure.nnz = static_cast<int>(matrix.number_nonzeros());
-      this->mumps_structure.job = MUMPSSolver::JOB_ANALYSIS;
-      // connect the local COO matrix with the pointers in the structure
+      // connect the local COO sparsity with the pointers in the structure
       this->mumps_structure.irn = this->COO_matrix.row_indices_pointer();
       this->mumps_structure.jcn = this->COO_matrix.column_indices_pointer();
       dmumps_c(&this->mumps_structure);
@@ -92,7 +95,7 @@ namespace uno {
       return this->dimension - this->number_zero_eigenvalues();
    }
 
-   void MUMPSSolver::save_matrix_to_local_format(const SymmetricMatrix<size_t, double>& matrix) {
+   void MUMPSSolver::save_sparsity_to_local_format(const SymmetricMatrix<size_t, double>& matrix) {
       // build the internal matrix representation
       this->COO_matrix.reset();
       for (const auto [row_index, column_index, element]: matrix) {
diff --git a/uno/solvers/MUMPS/MUMPSSolver.hpp b/uno/solvers/MUMPS/MUMPSSolver.hpp
index 83962b36..e79e3f61 100644
--- a/uno/solvers/MUMPS/MUMPSSolver.hpp
+++ b/uno/solvers/MUMPS/MUMPSSolver.hpp
@@ -40,7 +40,7 @@ namespace uno {
       static const int GENERAL_SYMMETRIC = 2;
 
       const size_t fortran_shift{1};
-      void save_matrix_to_local_format(const SymmetricMatrix<size_t, double>& row_index);
+      void save_sparsity_to_local_format(const SymmetricMatrix<size_t, double>& matrix);
    };
 } // namespace