Identifying the systems available linear solvers (With JuMP & Ipopt)

Domains Optimization (Mathematical)
#1

Hi all,

I’d like to be able to verify / list available linear solvers on any given machine before running a model. My stack is JuMP and Ipopt, which means I can do something like this:

import Ipopt; 
import JuMP; 
optimizer = JuMP.with_optimizer(Ipopt.Optimizer, linear_solver="ma97"); 
model = JuMP.Model(optimizer)

Here, optimizer is of type JuMP.OptimizerFactory, and seems to be lazily evaluated since at this point in time model has no issues even if one types garbage into the linear_solver string.
From my understanding, JuMP sends this factory info to Ipopt.jl, which uses MOI to then do a bunch of ccalls to run everything once JuMP.optimize!(model) is called, and not before.

By default Ipopt.jl packages the MUMPS linear solver, although there are more efficient ones out there: HSL MA97 as an example (as set in optimizer above).

In my specific case, I want to be able to prioritize certain linear solvers over others, depending on solvers available on a users system. By default, Ipopt will use MA27 if it’s compiled with HSL, but I’m needing MA97 if available and would fall back to MUMPS if not.

The only ways I can find to check if MA97 is extant at the moment is to run JuMP.optimize!(model) and wait for it to fail via the MOI ccall attempt via addOption at the JuMP level.

Which ultimately could be simplified to invoking something at the Ipopt.jl level by building some dummy IpoptProblem and capture the error:

prob = Ipopt.createProblem(1,[1.],[1.],1,[1.],[1.],1,1,sum,sum,sum,sum);
try 
    Ipopt.addOption(prob, "linear_solver", "ma97")
    println("HSL installed")
catch 
    println("Use MUMPS")
end

Is there a cleaner way to do this at all that I’ve just overlooked? This of course gets pretty messy if I want a list of installed solvers rather than just one check and a fallback.

Update: as a followup - I’ve realized this solution isn’t even enough. The addOption call will pass so long as it is possible to be accepted in Ipopt’s settings here. For example, I don’t have “pardiso” installed, but I won’t receive an error at this point if I try to test with it. Thus I need to actually call Ipopt.solveProblem(prob) as well and capture on that output instead.

#2

I’m not sure there is a good solution. Trying to solve the dummy model via the low-level interface is probably the easiest way forward.

#3

Thanks @odow. Just wanted to make sure I wasn’t doing something overtly complicated for no reason.

For completeness, here’s my final solution to this issue.

function linearSolver(solver_name::String = "ma97")
    prob = Ipopt.createProblem(1,[1.],[1.],1,[1.],[1.],1,1,sum,sum,sum,sum);
    Ipopt.addOption(prob, "sb", "yes");
    Ipopt.addOption(prob, "print_level", 0);
    # Initially, we must check that coinhsl is installed at all if we want a HSL solver.
    # If not, Ipopt will default to try and find any dynamically linked libhsl. If it
    # cannot find one it will hard panic and we can't capture that failure.
    if occursin("ma", solver_name)
        Ipopt.addOption(prob, "linear_solver", "ma27");
        try
            # No HSL was found on the system, we return with fallback.
            return runLinearSolverCheck(prob, solver_name)
        catch
            #Continue, we have access to at least some version of HSL
        end
    end
    try
        # Outer try will fail if solver string is not in the list of
        # possible Ipopt solvers.
        # For now that's ma27, ma57, ma77, ma86, ma97, pardiso, wsmp, mumps, custom
        Ipopt.addOption(prob, "linear_solver", solver_name);
        try
            # Inner try attempts to run the dummy program and will crash because
            # the dummy is malformed.
            # No error code will be returned if the solver is extant, and an Invalid_Option
            # if the library is not found.
            runLinearSolverCheck(prob, solver_name)
        catch
            # We can use the requested solver.
            solver_name
        end
    catch
        "mumps"
    end
end

function runLinearSolverCheck(prob::IpoptProblem, solver_name::String)
    result_code = Ipopt.solveProblem(prob);
    if Ipopt.ApplicationReturnStatus[result_code] == :Invalid_Option
        @info "Unable to set linear_solver = $(solver_name), defaulting to MUMPS."
        return "mumps"
    else
        error("Attempts to identify linear solvers on system returned unexpected results.");
    end
end

Linear solvers installed on my system: ma*, mumps. Missing pardiso, wsmp.

julia> linearSolver()
"ma97"

julia> linearSolver("ma27")
"ma27"

julia> linearSolver("pardiso")
┌ Warning: Unable to set linear_solver = pardiso, defaulting to MUMPS.
└ @ LS ~/testing/src/ls.jl:56
"mumps"

julia> linearSolver("mumps")
"mumps"
#4

This is a good question for the Ipopt mailing list. If Ipopt has an API for this, we can use it. If not, then your solution is the best that one can do.

1 Like
#5

Libbum, can I use your code? How it is licensed? I need to do something like this in a code I am writing.

#6

Yes, certainly! You can find everything I currently use in DICE.jl, which is MIT licensed.