Symmetric Indefinite Solvers (with block diagonal pivoting)

Specific Domains Numerics
, ,
1

Are there any modules that provide access to (sparse) symmetric indefinite solvers with Bunch-Kaufman-Parlett pivoting where I can access the L and D matrices, or would it be useful to build a module that does this? Right now, I have to do something like this.

As of now, is compiling IPOPT with the HSL linear solvers (for example, HSL_MA57) the easiest way to currently use the HSL routines? I saw a previous post build a custom IPOPT with these solvers, but if I only want to call the Bunch-Kaufman-Parlett LDLt factorization (MA57), can I access it through IPOPT, or is there a more direct way? Btw, as another alternative, Pardiso doesn’t yet allow you to access the L and D factors.

Thanks!

2

Have you seen bkfact?

help?> bkfact
search: bkfact bkfact!

  bkfact(A, uplo::Symbol=:U, symmetric::Bool=issymmetric(A), rook::Bool=false) -> BunchKaufman

  Compute the Bunch-Kaufman [^Bunch1977] factorization of a symmetric or
  Hermitian matrix A and return a BunchKaufman object. uplo indicates which
  triangle of matrix A to reference. If symmetric is true, A is assumed to be
  symmetric. If symmetric is false, A is assumed to be Hermitian. If rook is
  true, rook pivoting is used. If rook is false, rook pivoting is not used.
  The following functions are available for BunchKaufman objects: size, \,
  inv, issymmetric, ishermitian.

  | [^Bunch1977]
  |
  |  J R Bunch and L Kaufman, Some stable methods for calculating
  |  inertia and solving symmetric linear systems, Mathematics of
  |  Computation 31:137 (1977), 163-179. url.
3

Thanks, sorry, I wasn’t clear. I’m aware of bkfact and ldlfact, but neither allow me to access L and D separately (i.e., if I need to do a Hessian modification). I’ll clarify my question.

4

but neither allow me to access L and D

This is fixed on master where you can extract the factors with

julia> F = bkfact(A);

julia> F.U
10×10 UnitUpperTriangular{Float64,Array{Float64,2}}:
 1.0  -0.989314  1.76372    0.710251  -0.499019    …  -0.127294   0.546951  -0.374852
  ⋅    1.0       0.837403  -0.227568   0.715229        0.762732   0.216349   0.743333
  ⋅     ⋅        1.0       -0.958846  -0.0750025      -0.308083   0.115865  -0.00222697
  ⋅     ⋅         ⋅         1.0        0.00597047      0.876363  -0.129265   0.377775
  ⋅     ⋅         ⋅          ⋅         1.0             0.552132   0.909079   0.333753
  ⋅     ⋅         ⋅          ⋅          ⋅          …  -1.18999   -0.69973    0.0101069
  ⋅     ⋅         ⋅          ⋅          ⋅             -1.14554    0.269799   0.313673
  ⋅     ⋅         ⋅          ⋅          ⋅              1.0       -0.756874  -0.263642
  ⋅     ⋅         ⋅          ⋅          ⋅               ⋅         1.0        0.0
  ⋅     ⋅         ⋅          ⋅          ⋅               ⋅          ⋅         1.0

julia> F.D
10×10 Tridiagonal{Float64,Array{Float64,1}}:
 2.03091   0.0        ⋅         ⋅         ⋅     …   ⋅        ⋅         ⋅          ⋅
 0.0      -4.87095   0.0        ⋅         ⋅         ⋅        ⋅         ⋅          ⋅
  ⋅        0.0      -0.707841  0.0        ⋅         ⋅        ⋅         ⋅          ⋅
  ⋅         ⋅        0.0       1.13327   0.0        ⋅        ⋅         ⋅          ⋅
  ⋅         ⋅         ⋅        0.0      10.324      ⋅        ⋅         ⋅          ⋅
  ⋅         ⋅         ⋅         ⋅        0.0    …  0.0       ⋅         ⋅          ⋅
  ⋅         ⋅         ⋅         ⋅         ⋅        7.7207   0.0        ⋅          ⋅
  ⋅         ⋅         ⋅         ⋅         ⋅        0.0     -2.39402   0.0         ⋅
  ⋅         ⋅         ⋅         ⋅         ⋅         ⋅       0.0      -0.130385  -3.29071
  ⋅         ⋅         ⋅         ⋅         ⋅         ⋅        ⋅       -3.29071   -0.464019
3 Likes
5

Thanks, sorry I missed this. Does this do fill-reducing permutations / can it handle sparse matrices?

6

We only have a dense Bunch-Kaufman. The ldlfact works for sparse matrices and you can extract the factors but it only doesn’t fill-reducing permutations, i.e. no permutations during the numerical factorization.

1 Like
7

See HSL.jl. I wrote an interface to HSL_MA97 (essentially a multi-core version of HSL_MA57) that gives access to the factors and allows you to modify them. In the ma57 branch, there’s a similar interface to HSL_MA57 but a patch is required to access the factors. I’m not quite sure I’m allowed to distribute it (I should submit it to the HSL folks). In the mean time, I could send it privately. Feel free to open an issue there if anything doesn’t work right.

1 Like