feat: [linalg] add iterative solvers #994
Conversation
…r workspace sizes
|
Thank you @jalvesz . Nice approach. |
jalvesz
changed the title
There was a problem hiding this comment.
Thank you @jalvesz . Interesting approach that is very flexible. In most of my practical cases, the coefficient matrix cannot be computed explicitely, even if very sparse. Your DT approach will still allow me to use thses solvers.
I suggest to split this PR in two: 1 for the iterative solvers (ldlt, ssor,...) and 1 for the CG solvers.
Regarding the preconditioners, I think that proposing only identity/none and diagonal is good enough for a first PR. Additional preconditioners can be provided in following PRs.
I think that such an approach will help to review this PR.
| #:for k, t, s in R_KINDS_TYPES | ||
| type, public :: solver_workspace_${s}$ | ||
| ${t}$, allocatable :: tmp(:,:) | ||
| procedure(logger_sub_${s}$), pointer, nopass :: callback => null() |
There was a problem hiding this comment.
Could the default logger the logger provided by stdlib?
There was a problem hiding this comment.
Sure, I just wasn't sure what would be a good default logger, do you have a proposal?
There was a problem hiding this comment.
Actually, looking again to the code, I think that your approach will provide full flexibility to the user.
There was a problem hiding this comment.
@Simkern worked quite extensively on the logging system for LightKrylov, may be he has some suggestions. Note however that our logger, based on the stdlib one works for the whole package and reports way more info that what would be needed for the CG solver.
|
|
||
| `x`: 1-D array of `real(<kind>)` which serves as the input initial guess and the output solution. This argument is `intent(inout)`. | ||
|
|
||
| `tol`: scalar of type `real(<kind>)` specifying the requested tolerance with respect to the relative residual vector norm. This argument is `intent(in)`. |
There was a problem hiding this comment.
tolerance with respect to the relative residual vector norm:
Maybe for a follow-up PR: should we support different convergence criteria?
There was a problem hiding this comment.
Which criteria would you consider adding?
I guess a monitoring on the solution vector instead of the residual could also be considered. something like
There was a problem hiding this comment.
Yes, indeed. Here are a couple of others: https://doi.org/10.1186/s12711-021-00626-1 , related to the relative error in x.
Co-authored-by: Jeremie Vandenplas <[email protected]>
| !! The `linop` type is used to define the linear operator for the iterative solvers. | ||
| #:for k, t, s in R_KINDS_TYPES | ||
| type, public :: linop_${s}$_type | ||
| procedure(vector_sub_${s}$), nopass, pointer :: apply => null() |
There was a problem hiding this comment.
Thinking about other methods to add, it might be necessary to cover also applying an in-place conjugate transpose. Two ideas to enable this:
- the apply procedure takes as input also an "operation"
character(1)argument > ('N','T','H') - a secondary procedure pointer is added such as
apply_transpose
?
thought, there are variants that do not need the transpose, so not necessarily a must. For instance, instead of implementing BiCG, the BiCG-Stab method doesn't require applying the transpose. GMRES doesn't need it either...
There was a problem hiding this comment.
The choice we made in LightKrylov is to have two type-bound procedures, one matvec for rmatvec for scipy to define a LinearOperator.
Although CG, CMRES, BiCG-Stab and a handful of others do no use the transpose, CGNE and CGNR (conjugate gradient for the normal equations when solving large-scale linear least-squares problems) actually do. While maybe not necessary for this first PR, it may be a good thing to keep this in mind for future extensions.
There was a problem hiding this comment.
It might then be better to add the procedure in this PR even if keeping it as null pointer to start with.
Should we use the naming convention of SciPy as you did for LightKrylov (matvec+rmatvec)? or would it be more interesting to think about a different design?
There was a problem hiding this comment.
in e5fa984 I propose to keep the use of the op input argument to indicate No transpose, Transpose or Hermitian. That way there is just one procedure signature consistent with gemv and spmv.
There was a problem hiding this comment.
Pull Request Overview
This PR introduces iterative solvers to the linear algebra module, implementing the Conjugate Gradient (CG) and Preconditioned Conjugate Gradient (PCG) methods. The implementation provides both low-level kernel interfaces for maximum flexibility and high-level APIs for convenient use with dense and CSR matrices.
- Adds
stdlib_linalg_iterative_solversmodule with core types and interfaces - Implements CG and PCG solver algorithms in separate submodules
- Provides comprehensive test coverage and example programs
Reviewed Changes
Copilot reviewed 12 out of 12 changed files in this pull request and generated no comments.
Show a summary per file
| File | Description |
|---|---|
| src/stdlib_linalg_iterative_solvers.fypp | Main module defining public interfaces, types, and abstract interfaces for iterative solvers |
| src/stdlib_linalg_iterative_solvers_cg.fypp | Submodule implementing the Conjugate Gradient solver algorithm |
| src/stdlib_linalg_iterative_solvers_pcg.fypp | Submodule implementing the Preconditioned Conjugate Gradient solver algorithm |
| src/stdlib_sparse_constants.fypp | Removes duplicate constant definitions, now using stdlib_constants |
| test/linalg/test_linalg_solve_iterative.fypp | Test suite for both CG and PCG solvers |
| example/linalg/example_solve_cg.f90 | Simple example demonstrating CG solver usage |
| example/linalg/example_solve_pcg.f90 | Example showing PCG solver with Dirichlet boundary conditions |
| example/linalg/example_solve_custom.f90 | Advanced example showing custom solver implementation |
| doc/specs/stdlib_linalg_iterative_solvers.md | Comprehensive documentation for the new iterative solvers |
Tip: Customize your code reviews with copilot-instructions.md. Create the file or learn how to get started.
This daft PR introduces the base for iterative solvers.
Two methods are proposed:
Each method is made public with two public interface flavors:
solve_<method>_kernel: All arguments are mandatory (no optionals/no internal allocations), it contains the methods steps. The linear system (and preconditioner) is defined through a public DTlinopwhich enables to extend two key procedures:applyequivalent to a matrix-vector product andinnerequivalent to the dot_product. This is the interface recommended to extend the method when dealing with a matrix type not available instdlibor when working in distributed-memory frameworks for which the matvec, dot_product and factorization steps need to be adapted to account for parallel synchronization.solve_<method>: API with optional arguments, the linear system can be defined with dense orCSR_<>_typematrices. For the PCG, the following preconditioners are available:none(identity),jacobi(1/diagonal). It internally usessolve_<method>_kernel.