Skip to content
This repository was archived by the owner on May 15, 2025. It is now read-only.

Trying to improve the Klement-solver #22

Merged
merged 7 commits into from
Jan 4, 2023
Merged

Conversation

Deltadahl
Copy link
Contributor

Trying to solve #20 and #19.
I tried the approach described in #18

But I have not used the singularity test @YingboMa explained:
"Or F = lu(J), and the singularity check is then abs(F.factors[end, end]) < singular_tol"
Although that solution was faster I was afraid that the LU-decomposition lu(J) could result in an error if J was singular.

Therefore I read about good approaches to checking for a singular matrix and found that

cond(J) > tolerance_value

was a good approach (https://stackoverflow.com/questions/13145948/how-to-find-out-if-a-matrix-is-singular and https://discourse.julialang.org/t/checking-for-singularity-of-matrix/54746/4)

Is this OK, or do you want me to change something?

@codecov
Copy link

codecov bot commented Jan 4, 2023

Codecov Report

Merging #22 (8915554) into main (e681591) will decrease coverage by 1.85%.
The diff coverage is 93.18%.

@@            Coverage Diff             @@
##             main      #22      +/-   ##
==========================================
- Coverage   91.32%   89.47%   -1.86%     
==========================================
  Files           8        8              
  Lines         219      247      +28     
==========================================
+ Hits          200      221      +21     
- Misses         19       26       +7     
Impacted Files Coverage Δ
src/klement.jl 90.56% <91.89%> (+0.24%) ⬆️
src/broyden.jl 92.59% <100.00%> (ø)
src/utils.jl 100.00% <100.00%> (ø)
src/falsi.jl 84.09% <0.00%> (-11.37%) ⬇️

📣 We’re building smart automated test selection to slash your CI/CD build times. Learn more

@ChrisRackauckas
Copy link
Member

Although that solution was faster I was afraid that the LU-decomposition lu(J) could result in an error if J was singular.

Just do lu(A, check = false).

  lu(A, pivot = RowMaximum(); check = true) -> F::LU

  Compute the LU factorization of A.

  When check = true, an error is thrown if the decomposition fails. When check = false, responsibility for checking the decomposition's validity (via issuccess) lies with the user.

src/klement.jl Outdated
@@ -16,7 +16,7 @@ function SciMLBase.solve(prob::NonlinearProblem,
x = float(prob.u0)
fₙ = f(x)
T = eltype(x)
J = ArrayInterfaceCore.zeromatrix(x) + I
J = init_J(x)
Copy link
Member

Choose a reason for hiding this comment

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

Do this to Broyden too

@Deltadahl
Copy link
Contributor Author

I didn't manage to get the program non-allocating on static-arrays due to @ballocated lu(J, check = false) > 0.
To get the program non-allocating for scalars, I split the program into 2 cases.

Just tell me if I can improve anything,
And thx for all the good feedback, I'm learning a lot from it! :)

src/klement.jl Outdated
Comment on lines 92 to 97
F = lu(J, check = false)

iszero(fₙ) &&
return SciMLBase.build_solution(prob, alg, xₙ, fₙ;
retcode = ReturnCode.Success)
# Singularity test
if any(abs.(F.U[diagind(F.U)]) .< singular_tol)
J = init_J(xₙ)
F = lu(J, check = false)
Copy link
Member

Choose a reason for hiding this comment

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

This does one more LU factorization than is needed. Instead, move the factorization to the top

src/klement.jl Outdated
return SciMLBase.build_solution(prob, alg, xₙ, fₙ;
retcode = ReturnCode.Success)
# Singularity test
if any(abs.(F.U[diagind(F.U)]) .< singular_tol)
Copy link
Member

@YingboMa YingboMa Jan 4, 2023

Choose a reason for hiding this comment

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

This allocates unnecessarily.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I didn't get abs(F.factors[end, end]) to work for static arrays (F.factors gives error)

Copy link
Contributor Author

Choose a reason for hiding this comment

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

And for non-static arrays abs(F.factors[end, end]) worked fine if J was non-singular.
But it didn't work as I excpected when J was singular, i.e. it gave abs(F.factors[end, end])=3.0 and abs(F.factors[end-1, end-1])=0.0, resulting in the solver ending up with NaN as output.

Copy link
Member

Choose a reason for hiding this comment

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

Let's open up an issue on this. It may need to be specialized more, but for now we can go with this.

Sign up for free to subscribe to this conversation on GitHub. Already have an account? Sign in.
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants