# Problem with IncompleteLU.jl

**URL:** <https://discourse.julialang.org/t/problem-with-incompletelu-jl/90636>\
**Category:** Numerics\
**Created:** [November 22, 2022, 10:01am UTC](https://discourse.julialang.org/t/problem-with-incompletelu-jl/90636 "2022-11-22T10:01:13Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![nvenkov1](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nvenkov1/32/43679_2.png) [@nvenkov1](https://discourse.julialang.org/u/nvenkov1)\
**Post date:** [November 22, 2022, 10:01am UTC](https://discourse.julialang.org/t/problem-with-incompletelu-jl/90636/1 "2022-11-22T10:01:13Z")

</div>

I use the file transient.mtx from the SuiteSparse Matrix Collection, see [https://suitesparse-collection-website.herokuapp.com/MM/Freescale/transient.tar.gz](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:

```julia
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?

---

<div class="post-metadata">

**Author:** ![nvenkov1](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nvenkov1/32/43679_2.png) [@nvenkov1](https://discourse.julialang.org/u/nvenkov1)\
**Post date:** [November 22, 2022, 10:41am UTC](https://discourse.julialang.org/t/problem-with-incompletelu-jl/90636/2 "2022-11-22T10:41:50Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![nvenkov1](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nvenkov1/32/43679_2.png) [@nvenkov1](https://discourse.julialang.org/u/nvenkov1)\
**Post date:** [November 22, 2022, 1:39pm UTC](https://discourse.julialang.org/t/problem-with-incompletelu-jl/90636/3 "2022-11-22T13:39:00Z")

</div>

The following fix seems to work:

```julia
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.

---

<div class="post-metadata">

**Author:** ![amontoison](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/amontoison/32/218741_2.png) [@amontoison](https://discourse.julialang.org/u/amontoison)\
**Post date:** [November 22, 2022, 7:09pm UTC](https://discourse.julialang.org/t/problem-with-incompletelu-jl/90636/4 "2022-11-22T19:09:28Z")

</div>

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](https://juliasmoothoptimizers.github.io/Krylov.jl/dev/gpu/#Example-with-a-general-square-system)  
It’s a copy of the [`mc21`](https://www.hsl.rl.ac.uk/catalogue/mc21.html) algorithm for information.

---

<div class="post-metadata">

**Author:** ![amontoison](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/amontoison/32/218741_2.png) [@amontoison](https://discourse.julialang.org/u/amontoison)\
**Post date:** [November 22, 2022, 9:57pm UTC](https://discourse.julialang.org/t/problem-with-incompletelu-jl/90636/5 "2022-11-22T21:57:22Z")

</div>

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

```julia
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)

```