Merge pull request #73 from SciML/prec

Even further simplify preconditioner interface

Author: ChrisRackauckas

docs/src/basics/FAQ.md

@@ -9,8 +9,9 @@ IterativeSolvers.jl computes the norm after the application of the left precondt
 hack the system via the following formulation:
 
 ```julia
-Pl = LinearSolve.InvDiagonalPreconditioner(weights)
-Pr = LinearSolve.DiagonalPreconditioner(weights)
+using LinearSolve, LinearAlgebra
+Pl = LinearSolve.InvPreconditioner(Diagonal(weights))
+Pr = Diagonal(weights)
 
 A = rand(n,n)
 b = rand(n)
@@ -24,8 +25,9 @@ can use `ComposePreconditioner` to apply the preconditioner after the applicatio
 of the weights like as follows:
 
 ```julia
-Pl = ComposePreconitioner(LinearSolve.InvDiagonalPreconditioner(weights),realprec)
-Pr = LinearSolve.DiagonalPreconditioner(weights)
+using LinearSolve, LinearAlgebra
+Pl = ComposePreconitioner(LinearSolve.InvPreconditioner(Diagonal(weights),realprec))
+Pr = Diagonal(weights)
 
 A = rand(n,n)
 b = rand(n)

docs/src/basics/Preconditioners.md

@@ -42,8 +42,9 @@ multiplies by the same values. This is commonly used in the case where if, inste
 of random, `s` is an approximation to the eigenvalues of a system.
 
 ```julia
+using LinearSolve, LinearAlgebra
 s = rand(n)
-Pl = LinearSolve.DiagonalPreconditioner(s)
+Pl = Diagonal(s)
 
 A = rand(n,n)
 b = rand(n)
@@ -52,22 +53,27 @@ prob = LinearProblem(A,b)
 sol = solve(prob,IterativeSolvers_GMRES(),Pl=Pl)
 ```
 
-## Pre-Defined Preconditioners
-
-To simplify the usage of preconditioners, LinearSolve.jl comes with many standard
-preconditioners written to match the required interface.
-
-- `DiagonalPreconditioner(s::Union{Number,AbstractVector})`: the diagonal
-  preconditioner, defined as a diagonal matrix `Diagonal(s)`.
-- `InvDiagonalPreconditioner(s::Union{Number,AbstractVector})`: the diagonal
-  preconditioner, defined as a diagonal matrix `Diagonal(1./s)`.
-- `ComposePreconditioner(prec1,prec2)`: composes the preconditioners to apply
-  `prec1` before `prec2`.
-
 ## Preconditioner Interface
 
 To define a new preconditioner you define a Julia type which satisfies the
 following interface:
 
-- `Base.eltype(::Preconditioner)`
-- `LinearAlgebra.ldiv!(::AbstractVector,::Preconditioner,::AbstractVector)`
+- `Base.eltype(::Preconditioner)` (Required only for Krylov.jl)
+- `LinearAlgebra.ldiv!(::AbstractVector,::Preconditioner,::AbstractVector)` and
+  `LinearAlgebra.ldiv!(::Preconditioner,::AbstractVector)`
+
+## Curated List of Pre-Defined Preconditioners
+
+The following preconditioners are tested to match the interface of LinearSolve.jl.
+
+- `ComposePreconditioner(prec1,prec2)`: composes the preconditioners to apply
+  `prec1` before `prec2`.
+- `InvPreconditioner(prec)`: inverts `mul!` and `ldiv!` in a preconditioner
+  definition as a lazy inverse.
+- `LinearAlgera.Diagonal(s::Union{Number,AbstractVector})`: the lazy Diagonal
+  matrix type of Base.LinearAlgebra. Used for efficient construction of a
+  diagonal preconditioner.
+- Other `Base.LinearAlgera` types: all define the full Preconditioner interface.
+- [IncompleteLU.ilu](https://github.com/haampie/IncompleteLU.jl): an implementation
+  of the incomplete LU-factorization preconditioner. This requires `A` as a
+  `SparseMatrixCSC`.

docs/src/tutorials/caching_interface.md

@@ -21,6 +21,8 @@ means of solving and resolving linear systems. To do this with LinearSolve.jl,
 you simply `init` a cache, `solve`, replace `b`, and solve again. This looks like:
 
 ```julia
+using LinearSolve
+
 n = 4
 A = rand(n,n)
 b1 = rand(n); b2 = rand(n)

src/preconditioners.jl

@@ -1,56 +1,4 @@
-## Diagonal Preconditioners
-
-struct DiagonalPreconditioner{D}
-    diag::D
-end
-struct InvDiagonalPreconditioner{D}
-    diag::D
-end
-
-Base.eltype(A::Union{DiagonalPreconditioner,InvDiagonalPreconditioner}) = eltype(A.diag)
-
-function LinearAlgebra.ldiv!(A::DiagonalPreconditioner, x)
-    x .= x ./ A.diag
-end
-
-function LinearAlgebra.ldiv!(y, A::DiagonalPreconditioner, x)
-    y .= x ./ A.diag
-end
-
-#=
-function LinearAlgebra.ldiv!(y::Matrix, A::DiagonalPreconditioner, b::Matrix)
-    @inbounds @simd for j ∈ 1:size(y, 2)
-        for i ∈ 1:length(A.diag)
-            y[i,j] = b[i,j] / A.diag[i]
-        end
-    end
-    return y
-end
-=#
-
-function LinearAlgebra.ldiv!(A::InvDiagonalPreconditioner, x)
-    x .= x .* A.diag
-end
-
-function LinearAlgebra.ldiv!(y, A::InvDiagonalPreconditioner, x)
-    y .= x .* A.diag
-end
-
-#=
-function LinearAlgebra.ldiv!(y::Matrix, A::InvDiagonalPreconditioner, b::Matrix)
-    @inbounds @simd for j ∈ 1:size(y, 2)
-        for i ∈ 1:length(A.diag)
-            y[i,j] = b[i,j] * A.diag[i]
-        end
-    end
-    return y
-end
-=#
-
-LinearAlgebra.mul!(y, A::DiagonalPreconditioner, x) = LinearAlgebra.ldiv!(y, InvDiagonalPreconditioner(A.diag), x)
-LinearAlgebra.mul!(y, A::InvDiagonalPreconditioner, x) = LinearAlgebra.ldiv!(y, DiagonalPreconditioner(A.diag), x)
-
-## Compose Preconditioner
+# Tooling Preconditioners
 
 struct ComposePreconditioner{Ti,To}
     inner::Ti
@@ -78,4 +26,6 @@ struct InvPreconditioner{T}
 end
 
 Base.eltype(A::InvPreconditioner) = Base.eltype(A.P)
+LinearAlgebra.ldiv!(A::InvPreconditioner, x) = mul!(x, A.P, x)
+LinearAlgebra.ldiv!(y, A::InvPreconditioner, x) = mul!(y, A.P, x)
 LinearAlgebra.mul!(y, A::InvPreconditioner, x) = ldiv!(y, A.P, x)

test/runtests.jl

@@ -201,27 +201,9 @@ end
 end
 
 @testset "Preconditioners" begin
-    @testset "Scalar Diagonal Preconditioner" begin
-        s = rand()
-
-        x = rand(n,n)
-        y = rand(n,n)
-
-        Pl, Pr = LinearSolve.DiagonalPreconditioner(s),LinearSolve.InvDiagonalPreconditioner(s)
-
-        mul!(y, Pl, x); @test y ≈ s * x
-        mul!(y, Pr, x); @test y ≈ s \ x
-
-        y .= x; ldiv!(Pl, x); @test x ≈ s \ y
-        y .= x; ldiv!(Pr, x); @test x ≈ s * y
-
-        ldiv!(y, Pl, x); @test y ≈ s \ x
-        ldiv!(y, Pr, x); @test y ≈ s * x
-    end
-
     @testset "Vector Diagonal Preconditioner" begin
         s = rand(n)
-        Pl, Pr = LinearSolve.DiagonalPreconditioner(s),LinearSolve.InvDiagonalPreconditioner(s)
+        Pl, Pr = Diagonal(s),LinearSolve.InvPreconditioner(Diagonal(s))
 
         x = rand(n,n)
         y = rand(n,n)
@@ -237,21 +219,20 @@ end
     end
 
     @testset "ComposePreconditioenr" begin
-        s1 = rand()
-        s2 = rand()
+        s1 = rand(n)
+        s2 = rand(n)
 
         x = rand(n,n)
         y = rand(n,n)
 
-        P1 = LinearSolve.DiagonalPreconditioner(s1)
-        P2 = LinearSolve.DiagonalPreconditioner(s2)
+        P1 = Diagonal(s1)
+        P2 = Diagonal(s2)
 
         P  = LinearSolve.ComposePreconditioner(P1,P2)
 
         # ComposePreconditioner
         ldiv!(y, P, x);      @test y ≈ ldiv!(P2, ldiv!(P1, x))
         y .= x; ldiv!(P, x); @test x ≈ ldiv!(P2, ldiv!(P1, y))
-
     end
 end