Skip to content

Bugfix: numerical_setup! for PETScLinearSolver #104

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
wants to merge 5 commits into
base: master
Choose a base branch
from
Draft

Bugfix: numerical_setup! for PETScLinearSolver #104

Show file tree
Hide file tree
Changes from 4 commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
88e42e3
Bugfix: numerical_setup! now reuses the PETSc matrix
JordiManyer Oct 1, 2024
8832a02
Merge branch 'master' of github.com:gridap/GridapPETSc.jl into bugfix…
JordiManyer Oct 1, 2024
b0c90fc
Relaxed checks
JordiManyer Oct 1, 2024
b23197d
Added copy! between PETScMatrices
JordiManyer Oct 1, 2024
c0e94ae
Bugfix: aliasing
JordiManyer Oct 2, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
10 changes: 7 additions & 3 deletions src/PETScArrays.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -306,6 +306,13 @@ function Base.copy(a::PETScMatrix)
Init(v)
end

function Base.copy!(a::PETScMatrix,b::PETScMatrix)
if a !== b
@check_error_code PETSC.MatCopy(b.mat[],a.mat[],PETSC.SAME_NONZERO_PATTERN)
end
a
end

function Base.copy!(a::PETScMatrix,b::AbstractMatrix)
_copy!(a.mat[],b)
end
Expand Down Expand Up @@ -363,9 +370,6 @@ function _copy!(petscmat::Mat,mat::AbstractSparseMatrix)
@check_error_code PETSC.MatAssemblyEnd(petscmat, PETSC.MAT_FINAL_ASSEMBLY)
end




function Base.convert(::Type{PETScMatrix},a::PETScMatrix)
a
end
Expand Down
20 changes: 17 additions & 3 deletions src/PETScLinearSolvers.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -98,10 +98,24 @@ function Algebra.solve!(x::PVector,ns::PETScLinearSolverNS,b::PVector)
x
end

# NOTE:
# We previously threw away PETSc's matrix `ns.B`, and re-set the KSP object (commented code).
# This is not only unnecessary, but can also cause some issues with MUMPS.
# I am not completely sure why, but here are some notes on the issue:
# - It is matrix-dependent, and only happens in parallel (nprocs > 1).
# - It has to do with the re-use of the symmetric permutation created by MUMPS to
# find the pivots.
# I think it probably re-orders the matrix internally, and does not re-order it again when we swap it
# using `KSPSetOperators`. So when accessing the new (non-permuted) matrix using the old permutation,
# it throws a stack overflow error.
# In fact, when updating the SNES setups in the nonlinear solvers, we do not re-set the matrix,
# but use `copy!` instead.
function Algebra.numerical_setup!(ns::PETScLinearSolverNS,A::AbstractMatrix)
ns.A = A
ns.B = convert(PETScMatrix,A)
@check_error_code PETSC.KSPSetOperators(ns.ksp[],ns.B.mat[],ns.B.mat[])
# ns.A = A
# ns.B = convert(PETScMatrix,A)
# @check_error_code PETSC.KSPSetOperators(ns.ksp[],ns.B.mat[],ns.B.mat[])
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

are we sure that we do no longer need to call KSPSetOperators?

Copy link
Member Author

@JordiManyer JordiManyer Oct 2, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is what the manual says:

To solve successive linear systems that have different preconditioner matrices 
(i.e., the matrix elements and/or the matrix data structure change), the user 
must call KSPSetOperators and KSPSolve for each solve.

which would indicate we do. However, my tests point to this not being the case. This works, for instance:

solver = PETScLinearSolver(mykspsetup)
println("First call")
ns = numerical_setup(symbolic_setup(solver,A),A)
x = allocate_in_domain(A)
fill!(x,0.0)
solve!(x,ns,b)

println("Second call")
map(Ai -> rmul!(Ai,2.0), partition(A))
b *= 2.0
numerical_setup!(ns,A)
x2 = allocate_in_domain(A)
fill!(x2,0.0)
solve!(x2,ns,b)

e = norm(x-x2)
println("Error =", e)

@assert nnz(ns.A) == nnz(A) # This is weak, but it might catch some errors
copy!(ns.B,A)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In this copy! call, which are the possible types for A?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In theory, all supported types

@check_error_code PETSC.KSPSetUp(ns.ksp[])
ns
end
Loading