Reviewing the MA27 interface (in progress): fixed two bugs, which sol…

…ves the issue with byrd and hs015
cvanaret · Dec 1, 2024 · d319c89 · d319c89
1 parent 239b6c8
commit d319c89
Show file tree

Hide file tree

Showing 2 changed files with 14 additions and 12 deletions.
diff --git a/uno/solvers/MA27/MA27Solver.cpp b/uno/solvers/MA27/MA27Solver.cpp
@@ -111,11 +111,13 @@ namespace uno {
    MA27Solver::MA27Solver(size_t max_dimension, size_t max_number_nonzeros):
          DirectSymmetricIndefiniteLinearSolver<size_t, double>(max_dimension),
          nz_max(static_cast<int>(max_number_nonzeros)), n(static_cast<int>(max_dimension)), nnz(static_cast<int>(max_number_nonzeros)),
-         irn(max_number_nonzeros), icn(max_number_nonzeros), iw((2 * max_number_nonzeros + 3 * max_dimension + 1) * 6 / 5),
-         ikeep(3 * max_number_nonzeros), iw1(2 * max_dimension) {
-      iflag = 0;
-      // set the default values of the controlling parameters
+         irn(max_number_nonzeros), icn(max_number_nonzeros),
+         iw((2 * max_number_nonzeros + 3 * max_dimension + 1) * 6 / 5), // 20% more than 2*nnz + 3*n + 1
+         ikeep(3 * max_dimension), iw1(2 * max_dimension) {
+      // initialization: set the default values of the controlling parameters
       MA27ID(icntl.data(), cntl.data());
+      // a suitable pivot order is to be chosen automatically
+      iflag = 0;
       // suppress warning messages
       icntl[eICNTL::LP] = 0;
       icntl[eICNTL::MP] = 0;
@@ -144,14 +146,11 @@ namespace uno {
       MA27AD(&n, &nnz,                                   /* size info */
             irn.data(), icn.data(),                     /* matrix indices */
             iw.data(), &liw, ikeep.data(), iw1.data(),  /* solver workspace */
-            &nsteps, &iflag,
-            icntl.data(), cntl.data(),
-            info.data(), &ops);
+            &nsteps, &iflag, icntl.data(), cntl.data(), info.data(), &ops);
 
+      // resize the factor by at least INFO(5) (here, 50% more)
       factor.resize(static_cast<size_t>(3 * info[eINFO::NRLNEC] / 2));
 
-      std::copy(matrix.data_pointer(), matrix.data_pointer() + matrix.number_nonzeros(), factor.begin());
-
       assert(eIFLAG::SUCCESS == info[eINFO::IFLAG] && "MA27: the symbolic factorization failed");
       if (eIFLAG::SUCCESS != info[eINFO::IFLAG]) {
          WARNING << "MA27 has issued a warning: IFLAG = " << info[eINFO::IFLAG] << " additional info, IERROR = " << info[eINFO::IERROR] << '\n';
@@ -186,6 +185,9 @@ namespace uno {
       assert(matrix.dimension() <= iw1.capacity() && "MA27Solver: the dimension of the matrix is larger than the preallocated size");
       assert(nnz == static_cast<int>(matrix.number_nonzeros()) && "MA27Solver: the numbers of nonzeros do not match");
 
+      // initialize factor with the entries of the matrix. It will be modified by MA27BD
+      std::copy(matrix.data_pointer(), matrix.data_pointer() + matrix.number_nonzeros(), factor.begin());
+
       // numerical factorization
       int la = static_cast<int>(factor.size());
       int liw = static_cast<int>(iw.size());

diff --git a/uno/solvers/MA27/MA27Solver.hpp b/uno/solvers/MA27/MA27Solver.hpp
@@ -37,8 +37,8 @@ class MA27Solver
    int nz_max{};                    // maximal number of nonzeros entries
    int n{};			                // dimension of current factorisation (maximal value here is <= max_dimension)
    int nnz{};			            // number of nonzeros of current factorisation
-   std::array<int,30> icntl{};      // integer array of length 30; integer control values
-   std::array<double,5> cntl{};     // double array of length 5; double control values
+   std::array<int, 30> icntl{};      // integer array of length 30; integer control values
+   std::array<double, 5> cntl{};     // double array of length 5; double control values
 
    std::vector<int> irn{};          // row index of input
    std::vector<int> icn{};          // col index of input
@@ -48,7 +48,7 @@ class MA27Solver
    std::vector<int> iw1{};          // integer workspace array of length n
    int nsteps{};                    // integer, to be set by ma27
    int iflag{};                     // integer; 0 if pivot order chosen automatically; 1 if the pivot order set by ikeep
-   std::array<int,20> info{};       // integer array of length 20
+   std::array<int, 20> info{};       // integer array of length 20
    double ops{};                    // double, operations count
 
    std::vector<double> factor{};    // data array of length la;