Skip to content

Commit 463259b

Browse files
jalveszjalvesz
committed
start working on the documentation
1 parent 79dbbbe commit 463259b

File tree

3 files changed

+301
-35
lines changed

3 files changed

+301
-35
lines changed

‎doc/specs/stdlib_linalg_iterative_solvers.md

Lines changed: 249 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,249 @@
1+
---
2+
title: linalg_iterative_solvers
3+
---
4+
5+
# The `stdlib_linalg_iterative_solvers` module
6+
7+
[TOC]
8+
9+
## Introduction
10+
11+
The `stdlib_linalg_iterative_solvers` module provides base implementations for known iterative solver methods. Each method is exposed with two procedure flavors:
12+
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.
14+
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.
16+
17+
<!-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -->
18+
### The `linop` derived type
19+
20+
The `linop_<kind>` derive type is an auxiliary class enabling to abstract the definition of the linear system and the actual implementation of the solvers.
21+
22+
#### Type-bound procedures
23+
24+
The following type-bound procedures pointer enable customization of the solver:
25+
26+
##### `apply`
27+
28+
Proxy procedure for the matrix-vector product $y = alpha * M * x + beta * y$.
29+
30+
#### Syntax
31+
32+
`call ` [[stdlib_iterative_solvers(module):apply(interface)]] ` (x,y,alpha,beta)`
33+
34+
###### Class
35+
36+
Subroutine
37+
38+
###### Argument(s)
39+
40+
`x`: 1-D array of `real(<kind>)`. This argument is `intent(in)`.
41+
42+
`y`: 1-D array of `real(<kind>)`. This argument is `intent(inout)`.
43+
44+
`alpha`: scalar of `real(<kind>)`. This argument is `intent(in)`.
45+
46+
`beta`: scalar of `real(<kind>)`. This argument is `intent(in)`.
47+
48+
##### `inner_product`
49+
50+
Proxy procedure for the `dot_product`.
51+
52+
#### Syntax
53+
54+
`res = ` [[stdlib_iterative_solvers(module):inner_product(interface)]] ` (x,y)`
55+
56+
###### Class
57+
58+
Function
59+
60+
###### Argument(s)
61+
62+
`x`: 1-D array of `real(<kind>)`. This argument is `intent(in)`.
63+
64+
`y`: 1-D array of `real(<kind>)`. This argument is `intent(in)`.
65+
66+
###### Output value or Result value
67+
68+
The output is a scalar of `type` and `kind` same as to that of `x` and `y`.
69+
70+
<!-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -->
71+
### The `solver_workspace` derived type
72+
73+
The `solver_workspace_<kind>` derive type is an auxiliary class enabling to hold the data associated to the working arrays needed by the solvers to operate.
74+
75+
#### Type-bound procedures
76+
77+
- `callback`: null pointer procedure enabling to pass a callback at each iteration to check on the solvers status.
78+
79+
##### Class
80+
81+
Subroutine
82+
83+
##### Argument(s)
84+
85+
`x`: 1-D array of `real(<kind>)` type with the current state of the solution vector. This argument is `intent(in)` as it should not be modified by the callback.
86+
87+
`norm_sq`: scalar of `real(<kind>)` type representing the squared norm of the residual at the current iteration. This argument is `intent(in)`.
88+
89+
`iter`: scalar of `integer` type giving the current iteration counter. This argument is `intent(in)`.
90+
91+
<!-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -->
92+
### `solve_cg_kernel` subroutine
93+
94+
#### Description
95+
96+
Implements the Conjugate Gradient (CG) method for solving the linear system \( Ax = b \), where \( A \) is a symmetric positive-definite linear operator defined via the `linop` type. This is the core implementation, allowing flexibility for custom matrix types or parallel environments.
97+
98+
#### Syntax
99+
100+
`call ` [[stdlib_iterative_solvers(module):solve_cg_kernel(interface)]] ` (A, b, x, tol, maxiter, workspace)`
101+
102+
#### Status
103+
104+
Experimental
105+
106+
#### Class
107+
108+
Subroutine
109+
110+
#### Argument(s)
111+
112+
`A`: `class(linop_<kind>)` defining the linear operator. This argument is `intent(in)`.
113+
114+
`b`: 1-D array of `real(<kind>)` defining the loading conditions of the linear system. This argument is `intent(in)`.
115+
116+
`x`: 1-D array of `real(<kind>)` which serves as the input initial guess and the output solution. This argument is `intent(inout)`.
117+
118+
`tol`: scalar of type `real(<kind>)` specifying the requested tolerance with respect to the relative residual vector norm. This argument is `intent(in)`.
119+
120+
`maxiter`: scalar of type `integer` defining the maximum allowed number of iterations. This argument is `intent(in)`.
121+
122+
`workspace`: `type(solver_workspace_<kind>)` holding the work temporal array for the solver. This argument is `intent(inout)`.
123+
124+
<!-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -->
125+
### `solve_cg` subroutine
126+
127+
#### Description
128+
129+
Provides a user-friendly interface to the CG method for solving \( Ax = b \), supporting `dense` and `CSR_<kind>_type` matrices. It handles workspace allocation and optional parameters for customization.
130+
131+
#### Syntax
132+
133+
`call ` [[stdlib_iterative_solvers(module):solve_cg(interface)]] ` (A, b, x [, di, tol, maxiter, restart, workspace])`
134+
135+
#### Status
136+
137+
Experimental
138+
139+
#### Class
140+
141+
Subroutine
142+
143+
#### Argument(s)
144+
145+
`A`: `dense` or `CSR_<kind>_type` matrix defining the linear system. This argument is `intent(in)`.
146+
147+
`b`: 1-D array of `real(<kind>)` defining the loading conditions of the linear system. This argument is `intent(in)`.
148+
149+
`x`: 1-D array of `real(<kind>)` which serves as the input initial guess and the output solution. This argument is `intent(inout)`.
150+
151+
`di` (optional): 1-D mask array of type `logical(1)` 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)`.
152+
153+
`tol` (optional): scalar of type `real(<kind>)` specifying the requested tolerance with respect to the relative residual vector norm. If no value is given, a default value of `1.e-4` is set. This argument is `intent(in)`.
154+
155+
`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)`.
156+
157+
`workspace` (optional): `type(solver_workspace_<kind>)` 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)`.
158+
159+
#### Example
160+
161+
```fortran
162+
{!example/linalg/example_solve_cg.f90!}
163+
```
164+
165+
<!-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -->
166+
### `solve_pcg_kernel` subroutine
167+
168+
#### Description
169+
170+
Implements the Preconditionned Conjugate Gradient (PCG) method for solving the linear system \( Ax = b \), where \( A \) is a symmetric positive-definite linear operator defined via the `linop` type. This is the core implementation, allowing flexibility for custom matrix types or parallel environments.
171+
172+
#### Syntax
173+
174+
`call ` [[stdlib_iterative_solvers(module):solve_cg_kernel(interface)]] ` (A, M, b, x, tol, maxiter, workspace)`
175+
176+
#### Status
177+
178+
Experimental
179+
180+
#### Class
181+
182+
Subroutine
183+
184+
#### Argument(s)
185+
186+
`A`: `class(linop_<kind>)` defining the linear operator. This argument is `intent(in)`.
187+
188+
`M`: `class(linop_<kind>)` defining the preconditionner linear operator. This argument is `intent(in)`.
189+
190+
`b`: 1-D array of `real(<kind>)` defining the loading conditions of the linear system. This argument is `intent(in)`.
191+
192+
`x`: 1-D array of `real(<kind>)` which serves as the input initial guess and the output solution. This argument is `intent(inout)`.
193+
194+
`tol`: scalar of type `real(<kind>)` specifying the requested tolerance with respect to the relative residual vector norm. This argument is `intent(in)`.
195+
196+
`maxiter`: scalar of type `integer` defining the maximum allowed number of iterations. This argument is `intent(in)`.
197+
198+
`workspace`: `type(solver_workspace_<kind>)` holding the work temporal array for the solver. This argument is `intent(inout)`.
199+
200+
#### Example
201+
202+
```fortran
203+
{!example/linalg/example_solve_custom.f90!}
204+
```
205+
206+
<!-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -->
207+
### `solve_pcg` subroutine
208+
209+
#### Description
210+
211+
Provides a user-friendly interface to the PCG method for solving \( Ax = b \), supporting `dense` and `CSR_<kind>_type` matrices. It supports optional preconditioners and handles workspace allocation.
212+
213+
#### Syntax
214+
215+
`call ` [[stdlib_iterative_solvers(module):solve_pcg(interface)]] ` (A, b, x [, di, tol, maxiter, restart, precond, M, workspace])`
216+
217+
#### Status
218+
219+
Experimental
220+
221+
#### Class
222+
223+
Subroutine
224+
225+
#### Argument(s)
226+
227+
`A`: `dense` or `CSR_<kind>_type` matrix defining the linear system. This argument is `intent(in)`.
228+
229+
`b`: 1-D array of `real(<kind>)` defining the loading conditions of the linear system. This argument is `intent(in)`.
230+
231+
`x`: 1-D array of `real(<kind>)` which serves as the input initial guess and the output solution. This argument is `intent(inout)`.
232+
233+
`di` (optional): 1-D mask array of type `logical(1)` 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)`.
234+
235+
`tol` (optional): scalar of type `real(<kind>)` specifying the requested tolerance with respect to the relative residual vector norm. If no value is given, a default value of `1.e-4` is set. This argument is `intent(in)`.
236+
237+
`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)`.
238+
239+
`precond` (optional): scalar of type `integer` enabling to switch among the default preconditionners available. If no value is given, no preconditionning will be applied. This argument is `intent(in)`.
240+
241+
`M` (optional): `class(linop_<kind>)` defining a custom preconditionner linear operator. If given, `precond` will have no effect, a pointer is set to this custom preconditionner.
242+
243+
`workspace` (optional): `type(solver_workspace_<kind>)` 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)`.
244+
245+
#### Example
246+
247+
```fortran
248+
{!example/linalg/example_solve_pcg.f90!}
249+
```

‎src/stdlib_linalg_iterative_solvers.fypp

Lines changed: 49 additions & 32 deletions
Original file line numberDiff line numberDiff line change
@@ -3,7 +3,8 @@
33
#:set C_KINDS_TYPES = list(zip(CMPLX_KINDS, CMPLX_TYPES, CMPLX_SUFFIX))
44
#:set MATRIX_TYPES = ["dense", "CSR"]
55
#:set RANKS = range(1, 2+1)
6-
6+
!! The `stdlib_linalg_iterative_solvers` module provides interfaces for iterative solvers.
7+
!!
78
module stdlib_linalg_iterative_solvers
89
use stdlib_kinds
910
use stdlib_sparse
@@ -18,14 +19,20 @@ module stdlib_linalg_iterative_solvers
1819
end enum
1920
public :: size_wksp_cg, size_wksp_pcg
2021

21-
!> linear operator class for the iterative solvers
22+
!! version: experimental
23+
!!
24+
!! linop type holding the linear operator and its associated methods.
25+
!! The `linop` type is used to define the linear operator for the iterative solvers.
2226
#:for k, t, s in R_KINDS_TYPES
2327
type, public :: linop_${s}$
2428
procedure(vector_sub_${s}$), nopass, pointer :: apply => null()
2529
procedure(reduction_sub_${s}$), nopass, pointer :: inner_product => default_dot_${s}$
2630
end type
2731
#:endfor
2832

33+
!! version: experimental
34+
!!
35+
!! solver_workspace type holding temporal array data for the iterative solvers.
2936
#:for k, t, s in R_KINDS_TYPES
3037
type, public :: solver_workspace_${s}$
3138
${t}$, allocatable :: tmp(:,:)
@@ -57,15 +64,19 @@ module stdlib_linalg_iterative_solvers
5764
#:endfor
5865
end interface
5966

67+
!! version: experimental
68+
!!
69+
!! solve_cg_kernel interface for the conjugate gradient method.
70+
!! [Specifications](../page/specs/stdlib_linalg_iterative_solvers.html#solve_cg_kernel)
6071
interface solve_cg_kernel
6172
#:for k, t, s in R_KINDS_TYPES
6273
module subroutine solve_cg_kernel_${s}$(A,b,x,tol,maxiter,workspace)
63-
class(linop_${s}$), intent(in) :: A
64-
${t}$, intent(in) :: b(:)
65-
${t}$, intent(inout) :: x(:)
66-
${t}$, intent(in) :: tol
67-
integer, intent(in) :: maxiter
68-
type(solver_workspace_${s}$), intent(inout) :: workspace
74+
class(linop_${s}$), intent(in) :: A !! linear operator
75+
${t}$, intent(in) :: b(:) !! right-hand side vector
76+
${t}$, intent(inout) :: x(:) !! solution vector and initial guess
77+
${t}$, intent(in) :: tol !! tolerance for convergence
78+
integer, intent(in) :: maxiter !! maximum number of iterations
79+
type(solver_workspace_${s}$), intent(inout) :: workspace !! workspace for the solver
6980
end subroutine
7081
#:endfor
7182
end interface
@@ -75,34 +86,39 @@ module stdlib_linalg_iterative_solvers
7586
#:for matrix in MATRIX_TYPES
7687
#:for k, t, s in R_KINDS_TYPES
7788
module subroutine solve_cg_${matrix}$_${s}$(A,b,x,di,tol,maxiter,restart,workspace)
89+
!! linear operator matrix
7890
#:if matrix == "dense"
79-
${t}$, intent(in) :: A(:,:)
91+
${t}$, intent(in) :: A(:,:)
8092
#:else
8193
type(${matrix}$_${s}$_type), intent(in) :: A
8294
#:endif
83-
${t}$, intent(in) :: b(:)
84-
${t}$, intent(inout) :: x(:)
85-
${t}$, intent(in), optional :: tol
86-
logical(1), intent(in), optional, target :: di(:)
87-
integer, intent(in), optional :: maxiter
88-
logical, intent(in), optional :: restart
89-
type(solver_workspace_${s}$), optional, intent(inout), target :: workspace
95+
${t}$, intent(in) :: b(:) !! right-hand side vector
96+
${t}$, intent(inout) :: x(:) !! solution vector and initial guess
97+
${t}$, intent(in), optional :: tol !! tolerance for convergence
98+
logical(1), intent(in), optional, target :: di(:) !! dirichlet conditions mask
99+
integer, intent(in), optional :: maxiter !! maximum number of iterations
100+
logical, intent(in), optional :: restart !! restart flag
101+
type(solver_workspace_${s}$), optional, intent(inout), target :: workspace !! workspace for the solver
90102
end subroutine
91103
#:endfor
92104
#:endfor
93105
end interface
94106
public :: solve_cg
95107

108+
!! version: experimental
109+
!!
110+
!! solve_pcg_kernel interface for the preconditionned conjugate gradient method.
111+
!! [Specifications](../page/specs/stdlib_linalg_iterative_solvers.html#solve_pcg_kernel)
96112
interface solve_pcg_kernel
97113
#:for k, t, s in R_KINDS_TYPES
98114
module subroutine solve_pcg_kernel_${s}$(A,M,b,x,tol,maxiter,workspace)
99-
class(linop_${s}$), intent(in) :: A
100-
class(linop_${s}$), intent(in) :: M !> preconditionner
101-
${t}$, intent(in) :: b(:)
102-
${t}$, intent(inout) :: x(:)
103-
${t}$, intent(in) :: tol
104-
integer, intent(in) :: maxiter
105-
type(solver_workspace_${s}$), intent(inout) :: workspace
115+
class(linop_${s}$), intent(in) :: A !! linear operator
116+
class(linop_${s}$), intent(in) :: M !! preconditionner linear operator
117+
${t}$, intent(in) :: b(:) !! right-hand side vector
118+
${t}$, intent(inout) :: x(:) !! solution vector and initial guess
119+
${t}$, intent(in) :: tol !! tolerance for convergence
120+
integer, intent(in) :: maxiter !! maximum number of iterations
121+
type(solver_workspace_${s}$), intent(inout) :: workspace !! workspace for the solver
106122
end subroutine
107123
#:endfor
108124
end interface
@@ -112,20 +128,21 @@ module stdlib_linalg_iterative_solvers
112128
#:for matrix in MATRIX_TYPES
113129
#:for k, t, s in R_KINDS_TYPES
114130
module subroutine solve_pcg_${matrix}$_${s}$(A,b,x,di,tol,maxiter,restart,precond,M,workspace)
131+
!! linear operator matrix
115132
#:if matrix == "dense"
116133
${t}$, intent(in) :: A(:,:)
117134
#:else
118135
type(${matrix}$_${s}$_type), intent(in) :: A
119136
#:endif
120-
${t}$, intent(in) :: b(:)
121-
${t}$, intent(inout) :: x(:)
122-
${t}$, intent(in), optional :: tol
123-
logical(1), intent(in), optional, target :: di(:)
124-
integer, intent(in), optional :: maxiter
125-
logical, intent(in), optional :: restart
126-
integer, intent(in), optional :: precond !> preconditionner method flag
127-
class(linop_${s}$), optional , intent(in), target :: M !> preconditionner
128-
type(solver_workspace_${s}$), optional, intent(inout), target :: workspace
137+
${t}$, intent(in) :: b(:) !! right-hand side vector
138+
${t}$, intent(inout) :: x(:) !! solution vector and initial guess
139+
${t}$, intent(in), optional :: tol !! tolerance for convergence
140+
logical(1), intent(in), optional, target :: di(:) !! dirichlet conditions mask
141+
integer, intent(in), optional :: maxiter !! maximum number of iterations
142+
logical, intent(in), optional :: restart !! restart flag
143+
integer, intent(in), optional :: precond !! preconditionner method enumerator
144+
class(linop_${s}$), optional , intent(in), target :: M !! preconditionner linear operator
145+
type(solver_workspace_${s}$), optional, intent(inout), target :: workspace !! workspace for the solver
129146
end subroutine
130147
#:endfor
131148
#:endfor

0 commit comments

Comments
 (0)