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)
Pure 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
pure 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