implement bicgstab with copilot

Author: jalvesz

doc/specs/stdlib_linalg_iterative_solvers.md

@@ -248,4 +248,99 @@ 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_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_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 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 \( || 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
+
+```fortran
+{!example/linalg/example_solve_bicgstab.f90!}
 ```

example/linalg/CMakeLists.txt

@@ -41,6 +41,7 @@ ADD_EXAMPLE(get_norm)
 ADD_EXAMPLE(solve1)
 ADD_EXAMPLE(solve2)
 ADD_EXAMPLE(solve3)
+ADD_EXAMPLE(solve_bicgstab)
 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 

src/CMakeLists.txt

@@ -98,6 +98,7 @@ set(cppFiles
     stdlib_linalg_iterative_solvers.fypp
     stdlib_linalg_iterative_solvers_cg.fypp
     stdlib_linalg_iterative_solvers_pcg.fypp
+    stdlib_linalg_iterative_solvers_bicgstab.fypp
 )
 
 add_subdirectory(blas)

src/stdlib_linalg_iterative_solvers.fypp

@@ -15,9 +15,16 @@ module stdlib_linalg_iterative_solvers
     enum, bind(c)
         enumerator :: stdlib_size_wksp_cg = 3
         enumerator :: stdlib_size_wksp_pcg = 4
+        enumerator :: stdlib_size_wksp_bicgstab = 8
     end enum
-    public :: stdlib_size_wksp_cg, stdlib_size_wksp_pcg
+    public :: stdlib_size_wksp_cg, stdlib_size_wksp_pcg, stdlib_size_wksp_bicgstab
 
+    enum, bind(c)
+        enumerator :: pc_none = 0
+        enumerator :: pc_jacobi
+    end enum
+    public :: pc_none, pc_jacobi
+    
     !! version: experimental
     !!
     !! linop type holding the linear operator and its associated methods.
@@ -127,6 +134,26 @@ module stdlib_linalg_iterative_solvers
     end interface
     public :: stdlib_solve_pcg_kernel
 
+    !! version: experimental
+    !!
+    !! stdlib_solve_bicgstab_kernel interface for the biconjugate gradient stabilized method.
+    !! [Specifications](../page/specs/stdlib_linalg_iterative_solvers.html#stdlib_solve_bicgstab_kernel)
+    interface stdlib_solve_bicgstab_kernel
+        #:for k, t, s in R_KINDS_TYPES
+        module subroutine stdlib_solve_bicgstab_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) :: 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
+        #:endfor
+    end interface
+    public :: stdlib_solve_bicgstab_kernel
+
     interface stdlib_solve_pcg
         #:for matrix in MATRIX_TYPES
         #:for k, t, s in R_KINDS_TYPES
@@ -153,6 +180,32 @@ module stdlib_linalg_iterative_solvers
     end interface
     public :: stdlib_solve_pcg 
 
+    interface stdlib_solve_bicgstab
+        #:for matrix in MATRIX_TYPES
+        #:for k, t, s in R_KINDS_TYPES
+        module subroutine stdlib_solve_bicgstab_${matrix}$_${s}$(A,b,x,di,rtol,atol,maxiter,restart,precond,M,workspace)
+            !! 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(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
+            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
+        end subroutine
+        #:endfor
+        #:endfor
+    end interface
+    public :: stdlib_solve_bicgstab
+
 contains
 
     !------------------------------------------------------------------