preconditioners.f90 25 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
! preconditioners.f90 --- preconditioners for symmetric 5-point stencil system
!
! Copyright  (C)  2010  Florian Wilhelm <Florian.Wilhelm@kit.edu>
!
! Version: 1.0
! Keywords: scales ppm solver preconditioners ilu0 icc jacobi ssor
! Author: Florian Wilhelm <Florian.Wilhelm@kit.edu>
! Maintainer: Florian Wilhelm <Florian.Wilhelm@kit.edu>
! URL: https://www.dkrz.de/redmine/projects/show/scales-ppm
!
! Redistribution and use in source and binary forms, with or without
! modification, are  permitted provided that the following conditions are
! met:
!
! Redistributions of source code must retain the above copyright notice,
! this list of conditions and the following disclaimer.
!
! Redistributions in binary form must reproduce the above copyright
! notice, this list of conditions and the following disclaimer in the
! documentation and/or other materials provided with the distribution.
!
! Neither the name of the DKRZ GmbH nor the names of its contributors
! may be used to endorse or promote products derived from this software
! without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
! IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
! PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
! OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
! EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
! PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
! PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
! LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
! NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
!
! Commentary:
!
! Code:
!

43
MODULE preconditioners
44
    USE solver_public
45
    USE solver_internal
46
    USE mo_kind, ONLY: wp, dp
47
    USE ppm_base, ONLY: abort_ppm
48
    USE ppm_extents, ONLY: extent_start, extent_end, extent_size, extent_type => extent
49
50
51
52
53

    IMPLICIT NONE

    PRIVATE

Florian Wilhelm's avatar
Florian Wilhelm committed
54
    !-------------------------------------------------------------------
55
    ! Internal State Variables
56
57
58
    !-------------------------------------------------------------------
    ! Auxiliary variable for jacobi preconditioner
    REAL(wp), ALLOCATABLE :: jacobi_diag(:,:)
59
    ! Auxiliary variables for ilu0 preconditioner
60
    REAL(wp), ALLOCATABLE :: ilu0_diag(:,:)
61
62
    ! Auxiliary variables for icc(p) preconditioner
    REAL(wp), ALLOCATABLE :: ICC_C(:,:), ICC_W(:,:,:), ICC_S(:,:,:)
63
64
65
66

    !-------------------------------------------------------------------
    ! Module Function & Procedure Interface
    !-------------------------------------------------------------------
67
    PUBLIC :: NONE_PRECOND, JACOBI_PRECOND, ILU0_PRECOND, SSOR_PRECOND, ICC_PRECOND
68
    PUBLIC :: identity, prep_jacobi, jacobi, prep_ilu0, ilu0, ssor, prep_icc, icc, precond_prepared &
69
70
71
              , jacobi_precond_stencil, ilu0_precond_stencil, ssor_precond_stencil, icc_precond_stencil &
              , jacobi_precond_shifted_stencil, ilu0_precond_shifted_stencil, ssor_precond_shifted_stencil &
              , icc_precond_shifted_stencil
Florian Wilhelm's avatar
Florian Wilhelm committed
72

73
    CONTAINS
Florian Wilhelm's avatar
Florian Wilhelm committed
74

75
    ! Identity function, use this as preconditioner to get the non-preconditioned CG-method
Florian Wilhelm's avatar
Florian Wilhelm committed
76
    SUBROUTINE identity(r)
77
78
        REAL(wp), INTENT(INOUT) :: r(:,:)
    END SUBROUTINE identity
Florian Wilhelm's avatar
Florian Wilhelm committed
79

80
81
82
83
    ! Preparation for Jacobi preconditioner
    SUBROUTINE prep_jacobi()
        INTEGER :: ib, jb, il, jl, i, j
        REAL(wp), DIMENSION(:,:), POINTER :: uf, vf, ff
Florian Wilhelm's avatar
Florian Wilhelm committed
84

85
        ! Define some variables for the ranges of the fields
Florian Wilhelm's avatar
Florian Wilhelm committed
86
87
88
89
        ib = extent_start(stencil%extent(1))
        jb = extent_start(stencil%extent(2))
        il = extent_end(stencil%extent(1))
        jl = extent_end(stencil%extent(2))
90
91

        ALLOCATE(jacobi_diag(ib:il,jb:jl))
Florian Wilhelm's avatar
Florian Wilhelm committed
92

93
94
95
96
97
        DO j=jb,jl
            DO i=ib,il
                jacobi_diag(i,j) = 1._wp/stencil%central(i,j)
            ENDDO
        ENDDO
Florian Wilhelm's avatar
Florian Wilhelm committed
98

99
    END SUBROUTINE prep_jacobi
Florian Wilhelm's avatar
Florian Wilhelm committed
100

101
102
103
104
    ! Jacobi preconditioner
    SUBROUTINE jacobi(r)
        REAL(wp), INTENT(INOUT) :: r(:,:)
        INTEGER :: ib, jb, il, jl, i, j
Florian Wilhelm's avatar
Florian Wilhelm committed
105

106
        ! Define some variables for the ranges of the fields
Florian Wilhelm's avatar
Florian Wilhelm committed
107
108
109
110
111
        ib = extent_start(stencil%extent(1))
        jb = extent_start(stencil%extent(2))
        il = extent_end(stencil%extent(1))
        jl = extent_end(stencil%extent(2))

112
113
114
115
116
        DO j=jb,jl
            DO i=ib,il
                r(i,j) = r(i,j)*jacobi_diag(i,j)
            ENDDO
        ENDDO
Florian Wilhelm's avatar
Florian Wilhelm committed
117

118
    END SUBROUTINE jacobi
119
120
121
122
123

     ! Preparation of the ILU(0) preconditioner
    SUBROUTINE prep_ilu0()
        INTEGER :: ib, jb, il, jl, i, j
        REAL(wp), DIMENSION(:,:), POINTER :: uf, vf, ff
Florian Wilhelm's avatar
Florian Wilhelm committed
124

125
        ! Define some variables for the ranges of the fields
Florian Wilhelm's avatar
Florian Wilhelm committed
126
127
128
129
        ib = extent_start(stencil%extent(1))
        jb = extent_start(stencil%extent(2))
        il = extent_end(stencil%extent(1))
        jl = extent_end(stencil%extent(2))
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154

        ALLOCATE(ilu0_diag(ib:il,jb:jl))

        ! Define some aliases
        uf => stencil%zonal
        vf => stencil%meridional
        ff => stencil%central

        ! Calculate diagonal elements of ILU(0) factorization
        ilu0_diag = 0.
        ilu0_diag(ib,jb) = ff(ib,jb)
        DO i=ib+1,il
            ilu0_diag(i,jb) = ff(i,jb) - uf(i-1,jb)**2/ilu0_diag(i-1,jb)
        ENDDO
        DO j=jb+1,jl
            ilu0_diag(ib,j) = ff(ib,j) - vf(ib,j-1)**2/ilu0_diag(ib,j-1)
        ENDDO
        DO j=jb+1,jl
           DO i=ib+1,il
                ilu0_diag(i,j) = ff(i,j) - uf(i-1,j)**2/ilu0_diag(i-1,j) &
                                         - vf(i,j-1)**2/ilu0_diag(i,j-1)
           ENDDO
        ENDDO

    END SUBROUTINE prep_ilu0
Florian Wilhelm's avatar
Florian Wilhelm committed
155

156
157
158
159
160
    ! ILU(0) preconditioner
    SUBROUTINE ilu0(r)
        REAL(wp), INTENT(INOUT) :: r(:,:)
        INTEGER :: ib, jb, il, jl, i, j
        REAL(wp), DIMENSION(:,:), POINTER :: uf, vf, ff
Florian Wilhelm's avatar
Florian Wilhelm committed
161

162
        ! Define some variables for the ranges of the fields
Florian Wilhelm's avatar
Florian Wilhelm committed
163
164
165
166
167
        ib = extent_start(stencil%extent(1))
        jb = extent_start(stencil%extent(2))
        il = extent_end(stencil%extent(1))
        jl = extent_end(stencil%extent(2))

168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
        ! Define some aliases
        uf => stencil%zonal
        vf => stencil%meridional
        ff => stencil%central

        ! The overall system to solve is As = r <-> (LU)s = r

        ! 1.) Solve system Ly = r
        DO i=ib+1,il
            r(i,jb) = r(i,jb) + uf(i-1,jb)/ilu0_diag(i-1,jb)*r(i-1,jb)
        ENDDO
        DO j=jb+1,jl
            r(ib,j) = r(ib,j) + vf(ib,j-1)/ilu0_diag(ib,j-1)*r(ib,j-1)
        ENDDO
        DO j = jb+1,jl
            DO i = ib+1,il
                r(i,j) = r(i,j) + uf(i-1,j)/ilu0_diag(i-1,j)*r(i-1,j) &
                                + vf(i,j-1)/ilu0_diag(i,j-1)*r(i,j-1)
            ENDDO
        ENDDO

        ! 2.) Solve system Us = y
        r(il,jl) = 1._wp/ilu0_diag(il,jl)*r(il,jl)
        DO i=il-1,ib,-1
            r(i,jl) = 1._wp/ilu0_diag(i,jl)*( r(i,jl) + uf(i,jl)*r(i+1,jl) )
        ENDDO
        DO j=jl-1,jb,-1
            r(il,j) = 1._wp/ilu0_diag(il,j)*( r(il,j) + vf(il,j)*r(il,j+1) )
        ENDDO
        DO j = jl-1,jb,-1
            DO i = il-1,ib,-1
                r(i,j) = 1._wp/ilu0_diag(i,j)*( r(i,j) + uf(i,j)*r(i+1,j)   &
                                                       + vf(i,j)*r(i,j+1) )
            ENDDO
        ENDDO

    END SUBROUTINE ilu0
Florian Wilhelm's avatar
Florian Wilhelm committed
205

206
    ! Symmetric SOR preconditioner
207
    SUBROUTINE ssor(r)
208
        USE solver_config, ONLY: get_ssor_param
209
210
211
212
213
214
215
216
217
218

        REAL(wp), INTENT(INOUT) :: r(:,:)
        INTEGER :: i, j, ib, jb, il, jl
        REAL(wp):: ssorpar
        REAL(wp), DIMENSION(:,:), POINTER :: uf, vf, ff

        ! Get SSOR parameter from configuration
        ssorpar = get_ssor_param()

        ! Define some variables for the ranges of the fields
Florian Wilhelm's avatar
Florian Wilhelm committed
219
220
221
222
223
        ib = extent_start(stencil%extent(1))
        jb = extent_start(stencil%extent(2))
        il = extent_end(stencil%extent(1))
        jl = extent_end(stencil%extent(2))

224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
        ! Define some aliases
        uf => stencil%zonal
        vf => stencil%meridional
        ff => stencil%central

        !r(ib,jb) = r(ib,jb)
        DO i=ib+1,il
           r(i,jb) = r(i,jb) + ssorpar*uf(i-1,jb)/ff(i-1,jb)*r(i-1,jb)
        ENDDO
        DO j=jb+1,jl
           r(ib,j) = r(ib,j) + ssorpar*vf(ib,j-1)/ff(ib,j-1)*r(ib,j-1)
        ENDDO
        DO j=jb+1,jl
          DO i=ib+1,il
             r(i,j) = r(i,j) + ssorpar*( uf(i-1,j)/ff(i-1,j)*r(i-1,j)  &
                                       + vf(i,j-1)/ff(i,j-1)*r(i,j-1) )
          ENDDO
        ENDDO

        r(il,jl) = r(il,jl)/ff(il,jl)
        DO i=il-1,ib,-1
           r(i,jl) = ( r(i,jl) + ssorpar*uf(i,jl)*r(i+1,jl) )/ff(i,jl)
        ENDDO
        DO j=jl-1,jb,-1
           r(il,j) = ( r(il,j) + ssorpar*vf(il,j)*r(il,j+1) )/ff(il,j)
        ENDDO
        DO j = jl-1,jb,-1
           DO i = il-1,ib,-1
              r(i,j) = ( r(i,j) + ssorpar*( uf(i,j)*r(i+1,j)           &
                                          + vf(i,j)*r(i,j+1) ) )/ff(i,j)
           ENDDO
        ENDDO

    END SUBROUTINE ssor
Florian Wilhelm's avatar
Florian Wilhelm committed
258

259
260
261
262
263
264
265
    ! Calculates Matrix for Incomplete Cholesky Decomposition (ICC) with fillin p
    ! This routine makes use of a transformation from the barotropic stencil to a
    ! real matrix A with indices m_i, m_j. Variables s_i, s_j refer to the stencil
    ! coordinates, i, j are used according the actual context.
    ! The result L of L*L' ~= A is saved in ICC_C, ICC_W, ICC_S
    SUBROUTINE prep_icc(p)
        INTEGER, INTENT(IN) :: p  ! FILL-IN
266
        INTEGER :: ib, il, jb, jl, m_i, m_j, m_n, m_m, k
267
268
269
270
        REAL(wp) :: tmp
        REAL(wp), DIMENSION(:,:), POINTER :: uf, vf, ff

        ! Define some variables for the ranges of the fields
Florian Wilhelm's avatar
Florian Wilhelm committed
271
272
273
274
275
        ib = extent_start(stencil%extent(1))
        jb = extent_start(stencil%extent(2))
        il = extent_end(stencil%extent(1))
        jl = extent_end(stencil%extent(2))

276
277
278
279
280
281
282
        ! Define some aliases
        uf => stencil%zonal
        vf => stencil%meridional
        ff => stencil%central

        ! MAX(1,p-1) handles structurel exception for p=1
        ! ICC_C = Center, _W = West, _S = South in stencil
283
        ALLOCATE(ICC_C(ib:il,jb:jl), ICC_W(MAX(1,p-1),ib:il,jb:jl), ICC_S(1+p,ib:il,jb:jl))
284
285

        ! Initialize
286
        ICC_C = 1.
287
288
289
        ICC_W = 0.
        ICC_S = 0.

290
291
        m_n = extent_size(stencil%extent)    ! dim of barotropic Matrix, m_ = Matrix
        m_m = extent_size(stencil%extent(1)) ! m_m is distance between diagonal and outer diagonal
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314

        DO m_i = 1,m_n

            ! start treating diagonal
            ! L(i,i) = sqrt( A(i,i) - sum(L(i,1:i-1).^2) ) in OCTAVE
            tmp = get_value_of_A(m_i,m_i)

            ! outer diagonal
            DO k = MAX(1,m_i-m_m),m_i-m_m+p
                tmp = tmp - get_value_of_L(m_i,k)**2
            ENDDO

            ! secondary diagonals
            DO k = MAX(1,m_i-MAX(1,p-1)),m_i-1
                tmp = tmp - get_value_of_L(m_i,k)**2
            ENDDO

            tmp = SQRT(tmp)
            CALL set_value_of_L(m_i, m_i, tmp)
            ! end treating diagonal

            ! start treating secondary diagonals
            DO m_j = m_i+1,m_i+MAX(1,p-1)
315
                IF ( m_j > m_n ) EXIT
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
                ! L(j,i) = ( A(j,i) - sum(L(j,1:i-1).*L(i,1:i-1)) )/L(i,i) in OCTAVE
                tmp = get_value_of_A(m_j, m_i)

                ! outer diagonal
                DO k = MAX(1,m_i-m_m),m_i-m_m+p
                    tmp = tmp - get_value_of_L(m_j,k) * get_value_of_L(m_i,k)
                ENDDO

                ! secondary diagonals
                DO k = MAX(1,m_i-MAX(1,p-1)),m_i-1
                    tmp = tmp - get_value_of_L(m_j,k) * get_value_of_L(m_i,k)
                ENDDO

                tmp = tmp / get_value_of_L(m_i,m_i)
                CALL set_value_of_L(m_j, m_i, tmp)
            ENDDO
            ! end treating secondary diagonals

            ! start treating outer diagonals
            DO m_j = m_i+m_m-p,m_i+m_m
336
                IF ( m_j > m_n ) EXIT
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
                ! L(j,i) = ( A(j,i) - sum(L(j,1:i-1).*L(i,1:i-1)) )/L(i,i) in OCTAVE
                tmp = get_value_of_A(m_j, m_i)

                ! outer diagonal
                DO k = MAX(1,m_i-m_m),m_i-m_m+p
                    tmp = tmp - get_value_of_L(m_j,k) * get_value_of_L(m_i,k)
                ENDDO

                ! secondary diagonals
                DO k = MAX(1,m_i-MAX(1,p-1)),m_i-1
                    tmp = tmp - get_value_of_L(m_j,k) * get_value_of_L(m_i,k)
                ENDDO

                tmp = tmp / get_value_of_L(m_i,m_i)
                CALL set_value_of_L(m_j, m_i, tmp)
            ENDDO
            ! end treating outer diagonals

        ENDDO
Florian Wilhelm's avatar
Florian Wilhelm committed
356

357
358
        ! Calculate inverse of ICC_C to avoid division afterwards
        ICC_C = 1._wp/ICC_C
359
360
361
362
363
364

        CONTAINS

        ! Calculates stencil coordinates given a row in matrix A
        SUBROUTINE coord_by_row(row, i, j)
            INTEGER, INTENT(IN) :: row
Florian Wilhelm's avatar
Florian Wilhelm committed
365
            INTEGER, INTENT(OUT) :: i, j
366

367
            i = MOD(row - 1, m_m) + ib
368
            j = (row - 1) / m_m + jb
369
370
371
372
373
374

        END SUBROUTINE coord_by_row

        ! Get value in matrix L
        FUNCTION get_value_of_L(i,j) RESULT(res)
            INTEGER, INTENT(IN) :: i, j
375
            INTEGER :: d, s_i, s_j
376
377
378
379
380
381
382
383
384
385
386
387
            REAL(wp) :: res

            res = 0.

            ! diagonal
            IF ( i == j) THEN
                CALL coord_by_row(i, s_i, s_j)
                res = ICC_C(s_i, s_j)
                RETURN
            ENDIF

            ! secondary diagonal
388
            d = i - j
389
390
391
392
393
394
395
            IF ( d >= 1 .AND. d <= MAX(1,p-1) ) THEN
                CALL coord_by_row(i, s_i, s_j)
                res = ICC_W(d, s_i, s_j)
                RETURN
            ENDIF

            ! outer diagonal
396
397
            IF ( i - m_m > 0 ) THEN
                d = j - (i - m_m) + 1
398
399
400
401
402
403
404
405
406
407
408
409
410
                IF ( d >= 1 .AND. d <= p+1 ) THEN
                    CALL coord_by_row(i, s_i, s_j)
                    res = ICC_S(d, s_i, s_j)
                    RETURN
                ENDIF
            ENDIF

        END FUNCTION get_value_of_L

        ! Set value in matrix L
        SUBROUTINE set_value_of_L(i,j,v)
            INTEGER, INTENT(IN) :: i, j
            REAL(wp), INTENT(IN) :: v
411
            INTEGER :: d, s_i, s_j
412
413
414
415
416
417
418
419
420

            ! diagonal
            IF ( i == j) THEN
                CALL coord_by_row(i, s_i, s_j)
                ICC_C(s_i, s_j) = v
                RETURN
            ENDIF

            ! secondary diagonal
421
            d = i - j
422
423
424
425
426
427
428
            IF ( d >= 1 .AND. d <= MAX(1,p-1) ) THEN
                CALL coord_by_row(i, s_i, s_j)
                ICC_W(d, s_i, s_j) = v
                RETURN
            ENDIF

            ! outer diagonal
429
430
            IF ( i - m_m > 0 ) THEN
                d = j - (i - m_m) + 1
431
432
433
434
435
436
437
438
439
440
441
442
                IF ( d >= 1 .AND. d <= p+1 ) THEN
                    CALL coord_by_row(i, s_i, s_j)
                    ICC_S(d, s_i, s_j) = v
                    RETURN
                ENDIF
            ENDIF

        END SUBROUTINE set_value_of_L

        ! Get value of virtual barotropic matrix A
        FUNCTION get_value_of_A(i,j) RESULT(res)
            INTEGER, INTENT(IN) :: i, j
443
            INTEGER :: s_i, s_j
444
445
446
447
448
449
450
451
452
            REAL(wp) :: res

            res = 0.

            IF ( i == j) THEN ! diagonal
                CALL coord_by_row(i, s_i, s_j)
                res = ff(s_i, s_j)
            ELSEIF ( j-i == 1 ) THEN ! right secondary diagonal
                CALL coord_by_row(i, s_i, s_j)
453
                IF (s_i >= il) RETURN ! hit boundary
454
455
456
                res = -uf(s_i,s_j)
            ELSEIF ( j-i == -1 ) THEN ! left secondary diagonal
                CALL coord_by_row(i, s_i, s_j)
457
                IF (s_i <= ib) RETURN ! hit boundary
458
                res = -uf(s_i-1,s_j)
459
            ELSEIF ( j-i == m_m ) THEN ! right outer diagonal
460
                CALL coord_by_row(i, s_i, s_j)
461
                IF (s_j >= jl) RETURN ! hit boundary
462
                res = -vf(s_i,s_j)
463
            ELSEIF ( j-i == -m_m ) THEN ! left outer diagonal
464
465
466
467
468
469
470
471
                CALL coord_by_row(i, s_i, s_j)
                IF (s_j <= jb) RETURN ! hit boundary
                res = -vf(s_i,s_j-1)
            ENDIF

        END FUNCTION get_value_of_A

    END SUBROUTINE prep_icc
Florian Wilhelm's avatar
Florian Wilhelm committed
472

473
474
475
    ! ICC(p) preconditioner
    SUBROUTINE icc(r)
        REAL(wp), INTENT(INOUT) :: r(:,:)
476
        INTEGER :: ib, il, jb, jl, i, j, k, p, s_l, w_l
477
478
479
        REAL(wp), DIMENSION(:,:), POINTER :: uf, vf, ff

        ! Define some variables for the ranges of the fields
Florian Wilhelm's avatar
Florian Wilhelm committed
480
481
482
483
484
        ib = extent_start(stencil%extent(1))
        jb = extent_start(stencil%extent(2))
        il = extent_end(stencil%extent(1))
        jl = extent_end(stencil%extent(2))

485
486
487
488
489
        ! Define some aliases
        uf => stencil%zonal
        vf => stencil%meridional
        ff => stencil%central

490
491
492
        ! Set halos to zero
        CALL clear_halos(r, stencil%extent)

493
494
495
496
497
498
499
500
        ! get the fill-in p of ICC
        p = SIZE(ICC_S,1) - 1
        s_l = p + 1 ! length of south diagonal
        w_l = MAX(1,p-1) ! length of west diagonal

        ! The overall system to solve is As = r <-> (LL')s = r

        ! 1.) Solve system Ly = r
Florian Wilhelm's avatar
Florian Wilhelm committed
501

502
        ! first column
503
        r(ib,jb) = r(ib,jb)*ICC_C(ib,jb)
504
        DO i=ib+1,il
505
           r(i,jb) = ( r(i,jb) - ICC_W(1,i,jb)*r(i-1,jb) )*ICC_C(i,jb)
506
507
508
        ENDDO

        ! all other columns
509
510
        DO j = jb+1,jl
            DO i = ib,il
511
512
513
514
515
516
517
                DO k = 1,s_l
                    r(i,j) = r(i,j) - ICC_S(k,i,j)*r(i+k-1,j-1)
                ENDDO
                DO k = 1,w_l
                    r(i,j) = r(i,j) - ICC_W(k,i,j)*r(i-k,j)
                ENDDO
                r(i,j) = r(i,j)*ICC_C(i,j)
518
519
520
521
            ENDDO
        ENDDO

        ! 2.) Solve system L's = y
522
        ! Backward substitution is done on a column level instead of row-wise
523

524
525
526
527
        ! all but first columns
        DO j = jl,jb+1,-1
            DO i = il,ib,-1
                r(i,j) = r(i,j)*ICC_C(i,j)
Florian Wilhelm's avatar
Florian Wilhelm committed
528
529
                DO k = w_l,1,-1
                    r(i-k,j) = r(i-k,j) - r(i,j)*ICC_W(k,i,j) ! saxpy
530
531
532
533
                ENDDO
                DO k = s_l,1,-1
                    r(i+k-1,j-1) = r(i+k-1,j-1) - r(i,j)*ICC_S(k,i,j) ! saxpy
                ENDDO
534
            ENDDO
535
        ENDDO
Florian Wilhelm's avatar
Florian Wilhelm committed
536

537
        ! first column
538
        DO i = il,ib+1,-1
Florian Wilhelm's avatar
Florian Wilhelm committed
539
540
            r(i,jb) = r(i,jb)*ICC_C(i,jb)
            r(i-1,jb) = r(i-1,jb) - r(i,jb)*ICC_W(1,i,jb) ! saxpy
541
        ENDDO
542
        r(ib,jb) = r(ib,jb)*ICC_C(ib,jb)
543
544

    END SUBROUTINE icc
545

546
    ! Check if a certain preconditioner is already prepared
547
548
549
    FUNCTION precond_prepared(preconditioner) RESULT(ans)
        INTEGER, INTENT(IN) :: preconditioner
        LOGICAL :: ans
Florian Wilhelm's avatar
Florian Wilhelm committed
550

551
552
553
        SELECT CASE (preconditioner)
        CASE (NONE_PRECOND)
            ans = .TRUE.
554
555
        CASE (JACOBI_PRECOND)
            ans = ALLOCATED(jacobi_diag)
556
557
558
559
560
561
562
        CASE (ILU0_PRECOND)
            ans = ALLOCATED(ilu0_diag)
        CASE (SSOR_PRECOND)
            ans = config%ssor_param /= 0._wp
        CASE (ICC_PRECOND)
            ans = ALLOCATED(ICC_C) .AND. ALLOCATED(ICC_W) .AND. ALLOCATED(ICC_S)
        CASE DEFAULT
Thomas Jahns's avatar
Thomas Jahns committed
563
564
565
566
            CALL abort_ppm("Selected preconditioner [" &
                 // int2str(preconditioner) // "] does not exist!", &
                 __FILE__, &
                 __LINE__)
567
        END SELECT
Florian Wilhelm's avatar
Florian Wilhelm committed
568

569
    END FUNCTION precond_prepared
Florian Wilhelm's avatar
Florian Wilhelm committed
570

571
    ! Jacobi preconditioned matrix stencil operation
572
573
574
575
576
577
578
579
580
581
582
    SUBROUTINE jacobi_precond_stencil(field, res_field)
        USE solver_config, ONLY: apply_stencil

        REAL(wp), INTENT(IN) :: field(:,:)
        REAL(wp), INTENT(OUT) :: res_field(:,:)

        ! Apply stencil first
        IF (.NOT. precond_prepared(JACOBI_PRECOND)) THEN
            CALL prep_jacobi()
        ENDIF
        CALL apply_stencil(field, res_field)
Florian Wilhelm's avatar
Florian Wilhelm committed
583

584
585
        ! Apply Jacobi preconditioner
        CALL jacobi(res_field)
Florian Wilhelm's avatar
Florian Wilhelm committed
586

587
    END SUBROUTINE jacobi_precond_stencil
588
589
590
591
592
593
594
595
596
597

    ! Jacobi preconditioned matrix stencil operation shifted by a value
    SUBROUTINE jacobi_precond_shifted_stencil(field, res_field)
        USE solver_config, ONLY: apply_stencil

        REAL(wp), INTENT(IN) :: field(:,:)
        REAL(wp), INTENT(OUT) :: res_field(:,:)
        INTEGER :: i, j , ib, jb, il, jl

        ! Define some variables for the ranges of the fields
Florian Wilhelm's avatar
Florian Wilhelm committed
598
599
600
        ib = extent_start(stencil%extent(1))
        jb = extent_start(stencil%extent(2))
        il = extent_end(stencil%extent(1))
601
602
603
604
605
606
607
        jl = extent_end(stencil%extent(2))

        ! Apply stencil first
        IF (.NOT. precond_prepared(JACOBI_PRECOND)) THEN
            CALL prep_jacobi()
        ENDIF
        CALL apply_stencil(field, res_field)
Florian Wilhelm's avatar
Florian Wilhelm committed
608

609
610
611
612
613
        ! Apply Jacobi preconditioner
        CALL jacobi(res_field)

        ! Perform the shift
        FORALL(i = ib:il, j = jb:jl) res_field(i,j) = stencil%shift*field(i,j) - res_field(i,j)
Florian Wilhelm's avatar
Florian Wilhelm committed
614

615
    END SUBROUTINE jacobi_precond_shifted_stencil
Florian Wilhelm's avatar
Florian Wilhelm committed
616

617
    ! ILU(0) preconditioned matrix stencil operation
618
619
    SUBROUTINE ilu0_precond_stencil(field, res_field)
        USE solver_config, ONLY: apply_stencil
Florian Wilhelm's avatar
Florian Wilhelm committed
620

621
622
623
624
625
        REAL(wp), INTENT(IN) :: field(:,:)
        REAL(wp), INTENT(OUT) :: res_field(:,:)

        ! Apply stencil first
        CALL apply_stencil(field, res_field)
Florian Wilhelm's avatar
Florian Wilhelm committed
626

627
628
629
630
631
        ! Apply ILU0 preconditioner
        IF (.NOT. precond_prepared(ILU0_PRECOND)) THEN
            CALL prep_ilu0()
        ENDIF
        CALL ilu0(res_field)
Florian Wilhelm's avatar
Florian Wilhelm committed
632

633
    END SUBROUTINE ilu0_precond_stencil
634
635
636
637

    ! ILU(0) preconditioned matrix stencil operation shifted by a value
    SUBROUTINE ilu0_precond_shifted_stencil(field, res_field)
        USE solver_config, ONLY: apply_stencil
Florian Wilhelm's avatar
Florian Wilhelm committed
638

639
640
641
642
643
        REAL(wp), INTENT(IN) :: field(:,:)
        REAL(wp), INTENT(OUT) :: res_field(:,:)
        INTEGER :: i, j , ib, jb, il, jl

        ! Define some variables for the ranges of the fields
Florian Wilhelm's avatar
Florian Wilhelm committed
644
645
646
        ib = extent_start(stencil%extent(1))
        jb = extent_start(stencil%extent(2))
        il = extent_end(stencil%extent(1))
647
648
649
650
        jl = extent_end(stencil%extent(2))

        ! Apply stencil first
        CALL apply_stencil(field, res_field)
Florian Wilhelm's avatar
Florian Wilhelm committed
651

652
653
654
655
656
657
658
659
        ! Apply ILU0 preconditioner
        IF (.NOT. precond_prepared(ILU0_PRECOND)) THEN
            CALL prep_ilu0()
        ENDIF
        CALL ilu0(res_field)

        ! Perform the shift
        FORALL(i = ib:il, j = jb:jl) res_field(i,j) = stencil%shift*field(i,j) - res_field(i,j)
Florian Wilhelm's avatar
Florian Wilhelm committed
660

661
    END SUBROUTINE ilu0_precond_shifted_stencil
Florian Wilhelm's avatar
Florian Wilhelm committed
662

663
    ! SSOR preconditioned matrix stencil operation
664
665
    SUBROUTINE ssor_precond_stencil(field, res_field)
        USE solver_config, ONLY: apply_stencil
Florian Wilhelm's avatar
Florian Wilhelm committed
666

667
668
669
670
671
        REAL(wp), INTENT(IN) :: field(:,:)
        REAL(wp), INTENT(OUT) :: res_field(:,:)

        ! Apply stencil first
        CALL apply_stencil(field, res_field)
Florian Wilhelm's avatar
Florian Wilhelm committed
672

673
674
        ! Apply SSOR preconditioner
        IF (.NOT. precond_prepared(SSOR_PRECOND)) THEN
Thomas Jahns's avatar
Thomas Jahns committed
675
676
677
            CALL abort_ppm("No SSOR parameter provided!", &
                 __FILE__, &
                 __LINE__)
678
679
        ENDIF
        CALL ssor(res_field)
Florian Wilhelm's avatar
Florian Wilhelm committed
680

681
682
    END SUBROUTINE ssor_precond_stencil

683
684
685
686
687
688
689
690
691
    ! SSOR preconditioned matrix stencil operation shifted by a value
    SUBROUTINE ssor_precond_shifted_stencil(field, res_field)
        USE solver_config, ONLY: apply_stencil

        REAL(wp), INTENT(IN) :: field(:,:)
        REAL(wp), INTENT(OUT) :: res_field(:,:)
        INTEGER :: i, j , ib, jb, il, jl

        ! Define some variables for the ranges of the fields
Florian Wilhelm's avatar
Florian Wilhelm committed
692
693
694
        ib = extent_start(stencil%extent(1))
        jb = extent_start(stencil%extent(2))
        il = extent_end(stencil%extent(1))
695
        jl = extent_end(stencil%extent(2))
Florian Wilhelm's avatar
Florian Wilhelm committed
696

697
698
        ! Apply stencil first
        CALL apply_stencil(field, res_field)
Florian Wilhelm's avatar
Florian Wilhelm committed
699

700
701
        ! Apply SSOR preconditioner
        IF (.NOT. precond_prepared(SSOR_PRECOND)) THEN
Thomas Jahns's avatar
Thomas Jahns committed
702
703
704
            CALL abort_ppm("No SSOR parameter provided!", &
                 __FILE__, &
                 __LINE__)
705
706
707
708
709
        ENDIF
        CALL ssor(res_field)

        ! Perform the shift
        FORALL(i = ib:il, j = jb:jl) res_field(i,j) = stencil%shift*field(i,j) - res_field(i,j)
Florian Wilhelm's avatar
Florian Wilhelm committed
710

711
712
    END SUBROUTINE ssor_precond_shifted_stencil

713
    ! ICC(p) preconditioned matrix stencil operation
714
715
716
717
718
719
720
721
    SUBROUTINE icc_precond_stencil(field, res_field)
        USE solver_config, ONLY: apply_stencil

        REAL(wp), INTENT(IN) :: field(:,:)
        REAL(wp), INTENT(OUT) :: res_field(:,:)

        ! Apply stencil first
        CALL apply_stencil(field, res_field)
Florian Wilhelm's avatar
Florian Wilhelm committed
722

723
724
725
726
727
        ! Apply ICC preconditioner
        IF (.NOT. precond_prepared(ICC_PRECOND)) THEN
            CALL prep_icc(config%icc_param)
        ENDIF
        CALL icc(res_field)
Florian Wilhelm's avatar
Florian Wilhelm committed
728

729
730
    END SUBROUTINE icc_precond_stencil

731
732
733
    ! ICC(p) preconditioned matrix stencil operation shifted by a value
    SUBROUTINE icc_precond_shifted_stencil(field, res_field)
        USE solver_config, ONLY: apply_stencil
Florian Wilhelm's avatar
Florian Wilhelm committed
734

735
736
737
738
739
        REAL(wp), INTENT(IN) :: field(:,:)
        REAL(wp), INTENT(OUT) :: res_field(:,:)
        INTEGER :: i, j , ib, jb, il, jl

        ! Define some variables for the ranges of the fields
Florian Wilhelm's avatar
Florian Wilhelm committed
740
741
742
743
        ib = extent_start(stencil%extent(1))
        jb = extent_start(stencil%extent(2))
        il = extent_end(stencil%extent(1))
        jl = extent_end(stencil%extent(2))
744
745
746

        ! Apply stencil first
        CALL apply_stencil(field, res_field)
Florian Wilhelm's avatar
Florian Wilhelm committed
747

748
749
750
751
752
753
754
755
        ! Apply ICC preconditioner
        IF (.NOT. precond_prepared(ICC_PRECOND)) THEN
            CALL prep_icc(config%icc_param)
        ENDIF
        CALL icc(res_field)

        ! Perform the shift
        FORALL(i = ib:il, j = jb:jl) res_field(i,j) = stencil%shift*field(i,j) - res_field(i,j)
Florian Wilhelm's avatar
Florian Wilhelm committed
756

757
758
    END SUBROUTINE icc_precond_shifted_stencil

759
END MODULE preconditioners