Problem with IncompleteLU.jl

Specific Domains Numerics
1

I use the file transient.mtx from the SuiteSparse Matrix Collection, see https://suitesparse-collection-website.herokuapp.com/MM/Freescale/transient.tar.gz.
Then I try to define an incomplete LU factorization that I apply as a preconditioner as follows:

using MatrixMarket: mmread
using Random: seed!
using IncompleteLU: ilu

A = mmread("transient.mtx");
M = ilu(A);
seed!(1);
b = rand(A.n);
sum(isnan.(M \ b))

and the result I get is 9. What is wrong with IncompleteLU?

2

I just ran the same code again, and this time I got no NaN and I can successfully use M as a preconditioner for GMRES. The behavior of IncompleteLU.ilu does not seem to be deterministic.

3

The following fix seems to work:

using MatrixMarket: mmread
using Random: seed!
using IncompleteLU: ilu
using LinearAlgebra: I

A = mmread("transient.mtx");
dA = 10^-6
M = ilu(A+dA*I);
seed!(1);
b = rand(A.n);
sum(isnan.(M \ b))

Then I get zero NaN, and I can successfully use M in GMRES.

4

The transient matrix has a zero on its diagonal, you need to apply a permutation first to insure that you will not have zero pivots for the ILU factorization.
An other alternative is to add a small multiple of I.

If you have CUDA.jl, they have a routine for CPU matrices, it’s zfd and I use it before the computation of an ILU0 preconditioner: GPU support · Krylov.jl
It’s a copy of the mc21 algorithm for information.

1 Like
5

@nvenkov1
off-topic: If you use Suite Sparse Matrix Collection often, a Julia interface is available to easily download the matrices:

using SuiteSparseMatrixCollection, MatrixMarket, SparseArrays

ssmc = ssmc_db()
matrices = ssmc_matrices(ssmc, "Freescale", "transient")
paths = fetch_ssmc(matrices, format="MM")

path_transient = joinpath(paths[1], "transient.mtx")
M = MatrixMarket.mmread(path_transient)
1 Like