From ea1e8e3501c7e36d42a25f2d99aaae6afee44753 Mon Sep 17 00:00:00 2001 From: Manuel Schaich Date: Sat, 30 Nov 2024 11:22:51 +0000 Subject: [PATCH 01/12] ma27 integration --- uno/solvers/MA27/MA27Solver.cpp | 259 ++++++++++++++++++ uno/solvers/MA27/MA27Solver.hpp | 60 ++++ ...SymmetricIndefiniteLinearSolverFactory.cpp | 20 ++ 3 files changed, 339 insertions(+) create mode 100644 uno/solvers/MA27/MA27Solver.cpp create mode 100644 uno/solvers/MA27/MA27Solver.hpp diff --git a/uno/solvers/MA27/MA27Solver.cpp b/uno/solvers/MA27/MA27Solver.cpp new file mode 100644 index 00000000..b099bac0 --- /dev/null +++ b/uno/solvers/MA27/MA27Solver.cpp @@ -0,0 +1,259 @@ + +#include +#include "MA27Solver.hpp" +#include "linear_algebra/SymmetricMatrix.hpp" +#include "linear_algebra/Vector.hpp" +#include "tools/Logger.hpp" +#include "fortran_interface.h" + +#define MA27ID FC_GLOBAL(ma27id, MA27ID) +#define MA27AD FC_GLOBAL(ma27ad, MA27AD) +#define MA27BD FC_GLOBAL(ma27bd, MA27BD) +#define MA27CD FC_GLOBAL(ma27cd, MA27CD) + +extern "C" +{ + void MA27ID(int ICNTL[], double CNTL[]); + + void MA27AD(int *N, int *NZ, int IRN[], int ICN[], int IW[], int *LIW, int IKEEP[], int IW1[], + int *NSTEPS, int *IFLAG, int ICNTL[], double CNTL[], int INFO[], double *OPS); + + void MA27BD(int *N, int *NZ, int IRN[], int ICN[], double A[], int *LA, int IW[], int *LIW, + int IKEEP[], int *NSTEPS, int *MAXFRT, int IW1[], int ICNTL[], double CNTL[], int INFO[]); + + void MA27CD(int *N, double A[], int *LA, int IW[], int *LIW, double W[], int *MAXFRT, double RHS[], + int IW1[], int *NSTEPS, int ICNTL[], int INFO[]); +} + +namespace uno +{ +constexpr auto fortran_shift = 1; + +enum ICNTL { + LP = 0, // sed by the subroutines as the output stream for error messages. If it is set to zero these messages will be suppressed. The default value is 6. + MP, //used by the subroutines as the output stream for diagnostic printing and for warning messages. If it is set to zero then messages are suppressed. The default value is 6. + LDIAG, // used by the subroutines to control diagnostic printing. If ICNTL(3) is equal to zero (the default), no diagnostic printing will be produced, a value of 1 will print scalar parameters (both in argument lists and in the control and information arrays) and a few entries of array parameters on entry and successful exit from each subroutine while ICNTL(3) equal to 2 will print all parameter values on entry and successful exit. + /* The entries ICNTL(4) to ICNTL(25) are not of interest to the general user and are discussed more fully by Duff and Reid (AERE R-10533, 1982) under the internal names IOVFLO, NEMIN and IFRLVL + */ + IOVFLO, + NEMIN, + IFRLVL1, + IFRLVL2, + IFRLVL3, + IFRLVL4, + IFRLVL5, + IFRLVL6, + IFRLVL7, + IFRLVL8, + IFRLVL9, + IFRLVL10, + IFRLVL11, + IFRLVL12, + IFRLVL13, + IFRLVL14, + IFRLVL15, + IFRLVL16, + IFRLVL17, + IFRLVL18, + IFRLVL19, + IFRLVL20, + UNUSED_ICNTL1, + UNUSED_ICNTL2, + UNUSED_ICNTL3, + UNUSED_ICNTL4, + UNUSED_ICNTL5, +}; + +enum CNTL { + U = 0, // used by the subroutine to control numerical pivoting. Values greater than 0.5 are treated as 0.5 and less than –0.5 as –0.5. Its default value is 0.1. If U is positive, numerical pivoting will be performed. If U is non-positive, no pivoting will be performed, the subroutine will fail if a zero pivot is encountered, and a flag (see section 2.3) will be set if not all pivots are of the same sign; the factorization will continue after a sign change is detected if U is zero but will exit immediately if U is less than zero. If the system is definite, then setting U to zero will decrease the factorization time while still providing a stable decomposition. For problems requiring greater than average numerical care a higher value than the default would be advisable. + FRATIO, // given the default value of 1.0 by MA27I/ID. If MA27A/AD encounters a row of the reduced matrix with a proportion of entries greater than FRATIO, the row is treated as full. FRATIO is not altered by MA27. + PIVTOL, // given the default value of 0.0 by MA27I/ID. MA27B/BD will not accept an entry with absolute value less than PIVTOL as a 1×1 pivot or the off-diagonal entry of a 2×2 pivot. PIVTOL is not altered by MA27. + UNUSED_CNTL1, + UNUSED_CNTL2, +}; + +enum INFO { + IFLAG = 0, // is an error flag. A value of zero indicates that the subroutine has performed successfully. + IERROR, // provides supplementary information when there is an error. + NRLTOT, // gives the total amount of REAL words required for a successful completion of MA27B/BD without the need for data compression provided no numerical pivoting is performed. The actual amount required may be higher because of numerical pivoting, but probably not by more than 3%. + NIRTOT, // gives the total amount of INTEGER words required for a successful completion of MA27B/BD without the need for data compression provided no numerical pivoting is performed. The actual amount required may be higher because of numerical pivoting, but probably not by more than 3%. + NRLNEC, // gives the amount of REAL words required for successful completion of MA27B/BD allowing data compression (see NCMPBR returned in INFO(12)), again provided no numerical pivoting is performed. Numerical pivoting may cause a higher value to be required, but probably not by more than 3%. If storage was conserved by equivalencing IW(1) with IRN(1), NRLNEC and NIRNEC cannot be calculated exactly but instead an upper bound will be returned. Experience has shown that this can overestimate the exact values by 50% although the tightness of the bound is very problem dependent. For example, a tight bound will generally be obtained if there are many more entries in the factors than in the input matrix. + NIRNEC, // gives the amount of INTEGER words required for successful completion of MA27B/BD allowing data compression (see NCMPBR returned in INFO(12)), again provided no numerical pivoting is performed. Numerical pivoting may cause a higher value to be required, but probably not by more than 3%. If storage was conserved by equivalencing IW(1) with IRN(1), NRLNEC and NIRNEC cannot be calculated exactly but instead an upper bound will be returned. Experience has shown that this can overestimate the exact values by 50% although the tightness of the bound is very problem dependent. For example, a tight bound will generally be obtained if there are many more entries in the factors than in the input matrix. + NRLADU, // gives the number of REAL words required to hold the matrix factors if no numerical pivoting is performed by MA27B/BD. Numerical pivoting may change this slightly. + NIRADU, // gives the number of INTEGER words required to hold the matrix factors if no numerical pivoting is performed by MA27B/BD. Numerical pivoting may change this slightly. + NRLBDU, // give the amount of REAL words actually used to hold the factorization. + NIRBDU, // give the amount of INTEGER words actually used to hold the factorization. + NCMPA, // holds the number of compresses of the internal data structure performed by MA27A/AD. If this is high (say > 10), the performance of MA27A/AD may be improved by increasing the length of array IW. + NCMPBR, // holds the number of compresses of the real data structure required by the factorization. If either of these is high (say > 10), then the speed of the factorization may be increased by allocating more space to the arrays A as appropriate. + NCMPBI, // holds the number of compresses of the integer data structure required by the factorization. If either of these is high (say > 10), then the speed of the factorization may be increased by allocating more space to the arrays IW as appropriate. + NTWO, // gives the number of 2×2 pivots used during the factorization. + NEIG, // gives the number of negative eigenvalues of A. + UNUSED_INFO1, + UNUSED_INFO2, + UNUSED_INFO3, + UNUSED_INFO4, + UNUSED_INFO5, +}; + +enum IFLAG { + NSTEPS = -7, // Value of NSTEPS outside the range 1 ≤ NSTEPS ≤ N (MA27B/BD entry). + PIVOTSIGN = -6, // A change of sign of pivots has been detected when U was negative. INFO(2) is set to the pivot step at which the change was detected. (MA27B/BD entry only) + SINGULAR = -5, // Matrix is singular (MA27B/BD entry only). INFO(2) is set to the pivot step at which singularity was detected + INSUFFICIENTREAL = -4, // Failure due to insufficient space allocated to array A (MA27B/BD entry only). INFO(2) is set to a value that may suffice. + INSUFFICIENTINTEGER = -3, // Failure due to insufficient space allocated to array IW (MA27A/AD and MA27B/BD entries). INFO(2) is set to a value that may suffice. + NZOUTOFRANGE = -2, // Value of NZ out of range. NZ < 0. (MA27A/AD and MA27B/BD entries) + NOUTOFRANGE = -1, // Value of N out of range. N < 1. (MA27A/AD and MA27B/BD entries). + SUCCESS = 0, // Successful completion. + IDXOUTOFRANGE = 1, // ndex (in IRN or ICN) out of range. Action taken by subroutine is to ignore any such entries and continue (MA27A/AD and MA27B/BD entries). INFO(2) is set to the number of faulty entries. Details of the first ten are printed on unit ICNTL(2). + FALSEDEFINITENESS, // Pivots have different signs when factorizing a supposedly definite matrix (when the value of U in CNTL(1) is zero) (MA27B/BD entry only). INFO(2) is set to the number of sign changes. Note that this warning will overwrite an INFO(1)=1 warning. Details of the first ten are printed on unit ICNTL(2). + RANKDEFECT, // Matrix is rank deficient. In this case, a decomposition will still have been produced which will enable the subsequent solution of consistent equations (MA27B/BD entry only). INFO(2) will be set to the rank of the matrix. Note that this warning will overwrite an INFO(1)=1 or INFO(1)=2 warning. +}; + + + MA27Solver::MA27Solver(size_t max_dimension, size_t max_number_nonzeros) + : DirectSymmetricIndefiniteLinearSolver(max_dimension) + , nz_max(static_cast(max_number_nonzeros)) + , n(static_cast(max_dimension)) + , nnz(static_cast(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 + MA27ID(icntl.data(), cntl.data()); + // suppress warning messages + icntl[ICNTL::LP] = 0; + icntl[ICNTL::MP] = 0; + icntl[ICNTL::LDIAG] = 0; + } + + void MA27Solver::factorize(const SymmetricMatrix &matrix) + { + // // general factorization method: symbolic factorization and numerical factorization + do_symbolic_factorization(matrix); + do_numerical_factorization(matrix); + } + + void MA27Solver::do_symbolic_factorization(const SymmetricMatrix &matrix) + { + assert(matrix.dimension() <= iw1.capacity() && "MA27Solver: the dimension of the matrix is larger than the preallocated size"); + assert(matrix.number_nonzeros() <= irn.capacity() && + "MA27Solver: the number of nonzeros of the matrix is larger than the preallocated size"); + + // build the internal matrix representation + save_matrix_to_local_format(matrix); + + n = static_cast(matrix.dimension()); + nnz = static_cast(matrix.number_nonzeros()); + + // symbolic factorization + int liw = static_cast(iw.size()); + 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); + + factor.resize(static_cast(3 * info[INFO::NRLNEC] / 2)); + + std::copy(matrix.data_pointer(), matrix.data_pointer() + matrix.number_nonzeros(), factor.begin()); + + assert(IFLAG::SUCCESS == info[INFO::IFLAG] && "MA27: the symbolic factorization failed"); + if (IFLAG::SUCCESS != info[INFO::IFLAG] ) + { + WARNING << "MA27 has issued a warning: IFLAG = " << info[INFO::IFLAG] << " additional info, IERROR = " << info[INFO::IERROR] << '\n'; + } + } + + void MA27Solver::do_numerical_factorization([[maybe_unused]]const SymmetricMatrix &matrix) + { + assert(matrix.dimension() <= iw1.capacity() && "MA27Solver: the dimension of the matrix is larger than the preallocated size"); + assert(nnz == static_cast(matrix.number_nonzeros()) && "MA27Solver: the numbers of nonzeros do not match"); + + // numerical factorization + int la = static_cast(factor.size()); + int liw = static_cast(iw.size()); + MA27BD(&n, &nnz, irn.data(), icn.data(), factor.data(), &la, iw.data(), &liw, + ikeep.data(), &nsteps, &maxfrt, iw1.data(), icntl.data(), cntl.data(), info.data()); + + assert(IFLAG::SUCCESS == info[INFO::IFLAG] && "MA27: the numerical factorization failed"); + if (IFLAG::SUCCESS != info[INFO::IFLAG] ) + { + WARNING << "MA27 has issued a warning: IFLAG = " << info[INFO::IFLAG] << " additional info, IERROR = " << info[INFO::IERROR] << '\n'; + } + } + + void MA27Solver::solve_indefinite_system([[maybe_unused]]const SymmetricMatrix &matrix, const Vector &rhs, Vector &result) + { + // solve + std::vector w(maxfrt); // double workspace + int la = static_cast(factor.size()); + int liw = static_cast(iw.size()); + + result = rhs; + + MA27CD( &n, + factor.data(), + &la, + iw.data(), + &liw, + w.data(), + &maxfrt, + result.data(), + iw1.data(), + &nsteps, + icntl.data(), + info.data()); + + assert(info[INFO::IFLAG] == IFLAG::SUCCESS && "MA27: the solution failed"); + if (IFLAG::SUCCESS != info[INFO::IFLAG] ) + { + WARNING << "MA27 has issued a warning: IFLAG = " << info[INFO::IFLAG] << " additional info, IERROR = " << info[INFO::IERROR] << '\n'; + } + } + + std::tuple MA27Solver::get_inertia() const + { + // rank = number_positive_eigenvalues + number_negative_eigenvalues + // n = rank + number_zero_eigenvalues + const size_t rankA = rank(); + const size_t num_negative_eigenvalues = number_negative_eigenvalues(); + const size_t num_positive_eigenvalues = rankA - num_negative_eigenvalues; + const size_t num_zero_eigenvalues = static_cast(n) - rankA; + return std::make_tuple(num_positive_eigenvalues, num_negative_eigenvalues, num_zero_eigenvalues); + } + + size_t MA27Solver::number_negative_eigenvalues() const + { + return static_cast(info[INFO::NEIG]); + } + + bool MA27Solver::matrix_is_singular() const + { + return (info[INFO::IFLAG] == IFLAG::SINGULAR); + } + + size_t MA27Solver::rank() const + { + return info[INFO::IFLAG] == IFLAG::RANKDEFECT ? static_cast(info[INFO::IERROR]) : static_cast(n); + } + + void MA27Solver::save_matrix_to_local_format(const SymmetricMatrix &matrix) + { + // build the internal matrix representation + irn.clear(); + icn.clear(); + + for (const auto [row_index, column_index, element] : matrix) + { + irn.emplace_back(static_cast(row_index + fortran_shift)); + icn.emplace_back(static_cast(column_index + fortran_shift)); + } + } + +} // namespace uno \ No newline at end of file diff --git a/uno/solvers/MA27/MA27Solver.hpp b/uno/solvers/MA27/MA27Solver.hpp new file mode 100644 index 00000000..889ef650 --- /dev/null +++ b/uno/solvers/MA27/MA27Solver.hpp @@ -0,0 +1,60 @@ +#ifndef UNO_MA27SOLVER_H +#define UNO_MA27SOLVER_H + +#include +#include +#include "solvers/DirectSymmetricIndefiniteLinearSolver.hpp" + +namespace uno { +// forward declaration +template +class Vector; + + +class MA27Solver + : public DirectSymmetricIndefiniteLinearSolver +{ +public: + explicit MA27Solver(size_t max_dimension, size_t max_number_nonzeros); + ~MA27Solver() override = default; + + void factorize(const SymmetricMatrix& matrix) override; + void do_symbolic_factorization(const SymmetricMatrix& matrix) override; + void do_numerical_factorization(const SymmetricMatrix& matrix) override; + void solve_indefinite_system(const SymmetricMatrix& matrix, const Vector& rhs, Vector& result) override; + + + [[nodiscard]] std::tuple get_inertia() const override; + [[nodiscard]] size_t number_negative_eigenvalues() const override; + // [[nodiscard]] bool matrix_is_positive_definite() const override; + [[nodiscard]] bool matrix_is_singular() const override; + [[nodiscard]] size_t rank() const override; + +private: + 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 icntl{}; // integer array of length 30; integer control values + std::array cntl{}; // double array of length 5; double control values + + std::vector irn{}; // row index of input + std::vector icn{}; // col index of input + + std::vector iw{}; // integer workspace of length liw + std::vector ikeep{}; // integer array of 3*n; pivot sequence + std::vector 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 info{}; // integer array of length 20 + double ops{}; // double, operations count + + std::vector factor{}; // data array of length la; + int maxfrt{}; // integer, to be set by ma27 + + + // bool use_iterative_refinement{false}; // Not sure how to do this with ma27 + void save_matrix_to_local_format(const SymmetricMatrix& matrix); +}; + +} // namespace uno +#endif // UNO_MA27SOLVER_H \ No newline at end of file diff --git a/uno/solvers/SymmetricIndefiniteLinearSolverFactory.cpp b/uno/solvers/SymmetricIndefiniteLinearSolverFactory.cpp index edf79213..2059d201 100644 --- a/uno/solvers/SymmetricIndefiniteLinearSolverFactory.cpp +++ b/uno/solvers/SymmetricIndefiniteLinearSolverFactory.cpp @@ -12,6 +12,10 @@ #include "solvers/MA57/MA57Solver.hpp" #endif +#if defined(HAS_HSL) || defined(HAS_MA27) +#include "solvers/MA27/MA27Solver.hpp" +#endif + #ifdef HAS_HSL namespace uno { extern "C" { @@ -38,6 +42,17 @@ namespace uno { return std::make_unique(dimension, number_nonzeros); } #endif + +#if defined(HAS_HSL) || defined(HAS_MA27) + if (linear_solver_name == "MA27" +# ifdef HAS_HSL + && LIBHSL_isfunctional() +# endif + ) { + return std::make_unique(dimension, number_nonzeros); + } +#endif // HAS_HSL || HAS_MA27 + #ifdef HAS_MUMPS if (linear_solver_name == "MUMPS") { return std::make_unique(dimension, number_nonzeros); @@ -62,10 +77,15 @@ namespace uno { #ifdef HAS_HSL if (LIBHSL_isfunctional()) { solvers.emplace_back("MA57"); + solvers.emplace_back("MA27"); } #elif defined(HAS_MA57) solvers.emplace_back("MA57"); +#elif defined(HAS_MA27) + solvers.emplace_back("MA27"); #endif + + #ifdef HAS_MUMPS solvers.emplace_back("MUMPS"); #endif From 88552743987f7e1e6938e604d38c5c2b8b610fd9 Mon Sep 17 00:00:00 2001 From: Manuel Schaich Date: Sat, 30 Nov 2024 14:15:46 +0000 Subject: [PATCH 02/12] Allow refactorisation when workspace is insufficient. --- uno/solvers/MA27/MA27Solver.cpp | 32 ++++++++++++++++++++++++++++++++ uno/solvers/MA27/MA27Solver.hpp | 1 + 2 files changed, 33 insertions(+) diff --git a/uno/solvers/MA27/MA27Solver.cpp b/uno/solvers/MA27/MA27Solver.cpp index b099bac0..996746ac 100644 --- a/uno/solvers/MA27/MA27Solver.cpp +++ b/uno/solvers/MA27/MA27Solver.cpp @@ -170,6 +170,33 @@ enum IFLAG { } } + void MA27Solver::repeat_factorization_after_resizing() { + if (IFLAG::INSUFFICIENTINTEGER == info[INFO::IFLAG]) + { + INFO << "MA27: insufficient integer workspace, resizing and retrying. \n"; + // increase the size of iw + iw.resize(static_cast(info[INFO::IERROR])); + } + + if (IFLAG::INSUFFICIENTREAL == info[INFO::IFLAG]) + { + INFO << "MA27: insufficient real workspace, resizing and retrying. \n"; + // increase the size of factor + factor.resize(static_cast(info[INFO::IERROR])); + } + + int la = static_cast(factor.size()); + int liw = static_cast(iw.size()); + + MA27BD(&n, &nnz, irn.data(), icn.data(), factor.data(), &la, iw.data(), &liw, + ikeep.data(), &nsteps, &maxfrt, iw1.data(), icntl.data(), cntl.data(), info.data()); + + if (IFLAG::INSUFFICIENTINTEGER == info[INFO::IFLAG] || IFLAG::INSUFFICIENTREAL == info[INFO::IFLAG]) + { + repeat_factorization_after_resizing(); + } + } + void MA27Solver::do_numerical_factorization([[maybe_unused]]const SymmetricMatrix &matrix) { assert(matrix.dimension() <= iw1.capacity() && "MA27Solver: the dimension of the matrix is larger than the preallocated size"); @@ -181,6 +208,11 @@ enum IFLAG { MA27BD(&n, &nnz, irn.data(), icn.data(), factor.data(), &la, iw.data(), &liw, ikeep.data(), &nsteps, &maxfrt, iw1.data(), icntl.data(), cntl.data(), info.data()); + if (IFLAG::INSUFFICIENTINTEGER == info[INFO::IFLAG] || IFLAG::INSUFFICIENTREAL == info[INFO::IFLAG]) + { + repeat_factorization_after_resizing(); + } + assert(IFLAG::SUCCESS == info[INFO::IFLAG] && "MA27: the numerical factorization failed"); if (IFLAG::SUCCESS != info[INFO::IFLAG] ) { diff --git a/uno/solvers/MA27/MA27Solver.hpp b/uno/solvers/MA27/MA27Solver.hpp index 889ef650..7db2bd30 100644 --- a/uno/solvers/MA27/MA27Solver.hpp +++ b/uno/solvers/MA27/MA27Solver.hpp @@ -54,6 +54,7 @@ class MA27Solver // bool use_iterative_refinement{false}; // Not sure how to do this with ma27 void save_matrix_to_local_format(const SymmetricMatrix& matrix); + void repeat_factorization_after_resizing(); }; } // namespace uno From 3b6e5a4755ac5e19039e1fcca64e1d1fc10862e5 Mon Sep 17 00:00:00 2001 From: Manuel Schaich Date: Sat, 30 Nov 2024 14:48:24 +0000 Subject: [PATCH 03/12] INFO is a stream. Avoid aliasing. --- uno/solvers/MA27/MA27Solver.cpp | 81 +++++++++++++++++++++------------ 1 file changed, 52 insertions(+), 29 deletions(-) diff --git a/uno/solvers/MA27/MA27Solver.cpp b/uno/solvers/MA27/MA27Solver.cpp index 996746ac..f53b3217 100644 --- a/uno/solvers/MA27/MA27Solver.cpp +++ b/uno/solvers/MA27/MA27Solver.cpp @@ -27,9 +27,9 @@ extern "C" namespace uno { -constexpr auto fortran_shift = 1; -enum ICNTL { + +enum eICNTL { LP = 0, // sed by the subroutines as the output stream for error messages. If it is set to zero these messages will be suppressed. The default value is 6. MP, //used by the subroutines as the output stream for diagnostic printing and for warning messages. If it is set to zero then messages are suppressed. The default value is 6. LDIAG, // used by the subroutines to control diagnostic printing. If ICNTL(3) is equal to zero (the default), no diagnostic printing will be produced, a value of 1 will print scalar parameters (both in argument lists and in the control and information arrays) and a few entries of array parameters on entry and successful exit from each subroutine while ICNTL(3) equal to 2 will print all parameter values on entry and successful exit. @@ -64,7 +64,7 @@ enum ICNTL { UNUSED_ICNTL5, }; -enum CNTL { +enum eCNTL { U = 0, // used by the subroutine to control numerical pivoting. Values greater than 0.5 are treated as 0.5 and less than –0.5 as –0.5. Its default value is 0.1. If U is positive, numerical pivoting will be performed. If U is non-positive, no pivoting will be performed, the subroutine will fail if a zero pivot is encountered, and a flag (see section 2.3) will be set if not all pivots are of the same sign; the factorization will continue after a sign change is detected if U is zero but will exit immediately if U is less than zero. If the system is definite, then setting U to zero will decrease the factorization time while still providing a stable decomposition. For problems requiring greater than average numerical care a higher value than the default would be advisable. FRATIO, // given the default value of 1.0 by MA27I/ID. If MA27A/AD encounters a row of the reduced matrix with a proportion of entries greater than FRATIO, the row is treated as full. FRATIO is not altered by MA27. PIVTOL, // given the default value of 0.0 by MA27I/ID. MA27B/BD will not accept an entry with absolute value less than PIVTOL as a 1×1 pivot or the off-diagonal entry of a 2×2 pivot. PIVTOL is not altered by MA27. @@ -72,7 +72,7 @@ enum CNTL { UNUSED_CNTL2, }; -enum INFO { +enum eINFO { IFLAG = 0, // is an error flag. A value of zero indicates that the subroutine has performed successfully. IERROR, // provides supplementary information when there is an error. NRLTOT, // gives the total amount of REAL words required for a successful completion of MA27B/BD without the need for data compression provided no numerical pivoting is performed. The actual amount required may be higher because of numerical pivoting, but probably not by more than 3%. @@ -95,7 +95,7 @@ enum INFO { UNUSED_INFO5, }; -enum IFLAG { +enum eIFLAG { NSTEPS = -7, // Value of NSTEPS outside the range 1 ≤ NSTEPS ≤ N (MA27B/BD entry). PIVOTSIGN = -6, // A change of sign of pivots has been detected when U was negative. INFO(2) is set to the pivot step at which the change was detected. (MA27B/BD entry only) SINGULAR = -5, // Matrix is singular (MA27B/BD entry only). INFO(2) is set to the pivot step at which singularity was detected @@ -126,9 +126,9 @@ enum IFLAG { // set the default values of the controlling parameters MA27ID(icntl.data(), cntl.data()); // suppress warning messages - icntl[ICNTL::LP] = 0; - icntl[ICNTL::MP] = 0; - icntl[ICNTL::LDIAG] = 0; + icntl[eICNTL::LP] = 0; + icntl[eICNTL::MP] = 0; + icntl[eICNTL::LDIAG] = 0; } void MA27Solver::factorize(const SymmetricMatrix &matrix) @@ -159,30 +159,30 @@ enum IFLAG { icntl.data(), cntl.data(), info.data(), &ops); - factor.resize(static_cast(3 * info[INFO::NRLNEC] / 2)); + factor.resize(static_cast(3 * info[eINFO::NRLNEC] / 2)); std::copy(matrix.data_pointer(), matrix.data_pointer() + matrix.number_nonzeros(), factor.begin()); - assert(IFLAG::SUCCESS == info[INFO::IFLAG] && "MA27: the symbolic factorization failed"); - if (IFLAG::SUCCESS != info[INFO::IFLAG] ) + 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[INFO::IFLAG] << " additional info, IERROR = " << info[INFO::IERROR] << '\n'; + WARNING << "MA27 has issued a warning: IFLAG = " << info[eINFO::IFLAG] << " additional info, IERROR = " << info[eINFO::IERROR] << '\n'; } } void MA27Solver::repeat_factorization_after_resizing() { - if (IFLAG::INSUFFICIENTINTEGER == info[INFO::IFLAG]) + if (eIFLAG::INSUFFICIENTINTEGER == info[eINFO::IFLAG]) { INFO << "MA27: insufficient integer workspace, resizing and retrying. \n"; // increase the size of iw - iw.resize(static_cast(info[INFO::IERROR])); + iw.resize(static_cast(info[eINFO::IERROR])); } - if (IFLAG::INSUFFICIENTREAL == info[INFO::IFLAG]) + if (eIFLAG::INSUFFICIENTREAL == info[eINFO::IFLAG]) { INFO << "MA27: insufficient real workspace, resizing and retrying. \n"; // increase the size of factor - factor.resize(static_cast(info[INFO::IERROR])); + factor.resize(static_cast(info[eINFO::IERROR])); } int la = static_cast(factor.size()); @@ -191,7 +191,7 @@ enum IFLAG { MA27BD(&n, &nnz, irn.data(), icn.data(), factor.data(), &la, iw.data(), &liw, ikeep.data(), &nsteps, &maxfrt, iw1.data(), icntl.data(), cntl.data(), info.data()); - if (IFLAG::INSUFFICIENTINTEGER == info[INFO::IFLAG] || IFLAG::INSUFFICIENTREAL == info[INFO::IFLAG]) + if (eIFLAG::INSUFFICIENTINTEGER == info[eINFO::IFLAG] || eIFLAG::INSUFFICIENTREAL == info[eINFO::IFLAG]) { repeat_factorization_after_resizing(); } @@ -208,16 +208,39 @@ enum IFLAG { MA27BD(&n, &nnz, irn.data(), icn.data(), factor.data(), &la, iw.data(), &liw, ikeep.data(), &nsteps, &maxfrt, iw1.data(), icntl.data(), cntl.data(), info.data()); - if (IFLAG::INSUFFICIENTINTEGER == info[INFO::IFLAG] || IFLAG::INSUFFICIENTREAL == info[INFO::IFLAG]) + if (eIFLAG::INSUFFICIENTINTEGER == info[eINFO::IFLAG] || eIFLAG::INSUFFICIENTREAL == info[eINFO::IFLAG]) { repeat_factorization_after_resizing(); } - assert(IFLAG::SUCCESS == info[INFO::IFLAG] && "MA27: the numerical factorization failed"); - if (IFLAG::SUCCESS != info[INFO::IFLAG] ) + switch (info[eINFO::IFLAG]) { - WARNING << "MA27 has issued a warning: IFLAG = " << info[INFO::IFLAG] << " additional info, IERROR = " << info[INFO::IERROR] << '\n'; - } + case NSTEPS: + WARNING << "MA27BD: Value of NSTEPS outside the range 1 ≤ NSTEPS ≤ N" << '\n'; + break; + case PIVOTSIGN: + WARNING << "MA27BD: A change of sign of pivots has been detected when U was negative. Detected at pivot step " << info[eINFO::IERROR] << '\n'; + break; + case SINGULAR: + WARNING << "MA27BD: Matrix is singular. Singularity detected during pivot step "<< info[eINFO::IERROR] << '\n'; + break; + case NZOUTOFRANGE: + WARNING << "MA27BD: Value of NZ out of range. NZ < 0." << '\n'; + break; + case NOUTOFRANGE: + WARNING << "MA27BD: Value of N out of range. N < 1." << '\n'; + break; + case IDXOUTOFRANGE: + WARNING << "MA27BD: Index (in IRN or ICN) out of range. " << info[eINFO::IERROR] << " indices affected." << '\n'; + break; + case FALSEDEFINITENESS: + WARNING << "MA27BD: Matrix was supposed to be definite, but pivots have different signs when factorizing. Detected " << info[eINFO::IERROR] << " sign changes." << '\n'; + break; + case RANKDEFECT: + WARNING << "MA27BD: Matrix is rank deficient. Rank: " << info[eINFO::IERROR] << " whereas dimension " << n << '\n'; + break; + } + } void MA27Solver::solve_indefinite_system([[maybe_unused]]const SymmetricMatrix &matrix, const Vector &rhs, Vector &result) @@ -242,10 +265,10 @@ enum IFLAG { icntl.data(), info.data()); - assert(info[INFO::IFLAG] == IFLAG::SUCCESS && "MA27: the solution failed"); - if (IFLAG::SUCCESS != info[INFO::IFLAG] ) + assert(info[eINFO::IFLAG] == eIFLAG::SUCCESS && "MA27: the solution failed"); + if (eIFLAG::SUCCESS != info[eINFO::IFLAG] ) { - WARNING << "MA27 has issued a warning: IFLAG = " << info[INFO::IFLAG] << " additional info, IERROR = " << info[INFO::IERROR] << '\n'; + WARNING << "MA27 has issued a warning: IFLAG = " << info[eINFO::IFLAG] << " additional info, IERROR = " << info[eINFO::IERROR] << '\n'; } } @@ -262,17 +285,17 @@ enum IFLAG { size_t MA27Solver::number_negative_eigenvalues() const { - return static_cast(info[INFO::NEIG]); + return static_cast(info[eINFO::NEIG]); } bool MA27Solver::matrix_is_singular() const { - return (info[INFO::IFLAG] == IFLAG::SINGULAR); + return (info[eINFO::IFLAG] == eIFLAG::SINGULAR); } size_t MA27Solver::rank() const { - return info[INFO::IFLAG] == IFLAG::RANKDEFECT ? static_cast(info[INFO::IERROR]) : static_cast(n); + return info[eINFO::IFLAG] == eIFLAG::RANKDEFECT ? static_cast(info[eINFO::IERROR]) : static_cast(n); } void MA27Solver::save_matrix_to_local_format(const SymmetricMatrix &matrix) @@ -280,7 +303,7 @@ enum IFLAG { // build the internal matrix representation irn.clear(); icn.clear(); - + constexpr auto fortran_shift = 1; for (const auto [row_index, column_index, element] : matrix) { irn.emplace_back(static_cast(row_index + fortran_shift)); From 30588e2e5930ca0f4d7531c43b3d0631c0f7e6e9 Mon Sep 17 00:00:00 2001 From: Manuel Schaich Date: Sat, 30 Nov 2024 15:05:14 +0000 Subject: [PATCH 04/12] singular matrix detection. --- uno/solvers/MA27/MA27Solver.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/uno/solvers/MA27/MA27Solver.cpp b/uno/solvers/MA27/MA27Solver.cpp index f53b3217..f626ba94 100644 --- a/uno/solvers/MA27/MA27Solver.cpp +++ b/uno/solvers/MA27/MA27Solver.cpp @@ -290,7 +290,7 @@ enum eIFLAG { bool MA27Solver::matrix_is_singular() const { - return (info[eINFO::IFLAG] == eIFLAG::SINGULAR); + return (info[eINFO::IFLAG] == eIFLAG::SINGULAR || info[eINFO::IFLAG] == eIFLAG::RANKDEFECT); } size_t MA27Solver::rank() const From 9a0bf9310f5e1f975ace43a4a8256f83dc211151 Mon Sep 17 00:00:00 2001 From: Manuel Schaich Date: Sat, 30 Nov 2024 15:05:32 +0000 Subject: [PATCH 05/12] Unit tests for ma27 --- unotest/functional_tests/MA27SolverTests.cpp | 77 ++++++++++++++++++++ 1 file changed, 77 insertions(+) create mode 100644 unotest/functional_tests/MA27SolverTests.cpp diff --git a/unotest/functional_tests/MA27SolverTests.cpp b/unotest/functional_tests/MA27SolverTests.cpp new file mode 100644 index 00000000..f8e80da9 --- /dev/null +++ b/unotest/functional_tests/MA27SolverTests.cpp @@ -0,0 +1,77 @@ +// Copyright (c) 2018-2024 Charlie Vanaret +// Licensed under the MIT license. See LICENSE file in the project directory for details. + +#include +#include +#include "linear_algebra/SymmetricMatrix.hpp" +#include "solvers/MA27/MA27Solver.hpp" + +using namespace uno; + +TEST(MA27Solver, SystemSize5) { + const size_t n = 5; + const size_t nnz = 7; + SymmetricMatrix matrix(n, nnz, false, "COO"); + matrix.insert(2., 0, 0); + matrix.insert(3., 0, 1); + matrix.insert(4., 1, 2); + matrix.insert(6., 1, 4); + matrix.insert(1., 2, 2); + matrix.insert(5., 2, 3); + matrix.insert(1., 4, 4); + const Vector rhs{8., 45., 31., 15., 17.}; + Vector result(n); + result.fill(0.); + const std::array reference{1., 2., 3., 4., 5.}; + + MA27Solver solver(n, nnz); + solver.do_symbolic_factorization(matrix); + solver.do_numerical_factorization(matrix); + solver.solve_indefinite_system(matrix, rhs, result); + + for (size_t index: Range(n)) { + EXPECT_DOUBLE_EQ(result[index], reference[index]); + } +} + +TEST(MA27Solver, Inertia) { + const size_t n = 5; + const size_t nnz = 7; + SymmetricMatrix matrix(n, nnz, false, "COO"); + matrix.insert(2., 0, 0); + matrix.insert(3., 0, 1); + matrix.insert(4., 1, 2); + matrix.insert(6., 1, 4); + matrix.insert(1., 2, 2); + matrix.insert(5., 2, 3); + matrix.insert(1., 4, 4); + + MA27Solver solver(n, nnz); + solver.do_symbolic_factorization(matrix); + solver.do_numerical_factorization(matrix); + + const auto [number_positive, number_negative, number_zero] = solver.get_inertia(); + ASSERT_EQ(number_positive, 3); + ASSERT_EQ(number_negative, 2); + ASSERT_EQ(number_zero, 0); +} + +TEST(MA27Solver, SingularMatrix) { + const size_t n = 4; + const size_t nnz = 7; + // comes from hs015 solved with byrd preset + SymmetricMatrix matrix(n, nnz, false, "COO"); + matrix.insert( -0.0198, 0, 0); + matrix.insert(0.625075, 0, 0); + matrix.insert(-0.277512, 0, 1); + matrix.insert(-0.624975, 1, 1); + matrix.insert(0.625075, 1, 1); + matrix.insert(0., 2, 2); + matrix.insert(0., 3, 3); + MA27Solver solver(n, nnz); + solver.do_symbolic_factorization(matrix); + solver.do_numerical_factorization(matrix); + + // expected inertia (1, 1, 2) + ASSERT_TRUE(solver.matrix_is_singular()); +} From ea258b0f5a0972aeb8f493775ed678daf4fa8fc2 Mon Sep 17 00:00:00 2001 From: Manuel Schaich Date: Sat, 30 Nov 2024 15:10:33 +0000 Subject: [PATCH 06/12] Best guess for ma27 build --- CMakeLists.txt | 25 +++++++++++++++++++++++-- 1 file changed, 23 insertions(+), 2 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 4fbaa1f1..973f6f01 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,4 +1,4 @@ -# Copyright (c) 2018-2024 Charlie Vanaret +# Copyright (c) 2018-2024 Charlie Vanaret # Licensed under the MIT license. See LICENSE file in the project directory for details. cmake_minimum_required(VERSION 3.7) @@ -107,11 +107,32 @@ else() else() message(WARNING "Optional library MA57 was not found.") endif() + find_library(MA27 ma27) + if(MA27) + link_to_uno(ma27 ${MA27}) + else() + message(WARNING "Optional library MA27 was not found.") + endif() endif() -if(HSL OR MA57) +if(HSL) + list(APPEND UNO_SOURCE_FILES uno/solvers/MA57/MA57Solver.cpp) + list(APPEND UNO_SOURCE_FILES uno/solvers/MA27/MA27Solver.cpp) + + list(APPEND TESTS_UNO_SOURCE_FILES unotest/functional_tests/MA57SolverTests.cpp) + list(APPEND TESTS_UNO_SOURCE_FILES unotest/functional_tests/MA27SolverTests.cpp) + + find_package(BLAS REQUIRED) + list(APPEND LIBRARIES ${BLAS_LIBRARIES}) +elseif(MA57) list(APPEND UNO_SOURCE_FILES uno/solvers/MA57/MA57Solver.cpp) list(APPEND TESTS_UNO_SOURCE_FILES unotest/functional_tests/MA57SolverTests.cpp) + find_package(BLAS REQUIRED) + list(APPEND LIBRARIES ${BLAS_LIBRARIES}) +elseif(MA27) + list(APPEND UNO_SOURCE_FILES uno/solvers/MA27/MA27Solver.cpp) + list(APPEND TESTS_UNO_SOURCE_FILES unotest/functional_tests/MA27SolverTests.cpp) + find_package(BLAS REQUIRED) list(APPEND LIBRARIES ${BLAS_LIBRARIES}) endif() From e048dc2fb12dfd35f6404b32cac03a51ade2e0de Mon Sep 17 00:00:00 2001 From: Manuel Schaich Date: Sat, 30 Nov 2024 15:13:38 +0000 Subject: [PATCH 07/12] Typos etc. --- uno/solvers/MA27/MA27Solver.cpp | 42 ++++++++++++++++----------------- 1 file changed, 21 insertions(+), 21 deletions(-) diff --git a/uno/solvers/MA27/MA27Solver.cpp b/uno/solvers/MA27/MA27Solver.cpp index f626ba94..791100d6 100644 --- a/uno/solvers/MA27/MA27Solver.cpp +++ b/uno/solvers/MA27/MA27Solver.cpp @@ -30,9 +30,9 @@ namespace uno enum eICNTL { - LP = 0, // sed by the subroutines as the output stream for error messages. If it is set to zero these messages will be suppressed. The default value is 6. - MP, //used by the subroutines as the output stream for diagnostic printing and for warning messages. If it is set to zero then messages are suppressed. The default value is 6. - LDIAG, // used by the subroutines to control diagnostic printing. If ICNTL(3) is equal to zero (the default), no diagnostic printing will be produced, a value of 1 will print scalar parameters (both in argument lists and in the control and information arrays) and a few entries of array parameters on entry and successful exit from each subroutine while ICNTL(3) equal to 2 will print all parameter values on entry and successful exit. + LP = 0, // Used by the subroutines as the output stream for error messages. If it is set to zero these messages will be suppressed. The default value is 6. + MP, // Used by the subroutines as the output stream for diagnostic printing and for warning messages. If it is set to zero then messages are suppressed. The default value is 6. + LDIAG, // Used by the subroutines to control diagnostic printing. If ICNTL(3) is equal to zero (the default), no diagnostic printing will be produced, a value of 1 will print scalar parameters (both in argument lists and in the control and information arrays) and a few entries of array parameters on entry and successful exit from each subroutine while ICNTL(3) equal to 2 will print all parameter values on entry and successful exit. /* The entries ICNTL(4) to ICNTL(25) are not of interest to the general user and are discussed more fully by Duff and Reid (AERE R-10533, 1982) under the internal names IOVFLO, NEMIN and IFRLVL */ IOVFLO, @@ -65,29 +65,29 @@ enum eICNTL { }; enum eCNTL { - U = 0, // used by the subroutine to control numerical pivoting. Values greater than 0.5 are treated as 0.5 and less than –0.5 as –0.5. Its default value is 0.1. If U is positive, numerical pivoting will be performed. If U is non-positive, no pivoting will be performed, the subroutine will fail if a zero pivot is encountered, and a flag (see section 2.3) will be set if not all pivots are of the same sign; the factorization will continue after a sign change is detected if U is zero but will exit immediately if U is less than zero. If the system is definite, then setting U to zero will decrease the factorization time while still providing a stable decomposition. For problems requiring greater than average numerical care a higher value than the default would be advisable. - FRATIO, // given the default value of 1.0 by MA27I/ID. If MA27A/AD encounters a row of the reduced matrix with a proportion of entries greater than FRATIO, the row is treated as full. FRATIO is not altered by MA27. - PIVTOL, // given the default value of 0.0 by MA27I/ID. MA27B/BD will not accept an entry with absolute value less than PIVTOL as a 1×1 pivot or the off-diagonal entry of a 2×2 pivot. PIVTOL is not altered by MA27. + U = 0, // Used by the subroutine to control numerical pivoting. Values greater than 0.5 are treated as 0.5 and less than –0.5 as –0.5. Its default value is 0.1. If U is positive, numerical pivoting will be performed. If U is non-positive, no pivoting will be performed, the subroutine will fail if a zero pivot is encountered, and a flag (see section 2.3) will be set if not all pivots are of the same sign; the factorization will continue after a sign change is detected if U is zero but will exit immediately if U is less than zero. If the system is definite, then setting U to zero will decrease the factorization time while still providing a stable decomposition. For problems requiring greater than average numerical care a higher value than the default would be advisable. + FRATIO, // Given the default value of 1.0 by MA27I/ID. If MA27A/AD encounters a row of the reduced matrix with a proportion of entries greater than FRATIO, the row is treated as full. FRATIO is not altered by MA27. + PIVTOL, // Given the default value of 0.0 by MA27I/ID. MA27B/BD will not accept an entry with absolute value less than PIVTOL as a 1×1 pivot or the off-diagonal entry of a 2×2 pivot. PIVTOL is not altered by MA27. UNUSED_CNTL1, UNUSED_CNTL2, }; enum eINFO { - IFLAG = 0, // is an error flag. A value of zero indicates that the subroutine has performed successfully. - IERROR, // provides supplementary information when there is an error. - NRLTOT, // gives the total amount of REAL words required for a successful completion of MA27B/BD without the need for data compression provided no numerical pivoting is performed. The actual amount required may be higher because of numerical pivoting, but probably not by more than 3%. - NIRTOT, // gives the total amount of INTEGER words required for a successful completion of MA27B/BD without the need for data compression provided no numerical pivoting is performed. The actual amount required may be higher because of numerical pivoting, but probably not by more than 3%. - NRLNEC, // gives the amount of REAL words required for successful completion of MA27B/BD allowing data compression (see NCMPBR returned in INFO(12)), again provided no numerical pivoting is performed. Numerical pivoting may cause a higher value to be required, but probably not by more than 3%. If storage was conserved by equivalencing IW(1) with IRN(1), NRLNEC and NIRNEC cannot be calculated exactly but instead an upper bound will be returned. Experience has shown that this can overestimate the exact values by 50% although the tightness of the bound is very problem dependent. For example, a tight bound will generally be obtained if there are many more entries in the factors than in the input matrix. - NIRNEC, // gives the amount of INTEGER words required for successful completion of MA27B/BD allowing data compression (see NCMPBR returned in INFO(12)), again provided no numerical pivoting is performed. Numerical pivoting may cause a higher value to be required, but probably not by more than 3%. If storage was conserved by equivalencing IW(1) with IRN(1), NRLNEC and NIRNEC cannot be calculated exactly but instead an upper bound will be returned. Experience has shown that this can overestimate the exact values by 50% although the tightness of the bound is very problem dependent. For example, a tight bound will generally be obtained if there are many more entries in the factors than in the input matrix. - NRLADU, // gives the number of REAL words required to hold the matrix factors if no numerical pivoting is performed by MA27B/BD. Numerical pivoting may change this slightly. - NIRADU, // gives the number of INTEGER words required to hold the matrix factors if no numerical pivoting is performed by MA27B/BD. Numerical pivoting may change this slightly. - NRLBDU, // give the amount of REAL words actually used to hold the factorization. - NIRBDU, // give the amount of INTEGER words actually used to hold the factorization. - NCMPA, // holds the number of compresses of the internal data structure performed by MA27A/AD. If this is high (say > 10), the performance of MA27A/AD may be improved by increasing the length of array IW. - NCMPBR, // holds the number of compresses of the real data structure required by the factorization. If either of these is high (say > 10), then the speed of the factorization may be increased by allocating more space to the arrays A as appropriate. - NCMPBI, // holds the number of compresses of the integer data structure required by the factorization. If either of these is high (say > 10), then the speed of the factorization may be increased by allocating more space to the arrays IW as appropriate. - NTWO, // gives the number of 2×2 pivots used during the factorization. - NEIG, // gives the number of negative eigenvalues of A. + IFLAG = 0, // An error flag. A value of zero indicates that the subroutine has performed successfully. + IERROR, // Provides supplementary information when there is an error. + NRLTOT, // Gives the total amount of REAL words required for a successful completion of MA27B/BD without the need for data compression provided no numerical pivoting is performed. The actual amount required may be higher because of numerical pivoting, but probably not by more than 3%. + NIRTOT, // Gives the total amount of INTEGER words required for a successful completion of MA27B/BD without the need for data compression provided no numerical pivoting is performed. The actual amount required may be higher because of numerical pivoting, but probably not by more than 3%. + NRLNEC, // Gives the amount of REAL words required for successful completion of MA27B/BD allowing data compression (see NCMPBR returned in INFO(12)), again provided no numerical pivoting is performed. Numerical pivoting may cause a higher value to be required, but probably not by more than 3%. If storage was conserved by equivalencing IW(1) with IRN(1), NRLNEC and NIRNEC cannot be calculated exactly but instead an upper bound will be returned. Experience has shown that this can overestimate the exact values by 50% although the tightness of the bound is very problem dependent. For example, a tight bound will generally be obtained if there are many more entries in the factors than in the input matrix. + NIRNEC, // Gives the amount of INTEGER words required for successful completion of MA27B/BD allowing data compression (see NCMPBR returned in INFO(12)), again provided no numerical pivoting is performed. Numerical pivoting may cause a higher value to be required, but probably not by more than 3%. If storage was conserved by equivalencing IW(1) with IRN(1), NRLNEC and NIRNEC cannot be calculated exactly but instead an upper bound will be returned. Experience has shown that this can overestimate the exact values by 50% although the tightness of the bound is very problem dependent. For example, a tight bound will generally be obtained if there are many more entries in the factors than in the input matrix. + NRLADU, // Gives the number of REAL words required to hold the matrix factors if no numerical pivoting is performed by MA27B/BD. Numerical pivoting may change this slightly. + NIRADU, // Gives the number of INTEGER words required to hold the matrix factors if no numerical pivoting is performed by MA27B/BD. Numerical pivoting may change this slightly. + NRLBDU, // Gives the amount of REAL words actually used to hold the factorization. + NIRBDU, // Gives the amount of INTEGER words actually used to hold the factorization. + NCMPA, // Holds the number of compresses of the internal data structure performed by MA27A/AD. If this is high (say > 10), the performance of MA27A/AD may be improved by increasing the length of array IW. + NCMPBR, // Holds the number of compresses of the real data structure required by the factorization. If either of these is high (say > 10), then the speed of the factorization may be increased by allocating more space to the arrays A as appropriate. + NCMPBI, // Holds the number of compresses of the integer data structure required by the factorization. If either of these is high (say > 10), then the speed of the factorization may be increased by allocating more space to the arrays IW as appropriate. + NTWO, // Gives the number of 2×2 pivots used during the factorization. + NEIG, // Gives the number of negative eigenvalues of A. UNUSED_INFO1, UNUSED_INFO2, UNUSED_INFO3, From c6e48fe48670629309c1baac6f87773615bdae3a Mon Sep 17 00:00:00 2001 From: Manuel Schaich Date: Sat, 30 Nov 2024 15:33:06 +0000 Subject: [PATCH 08/12] Missing link between factors and matrix passed. --- uno/solvers/MA27/MA27Solver.cpp | 8 +++++--- uno/solvers/MA27/MA27Solver.hpp | 2 +- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/uno/solvers/MA27/MA27Solver.cpp b/uno/solvers/MA27/MA27Solver.cpp index 791100d6..2111ffed 100644 --- a/uno/solvers/MA27/MA27Solver.cpp +++ b/uno/solvers/MA27/MA27Solver.cpp @@ -170,7 +170,7 @@ enum eIFLAG { } } - void MA27Solver::repeat_factorization_after_resizing() { + void MA27Solver::repeat_factorization_after_resizing([[maybe_unused]]const SymmetricMatrix &matrix) { if (eIFLAG::INSUFFICIENTINTEGER == info[eINFO::IFLAG]) { INFO << "MA27: insufficient integer workspace, resizing and retrying. \n"; @@ -193,7 +193,7 @@ enum eIFLAG { if (eIFLAG::INSUFFICIENTINTEGER == info[eINFO::IFLAG] || eIFLAG::INSUFFICIENTREAL == info[eINFO::IFLAG]) { - repeat_factorization_after_resizing(); + repeat_factorization_after_resizing(matrix); } } @@ -210,7 +210,7 @@ enum eIFLAG { if (eIFLAG::INSUFFICIENTINTEGER == info[eINFO::IFLAG] || eIFLAG::INSUFFICIENTREAL == info[eINFO::IFLAG]) { - repeat_factorization_after_resizing(); + repeat_factorization_after_resizing(matrix); } switch (info[eINFO::IFLAG]) @@ -303,11 +303,13 @@ enum eIFLAG { // build the internal matrix representation irn.clear(); icn.clear(); + factor.clear(); constexpr auto fortran_shift = 1; for (const auto [row_index, column_index, element] : matrix) { irn.emplace_back(static_cast(row_index + fortran_shift)); icn.emplace_back(static_cast(column_index + fortran_shift)); + factor.emplace_back(element); } } diff --git a/uno/solvers/MA27/MA27Solver.hpp b/uno/solvers/MA27/MA27Solver.hpp index 7db2bd30..8c56b03e 100644 --- a/uno/solvers/MA27/MA27Solver.hpp +++ b/uno/solvers/MA27/MA27Solver.hpp @@ -54,7 +54,7 @@ class MA27Solver // bool use_iterative_refinement{false}; // Not sure how to do this with ma27 void save_matrix_to_local_format(const SymmetricMatrix& matrix); - void repeat_factorization_after_resizing(); + void repeat_factorization_after_resizing(const SymmetricMatrix &matrix); }; } // namespace uno From 454e2500fc18f3bcfd1764fa454cbbc164119bf7 Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Sat, 30 Nov 2024 19:59:12 +0100 Subject: [PATCH 09/12] Fixed CMakeLists.txt: we may want to link MA57 and MA27 simultaneously --- CMakeLists.txt | 14 +++----------- 1 file changed, 3 insertions(+), 11 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 973f6f01..423765fd 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -114,22 +114,14 @@ else() message(WARNING "Optional library MA27 was not found.") endif() endif() -if(HSL) - list(APPEND UNO_SOURCE_FILES uno/solvers/MA57/MA57Solver.cpp) - list(APPEND UNO_SOURCE_FILES uno/solvers/MA27/MA27Solver.cpp) - - list(APPEND TESTS_UNO_SOURCE_FILES unotest/functional_tests/MA57SolverTests.cpp) - list(APPEND TESTS_UNO_SOURCE_FILES unotest/functional_tests/MA27SolverTests.cpp) - - find_package(BLAS REQUIRED) - list(APPEND LIBRARIES ${BLAS_LIBRARIES}) -elseif(MA57) +if(HSL OR MA57) list(APPEND UNO_SOURCE_FILES uno/solvers/MA57/MA57Solver.cpp) list(APPEND TESTS_UNO_SOURCE_FILES unotest/functional_tests/MA57SolverTests.cpp) find_package(BLAS REQUIRED) list(APPEND LIBRARIES ${BLAS_LIBRARIES}) -elseif(MA27) +endif() +if(HSL OR MA27) list(APPEND UNO_SOURCE_FILES uno/solvers/MA27/MA27Solver.cpp) list(APPEND TESTS_UNO_SOURCE_FILES unotest/functional_tests/MA27SolverTests.cpp) From 1d36b4c093ce46238f45830f20e11c18555d5a2b Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Sat, 30 Nov 2024 20:07:42 +0100 Subject: [PATCH 10/12] Fixed SymmetricIndefiniteLinearSolverFactory.cpp: we may have MA57 and MA27 simultaneously --- .../SymmetricIndefiniteLinearSolverFactory.cpp | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/uno/solvers/SymmetricIndefiniteLinearSolverFactory.cpp b/uno/solvers/SymmetricIndefiniteLinearSolverFactory.cpp index 2059d201..1d83c324 100644 --- a/uno/solvers/SymmetricIndefiniteLinearSolverFactory.cpp +++ b/uno/solvers/SymmetricIndefiniteLinearSolverFactory.cpp @@ -35,9 +35,9 @@ namespace uno { [[maybe_unused]] const std::string& linear_solver_name = options.get_string("linear_solver"); #if defined(HAS_HSL) || defined(HAS_MA57) if (linear_solver_name == "MA57" -#ifdef HAS_HSL + #ifdef HAS_HSL && LIBHSL_isfunctional() -#endif + #endif ) { return std::make_unique(dimension, number_nonzeros); } @@ -45,9 +45,9 @@ namespace uno { #if defined(HAS_HSL) || defined(HAS_MA27) if (linear_solver_name == "MA27" -# ifdef HAS_HSL + # ifdef HAS_HSL && LIBHSL_isfunctional() -# endif + # endif ) { return std::make_unique(dimension, number_nonzeros); } @@ -79,16 +79,18 @@ namespace uno { solvers.emplace_back("MA57"); solvers.emplace_back("MA27"); } -#elif defined(HAS_MA57) +#else + #ifdef HAS_MA57 solvers.emplace_back("MA57"); -#elif defined(HAS_MA27) + #endif + #ifdef HAS_MA27 solvers.emplace_back("MA27"); + #endif #endif - #ifdef HAS_MUMPS solvers.emplace_back("MUMPS"); #endif return solvers; } -} // namespace \ No newline at end of file +} // namespace From c612b8cf658099f72fe802d1704613618ee37cc1 Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Sat, 30 Nov 2024 20:11:12 +0100 Subject: [PATCH 11/12] Added copyright in MA27Solver.hpp --- uno/solvers/MA27/MA27Solver.hpp | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/uno/solvers/MA27/MA27Solver.hpp b/uno/solvers/MA27/MA27Solver.hpp index 8c56b03e..149f68cb 100644 --- a/uno/solvers/MA27/MA27Solver.hpp +++ b/uno/solvers/MA27/MA27Solver.hpp @@ -1,3 +1,6 @@ +// Copyright (c) 2024 Manuel Schaich +// Licensed under the MIT license. See LICENSE file in the project directory for details. + #ifndef UNO_MA27SOLVER_H #define UNO_MA27SOLVER_H @@ -58,4 +61,4 @@ class MA27Solver }; } // namespace uno -#endif // UNO_MA27SOLVER_H \ No newline at end of file +#endif // UNO_MA27SOLVER_H From 27af512a897ec8878ae0e92d604b219e979f0fb7 Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Sat, 30 Nov 2024 20:11:45 +0100 Subject: [PATCH 12/12] Added copyright in MA27Solver.cpp --- uno/solvers/MA27/MA27Solver.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/uno/solvers/MA27/MA27Solver.cpp b/uno/solvers/MA27/MA27Solver.cpp index 2111ffed..61c34070 100644 --- a/uno/solvers/MA27/MA27Solver.cpp +++ b/uno/solvers/MA27/MA27Solver.cpp @@ -1,3 +1,5 @@ +// Copyright (c) 2024 Manuel Schaich +// Licensed under the MIT license. See LICENSE file in the project directory for details. #include #include "MA27Solver.hpp" @@ -313,4 +315,4 @@ enum eIFLAG { } } -} // namespace uno \ No newline at end of file +} // namespace uno