Skip to content

A Julia implementation of incomplete LU factorization with zero level of fill in.

License

Notifications You must be signed in to change notification settings

mcovalt/ILUZero.jl

Folders and files

NameName
Last commit message
Last commit date

Latest commit

 

History

56 Commits
 
 
 
 
 
 
 
 
 
 
 
 

Repository files navigation

ILUZero.jl

ILUZero.jl is a Julia implementation of incomplete LU factorization with zero level of fill-in. It allows for non-allocating updates of the factorization.

Requirements

  • Julia 1.0 and up

Installation

julia> ]
pkg> add ILUZero

Why use ILUZero.jl?

This package was created to exploit some specific properties of a problem I had:

  • the non-zero elements of the matrix were frequently updated
  • an approximation of the matrix inverse was needed each time the non-zero elements were changed
  • the sparsity structure of the matrix was never altered
  • the non-zero elements of the matrix weren't necessarily floating point numbers

The ILU(0) preconditioner was a great solution for this problem. The method is simple enough to write a performant version purely in Julia which allowed using Julia's flexible type system for the elements. Also, ILU(0) has constant memory requirements so the updating matrix could re-use previously allocated structures*. Win win (win win)!

*Hint: take a look at ConjugateGradients.jl if this type of preconditioner would be beneficial to your project. Combined, these packages can help reduce allocations in those hot paths.

How to use

julia> using ILUZero
  • LU = ilu0(A): Create a factorization based on a sparse matrix A
  • ilu0!(LU, A): Update factorization LU in-place based on a sparse matrix A. This assumes the original factorization was created with another sparse matrix with the exact same sparsity pattern as A. No check is made for this.
  • To solve for x in (LU)x=b, use the same methods as you typically would: \ or ldiv!(x, LU, b). See the docs for further information.
  • There's also:
    • Forward substitution: forward_substitution!(y, LU, b) solves L\b and stores the solution in y.
    • Backward substitution: backward_substitution!(x, LU, y) solves U\y and stores the solution in x.
    • Nonzero count: nnz(LU) will return the number of nonzero entries in LU.

Performance

julia> using ILUZero
julia> using BenchmarkTools, LinearAlgebra, SparseArrays
julia> A = sprand(1000, 1000, 5 / 1000) + 10I
julia> fact = @btime ilu0(A)
       119.475 μs (21 allocations: 165.31 KiB)
julia> updated_fact = @btime ilu0!($fact, $A)
       83.079 μs (0 allocations: 0 bytes)