JuliaSparse / SparseArrays.jl

SparseArrays.jl is a Julia stdlib
https://sparsearrays.juliasparse.org/
Other
89 stars 51 forks source link

Obtaining D matrix from LDLt factorization not working #113

Open lijas opened 3 years ago

lijas commented 3 years ago
julia> using LinearAlgebra

julia> using SparseArrays

julia> M = Symmetric(sparse(rand(4,4)));

julia> F = LinearAlgebra.ldlt(M)
SuiteSparse.CHOLMOD.Factor{Float64}
type:    LDLt
method:  simplicial
maxnnz:  16
nnz:     10
success: true

julia> sparse(F.D)
ERROR: getindex not defined for SuiteSparse.CHOLMOD.FactorComponent{Float64, :D}
Stacktrace:
  [1] error(::String, ::Type)
    @ Base ./error.jl:42
  [2] error_if_canonical_getindex(::IndexCartesian, ::SuiteSparse.CHOLMOD.FactorComponent{Float64, :D}, ::Int64, ::Int64)
    @ Base ./abstractarray.jl:1177
  [3] getindex
    @ ./abstractarray.jl:1163 [inlined]
  [4] _getindex
    @ ./abstractarray.jl:1208 [inlined]
  [5] getindex
    @ ./abstractarray.jl:1164 [inlined]
  [6] iterate
    @ ./abstractarray.jl:1090 [inlined]
  [7] iterate
    @ ./abstractarray.jl:1088 [inlined]
  [8] SparseMatrixCSC{Float64, Int64}(M::SuiteSparse.CHOLMOD.FactorComponent{Float64, :D})
    @ SparseArrays /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/SparseArrays/src/sparsematrix.jl:616
  [9] convert
    @ /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/SparseArrays/src/sparsematrix.jl:667 [inlined]
 [10] sparse(A::SuiteSparse.CHOLMOD.FactorComponent{Float64, :D})
    @ SparseArrays /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/SparseArrays/src/sparsematrix.jl:692
 [11] top-level scope
    @ REPL[5]:1

julia> sparse(F.L)
ERROR: SuiteSparse.CHOLMOD.CHOLMODException("sparse: supported only for :LD on LDLt factorizations")
Stacktrace:
 [1] sparse(FC::SuiteSparse.CHOLMOD.FactorComponent{Float64, :L})
   @ SuiteSparse.CHOLMOD /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/SuiteSparse/src/cholmod.jl:1098
 [2] top-level scope
   @ REPL[6]:1

julia> Diagonal(sparse(F.LD)) #Work around
4×4 Diagonal{Float64, SparseVector{Float64, Int64}}:
 0.435803   ⋅          ⋅         ⋅ 
  ⋅        0.612003    ⋅         ⋅ 
  ⋅         ⋅        -0.909383   ⋅ 
  ⋅         ⋅          ⋅        0.0107779

Version info:

julia> versioninfo()
Julia Version 1.6.0-rc1
Commit a58bdd9010 (2021-02-06 15:49 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: Intel(R) Core(TM) i7-8665U CPU @ 1.90GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-11.0.1 (ORCJIT, skylake)
ViralBShah commented 2 years ago

Seems like we don't have methods to convert these to regular sparse matrices. These are pointers to the internal opaque structures.

julia> F.D
SuiteSparse.CHOLMOD.FactorComponent{Float64, :D}
type:    LDLt
method:  simplicial
maxnnz:  16
nnz:     10
success: true

julia> F.L
SuiteSparse.CHOLMOD.FactorComponent{Float64, :L}
type:    LDLt
method:  simplicial
maxnnz:  16
nnz:     10
success: true
rayegun commented 2 years ago

Is it breaking to convert them?

ViralBShah commented 2 years ago

We don't want to convert them - just make sure that the conversion works when explicitly asked for.

ViralBShah commented 8 months ago

Noting, that you can do:

julia> sparse(F.LD)
4×4 SparseMatrixCSC{Float64, Int64} with 10 stored entries:
 0.855515    ⋅          ⋅        ⋅ 
 0.0135824  0.396216    ⋅        ⋅ 
 0.477432   2.28558   -1.6451    ⋅ 
 0.434555   2.34267    1.04879  0.076564