solvers.f90 18.4 KB
Newer Older
Florian Wilhelm's avatar
Florian Wilhelm committed
1
! solvers.f90 --- provides solvers ie. preconditioned CG, Chebyshev, Schwarz
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
!
! Copyright  (C)  2010  Florian Wilhelm <Florian.Wilhelm@kit.edu>
!
! Version: 1.0
! Keywords: scales ppm solver cg chebyshev schwarz
! 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 solvers
44
    USE solver_public
45
    USE solver_internal
46
    USE linear_algebra
Florian Wilhelm's avatar
Florian Wilhelm committed
47
    USE mo_kind, ONLY: wp, dp
48
    USE ppm_extents, ONLY: extent_start, extent_end, extent_size, extent_type => extent
49
    USE ppm_base, ONLY: abort_ppm
Florian Wilhelm's avatar
Florian Wilhelm committed
50

Florian Wilhelm's avatar
Florian Wilhelm committed
51
52
53
    IMPLICIT NONE

    PRIVATE
Florian Wilhelm's avatar
Florian Wilhelm committed
54

Florian Wilhelm's avatar
Florian Wilhelm committed
55
56
57
    !-------------------------------------------------------------------
    ! Module Function & Procedure Interface
    !-------------------------------------------------------------------
58
    PUBLIC :: CG_SOLVER, CHEBYSHEV_SOLVER
59
    PUBLIC :: cg_method, precond_cg_method, chebyshev_method, precond_chebyshev_method, schwarz_method
Florian Wilhelm's avatar
Florian Wilhelm committed
60

Florian Wilhelm's avatar
Florian Wilhelm committed
61
    CONTAINS
Florian Wilhelm's avatar
Florian Wilhelm committed
62

Florian Wilhelm's avatar
Florian Wilhelm committed
63
64
65
66
67
68
69
    ! Preconditioned CG-Method with A*x=b and startvalue x
    FUNCTION precond_cg_method(A, b, x, ext_x, precond, exchange, tol_opt, maxiter_opt) RESULT(kiter)
        REAL(wp), INTENT(IN) :: b(:,:)                  ! right-hand-side b
        REAL(wp), INTENT(INOUT) :: x(:,:)               ! startvalue and result
        TYPE(extent_type), INTENT(IN) :: ext_x(:)       ! extent of x
        REAL(wp), OPTIONAL, INTENT(IN) :: tol_opt       ! tolerance for residual
        INTEGER, OPTIONAL, INTENT(IN) :: maxiter_opt    ! maximum iterations
70
71
72
        INTEGER :: i, j, jb, ib, il, jl, ie, je, inner_size, kiter, maxiter
        REAL(wp) :: delta, delta0, delta_old, dAd, alpha, beta, tol

Florian Wilhelm's avatar
Florian Wilhelm committed
73
74
        ! Auxiliary fields for CG-Method
        REAL(wp), ALLOCATABLE :: d(:,:), r(:,:), Ad(:,:), s(:,:)
75

Florian Wilhelm's avatar
Florian Wilhelm committed
76
77
        INTERFACE
            ! Matrix-Vector multiplication of the linear system to solve
Florian Wilhelm's avatar
Florian Wilhelm committed
78
            SUBROUTINE A(field, res_field)
Florian Wilhelm's avatar
Florian Wilhelm committed
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
                USE mo_kind, ONLY: wp, dp
                REAL(wp), INTENT(IN) :: field(:,:)
                REAL(wp), INTENT(OUT) :: res_field(:,:)
            END SUBROUTINE
            ! Preconditioner
            SUBROUTINE precond(r)
                USE mo_kind, ONLY: wp, dp
                REAL(wp), INTENT(INOUT) :: r(:,:)
            END SUBROUTINE
            ! Function to exchange boundaries if necessary
            SUBROUTINE exchange(a0, text)
                USE mo_kind, ONLY: wp, dp
                REAL(wp), INTENT(INOUT) :: a0(:,:)
                CHARACTER (LEN=*), INTENT(IN), OPTIONAL :: text
            END SUBROUTINE exchange
        END INTERFACE

        ! Check and set optional arguments
        IF ( PRESENT(maxiter_opt) ) THEN
            maxiter = maxiter_opt
        ELSE
            maxiter = config%maxiter
        ENDIF
        IF ( PRESENT(tol_opt) ) THEN
            tol = tol_opt**2 ! **2, because we compare later the square of a norm
        ELSE
            tol = config%tol**2
        ENDIF
Florian Wilhelm's avatar
Florian Wilhelm committed
107

Florian Wilhelm's avatar
Florian Wilhelm committed
108
109
110
        ! Define some variables for the ranges of the fields
        ie = SIZE(x,1)
        je = SIZE(x,2)
Florian Wilhelm's avatar
Florian Wilhelm committed
111
112
113
114
        ib = extent_start(ext_x(1))
        jb = extent_start(ext_x(2))
        il = extent_end(ext_x(1))
        jl = extent_end(ext_x(2))
115
        inner_size = extent_size(ext_x)
Florian Wilhelm's avatar
Florian Wilhelm committed
116

Florian Wilhelm's avatar
Florian Wilhelm committed
117
        ALLOCATE(d(ie,je), r(ie,je), Ad(ie,je), s(ie,je))
118

119
120
        ! Initialize boundary/halos of d
        CALL clear_halos(d, ext_x)
Florian Wilhelm's avatar
Florian Wilhelm committed
121
122
123
124

        ! Initialise residual and direction
        CALL A(x, Ad)
        r(ib:il,jb:jl) = b(ib:il,jb:jl) - Ad(ib:il,jb:jl)
Florian Wilhelm's avatar
Florian Wilhelm committed
125

Florian Wilhelm's avatar
Florian Wilhelm committed
126
        ! Apply Preconditioner
Florian Wilhelm's avatar
Florian Wilhelm committed
127
        d(ib:il,jb:jl) = r(ib:il,jb:jl)
Florian Wilhelm's avatar
Florian Wilhelm committed
128
129
130
        CALL precond(d)

        ! Calculate scalar product of residual and direction
131
        delta = arr_dotproduct(r(ib:il,jb:jl), d(ib:il,jb:jl))
Florian Wilhelm's avatar
Florian Wilhelm committed
132
133
134

        delta0 = delta
        kiter = 0
Florian Wilhelm's avatar
Florian Wilhelm committed
135

Florian Wilhelm's avatar
Florian Wilhelm committed
136
137
138
139
140
141
        DO WHILE (kiter < maxiter .AND. delta > tol*delta0)
            IF (config%do_exchange) CALL exchange(d,'precond_cg_method 1')
            kiter = kiter + 1

            ! Initialise auxiliary variables
            CALL A(d, Ad)
Florian Wilhelm's avatar
Florian Wilhelm committed
142

143
            dAd = arr_dotproduct( d(ib:il,jb:jl), Ad(ib:il,jb:jl) )
Florian Wilhelm's avatar
Florian Wilhelm committed
144

145
            alpha = delta/dAd
Florian Wilhelm's avatar
Florian Wilhelm committed
146

Florian Wilhelm's avatar
Florian Wilhelm committed
147
            ! One step of preconditioned CG-method
Florian Wilhelm's avatar
Florian Wilhelm committed
148
#ifndef BLAS
Florian Wilhelm's avatar
Florian Wilhelm committed
149
150
151
            x(ib:il,jb:jl) = x(ib:il,jb:jl) + alpha*d(ib:il,jb:jl)
            r(ib:il,jb:jl) = r(ib:il,jb:jl) - alpha*Ad(ib:il,jb:jl)
#else
152
            CALL DAXPY(inner_size,alpha,d(ib:il,jb:jl),1,x(ib:il,jb:jl),1)
Florian Wilhelm's avatar
Florian Wilhelm committed
153
154
            CALL DAXPY(inner_size,-alpha,Ad(ib:il,jb:jl),1,r(ib:il,jb:jl),1)
#endif
Florian Wilhelm's avatar
Florian Wilhelm committed
155
156
157
158
159
160
161

            ! Apply Preconditioner
            s(ib:il,jb:jl) = r(ib:il,jb:jl)
            CALL precond(s)

            ! Calculate scalar product of residual and direction
            delta_old = delta
162
163
            delta = arr_dotproduct( r(ib:il,jb:jl), s(ib:il,jb:jl) )

Florian Wilhelm's avatar
Florian Wilhelm committed
164
            beta = delta/delta_old
Florian Wilhelm's avatar
Florian Wilhelm committed
165
166
167

            !write(*,*) 'PCG iteration: ', kiter, ' current delta: ', delta

Florian Wilhelm's avatar
Florian Wilhelm committed
168
            ! Compute new direction d
Florian Wilhelm's avatar
Florian Wilhelm committed
169
#ifndef BLAS
Florian Wilhelm's avatar
Florian Wilhelm committed
170
171
            d(ib:il,jb:jl) = s(ib:il,jb:jl) + beta*d(ib:il,jb:jl)
#else
172
173
            CALL DSCAL(inner_size,beta,d(ib:il,jb:jl),1)
            CALL DAXPY(inner_size,1._wp,s(ib:il,jb:jl),1,d(ib:il,jb:jl),1)
Florian Wilhelm's avatar
Florian Wilhelm committed
174
175
176
#endif

        ENDDO
Florian Wilhelm's avatar
Florian Wilhelm committed
177

Florian Wilhelm's avatar
Florian Wilhelm committed
178
        DEALLOCATE(d, r, s, Ad)
179

180
181
        ! Communicate x (only d was used the whole time)
        IF (config%do_exchange) CALL exchange(x, 'precond_cg_method 2')
Florian Wilhelm's avatar
Florian Wilhelm committed
182

Florian Wilhelm's avatar
Florian Wilhelm committed
183
    END FUNCTION precond_cg_method
Florian Wilhelm's avatar
Florian Wilhelm committed
184

Florian Wilhelm's avatar
Florian Wilhelm committed
185
186
    ! CG-Method with A*x=b and startvalue x
    FUNCTION cg_method(A, b, x, ext_x, exchange, tol_opt, maxiter_opt) RESULT(kiter)
187
        USE preconditioners, ONLY: identity
Florian Wilhelm's avatar
Florian Wilhelm committed
188

Florian Wilhelm's avatar
Florian Wilhelm committed
189
190
191
192
193
194
        REAL(wp), INTENT(IN) :: b(:,:)                  ! right-hand-side b
        REAL(wp), INTENT(INOUT) :: x(:,:)               ! startvalue and result
        TYPE(extent_type), INTENT(IN) :: ext_x(:)       ! extent of x
        REAL(wp), OPTIONAL, INTENT(IN) :: tol_opt       ! tolerance for residual
        INTEGER, OPTIONAL, INTENT(IN) :: maxiter_opt    ! maximum iterations
        INTEGER :: kiter                                ! number of iterations
Florian Wilhelm's avatar
Florian Wilhelm committed
195

Florian Wilhelm's avatar
Florian Wilhelm committed
196
197
        INTERFACE
            ! Matrix-Vector multiplication of the linear system to solve
Florian Wilhelm's avatar
Florian Wilhelm committed
198
            SUBROUTINE A(field, res_field)
Florian Wilhelm's avatar
Florian Wilhelm committed
199
200
201
202
203
204
205
206
207
208
209
                USE mo_kind, ONLY: wp, dp
                REAL(wp), INTENT(IN) :: field(:,:)
                REAL(wp), INTENT(OUT) :: res_field(:,:)
            END SUBROUTINE
            ! Function to exchange boundaries if necessary
            SUBROUTINE exchange(a0, text)
                USE mo_kind, ONLY: wp, dp
                REAL(wp), INTENT(INOUT) :: a0(:,:)
                CHARACTER (LEN=*), INTENT(IN), OPTIONAL :: text
            END SUBROUTINE exchange
        END INTERFACE
Florian Wilhelm's avatar
Florian Wilhelm committed
210

Florian Wilhelm's avatar
Florian Wilhelm committed
211
212
213
214
215
216
217
218
219
        IF ( PRESENT(tol_opt) ) THEN
            IF ( PRESENT(maxiter_opt) ) THEN
                kiter = precond_cg_method(A, b, x, ext_x, identity, exchange, tol_opt, maxiter_opt)
            ELSE
                kiter = precond_cg_method(A, b, x, ext_x, identity, exchange, tol_opt)
            ENDIF
        ELSE
            kiter = precond_cg_method(A, b, x, ext_x, identity, exchange)
        ENDIF
Florian Wilhelm's avatar
Florian Wilhelm committed
220

Florian Wilhelm's avatar
Florian Wilhelm committed
221
    END FUNCTION cg_method
Florian Wilhelm's avatar
Florian Wilhelm committed
222

223
    ! Additive Schwarz method which uses a local solver for each partition
Florian Wilhelm's avatar
Florian Wilhelm committed
224
225
226
227
228
229
    FUNCTION schwarz_method(A, b, x, ext_x, local_solver, exchange, tol_opt, maxiter_opt) RESULT(kiter)
        REAL(wp), INTENT(IN) :: b(:,:)                  ! right-hand-side b
        REAL(wp), INTENT(INOUT) :: x(:,:)               ! startvalue and result
        TYPE(extent_type), INTENT(IN) :: ext_x(:)       ! extent of x
        REAL(wp), OPTIONAL, INTENT(IN) :: tol_opt       ! tolerance for residual
        INTEGER, OPTIONAL, INTENT(IN) :: maxiter_opt    ! maximum iterations
230
231
        INTEGER :: kiter, maxiter, iiter, ie, je, ierror
        !INTEGER :: my_rank
232
        REAL(wp) :: tol, numer, denom, local_relres
Florian Wilhelm's avatar
Florian Wilhelm committed
233

Florian Wilhelm's avatar
Florian Wilhelm committed
234
        INTERFACE
235
            ! Matrix-Vector multiplication of the linear system to solve
Florian Wilhelm's avatar
Florian Wilhelm committed
236
237
238
239
240
            SUBROUTINE A(field, res_field)
                USE mo_kind, ONLY: wp, dp
                REAL(wp), INTENT(IN) :: field(:,:)
                REAL(wp), INTENT(OUT) :: res_field(:,:)
            END SUBROUTINE
241
            ! The solver used by Schwarz method
242
            FUNCTION local_solver(A, b, x, ext_x, exchange, tol_opt, maxiter_opt) RESULT(kiter)
Florian Wilhelm's avatar
Florian Wilhelm committed
243
                USE mo_kind, ONLY: wp, dp
244
                USE ppm_extents, ONLY: extent_type => extent
Florian Wilhelm's avatar
Florian Wilhelm committed
245

246
247
248
249
250
251
                REAL(wp), INTENT(IN) :: b(:,:)                  ! right-hand-side b
                REAL(wp), INTENT(INOUT) :: x(:,:)               ! startvalue and result
                TYPE(extent_type), INTENT(IN) :: ext_x(:)       ! extent of x
                REAL(wp), OPTIONAL, INTENT(IN) :: tol_opt       ! tolerance for residual
                INTEGER, OPTIONAL, INTENT(IN) :: maxiter_opt    ! maximum iterations
                INTEGER :: kiter                                ! needed iterations
Florian Wilhelm's avatar
Florian Wilhelm committed
252
                INTERFACE
253
                    ! Matrix-Vector multiplication of the linear system to solve
Florian Wilhelm's avatar
Florian Wilhelm committed
254
                    SUBROUTINE A(field, res_field)
Florian Wilhelm's avatar
Florian Wilhelm committed
255
256
257
258
                        USE mo_kind, ONLY: wp, dp
                        REAL(wp), INTENT(IN) :: field(:,:)
                        REAL(wp), INTENT(OUT) :: res_field(:,:)
                    END SUBROUTINE
259
260
261
262
263
264
                    ! Function to exchange boundaries if necessary
                    SUBROUTINE exchange(a0, text)
                        USE mo_kind, ONLY: wp, dp
                        REAL(wp), INTENT(INOUT) :: a0(:,:)
                        CHARACTER (LEN=*), INTENT(IN), OPTIONAL :: text
                    END SUBROUTINE exchange
Florian Wilhelm's avatar
Florian Wilhelm committed
265
266
                END INTERFACE
            END FUNCTION
267
            ! Function to exchange boundaries if necessary
268
            SUBROUTINE exchange(a0, text)
269
                USE mo_kind, ONLY: wp, dp
270
271
272
                REAL(wp), INTENT(INOUT) :: a0(:,:)
                CHARACTER (LEN=*), INTENT(IN), OPTIONAL :: text
            END SUBROUTINE exchange
Florian Wilhelm's avatar
Florian Wilhelm committed
273
274
        END INTERFACE

Florian Wilhelm's avatar
Florian Wilhelm committed
275
276
277
278
279
280
281
282
283
284
        IF ( PRESENT(maxiter_opt) ) THEN
            maxiter = maxiter_opt
        ELSE
            maxiter = config%schwarz_maxiter
        ENDIF
        IF ( PRESENT(tol_opt) ) THEN
            tol = tol_opt
        ELSE
            tol = config%schwarz_tol
        ENDIF
Florian Wilhelm's avatar
Florian Wilhelm committed
285
286

        kiter = 0
Florian Wilhelm's avatar
Florian Wilhelm committed
287
288
289
        ie = SIZE(b,1)
        je = SIZE(b,2)

290
291
        denom = calc_abs_res(A, b, x, ext_x, .TRUE.)
        numer = denom
Florian Wilhelm's avatar
Florian Wilhelm committed
292

293
        !        CALL MPI_Comm_rank(mpi_comm_world, my_rank, ierror)
Florian Wilhelm's avatar
Florian Wilhelm committed
294

295
        DO WHILE (kiter < maxiter .AND. numer > tol*denom)
Florian Wilhelm's avatar
Florian Wilhelm committed
296
            kiter = kiter + 1
Florian Wilhelm's avatar
Florian Wilhelm committed
297

298
299
            ! Solve local problem according to solver as configured
            iiter = local_solver(A, b, x, ext_x, exchange, config%tol, config%maxiter)
Florian Wilhelm's avatar
Florian Wilhelm committed
300

301
            CALL exchange(x, 'schwarz_method 1')
Florian Wilhelm's avatar
Florian Wilhelm committed
302

303
304
!~             local_relres = calc_rel_res(A, b, x, ext_x, .FALSE.)
!~             WRITE(*,*) "CPU: ", my_rank, "kiter: ", kiter, "iiter: ", iiter, "local_relres: ", local_relres
Florian Wilhelm's avatar
Florian Wilhelm committed
305
306
307
!~             IF( my_rank == 0 ) WRITE(*,*) "Global relres", kiter, ":", numer/denom

            IF( kiter >= config%schwarz_miniter .AND. MOD(kiter, config%schwarz_checkrate) == 0 ) THEN
308
309
                numer = calc_abs_res(A, b, x, ext_x, .TRUE.)
            ENDIF
Florian Wilhelm's avatar
Florian Wilhelm committed
310
        END DO
Florian Wilhelm's avatar
Florian Wilhelm committed
311

Florian Wilhelm's avatar
Florian Wilhelm committed
312
    END FUNCTION schwarz_method
Florian Wilhelm's avatar
Florian Wilhelm committed
313

314
    ! Non-preconditioned Chebyshev method
Florian Wilhelm's avatar
Florian Wilhelm committed
315
    FUNCTION chebyshev_method(A, b, x, ext_x, lambda_min, lambda_max, exchange, tol_opt, maxiter_opt) RESULT(kiter)
316
        USE preconditioners, ONLY: identity
Florian Wilhelm's avatar
Florian Wilhelm committed
317

318
319
320
321
322
323
324
        REAL(wp), INTENT(IN) :: b(:,:)                  ! right-hand-side b
        REAL(wp), INTENT(INOUT) :: x(:,:)               ! startvalue and result
        TYPE(extent_type), INTENT(IN) :: ext_x(:)       ! extent of x
        REAL(wp), INTENT(IN) :: lambda_min, lambda_max  ! smallest and largest absolute eigenvalue
        REAL(wp), OPTIONAL, INTENT(IN) :: tol_opt       ! tolerance for residual
        INTEGER, OPTIONAL, INTENT(IN) :: maxiter_opt    ! maximum iterations
        INTEGER :: kiter                                ! number of iterations
Florian Wilhelm's avatar
Florian Wilhelm committed
325

326
327
        INTERFACE
            ! Matrix-Vector multiplication of the linear system to solve
Florian Wilhelm's avatar
Florian Wilhelm committed
328
            SUBROUTINE A(field, res_field)
329
                USE mo_kind, ONLY: wp, dp
330
331
332
333
334
                REAL(wp), INTENT(IN) :: field(:,:)
                REAL(wp), INTENT(OUT) :: res_field(:,:)
            END SUBROUTINE
            ! Function to exchange boundaries if necessary
            SUBROUTINE exchange(a0, text)
335
                USE mo_kind, ONLY: wp, dp
336
337
338
339
                REAL(wp), INTENT(INOUT) :: a0(:,:)
                CHARACTER (LEN=*), INTENT(IN), OPTIONAL :: text
            END SUBROUTINE exchange
        END INTERFACE
Florian Wilhelm's avatar
Florian Wilhelm committed
340

341
342
343
344
345
346
347
348
349
350
351
        IF ( PRESENT(tol_opt) ) THEN
            IF ( PRESENT(maxiter_opt) ) THEN
                kiter = precond_chebyshev_method(A, b, x, ext_x, lambda_min, lambda_max, identity, exchange, tol_opt, maxiter_opt)
            ELSE
                kiter = precond_chebyshev_method(A, b, x, ext_x, lambda_min, lambda_max, identity, exchange, tol_opt)
            ENDIF
        ELSE
            kiter = precond_chebyshev_method(A, b, x, ext_x, lambda_min, lambda_max, identity, exchange)
        ENDIF

    END FUNCTION chebyshev_method
Florian Wilhelm's avatar
Florian Wilhelm committed
352

353
    ! Preconditioned Chebyshev acceleration according to Saad "Iterative Methods for Sparse Linear Systems", Algorithm 12.1
Florian Wilhelm's avatar
Florian Wilhelm committed
354
    FUNCTION precond_chebyshev_method(A, b, x, ext_x, lambda_min, lambda_max, precond, exchange, tol_opt, maxiter_opt) RESULT(kiter)
Florian Wilhelm's avatar
Florian Wilhelm committed
355
356
357
        REAL(wp), INTENT(IN) :: b(:,:)                  ! right-hand-side b
        REAL(wp), INTENT(INOUT) :: x(:,:)               ! startvalue and result
        TYPE(extent_type), INTENT(IN) :: ext_x(:)       ! extent of x
358
        REAL(wp), INTENT(IN) :: lambda_min, lambda_max  ! smallest and largest absolute eigenvalue
Florian Wilhelm's avatar
Florian Wilhelm committed
359
360
        REAL(wp), OPTIONAL, INTENT(IN) :: tol_opt       ! tolerance for residual
        INTEGER, OPTIONAL, INTENT(IN) :: maxiter_opt    ! maximum iterations
361
        INTEGER :: i, j, jb, ib, il, jl, ie, je, kiter, maxiter
362
        REAL(wp) :: rho_old, rho_new, theta, delta, sigma, tol, denom, numer
Florian Wilhelm's avatar
Florian Wilhelm committed
363
        ! Auxiliary fields for CG-Method
364
        REAL(wp), ALLOCATABLE :: r(:,:), Ax(:,:), d(:,:)
Florian Wilhelm's avatar
Florian Wilhelm committed
365

Florian Wilhelm's avatar
Florian Wilhelm committed
366
367
        INTERFACE
            ! Matrix-Vector multiplication of the linear system to solve
Florian Wilhelm's avatar
Florian Wilhelm committed
368
            SUBROUTINE A(field, res_field)
Florian Wilhelm's avatar
Florian Wilhelm committed
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
                USE mo_kind, ONLY: wp, dp
                REAL(wp), INTENT(IN) :: field(:,:)
                REAL(wp), INTENT(OUT) :: res_field(:,:)
            END SUBROUTINE
            ! Preconditioner
            SUBROUTINE precond(r)
                USE mo_kind, ONLY: wp, dp
                REAL(wp), INTENT(INOUT) :: r(:,:)
            END SUBROUTINE
            ! Function to exchange boundaries if necessary
            SUBROUTINE exchange(a0, text)
                USE mo_kind, ONLY: wp, dp
                REAL(wp), INTENT(INOUT) :: a0(:,:)
                CHARACTER (LEN=*), INTENT(IN), OPTIONAL :: text
            END SUBROUTINE exchange
        END INTERFACE

        ! Check and set optional arguments
        IF ( PRESENT(maxiter_opt) ) THEN
            maxiter = maxiter_opt
        ELSE
            maxiter = config%maxiter
        ENDIF
        IF ( PRESENT(tol_opt) ) THEN
            tol = tol_opt**2 ! **2, because we compare later the square of a norm
        ELSE
            tol = config%tol**2
        ENDIF

398
399
        ie = SIZE(x,1)
        je = SIZE(x,2)
Florian Wilhelm's avatar
Florian Wilhelm committed
400
401
402
403
404
        ib = extent_start(ext_x(1))
        jb = extent_start(ext_x(2))
        il = extent_end(ext_x(1))
        jl = extent_end(ext_x(2))

405
406
        theta = (lambda_max + lambda_min) / 2._wp
        delta = (lambda_max - lambda_min) / 2._wp
Florian Wilhelm's avatar
Florian Wilhelm committed
407

408
409
410
411
412
413
414
415
416
417
418
419
        ! Check for trivial solution
        IF (delta == 0._wp) THEN
            IF (lambda_max == 0._wp) THEN
                x(ib:il,jb:jl) = 0._wp
            ELSE
                x(ib:il,jb:jl) = b(ib:il,jb:jl)/lambda_max
            ENDIF
            RETURN
        ENDIF

        ALLOCATE(r(ie,je), Ax(ie,je), d(ie,je))

420
        ! Residual and direction for first iteration
Florian Wilhelm's avatar
Florian Wilhelm committed
421
422
        CALL A(x, Ax)
        r(ib:il,jb:jl) = b(ib:il,jb:jl) - Ax(ib:il,jb:jl)
423
424
        CALL precond(r)
        d(ib:il,jb:jl) = 1._wp/theta*r(ib:il,jb:jl)
Florian Wilhelm's avatar
Florian Wilhelm committed
425

426
        ! Denominator for relativ residual
Florian Wilhelm's avatar
Florian Wilhelm committed
427
        denom = arr_dotproduct(r(ib:il,jb:jl))
Florian Wilhelm's avatar
Florian Wilhelm committed
428
        numer = denom
Florian Wilhelm's avatar
Florian Wilhelm committed
429

430
        kiter = 0
431
432
        sigma = theta/delta
        rho_old = 1._wp/sigma
Florian Wilhelm's avatar
Florian Wilhelm committed
433

Florian Wilhelm's avatar
Florian Wilhelm committed
434
        DO WHILE (kiter < maxiter .AND. numer > tol*denom)
435
            kiter = kiter + 1
Florian Wilhelm's avatar
Florian Wilhelm committed
436

437
438
            ! x_new = x_old + d_old
            x(ib:il,jb:jl) = x(ib:il,jb:jl) + d(ib:il,jb:jl)
Florian Wilhelm's avatar
Florian Wilhelm committed
439
            IF (config%do_exchange) CALL exchange(x, 'chebyshev 1')
Florian Wilhelm's avatar
Florian Wilhelm committed
440

441
            ! r_new = r_old - A*d_old = b - A*x_new
Florian Wilhelm's avatar
Florian Wilhelm committed
442
            CALL A(x, Ax)
443
444
            r(ib:il,jb:jl) = b(ib:il,jb:jl) - Ax(ib:il,jb:jl) ! == r - Ad in Saad Algorithm
            CALL precond(r)
Florian Wilhelm's avatar
Florian Wilhelm committed
445

446
447
            ! update rho
            rho_new = 1._wp/(2._wp*sigma - rho_old)
Florian Wilhelm's avatar
Florian Wilhelm committed
448

449
450
451
            ! d_new = rho_new*rho_old * d_old + (2*rho_new)/delta * r_new
            d(ib:il,jb:jl) = rho_new*rho_old*d(ib:il,jb:jl) + (2*rho_old)/delta*r(ib:il,jb:jl)
            rho_old = rho_new
Florian Wilhelm's avatar
Florian Wilhelm committed
452

453
            ! calculate current residual
Florian Wilhelm's avatar
Florian Wilhelm committed
454
            IF( kiter >= config%cheby_miniter .AND. MOD(kiter, config%cheby_checkrate) == 0 ) THEN
455
                numer = arr_dotproduct(r(ib:il,jb:jl))
Florian Wilhelm's avatar
Florian Wilhelm committed
456
            ENDIF
Florian Wilhelm's avatar
Florian Wilhelm committed
457

458
            !WRITE(*,*) 'CHEBY: kiter, ', kiter, 'rel_res', numer/denom
Florian Wilhelm's avatar
Florian Wilhelm committed
459
        ENDDO
Florian Wilhelm's avatar
Florian Wilhelm committed
460

461
        DEALLOCATE(r, Ax, d)
462

463
    END FUNCTION precond_chebyshev_method
Florian Wilhelm's avatar
Florian Wilhelm committed
464

465
END MODULE solvers