diff --git a/doc/specs/stdlib_linalg_iterative_solvers.md b/doc/specs/stdlib_linalg_iterative_solvers.md index 6a0b69a95..949f55a44 100644 --- a/doc/specs/stdlib_linalg_iterative_solvers.md +++ b/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__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__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_` which proposes an off-the-shelf ready to use interface for `dense` and `CSR__type` matrices for all `real` kinds. +* A `stdlib_solve_` which proposes an off-the-shelf ready to use interface for `dense` and `CSR__type` matrices for all `real` kinds. -### The `linop` derived type +### The `stdlib_linop` derived type The `stdlib_linop__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__type)` defining the linear operator. This argument is `intent(in)`. + +`M`: `class(stdlib_linop__type)` defining the preconditioner linear operator. This argument is `intent(in)`. + +`b`: 1-D array of `real()` defining the loading conditions of the linear system. This argument is `intent(in)`. + +`x`: 1-D array of `real()` which serves as the input initial guess and the output solution. This argument is `intent(inout)`. + +`rtol` and `atol`: scalars of type `real()` 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__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__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__type` matrix defining the linear system. This argument is `intent(in)`. + +`b`: 1-D array of `real()` defining the loading conditions of the linear system. This argument is `intent(in)`. + +`x`: 1-D array of `real()` 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()` 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._)`. 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__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__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!} ``` \ No newline at end of file diff --git a/example/linalg/CMakeLists.txt b/example/linalg/CMakeLists.txt index 2aa9ca2e8..134999a8b 100644 --- a/example/linalg/CMakeLists.txt +++ b/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) diff --git a/example/linalg/example_solve_bicgstab.f90 b/example/linalg/example_solve_bicgstab.f90 new file mode 100644 index 000000000..7fdb62070 --- /dev/null +++ b/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 diff --git a/example/linalg/example_solve_bicgstab_wilkinson.f90 b/example/linalg/example_solve_bicgstab_wilkinson.f90 new file mode 100644 index 000000000..7e1c54ddc --- /dev/null +++ b/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 \ No newline at end of file diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 9280d3397..6b255ae9b 100644 --- a/src/CMakeLists.txt +++ b/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 diff --git a/src/stdlib_linalg_iterative_solvers.fypp b/src/stdlib_linalg_iterative_solvers.fypp index 70e412732..1b216a588 100644 --- a/src/stdlib_linalg_iterative_solvers.fypp +++ b/src/stdlib_linalg_iterative_solvers.fypp @@ -15,11 +15,20 @@ 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 !! + !! [Specifications](../page/specs/stdlib_linalg_iterative_solvers.html#stdlib_linop) + !! !! linop type holding the linear operator and its associated methods. !! The `linop` type is used to define the linear operator for the iterative solvers. #:for k, t, s in R_KINDS_TYPES @@ -31,6 +40,8 @@ module stdlib_linalg_iterative_solvers !! version: experimental !! + !! [Specifications](../page/specs/stdlib_linalg_iterative_solvers.html#stdlib_solver_workspace) + !! !! solver_workspace type holding temporal array data for the iterative solvers. #:for k, t, s in R_KINDS_TYPES type, public :: stdlib_solver_workspace_${s}$_type @@ -82,7 +93,10 @@ module stdlib_linalg_iterative_solvers #:endfor end interface public :: stdlib_solve_cg_kernel - + + !! version: experimental + !! + !! [Specifications](../page/specs/stdlib_linalg_iterative_solvers.html#stdlib_solve_cg) interface stdlib_solve_cg #:for matrix in MATRIX_TYPES #:for k, t, s in R_KINDS_TYPES @@ -97,7 +111,7 @@ module stdlib_linalg_iterative_solvers ${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 + logical(int8), intent(in), optional, target :: di(:) !! dirichlet conditions mask integer, intent(in), optional :: maxiter !! maximum number of iterations logical, intent(in), optional :: restart !! restart flag type(stdlib_solver_workspace_${s}$_type), optional, intent(inout), target :: workspace !! workspace for the solver @@ -127,6 +141,29 @@ 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 + + !! version: experimental + !! + !! [Specifications](../page/specs/stdlib_linalg_iterative_solvers.html#stdlib_solve_pcg) interface stdlib_solve_pcg #:for matrix in MATRIX_TYPES #:for k, t, s in R_KINDS_TYPES @@ -141,7 +178,7 @@ module stdlib_linalg_iterative_solvers ${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 + logical(int8), 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 @@ -153,6 +190,35 @@ module stdlib_linalg_iterative_solvers end interface public :: stdlib_solve_pcg + !! version: experimental + !! + !! [Specifications](../page/specs/stdlib_linalg_iterative_solvers.html#stdlib_solve_bicgstab) + 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(int8), 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 !------------------------------------------------------------------ diff --git a/src/stdlib_linalg_iterative_solvers_bicgstab.fypp b/src/stdlib_linalg_iterative_solvers_bicgstab.fypp new file mode 100644 index 000000000..73d3e7a09 --- /dev/null +++ b/src/stdlib_linalg_iterative_solvers_bicgstab.fypp @@ -0,0 +1,252 @@ +#:include "common.fypp" +#:set R_KINDS_TYPES = list(zip(REAL_KINDS, REAL_TYPES, REAL_SUFFIX)) +#:set C_KINDS_TYPES = list(zip(CMPLX_KINDS, CMPLX_TYPES, CMPLX_SUFFIX)) +#:set MATRIX_TYPES = ["dense", "CSR"] +#:set RANKS = range(1, 2+1) + +submodule(stdlib_linalg_iterative_solvers) stdlib_linalg_iterative_bicgstab + use stdlib_kinds + use stdlib_sparse + use stdlib_constants + use stdlib_linalg_iterative_solvers + use stdlib_optval, only: optval + implicit none + +contains + + #: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 + class(stdlib_linop_${s}$_type), intent(in) :: M + ${t}$, intent(in) :: b(:), rtol, atol + ${t}$, intent(inout) :: x(:) + integer, intent(in) :: maxiter + type(stdlib_solver_workspace_${s}$_type), intent(inout) :: workspace + !------------------------- + integer :: iter + ${t}$ :: norm_sq, norm_sq0, tolsq + ${t}$ :: rho, rho_prev, alpha, beta, omega, rv + ${t}$, parameter :: rhotol = epsilon(one_${s}$)**2 + ${t}$, parameter :: omegatol = epsilon(one_${s}$)**2 + !------------------------- + iter = 0 + associate( R => workspace%tmp(:,1), & + Rt => workspace%tmp(:,2), & + P => workspace%tmp(:,3), & + Pt => workspace%tmp(:,4), & + V => workspace%tmp(:,5), & + S => workspace%tmp(:,6), & + St => workspace%tmp(:,7), & + T => workspace%tmp(:,8)) + + norm_sq = A%inner_product( b, b ) + norm_sq0 = norm_sq + if(associated(workspace%callback)) call workspace%callback(x, norm_sq, iter) + + if ( norm_sq0 > zero_${s}$ ) then + + ! Compute initial residual: r = b - A*x + R = B + call A%matvec(X, R, alpha= -one_${s}$, beta=one_${s}$, op='N') ! R = B - A*X + + ! Choose arbitrary Rt (often Rt = r0) + Rt = R + + tolsq = max(rtol*rtol * norm_sq0, atol*atol) + + rho_prev = one_${s}$ + alpha = one_${s}$ + omega = one_${s}$ + + do while ( (iter < maxiter) .AND. (norm_sq >= tolsq) ) + + rho = A%inner_product( Rt, R ) + + ! Check for rho breakdown + if (abs(rho) < rhotol) exit + + if (iter > 0) then + ! Check for omega breakdown + if (abs(omega) < omegatol) exit + + beta = (rho / rho_prev) * (alpha / omega) + P = R + beta * (P - omega * V) + else + P = R + end if + + ! Preconditioned BiCGSTAB step + call M%matvec(P, Pt, alpha=one_${s}$, beta=zero_${s}$, op='N') ! Pt = M^{-1}*P + call A%matvec(Pt, V, alpha=one_${s}$, beta=zero_${s}$, op='N') ! V = A*Pt + + rv = A%inner_product( Rt, V ) + if (abs(rv) < epsilon(one_${s}$)) exit ! rv breakdown + + alpha = rho / rv + + ! Update residual: s = r - alpha*v + S = R - alpha * V + + ! Check if s is small enough + norm_sq = A%inner_product( S, S ) + if (norm_sq < tolsq) then + X = X + alpha * Pt + exit + end if + + ! Preconditioned update for t = A * M^{-1} * s + call M%matvec(S, St, alpha=one_${s}$, beta=zero_${s}$, op='N') ! St = M^{-1}*S + call A%matvec(St, T, alpha=one_${s}$, beta=zero_${s}$, op='N') ! T = A*St + + ! Compute omega + omega = A%inner_product( T, S ) / A%inner_product( T, T ) + + ! Update solution and residual + X = X + alpha * Pt + omega * St + R = S - omega * T + + norm_sq = A%inner_product( R, R ) + rho_prev = rho + iter = iter + 1 + if(associated(workspace%callback)) call workspace%callback(x, norm_sq, iter) + end do + end if + end associate + end subroutine + #:endfor + + #: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) + #:if matrix == "dense" + use stdlib_linalg, only: diag + ${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 :: rtol, atol + logical(int8), intent(in), optional, target :: di(:) + integer, intent(in), optional :: maxiter + logical, intent(in), optional :: restart + integer, intent(in), optional :: precond + class(stdlib_linop_${s}$_type), optional , intent(in), target :: M + type(stdlib_solver_workspace_${s}$_type), optional, intent(inout), target :: workspace + !------------------------- + type(stdlib_linop_${s}$_type) :: op + type(stdlib_linop_${s}$_type), pointer :: M_ => null() + type(stdlib_solver_workspace_${s}$_type), pointer :: workspace_ + integer :: n, maxiter_ + ${t}$ :: rtol_, atol_ + logical :: restart_ + logical(int8), pointer :: di_(:) + !------------------------- + ! working data for preconditioner + integer :: precond_ + ${t}$, allocatable :: diagonal(:) + + !------------------------- + n = size(b) + maxiter_ = optval(x=maxiter, default=n) + restart_ = optval(x=restart, default=.true.) + rtol_ = optval(x=rtol, default=1.e-5_${s}$) + atol_ = optval(x=atol, default=epsilon(one_${s}$)) + precond_ = optval(x=precond, default=pc_none) + !------------------------- + ! internal memory setup + ! preconditioner + if(present(M)) then + M_ => M + else + allocate( M_ ) + allocate(diagonal(n),source=zero_${s}$) + + select case(precond_) + case(pc_jacobi) + #:if matrix == "dense" + diagonal = diag(A) + #:else + call diag(A,diagonal) + #:endif + M_%matvec => precond_jacobi + case default + M_%matvec => precond_none + end select + where(abs(diagonal)>epsilon(zero_${s}$)) diagonal = one_${s}$/diagonal + end if + ! matvec for the operator + op%matvec => matvec + + ! direchlet boundary conditions mask + if(present(di))then + di_ => di + else + allocate(di_(n),source=.false._int8) + end if + + ! workspace for the solver + if(present(workspace)) then + workspace_ => workspace + else + allocate( workspace_ ) + end if + if(.not.allocated(workspace_%tmp)) allocate( workspace_%tmp(n,stdlib_size_wksp_bicgstab) , source = zero_${s}$ ) + !------------------------- + ! 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_bicgstab_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 + M_ => null() + workspace_ => null() + contains + + subroutine matvec(x,y,alpha,beta,op) + #:if matrix == "dense" + use stdlib_linalg_blas, only: gemv + #:endif + ${t}$, intent(in) :: x(:) + ${t}$, intent(inout) :: y(:) + ${t}$, intent(in) :: alpha + ${t}$, intent(in) :: beta + character(1), intent(in) :: op + #:if matrix == "dense" + call gemv(op,m=size(A,1),n=size(A,2),alpha=alpha,a=A,lda=size(A,1),x=x,incx=1,beta=beta,y=y,incy=1) + #:else + call spmv( A , x, y , alpha, beta , op) + #:endif + y = merge( zero_${s}$, y, di_ ) + end subroutine + + subroutine precond_none(x,y,alpha,beta,op) + ${t}$, intent(in) :: x(:) + ${t}$, intent(inout) :: y(:) + ${t}$, intent(in) :: alpha + ${t}$, intent(in) :: beta + character(1), intent(in) :: op + y = merge( zero_${s}$, x, di_ ) + end subroutine + subroutine precond_jacobi(x,y,alpha,beta,op) + ${t}$, intent(in) :: x(:) + ${t}$, intent(inout) :: y(:) + ${t}$, intent(in) :: alpha + ${t}$, intent(in) :: beta + character(1), intent(in) :: op + y = merge( zero_${s}$, diagonal * x, di_ ) ! inverted diagonal + end subroutine + end subroutine + + #:endfor + #:endfor + +end submodule stdlib_linalg_iterative_bicgstab diff --git a/src/stdlib_linalg_iterative_solvers_pcg.fypp b/src/stdlib_linalg_iterative_solvers_pcg.fypp index 312525ae6..d75ea8d6f 100644 --- a/src/stdlib_linalg_iterative_solvers_pcg.fypp +++ b/src/stdlib_linalg_iterative_solvers_pcg.fypp @@ -12,11 +12,6 @@ submodule(stdlib_linalg_iterative_solvers) stdlib_linalg_iterative_pcg use stdlib_optval, only: optval implicit none - enum, bind(c) - enumerator :: pc_none = 0 - enumerator :: pc_jacobi - end enum - contains #:for k, t, s in R_KINDS_TYPES diff --git a/test/linalg/test_linalg_solve_iterative.fypp b/test/linalg/test_linalg_solve_iterative.fypp index 2e4bd4098..8da08fcb8 100644 --- a/test/linalg/test_linalg_solve_iterative.fypp +++ b/test/linalg/test_linalg_solve_iterative.fypp @@ -22,7 +22,9 @@ module test_linalg_solve_iterative allocate(tests(0)) tests = [ new_unittest("stdlib_solve_cg",test_stdlib_solve_cg), & - new_unittest("stdlib_solve_pcg",test_stdlib_solve_pcg) ] + new_unittest("stdlib_solve_pcg",test_stdlib_solve_pcg), & + new_unittest("stdlib_solve_bicgstab",test_stdlib_solve_bicgstab), & + new_unittest("stdlib_solve_bicgstab_nonsymmetric",test_stdlib_solve_bicgstab_nonsymmetric) ] end subroutine test_linear_systems @@ -78,6 +80,122 @@ module test_linalg_solve_iterative #:endfor end subroutine test_stdlib_solve_pcg + subroutine test_stdlib_solve_bicgstab(error) + type(error_type), allocatable, intent(out) :: error + + #:for k, t, s in R_KINDS_TYPES + ! Test 1: Simple non-symmetric matrix (same as SciPy example) + block + ${t}$, parameter :: tol = 1000*epsilon(0.0_${k}$) + ${t}$ :: A(4,4) = reshape([${t}$ :: 4, 2, 0, 1, & + 3, 0, 0, 2, & + 0, 1, 1, 1, & + 0, 2, 1, 0], [4,4]) + ${t}$ :: x(4), load(4), xref(4) + + ! Reference solution computed with high precision + xref = [12.5_${k}$, -17._${k}$, 23.5_${k}$, -24.5_${k}$] + x = 0.0_${k}$ ! initial guess + load = [-1.0_${k}$, -0.5_${k}$, -1.0_${k}$, 2.0_${k}$] ! load vector + + call stdlib_solve_bicgstab(A, load, x, rtol=1.e-10_${k}$) + + call check(error, norm2(x-xref)