Skip to content

Commit

Permalink
Merge pull request #63 from SciML/prec
Browse files Browse the repository at this point in the history
Change scaling_preconditioner to a proper set of DiagonalPreconditioner
  • Loading branch information
ChrisRackauckas authored Dec 16, 2021
2 parents 7b8bbbf + 2bb64f9 commit ad7c877
Show file tree
Hide file tree
Showing 3 changed files with 62 additions and 9 deletions.
1 change: 1 addition & 0 deletions docs/make.jl
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,7 @@ makedocs(
"Basics" => Any[
"basics/LinearProblem.md",
"basics/CachingAPI.md",
"basics/Preconditioners.md",
"basics/FAQ.md"
],
"Solvers" => Any[
Expand Down
58 changes: 55 additions & 3 deletions src/preconditioners.jl
Original file line number Diff line number Diff line change
@@ -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
Expand Down
12 changes: 6 additions & 6 deletions test/runtests.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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)
Expand All @@ -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)
Expand Down

0 comments on commit ad7c877

Please sign in to comment.