Merge pull request #63 from SciML/prec

Change scaling_preconditioner to a proper set of DiagonalPreconditioner

Author: ChrisRackauckas

docs/make.jl

@@ -20,6 +20,7 @@ makedocs(
         "Basics" => Any[
             "basics/LinearProblem.md",
             "basics/CachingAPI.md",
+            "basics/Preconditioners.md",
             "basics/FAQ.md"
         ],
         "Solvers" => Any[

src/preconditioners.jl

@@ -1,7 +1,59 @@
-## Preconditioners
+## Diagonal Preconditioners
 
-scaling_preconditioner(s::Number) = I * (1/s), I * s
-scaling_preconditioner(s::AbstractVector) = Diagonal(inv.(s)),Diagonal(s)
+struct DiagonalPreconditioner{D}
+    diag::D
+end
+struct InvDiagonalPreconditioner{D}
+    diag::D
+end
+
+Base.eltype(A::Union{DiagonalPreconditioner,InvDiagonalPreconditioner}) = eltype(A.diag)
+Base.adjoint(A::Union{DiagonalPreconditioner,InvDiagonalPreconditioner}) = A
+Base.inv(A::DiagonalPreconditioner) = InvDiagonalPreconditioner(A.diag)
+Base.inv(A::InvDiagonalPreconditioner) = DiagonalPreconditioner(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
 
 struct ComposePreconditioner{Ti,To}
     inner::Ti

test/runtests.jl

@@ -197,13 +197,13 @@ end
 end
 
 @testset "Preconditioners" begin
-    @testset "scaling_preconditioner" begin
+    @testset "Scalar Diagonal Preconditioner" begin
         s = rand()
 
         x = rand(n,n)
         y = rand(n,n)
 
-        Pl, Pr = LinearSolve.scaling_preconditioner(1/s)
+        Pl, Pr = LinearSolve.DiagonalPreconditioner(s),LinearSolve.InvDiagonalPreconditioner(s)
 
         mul!(y, Pl, x); @test y ≈ s * x
         mul!(y, Pr, x); @test y ≈ s \ x
@@ -215,9 +215,9 @@ end
         ldiv!(y, Pr, x); @test y ≈ s * x
     end
 
-    @testset "vector scaling_preconditioner" begin
+    @testset "Vector Diagonal Preconditioner" begin
         s = rand(n)
-        Pl, Pr = LinearSolve.scaling_preconditioner(1 ./ s)
+        Pl, Pr = LinearSolve.DiagonalPreconditioner(s),LinearSolve.InvDiagonalPreconditioner(s)
 
         x = rand(n,n)
         y = rand(n,n)
@@ -239,8 +239,8 @@ end
         x = rand(n,n)
         y = rand(n,n)
 
-        P1, _ = LinearSolve.scaling_preconditioner(s1)
-        P2, _ = LinearSolve.scaling_preconditioner(s2)
+        P1 = LinearSolve.DiagonalPreconditioner(s1)
+        P2 = LinearSolve.DiagonalPreconditioner(s2)
 
         P  = LinearSolve.ComposePreconditioner(P1,P2)
         Pi = LinearSolve.InvComposePreconditioner(P)