Add the code of IncompleteLU.jl in KrylovPreconditioners.jl

src/KrylovPreconditioners.jl

@@ -43,6 +43,7 @@ export TriangularOperator
 include("ic0.jl")
 include("ilu0.jl")
 include("blockjacobi.jl")
+include("ilu/IncompleteLU.jl")
 
 # Scaling
 include("scaling.jl")

src/ilu/IncompleteLU.jl

@@ -0,0 +1,15 @@
+using LinearAlgebra: Factorization, AdjointFactorization, LowerTriangular, UnitLowerTriangular, UpperTriangular
+using SparseArrays
+using Base: @propagate_inbounds
+
+struct ILUFactorization{Tv,Ti} <: Factorization{Tv}
+    L::SparseMatrixCSC{Tv,Ti}
+    U::SparseMatrixCSC{Tv,Ti}
+end
+
+include("sorted_set.jl")
+include("linked_list.jl")
+include("sparse_vector_accumulator.jl")
+include("insertion_sort_update_vector.jl")
+include("application.jl")
+include("crout_ilu.jl")

src/ilu/application.jl

@@ -0,0 +1,221 @@
+import SparseArrays: nnz
+import LinearAlgebra: ldiv!
+import Base.\
+
+export forward_substitution!, backward_substitution!
+export adjoint_forward_substitution!, adjoint_backward_substitution!
+
+"""
+Returns the number of nonzeros of the `L` and `U`
+factor combined.
+
+Excludes the unit diagonal of the `L` factor,
+which is not stored.
+"""
+nnz(F::ILUFactorization) = nnz(F.L) + nnz(F.U)
+
+function ldiv!(F::ILUFactorization, y::AbstractVecOrMat)
+    forward_substitution!(F, y)
+    backward_substitution!(F, y)
+end
+
+function ldiv!(F::AdjointFactorization{<:Any,<:ILUFactorization}, y::AbstractVecOrMat)
+    adjoint_forward_substitution!(F.parent, y)
+    adjoint_backward_substitution!(F.parent, y)
+end
+
+function ldiv!(y::AbstractVector, F::ILUFactorization, x::AbstractVector)
+    y .= x
+    ldiv!(F, y)
+end
+
+function ldiv!(y::AbstractVector, F::AdjointFactorization{<:Any,<:ILUFactorization}, x::AbstractVector)
+    y .= x
+    ldiv!(F, y)
+end
+
+function ldiv!(y::AbstractMatrix, F::ILUFactorization, x::AbstractMatrix)
+    y .= x
+    ldiv!(F, y)
+end
+
+function ldiv!(y::AbstractMatrix, F::AdjointFactorization{<:Any,<:ILUFactorization}, x::AbstractMatrix)
+    y .= x
+    ldiv!(F, y)
+end
+
+(\)(F::ILUFactorization, y::AbstractVecOrMat) = ldiv!(F, copy(y))
+(\)(F::AdjointFactorization{<:Any,<:ILUFactorization}, y::AbstractVecOrMat) = ldiv!(F, copy(y))
+
+"""
+Applies in-place backward substitution with the U factor of F, under the assumptions:
+
+1. U is stored transposed / row-wise
+2. U has no lower-triangular elements stored
+3. U has (nonzero) diagonal elements stored.
+"""
+function backward_substitution!(F::ILUFactorization, y::AbstractVector)
+    U = F.U
+    @inbounds for col = U.n : -1 : 1
+
+        # Substitutions
+        for idx = U.colptr[col + 1] - 1 : -1 : U.colptr[col] + 1
+            y[col] -= U.nzval[idx] * y[U.rowval[idx]]
+        end
+
+        # Final value for y[col]
+        y[col] /= U.nzval[U.colptr[col]]
+    end
+
+    y
+end
+
+function backward_substitution!(F::ILUFactorization, y::AbstractMatrix)
+    U = F.U
+    p = size(y, 2)
+
+    @inbounds for c = 1 : p
+        @inbounds for col = U.n : -1 : 1
+
+            # Substitutions
+            for idx = U.colptr[col + 1] - 1 : -1 : U.colptr[col] + 1
+                y[col,c] -= U.nzval[idx] * y[U.rowval[idx],c]
+            end
+
+            # Final value for y[col,c]
+            y[col,c] /= U.nzval[U.colptr[col]]
+        end
+    end
+
+    y
+end
+
+function backward_substitution!(v::AbstractVector, F::ILUFactorization, y::AbstractVector)
+    v .= y
+    backward_substitution!(F, v)
+end
+
+function backward_substitution!(v::AbstractMatrix, F::ILUFactorization, y::AbstractMatrix)
+    v .= y
+    backward_substitution!(F, v)
+end
+
+function adjoint_backward_substitution!(F::ILUFactorization, y::AbstractVector)
+    L = F.L
+    @inbounds for col = L.n - 1 : -1 : 1
+        # Substitutions
+        for idx = L.colptr[col + 1] - 1 : -1 : L.colptr[col]
+            y[col] -= L.nzval[idx] * y[L.rowval[idx]]
+        end
+    end
+
+    y
+end
+
+function adjoint_backward_substitution!(F::ILUFactorization, y::AbstractMatrix)
+    L = F.L
+    p = size(y, 2)
+    @inbounds for c = 1 : p
+        @inbounds for col = L.n - 1 : -1 : 1
+            # Substitutions
+            for idx = L.colptr[col + 1] - 1 : -1 : L.colptr[col]
+                y[col,c] -= L.nzval[idx] * y[L.rowval[idx],c]
+            end
+        end
+    end
+
+    y
+end
+
+function adjoint_backward_substitution!(v::AbstractVector, F::ILUFactorization, y::AbstractVector)
+    v .= y
+    adjoint_backward_substitution!(F, v)
+end
+
+function adjoint_backward_substitution!(v::AbstractMatrix, F::ILUFactorization, y::AbstractMatrix)
+    v .= y
+    adjoint_backward_substitution!(F, v)
+end
+
+"""
+Applies in-place forward substitution with the L factor of F, under the assumptions:
+
+1. L is stored column-wise (unlike U)
+2. L has no upper triangular elements
+3. L has *no* diagonal elements
+"""
+function forward_substitution!(F::ILUFactorization, y::AbstractVector)
+    L = F.L
+    @inbounds for col = 1 : L.n - 1
+        for idx = L.colptr[col] : L.colptr[col + 1] - 1
+            y[L.rowval[idx]] -= L.nzval[idx] * y[col]
+        end
+    end
+
+    y
+end
+
+function forward_substitution!(F::ILUFactorization, y::AbstractMatrix)
+    L = F.L
+    p = size(y, 2)
+    @inbounds for c = 1 : p
+        @inbounds for col = 1 : L.n - 1
+            for idx = L.colptr[col] : L.colptr[col + 1] - 1
+                y[L.rowval[idx],c] -= L.nzval[idx] * y[col,c]
+            end
+        end
+    end
+
+    y
+end
+
+function forward_substitution!(v::AbstractVector, F::ILUFactorization, y::AbstractVector)
+    v .= y
+    forward_substitution!(F, v)
+end
+
+function forward_substitution!(v::AbstractMatrix, F::ILUFactorization, y::AbstractMatrix)
+    v .= y
+    forward_substitution!(F, v)
+end
+
+function adjoint_forward_substitution!(F::ILUFactorization, y::AbstractVector)
+    U = F.U
+    @inbounds for col = 1 : U.n
+        # Final value for y[col]
+        y[col] /= U.nzval[U.colptr[col]]
+
+        for idx = U.colptr[col] + 1 : U.colptr[col + 1] - 1
+            y[U.rowval[idx]] -= U.nzval[idx] * y[col]
+        end
+    end
+
+    y
+end
+
+function adjoint_forward_substitution!(F::ILUFactorization, y::AbstractMatrix)
+    U = F.U
+    p = size(y, 2)
+    @inbounds for c = 1 : p
+        @inbounds for col = 1 : U.n
+            # Final value for y[col,c]
+            y[col,c] /= U.nzval[U.colptr[col]]
+
+            for idx = U.colptr[col] + 1 : U.colptr[col + 1] - 1
+                y[U.rowval[idx],c] -= U.nzval[idx] * y[col,c]
+            end
+        end
+    end
+
+    y
+end
+
+function adjoint_forward_substitution!(v::AbstractVector, F::ILUFactorization, y::AbstractVector)
+    v .= y
+    adjoint_forward_substitution!(F, v)
+end
+
+function adjoint_forward_substitution!(v::AbstractMatrix, F::ILUFactorization, y::AbstractMatrix)
+    v .= y
+    adjoint_forward_substitution!(F, v)
+end

src/ilu/crout_ilu.jl

@@ -0,0 +1,118 @@
+export ilu
+
+struct ILUFactorization{Tv,Ti} <: Factorization{Tv}
+    L::SparseMatrixCSC{Tv,Ti}
+    U::SparseMatrixCSC{Tv,Ti}
+end
+
+function lutype(T::Type)
+    UT = typeof(oneunit(T) - oneunit(T) * (oneunit(T) / (oneunit(T) + zero(T))))
+    LT = typeof(oneunit(UT) / oneunit(UT))
+    S = promote_type(T, LT, UT)
+end
+
+function ilu(A::SparseMatrixCSC{ATv,Ti}; τ = 1e-3) where {ATv,Ti}
+    n = size(A, 1)
+
+    Tv = lutype(ATv)
+
+    L = spzeros(Tv, Ti, n, n)
+    U = spzeros(Tv, Ti, n, n)
+
+    U_row = SparseVectorAccumulator{Tv,Ti}(n)
+    L_col = SparseVectorAccumulator{Tv,Ti}(n)
+
+    A_reader = RowReader(A)
+    L_reader = RowReader(L, Val{false})
+    U_reader = RowReader(U, Val{false})
+
+    @inbounds for k = Ti(1) : Ti(n)
+
+        ##
+        ## Copy the new row into U_row and the new column into L_col
+        ##
+
+        col::Int = first_in_row(A_reader, k)
+
+        while is_column(col)
+            add!(U_row, nzval(A_reader, col), col)
+            next_col = next_column(A_reader, col)
+            next_row!(A_reader, col)
+
+            # Check if the next nonzero in this column
+            # is still above the diagonal
+            if has_next_nonzero(A_reader, col) && nzrow(A_reader, col) ≤ col
+                enqueue_next_nonzero!(A_reader, col)
+            end
+
+            col = next_col
+        end
+
+        # Copy the remaining part of the column into L_col
+        axpy!(one(Tv), A, k, nzidx(A_reader, k), L_col)
+
+        ##
+        ## Combine the vectors:
+        ##
+
+        # U_row[k:n] -= L[k,i] * U[i,k:n] for i = 1 : k - 1
+        col = first_in_row(L_reader, k)
+
+        while is_column(col)
+            axpy!(-nzval(L_reader, col), U, col, nzidx(U_reader, col), U_row)
+
+            next_col = next_column(L_reader, col)
+            next_row!(L_reader, col)
+
+            if has_next_nonzero(L_reader, col)
+                enqueue_next_nonzero!(L_reader, col)
+            end
+
+            col = next_col
+        end
+
+        # Nothing is happening here when k = n, maybe remove?
+        # L_col[k+1:n] -= U[i,k] * L[i,k+1:n] for i = 1 : k - 1
+        if k < n
+            col = first_in_row(U_reader, k)
+
+            while is_column(col)
+                axpy!(-nzval(U_reader, col), L, col, nzidx(L_reader, col), L_col)
+
+                next_col = next_column(U_reader, col)
+                next_row!(U_reader, col)
+
+                if has_next_nonzero(U_reader, col)
+                    enqueue_next_nonzero!(U_reader, col)
+                end
+
+                col = next_col
+            end
+        end
+
+        ## 
+        ## Apply a drop rule
+        ##
+
+        U_diag_element = U_row.nzval[k]
+        # U_diag_element = U_row.values[k]
+
+        # Append the columns        
+        append_col!(U, U_row, k, τ)
+        append_col!(L, L_col, k, τ, inv(U_diag_element))
+
+        # Add the new row and column to U_nonzero_col, L_nonzero_row, U_first, L_first
+        # (First index *after* the diagonal)
+        U_reader.next_in_column[k] = U.colptr[k] + 1
+        if U.colptr[k] < U.colptr[k + 1] - 1
+            enqueue_next_nonzero!(U_reader, k)
+        end
+
+        L_reader.next_in_column[k] = L.colptr[k]
+        if L.colptr[k] < L.colptr[k + 1]
+            enqueue_next_nonzero!(L_reader, k)
+        end
+    end
+
+    return ILUFactorization(L, U)
+end