CODEWITHSUNDEEP
julialinear-algebranumerical-methods
Smoop94
Smoop94

Reputation: 73

Solve Complex Symmetric Sparse Linear Systems in Julia

I am building a numerical code that requires solving complex, structurally symmetric linear systems.

An example of system I am tying to solve is the following:

A = [20.0          0        0          0
    0          0        0 -10.6458im
    0          0 8.333333          0
    0 -10.6458im        0          1]

b = [0.0, 6.473417506631299,  0.0, 0.0]

x = A\b

The expected result is approximately:

x = [0, 0.0571185, 0, 0.6080724329436303im]

Since I need to work with sparse systems I built the matrix A as sparse (sparseCSC format) by storing only the lower-triangular part. I then made it symmetric with the call:

A = Symmetric(A,:L)

The type returned by typeof is:

LinearAlgebra.Symmetric{ComplexF64, SparseArrays.SparseMatrixCSC{ComplexF64, Int64}}

Everything seems OK until I need to solve a linear system that involves said matrix (using x = A\b). The following error appears:

LoadError: MethodError: no method matching lu!(::SparseArrays.SparseMatrixCSC{ComplexF64, Int64}, ::Val{true}; check=true)
Closest candidates are:
  lu!(::StridedMatrix{T}, ::Union{Val{true}, Val{false}}; check) where T<:Union{Float32, 
Float64, ComplexF32, ComplexF64} at 
C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\LinearAlgebra\src\lu.jl:79
  lu!(::Union{LinearAlgebra.Hermitian{T, S}, LinearAlgebra.Symmetric{T, S}} where {T, S}, ::Union{Val{true}, Val{false}}; check) at 
C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\LinearAlgebra\src\lu.jl:88
  lu!(::StridedMatrix{T} where T, ::Union{Val{true}, Val{false}}; check) at 
C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\LinearAlgebra\src\lu.jl:130

Has anyone experienced this problem? Is there a more correct way to solve complex symmetric sparse systems in Julia? (I am on Julia 1.6.1 on Windows 10).

Upvotes: 1

Views: 237

Answers (1)

Przemyslaw Szufel
Przemyslaw Szufel

Reputation: 42264

This seems like a bug very almost identical to this one: https://github.com/JuliaLang/SuiteSparse.jl/issues/18

I believe that as of today the only option might be to drop the symmetry information when solving A for b:

julia> sparse(A)\b
4-element Vector{ComplexF64}:
                 0.0 + 0.0im
 0.05711852871025479 + 0.0im
                 0.0 + 0.0im
                -0.0 + 0.6080724329436303im

Upvotes: 1

Related Questions