Skip to content

Commit acefaaf

Browse files
jalveszjalvesz
committed
review csr factorization
1 parent c596ac0 commit acefaaf

File tree

2 files changed

+47
-71
lines changed

2 files changed

+47
-71
lines changed

‎src/stdlib_linalg_iterative_aux.fypp

Lines changed: 39 additions & 63 deletions
Original file line numberDiff line numberDiff line change
@@ -137,8 +137,8 @@ contains
137137
module subroutine factorize_ssor_dense_${s}$(A, w, L, D)
138138
${t}$, intent(in) :: A(:,:)
139139
${t}$, intent(in) :: w
140-
${t}$, intent(out) :: L(:,:)
141-
${t}$, intent(out) :: D(:)
140+
${t}$, intent(inout) :: L(:,:)
141+
${t}$, intent(inout) :: D(:)
142142

143143
integer(ilp) :: i, j, n
144144
${t}$ :: inv_w
@@ -164,8 +164,8 @@ contains
164164
module subroutine factorize_ssor_csr_${s}$(A, w, L, D)
165165
type(CSR_${s}$_type), intent(in) :: A
166166
${t}$, intent(in) :: w
167-
type(CSR_${s}$_type), intent(out) :: L
168-
${t}$, intent(out) :: D(:)
167+
type(CSR_${s}$_type), intent(inout) :: L
168+
${t}$, intent(inout) :: D(:)
169169

170170
integer(ilp) :: i, j, n
171171
${t}$ :: inv_w
@@ -204,8 +204,8 @@ contains
204204
#:for k, t, s in R_KINDS_TYPES
205205
module subroutine factorize_ldlt_dense_${s}$(A, L, D)
206206
${t}$, intent(in) :: A(:,:)
207-
${t}$, intent(out) :: L(:,:)
208-
${t}$, intent(out) :: D(:)
207+
${t}$, intent(inout) :: L(:,:)
208+
${t}$, intent(inout) :: D(:)
209209

210210
${t}$:: auxsum
211211
integer(ilp) :: i, j, k, n
@@ -244,67 +244,43 @@ contains
244244
#:for k, t, s in R_KINDS_TYPES
245245
module subroutine factorize_ldlt_csr_${s}$(A, L, D)
246246
type(CSR_${s}$_type), intent(in) :: A
247-
type(CSR_${s}$_type), intent(out) :: L
248-
${t}$, intent(out) :: D(:)
247+
type(CSR_${s}$_type), intent(inout) :: L
248+
${t}$, intent(inout) :: D(:)
249249

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}$)
250+
${t}$:: aux_diag, nondiag
251+
integer(ilp) :: i, j, k, m, n, p, q, r, ad1, ad2
258252

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
253+
D = zero_${s}$
254+
L%data = zero_${s}$
271255

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
256+
select case(A%storage)
257+
case(sparse_lower)
258+
do i = 1, A%nrows
259+
ad1 = A%rowptr(i)
260+
ad2 = A%rowptr(i+1)-1
261+
L%data(ad2) = one_${s}$
262+
aux_diag = zero_${s}$
263+
do m = ad1, ad2-1
264+
j = A%col(m)
265+
nondiag = A%data(m)
266+
q = ad1
267+
do p = A%rowptr(j), A%rowptr(j+1)-2
268+
k = A%col(p)
269+
do r = q, m-1
270+
n = A%col(r)
271+
if( k == n ) then
272+
nondiag = nondiag - L%data(r)*L%data(p)*D(k) ! Aij - Lik Ljk Dk
273+
q = r + 1
274+
end if
275+
end do
276+
end do
277+
L%data(m) = nondiag / D(j) ! Lij = (Aij - Lik Ljk Dk)/Dj
278+
aux_diag = aux_diag + (L%data(m) * nondiag)
279+
end do
280+
D(i) = A%data(ad2) - aux_diag ! Dj = Ajj - Ljk^2 Dk
281+
if (D(i) <= zero_${s}$) D(i)= one_${s}$
306282
end do
307-
end do
283+
end select
308284
end subroutine
309285

310286
#:endfor

‎src/stdlib_linalg_iterative_solvers.fypp

Lines changed: 8 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -171,13 +171,13 @@ module stdlib_linalg_iterative_solvers
171171
#:for k, t, s in R_KINDS_TYPES
172172
module subroutine factorize_ldlt_dense_${s}$(A, L, D)
173173
${t}$, intent(in) :: A(:,:)
174-
${t}$, intent(out) :: L(:,:)
175-
${t}$, intent(out) :: D(:)
174+
${t}$, intent(inout) :: L(:,:)
175+
${t}$, intent(inout) :: D(:)
176176
end subroutine
177177
module subroutine factorize_ldlt_csr_${s}$(A, L, D)
178178
type(CSR_${s}$_type), intent(in) :: A
179-
type(CSR_${s}$_type), intent(out) :: L
180-
${t}$, intent(out) :: D(:)
179+
type(CSR_${s}$_type), intent(inout) :: L
180+
${t}$, intent(inout) :: D(:)
181181
end subroutine
182182
#:endfor
183183
end interface
@@ -188,14 +188,14 @@ module stdlib_linalg_iterative_solvers
188188
module subroutine factorize_ssor_dense_${s}$(A, w, L, D)
189189
${t}$, intent(in) :: A(:,:)
190190
${t}$, intent(in) :: w
191-
${t}$, intent(out) :: L(:,:)
192-
${t}$, intent(out) :: D(:)
191+
${t}$, intent(inout) :: L(:,:)
192+
${t}$, intent(inout) :: D(:)
193193
end subroutine
194194
module subroutine factorize_ssor_csr_${s}$(A, w, L, D)
195195
type(CSR_${s}$_type), intent(in) :: A
196196
${t}$, intent(in) :: w
197-
type(CSR_${s}$_type), intent(out) :: L
198-
${t}$, intent(out) :: D(:)
197+
type(CSR_${s}$_type), intent(inout) :: L
198+
${t}$, intent(inout) :: D(:)
199199
end subroutine
200200
#:endfor
201201
end interface

0 commit comments

Comments
 (0)