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

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 68.84%. Comparing base (4c8521d) to head (825fd51).

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #1186      +/-   ##
==========================================
+ Coverage   68.81%   68.84%   +0.02%     
==========================================
  Files         408      408              
  Lines       13726    13726              
  Branches     1552     1552              
==========================================
+ Hits         9446     9450       +4     
+ Misses       4280     4276       -4     

☔ View full report in Codecov by Sentry.
📢 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.

[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.