feat (linalg_iterative) : add gmres solver to the iterative methods library

[jalvesz]

This PR proposes adding the generalized minimal residual method (GMRES) to the linalg_iterative module.

It was drafted with the help of CoPilot and received the help of @AlexanderGSC to add the Modified Gram Schmidt scheme for increased robustness.

• Public core kernel implementation
• Standard interface for dense and CSR matrix types
• Unit test
• Examples
• Documentation

[codecov]

Codecov Report

❌ Patch coverage is 0% with 10 lines in your changes missing coverage. Please review.
✅ Project coverage is 68.79%. Comparing base (4c8521d) to head (0069cd3).

Files with missing lines	Patch %	Lines
example/linalg/example_solve_gmres.f90	0.00%	10 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #1186      +/-   ##
==========================================
- Coverage   68.81%   68.79%   -0.03%     
==========================================
  Files         408      409       +1     
  Lines       13726    13736      +10     
  Branches     1552     1552              
==========================================
+ Hits         9446     9450       +4     
- Misses       4280     4286       +6     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[AlexanderGSC]

Thanks @jalvesz!

If more tests or documentation need to be added, or if the maintainers think it’s necessary to make any changes, I’d be happy to help.

[jvdp1]

Thank you for this addition. Here are some suggestions.

[jvdp1]

is it on purpose the same as stdlib_size_wksp_cg?

[jalvesz]

It is not on purpose the as stdlib_size_wksp_cg. The workspace for the GMRES algorithm is bit more tricky than for the other gradient descent methods as not only are there reference buffer vectors with the size of the problem but also there are arrays depending on the kdim (number of internal iterations per restart cycle).

The current workspace design :
r(:) → residual (1 vector)
w(:) → Arnoldi work vector (1 vector)
v(:,1:kdim+1) → Krylov basis (kdim+1 vectors)
z(:,1:kdim) → preconditioned basis (kdim vectors)

Actually the last vector z could be reduced to z(:) as the preconditionner is applied per each j inner iteration.

I'll review this and try to propose a more lean version in terms of internal storage size. This won't change the fact that the workspace depends not only on n (problem size) but also on the selected number of internal iterations kdim.

Maybe @AlexanderGSC you have some suggestions here?

[AlexanderGSC]

Hi @jvdp1, @jalvesz!

@jalvesz is right; the size of the workspace depends on the size of the restart (kdim), and there isn’t much that can be done about it. Storing only the last vector from the preconditioner saves a lot of memory – well spotted!

I can think of a few ideas for keeping the interface unchanged, each with its own pros and cons:

• move the allocation of v inside the kernel:
  pros: keeps the current interface unchanged.
  cons: I don’t know if it’s allowed to reserve memory inside the kernel. At the moment there are variables such as the Hessemberg matrix, the rotation vectors, etc. that are allocated inside the kernel, but they are small, depending only on kdim, which by default is 30, meaning very little memory is used.

• Explicit workspace handling:
  Pros: no memory allocation within the kernel.
  Cons: we would need to document that the workspace size is 3+(kdim+1). And we would need to add a check in the code to ensure this is the case. Furthermore, in the interface, the value would not reflect the workspace size.

• Implicit kdim handling: the workspace size is checked, and the kdim size is deduced implicitly. If a value of 33 is set in the interface, gmres would run with a default kdim=30.
  Pros: allows running ‘by default’ or with ‘fine-tuning’ without allocating memory in the kernel.
  Cons: the way the workspace is managed becomes a little odd and less straightforward. Restart is a parameter that every version of gmres has.

In my opinion, if it isn’t a problem, I would move the v memory allocation into the kernel. If that is a problem, then we could explore other options.

[jalvesz]

I gave another round to the proposal. I added a compact option to allow switching between a memory efficient or a speed oriented implementation versions. Let me know your thoughts on this.

[AlexanderGSC]

Hello @jalvesz,

What you’re suggesting is much better; it’s a definitive solution that would allow for much greater flexibility and provide a clean interface for any method requiring a variable number of vectors.

Is the Hilbert test causing problems again in single precision? If so, I think the best thing is to remove it, or if necessary, I’ll modify it so that the test is skipped entirely in single precision. Adding the interface is much more important. There’ll be time to fix the test later on.

Great work!

[jalvesz]

Yes, the Hilbert test in single precision is causing problems with two of the intel jobs, don't know what is best, if just remove it or try with a jacobi preconditioner ?

[AlexanderGSC]

I don’t think using a preconditioner will make any difference.

The fact that it fails with Intel compilers suggests to me that they may be using aggressive instruction fusion optimisations that cause some numerical drift. The condition number of the Hilbert's matrix must be of the order of $10^{12}$, which is impossible for SP.

The simplest solution: reduce the size of the Hilbert matrix to something manageable for SP. With $n=4$ it should work.

If you agree, I can push the change to update the test to $n=4$ directly so we can see if the Intel CI jobs finally turn green.

[jalvesz]

Indeed! I'm thinking that another option would be to check the condition number estimate and pick a size such that the condition is below some critical threshold ? Using
$$\kappa \approx \frac{(1+\sqrt{2})^{4*N}}{\sqrt{N}}$$
a limit size could be precomputed at compile time

... or simply just do the test for kind >sp with fixed size.

[jvdp1]

should it be !> as the docstring precedes the declaration?

Suggested change

	!! linear operator matrix
	!> linear operator matrix

[jalvesz]

Oh yes it could. In that case all the docstring characters in this file should be changed.