From 88e42e3cfda0ffce1ad18fff7d4200e463f49b99 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Tue, 1 Oct 2024 16:55:45 +1000 Subject: [PATCH 1/4] Bugfix: numerical_setup! now reuses the PETSc matrix --- src/PETScLinearSolvers.jl | 24 +++++++++++++++++++++--- 1 file changed, 21 insertions(+), 3 deletions(-) diff --git a/src/PETScLinearSolvers.jl b/src/PETScLinearSolvers.jl index 291ad40..f7627c6 100644 --- a/src/PETScLinearSolvers.jl +++ b/src/PETScLinearSolvers.jl @@ -98,10 +98,28 @@ 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. +# +# TODO: +# I've added a check to ensure the input matrix `A` is always teh initial one. +# Is this strictly necessary? Could we check the sparsity pattern 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[]) + @assert ns.A === A + copy!(ns.B,A) @check_error_code PETSC.KSPSetUp(ns.ksp[]) ns end From b0c90fc763913cd2d506db9147d60412d2ce769b Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Wed, 2 Oct 2024 09:00:54 +1000 Subject: [PATCH 2/4] Relaxed checks --- src/PETScLinearSolvers.jl | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/src/PETScLinearSolvers.jl b/src/PETScLinearSolvers.jl index f7627c6..6596425 100644 --- a/src/PETScLinearSolvers.jl +++ b/src/PETScLinearSolvers.jl @@ -110,15 +110,11 @@ end # 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. -# -# TODO: -# I've added a check to ensure the input matrix `A` is always teh initial one. -# Is this strictly necessary? Could we check the sparsity pattern 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[]) - @assert ns.A === A + @assert nnz(ns.A) == nnz(A) # This is weak, but it might catch some errors copy!(ns.B,A) @check_error_code PETSC.KSPSetUp(ns.ksp[]) ns From b23197df67030cb7bb5c8d5caa9ff8ff43c67143 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Wed, 2 Oct 2024 09:48:16 +1000 Subject: [PATCH 3/4] Added copy! between PETScMatrices --- src/PETScArrays.jl | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/src/PETScArrays.jl b/src/PETScArrays.jl index 1f455ba..6fb6e53 100644 --- a/src/PETScArrays.jl +++ b/src/PETScArrays.jl @@ -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 @@ -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 From c0e94aecd201f0c844fa7f64efa87b0357649555 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Wed, 2 Oct 2024 12:58:41 +1000 Subject: [PATCH 4/4] Bugfix: aliasing --- src/PETScLinearSolvers.jl | 27 ++++++++++++++++++++------- 1 file changed, 20 insertions(+), 7 deletions(-) diff --git a/src/PETScLinearSolvers.jl b/src/PETScLinearSolvers.jl index 6596425..76e18e3 100644 --- a/src/PETScLinearSolvers.jl +++ b/src/PETScLinearSolvers.jl @@ -110,12 +110,25 @@ end # 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[]) +# @assert nnz(ns.A) == nnz(A) # This is weak, but it might catch some errors +# copy!(ns.B,A) +# @check_error_code PETSC.KSPSetUp(ns.ksp[]) +# ns +#end + 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[]) - @assert nnz(ns.A) == nnz(A) # This is weak, but it might catch some errors - copy!(ns.B,A) - @check_error_code PETSC.KSPSetUp(ns.ksp[]) + if ns.A === A + copy!(ns.B,A) + @check_error_code PETSC.KSPSetUp(ns.ksp[]) + else + ns.A = A + ns.B = convert(PETScMatrix,A) + @check_error_code PETSC.KSPSetOperators(ns.ksp[],ns.B.mat[],ns.B.mat[]) + @check_error_code PETSC.KSPSetUp(ns.ksp[]) + end ns -end +end \ No newline at end of file