linear_algebra.f90 7.61 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
! linear_algebra.f90 --- basic tools from linear algebra
!
! Copyright  (C)  2010  Florian Wilhelm <Florian.Wilhelm@kit.edu>
!
! Version: 1.0
! Keywords: scales ppm solver linear algebra residual
! 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 linear_algebra
44
    USE solver_internal, ONLY: config
45
    USE mo_kind, ONLY: wp, dp
46
47
    USE ppm_extents, ONLY: extent_type => extent, extent_start, &
         extent_end, extent_size
48
    USE ppm_base, ONLY: abort_ppm
Florian Wilhelm's avatar
Florian Wilhelm committed
49

50
51
52
    IMPLICIT NONE

    PRIVATE
Florian Wilhelm's avatar
Florian Wilhelm committed
53

54
55
56
    !-------------------------------------------------------------------
    ! Module Function & Procedure Interface
    !-------------------------------------------------------------------
Florian Wilhelm's avatar
Florian Wilhelm committed
57

58
59
60
    INTERFACE arr_dotproduct
        MODULE PROCEDURE arr_dotproduct1
        MODULE PROCEDURE arr_dotproduct2
61
62
    END INTERFACE

63
    PUBLIC:: calc_abs_res, calc_rel_res, arr_dotproduct, global_sum, arr_norm_2
Florian Wilhelm's avatar
Florian Wilhelm committed
64

65
    CONTAINS
Florian Wilhelm's avatar
Florian Wilhelm committed
66

67
     ! Array-wise dotproduct for one matrix with itself
68
    FUNCTION arr_dotproduct1(x, global_opt) RESULT(ans)
69
70
71
72
73
74
75
76
77
78
79
        REAL(wp), INTENT(IN) :: x(:,:)
        LOGICAL, INTENT(IN), OPTIONAL :: global_opt
        REAL(wp) :: ans
        LOGICAL :: global

        IF ( PRESENT(global_opt) ) THEN
            global = global_opt
        ELSE
            global = config%do_exchange
        ENDIF

Florian Wilhelm's avatar
Florian Wilhelm committed
80
        ans = SUM(x**2)
81
82

        IF (global) ans = global_sum(ans)
Florian Wilhelm's avatar
Florian Wilhelm committed
83

84
    END FUNCTION arr_dotproduct1
Florian Wilhelm's avatar
Florian Wilhelm committed
85

86
    ! Array-wise dotproduct for two different matrices
87
88
89
    FUNCTION arr_dotproduct2(x, y, global_opt) RESULT(ans)
#ifdef BLAS
        REAL(wp) :: DDOT
Florian Wilhelm's avatar
Florian Wilhelm committed
90
#endif
91
        REAL(wp), INTENT(IN) :: x(:,:), y(:,:)
92
93
        LOGICAL, INTENT(IN), OPTIONAL :: global_opt
        REAL(wp) :: ans
94
        REAL(wp), ALLOCATABLE :: xy(:)
95
        LOGICAL :: global
96
        INTEGER :: j, ie, je
97

98
99
100
101
102
103
        IF ( PRESENT(global_opt) ) THEN
            global = global_opt
        ELSE
            global = config%do_exchange
        ENDIF

104
105
        ie = SIZE(x,1)
        je = SIZE(x,2)
Florian Wilhelm's avatar
Florian Wilhelm committed
106

107
#ifndef BLAS
108
        ALLOCATE(xy(je))
109

110
111
        DO j=1,je
            xy(j) = DOT_PRODUCT(x(:,j), y(:,j))
112
        ENDDO
Florian Wilhelm's avatar
Florian Wilhelm committed
113

114
        ans = SUM(xy)
115
        DEALLOCATE(xy)
Florian Wilhelm's avatar
Florian Wilhelm committed
116
117
#else
        ans = DDOT(ie*je, x, 1, y, 1)
118
#endif
119
        IF (global) ans = global_sum(ans)
Florian Wilhelm's avatar
Florian Wilhelm committed
120

121
    END FUNCTION arr_dotproduct2
Florian Wilhelm's avatar
Florian Wilhelm committed
122

123
    ! Array-wise 2-norm
124
    FUNCTION arr_norm_2(x, global_opt) RESULT(ans)
125
        REAL(wp), INTENT(IN) :: x(:,:)
126
127
128
        LOGICAL, INTENT(IN), OPTIONAL :: global_opt
        REAL(wp) :: ans
        LOGICAL :: global
129

130
131
132
133
134
135
        IF ( PRESENT(global_opt) ) THEN
            global = global_opt
        ELSE
            global = config%do_exchange
        ENDIF

136
        ans = SQRT(arr_dotproduct(x, global))
137

138
    END FUNCTION arr_norm_2
Florian Wilhelm's avatar
Florian Wilhelm committed
139

140
    ! Sums a variable globally up
141
    FUNCTION global_sum(summand, comm_opt) RESULT(all_sum)
142
#ifdef USE_MPI
143
        USE solver_internal, ONLY: mpi_real_dp, mpi_comm_world, MPI_SUM
144
#endif
145

146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
        REAL(wp), INTENT(in) :: summand
        INTEGER, OPTIONAL, INTENT(in) :: comm_opt
        REAL(wp) :: all_sum
#ifdef USE_MPI
        INTEGER :: comm, ierror

        IF (PRESENT(comm_opt)) THEN
            comm = comm_opt
        ELSE
            comm = mpi_comm_world
        ENDIF

        CALL MPI_ALLREDUCE (summand, all_sum, 1, mpi_real_dp, MPI_SUM, comm, ierror)
#else
        all_sum = summand
#endif

163
    END FUNCTION global_sum
164
165
166
167
168
169
170
171
172
173

    ! Calculate absolute residual of system
    ! If global_opt is .TRUE. then calculate the global absolute residual
    FUNCTION calc_abs_res(A, b, x, ext_x, global_opt) RESULT(abs_res)
        REAL(wp), INTENT(IN) :: x(:,:)
        REAL(wp), INTENT(IN) :: b(:,:)
        TYPE(extent_type), INTENT(IN) :: ext_x(:)      ! extent of x
        LOGICAL, INTENT(IN), OPTIONAL :: global_opt
        REAL(wp), ALLOCATABLE :: r(:,:), Ax(:,:)
        INTEGER :: ie, je, ib, jb, il, jl
174
        REAL(wp) :: abs_res
175
        LOGICAL :: global
Florian Wilhelm's avatar
Florian Wilhelm committed
176

177
178
        INTERFACE
            SUBROUTINE A(field, res_field)
179
                USE mo_kind, ONLY: wp, dp
180
181
182
183
                REAL(wp), INTENT(IN) :: field(:,:)
                REAL(wp), INTENT(OUT) :: res_field(:,:)
            END SUBROUTINE
        END INTERFACE
Florian Wilhelm's avatar
Florian Wilhelm committed
184

185
186
187
188
189
        IF ( PRESENT(global_opt) ) THEN
            global = global_opt
        ELSE
            global = config%do_exchange
        ENDIF
Florian Wilhelm's avatar
Florian Wilhelm committed
190

191
192
        ie = SIZE(x,1)
        je = SIZE(x,2)
Florian Wilhelm's avatar
Florian Wilhelm committed
193
194
195
196
197
        ib = extent_start(ext_x(1))
        jb = extent_start(ext_x(2))
        il = extent_end(ext_x(1))
        jl = extent_end(ext_x(2))

198
        ALLOCATE(r(ie,je), Ax(ie,je))
Florian Wilhelm's avatar
Florian Wilhelm committed
199

200
        CALL A(x, Ax)
201
        r(ib:il,jb:jl) = b(ib:il,jb:jl) - Ax(ib:il,jb:jl)
Florian Wilhelm's avatar
Florian Wilhelm committed
202

203
        abs_res = arr_norm_2(r(ib:il,jb:jl), global)
Florian Wilhelm's avatar
Florian Wilhelm committed
204

205
        DEALLOCATE(Ax, r)
Florian Wilhelm's avatar
Florian Wilhelm committed
206

207
    END FUNCTION calc_abs_res
Florian Wilhelm's avatar
Florian Wilhelm committed
208

209
210
    ! Calculate relative residual of system
    ! If global_opt is .TRUE. then calculate the global relative residual
211
    FUNCTION calc_rel_res(A, b, x, ext_x, global_opt, x0_opt) RESULT(rel_res)
212
213
214
215
216
217
218
        REAL(wp), INTENT(IN) :: x(:,:)
        REAL(wp), INTENT(IN) :: b(:,:)
        TYPE(extent_type), INTENT(IN) :: ext_x(:)      ! extent of x
        LOGICAL, INTENT(IN), OPTIONAL :: global_opt
        REAL(wp), INTENT(IN), OPTIONAL :: x0_opt(:,:)
        REAL(wp), ALLOCATABLE :: r(:,:), Ax(:,:), x0(:,:)
        INTEGER :: ie, je
219
        REAL(wp) :: rel_res, numer, denom
Florian Wilhelm's avatar
Florian Wilhelm committed
220

221
222
        INTERFACE
            SUBROUTINE A(field, res_field)
223
                USE mo_kind, ONLY: wp, dp
224
225
226
227
                REAL(wp), INTENT(IN) :: field(:,:)
                REAL(wp), INTENT(OUT) :: res_field(:,:)
            END SUBROUTINE
        END INTERFACE
Florian Wilhelm's avatar
Florian Wilhelm committed
228

229
230
        IF ( PRESENT(x0_opt) ) THEN
            x0 = x0_opt
Florian Wilhelm's avatar
Florian Wilhelm committed
231
        ELSE
232
233
234
235
236
            ie = SIZE(x,1)
            je = SIZE(x,2)
            ALLOCATE(x0(ie,je))
            x0 = 0._wp
        ENDIF
Florian Wilhelm's avatar
Florian Wilhelm committed
237

238
239
240
241
242
243
244
        IF ( PRESENT(global_opt) ) THEN
            numer = calc_abs_res(A, b, x, ext_x, global_opt)
            denom = calc_abs_res(A, b, x0, ext_x, global_opt)
        ELSE
            numer = calc_abs_res(A, b, x, ext_x)
            denom = calc_abs_res(A, b, x0, ext_x)
        ENDIF
Florian Wilhelm's avatar
Florian Wilhelm committed
245

246
        rel_res = numer/denom
Florian Wilhelm's avatar
Florian Wilhelm committed
247

248
249
250
251
252
253
254
        IF (.NOT. PRESENT(x0_opt)) THEN
            DEALLOCATE(x0)
        ENDIF

    END FUNCTION calc_rel_res

END MODULE linear_algebra