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

* start gmres draft

* add comments

* remove kdim from test

* use dot_product

* gmres: add Hilbert matrix test

* linalg: add MGS reorthogonalization to GMRES

* fix: rename loop variable to avoid fypp collision

* fix: hilbert sp tolerance test

* fix: hilbert sp tolerance test

* fix: hilbert dp tolerance test

* fix: hilbert dp tolerance test for Intel 2024 1 cmake

* Update test/linalg/test_linalg_solve_iterative.fypp

Co-authored-by: Jeremie Vandenplas <jeremie.vandenplas@gmail.com>

* Update test/linalg/test_linalg_solve_iterative.fypp

Co-authored-by: Jeremie Vandenplas <jeremie.vandenplas@gmail.com>

* Update test/linalg/test_linalg_solve_iterative.fypp

Co-authored-by: Jeremie Vandenplas <jeremie.vandenplas@gmail.com>

* Update test/linalg/test_linalg_solve_iterative.fypp

Co-authored-by: Jeremie Vandenplas <jeremie.vandenplas@gmail.com>

* Update test/linalg/test_linalg_solve_iterative.fypp

Co-authored-by: Jeremie Vandenplas <jeremie.vandenplas@gmail.com>

* fix declarations

* make gmres implementation flexible by allowing switching between memory efficiency or runtime speed

* add example

* convert static enum to integer helper function to determine the number of colums of the gmres workspace

* simplify iterations logic

* replace manual given rotations and the triangular solve by lapack kernels

* try removing the hilbert test for sp as it is meaningless for the selected size

* Refactor coverage report generation in workflow

Test proposal by copilot to fix coverage report issues. Testing in an actual PR exhibiting coverage issues because of the examples.

* Restrict fpm-deployment coverage upload to the filtered `src/` report (#13)

* Initial plan

* Disable Codecov auto-search in coverage job

---------

Co-authored-by: copilot-swe-agent[bot] <198982749+Copilot@users.noreply.github.com>

* add review corrections

* parametrize system size

---------

Co-authored-by: AlexanderGSC <alexander.gsc@gmail.com>
Co-authored-by: Jeremie Vandenplas <jeremie.vandenplas@gmail.com>
Co-authored-by: Copilot <198982749+Copilot@users.noreply.github.com>

Author: 4 people

.github/workflows/fpm-deployment.yml

@@ -114,20 +114,25 @@ jobs:
 
       - name: Create coverage report
         run: |
+          set -euo pipefail
           mkdir -p ${{ env.COV_DIR }}
           mv ./build/gfortran_*/*/* ${{ env.COV_DIR }}
           lcov --capture --initial --base-directory . --directory ${{ env.COV_DIR }} --output-file ${{ env.COV_DIR }}/coverage.base --ignore-errors mismatch
           lcov --capture           --base-directory . --directory ${{ env.COV_DIR }} --output-file ${{ env.COV_DIR }}/coverage.capture --ignore-errors mismatch
-          lcov --add-tracefile ${{ env.COV_DIR }}/coverage.base --add-tracefile ${{ env.COV_DIR }}/coverage.capture --output-file ${{ env.COV_DIR }}/coverage.info --ignore-errors mismatch
-          # Exclude example files from coverage: use both absolute-path pattern (*/example/*)
-          # and relative-path pattern (example/**) because lcov 2.x requires ** for recursive matching
-          # and the fpm build stores paths relative to the project root.
-          lcov --remove ${{ env.COV_DIR }}/coverage.info '*/example/*' 'example/**' --output-file ${{ env.COV_DIR }}/coverage.info --ignore-errors unused,mismatch,inconsistent
+          lcov --add-tracefile ${{ env.COV_DIR }}/coverage.base --add-tracefile ${{ env.COV_DIR }}/coverage.capture --output-file ${{ env.COV_DIR }}/coverage.full.info --ignore-errors mismatch
+          # Keep only src files in final coverage report to avoid example/ and test/ sources affecting the metric.
+          lcov --extract ${{ env.COV_DIR }}/coverage.full.info '*/src/*' 'src/**' --output-file ${{ env.COV_DIR }}/coverage.info --ignore-errors unused,mismatch,inconsistent
+          echo "Files included in final coverage report:"
+          grep '^SF:' ${{ env.COV_DIR }}/coverage.info | sed 's#^SF:##' | sort -u
         env:
             COV_DIR: build/coverage
 
       - name: Upload coverage report
         uses: codecov/codecov-action@v4
         with:
           token: ${{ secrets.CODECOV_TOKEN }}
+          # Upload only the lcov report rooted at src/; otherwise Codecov also
+          # discovers generated *.gcov files under example/ and treats them as
+          # primary analyzed sources.
+          disable_search: true
           files: build/coverage/coverage.info

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
+    integer, parameter :: n = 4
+    real(dp), parameter :: A(n,n) = 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], [n,n])
+    real(dp), parameter :: x0(n) = [0.2_dp, -0.1_dp, 0.3_dp, 0.0_dp]
+    real(dp), parameter :: xref(n) = [1.0_dp, 2.0_dp, -1.0_dp, 0.5_dp]
+    real(dp) :: rhs(n), x(n)
+
+    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
+    !------------------------------------------------------------------
+    elemental 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