Linear solvers for Ipopt

[jbcaillau]

@PierreMartinon @gergaud @ocots Many thanks @amontoison for your input on this. To summarise:

• we should add to our direct solve an option allowing to change the linear solver used by Ipopt: the default is mumps, but there (at least) the other options below: ma27 and ma57 from HSL¹, and spral from STFC/RAL:

stats = ipopt(nlp, linear_solver="mumps")
stats = ipopt(nlp, linear_solver="spral")

using HSL_jll
stats = ipopt(nlp, linear_solver="ma27")
stats = ipopt(nlp, linear_solver="ma57")

• we should default to spral instead of mumps in our code, provided that Julia $\geq$ 1.9 is used and that the user defined environnement variables below are set:

export OMP_CANCELLATION=TRUE
export OMP_PROC_BIND=TRUE

Extra

Users on Intel processors (and Julia $\geq$ 1.9 that has libblastrampoline) can change their OpenBLAS backend to the Intel one:

using MKL  # Replace OpenBLAS by Intel MKL

To display the current backend, do

import LinearAlgebra
LinearAlgebra.BLAS.lbt_get_config()

See also

Performance tips

Footnotes

1. free academic license available on request ↩

[PierreMartinon]

Agreed. Would be interesting to check with some benchmarking.

[amontoison]

@PierreMartinon Do you use the function ipopt only in solve.jl?
I can open a PR to easily support HSL solvers.

[ocots]

@amontoison I think we use Ipopt only here

CTDirect.jl/src/solve.jl

Line 77 in b0b6921

function solve(docp::DOCP;

We can see this here: https://github.com/search?q=repo%3Acontrol-toolbox%2FCTDirect.jl%20ipopt&type=code

[jbcaillau]

@PierreMartinon Do you use the function ipopt only in solve.jl? I can open a PR to easily support HSL solvers.

@amontoison YES. PR most welcome on this. All the more as there is currently an internship with @ocots @frapac that will benefit a lot from being able to easily try and change linear solvers.

[amontoison]

Done in #113 @jbcaillau 😉

[jbcaillau]

@ocots @0Yassine0 @amontoison guys, I need a serious update, here:

• the initial post above goes back to Sep. 2023: are the information still up to date?
• after Support HSL solvers #113 is it really as simple as using HSL? what about importing HSL_jll? licensing...? not clear to me
• what about using spral instead of mumps? (I had noticed Ipopt was compiled with it for >= 1.9)

[amontoison]

@ocots @0Yassine0 @amontoison guys, I need a serious update, here:

• the initial post above goes back to Sep. 2023: are the information still up to date?

Yes, the information are still correct.
For the BLAS / LAPACK backends, I recently updated the README of Ipopt to explain how to use other backends (like BLIS) and that I recommend a sequential BLAS / LAPACK with MA86/MA97.

• after Support HSL solvers #113 is it really as simple as using HSL? what about importing HSL_jll licensing...? not clear to me

HSL_jll.jl has two versions, the official one that is available on the website of RAL and a dummy one that contains one C routine and that is available in the Julia registry.
Because we have the dummy version, HSL_jll.jl is an official Julia package and we can use it as any depenendency of a Julia package (like HSL.jl).
Thanks to the routine available in the dummy version, we can determine if the dummy or official HSL_jll.jl is installed. The C routine is LIBHSL_isfunctional and returns a boolean.
In HSL.jl, I did a Julia wrapper to easily check that at the Julia level:
https://github.com/JuliaSmoothOptimizers/HSL.jl/blob/main/src/C/libhsl.jl#L1
and also take care of loading a BLAS / LAPACK backend for the user if he didn't loaded one before (like MKL or AppleAccelerate):
https://github.com/JuliaSmoothOptimizers/HSL.jl/blob/main/src/HSL.jl#L16

Conclusion: If the user has installed the official HSL_jll.jl, we automatically use MA57 here otherwise we use MUMPS as a fallback but everything is transparent for the user because we can check which version of HSL_jll.jl is installed.

• what about using spral instead of mumps? (I had noticed Ipopt was compiled with it for >= 1.9)

SPRAL requires two environnement variables related to OpenMP. Because CTDirect.jl requires at least Julia 1.9, I don't check the Julia version but I ensure that the two environnement variables are well defined if the user wants to use it.
I think MUMPS is more efficient than SPRAL so I don't dispatch to it by default if we can't use the HSL solvers but we need some benchmarks to verify that.

Note that I compiled MUMPS 5.7.2 yesterday with Yggdrasil. I wait a new release of SPRAL and I will recompile Ipopt with the updated linear solvers.

[jbcaillau]

🙏🏽 @amontoison now I get it. Actually, you probably already the dummy routine trick, but did not remembered it so was not able to figure where licensing was coming into play. having a try (before forgetting it again 🤞🏽)

[jbcaillau]

@ocots @PierreMartinon @0Yassine0 As discussed now with @amontoison 🙏🏽: for sparse problems, typically for constraints resulting from discretising $\dot{x} = f(x,u)$ according to sth like

x[i+1] == x[i] + h[i] / 2 * ( f(x[i], u[i]) + f(x[i]+1, u[i+1]) )

one has very sparse Hessians and @amontoison suggests that ma27 should perform better than ma57. TBT.

[jbcaillau]

@0Yassine0 @ocots I installed HSL lib on a linux server, dev-ed HSL_jll.jl on the corresponding path but still git the following error (falling back to MUMPS, then):

Can you please share how you did the install end to end?

[0Yassine0]

@0Yassine0 @ocots I installed HSL lib on a linux server, dev-ed HSL_jll.jl on the corresponding path but still git the following error (falling back to MUMPS, then):

Can you please share how you did the install end to end?

I didn't use HSLpackage from @JuliaSmoothOptimizers for the matter.
The installation steps are explained here in details.
To use it with OptimalControl:
import HSL_jll
sol = OptimalControl.solve(Model,linear_solver="maXX" , hsllib=HSL_jll.libhsl_path)

[amontoison]

@0Yassine0 @ocots I installed HSL lib on a linux server, dev-ed HSL_jll.jl on the corresponding path but still git the following error (falling back to MUMPS, then):

Can you please share how you did the install end to end?

@jbcaillau You have an old version of HSL_jll.jl (the first one!!), please download the last one.
We renamed JULIAHSL into LIBHSL and in your version you only have the routine JULIAHSL_isfunctional.
You can directly download the new version of HSL_jll if you go to my order on the HSL website.

[jbcaillau]

@0Yassine0 @ocots I installed HSL lib on a linux server, dev-ed HSL_jll.jl on the corresponding path but still git the following error (falling back to MUMPS, then):
Can you please share how you did the install end to end?

@jbcaillau You have an old version of HSL_jll.jl (the first one!!), please download the last one. We renamed JULIAHSL into LIBHSL and in your version you only have the routine JULIAHSL_isfunctional. You can directly download the new version of HSL_jll if you go to my order on the HSL website.

@amontoison outdated me 🙃

MKL + HSL set 🙏🏽

[PierreMartinon]

Started a basic benchmarking script. SPRAL does seems slower than MUMPS

julia> include("test/benchmark.jl")
Profile: linear_solver=mumps

Settings: tol=1e-08 grid_size=100 precompile=true

Precompilation step
goddard param 

Benchmark step
goddard              completed in   1.81 s
param                completed in   2.74 s

Total time (s):   4.56 

julia> include("test/benchmark.jl")
Profile: linear_solver=spral

Settings: tol=1e-08 grid_size=100 precompile=true

Precompilation step
goddard param 

Benchmark step
goddard              completed in   9.69 s
param                completed in  12.20 s

Total time (s):  21.90 

[jbcaillau]

@PierreMartinon nice. what about ma57 and ma27 (the last might be better for our problems according to @amontoison)?

[jbcaillau]

@PierreMartinon I suggest to

• go back to MUMPS as default (same as Ipopt default), even if HSL is installed; for instance, our Goddard tuto does not work properly (small change needed but still...) with ma57
• AFAIK passing symbols (:mumps...) instead of strings is more (we can also keep the string way for compatibility, not sure where it already appears...)

[ocots]

I am ok with symbols. We can provide with the description.

[amontoison]

@jbcaillau What is the issue with MA57?
Is it an error in the linear solver?

[jbcaillau]

@amontoison just a slightly deprecated / less precise convergence for this particular case that is not compatible with some further solution analysis.

[amontoison]

Did you tune the tolerance of the optimization solver?
I don't think that the issue is the linear solver.

[jbcaillau]

only difference between the two runs = linear solver. same params otherwise (tols, etc.)

as we cannot be sure that a user will have HSL installed, it is safer to ensure that the default (= MUMPS) for Ipopt gives the requested precision for post computation. minor point, though.

[PierreMartinon]

I am ok with symbols. We can provide with the description.

ok for symbols as well (should we set this as a general rule ?) Equality tests may have to be updated for instance in https://github.com/control-toolbox/CTDirect.jl/blob/a9ae483a57382d44dba0612af914d78275fd6d7a/ext/CTSolveExt.jl

As for the default choice, this is a 1-line change here, so pretty straightforward

CTDirect.jl/src/CTDirect.jl

Line 21 in a9ae483

const __linear_solver() = "ma57"