Skip to content

Commit 716b3c5

Browse files
jalveszjalvesz
committed
start iterative solvers
1 parent 69eaa20 commit 716b3c5

File tree

4 files changed

+185
-10
lines changed

4 files changed

+185
-10
lines changed

‎src/CMakeLists.txt

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -92,6 +92,8 @@ set(cppFiles
9292

9393
stdlib_linalg_blas.fypp
9494
stdlib_linalg_lapack.fypp
95+
stdlib_linalg_iterative_solvers.fypp
96+
stdlib_linalg_iterative_solvers_cg.fypp
9597
)
9698

9799
add_subdirectory(blas)

‎src/stdlib_linalg_iterative_solvers.fypp

Lines changed: 77 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,77 @@
1+
#:include "common.fypp"
2+
#:set R_KINDS_TYPES = list(zip(REAL_KINDS, REAL_TYPES, REAL_SUFFIX))
3+
#:set C_KINDS_TYPES = list(zip(CMPLX_KINDS, CMPLX_TYPES, CMPLX_SUFFIX))
4+
#:set MATRIX_TYPES = ["dense", "CSR"]
5+
#:set RANKS = range(1, 2+1)
6+
7+
module stdlib_linalg_iterative_solvers
8+
use stdlib_kinds
9+
use stdlib_sparse
10+
implicit none
11+
private
12+
13+
#:for k, t, s in R_KINDS_TYPES
14+
type, public :: linop_${s}$
15+
procedure(vector_sub_${s}$), nopass, pointer :: matvec => null()
16+
procedure(vector_sub_${s}$), nopass, pointer :: preconditionner => null()
17+
procedure(reduction_sub_${s}$), nopass, pointer :: inner_product => null()
18+
end type
19+
#:endfor
20+
21+
#:for k, t, s in R_KINDS_TYPES
22+
type, abstract, public :: solver_workspace_${s}$
23+
end type
24+
type, public, extends(solver_workspace_${s}$) :: cg_workspace_${s}$
25+
${t}$, allocatable :: r(:), p(:), Ap(:)
26+
end type
27+
28+
#:endfor
29+
30+
abstract interface
31+
#:for k, t, s in R_KINDS_TYPES
32+
subroutine vector_sub_${s}$(x,y)
33+
import :: ${k}$
34+
${t}$, intent(in) :: x(:)
35+
${t}$, intent(inout) :: y(:)
36+
end subroutine
37+
pure ${t}$ function reduction_sub_${s}$(x,y) result(r)
38+
import :: ${k}$
39+
${t}$, intent(in) :: x(:)
40+
${t}$, intent(in) :: y(:)
41+
end function
42+
#:endfor
43+
end interface
44+
45+
interface solve_cg_generic
46+
#:for k, t, s in R_KINDS_TYPES
47+
module subroutine solve_cg_generic_${s}$(A,b,x,tol,maxiter,workspace)
48+
class(linop_${s}$), intent(in) :: A
49+
${t}$, intent(in) :: b(:), tol
50+
${t}$, intent(inout) :: x(:)
51+
integer, intent(in) :: maxiter
52+
type(cg_workspace_${s}$), intent(inout) :: workspace
53+
end subroutine
54+
#:endfor
55+
end interface
56+
public :: solve_cg_generic
57+
58+
interface solve_cg
59+
#:for matrix in MATRIX_TYPES
60+
#:for k, t, s in R_KINDS_TYPES
61+
module subroutine solve_cg_${matrix}$_${s}$(A,b,x,tol,maxiter,workspace)
62+
#:if matrix == "dense"
63+
${t}$, intent(in) :: A(:,:)
64+
#:else
65+
type(${matrix}$_${s}$_type), intent(in) :: A
66+
#:endif
67+
${t}$, intent(in) :: b(:), tol
68+
${t}$, intent(inout) :: x(:)
69+
integer, intent(in), optional :: maxiter
70+
type(cg_workspace_${s}$), optional, intent(inout), target :: workspace
71+
end subroutine
72+
#:endfor
73+
#:endfor
74+
end interface
75+
public :: solve_cg
76+
77+
end module stdlib_linalg_iterative_solvers

‎src/stdlib_linalg_iterative_solvers_cg.fypp

Lines changed: 105 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,105 @@
1+
#:include "common.fypp"
2+
#:set R_KINDS_TYPES = list(zip(REAL_KINDS, REAL_TYPES, REAL_SUFFIX))
3+
#:set C_KINDS_TYPES = list(zip(CMPLX_KINDS, CMPLX_TYPES, CMPLX_SUFFIX))
4+
#:set MATRIX_TYPES = ["dense", "CSR"]
5+
#:set RANKS = range(1, 2+1)
6+
7+
submodule(stdlib_linalg_iterative_solvers) stdlib_linalg_iterative_cg
8+
use stdlib_kinds
9+
use stdlib_sparse
10+
use stdlib_constants
11+
use stdlib_intrinsics, only: dot_product => stdlib_dot_product
12+
use stdlib_linalg_iterative_solvers
13+
implicit none
14+
15+
contains
16+
17+
#:for k, t, s in R_KINDS_TYPES
18+
module subroutine solve_cg_generic_${s}$(A,b,x,tol,maxiter,workspace)
19+
class(linop_${s}$), intent(in) :: A
20+
${t}$, intent(in) :: b(:), tol
21+
${t}$, intent(inout) :: x(:)
22+
integer, intent(in) :: maxiter
23+
type(cg_workspace_${s}$), intent(inout) :: workspace
24+
!-------------------------
25+
integer :: iter
26+
${t}$ :: rtr, rtrold, alpha, beta, norm0_sq
27+
!-------------------------
28+
associate( p => workspace%p, &
29+
r => workspace%r, &
30+
Ap => workspace%Ap)
31+
x = zero_${s}$
32+
rtr = A%inner_product(r, r)
33+
norm0_sq = A%inner_product(b, b)
34+
p = b
35+
beta = zero_${s}$
36+
iter = 1
37+
do while( rtr > tol**2 * norm0_sq .and. iter < maxiter)
38+
p = r + beta * p
39+
call A%matvec(p,Ap)
40+
alpha = rtr / A%inner_product(p, Ap)
41+
x = x + alpha * p
42+
r = r - alpha * Ap
43+
rtrold = rtr
44+
rtr = A%inner_product(r, r)
45+
beta = rtr / rtrold
46+
iter = iter + 1
47+
end do
48+
end associate
49+
end subroutine
50+
#:endfor
51+
52+
#:for matrix in MATRIX_TYPES
53+
#:for k, t, s in R_KINDS_TYPES
54+
module subroutine solve_cg_${matrix}$_${s}$(A,b,x,tol,maxiter,workspace)
55+
#:if matrix == "dense"
56+
${t}$, intent(in) :: A(:,:)
57+
#:else
58+
type(${matrix}$_${s}$_type), intent(in) :: A
59+
#:endif
60+
${t}$, intent(in) :: b(:), tol
61+
${t}$, intent(inout) :: x(:)
62+
integer, intent(in), optional :: maxiter
63+
type(cg_workspace_${s}$), optional, intent(inout), target :: workspace
64+
!-------------------------
65+
type(linop_${s}$) :: op
66+
type(cg_workspace_${s}$), pointer :: workspace_
67+
integer :: n, maxiter_
68+
!-------------------------
69+
n = size(b)
70+
op%matvec => default_matvec
71+
op%inner_product => default_dot
72+
73+
maxiter_ = n
74+
if(present(maxiter)) maxiter_ = maxiter
75+
if(present(workspace)) then
76+
workspace_ => workspace
77+
else
78+
allocate( workspace_%r(n), &
79+
workspace_%p(n), &
80+
workspace_%Ap(n))
81+
82+
end if
83+
call solve_cg_generic(op,b,x,tol,maxiter_,workspace_)
84+
contains
85+
86+
subroutine default_matvec(x,y)
87+
${t}$, intent(in) :: x(:)
88+
${t}$, intent(inout) :: y(:)
89+
#:if matrix == "dense"
90+
y = matmul(A,x)
91+
#:else
92+
call spmv( A , x, y )
93+
#:endif
94+
end subroutine
95+
pure ${t}$ function default_dot(x,y) result(r)
96+
${t}$, intent(in) :: x(:)
97+
${t}$, intent(in) :: y(:)
98+
r = dot_product(x,y)
99+
end function
100+
end subroutine
101+
102+
#:endfor
103+
#:endfor
104+
105+
end submodule stdlib_linalg_iterative_cg

‎src/stdlib_sparse_constants.fypp

Lines changed: 1 addition & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -3,7 +3,7 @@
33
#:set C_KINDS_TYPES = list(zip(CMPLX_KINDS, CMPLX_TYPES, CMPLX_SUFFIX))
44
module stdlib_sparse_constants
55
use stdlib_kinds, only: int8, int16, int32, int64, sp, dp, xdp, qp
6-
6+
use stdlib_constants
77
implicit none
88
public
99

@@ -20,13 +20,4 @@ module stdlib_sparse_constants
2020
! Integer size support for ILP64 builds should be done here
2121
integer, parameter :: ilp = int32
2222

23-
#:for k1, t1, s1 in (R_KINDS_TYPES)
24-
${t1}$, parameter :: zero_${s1}$ = 0._${k1}$
25-
${t1}$, parameter :: one_${s1}$ = 1._${k1}$
26-
#:endfor
27-
#:for k1, t1, s1 in (C_KINDS_TYPES)
28-
${t1}$, parameter :: zero_${s1}$ = (0._${k1}$,0._${k1}$)
29-
${t1}$, parameter :: one_${s1}$ = (1._${k1}$,1._${k1}$)
30-
#:endfor
31-
3223
end module stdlib_sparse_constants

0 commit comments

Comments
 (0)