From 3dc40884d48517ea165acf6d381cae77833af615 Mon Sep 17 00:00:00 2001 From: jalvesz Date: Sat, 27 Sep 2025 13:20:42 +0200 Subject: [PATCH 1/8] implement bicgstab with copilot --- doc/specs/stdlib_linalg_iterative_solvers.md | 95 +++++++ example/linalg/CMakeLists.txt | 1 + example/linalg/example_solve_bicgstab.f90 | 41 +++ src/CMakeLists.txt | 1 + src/stdlib_linalg_iterative_solvers.fypp | 55 +++- ...lib_linalg_iterative_solvers_bicgstab.fypp | 252 ++++++++++++++++++ src/stdlib_linalg_iterative_solvers_pcg.fypp | 5 - test/linalg/test_linalg_solve_iterative.fypp | 120 ++++++++- 8 files changed, 563 insertions(+), 7 deletions(-) create mode 100644 example/linalg/example_solve_bicgstab.f90 create mode 100644 src/stdlib_linalg_iterative_solvers_bicgstab.fypp diff --git a/doc/specs/stdlib_linalg_iterative_solvers.md b/doc/specs/stdlib_linalg_iterative_solvers.md index 6a0b69a95..8672bbb5f 100644 --- a/doc/specs/stdlib_linalg_iterative_solvers.md +++ b/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__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_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 conditions 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 + +```fortran +{!example/linalg/example_solve_bicgstab.f90!} ``` \ No newline at end of file diff --git a/example/linalg/CMakeLists.txt b/example/linalg/CMakeLists.txt index 2aa9ca2e8..70f6d252b 100644 --- a/example/linalg/CMakeLists.txt +++ b/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) 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/src/CMakeLists.txt b/src/CMakeLists.txt index 5ff2342fc..4af2be566 100644 --- a/src/CMakeLists.txt +++ b/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) diff --git a/src/stdlib_linalg_iterative_solvers.fypp b/src/stdlib_linalg_iterative_solvers.fypp index 70e412732..4bff77dbd 100644 --- a/src/stdlib_linalg_iterative_solvers.fypp +++ b/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 !------------------------------------------------------------------ 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) Date: Sat, 27 Sep 2025 14:08:20 +0200 Subject: [PATCH 2/8] replace logical size indicator --- src/stdlib_linalg_iterative_solvers.fypp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/stdlib_linalg_iterative_solvers.fypp b/src/stdlib_linalg_iterative_solvers.fypp index 4bff77dbd..0401830d3 100644 --- a/src/stdlib_linalg_iterative_solvers.fypp +++ b/src/stdlib_linalg_iterative_solvers.fypp @@ -104,7 +104,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 @@ -168,7 +168,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 @@ -194,7 +194,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 From a3a3bdde3987fb5166adce002b45a8db78951670 Mon Sep 17 00:00:00 2001 From: jalvesz Date: Sun, 28 Sep 2025 12:04:34 +0200 Subject: [PATCH 3/8] add example with wilkilson matrix --- example/linalg/CMakeLists.txt | 1 + .../example_solve_bicgstab_wilkinson.f90 | 52 +++++++++++++++++++ 2 files changed, 53 insertions(+) create mode 100644 example/linalg/example_solve_bicgstab_wilkinson.f90 diff --git a/example/linalg/CMakeLists.txt b/example/linalg/CMakeLists.txt index 70f6d252b..134999a8b 100644 --- a/example/linalg/CMakeLists.txt +++ b/example/linalg/CMakeLists.txt @@ -42,6 +42,7 @@ 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_wilkinson.f90 b/example/linalg/example_solve_bicgstab_wilkinson.f90 new file mode 100644 index 000000000..c93f025aa --- /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) = [( dble(abs(i)), 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 From 938c0dec26e5c68c77cee39b0472df0cc81f3567 Mon Sep 17 00:00:00 2001 From: jalvesz Date: Sun, 28 Sep 2025 12:35:02 +0200 Subject: [PATCH 4/8] solve conflict with master branch --- cmake/stdlib.cmake | 65 ++++++++++++++++++- src/CMakeLists.txt | 83 ++++-------------------- src/blas/CMakeLists.txt | 41 ++++++------ src/lapack/CMakeLists.txt | 129 ++++++++++++++++++++++---------------- 4 files changed, 172 insertions(+), 146 deletions(-) diff --git a/cmake/stdlib.cmake b/cmake/stdlib.cmake index bdaf87d1f..ef71c6371 100644 --- a/cmake/stdlib.cmake +++ b/cmake/stdlib.cmake @@ -12,7 +12,8 @@ function(preprocess preproc preprocopts srcext trgext srcfiles trgfiles) set(_trgfiles) foreach(srcfile IN LISTS srcfiles) - string(REGEX REPLACE "\\.${srcext}$" ".${trgext}" trgfile ${srcfile}) + get_filename_component(filename ${srcfile} NAME) + string(REGEX REPLACE "\\.${srcext}$" ".${trgext}" trgfile ${filename}) add_custom_command( OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/${trgfile} COMMAND ${preproc} ${preprocopts} ${CMAKE_CURRENT_SOURCE_DIR}/${srcfile} ${CMAKE_CURRENT_BINARY_DIR}/${trgfile} @@ -47,3 +48,65 @@ function (fypp_f90pp fyppopts fyppfiles F90files) set(${F90files} ${_F90files} PARENT_SCOPE) endfunction() +# Helper function to configure stdlib targets +# +# It preprocesses the given fypp and fypp+cpp files, combines them with the +# regular Fortran files, and creates a library target with the given name. +# Args: +# target_name [in]: Name of the library target to create +# regular_sources_var [in]: Regular Fortran sources +# fypp_files_var [in]: Sources to be preprocessed with fypp +# cpp_files_var [in]: Sources to be preprocessed with fypp and cpp +# +function(configure_stdlib_target target_name regular_sources_var fypp_files_var cpp_files_var) + #### Pre-process: .fpp -> .f90 via Fypp + fypp_f90("${fyppFlags}" "${${fypp_files_var}}" ${target_name}_fypp_outFiles) + #### Pre-process: .fypp -> .F90 via Fypp (for C preprocessor directives) + fypp_f90pp("${fyppFlags}" "${${cpp_files_var}}" ${target_name}_cpp_outFiles) + + list(APPEND all_sources ${${target_name}_fypp_outFiles}) + list(APPEND all_sources ${${target_name}_cpp_outFiles}) + list(APPEND all_sources ${${regular_sources_var}}) + + add_library(${target_name} ${all_sources}) + add_library(${PROJECT_NAME}::${target_name} ALIAS ${target_name}) + + set_target_properties( + ${target_name} + PROPERTIES + POSITION_INDEPENDENT_CODE ON + WINDOWS_EXPORT_ALL_SYMBOLS ON + ) + + if(CMAKE_Fortran_COMPILER_ID STREQUAL GNU AND CMAKE_Fortran_COMPILER_VERSION VERSION_LESS 10.0) + target_compile_options( + ${target_name} + PRIVATE + $<$:-fno-range-check> + ) + endif() + + set(LIB_MOD_DIR ${CMAKE_CURRENT_BINARY_DIR}/mod_files/${target_name}/) + set(INSTALL_MOD_DIR "${CMAKE_INSTALL_MODULEDIR}/${target_name}") + # We need the module directory before we finish the configure stage since the + # build interface might resolve before the module directory is generated by CMake + if(NOT EXISTS "${LIB_MOD_DIR}") + file(MAKE_DIRECTORY "${LIB_MOD_DIR}") + endif() + + set_target_properties(${target_name} PROPERTIES + Fortran_MODULE_DIRECTORY ${LIB_MOD_DIR} + ) + target_include_directories(${target_name} PUBLIC + $ + $ + ) + + install(TARGETS ${target_name} + EXPORT ${PROJECT_NAME}-targets + RUNTIME DESTINATION "${CMAKE_INSTALL_BINDIR}" + ARCHIVE DESTINATION "${CMAKE_INSTALL_LIBDIR}" + LIBRARY DESTINATION "${CMAKE_INSTALL_LIBDIR}" + ) + install(DIRECTORY ${LIB_MOD_DIR} DESTINATION "${INSTALL_MOD_DIR}") +endfunction() \ No newline at end of file diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 4af2be566..6b255ae9b 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,6 +1,6 @@ -#### Pre-process: .fpp -> .f90 via Fypp +add_subdirectory(blas) +add_subdirectory(lapack) -# Create a list of the files to be preprocessed set(fppFiles stdlib_ascii.fypp stdlib_bitsets.fypp @@ -36,6 +36,10 @@ set(fppFiles stdlib_linalg_determinant.fypp 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 stdlib_linalg_norms.fypp stdlib_linalg_state.fypp @@ -89,25 +93,8 @@ set(fppFiles stdlib_strings.fypp stdlib_version.fypp ) - -# Preprocessed files to contain preprocessor directives -> .F90 -set(cppFiles - stdlib_linalg_constants.fypp - stdlib_linalg_blas.fypp - stdlib_linalg_lapack.fypp - 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) -add_subdirectory(lapack) - -fypp_f90("${fyppFlags}" "${fppFiles}" outFiles) -fypp_f90pp("${fyppFlags}" "${cppFiles}" outPreprocFiles) - -set(SRC +set(cppFiles stdlib_linalg_constants.fypp) +set(f90Files stdlib_ansi.f90 stdlib_ansi_operator.f90 stdlib_ansi_to_string.f90 @@ -128,57 +115,9 @@ set(SRC stdlib_specialfunctions_legendre.f90 stdlib_quadrature_gauss.f90 stdlib_stringlist_type.f90 - ${outFiles} - ${outPreprocFiles} + $,f18estop.f90,f08estop.f90> ) -add_library(${PROJECT_NAME} ${SRC}) - -# Link to BLAS and LAPACK -if(BLAS_FOUND AND LAPACK_FOUND) - target_link_libraries(${PROJECT_NAME} "BLAS::BLAS") - target_link_libraries(${PROJECT_NAME} "LAPACK::LAPACK") -endif() - -set_target_properties( - ${PROJECT_NAME} - PROPERTIES - POSITION_INDEPENDENT_CODE ON - WINDOWS_EXPORT_ALL_SYMBOLS ON -) +configure_stdlib_target(${PROJECT_NAME} f90Files fppFiles cppFiles) -if(CMAKE_Fortran_COMPILER_ID STREQUAL GNU AND CMAKE_Fortran_COMPILER_VERSION VERSION_LESS 10.0) - target_compile_options( - ${PROJECT_NAME} - PRIVATE - $<$:-fno-range-check> - ) -endif() - -set(LIB_MOD_DIR ${CMAKE_CURRENT_BINARY_DIR}/mod_files/) -# We need the module directory before we finish the configure stage since the -# build interface might resolve before the module directory is generated by CMake -if(NOT EXISTS "${LIB_MOD_DIR}") - file(MAKE_DIRECTORY "${LIB_MOD_DIR}") -endif() - -set_target_properties(${PROJECT_NAME} PROPERTIES - Fortran_MODULE_DIRECTORY ${LIB_MOD_DIR}) -target_include_directories(${PROJECT_NAME} PUBLIC - $ - $ -) - -if(f18errorstop) - target_sources(${PROJECT_NAME} PRIVATE f18estop.f90) -else() - target_sources(${PROJECT_NAME} PRIVATE f08estop.f90) -endif() - -install(TARGETS ${PROJECT_NAME} - EXPORT ${PROJECT_NAME}-targets - RUNTIME DESTINATION "${CMAKE_INSTALL_BINDIR}" - ARCHIVE DESTINATION "${CMAKE_INSTALL_LIBDIR}" - LIBRARY DESTINATION "${CMAKE_INSTALL_LIBDIR}" -) -install(DIRECTORY ${LIB_MOD_DIR} DESTINATION "${CMAKE_INSTALL_MODULEDIR}") +target_link_libraries(${PROJECT_NAME} PUBLIC blas lapack) diff --git a/src/blas/CMakeLists.txt b/src/blas/CMakeLists.txt index 5099b909e..c00e13776 100644 --- a/src/blas/CMakeLists.txt +++ b/src/blas/CMakeLists.txt @@ -1,20 +1,25 @@ -#### - -set(dir "${CMAKE_CURRENT_SOURCE_DIR}") - -list(APPEND fppFiles - blas/stdlib_blas_constants.fypp - blas/stdlib_blas.fypp - blas/stdlib_blas_level1.fypp - blas/stdlib_blas_level2_ban.fypp - blas/stdlib_blas_level2_gen.fypp - blas/stdlib_blas_level2_pac.fypp - blas/stdlib_blas_level2_sym.fypp - blas/stdlib_blas_level2_tri.fypp - blas/stdlib_blas_level3_gen.fypp - blas/stdlib_blas_level3_sym.fypp - blas/stdlib_blas_level3_tri.fypp - blas/stdlib_linalg_blas_aux.fypp +set(blas_fppFiles + ../stdlib_kinds.fypp + stdlib_blas_constants.fypp + stdlib_blas.fypp + stdlib_blas_level1.fypp + stdlib_blas_level2_ban.fypp + stdlib_blas_level2_gen.fypp + stdlib_blas_level2_pac.fypp + stdlib_blas_level2_sym.fypp + stdlib_blas_level2_tri.fypp + stdlib_blas_level3_gen.fypp + stdlib_blas_level3_sym.fypp + stdlib_blas_level3_tri.fypp + stdlib_linalg_blas_aux.fypp ) +set(blas_cppFiles + ../stdlib_linalg_constants.fypp + ../stdlib_linalg_blas.fypp +) + +configure_stdlib_target(blas "" blas_fppFiles blas_cppFiles) -set(fppFiles "${fppFiles}" PARENT_SCOPE) +if(BLAS_FOUND) + target_link_libraries(blas PUBLIC "BLAS::BLAS") +endif() \ No newline at end of file diff --git a/src/lapack/CMakeLists.txt b/src/lapack/CMakeLists.txt index a8f3445ce..994a9e513 100644 --- a/src/lapack/CMakeLists.txt +++ b/src/lapack/CMakeLists.txt @@ -1,58 +1,77 @@ -list(APPEND fppFiles - lapack/stdlib_lapack_base.fypp - lapack/stdlib_lapack_solve.fypp - lapack/stdlib_lapack_others.fypp - lapack/stdlib_lapack_orthogonal_factors.fypp - lapack/stdlib_lapack_eig_svd_lsq.fypp - lapack/stdlib_linalg_lapack_aux.fypp - lapack/stdlib_lapack_auxiliary.fypp - lapack/stdlib_lapack_blas_like_base.fypp - lapack/stdlib_lapack_blas_like_l1.fypp - lapack/stdlib_lapack_blas_like_l2.fypp - lapack/stdlib_lapack_blas_like_l3.fypp - lapack/stdlib_lapack_blas_like_mnorm.fypp - lapack/stdlib_lapack_blas_like_scalar.fypp - lapack/stdlib_lapack_cosine_sine.fypp - lapack/stdlib_lapack_cosine_sine2.fypp - lapack/stdlib_lapack_eigv_comp.fypp - lapack/stdlib_lapack_eigv_comp2.fypp - lapack/stdlib_lapack_eigv_gen.fypp - lapack/stdlib_lapack_eigv_gen2.fypp - lapack/stdlib_lapack_eigv_gen3.fypp - lapack/stdlib_lapack_eigv_std_driver.fypp - lapack/stdlib_lapack_eigv_svd_bidiag_dc.fypp - lapack/stdlib_lapack_eigv_svd_drivers.fypp - lapack/stdlib_lapack_eigv_svd_drivers2.fypp - lapack/stdlib_lapack_eigv_svd_drivers3.fypp - lapack/stdlib_lapack_eigv_sym_comp.fypp - lapack/stdlib_lapack_eigv_sym.fypp - lapack/stdlib_lapack_eigv_tridiag.fypp - lapack/stdlib_lapack_eigv_tridiag2.fypp - lapack/stdlib_lapack_eigv_tridiag3.fypp - lapack/stdlib_lapack_givens_jacobi_rot.fypp - lapack/stdlib_lapack_householder_reflectors.fypp - lapack/stdlib_lapack_lsq.fypp - lapack/stdlib_lapack_lsq_aux.fypp - lapack/stdlib_lapack_lsq_constrained.fypp - lapack/stdlib_lapack_orthogonal_factors_ql.fypp - lapack/stdlib_lapack_orthogonal_factors_qr.fypp - lapack/stdlib_lapack_orthogonal_factors_rz.fypp - lapack/stdlib_lapack_others_sm.fypp - lapack/stdlib_lapack_solve_aux.fypp - lapack/stdlib_lapack_solve_chol_comp.fypp - lapack/stdlib_lapack_solve_chol.fypp - lapack/stdlib_lapack_solve_ldl_comp.fypp - lapack/stdlib_lapack_solve_ldl_comp2.fypp - lapack/stdlib_lapack_solve_ldl_comp3.fypp - lapack/stdlib_lapack_solve_ldl_comp4.fypp - lapack/stdlib_lapack_solve_ldl.fypp - lapack/stdlib_lapack_solve_lu_comp.fypp - lapack/stdlib_lapack_solve_lu.fypp - lapack/stdlib_lapack_solve_tri_comp.fypp - lapack/stdlib_lapack_svd_bidiag_qr.fypp - lapack/stdlib_lapack_svd_comp.fypp - lapack/stdlib_lapack_svd_comp2.fypp +set(lapack_fppFiles + ../stdlib_kinds.fypp + ../stdlib_linalg_state.fypp + ../stdlib_error.fypp + ../stdlib_optval.fypp + ../stdlib_io.fypp + ../stdlib_ascii.fypp + ../stdlib_string_type.fypp + stdlib_lapack_base.fypp + stdlib_lapack_solve.fypp + stdlib_lapack_others.fypp + stdlib_lapack_orthogonal_factors.fypp + stdlib_lapack_eig_svd_lsq.fypp + stdlib_linalg_lapack_aux.fypp + stdlib_lapack_auxiliary.fypp + stdlib_lapack_blas_like_base.fypp + stdlib_lapack_blas_like_l1.fypp + stdlib_lapack_blas_like_l2.fypp + stdlib_lapack_blas_like_l3.fypp + stdlib_lapack_blas_like_mnorm.fypp + stdlib_lapack_blas_like_scalar.fypp + stdlib_lapack_cosine_sine.fypp + stdlib_lapack_cosine_sine2.fypp + stdlib_lapack_eigv_comp.fypp + stdlib_lapack_eigv_comp2.fypp + stdlib_lapack_eigv_gen.fypp + stdlib_lapack_eigv_gen2.fypp + stdlib_lapack_eigv_gen3.fypp + stdlib_lapack_eigv_std_driver.fypp + stdlib_lapack_eigv_svd_bidiag_dc.fypp + stdlib_lapack_eigv_svd_drivers.fypp + stdlib_lapack_eigv_svd_drivers2.fypp + stdlib_lapack_eigv_svd_drivers3.fypp + stdlib_lapack_eigv_sym_comp.fypp + stdlib_lapack_eigv_sym.fypp + stdlib_lapack_eigv_tridiag.fypp + stdlib_lapack_eigv_tridiag2.fypp + stdlib_lapack_eigv_tridiag3.fypp + stdlib_lapack_givens_jacobi_rot.fypp + stdlib_lapack_householder_reflectors.fypp + stdlib_lapack_lsq.fypp + stdlib_lapack_lsq_aux.fypp + stdlib_lapack_lsq_constrained.fypp + stdlib_lapack_orthogonal_factors_ql.fypp + stdlib_lapack_orthogonal_factors_qr.fypp + stdlib_lapack_orthogonal_factors_rz.fypp + stdlib_lapack_others_sm.fypp + stdlib_lapack_solve_aux.fypp + stdlib_lapack_solve_chol_comp.fypp + stdlib_lapack_solve_chol.fypp + stdlib_lapack_solve_ldl_comp.fypp + stdlib_lapack_solve_ldl_comp2.fypp + stdlib_lapack_solve_ldl_comp3.fypp + stdlib_lapack_solve_ldl_comp4.fypp + stdlib_lapack_solve_ldl.fypp + stdlib_lapack_solve_lu_comp.fypp + stdlib_lapack_solve_lu.fypp + stdlib_lapack_solve_tri_comp.fypp + stdlib_lapack_svd_bidiag_qr.fypp + stdlib_lapack_svd_comp.fypp + stdlib_lapack_svd_comp2.fypp ) +set(lapack_cppFiles + ../stdlib_linalg_constants.fypp + ../stdlib_linalg_lapack.fypp +) +set(lapack_f90Files + $,../f18estop.f90,../f08estop.f90> +) + +configure_stdlib_target(lapack lapack_f90Files lapack_fppFiles lapack_cppFiles) -set(fppFiles "${fppFiles}" PARENT_SCOPE) +if(LAPACK_FOUND) + target_link_libraries(lapack PUBLIC "LAPACK::LAPACK") +endif() +target_link_libraries(lapack PUBLIC blas) \ No newline at end of file From 2d08a292a7ef18f99649d94809c9596962a35521 Mon Sep 17 00:00:00 2001 From: jalvesz Date: Sun, 28 Sep 2025 12:45:50 +0200 Subject: [PATCH 5/8] revert space --- cmake/stdlib.cmake | 1 - 1 file changed, 1 deletion(-) diff --git a/cmake/stdlib.cmake b/cmake/stdlib.cmake index e47805db7..0ec86e299 100644 --- a/cmake/stdlib.cmake +++ b/cmake/stdlib.cmake @@ -109,5 +109,4 @@ function(configure_stdlib_target target_name regular_sources_var fypp_files_var LIBRARY DESTINATION "${CMAKE_INSTALL_LIBDIR}" ) install(DIRECTORY ${LIB_MOD_DIR} DESTINATION "${INSTALL_MOD_DIR}") - endfunction() From 56345f9fccb1f9ae03ae5b50f16d7fbc000a4b07 Mon Sep 17 00:00:00 2001 From: jalvesz Date: Sun, 28 Sep 2025 13:55:29 +0200 Subject: [PATCH 6/8] fix documentation links and descriptions --- doc/specs/stdlib_linalg_iterative_solvers.md | 31 ++++++++++++-------- src/stdlib_linalg_iterative_solvers.fypp | 15 +++++++++- 2 files changed, 32 insertions(+), 14 deletions(-) diff --git a/doc/specs/stdlib_linalg_iterative_solvers.md b/doc/specs/stdlib_linalg_iterative_solvers.md index 8672bbb5f..3f68cef22 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 @@ -259,7 +259,7 @@ Implements the Biconjugate Gradient Stabilized (BiCGSTAB) method for solving the #### Syntax -`call ` [[stdlib_iterative_solvers(module):stdlib_solve_bicgstab_kernel(interface)]] ` (A, M, b, x, rtol, atol, maxiter, workspace)` +`call ` [[stdlib_linalg_iterative_solvers(module):stdlib_solve_bicgstab_kernel(interface)]] ` (A, M, b, x, rtol, atol, maxiter, workspace)` #### Status @@ -298,7 +298,7 @@ Provides a user-friendly interface to the BiCGSTAB method for solving \( Ax = b #### Syntax -`call ` [[stdlib_iterative_solvers(module):stdlib_solve_bicgstab(interface)]] ` (A, b, x [, di, rtol, atol, maxiter, restart, precond, M, workspace])` +`call ` [[stdlib_linalg_iterative_solvers(module):stdlib_solve_bicgstab(interface)]] ` (A, b, x [, di, rtol, atol, maxiter, restart, precond, M, workspace])` #### Status @@ -339,8 +339,13 @@ BiCGSTAB is particularly effective for: The method uses 8 auxiliary vectors internally, requiring more memory than simpler methods but often providing better stability and convergence properties. -#### Example +#### 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/src/stdlib_linalg_iterative_solvers.fypp b/src/stdlib_linalg_iterative_solvers.fypp index 0401830d3..1b216a588 100644 --- a/src/stdlib_linalg_iterative_solvers.fypp +++ b/src/stdlib_linalg_iterative_solvers.fypp @@ -27,6 +27,8 @@ module stdlib_linalg_iterative_solvers !! 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 @@ -38,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 @@ -89,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 @@ -154,6 +161,9 @@ module stdlib_linalg_iterative_solvers 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 @@ -180,6 +190,9 @@ 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 From ef5611a8a136567d51e1d065cbf3cec2a1c8a8da Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jos=C3=A9=20Alves?= <102541118+jalvesz@users.noreply.github.com> Date: Mon, 29 Sep 2025 22:10:04 +0200 Subject: [PATCH 7/8] Update example/linalg/example_solve_bicgstab_wilkinson.f90 Co-authored-by: Jeremie Vandenplas --- example/linalg/example_solve_bicgstab_wilkinson.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/example/linalg/example_solve_bicgstab_wilkinson.f90 b/example/linalg/example_solve_bicgstab_wilkinson.f90 index c93f025aa..7e1c54ddc 100644 --- a/example/linalg/example_solve_bicgstab_wilkinson.f90 +++ b/example/linalg/example_solve_bicgstab_wilkinson.f90 @@ -15,7 +15,7 @@ program example_solve_bicgstab_wilkinson ! 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) = [( dble(abs(i)), i=-(n-1)/2, (n-1)/2)] + 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] From 8a5f861d0427df25f1d3e7ba2bdf8d546c061985 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jos=C3=A9=20Alves?= <102541118+jalvesz@users.noreply.github.com> Date: Sun, 5 Oct 2025 18:43:36 +0200 Subject: [PATCH 8/8] Update doc/specs/stdlib_linalg_iterative_solvers.md Co-authored-by: Jeremie Vandenplas --- doc/specs/stdlib_linalg_iterative_solvers.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/specs/stdlib_linalg_iterative_solvers.md b/doc/specs/stdlib_linalg_iterative_solvers.md index 3f68cef22..949f55a44 100644 --- a/doc/specs/stdlib_linalg_iterative_solvers.md +++ b/doc/specs/stdlib_linalg_iterative_solvers.md @@ -316,7 +316,7 @@ Subroutine `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 conditions values should be stored in the `b` load array. This argument is `intent(in)`. +`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)`.