stdlib_linalg_iterative_solvers moduleThe stdlib_linalg_iterative_solvers module provides base implementations for known iterative solver methods. Each method is exposed with two procedure flavors:
A stdlib_solve_<method>_kernel which holds the method's base implementation. The linear system argument is defined through a stdlib_linop derived type which enables extending the method for implicit or unknown (by stdlib) matrices or to complex scenarios involving distributed parallelism for which the user shall extend the inner_product and/or matrix-vector product to account for parallel syncrhonization.
A stdlib_solve_<method> which proposes an off-the-shelf ready to use interface for dense and CSR_<kind>_type matrices for all real kinds.
stdlib_linop derived typeThe stdlib_linop_<kind>_type derive type is an auxiliary class enabling to abstract the definition of the linear system and the actual implementation of the solvers.
The following type-bound procedure pointers enable customization of the solver:
matvecProxy procedure for the matrix-vector product .
call stdlib_linop_dp_type %matvec(x,y,alpha,beta,op)
Subroutine
x: 1-D array of real(<kind>). This argument is intent(in).
y: 1-D array of real(<kind>). This argument is intent(inout).
alpha: scalar of real(<kind>). This argument is intent(in).
beta: scalar of real(<kind>). This argument is intent(in).
op: character(1) scalar which can be have any of the following values: N (no transpose), T (transpose) or H (conjugate transpose). This argument is intent(in).
inner_productProxy procedure for the dot_product.
res = stdlib_linop_dp_type %inner_product(x,y)
Function
x: 1-D array of real(<kind>). This argument is intent(in).
y: 1-D array of real(<kind>). This argument is intent(in).
The output is a scalar of type and kind same as to that of x and y.
solver_workspace derived typeThe stdlib_solver_workspace_<kind>_type derive type is an auxiliary class enabling to hold the data associated to the working arrays needed by the solvers to operate.
callback: null pointer procedure enabling to pass a callback at each iteration to check on the solvers status.Subroutine
x: 1-D array of real(<kind>) type with the current state of the solution vector. This argument is intent(in) as it should not be modified by the callback.
norm_sq: scalar of real(<kind>) type representing the squared norm of the residual at the current iteration. This argument is intent(in).
iter: scalar of integer type giving the current iteration counter. This argument is intent(in).
stdlib_solve_cg_kernel subroutineImplements the Conjugate Gradient (CG) method for solving the linear system , where is a symmetric positive-definite linear operator defined via the stdlib_linop type. This is the core implementation, allowing flexibility for custom matrix types or parallel environments.
call stdlib_solve_cg_kernel (A, b, x, tol, maxiter, workspace)
Experimental
Subroutine
A: class(stdlib_linop_<kind>_type) defining the 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 . These arguments are intent(in).
maxiter: scalar of type integer defining the maximum allowed number of iterations. 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).
solve_cg subroutineProvides a user-friendly interface to the CG method for solving , supporting dense and CSR_<kind>_type matrices. It handles workspace allocation and optional parameters for customization.
call stdlib_solve_cg (A, b, x [, di, rtol, atol, maxiter, restart, workspace])
Experimental
Subroutine
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 right-hand-side (or loading) of the linear system. This argument is intent(in).
x: 1-D array of real(<kind>) which serves as the input initial guess and as 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 conditions 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 . Defaults 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).
workspace (optional): scalar derived type of type(stdlib_solver_workspace_<kind>_type) holding the work array for the solver. If the user passes its own workspace, then a pointer is set internally to it. Otherwise, memory will be internally allocated and deallocated before exiting the procedure. This argument is intent(inout).
program example_solve_cg
use stdlib_kinds, only: int8, dp
use stdlib_linalg_iterative_solvers, only: stdlib_solve_cg
real(dp) :: matrix(2,2)
real(dp) :: x(2), rhs(2)
matrix = reshape( [4, 1,&
1, 3] , [2,2])
x = dble( [2,1] ) !> initial guess
rhs = dble( [1,2] ) !> rhs vector
call stdlib_solve_cg(matrix, rhs, x, restart=.false.)
print *, x !> solution: [0.0909, 0.6364]
end program
stdlib_solve_pcg_kernel subroutineImplements the Preconditioned Conjugate Gradient (PCG) method for solving the linear system , where is a symmetric positive-definite linear operator defined via the stdlib_linop type. This is the core implementation, allowing flexibility for custom matrix types or parallel environments.
call stdlib_solve_cg_kernel (A, M, b, x, tol, maxiter, workspace)
Experimental
Subroutine
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 (optional): scalars of type real(<kind>) specifying the convergence test. For convergence, the following criterion is used . These arguments are intent(in).
maxiter: scalar of type integer defining the maximum allowed number of iterations. 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).
module custom_solver
use stdlib_kinds, only: int8, dp
use stdlib_sparse, only: CSR_dp_type, spmv, diag
use stdlib_linalg_iterative_solvers, only: stdlib_linop_dp_type, &
stdlib_solver_workspace_dp_type, &
stdlib_solve_pcg_kernel, &
stdlib_size_wksp_pcg
use stdlib_optval, only: optval
implicit none
private
public :: stdlib_solve_pcg_custom
contains
subroutine stdlib_solve_pcg_custom(A,b,x,di,rtol,atol,maxiter,restart,workspace)
type(CSR_dp_type), intent(in) :: A
real(dp), intent(in) :: b(:)
real(dp), intent(inout) :: x(:)
real(dp), intent(in), optional :: rtol
real(dp), intent(in), optional :: atol
logical(int8), intent(in), optional, target :: di(:)
integer, intent(in), optional :: maxiter
logical, intent(in), optional :: restart
type(stdlib_solver_workspace_dp_type), optional, intent(inout), target :: workspace
!-------------------------
type(stdlib_linop_dp_type) :: op
type(stdlib_linop_dp_type) :: M
type(stdlib_solver_workspace_dp_type), pointer :: workspace_
integer :: n, maxiter_
real(dp) :: rtol_, atol_
logical :: restart_
logical(int8), pointer :: di_(:)
real(dp), allocatable :: diagonal(:)
real(dp) :: norm_sq0
!-------------------------
n = size(b)
maxiter_ = optval(x=maxiter, default=n)
restart_ = optval(x=restart, default=.true.)
rtol_ = optval(x=rtol, default=1.e-4_dp)
atol_ = optval(x=atol, default=0._dp)
norm_sq0 = 0._dp
!-------------------------
! internal memory setup
op%matvec => my_matvec
op%inner_product => my_dot
M%matvec => my_jacobi_preconditioner
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)) allocate( workspace_%tmp(n,stdlib_size_wksp_pcg) , source = 0._dp )
workspace_%callback => my_logger
!-------------------------
! Jacobi preconditioner factorization
call diag(A,diagonal)
where(abs(diagonal)>epsilon(0._dp)) diagonal = 1._dp/diagonal
!-------------------------
! main call to the solver
call stdlib_solve_pcg_kernel(op,M,b,x,rtol_,atol_,maxiter_,workspace_)
!-------------------------
! internal memory cleanup
if(.not.present(di)) deallocate(di_)
di_ => null()
if(.not.present(workspace)) then
deallocate( workspace_%tmp )
deallocate( workspace_ )
end if
workspace_ => null()
contains
subroutine my_matvec(x,y,alpha,beta,op)
real(dp), intent(in) :: x(:)
real(dp), intent(inout) :: y(:)
real(dp), intent(in) :: alpha
real(dp), intent(in) :: beta
character(1), intent(in) :: op
call spmv( A , x, y , alpha, beta , op)
y = merge( 0._dp, y, di_ )
end subroutine
subroutine my_jacobi_preconditioner(x,y,alpha,beta,op)
real(dp), intent(in) :: x(:)
real(dp), intent(inout) :: y(:)
real(dp), intent(in) :: alpha
real(dp), intent(in) :: beta
character(1), intent(in) :: op
y = merge( 0._dp, diagonal * x , di_ )
end subroutine
real(dp) function my_dot(x,y) result(r)
real(dp), intent(in) :: x(:)
real(dp), intent(in) :: y(:)
r = dot_product(x,y)
end function
subroutine my_logger(x,norm_sq,iter)
real(dp), intent(in) :: x(:)
real(dp), intent(in) :: norm_sq
integer, intent(in) :: iter
if(iter == 0) norm_sq0 = norm_sq
print *, "Iteration: ", iter, " Residual: ", sqrt(norm_sq), " Relative: ", sqrt(norm_sq)/sqrt(norm_sq0)
end subroutine
end subroutine
end module custom_solver
program example_solve_custom
use custom_solver
use stdlib_kinds, only: int8, dp
use stdlib_sparse, only: CSR_dp_type, COO_dp_type, dense2coo, coo2csr
implicit none
type(CSR_dp_type) :: laplacian_csr
type(COO_dp_type) :: COO
real(dp) :: laplacian(5,5)
real(dp) :: x(5), rhs(5)
logical(int8) :: dirichlet(5)
laplacian = reshape( [1, -1, 0, 0, 0,&
-1, 2, -1, 0, 0,&
0, -1, 2, -1, 0,&
0, 0, -1, 2, -1,&
0, 0, 0, -1, 1] , [5,5])
call dense2coo(laplacian,COO)
call coo2csr(COO,laplacian_csr)
x = 0._dp
rhs = dble( [0,0,5,0,0] )
dirichlet = .false._int8
dirichlet([1,5]) = .true._int8
call stdlib_solve_pcg_custom(laplacian_csr, rhs, x, rtol=1.d-6, di=dirichlet)
print *, x !> solution: [0.0, 2.5, 5.0, 2.5, 0.0]
end program example_solve_custom
stdlib_solve_pcg subroutineProvides a user-friendly interface to the PCG method for solving , supporting dense and CSR_<kind>_type matrices. It supports optional preconditioners and handles workspace allocation.
call stdlib_solve_pcg (A, b, x [, di, tol, maxiter, restart, precond, M, workspace])
Experimental
Subroutine
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 conditions 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 . Defaults 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).
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 preconditionning 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, a pointer is set to this custom preconditioner.
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).
program example_solve_pcg
use stdlib_kinds, only: int8, dp
use stdlib_sparse
use stdlib_linalg_iterative_solvers, only: stdlib_solve_pcg
type(CSR_dp_type) :: laplacian_csr
type(COO_dp_type) :: COO
real(dp) :: laplacian(5,5)
real(dp) :: x(5), rhs(5)
logical(int8) :: dirichlet(5)
laplacian = reshape( [1, -1, 0, 0, 0,&
-1, 2, -1, 0, 0,&
0, -1, 2, -1, 0,&
0, 0, -1, 2, -1,&
0, 0, 0, -1, 1] , [5,5])
call dense2coo(laplacian,COO)
call coo2csr(COO,laplacian_csr)
x = 0._dp
rhs = real( [0,0,5,0,0], kind=dp )
dirichlet = .false._int8
dirichlet([1,5]) = .true._int8
call stdlib_solve_pcg(laplacian, rhs, x, rtol=1.d-6, di=dirichlet)
print *, x !> solution: [0.0, 2.5, 5.0, 2.5, 0.0]
x = 0._dp
call stdlib_solve_pcg(laplacian_csr, rhs, x, rtol=1.d-6, di=dirichlet)
print *, x !> solution: [0.0, 2.5, 5.0, 2.5, 0.0]
end program
stdlib_solve_bicgstab_kernel subroutineImplements the Biconjugate Gradient Stabilized (BiCGSTAB) method for solving the linear system , where is a general (non-symmetric) linear operator defined via the stdlib_linop type. BiCGSTAB is particularly suitable for solving non-symmetric linear systems and provides better stability than the basic BiCG method. This is the core implementation, allowing flexibility for custom matrix types or parallel environments.
call stdlib_solve_bicgstab_kernel (A, M, b, x, rtol, atol, maxiter, workspace)
Experimental
Subroutine
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 . These arguments are intent(in).
maxiter: scalar of type integer defining the maximum allowed number of iterations. 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).
The BiCGSTAB method requires 8 auxiliary vectors in its workspace, making it more memory-intensive than CG or PCG methods. However, it can handle general non-symmetric matrices and often converges faster than BiCG for many problems.
stdlib_solve_bicgstab subroutineProvides a user-friendly interface to the BiCGSTAB method for solving , supporting dense and CSR_<kind>_type matrices. BiCGSTAB is suitable for general (non-symmetric) linear systems and supports optional preconditioners for improved convergence. It handles workspace allocation and optional parameters for customization.
call stdlib_solve_bicgstab (A, b, x [, di, rtol, atol, maxiter, restart, precond, M, workspace])
Experimental
Subroutine
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 . 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).
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).
BiCGSTAB is particularly effective for: - Non-symmetric linear systems - Systems where CG cannot be applied - Cases where BiCG suffers from irregular convergence
The method uses 8 auxiliary vectors internally, requiring more memory than simpler methods but often providing better stability and convergence properties.
program example_solve_bicgstab
use stdlib_kinds, only: dp
use stdlib_linalg_iterative_solvers
implicit none
integer, parameter :: n = 4
real(dp) :: A(n,n), b(n), x(n)
integer :: i
! Example matrix (same as SciPy documentation)
A = reshape([4.0_dp, 2.0_dp, 0.0_dp, 1.0_dp, &
3.0_dp, 0.0_dp, 0.0_dp, 2.0_dp, &
0.0_dp, 1.0_dp, 1.0_dp, 1.0_dp, &
0.0_dp, 2.0_dp, 1.0_dp, 0.0_dp], [n,n])
b = [-1.0_dp, -0.5_dp, -1.0_dp, 2.0_dp]
x = 0.0_dp ! Initial guess
print *, 'Solving Ax = b using BiCGSTAB method:'
print *, 'Matrix A:'
do i = 1, n
print '(4F8.2)', A(i,:)
end do
print *, 'Right-hand side b:'
print '(4F8.2)', b
! Solve using BiCGSTAB
call stdlib_solve_bicgstab(A, b, x, rtol=1e-10_dp, atol=1e-12_dp)
print *, 'Solution x:'
print '(4F10.6)', x
! Verify solution
print *, 'Verification A*x:'
print '(4F10.6)', matmul(A, x)
print *, 'Residual ||b - A*x||:'
print *, norm2(b - matmul(A, x))
end program
program example_solve_bicgstab_wilkinson
use stdlib_kinds, only: dp
use stdlib_linalg_iterative_solvers
use stdlib_sparse
use stdlib_sparse_spmv
implicit none
integer, parameter :: n = 21
type(COO_dp_type) :: COO
type(CSR_dp_type) :: A
type(stdlib_solver_workspace_dp_type) :: workspace
real(dp) :: b(n), x(n), norm_sq0
integer :: i
! Construct the Wilkinson's matrix in COO format
! https://en.wikipedia.org/wiki/Wilkinson_matrix
call COO%malloc(n, n, n + 2*(n-1))
COO%data(1:n) = [( real(abs(i), kind=dp), i=-(n-1)/2, (n-1)/2)]
COO%data(n+1:) = 1.0_dp
do i = 1, n
COO%index(1:2, i) = [i,i]
end do
do i = 1, n-1
COO%index(1:2, n+i) = [i,i+1]
COO%index(1:2, n+n-1+i) = [i+1,i]
end do
call coo2ordered(COO,sort_data=.true.)
! Convert COO to CSR format
call coo2csr(COO, A)
! Set up the right-hand side for the solution to be ones
b = 0.0_dp
x = 1.0_dp
call spmv(A, x, b) ! b = A * 1
x = 0.0_dp ! initial guess
! Solve the system using BiCGSTAB
workspace%callback => my_logger
call stdlib_solve_bicgstab(A, b, x, rtol=1e-12_dp, maxiter=100, workspace=workspace)
contains
subroutine my_logger(x,norm_sq,iter)
real(dp), intent(in) :: x(:)
real(dp), intent(in) :: norm_sq
integer, intent(in) :: iter
if(iter == 0) norm_sq0 = norm_sq
print *, "Iteration: ", iter, " Residual: ", norm_sq, " Relative: ", norm_sq/norm_sq0
end subroutine
end program example_solve_bicgstab_wilkinson
stdlib_solve_gmres_kernel subroutineImplements the restarted Generalized Minimal Residual (GMRES) method for solving the linear system , where 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.
call stdlib_solve_gmres_kernel (A, M, b, x, rtol, atol, maxiter, kdim, workspace, compact)
Experimental
Subroutine
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 . 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).
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 subroutineProvides a user-friendly interface to the restarted GMRES method for solving , 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.
call stdlib_solve_gmres (A, b, x [, di, rtol, atol, maxiter, restart, kdim, precond, M, workspace, compact])
Experimental
Subroutine
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 . 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).
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.
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