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