feat: [linalg] add iterative solvers

doc/specs/stdlib_linalg_iterative_solvers.md

@@ -95,7 +95,7 @@ Subroutine
 
 #### Description
 
-Implements the Conjugate Gradient (CG) method for solving the linear system \( Ax = b \), where \( A \) is a symmetric positive-definite linear operator defined via the `linop` type. This is the core implementation, allowing flexibility for custom matrix types or parallel environments.
+Implements the Conjugate Gradient (CG) method for solving the linear system \( Ax = b \), where \( A \) 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.
 
 #### Syntax
 
@@ -117,7 +117,7 @@ Subroutine
 
 `x`: 1-D array of `real(<kind>)` which serves as the input initial guess and the output solution. This argument is `intent(inout)`.
 
-`tol`: scalar of type `real(<kind>)` specifying the requested tolerance with respect to the relative residual vector norm. This argument is `intent(in)`.
+`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)`.
 
@@ -132,7 +132,7 @@ Provides a user-friendly interface to the CG method for solving \( Ax = b \), su
 
 #### Syntax
 
-`call ` [[stdlib_iterative_solvers(module):solve_cg(interface)]] ` (A, b, x [, di, tol, maxiter, restart, workspace])`
+`call ` [[stdlib_iterative_solvers(module):solve_cg(interface)]] ` (A, b, x [, di, rtol, atol, maxiter, restart, workspace])`
 
 #### Status
 
@@ -152,7 +152,7 @@ Subroutine
 
 `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)`.
 
-`tol` (optional): scalar of type `real(<kind>)` specifying the requested tolerance with respect to the relative residual vector norm. If no value is given, a default value of `1.e-4` is set. 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 ) \). Defaults values are `rtol=1.e-4` and `atol=0`. 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)`.
 
@@ -169,7 +169,7 @@ Subroutine
 
 #### Description
 
-Implements the Preconditioned Conjugate Gradient (PCG) method for solving the linear system \( Ax = b \), where \( A \) is a symmetric positive-definite linear operator defined via the `linop` type. This is the core implementation, allowing flexibility for custom matrix types or parallel environments.
+Implements the Preconditioned Conjugate Gradient (PCG) method for solving the linear system \( Ax = b \), where \( A \) 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.
 
 #### Syntax
 
@@ -193,7 +193,7 @@ Subroutine
 
 `x`: 1-D array of `real(<kind>)` which serves as the input initial guess and the output solution. This argument is `intent(inout)`.
 
-`tol`: scalar of type `real(<kind>)` specifying the requested tolerance with respect to the relative residual vector norm. 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 ) \). These arguments are `intent(in)`.
 
 `maxiter`: scalar of type `integer` defining the maximum allowed number of iterations. This argument is `intent(in)`.
 
@@ -234,7 +234,7 @@ Subroutine
 
 `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)`.
 
-`tol` (optional): scalar of type `real(<kind>)` specifying the requested tolerance with respect to the relative residual vector norm. If no value is given, a default value of `1.e-4` is set. 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 ) \). Defaults values are `rtol=1.e-4` and `atol=0`. 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)`.
 

example/linalg/example_solve_custom.f90

@@ -12,11 +12,12 @@ module custom_solver
 
 contains
 
-    subroutine stdlib_solve_pcg_custom(A,b,x,di,tol,maxiter,restart,workspace)
+    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 :: tol
+        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
@@ -26,7 +27,7 @@ subroutine stdlib_solve_pcg_custom(A,b,x,di,tol,maxiter,restart,workspace)
         type(stdlib_linop_dp_type) :: M
         type(stdlib_solver_workspace_dp_type), pointer :: workspace_
         integer :: n, maxiter_
-        real(dp) :: tol_
+        real(dp) :: rtol_, atol_
         logical :: restart_
         logical(int8), pointer :: di_(:)
         real(dp), allocatable :: diagonal(:)
@@ -35,7 +36,8 @@ subroutine stdlib_solve_pcg_custom(A,b,x,di,tol,maxiter,restart,workspace)
         n = size(b)
         maxiter_ = optval(x=maxiter, default=n)
         restart_ = optval(x=restart, default=.true.)
-        tol_     = optval(x=tol,     default=1.e-4_dp)
+        rtol_    = optval(x=rtol,    default=1.e-4_dp)
+        atol_    = optval(x=atol,    default=0._dp)
         norm_sq0 = 0._dp
         !-------------------------
         ! internal memory setup
@@ -61,7 +63,7 @@ subroutine stdlib_solve_pcg_custom(A,b,x,di,tol,maxiter,restart,workspace)
         where(abs(diagonal)>epsilon(0._dp)) diagonal = 1._dp/diagonal
         !-------------------------
         ! main call to the solver
-        call stdlib_solve_pcg_kernel(op,M,b,x,tol_,maxiter_,workspace_)
+        call stdlib_solve_pcg_kernel(op,M,b,x,rtol_,atol_,maxiter_,workspace_)
 
         !-------------------------
         ! internal memory cleanup
@@ -135,7 +137,7 @@ program example_solve_custom
     dirichlet = .false._int8
     dirichlet([1,5]) = .true._int8
 
-    call stdlib_solve_pcg_custom(laplacian_csr, rhs, x, tol=1.d-6, di=dirichlet)
+    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

example/linalg/example_solve_pcg.f90

@@ -23,10 +23,10 @@ program example_solve_pcg
     dirichlet = .false._int8
     dirichlet([1,5]) = .true._int8
 
-    call stdlib_solve_pcg(laplacian, rhs, x, tol=1.d-6, di=dirichlet)
+    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, tol=1.d-6, di=dirichlet)
+    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

src/stdlib_linalg_iterative_solvers.fypp

@@ -70,11 +70,12 @@ module stdlib_linalg_iterative_solvers
     !! [Specifications](../page/specs/stdlib_linalg_iterative_solvers.html#stdlib_solve_cg_kernel)
     interface stdlib_solve_cg_kernel
         #:for k, t, s in R_KINDS_TYPES
-        module subroutine stdlib_solve_cg_kernel_${s}$(A,b,x,tol,maxiter,workspace)
+        module subroutine stdlib_solve_cg_kernel_${s}$(A,b,x,rtol,atol,maxiter,workspace)
             class(stdlib_linop_${s}$_type), intent(in) :: A !! linear operator
             ${t}$, intent(in) :: b(:) !! right-hand side vector
             ${t}$, intent(inout) :: x(:) !! solution vector and initial guess
-            ${t}$, intent(in) :: tol !! tolerance for convergence
+            ${t}$, intent(in) :: rtol !! relative tolerance for convergence
+            ${t}$, intent(in) :: atol !! absolut tolerance for convergence
             integer, intent(in) :: maxiter !! maximum number of iterations
             type(stdlib_solver_workspace_${s}$_type), intent(inout) :: workspace !! workspace for the solver
         end subroutine
@@ -85,7 +86,7 @@ module stdlib_linalg_iterative_solvers
     interface stdlib_solve_cg
         #:for matrix in MATRIX_TYPES
         #:for k, t, s in R_KINDS_TYPES
-        module subroutine stdlib_solve_cg_${matrix}$_${s}$(A,b,x,di,tol,maxiter,restart,workspace)
+        module subroutine stdlib_solve_cg_${matrix}$_${s}$(A,b,x,di,rtol,atol,maxiter,restart,workspace)
             !! linear operator matrix
             #:if matrix == "dense"
             ${t}$, intent(in) :: A(:,:) 
@@ -94,7 +95,8 @@ module stdlib_linalg_iterative_solvers
             #:endif
             ${t}$, intent(in) :: b(:) !! right-hand side vector
             ${t}$, intent(inout) :: x(:) !! solution vector and initial guess
-            ${t}$, intent(in), optional :: tol !! tolerance for convergence
+            ${t}$, intent(in), optional :: rtol !! relative tolerance for convergence
+            ${t}$, intent(in), optional :: atol !! absolute tolerance for convergence
             logical(1), intent(in), optional, target  :: di(:) !! dirichlet conditions mask
             integer, intent(in), optional :: maxiter !! maximum number of iterations
             logical, intent(in), optional :: restart !! restart flag
@@ -111,12 +113,13 @@ module stdlib_linalg_iterative_solvers
     !! [Specifications](../page/specs/stdlib_linalg_iterative_solvers.html#stdlib_solve_pcg_kernel)
     interface stdlib_solve_pcg_kernel
         #:for k, t, s in R_KINDS_TYPES
-        module subroutine stdlib_solve_pcg_kernel_${s}$(A,M,b,x,tol,maxiter,workspace)
+        module subroutine stdlib_solve_pcg_kernel_${s}$(A,M,b,x,rtol,atol,maxiter,workspace)
             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) :: tol !! tolerance for convergence
+            ${t}$, intent(in) :: rtol !! relative tolerance for convergence
+            ${t}$, intent(in) :: atol !! absolute tolerance for convergence
             integer, intent(in) :: maxiter !! maximum number of iterations
             type(stdlib_solver_workspace_${s}$_type), intent(inout) :: workspace !! workspace for the solver
         end subroutine
@@ -127,7 +130,7 @@ module stdlib_linalg_iterative_solvers
     interface stdlib_solve_pcg
         #:for matrix in MATRIX_TYPES
         #:for k, t, s in R_KINDS_TYPES
-        module subroutine stdlib_solve_pcg_${matrix}$_${s}$(A,b,x,di,tol,maxiter,restart,precond,M,workspace)
+        module subroutine stdlib_solve_pcg_${matrix}$_${s}$(A,b,x,di,rtol,atol,maxiter,restart,precond,M,workspace)
             !! linear operator matrix
             #:if matrix == "dense"
             ${t}$, intent(in) :: A(:,:)
@@ -136,7 +139,8 @@ module stdlib_linalg_iterative_solvers
             #:endif
             ${t}$, intent(in) :: b(:) !! right-hand side vector
             ${t}$, intent(inout) :: x(:) !! solution vector and initial guess
-            ${t}$, intent(in), optional :: tol !! tolerance for convergence
+            ${t}$, intent(in), optional :: rtol !! relative tolerance for convergence
+            ${t}$, intent(in), optional :: atol !! absolute tolerance for convergence
             logical(1), intent(in), optional, target  :: di(:) !! dirichlet conditions mask
             integer, intent(in), optional  :: maxiter !! maximum number of iterations
             logical, intent(in), optional :: restart !! restart flag

src/stdlib_linalg_iterative_solvers_cg.fypp

@@ -15,9 +15,9 @@ submodule(stdlib_linalg_iterative_solvers) stdlib_linalg_iterative_cg
 contains
 
     #:for k, t, s in R_KINDS_TYPES
-    module subroutine stdlib_solve_cg_kernel_${s}$(A,b,x,tol,maxiter,workspace)
+    module subroutine stdlib_solve_cg_kernel_${s}$(A,b,x,rtol,atol,maxiter,workspace)
         class(stdlib_linop_${s}$_type), intent(in) :: A
-        ${t}$, intent(in) :: b(:), tol
+        ${t}$, intent(in) :: b(:), rtol, atol
         ${t}$, intent(inout) :: x(:)
         integer, intent(in) :: maxiter
         type(stdlib_solver_workspace_${s}$_type), intent(inout) :: workspace
@@ -40,10 +40,10 @@ contains
 
             P = R
 
-            tolsq = tol*tol
+            tolsq = max(rtol*rtol * norm_sq0, atol*atol)
             beta = zero_${s}$
             if(associated(workspace%callback)) call workspace%callback(x, norm_sq, iter)
-            do while( norm_sq > tolsq * norm_sq0 .and. iter < maxiter)
+            do while( norm_sq >= tolsq .and. iter < maxiter)
                 call A%matvec(P,Ap, alpha= one_${s}$, beta=zero_${s}$, op='N') ! Ap = A*P
 
                 alpha = norm_sq / A%inner_product(P, Ap)
@@ -67,15 +67,15 @@ contains
 
     #:for matrix in MATRIX_TYPES
     #:for k, t, s in R_KINDS_TYPES
-    module subroutine stdlib_solve_cg_${matrix}$_${s}$(A,b,x,di,tol,maxiter,restart,workspace)
+    module subroutine stdlib_solve_cg_${matrix}$_${s}$(A,b,x,di,rtol,atol,maxiter,restart,workspace)
         #:if matrix == "dense"
         ${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 :: tol
+        ${t}$, intent(in), optional :: rtol, atol
         logical(int8), intent(in), optional, target  :: di(:)
         integer, intent(in), optional :: maxiter
         logical, intent(in), optional  :: restart
@@ -84,14 +84,15 @@ contains
         type(stdlib_linop_${s}$_type) :: op
         type(stdlib_solver_workspace_${s}$_type), pointer :: workspace_
         integer :: n, maxiter_
-        ${t}$ :: tol_
+        ${t}$ :: rtol_, atol_
         logical :: restart_
         logical(int8), pointer :: di_(:)
         !-------------------------
         n = size(b)
         maxiter_ = optval(x=maxiter, default=n)
         restart_ = optval(x=restart, default=.true.)
-        tol_     = optval(x=tol,     default=1.e-4_${s}$)
+        rtol_    = optval(x=rtol,    default=1.e-4_${s}$)
+        atol_    = optval(x=atol,    default=zero_${s}$)
 
         !-------------------------
         ! internal memory setup
@@ -113,7 +114,7 @@ contains
         ! main call to the solver
         if(restart_) x = zero_${s}$
         x = merge( b, x, di_ ) ! copy dirichlet load conditions encoded in B and indicated by di
-        call stdlib_solve_cg_kernel(op,b,x,tol_,maxiter_,workspace_)
+        call stdlib_solve_cg_kernel(op,b,x,rtol_,atol_,maxiter_,workspace_)
 
         !-------------------------
         ! internal memory cleanup

src/stdlib_linalg_iterative_solvers_pcg.fypp

@@ -20,10 +20,10 @@ submodule(stdlib_linalg_iterative_solvers) stdlib_linalg_iterative_pcg
 contains
 
     #:for k, t, s in R_KINDS_TYPES
-    module subroutine stdlib_solve_pcg_kernel_${s}$(A,M,b,x,tol,maxiter,workspace)
+    module subroutine stdlib_solve_pcg_kernel_${s}$(A,M,b,x,rtol,atol,maxiter,workspace)
         class(stdlib_linop_${s}$_type), intent(in) :: A
         class(stdlib_linop_${s}$_type), intent(in) :: M 
-        ${t}$, intent(in) :: b(:), tol
+        ${t}$, intent(in) :: b(:), rtol, atol
         ${t}$, intent(inout) :: x(:)
         integer, intent(in) :: maxiter
         type(stdlib_solver_workspace_${s}$_type), intent(inout) :: workspace
@@ -48,11 +48,11 @@ contains
 
             call M%matvec(R,P, alpha= one_${s}$, beta=zero_${s}$, op='N') ! P = M^{-1}*R
 
-            tolsq = tol*tol
+            tolsq = max(rtol*rtol * norm_sq0, atol*atol)
 
             zr1 = zero_${s}$
             zr2 = one_${s}$
-            do while ( (iter < maxiter) .AND. (norm_sq > tolsq * norm_sq0) )
+            do while ( (iter < maxiter) .AND. (norm_sq >= tolsq) )
 
                 call M%matvec(R,S, alpha= one_${s}$, beta=zero_${s}$, op='N') ! S = M^{-1}*R
                 zr2 = A%inner_product( R, S )
@@ -84,7 +84,7 @@ contains
 
     #:for matrix in MATRIX_TYPES
     #:for k, t, s in R_KINDS_TYPES
-    module subroutine stdlib_solve_pcg_${matrix}$_${s}$(A,b,x,di,tol,maxiter,restart,precond,M,workspace)
+    module subroutine stdlib_solve_pcg_${matrix}$_${s}$(A,b,x,di,rtol,atol,maxiter,restart,precond,M,workspace)
         #:if matrix == "dense"
         use stdlib_linalg, only: diag
         ${t}$, intent(in) :: A(:,:)
@@ -93,7 +93,7 @@ contains
         #:endif
         ${t}$, intent(in) :: b(:)
         ${t}$, intent(inout) :: x(:)
-        ${t}$, intent(in), optional :: tol
+        ${t}$, intent(in), optional :: rtol, atol
         logical(int8), intent(in), optional, target  :: di(:)
         integer, intent(in), optional  :: maxiter
         logical, intent(in), optional  :: restart
@@ -105,7 +105,7 @@ contains
         type(stdlib_linop_${s}$_type), pointer :: M_ => null()
         type(stdlib_solver_workspace_${s}$_type), pointer :: workspace_
         integer :: n, maxiter_
-        ${t}$ :: tol_
+        ${t}$ :: rtol_, atol_
         logical :: restart_
         logical(int8), pointer :: di_(:)
         !-------------------------
@@ -117,7 +117,8 @@ contains
         n = size(b)
         maxiter_ = optval(x=maxiter, default=n)
         restart_ = optval(x=restart, default=.true.)
-        tol_     = optval(x=tol,     default=1.e-4_${s}$)
+        rtol_    = optval(x=rtol,    default=1.e-4_${s}$)
+        atol_    = optval(x=atol,    default=zero_${s}$)
         precond_ = optval(x=precond, default=pc_none)
         !-------------------------
         ! internal memory setup
@@ -162,7 +163,7 @@ contains
         ! main call to the solver
         if(restart_) x = zero_${s}$
         x = merge( b, x, di_ ) ! copy dirichlet load conditions encoded in B and indicated by di
-        call stdlib_solve_pcg_kernel(op,M_,b,x,tol_,maxiter_,workspace_)
+        call stdlib_solve_pcg_kernel(op,M_,b,x,rtol_,atol_,maxiter_,workspace_)
 
         !-------------------------
         ! internal memory cleanup

test/linalg/test_linalg_solve_iterative.fypp

@@ -69,7 +69,7 @@ module test_linalg_solve_iterative
         dirichlet = .false._int8
         dirichlet([1,5]) = .true._int8
 
-        call stdlib_solve_pcg(A, load, x, di=dirichlet, tol=1.e-6_${k}$)
+        call stdlib_solve_pcg(A, load, x, di=dirichlet, rtol=1.e-6_${k}$)
 
         call check(error, norm2(x-xref)<1.e-6_${k}$*norm2(xref), 'error in preconditionned conjugate gradient solver')
         if (allocated(error)) return