Skip to content

Commit

Permalink
Handle the "maybe extend" behaviour of extendMatrix() more
Browse files Browse the repository at this point in the history
appropriately.

Currently, if the given SDDM has a very small diagonal, the extension
does not happen, but the RHS is augmented anyway. This causes a crash.
  • Loading branch information
pratyai committed Feb 14, 2024
1 parent c18861f commit 530d602
Showing 1 changed file with 17 additions and 11 deletions.
28 changes: 17 additions & 11 deletions src/solverInterface.jl
Original file line number Diff line number Diff line change
Expand Up @@ -323,25 +323,31 @@ Uses a `lapSolver` to solve systems of linear equations in sddm matrices.
"""
function sddmWrapLap(lapSolver, sddm::AbstractArray; tol::Real=1e-6, maxits=Inf, maxtime=Inf, verbose=false, pcgIts=Int[], params...)

# Make a new adj matrix, a1, with an extra entry at the end.
a, dVal, dExcess = adjValAndExcess(sddm)
a1 = extendMatrix(a,dVal, dExcess)
F = lapSolver(a1; tol=tol, maxits=maxits, maxtime=maxtime, verbose=verbose, pcgIts=pcgIts, params...)

# make a function that solves the extended system, modulo the last entry
tol_=tol
maxits_=maxits
maxtime_=maxtime
verbose_=verbose
pcgIts_=pcgIts

f = function(b; tol=tol_, maxits=maxits_, maxtime=maxtime_, verbose=verbose_, pcgIts=pcgIts_)
xaug = F([b; -sum(b)], tol=tol, maxits=maxits, maxtime=maxtime, verbose=verbose, pcgIts=pcgIts)
xaug = xaug .- xaug[end]
return xaug[1:a.n]
end
# Make a new adj matrix, a1, with an extra entry at the end.
a, dVal, dExcess = adjValAndExcess(sddm)
a1 = extendMatrix(a,dVal, dExcess)

return f
# `extendMatrix()` may or may not have actually made an extension.
if size(a1) == size(sddm)
# If it didn't, then just treat the SDDM as a regular Laplacian.
local adjm = adj(sddm)[1]
return lapSolver(adjm; tol=tol, maxits=maxits, maxtime=maxtime, verbose=verbose, pcgIts=pcgIts, params...)
else
F = lapSolver(a1; tol=tol, maxits=maxits, maxtime=maxtime, verbose=verbose, pcgIts=pcgIts, params...)
f = function(b; tol=tol_, maxits=maxits_, maxtime=maxtime_, verbose=verbose_, pcgIts=pcgIts_)
xaug = F([b; -sum(b)], tol=tol, maxits=maxits, maxtime=maxtime, verbose=verbose, pcgIts=pcgIts)
xaug = xaug .- xaug[end]
return xaug[1:a.n]
end
return f
end

end

Expand Down

0 comments on commit 530d602

Please sign in to comment.