diff --git a/doc/specs/stdlib_linalg_iterative_solvers.md b/doc/specs/stdlib_linalg_iterative_solvers.md index 6cebb61aa..f2c9173b9 100644 --- a/doc/specs/stdlib_linalg_iterative_solvers.md +++ b/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__type)` defining the linear operator. This argument is `intent(in)`. + +`M`: `class(stdlib_linop__type)` defining the preconditioner linear operator. This argument is `intent(in)`. + +`b`: 1-D array of `real()` defining the loading conditions of the linear system. This argument is `intent(in)`. + +`x`: 1-D array of `real()` which serves as the input initial guess and the output solution. This argument is `intent(inout)`. + +`rtol` and `atol`: scalars of type `real()` 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__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__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__type` matrix defining the linear system. This argument is `intent(in)`. + +`b`: 1-D array of `real()` defining the loading conditions of the linear system. This argument is `intent(in)`. + +`x`: 1-D array of `real()` 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()` 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._)`. 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__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__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!} ``` \ No newline at end of file diff --git a/example/linalg/CMakeLists.txt b/example/linalg/CMakeLists.txt index dd65b48e4..ae2e0f6f3 100644 --- a/example/linalg/CMakeLists.txt +++ b/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() diff --git a/example/linalg/example_solve_gmres.f90 b/example/linalg/example_solve_gmres.f90 new file mode 100644 index 000000000..02ea8ad12 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/src/linalg_iterative/CMakeLists.txt b/src/linalg_iterative/CMakeLists.txt index c91188eba..adea11ca4 100644 --- a/src/linalg_iterative/CMakeLists.txt +++ b/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) diff --git a/src/linalg_iterative/stdlib_linalg_iterative_solvers.fypp b/src/linalg_iterative/stdlib_linalg_iterative_solvers.fypp index 6c7725f3f..9bc1b6f8f 100644 --- a/src/linalg_iterative/stdlib_linalg_iterative_solvers.fypp +++ b/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 diff --git a/src/linalg_iterative/stdlib_linalg_iterative_solvers_gmres.fypp b/src/linalg_iterative/stdlib_linalg_iterative_solvers_gmres.fypp new file mode 100644 index 000000000..31c4ce79e --- /dev/null +++ b/src/linalg_iterative/stdlib_linalg_iterative_solvers_gmres.fypp @@ -0,0 +1,304 @@ +#:include "common.fypp" +#:set R_KINDS_TYPES = list(zip(REAL_KINDS, REAL_TYPES, REAL_SUFFIX)) +#:set C_KINDS_TYPES = list(zip(CMPLX_KINDS, CMPLX_TYPES, CMPLX_SUFFIX)) +#:set MATRIX_TYPES = ["dense", "CSR"] +#:set RANKS = range(1, 2+1) + +submodule(stdlib_linalg_iterative_solvers) stdlib_linalg_iterative_gmres + use stdlib_kinds + use stdlib_sparse + use stdlib_constants + use stdlib_optval, only: optval + use stdlib_linalg_blas, only: gemv + implicit none + +contains + + #: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 + class(stdlib_linop_${s}$_type), intent(in) :: M + ${t}$, intent(in) :: b(:), rtol, atol + ${t}$, intent(inout) :: x(:) + integer, intent(in) :: maxiter, kdim + type(stdlib_solver_workspace_${s}$_type), intent(inout) :: workspace + logical, intent(in) :: compact + integer :: i, iter, j, j_final, iorth, jz, zbase + ${t}$ :: beta, denom, hnext, htmp, norm_sq, norm_sq0, temp, tolsq + ${t}$, allocatable :: cs(:), g(:), h(:,:), sn(:), y(:) + + allocate(h(kdim+1, kdim), cs(kdim), sn(kdim), g(kdim+1), y(kdim) ) + zbase = kdim + 4 + + associate( r => workspace%tmp(:,1), & + w => workspace%tmp(:,2), & + v => workspace%tmp(:,3:kdim+3), & + z => workspace%tmp(:,zbase:zbase+merge(0,kdim-1,compact)) ) + + iter = 0 + ! Initialize convergence targets from the right-hand side norm. + norm_sq0 = A%inner_product(b, b) + tolsq = max(rtol*rtol*norm_sq0, atol*atol) + + ! Form the initial residual and report the starting iterate. + r = b + call A%matvec(x, r, alpha=-one_${s}$, beta=one_${s}$, op='N') + norm_sq = A%inner_product(r, r) + if (associated(workspace%callback)) call workspace%callback(x, norm_sq, iter) + + if (norm_sq <= tolsq) return + + do while (iter < maxiter .and. norm_sq >= tolsq) + iter = iter + 1 + ! Start a new GMRES cycle from the current residual. + beta = sqrt(max(norm_sq, zero_${s}$)) + if (beta <= epsilon(one_${s}$)) exit + + h = zero_${s}$ + cs = zero_${s}$ + sn = zero_${s}$ + g = zero_${s}$ + y = zero_${s}$ + + ! Initialize the Krylov basis and least-squares right-hand side. + v(:,1) = r / beta + g(1) = beta + + j_final = 0 + do j = 1, kdim + ! Run Arnoldi with the preconditioned basis vector. + jz = merge(1, j, compact) + call M%matvec(v(:,j), z(:,jz), alpha=one_${s}$, beta=zero_${s}$, op='N') + call A%matvec(z(:,jz), w, alpha=one_${s}$, beta=zero_${s}$, op='N') + + ! Modified Gram Schmidt (MGSR) + do iorth = 1, 2 ! reorthogonalization + do i = 1, j + htmp = A%inner_product(v(:,i), w) + h(i,j) = h(i,j) + htmp + w = w - htmp*v(:,i) + end do + end do + + hnext = sqrt(max(A%inner_product(w, w), zero_${s}$)) + h(j+1,j) = hnext + if (hnext > epsilon(one_${s}$)) then + v(:,j+1) = w / hnext + else + v(:,j+1) = zero_${s}$ + end if + + ! Apply the previously accumulated Givens rotations. + do i = 1, j - 1 + temp = cs(i) * h(i,j) + sn(i) * h(i+1,j) + h(i+1,j) = -sn(i) * h(i,j) + cs(i) * h(i+1,j) + h(i,j) = temp + end do + + ! Build and apply the next Givens rotation. + denom = sqrt(h(j,j) * h(j,j) + h(j+1,j) * h(j+1,j)) + if (denom > epsilon(one_${s}$)) then + cs(j) = h(j,j) / denom + sn(j) = h(j+1,j) / denom + else + cs(j) = one_${s}$ + sn(j) = zero_${s}$ + end if + + temp = cs(j) * h(j,j) + sn(j) * h(j+1,j) + h(j+1,j) = -sn(j) * h(j,j) + cs(j) * h(j+1,j) + h(j,j) = temp + + temp = cs(j) * g(j) + sn(j) * g(j+1) + g(j+1) = -sn(j) * g(j) + cs(j) * g(j+1) + g(j) = temp + + ! Cheap residual-norm estimate; no solution rebuild needed. + norm_sq = g(j+1) * g(j+1) + j_final = j + + if (norm_sq < tolsq .or. hnext <= epsilon(one_${s}$)) exit + end do + + ! Cycle-end update from the least-squares correction. + if (j_final > 0) then + call upper_triangular_solve(h, g, y, j_final) + if (compact) then + call gemv('N', m=size(x), n=j_final, alpha=one_${s}$, & + a=v, lda=size(v,1), & + x=y, incx=1, & + beta=zero_${s}$, y=w, incy=1) + call M%matvec(w, z(:,1), alpha=one_${s}$, beta=zero_${s}$, op='N') + x = x + z(:,1) + else + call gemv('N', m=size(x), n=j_final, alpha=one_${s}$, & + a=z, lda=size(z,1), & + x=y, incx=1, & + beta=one_${s}$, y=x, incy=1) + end if + end if + + ! Refresh the true residual so the convergence test per restart cycle and the logged + ! value use the true ||b - A*x||, not the Hessenberg estimate. + r = b + call A%matvec(x, r, alpha=-one_${s}$, beta=one_${s}$, op='N') + norm_sq = A%inner_product(r, r) + + if (associated(workspace%callback)) call workspace%callback(x, norm_sq, iter) + end do + end associate + + contains + + subroutine upper_triangular_solve(h, g, y, n) + ${t}$, intent(in) :: h(:,:), g(:) + ${t}$, intent(inout) :: y(:) + integer, intent(in) :: n + integer :: row + + y(1:n) = g(1:n) + do row = n, 1, -1 + if (row < n) y(row) = y(row) - dot_product(h(row,row+1:n) , y(row+1:n)) + if (abs(h(row,row)) > epsilon(one_${s}$)) then + y(row) = y(row) / h(row,row) + else + y(row) = zero_${s}$ + end if + end do + end subroutine + end subroutine + #:endfor + + #: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) + #:if matrix == "dense" + use stdlib_linalg, only: diag + ${t}$, intent(in) :: A(:,:) + #:else + type(${matrix}$_${s}$_type), intent(in) :: A + #:endif + ${t}$, intent(in) :: b(:) + ${t}$, intent(inout) :: x(:) + ${t}$, intent(in), optional :: rtol, atol + logical(int8), intent(in), optional, target :: di(:) + integer, intent(in), optional :: maxiter, kdim + logical, intent(in), optional :: restart + integer, intent(in), optional :: precond + class(stdlib_linop_${s}$_type), optional, intent(in), target :: M + type(stdlib_solver_workspace_${s}$_type), optional, intent(inout), target :: workspace + logical, intent(in), optional :: compact + type(stdlib_linop_${s}$_type) :: op + type(stdlib_linop_${s}$_type), pointer :: M_ => null() + type(stdlib_solver_workspace_${s}$_type), pointer :: workspace_ + integer :: kdim_, maxiter_, n, ncols, precond_ + ${t}$ :: rtol_, atol_ + logical :: compact_, restart_ + logical(int8), pointer :: di_(:) + ${t}$, allocatable :: diagonal(:) + + n = size(b) + maxiter_ = optval(x=maxiter, default=n) + kdim_ = max(1, min(optval(x=kdim, default=min(30, n)), n)) + restart_ = optval(x=restart, default=.true.) + compact_ = optval(x=compact, default=.true.) + rtol_ = optval(x=rtol, default=1.e-5_${s}$) + atol_ = optval(x=atol, default=epsilon(one_${s}$)) + precond_ = optval(x=precond, default=pc_none) + ncols = stdlib_size_wksp_gmres(kdim_,compact_) + + if (present(M)) then + M_ => M + else + allocate(M_) + allocate(diagonal(n), source=zero_${s}$) + + select case(precond_) + case(pc_jacobi) + #:if matrix == "dense" + diagonal = diag(A) + #:else + call diag(A, diagonal) + #:endif + M_%matvec => precond_jacobi + case default + M_%matvec => precond_none + end select + where(abs(diagonal) > epsilon(zero_${s}$)) diagonal = one_${s}$ / diagonal + end if + + op%matvec => matvec + + if (present(di)) then + di_ => di + else + allocate(di_(n), source=.false._int8) + end if + + if (present(workspace)) then + workspace_ => workspace + else + allocate(workspace_) + end if + if (.not.allocated(workspace_%tmp)) then + allocate(workspace_%tmp(n, ncols), source=zero_${s}$) + else if (size(workspace_%tmp,1) /= n .or. size(workspace_%tmp,2) < ncols) then + deallocate(workspace_%tmp) + allocate(workspace_%tmp(n, ncols), source=zero_${s}$) + end if + + if (restart_) x = zero_${s}$ + x = merge(b, x, di_) + call stdlib_solve_gmres_kernel(op, M_, b, x, rtol_, atol_, maxiter_, kdim_, workspace_, compact=compact_) + + if (.not.present(di)) deallocate(di_) + di_ => null() + + if (.not.present(workspace)) then + deallocate(workspace_%tmp) + deallocate(workspace_) + end if + M_ => null() + workspace_ => null() + + contains + + subroutine matvec(x,y,alpha,beta,op) + #:if matrix == "dense" + use stdlib_linalg_blas, only: gemv + #:endif + ${t}$, intent(in) :: x(:) + ${t}$, intent(inout) :: y(:) + ${t}$, intent(in) :: alpha + ${t}$, intent(in) :: beta + character(1), intent(in) :: op + #:if matrix == "dense" + call gemv(op, m=size(A,1), n=size(A,2), alpha=alpha, a=A, lda=size(A,1), x=x, incx=1, beta=beta, y=y, incy=1) + #:else + call spmv(A, x, y, alpha, beta, op) + #:endif + y = merge(zero_${s}$, y, di_) + end subroutine + + subroutine precond_none(x,y,alpha,beta,op) + ${t}$, intent(in) :: x(:) + ${t}$, intent(inout) :: y(:) + ${t}$, intent(in) :: alpha + ${t}$, intent(in) :: beta + character(1), intent(in) :: op + y = merge(zero_${s}$, x, di_) + end subroutine + + subroutine precond_jacobi(x,y,alpha,beta,op) + ${t}$, intent(in) :: x(:) + ${t}$, intent(inout) :: y(:) + ${t}$, intent(in) :: alpha + ${t}$, intent(in) :: beta + character(1), intent(in) :: op + y = merge(zero_${s}$, diagonal * x, di_) + end subroutine + end subroutine + #:endfor + #:endfor + +end submodule stdlib_linalg_iterative_gmres \ No newline at end of file diff --git a/test/linalg/test_linalg_solve_iterative.fypp b/test/linalg/test_linalg_solve_iterative.fypp index 8da08fcb8..e3c47b32c 100644 --- a/test/linalg/test_linalg_solve_iterative.fypp +++ b/test/linalg/test_linalg_solve_iterative.fypp @@ -23,6 +23,8 @@ module test_linalg_solve_iterative tests = [ new_unittest("stdlib_solve_cg",test_stdlib_solve_cg), & new_unittest("stdlib_solve_pcg",test_stdlib_solve_pcg), & + new_unittest("stdlib_solve_gmres",test_stdlib_solve_gmres), & + new_unittest("stdlib_solve_gmres_hilbert",test_stdlib_solve_gmres_hilbert), & new_unittest("stdlib_solve_bicgstab",test_stdlib_solve_bicgstab), & new_unittest("stdlib_solve_bicgstab_nonsymmetric",test_stdlib_solve_bicgstab_nonsymmetric) ] @@ -80,6 +82,101 @@ module test_linalg_solve_iterative #:endfor end subroutine test_stdlib_solve_pcg + subroutine test_stdlib_solve_gmres(error) + type(error_type), allocatable, intent(out) :: error + + #:for k, t, s in R_KINDS_TYPES + block + ${t}$, parameter :: tol = 1000*epsilon(0.0_${k}$) + ${t}$, parameter :: A(4,4) = reshape([${t}$ :: 5, 1, 2, 0, & + 0, 4, -1, 1, & + 1, 0, 3, 2, & + 2, -1, 0, 6], [4,4]) + ${t}$ :: x(4), load(4) + ${t}$, parameter :: xref(*) = [1.0_${k}$, 2.0_${k}$, -1.0_${k}$, 0.5_${k}$] + load = matmul(A, xref) + x = 0.0_${k}$ + + call stdlib_solve_gmres(A, load, x, rtol=tol, maxiter=12) + + call check(error, norm2(x-xref)