Skip to content

Commit f027017

Browse files
jalveszjalvesz
committed
add factorizations
1 parent 05c076b commit f027017

File tree

3 files changed

+228
-8
lines changed

3 files changed

+228
-8
lines changed

‎src/stdlib_linalg_iterative_aux.fypp

Lines changed: 179 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -130,4 +130,183 @@ contains
130130
#:endfor
131131

132132

133+
!==============================================================
134+
! Factorization for preconditioners
135+
!==============================================================
136+
#:for k, t, s in R_KINDS_TYPES
137+
module subroutine factorize_ssor_dense_${s}$(A, w, L, D)
138+
${t}$, intent(in) :: A(:,:)
139+
${t}$, intent(in) :: w
140+
${t}$, intent(out) :: L(:,:)
141+
${t}$, intent(out) :: D(:)
142+
143+
integer(ilp) :: i, j, n
144+
${t}$ :: inv_w
145+
146+
n = size(A,dim=1)
147+
inv_w = 1._${s}$ / w
148+
do j = 1, n
149+
D(j) = A(j,j) * (2._${s}$-w)*inv_w
150+
do i = 1, j-1
151+
L(i,j) = 0._${s}$
152+
end do
153+
L(j,j) = L(j,j) * inv_w
154+
do i = j+1, n
155+
L(i,j) = A(i,j)
156+
end do
157+
end do
158+
159+
end subroutine
160+
161+
#:endfor
162+
163+
#:for k, t, s in R_KINDS_TYPES
164+
module subroutine factorize_ssor_csr_${s}$(A, w, L, D)
165+
type(CSR_${s}$_type), intent(in) :: A
166+
${t}$, intent(in) :: w
167+
type(CSR_${s}$_type), intent(out) :: L
168+
${t}$, intent(out) :: D(:)
169+
170+
integer(ilp) :: i, j, n
171+
${t}$ :: inv_w
172+
173+
inv_w = 1._${s}$ / w
174+
L%data = A%data
175+
select case(A%storage)
176+
case(sparse_lower)
177+
do i = 1, A%nrows
178+
D(i) = A%data( A%rowptr(i+1)-1 ) * (2._${s}$-w)*inv_w
179+
L%data( A%rowptr(i+1)-1 ) = L%data( A%rowptr(i+1)-1 ) * inv_w
180+
end do
181+
case(sparse_upper)
182+
do i = 1, A%nrows
183+
D(i) = A%data( A%rowptr(i) ) * (2._${s}$-w)*inv_w
184+
L%data( A%rowptr(i) ) = L%data( A%rowptr(i) ) * inv_w
185+
end do
186+
case(sparse_full)
187+
do i = 1, A%nrows
188+
do j = A%rowptr(i), A%rowptr(i+1)-1
189+
if( A%col(j) == i ) then
190+
D(i) = A%data(j) * (2._${s}$-w)*inv_w
191+
L%data(j) = L%data(j) * inv_w
192+
exit
193+
end if
194+
end do
195+
end do
196+
end select
197+
198+
end subroutine
199+
200+
#:endfor
201+
202+
!> Bunch-Kaufman factorization of a symmetric positive definite matrix A.
203+
!> The matrix A is assumed to be symmetric and positive definite.
204+
#:for k, t, s in R_KINDS_TYPES
205+
module subroutine factorize_ldlt_dense_${s}$(A, L, D)
206+
${t}$, intent(in) :: A(:,:)
207+
${t}$, intent(out) :: L(:,:)
208+
${t}$, intent(out) :: D(:)
209+
210+
${t}$:: auxsum
211+
integer(ilp) :: i, j, k, n
212+
213+
! Initialize L to zero
214+
n = size(A,dim=1)
215+
L = zero_${s}$
216+
217+
do i = 1, n
218+
! Compute D(i)
219+
auxsum = A(i, i)
220+
do k = 1, i - 1
221+
auxsum = auxsum - (L(i, k) ** 2) * D(k)
222+
end do
223+
D(i) = auxsum
224+
225+
! Check for positive definiteness
226+
if (D(i) <= zero_${s}$ ) D(i) = one_${s}$
227+
228+
! Set unit diagonal
229+
L(i, i) = one_${s}$
230+
231+
! Compute L(j, i) for j = i+1 to n
232+
do j = i + 1, n
233+
auxsum = A(j, i)
234+
do k = 1, i - 1
235+
auxsum = auxsum - L(j, k) * L(i, k) * D(k)
236+
end do
237+
L(j, i) = auxsum / D(i)
238+
end do
239+
end do
240+
end subroutine
241+
242+
#:endfor
243+
244+
#:for k, t, s in R_KINDS_TYPES
245+
module subroutine factorize_ldlt_csr_${s}$(A, L, D)
246+
type(CSR_${s}$_type), intent(in) :: A
247+
type(CSR_${s}$_type), intent(out) :: L
248+
${t}$, intent(out) :: D(:)
249+
250+
${t}$:: auxsum
251+
integer(ilp) :: i, j, k, n, p, q
252+
253+
integer, allocatable :: col_to_pos(:)
254+
${t}$, allocatable :: temp(:)
255+
256+
allocate(col_to_pos(n), source=-1)
257+
allocate(temp(n),source=zero_${s}$)
258+
259+
n = A%nrows
260+
do i = 1, n
261+
! Zero the temp workspace
262+
temp = zero_${s}$
263+
264+
! Load row i into temp (lower triangle only)
265+
do p = A%rowptr(i), A%rowptr(i+1) - 1
266+
j = A%col(p)
267+
if (j > i) cycle
268+
temp(j) = A%data(p)
269+
col_to_pos(j) = p ! map column to position in L%data
270+
end do
271+
272+
! Compute L(i, k) * D(k) * L(i, k) contributions
273+
do p = A%rowptr(i), A%rowptr(i+1) - 1
274+
k = A%col(p)
275+
if (k >= i) exit
276+
277+
auxsum = temp(k)
278+
279+
! Compute L(i, 1:k-1) ⋅ D(1:k-1) ⋅ L(k, 1:k-1)
280+
do q = A%rowptr(k), A%rowptr(k+1) - 1
281+
j = A%col(q)
282+
if (j >= k) exit
283+
auxsum = auxsum - L%data(col_to_pos(j)) * D(j) * L%data(q)
284+
end do
285+
286+
L%data(p) = auxsum / D(k)
287+
temp(k) = zero_${s}$ ! Clear temp slot
288+
end do
289+
290+
! Compute D(i)
291+
auxsum = temp(i)
292+
do p = A%rowptr(i), A%rowptr(i+1) - 1
293+
k = A%col(p)
294+
if (k >= i) exit
295+
auxsum = auxsum - L%data(p)**2 * D(k)
296+
end do
297+
298+
D(i) = merge( auxsum, one_${s}$ , auxsum > zero_${s}$)
299+
300+
! Set L(i, i) = 1.0
301+
L%data(col_to_pos(i)) = one_${s}$
302+
303+
! Reset col_to_pos map
304+
do p = A%rowptr(i), A%rowptr(i+1) - 1
305+
col_to_pos(A%col(p)) = -1
306+
end do
307+
end do
308+
end subroutine
309+
310+
#:endfor
311+
133312
end submodule stdlib_linalg_iterative_aux

‎src/stdlib_linalg_iterative_solvers.fypp

Lines changed: 34 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -167,6 +167,40 @@ module stdlib_linalg_iterative_solvers
167167
end interface
168168
public :: solve_backward_triangular
169169

170+
interface factorize_ldlt
171+
#:for k, t, s in R_KINDS_TYPES
172+
module subroutine factorize_ldlt_dense_${s}$(A, L, D)
173+
${t}$, intent(in) :: A(:,:)
174+
${t}$, intent(out) :: L(:,:)
175+
${t}$, intent(out) :: D(:)
176+
end subroutine
177+
module subroutine factorize_ldlt_csr_${s}$(A, L, D)
178+
type(CSR_${s}$_type), intent(in) :: A
179+
type(CSR_${s}$_type), intent(out) :: L
180+
${t}$, intent(out) :: D(:)
181+
end subroutine
182+
#:endfor
183+
end interface
184+
public :: factorize_ldlt
185+
186+
interface factorize_ssor
187+
#:for k, t, s in R_KINDS_TYPES
188+
module subroutine factorize_ssor_dense_${s}$(A, w, L, D)
189+
${t}$, intent(in) :: A(:,:)
190+
${t}$, intent(in) :: w
191+
${t}$, intent(out) :: L(:,:)
192+
${t}$, intent(out) :: D(:)
193+
end subroutine
194+
module subroutine factorize_ssor_csr_${s}$(A, w, L, D)
195+
type(CSR_${s}$_type), intent(in) :: A
196+
${t}$, intent(in) :: w
197+
type(CSR_${s}$_type), intent(out) :: L
198+
${t}$, intent(out) :: D(:)
199+
end subroutine
200+
#:endfor
201+
end interface
202+
public :: factorize_ssor
203+
170204

171205
contains
172206

‎src/stdlib_linalg_iterative_solvers_pccg.fypp

Lines changed: 15 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -7,7 +7,6 @@
77
submodule(stdlib_linalg_iterative_solvers) stdlib_linalg_iterative_pccg
88
use stdlib_kinds
99
use stdlib_sparse
10-
use stdlib_linalg, only: diag
1110
use stdlib_constants
1211
use stdlib_linalg_iterative_solvers
1312
implicit none
@@ -16,7 +15,7 @@ submodule(stdlib_linalg_iterative_solvers) stdlib_linalg_iterative_pccg
1615
enumerator :: none = 0
1716
enumerator :: jacobi
1817
enumerator :: ssor
19-
enumerator :: ichol
18+
enumerator :: ldlt
2019
end enum
2120

2221
contains
@@ -88,6 +87,7 @@ contains
8887
#:for k, t, s in R_KINDS_TYPES
8988
module subroutine solve_pccg_${matrix}$_${s}$(A,b,x,di,tol,maxiter,restart,precond,M,workspace)
9089
#:if matrix == "dense"
90+
use stdlib_linalg, only: diag
9191
${t}$, intent(in) :: A(:,:)
9292
#:else
9393
type(${matrix}$_${s}$_type), intent(in) :: A
@@ -126,43 +126,50 @@ contains
126126
precond_ = none ; if(present(precond)) precond_ = precond
127127
!-------------------------
128128
! internal memory setup
129-
op%apply => matvec
129+
! Preconditionner
130130
if(present(M)) then
131131
M_ => M
132132
else
133133
allocate( M_ )
134+
allocate(diagonal(n),source=zero_${s}$)
135+
134136
select case(precond_)
135137
case(jacobi)
136138
#:if matrix == "dense"
137139
diagonal = diag(A)
138140
#:else
139141
call diag(A,diagonal)
140142
#:endif
141-
where(abs(diagonal)>epsilon(0.d0)) diagonal = 1._dp/diagonal
142143
M_%apply => precond_jacobi
143144
case(ssor)
144145
#:if matrix == "dense"
145146
diagonal = diag(A)
146147
#:else
147148
call diag(A,diagonal)
148149
#:endif
149-
where(abs(diagonal)>epsilon(0.d0)) diagonal = 1._dp/diagonal
150-
! TODO: Extract the lower triangular part of A in L
150+
L = A !> copy A structure to L
151+
call factorize_ssor( A , one_${s}$ , L, diagonal )
151152
M_%apply => precond_ldl
152-
case(ichol)
153-
! TODO: Add LDL^t factorization
153+
case(ldlt)
154+
L = A !> copy A structure to L
155+
call factorize_ldlt( A , L, diagonal )
154156
M_%apply => precond_ldl
155157
case default
156158
M_%apply => precond_none
157159
end select
160+
where(abs(diagonal)>epsilon(zero_${s}$)) diagonal = one_${s}$/diagonal
158161
end if
162+
! matvec for the operator
163+
op%apply => matvec
159164

165+
! direchlet boundary conditions mask
160166
if(present(di))then
161167
di_ => di
162168
else
163169
allocate(di_(n),source=.false._1)
164170
end if
165171

172+
! workspace for the solver
166173
if(present(workspace)) then
167174
workspace_ => workspace
168175
else

0 commit comments

Comments
 (0)