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

doc/specs/stdlib_linalg_iterative_solvers.md

@@ -348,4 +348,102 @@ The method uses 8 auxiliary vectors internally, requiring more memory than simpl
 #### Example 2
 ```fortran
 {!example/linalg/example_solve_bicgstab_wilkinson.f90!}
+```
+
+<!-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -->
+### `stdlib_solve_gmres_kernel` subroutine
+
+#### Description
+
+Implements the restarted Generalized Minimal Residual (GMRES) method for solving the linear system \( Ax = b \), where \( A \) is a general linear operator defined via the `stdlib_linop` type. This is the core implementation, allowing custom matrix types, custom preconditioners, and user-managed workspace.
+
+#### Syntax
+
+`call ` [[stdlib_linalg_iterative_solvers(module):stdlib_solve_gmres_kernel(interface)]] ` (A, M, b, x, rtol, atol, maxiter, kdim, workspace, compact)`
+
+#### Status
+
+Experimental
+
+#### Class
+
+Subroutine
+
+#### Argument(s)
+
+`A`: `class(stdlib_linop_<kind>_type)` defining the linear operator. This argument is `intent(in)`.
+
+`M`: `class(stdlib_linop_<kind>_type)` defining the preconditioner linear operator. This argument is `intent(in)`.
+
+`b`: 1-D array of `real(<kind>)` defining the loading conditions of the linear system. This argument is `intent(in)`.
+
+`x`: 1-D array of `real(<kind>)` which serves as the input initial guess and the output solution. This argument is `intent(inout)`.
+
+`rtol` and `atol`: scalars of type `real(<kind>)` specifying the convergence test. For convergence, the following criterion is used \( || b - Ax ||^2 <= max(rtol^2 * || b ||^2 , atol^2 ) \). These arguments are `intent(in)`.
+
+`maxiter`: scalar of type `integer` defining the maximum allowed number of iterations. This argument is `intent(in)`.
+
+`kdim`: scalar of type `integer` defining the Krylov subspace size before restart. This argument is `intent(in)`.
+
+`workspace`: scalar derived type of `type(stdlib_solver_workspace_<kind>_type)` holding the work array for the solver. This argument is `intent(inout)`.
+
+`compact`: scalar of type `logical`. If true, GMRES stores only one preconditioned basis vector at a time and rebuilds the cycle-end correction through the preconditioner, reducing memory usage. If false, GMRES caches the full preconditioned basis for a faster cycle-end update at the cost of additional memory. Default is `.true.`. This argument is `intent(in)`.
+
+#### Note
+
+GMRES trades increasing memory and orthogonalization cost for robust convergence on general non-symmetric systems. The `compact` option allows choosing between a lower-memory layout and a slightly faster restart update.
+
+<!-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -->
+### `stdlib_solve_gmres` subroutine
+
+#### Description
+
+Provides a user-friendly interface to the restarted GMRES method for solving \( Ax = b \), supporting `dense` and `CSR_<kind>_type` matrices. GMRES is suitable for general non-symmetric systems and supports optional preconditioning, workspace reuse, and a storage-mode choice for balancing memory use against cycle-end update speed.
+
+#### Syntax
+
+`call ` [[stdlib_linalg_iterative_solvers(module):stdlib_solve_gmres(interface)]] ` (A, b, x [, di, rtol, atol, maxiter, restart, kdim, precond, M, workspace, compact])`
+
+#### Status
+
+Experimental
+
+#### Class
+
+Subroutine
+
+#### Argument(s)
+
+`A`: `dense` or `CSR_<kind>_type` matrix defining the linear system. This argument is `intent(in)`.
+
+`b`: 1-D array of `real(<kind>)` defining the loading conditions of the linear system. This argument is `intent(in)`.
+
+`x`: 1-D array of `real(<kind>)` which serves as the input initial guess and the output solution. This argument is `intent(inout)`.
+
+`di` (optional): 1-D mask array of type `logical(int8)` defining the degrees of freedom subject to dirichlet boundary conditions. The actual boundary condition values should be stored in the `b` load array. This argument is `intent(in)`.
+
+`rtol` and `atol` (optional): scalars of type `real(<kind>)` specifying the convergence test. For convergence, the following criterion is used \( || b - Ax ||^2 <= max(rtol^2 * || b ||^2 , atol^2 ) \). Default values are `rtol=1.e-5` and `atol=epsilon(1._<kind>)`. These arguments are `intent(in)`.
+
+`maxiter` (optional): scalar of type `integer` defining the maximum allowed number of iterations. If no value is given, a default of `N` is set, where `N = size(b)`. This argument is `intent(in)`.
+
+`restart` (optional): scalar of type `logical` indicating whether to restart the iteration with zero initial guess. Default is `.true.`. This argument is `intent(in)`.
+
+`kdim` (optional): scalar of type `integer` defining the Krylov subspace size before restart. If no value is given, a default of `min(30, N)` is used. This argument is `intent(in)`.
+
+`precond` (optional): scalar of type `integer` enabling to switch among the default preconditioners available with the following enum (`pc_none`, `pc_jacobi`). If no value is given, no preconditioning will be applied. This argument is `intent(in)`.
+
+`M` (optional): scalar derived type of `class(stdlib_linop_<kind>_type)` defining a custom preconditioner linear operator. If given, `precond` will have no effect, and a pointer is set to this custom preconditioner. This argument is `intent(in)`.
+
+`workspace` (optional): scalar derived type of `type(stdlib_solver_workspace_<kind>_type)` holding the work temporal array for the solver. If the user passes its own `workspace`, then internally a pointer is set to it, otherwise, memory will be internally allocated and deallocated before exiting the procedure. This argument is `intent(inout)`.
+
+`compact` (optional): scalar of type `logical`. If true, GMRES stores only one preconditioned basis vector per iteration cycle and minimizes memory use. If false, GMRES caches the full preconditioned basis to speed up the cycle-end update. Default is `.true.`. This argument is `intent(in)`.
+
+#### Note
+
+GMRES is especially useful for general non-symmetric systems where CG and PCG do not apply. Its main cost is the Krylov basis storage and orthogonalization work within each restart cycle. With `compact=.true.`, the solver reduces workspace requirements significantly. With `compact=.false.`, it retains the previous cached-basis update strategy.
+
+#### Example
+
+```fortran
+{!example/linalg/example_solve_gmres.f90!}
 ```

example/linalg/CMakeLists.txt

@@ -48,6 +48,7 @@ if (STDLIB_LINALG_ITERATIVE)
   ADD_EXAMPLE(solve_bicgstab)
   ADD_EXAMPLE(solve_bicgstab_wilkinson)
   ADD_EXAMPLE(solve_cg)
+  ADD_EXAMPLE(solve_gmres)
   ADD_EXAMPLE(solve_pcg)
   ADD_EXAMPLE(solve_custom)
 endif()

example/linalg/example_solve_gmres.f90

@@ -0,0 +1,24 @@
+program example_solve_gmres
+    use stdlib_kinds, only: dp
+    use stdlib_linalg_iterative_solvers, only: stdlib_solve_gmres, pc_jacobi
+    implicit none
+
+    real(dp), parameter :: A(4,4) = reshape([5.0_dp,  1.0_dp,  2.0_dp, 0.0_dp, &
+                                             0.0_dp,  4.0_dp, -1.0_dp, 1.0_dp, &
+                                             1.0_dp,  0.0_dp,  3.0_dp, 2.0_dp, &
+                                             2.0_dp, -1.0_dp,  0.0_dp, 6.0_dp], [4,4])
+    real(dp), parameter :: x0(4) = [0.2_dp, -0.1_dp, 0.3_dp, 0.0_dp]
+    real(dp), parameter :: xref(4) = [1.0_dp, 2.0_dp, -1.0_dp, 0.5_dp]
+    real(dp) :: rhs(4), x(4)
+
+    rhs = matmul(A, xref)
+
+    ! GMRES lets the user tune the restart size and the storage mode.
+    x = x0
+    call stdlib_solve_gmres(A, rhs, x, restart=.false., kdim=2, maxiter=12, precond=pc_jacobi)
+    print *, 'compact:', x
+
+    x = x0
+    call stdlib_solve_gmres(A, rhs, x, restart=.false., kdim=2, maxiter=12, precond=pc_jacobi, compact=.false.)
+    print *, 'cached :', x
+end program example_solve_gmres

src/linalg_iterative/CMakeLists.txt

@@ -1,6 +1,7 @@
 set(linalg_iterative_fppFiles
     stdlib_linalg_iterative_solvers_bicgstab.fypp
     stdlib_linalg_iterative_solvers_cg.fypp
+    stdlib_linalg_iterative_solvers_gmres.fypp
     stdlib_linalg_iterative_solvers.fypp
     stdlib_linalg_iterative_solvers_pcg.fypp
     )
@@ -13,4 +14,8 @@ set(linalg_iterative_f90Files
 
 configure_stdlib_target(${PROJECT_NAME}_linalg_iterative linalg_iterative_f90Files linalg_iterative_fppFiles linalg_iterative_cppFiles)
 
-target_link_libraries(${PROJECT_NAME}_linalg_iterative PUBLIC ${PROJECT_NAME}_linalg ${PROJECT_NAME}_sparse)
+target_link_libraries(${PROJECT_NAME}_linalg_iterative 
+                        PUBLIC 
+                        ${PROJECT_NAME}_blas
+                        ${PROJECT_NAME}_linalg 
+                        ${PROJECT_NAME}_sparse)

src/linalg_iterative/stdlib_linalg_iterative_solvers.fypp

@@ -17,7 +17,7 @@ module stdlib_linalg_iterative_solvers
         enumerator :: stdlib_size_wksp_pcg = 4
         enumerator :: stdlib_size_wksp_bicgstab = 8
     end enum
-    public :: stdlib_size_wksp_cg, stdlib_size_wksp_pcg, stdlib_size_wksp_bicgstab
+    public :: stdlib_size_wksp_cg, stdlib_size_wksp_pcg, stdlib_size_wksp_bicgstab, stdlib_size_wksp_gmres
 
     enum, bind(c)
         enumerator :: pc_none = 0
@@ -161,6 +161,28 @@ module stdlib_linalg_iterative_solvers
     end interface
     public :: stdlib_solve_bicgstab_kernel
 
+    !! version: experimental
+    !!
+    !! stdlib_solve_gmres_kernel interface for the restarted generalized minimal residual method.
+    !! [Specifications](../page/specs/stdlib_linalg_iterative_solvers.html#stdlib_solve_gmres_kernel)
+    interface stdlib_solve_gmres_kernel
+        #:for k, t, s in R_KINDS_TYPES
+        module subroutine stdlib_solve_gmres_kernel_${s}$(A,M,b,x,rtol,atol,maxiter,kdim,workspace,compact)
+            class(stdlib_linop_${s}$_type), intent(in) :: A !! linear operator
+            class(stdlib_linop_${s}$_type), intent(in) :: M !! preconditioner linear operator
+            ${t}$, intent(in) :: b(:) !! right-hand side vector
+            ${t}$, intent(inout) :: x(:) !! solution vector and initial guess
+            ${t}$, intent(in) :: rtol !! relative tolerance for convergence
+            ${t}$, intent(in) :: atol !! absolute tolerance for convergence
+            integer, intent(in) :: maxiter !! maximum number of iterations
+            integer, intent(in) :: kdim !! Krylov subspace size before restart
+            type(stdlib_solver_workspace_${s}$_type), intent(inout) :: workspace !! workspace for the solver
+            logical, intent(in) :: compact !! keep only one preconditioned basis vector and rebuild the cycle update at restart
+        end subroutine
+        #:endfor
+    end interface
+    public :: stdlib_solve_gmres_kernel
+
     !! version: experimental
     !!
     !! [Specifications](../page/specs/stdlib_linalg_iterative_solvers.html#stdlib_solve_pcg)
@@ -219,7 +241,47 @@ module stdlib_linalg_iterative_solvers
     end interface
     public :: stdlib_solve_bicgstab
 
+    !! version: experimental
+    !!
+    !! [Specifications](../page/specs/stdlib_linalg_iterative_solvers.html#stdlib_solve_gmres)
+    interface stdlib_solve_gmres
+        #:for matrix in MATRIX_TYPES
+        #:for k, t, s in R_KINDS_TYPES
+        module subroutine stdlib_solve_gmres_${matrix}$_${s}$(A,b,x,di,rtol,atol,maxiter,restart,kdim,precond,M,workspace,compact)
+            !! linear operator matrix
+            #:if matrix == "dense"
+            ${t}$, intent(in) :: A(:,:)
+            #:else
+            type(${matrix}$_${s}$_type), intent(in) :: A
+            #:endif
+            ${t}$, intent(in) :: b(:) !! right-hand side vector
+            ${t}$, intent(inout) :: x(:) !! solution vector and initial guess
+            ${t}$, intent(in), optional :: rtol !! relative tolerance for convergence
+            ${t}$, intent(in), optional :: atol !! absolute tolerance for convergence
+            logical(int8), intent(in), optional, target :: di(:) !! dirichlet conditions mask
+            integer, intent(in), optional :: maxiter !! maximum number of iterations
+            logical, intent(in), optional :: restart !! restart flag
+            integer, intent(in), optional :: kdim !! Krylov subspace size before restart
+            integer, intent(in), optional :: precond !! preconditioner method enumerator
+            class(stdlib_linop_${s}$_type), optional, intent(in), target :: M !! preconditioner linear operator
+            type(stdlib_solver_workspace_${s}$_type), optional, intent(inout), target :: workspace !! workspace for the solver
+            logical, intent(in), optional :: compact !! if true, favor lower memory use over cycle-end update speed
+        end subroutine
+        #:endfor
+        #:endfor
+    end interface
+    public :: stdlib_solve_gmres
+
 contains
+
+    !------------------------------------------------------------------
+    ! helpers
+    !------------------------------------------------------------------
+    integer(int32) function stdlib_size_wksp_gmres(kdim,compact)
+        integer(int32), intent(in) :: kdim
+        logical, intent(in) :: compact
+        stdlib_size_wksp_gmres = 3 + kdim + merge(1, kdim, compact)
+    end function
 
     !------------------------------------------------------------------
     ! defaults