Skip to content

Navigation Menu

Sign in
Appearance settings

Search code, repositories, users, issues, pull requests...

Provide feedback

We read every piece of feedback, and take your input very seriously.

Saved searches

Use saved searches to filter your results more quickly

Sign up
Appearance settings

Commit eb46153

Browse files
jalveszjvdp1
andauthored
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 <jeremie.vandenplas@gmail.com> * Update doc/specs/stdlib_linalg_iterative_solvers.md Co-authored-by: Jeremie Vandenplas <jeremie.vandenplas@gmail.com> --------- Co-authored-by: Jeremie Vandenplas <jeremie.vandenplas@gmail.com>
1 parent 80a8164 commit eb46153

9 files changed

+647
-20
lines changed

‎doc/specs/stdlib_linalg_iterative_solvers.md‎

Lines changed: 110 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -10,12 +10,12 @@ title: linalg_iterative_solvers
1010

1111
The `stdlib_linalg_iterative_solvers` module provides base implementations for known iterative solver methods. Each method is exposed with two procedure flavors:
1212

13-
* 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.
13+
* 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.
1414

15-
* A `solve_<method>` which proposes an off-the-shelf ready to use interface for `dense` and `CSR_<kind>_type` matrices for all `real` kinds.
15+
* 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.
1616

1717
<!-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -->
18-
### The `linop` derived type
18+
### The `stdlib_linop` derived type
1919

2020
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.
2121

@@ -25,11 +25,11 @@ The following type-bound procedure pointers enable customization of the solver:
2525

2626
##### `matvec`
2727

28-
Proxy procedure for the matrix-vector product $y = alpha * op(M) * x + beta * y$.
28+
Proxy procedure for the matrix-vector product \(y = alpha * op(M) * x + beta * y\).
2929

3030
#### Syntax
3131

32-
`call ` [[stdlib_iterative_solvers(module):matvec(interface)]] `(x,y,alpha,beta,op)`
32+
`call ` [[stdlib_linalg_iterative_solvers(module):stdlib_linop_dp_type(type)]] `%matvec(x,y,alpha,beta,op)`
3333

3434
###### Class
3535

@@ -53,7 +53,7 @@ Proxy procedure for the `dot_product`.
5353

5454
#### Syntax
5555

56-
`res = ` [[stdlib_iterative_solvers(module):inner_product(interface)]] `(x,y)`
56+
`res = ` [[stdlib_linalg_iterative_solvers(module):stdlib_linop_dp_type(type)]] `%inner_product(x,y)`
5757

5858
###### Class
5959

@@ -99,7 +99,7 @@ Implements the Conjugate Gradient (CG) method for solving the linear system \( A
9999

100100
#### Syntax
101101

102-
`call ` [[stdlib_iterative_solvers(module):stdlib_solve_cg_kernel(interface)]] ` (A, b, x, tol, maxiter, workspace)`
102+
`call ` [[stdlib_linalg_iterative_solvers(module):stdlib_solve_cg_kernel(interface)]] ` (A, b, x, tol, maxiter, workspace)`
103103

104104
#### Status
105105

@@ -132,7 +132,7 @@ Provides a user-friendly interface to the CG method for solving \( Ax = b \), su
132132

133133
#### Syntax
134134

135-
`call ` [[stdlib_iterative_solvers(module):solve_cg(interface)]] ` (A, b, x [, di, rtol, atol, maxiter, restart, workspace])`
135+
`call ` [[stdlib_linalg_iterative_solvers(module):stdlib_solve_cg(interface)]] ` (A, b, x [, di, rtol, atol, maxiter, restart, workspace])`
136136

137137
#### Status
138138

@@ -173,7 +173,7 @@ Implements the Preconditioned Conjugate Gradient (PCG) method for solving the li
173173

174174
#### Syntax
175175

176-
`call ` [[stdlib_iterative_solvers(module):stdlib_solve_cg_kernel(interface)]] ` (A, M, b, x, tol, maxiter, workspace)`
176+
`call ` [[stdlib_linalg_iterative_solvers(module):stdlib_solve_cg_kernel(interface)]] ` (A, M, b, x, tol, maxiter, workspace)`
177177

178178
#### Status
179179

@@ -214,7 +214,7 @@ Provides a user-friendly interface to the PCG method for solving \( Ax = b \), s
214214

215215
#### Syntax
216216

217-
`call ` [[stdlib_iterative_solvers(module):stdlib_solve_pcg(interface)]] ` (A, b, x [, di, tol, maxiter, restart, precond, M, workspace])`
217+
`call ` [[stdlib_linalg_iterative_solvers(module):stdlib_solve_pcg(interface)]] ` (A, b, x [, di, tol, maxiter, restart, precond, M, workspace])`
218218

219219
#### Status
220220

@@ -248,4 +248,104 @@ Subroutine
248248

249249
```fortran
250250
{!example/linalg/example_solve_pcg.f90!}
251+
```
252+
253+
<!-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -->
254+
### `stdlib_solve_bicgstab_kernel` subroutine
255+
256+
#### Description
257+
258+
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.
259+
260+
#### Syntax
261+
262+
`call ` [[stdlib_linalg_iterative_solvers(module):stdlib_solve_bicgstab_kernel(interface)]] ` (A, M, b, x, rtol, atol, maxiter, workspace)`
263+
264+
#### Status
265+
266+
Experimental
267+
268+
#### Class
269+
270+
Subroutine
271+
272+
#### Argument(s)
273+
274+
`A`: `class(stdlib_linop_<kind>_type)` defining the linear operator. This argument is `intent(in)`.
275+
276+
`M`: `class(stdlib_linop_<kind>_type)` defining the preconditioner linear operator. This argument is `intent(in)`.
277+
278+
`b`: 1-D array of `real(<kind>)` defining the loading conditions of the linear system. This argument is `intent(in)`.
279+
280+
`x`: 1-D array of `real(<kind>)` which serves as the input initial guess and the output solution. This argument is `intent(inout)`.
281+
282+
`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)`.
283+
284+
`maxiter`: scalar of type `integer` defining the maximum allowed number of iterations. This argument is `intent(in)`.
285+
286+
`workspace`: scalar derived type of `type(stdlib_solver_workspace_<kind>_type)` holding the work array for the solver. This argument is `intent(inout)`.
287+
288+
#### Note
289+
290+
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.
291+
292+
<!-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -->
293+
### `stdlib_solve_bicgstab` subroutine
294+
295+
#### Description
296+
297+
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.
298+
299+
#### Syntax
300+
301+
`call ` [[stdlib_linalg_iterative_solvers(module):stdlib_solve_bicgstab(interface)]] ` (A, b, x [, di, rtol, atol, maxiter, restart, precond, M, workspace])`
302+
303+
#### Status
304+
305+
Experimental
306+
307+
#### Class
308+
309+
Subroutine
310+
311+
#### Argument(s)
312+
313+
`A`: `dense` or `CSR_<kind>_type` matrix defining the linear system. This argument is `intent(in)`.
314+
315+
`b`: 1-D array of `real(<kind>)` defining the loading conditions of the linear system. This argument is `intent(in)`.
316+
317+
`x`: 1-D array of `real(<kind>)` which serves as the input initial guess and the output solution. This argument is `intent(inout)`.
318+
319+
`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)`.
320+
321+
`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)`.
322+
323+
`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)`.
324+
325+
`restart` (optional): scalar of type `logical` indicating whether to restart the iteration with zero initial guess. Default is `.true.`. This argument is `intent(in)`.
326+
327+
`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)`.
328+
329+
`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)`.
330+
331+
`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)`.
332+
333+
#### Note
334+
335+
BiCGSTAB is particularly effective for:
336+
- Non-symmetric linear systems
337+
- Systems where CG cannot be applied
338+
- Cases where BiCG suffers from irregular convergence
339+
340+
The method uses 8 auxiliary vectors internally, requiring more memory than simpler methods but often providing better stability and convergence properties.
341+
342+
#### Example 1
343+
344+
```fortran
345+
{!example/linalg/example_solve_bicgstab.f90!}
346+
```
347+
348+
#### Example 2
349+
```fortran
350+
{!example/linalg/example_solve_bicgstab_wilkinson.f90!}
251351
```

‎example/linalg/CMakeLists.txt‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -41,6 +41,8 @@ ADD_EXAMPLE(get_norm)
4141
ADD_EXAMPLE(solve1)
4242
ADD_EXAMPLE(solve2)
4343
ADD_EXAMPLE(solve3)
44+
ADD_EXAMPLE(solve_bicgstab)
45+
ADD_EXAMPLE(solve_bicgstab_wilkinson)
4446
ADD_EXAMPLE(solve_cg)
4547
ADD_EXAMPLE(solve_pcg)
4648
ADD_EXAMPLE(solve_custom)
Lines changed: 41 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,41 @@
1+
program example_solve_bicgstab
2+
use stdlib_kinds, only: dp
3+
use stdlib_linalg_iterative_solvers
4+
implicit none
5+
6+
integer, parameter :: n = 4
7+
real(dp) :: A(n,n), b(n), x(n)
8+
integer :: i
9+
10+
! Example matrix (same as SciPy documentation)
11+
A = reshape([4.0_dp, 2.0_dp, 0.0_dp, 1.0_dp, &
12+
3.0_dp, 0.0_dp, 0.0_dp, 2.0_dp, &
13+
0.0_dp, 1.0_dp, 1.0_dp, 1.0_dp, &
14+
0.0_dp, 2.0_dp, 1.0_dp, 0.0_dp], [n,n])
15+
16+
b = [-1.0_dp, -0.5_dp, -1.0_dp, 2.0_dp]
17+
18+
x = 0.0_dp ! Initial guess
19+
20+
print *, 'Solving Ax = b using BiCGSTAB method:'
21+
print *, 'Matrix A:'
22+
do i = 1, n
23+
print '(4F8.2)', A(i,:)
24+
end do
25+
print *, 'Right-hand side b:'
26+
print '(4F8.2)', b
27+
28+
! Solve using BiCGSTAB
29+
call stdlib_solve_bicgstab(A, b, x, rtol=1e-10_dp, atol=1e-12_dp)
30+
31+
print *, 'Solution x:'
32+
print '(4F10.6)', x
33+
34+
! Verify solution
35+
print *, 'Verification A*x:'
36+
print '(4F10.6)', matmul(A, x)
37+
38+
print *, 'Residual ||b - A*x||:'
39+
print *, norm2(b - matmul(A, x))
40+
41+
end program
Lines changed: 52 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,52 @@
1+
program example_solve_bicgstab_wilkinson
2+
use stdlib_kinds, only: dp
3+
use stdlib_linalg_iterative_solvers
4+
use stdlib_sparse
5+
use stdlib_sparse_spmv
6+
implicit none
7+
8+
integer, parameter :: n = 21
9+
type(COO_dp_type) :: COO
10+
type(CSR_dp_type) :: A
11+
type(stdlib_solver_workspace_dp_type) :: workspace
12+
real(dp) :: b(n), x(n), norm_sq0
13+
integer :: i
14+
15+
! Construct the Wilkinson's matrix in COO format
16+
! https://en.wikipedia.org/wiki/Wilkinson_matrix
17+
call COO%malloc(n, n, n + 2*(n-1))
18+
COO%data(1:n) = [( real(abs(i), kind=dp), i=-(n-1)/2, (n-1)/2)]
19+
COO%data(n+1:) = 1.0_dp
20+
do i = 1, n
21+
COO%index(1:2, i) = [i,i]
22+
end do
23+
do i = 1, n-1
24+
COO%index(1:2, n+i) = [i,i+1]
25+
COO%index(1:2, n+n-1+i) = [i+1,i]
26+
end do
27+
call coo2ordered(COO,sort_data=.true.)
28+
29+
! Convert COO to CSR format
30+
call coo2csr(COO, A)
31+
32+
! Set up the right-hand side for the solution to be ones
33+
b = 0.0_dp
34+
x = 1.0_dp
35+
call spmv(A, x, b) ! b = A * 1
36+
x = 0.0_dp ! initial guess
37+
38+
! Solve the system using BiCGSTAB
39+
workspace%callback => my_logger
40+
call stdlib_solve_bicgstab(A, b, x, rtol=1e-12_dp, maxiter=100, workspace=workspace)
41+
42+
contains
43+
44+
subroutine my_logger(x,norm_sq,iter)
45+
real(dp), intent(in) :: x(:)
46+
real(dp), intent(in) :: norm_sq
47+
integer, intent(in) :: iter
48+
if(iter == 0) norm_sq0 = norm_sq
49+
print *, "Iteration: ", iter, " Residual: ", norm_sq, " Relative: ", norm_sq/norm_sq0
50+
end subroutine
51+
52+
end program example_solve_bicgstab_wilkinson

‎src/CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -37,6 +37,7 @@ set(fppFiles
3737
stdlib_linalg_qr.fypp
3838
stdlib_linalg_inverse.fypp
3939
stdlib_linalg_iterative_solvers.fypp
40+
stdlib_linalg_iterative_solvers_bicgstab.fypp
4041
stdlib_linalg_iterative_solvers_cg.fypp
4142
stdlib_linalg_iterative_solvers_pcg.fypp
4243
stdlib_linalg_pinv.fypp

0 commit comments

Comments
(0)

AltStyle によって変換されたページ (->オリジナル) /