Skip to content

Commit 5e15a33

Browse files
jalveszjalvesz
committed
add forward/backward solvers for preconditionning
1 parent 379cd81 commit 5e15a33

4 files changed

+180
-3
lines changed

‎src/stdlib_linalg_iterative_aux.fypp

Lines changed: 137 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,137 @@
1+
#:include "common.fypp"
2+
#:set R_KINDS_TYPES = list(zip(REAL_KINDS, REAL_TYPES, REAL_SUFFIX))
3+
#:set MATRIX_TYPES = ["dense", "CSR"]
4+
#:set RANKS = range(1, 2+1)
5+
6+
submodule(stdlib_linalg_iterative_solvers) stdlib_linalg_iterative_aux
7+
use stdlib_kinds
8+
use stdlib_sparse
9+
use stdlib_constants
10+
use stdlib_linalg_iterative_solvers
11+
use stdlib_constants
12+
implicit none
13+
14+
integer, parameter :: ilp = int32
15+
16+
contains
17+
18+
#:for k, t, s in R_KINDS_TYPES
19+
module subroutine solve_forward_triangular_dense_${s}$(L,b,x)
20+
${t}$, intent(in) :: L(:,:)
21+
${t}$, intent(in) :: b(:)
22+
${t}$, intent(inout) :: x(:)
23+
${t}$ :: aux
24+
25+
integer(ilp) :: i, j, m
26+
27+
do i = 1, size(L,dim=1)
28+
aux = zero_${s}$
29+
do j = 1, i-1
30+
aux = aux + L(i,j)*x(j)
31+
end do
32+
x(i) = b(i) - aux
33+
end do
34+
end subroutine
35+
36+
module subroutine solve_backward_triangular_dense_${s}$(U,b,x)
37+
${t}$, intent(in) :: U(:,:)
38+
${t}$, intent(in) :: b(:)
39+
${t}$, intent(inout) :: x(:)
40+
${t}$ :: aux
41+
${t}$ :: baux(size(x))
42+
43+
integer(ilp) :: i, j, m
44+
45+
baux = zero_${s}$
46+
do i = size(U,dim=1), 1, -1
47+
x(i) = x(i) - baux(i)
48+
do j = i+1, size(U,dim=2)
49+
baux(j) = baux(j) + U(i,j)*x(i)
50+
end do
51+
end do
52+
end subroutine
53+
54+
#:endfor
55+
56+
#:for k, t, s in R_KINDS_TYPES
57+
module subroutine solve_forward_triangular_csr_${s}$(L,b,x)
58+
type(CSR_${s}$_type), intent(in) :: L
59+
${t}$, intent(in) :: b(:)
60+
${t}$, intent(inout) :: x(:)
61+
${t}$ :: aux
62+
63+
integer(ilp) :: i, j, m
64+
select case (L%storage)
65+
case(sparse_full)
66+
do i = 1, L%nrows
67+
aux = zero_${s}$
68+
do m = L%rowptr(i), L%rowptr(i+1)-1
69+
j = L%col(m)
70+
if(j>i) cycle !> skip upper part of the matrix
71+
aux = aux + L%data(m)*x(j)
72+
end do
73+
x(i) = b(i) - aux
74+
end do
75+
case(sparse_lower)
76+
do i = 1, L%nrows
77+
aux = zero_${s}$
78+
do m = L%rowptr(i), L%rowptr(i+1)-2
79+
j = L%col(m)
80+
aux = aux + L%data(m)*x(j)
81+
end do
82+
x(i) = b(i) - aux
83+
end do
84+
case(sparse_upper) !> treates as lower triangular (thus transpose)
85+
do i = 1, L%nrows
86+
x(i) = b(i)
87+
do m = L%rowptr(i)+1, L%rowptr(i+1)-1
88+
j = L%col(m)
89+
x(j) = x(j) - L%data(m)*x(i)
90+
end do
91+
end do
92+
end select
93+
end subroutine
94+
95+
module subroutine solve_backward_triangular_csr_${s}$(U,b,x)
96+
type(CSR_${s}$_type), intent(in) :: U
97+
${t}$, intent(in) :: b(:)
98+
${t}$, intent(inout) :: x(:)
99+
${t}$ :: aux
100+
${t}$ :: baux(size(x))
101+
102+
integer(ilp) :: i, j, m
103+
104+
baux = zero_${s}$
105+
select case (U%storage)
106+
case(sparse_full)
107+
do i = U%nrows, 1, -1
108+
x(i) = b(i) - baux(i)
109+
do m = U%rowptr(i), U%rowptr(i+1)-1
110+
j = U%col(m)
111+
if(j<i) cycle !> skip lower part of the matrix
112+
baux(j) = baux(j) + U%data(m)*x(i)
113+
end do
114+
end do
115+
case(sparse_lower)
116+
do i = U%nrows, 1, -1
117+
x(i) = b(i) - baux(i)
118+
do m = U%rowptr(i), U%rowptr(i+1)-2
119+
j = U%col(m)
120+
baux(j) = baux(j) + U%data(m)*x(i)
121+
end do
122+
end do
123+
case(sparse_upper)
124+
do i = U%nrows, 1, -1
125+
x(i) = b(i)
126+
do m = U%rowptr(i)+1, U%rowptr(i+1)-1
127+
j = U%col(m)
128+
x(i) = x(i) - U%data(m)*x(j)
129+
end do
130+
end do
131+
end select
132+
end subroutine
133+
134+
#:endfor
135+
136+
137+
end submodule stdlib_linalg_iterative_aux

‎src/stdlib_linalg_iterative_solvers.fypp

Lines changed: 36 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
#:include "common.fypp"
22
#:set R_KINDS_TYPES = list(zip(REAL_KINDS, REAL_TYPES, REAL_SUFFIX))
33
#:set C_KINDS_TYPES = list(zip(CMPLX_KINDS, CMPLX_TYPES, CMPLX_SUFFIX))
4-
#:set MATRIX_TYPES = ["dense", "CSR", "CSC", "ELL", "SELLC"]
4+
#:set MATRIX_TYPES = ["dense", "CSR"]
55
#:set RANKS = range(1, 2+1)
66

77
module stdlib_linalg_iterative_solvers
@@ -131,6 +131,41 @@ module stdlib_linalg_iterative_solvers
131131
end interface
132132
public :: solve_pccg
133133

134+
interface solve_forward_triangular
135+
!> Solve the system L x = b, where L is a strictly lower triangular matrix (diagonal assumed = 1).
136+
#:for k, t, s in R_KINDS_TYPES
137+
module subroutine solve_forward_triangular_dense_${s}$(L,b,x)
138+
${t}$, intent(in) :: L(:,:)
139+
${t}$, intent(in) :: b(:)
140+
${t}$, intent(inout) :: x(:)
141+
end subroutine
142+
module subroutine solve_forward_triangular_csr_${s}$(L,b,x)
143+
type(CSR_${s}$_type), intent(in) :: L
144+
${t}$, intent(in) :: b(:)
145+
${t}$, intent(inout) :: x(:)
146+
end subroutine
147+
#:endfor
148+
end interface
149+
public :: solve_forward_triangular
150+
151+
152+
interface solve_backward_triangular
153+
!> Solve the system U x = b, where U is a strictly upper triangular matrix (diagonal assumed = 1).
154+
#:for k, t, s in R_KINDS_TYPES
155+
module subroutine solve_backward_triangular_dense_${s}$(U,b,x)
156+
${t}$, intent(in) :: U(:,:)
157+
${t}$, intent(in) :: b(:)
158+
${t}$, intent(inout) :: x(:)
159+
end subroutine
160+
module subroutine solve_backward_triangular_csr_${s}$(U,b,x)
161+
type(CSR_${s}$_type), intent(in) :: U
162+
${t}$, intent(in) :: b(:)
163+
${t}$, intent(inout) :: x(:)
164+
end subroutine
165+
#:endfor
166+
end interface
167+
public :: solve_backward_triangular
168+
134169

135170
contains
136171

‎src/stdlib_linalg_iterative_solvers_cg.fypp

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
#:include "common.fypp"
22
#:set R_KINDS_TYPES = list(zip(REAL_KINDS, REAL_TYPES, REAL_SUFFIX))
33
#:set C_KINDS_TYPES = list(zip(CMPLX_KINDS, CMPLX_TYPES, CMPLX_SUFFIX))
4-
#:set MATRIX_TYPES = ["dense", "CSR", "CSC", "ELL", "SELLC"]
4+
#:set MATRIX_TYPES = ["dense", "CSR"]
55
#:set RANKS = range(1, 2+1)
66

77
submodule(stdlib_linalg_iterative_solvers) stdlib_linalg_iterative_cg

‎src/stdlib_linalg_iterative_solvers_pccg.fypp

Lines changed: 6 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
#:include "common.fypp"
22
#:set R_KINDS_TYPES = list(zip(REAL_KINDS, REAL_TYPES, REAL_SUFFIX))
33
#:set C_KINDS_TYPES = list(zip(CMPLX_KINDS, CMPLX_TYPES, CMPLX_SUFFIX))
4-
#:set MATRIX_TYPES = ["dense", "CSR", "CSC", "ELL", "SELLC"]
4+
#:set MATRIX_TYPES = ["dense", "CSR"]
55
#:set RANKS = range(1, 2+1)
66

77
submodule(stdlib_linalg_iterative_solvers) stdlib_linalg_iterative_pccg
@@ -115,6 +115,11 @@ contains
115115
allocate( M_ )
116116
M_%apply => default_preconditionner
117117
end if
118+
!> TODO: factorize the preconditionner
119+
!> Jacobi: M = diag(A)^{-1}
120+
!> SOR: M = (w/(2-w)) (D/w+L) D^{-1} (D/w+U) ; A = D + L + U
121+
!> ILU: M = ( LU )^{-1} ; A = LU
122+
!> ICHOL: M = (L*D*U)^{-1} ; A = L*D*U
118123
if(present(di))then
119124
di_ => di
120125
else

0 commit comments

Comments
 (0)