feat : [linalg] add bi-conjugate gradient stabilized method (#1034)

* implement bicgstab with copilot

* replace logical size indicator

* add example with wilkilson matrix

* solve conflict with master branch

* revert space

* fix documentation links and descriptions

* Update example/linalg/example_solve_bicgstab_wilkinson.f90

Co-authored-by: Jeremie Vandenplas <[email protected]>

* Update doc/specs/stdlib_linalg_iterative_solvers.md

Co-authored-by: Jeremie Vandenplas <[email protected]>

---------

Co-authored-by: Jeremie Vandenplas <[email protected]>

Author: jalvesz

doc/specs/stdlib_linalg_iterative_solvers.md

@@ -10,12 +10,12 @@ title: linalg_iterative_solvers
 
 The `stdlib_linalg_iterative_solvers` module provides base implementations for known iterative solver methods. Each method is exposed with two procedure flavors: 
 
-* A `solve_<method>_kernel` which holds the method's base implementation. The linear system argument is defined through a `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>_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 `solve_<method>` which proposes an off-the-shelf ready to use interface for `dense` and `CSR_<kind>_type` matrices for all `real` kinds.
+* 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.
 
 <!-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -->
-### The `linop` derived type
+### The `stdlib_linop` derived type
 
 The `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.
 
@@ -25,11 +25,11 @@ The following type-bound procedure pointers enable customization of the solver:
 
 ##### `matvec`
 
-Proxy procedure for the matrix-vector product $y = alpha * op(M) * x + beta * y$.
+Proxy procedure for the matrix-vector product \( y = alpha * op(M) * x + beta * y \).
 
 #### Syntax
 
-`call ` [[stdlib_iterative_solvers(module):matvec(interface)]] ` (x,y,alpha,beta,op)`
+`call ` [[stdlib_linalg_iterative_solvers(module):stdlib_linop_dp_type(type)]] `%matvec(x,y,alpha,beta,op)`
 
 ###### Class
 
@@ -53,7 +53,7 @@ Proxy procedure for the `dot_product`.
 
 #### Syntax
 
-`res = ` [[stdlib_iterative_solvers(module):inner_product(interface)]] ` (x,y)`
+`res = ` [[stdlib_linalg_iterative_solvers(module):stdlib_linop_dp_type(type)]] `%inner_product(x,y)`
 
 ###### Class
 
@@ -99,7 +99,7 @@ Implements the Conjugate Gradient (CG) method for solving the linear system \( A
 
 #### Syntax
 
-`call ` [[stdlib_iterative_solvers(module):stdlib_solve_cg_kernel(interface)]] ` (A, b, x, tol, maxiter, workspace)`
+`call ` [[stdlib_linalg_iterative_solvers(module):stdlib_solve_cg_kernel(interface)]] ` (A, b, x, tol, maxiter, workspace)`
 
 #### Status
 
@@ -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, rtol, atol, maxiter, restart, workspace])`
+`call ` [[stdlib_linalg_iterative_solvers(module):stdlib_solve_cg(interface)]] ` (A, b, x [, di, rtol, atol, maxiter, restart, workspace])`
 
 #### Status
 
@@ -173,7 +173,7 @@ Implements the Preconditioned Conjugate Gradient (PCG) method for solving the li
 
 #### Syntax
 
-`call ` [[stdlib_iterative_solvers(module):stdlib_solve_cg_kernel(interface)]] ` (A, M, b, x, tol, maxiter, workspace)`
+`call ` [[stdlib_linalg_iterative_solvers(module):stdlib_solve_cg_kernel(interface)]] ` (A, M, b, x, tol, maxiter, workspace)`
 
 #### Status
 
@@ -214,7 +214,7 @@ Provides a user-friendly interface to the PCG method for solving \( Ax = b \), s
 
 #### Syntax
 
-`call ` [[stdlib_iterative_solvers(module):stdlib_solve_pcg(interface)]] ` (A, b, x [, di, tol, maxiter, restart, precond, M, workspace])`
+`call ` [[stdlib_linalg_iterative_solvers(module):stdlib_solve_pcg(interface)]] ` (A, b, x [, di, tol, maxiter, restart, precond, M, workspace])`
 
 #### Status
 
@@ -248,4 +248,104 @@ Subroutine
 
 ```fortran
 {!example/linalg/example_solve_pcg.f90!}
+```
+
+<!-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -->
+### `stdlib_solve_bicgstab_kernel` subroutine
+
+#### Description
+
+Implements the Biconjugate Gradient Stabilized (BiCGSTAB) method for solving the linear system \( Ax = b \), where \( A \) 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.
+
+#### Syntax
+
+`call ` [[stdlib_linalg_iterative_solvers(module):stdlib_solve_bicgstab_kernel(interface)]] ` (A, M, b, x, rtol, atol, maxiter, workspace)`
+
+#### 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)`.
+
+`workspace`: scalar derived type of `type(stdlib_solver_workspace_<kind>_type)` holding the work array for the solver. This argument is `intent(inout)`.
+
+#### Note
+
+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` subroutine
+
+#### Description
+
+Provides a user-friendly interface to the BiCGSTAB method for solving \( Ax = b \), 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.
+
+#### Syntax
+
+`call ` [[stdlib_linalg_iterative_solvers(module):stdlib_solve_bicgstab(interface)]] ` (A, b, x [, di, rtol, atol, maxiter, restart, precond, M, workspace])`
+
+#### 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)`.
+
+`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)`.
+
+#### Note
+
+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.
+
+#### Example 1
+
+```fortran
+{!example/linalg/example_solve_bicgstab.f90!}
+```
+
+#### Example 2
+```fortran
+{!example/linalg/example_solve_bicgstab_wilkinson.f90!}
 ```

example/linalg/CMakeLists.txt

@@ -41,6 +41,8 @@ ADD_EXAMPLE(get_norm)
 ADD_EXAMPLE(solve1)
 ADD_EXAMPLE(solve2)
 ADD_EXAMPLE(solve3)
+ADD_EXAMPLE(solve_bicgstab)
+ADD_EXAMPLE(solve_bicgstab_wilkinson)
 ADD_EXAMPLE(solve_cg)
 ADD_EXAMPLE(solve_pcg)
 ADD_EXAMPLE(solve_custom)

example/linalg/example_solve_bicgstab.f90

@@ -0,0 +1,41 @@
+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 

example/linalg/example_solve_bicgstab_wilkinson.f90

@@ -0,0 +1,52 @@
+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

src/CMakeLists.txt

@@ -37,6 +37,7 @@ set(fppFiles
     stdlib_linalg_qr.fypp
     stdlib_linalg_inverse.fypp
     stdlib_linalg_iterative_solvers.fypp
+    stdlib_linalg_iterative_solvers_bicgstab.fypp
     stdlib_linalg_iterative_solvers_cg.fypp
     stdlib_linalg_iterative_solvers_pcg.fypp
     stdlib_linalg_pinv.fypp