MKLPardisoKKTSolver have inconsistent behavior vs default KKTSolver

[Sage0614]

Hi,

I find out a corner case where using MKLPardisoKKTSolver have different behavior vs. base solver

`
using LinearAlgebra, SparseArrays, Random, COSMO, JuMP, Test
using Pardiso

rng = Random.MersenneTwister(1)
k = 200; # number of factors
n = 2500; # number of assets
D = spdiagm(0 => rand(rng, n) .* sqrt(k))
F = sprandn(rng, n, k, 0.5); # factor loading matrix
μ = (3 .+ 9. * rand(rng, n)) / 100. # expected returns between 3% - 12%
γ = 1.0; # risk aversion parameter
d = 1 # we are starting from all cash
x0 = zeros(n);

model = JuMP.Model(optimizer_with_attributes(COSMO.Optimizer,
"rho" => 3.0,
"alpha" =>1.0,
"eps_abs" => 1e-6,
"eps_rel" => 1e-5,
"accelerator" => with_options(AndersonAccelerator, mem = 15) ));

set_optimizer_attribute(model,"kkt_solver",with_options(MKLPardisoKKTSolver))

Mt = [D.^0.5; F']
a = 1e-3
b = 1e-1

@variable(model, x[1:n]);
@variable(model, y[1:k]);
@variable(model, t[1:n]);
@variable(model, s[1:n]);

@objective(model, Min, x' * D * x + y' * y - 1/γ * μ' * x);

@constraint(model, y .== F' * x);
@constraint(model,[1;x] in MOI.NormOneCone(1+n));

@constraint(model, sum(x) + a * sum(s) + b * sum(t) == d + sum(x0) );
@constraint(model, [i = 1:n], x[i] - x0[i] <= s[i]); # model the absolute value with slack variable s
@constraint(model, [i = 1:n], x0[i] - x[i] <= s[i]);
@constraint(model, [i = 1:n], [t[i], 1, x[i] - x0[i]] in MOI.PowerCone(2/3));
JuMP.optimize!(model)
`
The result converged after 350 iters. however, if you start from the converged solution, and using default QDLDL solver, the program knows the solution is optimal and stopped after 1 iter:

`

      COSMO v0.8.6 - A Quadratic Objective Conic Solver
                     Michael Garstka
            University of Oxford, 2017 - 2022

---

Problem: x ∈ R^{10200},
constraints: A ∈ R^{17702x10200} (285005 nnz),
matrix size to factor: 27902x27902,
Floating-point precision: Float64
Sets: Nonnegatives of dim: 10001
ZeroSet of dim: 201
PowerCone of dim: 3
PowerCone of dim: 3
PowerCone of dim: 3
... and 2497 more
Settings: ϵ_abs = 1.0e-06, ϵ_rel = 1.0e-05,
ϵ_prim_inf = 1.0e-04, ϵ_dual_inf = 1.0e-04,
ρ = 3, σ = 1e-06, α = 1,
max_iter = 5000,
scaling iter = 10 (on),
check termination every 25 iter,
check infeasibility every 40 iter,
KKT system solver: QDLDL
Acc: Anderson Type2{QRDecomp},
Memory size = 15, RestartedMemory,
Safeguarded: true, tol: 2.0
Setup Time: 1.54ms

Iter: Objective: Primal Res: Dual Res: Rho:
1 -9.9443e-02 1.0155e-05 1.2928e-06 3.0000e+00

---

Results
Status: Solved
Iterations: 1
Optimal objective: -0.09944
Runtime: 0.022s (22.0ms)
`

if you are using the MKLPardiso solver, turns out it stalls there without any improvement until time out. truncated output below:

`

---

      COSMO v0.8.6 - A Quadratic Objective Conic Solver
                     Michael Garstka
            University of Oxford, 2017 - 2022

---

Problem: x ∈ R^{10200},
constraints: A ∈ R^{17702x10200} (285005 nnz),
matrix size to factor: 27902x27902,
Floating-point precision: Float64
Sets: Nonnegatives of dim: 10001
ZeroSet of dim: 201
PowerCone of dim: 3
PowerCone of dim: 3
PowerCone of dim: 3
... and 2497 more
Settings: ϵ_abs = 1.0e-06, ϵ_rel = 1.0e-05,
ϵ_prim_inf = 1.0e-04, ϵ_dual_inf = 1.0e-04,
ρ = 3, σ = 1e-06, α = 1,
max_iter = 5000,
scaling iter = 10 (on),
check termination every 25 iter,
check infeasibility every 40 iter,
KKT system solver: MKL Pardiso
Acc: Anderson Type2{QRDecomp},
Memory size = 15, RestartedMemory,
Safeguarded: true, tol: 2.0
Setup Time: 1.34ms

Iter: Objective: Primal Res: Dual Res: Rho:
1 0.0000e+00 1.0000e+00 1.1992e-01 3.0000e+00
25 0.0000e+00 1.0000e+00 1.1992e-01 3.0000e+00
50 0.0000e+00 1.0000e+00 1.1992e-01 3.0000e+00
75 0.0000e+00 1.0000e+00 1.1992e-01 3.0000e+00
100 0.0000e+00 1.0000e+00 1.1992e-01 3.0000e+00
125 0.0000e+00 1.0000e+00 1.1992e-01 3.0000e+00
150 0.0000e+00 1.0000e+00 1.1992e-01 3.0000e+00
175 0.0000e+00 1.0000e+00 1.1992e-01 3.0000e+00
200 0.0000e+00 1.0000e+00 1.1992e-01 3.0000e+00
225 0.0000e+00 1.0000e+00 1.1992e-01 3.0000e+00
250 0.0000e+00 1.0000e+00 1.1992e-01 3.0000e+00
275 0.0000e+00 1.0000e+00 1.1992e-01 3.0000e+00
300 0.0000e+00 1.0000e+00 1.1992e-01 3.0000e+00
325 0.0000e+00 1.0000e+00 1.1992e-01 3.0000e+00
350 0.0000e+00 1.0000e+00 1.1992e-01 3.0000e+00
375 0.0000e+00 1.0000e+00 1.1992e-01 3.0000e+00
400 0.0000e+00 1.0000e+00 1.1992e-01 3.0000e+00
425 0.0000e+00 1.0000e+00 1.1992e-01 3.0000e+00
450 0.0000e+00 1.0000e+00 1.1992e-01 3.0000e+00
475 0.0000e+00 1.0000e+00 1.1992e-01 3.0000e+00
500 0.0000e+00 1.0000e+00 1.1992e-01 3.0000e+00
`
Thanks.

[migarstka]

Thanks for reporting. This could be a problem with the way the warm starting is implemented with MKLPardiso. Will have a look at this.

[oisinfaust]

I also ran into this problem. I think it's not related to warm starting, but to the call of free_memory! at the very end of the optimize! method in solver.jl. The only effect of this is to release the memory of a Pardiso solver if one is being used.
Perhaps it would be better to remove this line, and attach free_memory! as a finalizer to all new Workspace objects at construction?