Arbitary precision for PSD Cone

[migarstka]

Is your feature request related to a problem? Please describe.

The positive semidefinite cone does not yet allow arbitrary precision types. However, it seems we now have a blueprint in Clarabel.jl on how this can be done. We can probably port the logic to COSMO's projection method.

[jeremypparker]

I would also like this. Currently the issue is that COSMO uses raw LAPACK syevr for eigendecomposition. It seems straightforward to replace this with a pure Julia alternative but presumably there was some reason for not doing this in the first place?

[goulart-paul]

My recollection is that we use a call to the lapack syevr directly since there is no obvious way of doing eigen-decomposition in Julia otherwise without allocating. Calling either of the julia eigen or eigen! functions eventually lands (I think) here:

https://github.com/JuliaLang/LinearAlgebra.jl/blob/9bc292de83f6178b4a92688da76a13e71faa3c58/src/lapack.jl#L5375

which allocates vectors of float and int work vectors on every call. For COSMO this is not ideal since the operations on cones, particularly projections to the PSD cone, dominate the computation time.

In Clarabel.jl we use GenericLinearAlgebra.jl for eigen decomposition of non-BLAS floating point types, which works fine. We just call the high level decomposition functions there though and don't worry about allocating, since in that case the decomposition is only a small fraction of the overall solve time anyway and the allocations don't really matter. I don't know how easy it would be to reach into the internals of GenericLinearAlgebra and preallocate workspace data once upfront for it. It's presumably possible at least.

[jeremypparker]

I note that very recently symmetric eigen! was added to GenericLinearAlgebra: JuliaLinearAlgebra/GenericLinearAlgebra.jl#145
I haven't investigated whether it allocates anything but it's at least partly in-place

As a user, I expect COSMO to be slow when using arbitrary precision. Therefore I maybe don't mind if there are lots of allocations in this case. Obviously we want to keep the performant version for Float64 and Float32, but perhaps those can be specialisations of _syevr! and the general version calls GenericLinearAlgebra?

[araujoms]

I'm sorry, the only thing I added to GenericLinearAlgebra was the missing dispatch for Symmetric{<:Real}, the algorithm itself was already there, it just required Hermitian. In any case, it is easy to call the internal GenericLinearAlgebra functions with pre-allocated workspace. It just doesn't make much of a difference, it still allocates quite a lot.

Also, there is already an open PR adding support for arbitrary precision types: #136