← Back to branch summary

~nickpapior/siesta/mixing

« back to all changes in this revision

Viewing changes to Src/SiestaXC/array.F90

  • Committer: Nick Papior
  • Date: 2017-09-08 07:49:32 UTC
  • mfrom: (591.1.41 trunk)
  • Revision ID: nickpapior@gmail.com-20170908074932-l6hiee82h1vjnn2h
Merged trunk-r594-632

Show diffs side-by-side

added added

removed removed

Lines of Context:
Src/SiestaXC/array.F90
 
1
!!@LICENSE
 
2
!
 
3
!******************************************************************************
 
4
! MODULE array
 
5
! Provides utilities for copying data from one array shape to another
 
6
! array shape.
 
7
! Written by Nick R. Papior, Feb.2017
 
8
!******************************************************************************
 
9
!
 
10
!   PUBLIC procedures available from this module:
 
11
! array_copy    : Copies an array from one dimension to another array with another dimension
 
12
! array_add     : Adds an array from one dimension to another array with another dimension
 
13
!
 
14
!   PUBLIC parameters, types, and variables available from this module:
 
15
! none
 
16
!
 
17
! The intent of this module is to bypass the use of the intrinsic reshape
 
18
! function.
 
19
! Sadly many compilers implement the reshape with a guard on pointers such
 
20
! that a temporary array is created upon reshaping. For very large segments
 
21
! this turns out to be a substantial memory consumption on the heap which
 
22
! may easily be circumvented by doing explicit copying.
 
23
! These routines are extremely simple and does nothing but copy data from
 
24
! one shape to another.
 
25
!
 
26
! Currently only these copies and adds are implemented:
 
27
!
 
28
!   1D -> 2D [integer, real(sp), real(dp)]
 
29
!   1D -> 3D [integer, real(sp), real(dp)]
 
30
!   2D -> 1D [integer, real(sp), real(dp)]
 
31
!   3D -> 1D [integer, real(sp), real(dp)]
 
32
 
 
33
module m_array
 
34
 
 
35
  use precision, only: sp, dp
 
36
  use sys,       only: die       ! Termination routine
 
37
 
 
38
  private
 
39
 
 
40
  public :: array_copy
 
41
  public :: array_add
 
42
 
 
43
  interface array_copy
 
44
     module procedure ac_1d_2d_ip, ac_1d_2d_sp, ac_1d_2d_dp
 
45
     module procedure ac_1d_3d_ip, ac_1d_3d_sp, ac_1d_3d_dp
 
46
     module procedure ac_2d_1d_ip, ac_2d_1d_sp, ac_2d_1d_dp
 
47
     module procedure ac_3d_1d_ip, ac_3d_1d_sp, ac_3d_1d_dp
 
48
     module procedure ac_4d_1d_ip, ac_4d_1d_sp, ac_4d_1d_dp
 
49
  end interface array_copy
 
50
 
 
51
  interface array_add
 
52
     module procedure aa_1d_2d_ip, aa_1d_2d_sp, aa_1d_2d_dp
 
53
     module procedure aa_1d_3d_ip, aa_1d_3d_sp, aa_1d_3d_dp
 
54
     module procedure aa_1d_4d_ip, aa_1d_4d_sp, aa_1d_4d_dp
 
55
     module procedure aa_2d_1d_ip, aa_2d_1d_sp, aa_2d_1d_dp
 
56
     module procedure aa_3d_1d_ip, aa_3d_1d_sp, aa_3d_1d_dp
 
57
     module procedure aa_4d_1d_ip, aa_4d_1d_sp, aa_4d_1d_dp
 
58
  end interface array_add
 
59
 
 
60
contains
 
61
 
 
62
  ! Copies a 1D array to a 2D array
 
63
  subroutine ac_1d_2d_ip(from_1D, to_1D, in1D, from_2D, to_2D, out2D)
 
64
    integer, intent(in) :: from_1D, to_1D
 
65
    integer, intent(in) :: in1D(:)
 
66
    integer, intent(in) :: from_2D(2), to_2D(2)
 
67
    integer, intent(inout) :: out2D(:,:)
 
68
 
 
69
    ! Local variables for copying
 
70
    integer :: i1D, i2D, j2D
 
71
 
 
72
    i2D = from_2D(1)
 
73
    j2D = from_2D(2)
 
74
    do i1D = from_1D, to_1D
 
75
       out2D(i2D, j2D) = in1D(i1D)
 
76
       i2D = i2D + 1
 
77
       if ( i2D > to_2D(1) ) then
 
78
          i2D = from_2D(1)
 
79
          j2D = j2D + 1
 
80
       end if
 
81
    end do
 
82
 
 
83
    if ( i2D /= from_2D(1) ) &
 
84
         call die("integer: 1D->2D failed (i)")
 
85
    if ( j2D <= to_2D(2) ) &
 
86
         call die("integer: 1D->2D failed (j)")
 
87
 
 
88
  end subroutine ac_1d_2d_ip
 
89
  subroutine ac_1d_2d_sp(from_1D, to_1D, in1D, from_2D, to_2D, out2D)
 
90
    integer, intent(in) :: from_1D, to_1D
 
91
    real(sp), intent(in) :: in1D(:)
 
92
    integer, intent(in) :: from_2D(2), to_2D(2)
 
93
    real(sp), intent(inout) :: out2D(:,:)
 
94
 
 
95
    ! Local variables for copying
 
96
    integer :: i1D, i2D, j2D
 
97
 
 
98
    i2D = from_2D(1)
 
99
    j2D = from_2D(2)
 
100
    do i1D = from_1D, to_1D
 
101
       out2D(i2D, j2D) = in1D(i1D)
 
102
       i2D = i2D + 1
 
103
       if ( i2D > to_2D(1) ) then
 
104
          i2D = from_2D(1)
 
105
          j2D = j2D + 1
 
106
       end if
 
107
    end do
 
108
 
 
109
    if ( i2D /= from_2D(1) ) &
 
110
         call die("real: 1D->2D failed (i)")
 
111
    if ( j2D <= to_2D(2) ) &
 
112
         call die("real: 1D->2D failed (j)")
 
113
 
 
114
  end subroutine ac_1d_2d_sp
 
115
  subroutine ac_1d_2d_dp(from_1D, to_1D, in1D, from_2D, to_2D, out2D)
 
116
    integer, intent(in) :: from_1D, to_1D
 
117
    real(dp), intent(in) :: in1D(:)
 
118
    integer, intent(in) :: from_2D(2), to_2D(2)
 
119
    real(dp), intent(inout) :: out2D(:,:)
 
120
 
 
121
    ! Local variables for copying
 
122
    integer :: i1D, i2D, j2D
 
123
 
 
124
    i2D = from_2D(1)
 
125
    j2D = from_2D(2)
 
126
    do i1D = from_1D, to_1D
 
127
       out2D(i2D, j2D) = in1D(i1D)
 
128
       i2D = i2D + 1
 
129
       if ( i2D > to_2D(1) ) then
 
130
          i2D = from_2D(1)
 
131
          j2D = j2D + 1
 
132
       end if
 
133
    end do
 
134
 
 
135
    if ( i2D /= from_2D(1) ) &
 
136
         call die("double: 1D->2D failed (i)")
 
137
    if ( j2D <= to_2D(2) ) &
 
138
         call die("double: 1D->2D failed (j)")
 
139
 
 
140
  end subroutine ac_1d_2d_dp
 
141
 
 
142
  ! Copies a 1D array to a 3D array
 
143
  subroutine ac_1d_3d_ip(from_1D, to_1D, in1D, from_3D, to_3D, out3D)
 
144
    integer, intent(in) :: from_1D, to_1D
 
145
    integer, intent(in) :: in1D(:)
 
146
    integer, intent(in) :: from_3D(3), to_3D(3)
 
147
    integer, intent(inout) :: out3D(:,:,:)
 
148
 
 
149
    ! Local variables for copying
 
150
    integer :: i1D, i3D, j3D, k3D
 
151
 
 
152
    i3D = from_3D(1)
 
153
    j3D = from_3D(2)
 
154
    k3D = from_3D(3)
 
155
    do i1D = from_1D, to_1D
 
156
       out3D(i3D, j3D, k3D) = in1D(i1D)
 
157
       i3D = i3D + 1
 
158
       if ( i3D > to_3D(1) ) then
 
159
          i3D = from_3D(1)
 
160
          j3D = j3D + 1
 
161
       end if
 
162
       if ( j3D > to_3D(2) ) then
 
163
          j3D = from_3D(2)
 
164
          k3D = k3D + 1
 
165
       end if
 
166
    end do
 
167
 
 
168
    if ( i3D /= from_3D(1) ) &
 
169
         call die("integer: 1D->3D failed (i)")
 
170
    if ( j3D /= from_3D(2) ) &
 
171
         call die("integer: 1D->3D failed (j)")
 
172
    if ( k3D <= to_3D(3) ) &
 
173
         call die("integer: 1D->3D failed (k)")
 
174
 
 
175
  end subroutine ac_1d_3d_ip
 
176
  subroutine ac_1d_3d_sp(from_1D, to_1D, in1D, from_3D, to_3D, out3D)
 
177
    integer, intent(in) :: from_1D, to_1D
 
178
    real(sp), intent(in) :: in1D(:)
 
179
    integer, intent(in) :: from_3D(3), to_3D(3)
 
180
    real(sp), intent(inout) :: out3D(:,:,:)
 
181
 
 
182
    ! Local variables for copying
 
183
    integer :: i1D, i3D, j3D, k3D
 
184
 
 
185
    i3D = from_3D(1)
 
186
    j3D = from_3D(2)
 
187
    k3D = from_3D(3)
 
188
    do i1D = from_1D, to_1D
 
189
       out3D(i3D, j3D, k3D) = in1D(i1D)
 
190
       i3D = i3D + 1
 
191
       if ( i3D > to_3D(1) ) then
 
192
          i3D = from_3D(1)
 
193
          j3D = j3D + 1
 
194
       end if
 
195
       if ( j3D > to_3D(2) ) then
 
196
          j3D = from_3D(2)
 
197
          k3D = k3D + 1
 
198
       end if
 
199
    end do
 
200
 
 
201
    if ( i3D /= from_3D(1) ) &
 
202
         call die("real: 1D->3D failed (i)")
 
203
    if ( j3D /= from_3D(2) ) &
 
204
         call die("real: 1D->3D failed (j)")
 
205
    if ( k3D <= to_3D(3) ) &
 
206
         call die("real: 1D->3D failed (k)")
 
207
 
 
208
  end subroutine ac_1d_3d_sp
 
209
  subroutine ac_1d_3d_dp(from_1D, to_1D, in1D, from_3D, to_3D, out3D)
 
210
    integer, intent(in) :: from_1D, to_1D
 
211
    real(dp), intent(in) :: in1D(:)
 
212
    integer, intent(in) :: from_3D(3), to_3D(3)
 
213
    real(dp), intent(inout) :: out3D(:,:,:)
 
214
 
 
215
    ! Local variables for copying
 
216
    integer :: i1D, i3D, j3D, k3D
 
217
 
 
218
    i3D = from_3D(1)
 
219
    j3D = from_3D(2)
 
220
    k3D = from_3D(3)
 
221
    do i1D = from_1D, to_1D
 
222
       out3D(i3D, j3D, k3D) = in1D(i1D)
 
223
       i3D = i3D + 1
 
224
       if ( i3D > to_3D(1) ) then
 
225
          i3D = from_3D(1)
 
226
          j3D = j3D + 1
 
227
       end if
 
228
       if ( j3D > to_3D(2) ) then
 
229
          j3D = from_3D(2)
 
230
          k3D = k3D + 1
 
231
       end if
 
232
    end do
 
233
 
 
234
    if ( i3D /= from_3D(1) ) &
 
235
         call die("double: 1D->3D failed (i)")
 
236
    if ( j3D /= from_3D(2) ) &
 
237
         call die("double: 1D->3D failed (j)")
 
238
    if ( k3D <= to_3D(3) ) &
 
239
         call die("double: 1D->3D failed (k)")
 
240
 
 
241
  end subroutine ac_1d_3d_dp
 
242
 
 
243
  ! Copies a 2D array to a 1D array
 
244
  subroutine ac_2d_1d_ip(from_2D, to_2D, in2D, from_1D, to_1D, out1D)
 
245
    integer, intent(in) :: from_2D(2), to_2D(2)
 
246
    integer, intent(in) :: in2D(:,:)
 
247
    integer, intent(in) :: from_1D, to_1D
 
248
    integer, intent(inout) :: out1D(:)
 
249
 
 
250
    ! Local variables for copying
 
251
    integer :: i2D, j2D, i1D
 
252
 
 
253
    i1D = from_1D
 
254
    do j2D = from_2D(2), to_2D(2)
 
255
       do i2D = from_2D(1), to_2D(1)
 
256
          out1D(i1D) = in2D(i2D, j2D)
 
257
          i1D = i1D + 1
 
258
       end do
 
259
    end do
 
260
 
 
261
    if ( i1D <= to_1D ) &
 
262
         call die("integer: 2D->1D failed")
 
263
 
 
264
  end subroutine ac_2d_1d_ip
 
265
  subroutine ac_2d_1d_sp(from_2D, to_2D, in2D, from_1D, to_1D, out1D)
 
266
    integer, intent(in) :: from_2D(2), to_2D(2)
 
267
    real(sp), intent(in) :: in2D(:,:)
 
268
    integer, intent(in) :: from_1D, to_1D
 
269
    real(sp), intent(inout) :: out1D(:)
 
270
 
 
271
    ! Local variables for copying
 
272
    integer :: i2D, j2D, i1D
 
273
 
 
274
    i1D = from_1D
 
275
    do j2D = from_2D(2), to_2D(2)
 
276
       do i2D = from_2D(1), to_2D(1)
 
277
          out1D(i1D) = in2D(i2D, j2D)
 
278
          i1D = i1D + 1
 
279
       end do
 
280
    end do
 
281
 
 
282
    if ( i1D <= to_1D ) &
 
283
         call die("real: 2D->1D failed")
 
284
 
 
285
  end subroutine ac_2d_1d_sp
 
286
  subroutine ac_2d_1d_dp(from_2D, to_2D, in2D, from_1D, to_1D, out1D)
 
287
    integer, intent(in) :: from_2D(2), to_2D(2)
 
288
    real(dp), intent(in) :: in2D(:,:)
 
289
    integer, intent(in) :: from_1D, to_1D
 
290
    real(dp), intent(inout) :: out1D(:)
 
291
 
 
292
    ! Local variables for copying
 
293
    integer :: i2D, j2D, i1D
 
294
 
 
295
    i1D = from_1D
 
296
    do j2D = from_2D(2), to_2D(2)
 
297
       do i2D = from_2D(1), to_2D(1)
 
298
          out1D(i1D) = in2D(i2D, j2D)
 
299
          i1D = i1D + 1
 
300
       end do
 
301
    end do
 
302
 
 
303
    if ( i1D <= to_1D ) &
 
304
         call die("double: 2D->1D failed")
 
305
 
 
306
  end subroutine ac_2d_1d_dp
 
307
  
 
308
  ! Copies a 3D array to a 1D array
 
309
  subroutine ac_3d_1d_ip(from_3D, to_3D, in3D, from_1D, to_1D, out1D)
 
310
    integer, intent(in) :: from_3D(3), to_3D(3)
 
311
    integer, intent(in) :: in3D(:,:,:)
 
312
    integer, intent(in) :: from_1D, to_1D
 
313
    integer, intent(inout) :: out1D(:)
 
314
 
 
315
    ! Local variables for copying
 
316
    integer :: i3D, j3D, k3D, i1D
 
317
 
 
318
    i1D = from_1D
 
319
    do k3D = from_3D(3), to_3D(3)
 
320
       do j3D = from_3D(2), to_3D(2)
 
321
          do i3D = from_3D(1), to_3D(1)
 
322
             out1D(i1D) = in3D(i3D, j3D, k3D)
 
323
             i1D = i1D + 1
 
324
          end do
 
325
       end do
 
326
    end do
 
327
 
 
328
    if ( i1D <= to_1D ) &
 
329
         call die("integer: 3D->1D failed")
 
330
 
 
331
  end subroutine ac_3d_1d_ip
 
332
  subroutine ac_3d_1d_sp(from_3D, to_3D, in3D, from_1D, to_1D, out1D)
 
333
    integer, intent(in) :: from_3D(3), to_3D(3)
 
334
    real(sp), intent(in) :: in3D(:,:,:)
 
335
    integer, intent(in) :: from_1D, to_1D
 
336
    real(sp), intent(inout) :: out1D(:)
 
337
 
 
338
    ! Local variables for copying
 
339
    integer :: i3D, j3D, k3D, i1D
 
340
 
 
341
    i1D = from_1D
 
342
    do k3D = from_3D(3), to_3D(3)
 
343
       do j3D = from_3D(2), to_3D(2)
 
344
          do i3D = from_3D(1), to_3D(1)
 
345
             out1D(i1D) = in3D(i3D, j3D, k3D)
 
346
             i1D = i1D + 1
 
347
          end do
 
348
       end do
 
349
    end do
 
350
 
 
351
    if ( i1D <= to_1D ) &
 
352
         call die("real: 3D->1D failed")
 
353
 
 
354
  end subroutine ac_3d_1d_sp
 
355
  subroutine ac_3d_1d_dp(from_3D, to_3D, in3D, from_1D, to_1D, out1D)
 
356
    integer, intent(in) :: from_3D(3), to_3D(3)
 
357
    real(dp), intent(in) :: in3D(:,:,:)
 
358
    integer, intent(in) :: from_1D, to_1D
 
359
    real(dp), intent(inout) :: out1D(:)
 
360
 
 
361
    ! Local variables for copying
 
362
    integer :: i3D, j3D, k3D, i1D
 
363
 
 
364
    i1D = from_1D
 
365
    do k3D = from_3D(3), to_3D(3)
 
366
       do j3D = from_3D(2), to_3D(2)
 
367
          do i3D = from_3D(1), to_3D(1)
 
368
             out1D(i1D) = in3D(i3D, j3D, k3D)
 
369
             i1D = i1D + 1
 
370
          end do
 
371
       end do
 
372
    end do
 
373
 
 
374
    if ( i1D <= to_1D ) &
 
375
         call die("double: 3D->1D failed")
 
376
 
 
377
  end subroutine ac_3d_1d_dp
 
378
 
 
379
  ! Copies a 4D array to a 1D array
 
380
  subroutine ac_4d_1d_ip(from_4D, to_4D, in4D, from_1D, to_1D, out1D)
 
381
    integer, intent(in) :: from_4D(4), to_4D(4)
 
382
    integer, intent(in) :: in4D(:,:,:,:)
 
383
    integer, intent(in) :: from_1D, to_1D
 
384
    integer, intent(inout) :: out1D(:)
 
385
 
 
386
    ! Local variables for copying
 
387
    integer :: i4D, j4D, k4D, m4D, i1D
 
388
 
 
389
    i1D = from_1D
 
390
    do m4D = from_4D(4), to_4D(4)
 
391
       do k4D = from_4D(3), to_4D(3)
 
392
          do j4D = from_4D(2), to_4D(2)
 
393
             do i4D = from_4D(1), to_4D(1)
 
394
                out1D(i1D) = in4D(i4D, j4D, k4D, m4D)
 
395
                i1D = i1D + 1
 
396
             end do
 
397
          end do
 
398
       end do
 
399
    end do
 
400
 
 
401
    if ( i1D <= to_1D ) &
 
402
         call die("integer: 4D->1D failed")
 
403
 
 
404
  end subroutine ac_4d_1d_ip
 
405
  subroutine ac_4d_1d_sp(from_4D, to_4D, in4D, from_1D, to_1D, out1D)
 
406
    integer, intent(in) :: from_4D(4), to_4D(4)
 
407
    real(sp), intent(in) :: in4D(:,:,:,:)
 
408
    integer, intent(in) :: from_1D, to_1D
 
409
    real(sp), intent(inout) :: out1D(:)
 
410
 
 
411
    ! Local variables for copying
 
412
    integer :: i4D, j4D, k4D, m4D, i1D
 
413
 
 
414
    i1D = from_1D
 
415
    do m4D = from_4D(4), to_4D(4)
 
416
       do k4D = from_4D(3), to_4D(3)
 
417
          do j4D = from_4D(2), to_4D(2)
 
418
             do i4D = from_4D(1), to_4D(1)
 
419
                out1D(i1D) = in4D(i4D, j4D, k4D, m4D)
 
420
                i1D = i1D + 1
 
421
             end do
 
422
          end do
 
423
       end do
 
424
    end do
 
425
    
 
426
    if ( i1D <= to_1D ) &
 
427
         call die("real: 4D->1D failed")
 
428
 
 
429
  end subroutine ac_4d_1d_sp
 
430
  subroutine ac_4d_1d_dp(from_4D, to_4D, in4D, from_1D, to_1D, out1D)
 
431
    integer, intent(in) :: from_4D(4), to_4D(4)
 
432
    real(dp), intent(in) :: in4D(:,:,:,:)
 
433
    integer, intent(in) :: from_1D, to_1D
 
434
    real(dp), intent(inout) :: out1D(:)
 
435
 
 
436
    ! Local variables for copying
 
437
    integer :: i4D, j4D, k4D, m4D, i1D
 
438
 
 
439
    i1D = from_1D
 
440
    do m4D = from_4D(4), to_4D(4)
 
441
       do k4D = from_4D(3), to_4D(3)
 
442
          do j4D = from_4D(2), to_4D(2)
 
443
             do i4D = from_4D(1), to_4D(1)
 
444
                out1D(i1D) = in4D(i4D, j4D, k4D, m4D)
 
445
                i1D = i1D + 1
 
446
             end do
 
447
          end do
 
448
       end do
 
449
    end do
 
450
    
 
451
    if ( i1D <= to_1D ) &
 
452
         call die("double: 4D->1D failed")
 
453
 
 
454
  end subroutine ac_4d_1d_dp
 
455
  
 
456
 
 
457
  ! Adds a 1D array to a 2D array
 
458
  subroutine aa_1d_2d_ip(from_1D, to_1D, a1D, from_2D, to_2D, out2D)
 
459
    integer, intent(in) :: from_1D, to_1D
 
460
    integer, intent(in) :: a1D(:)
 
461
    integer, intent(in) :: from_2D(2), to_2D(2)
 
462
    integer, intent(inout) :: out2D(:,:)
 
463
 
 
464
    ! Local variables for copying
 
465
    integer :: i1D, i2D, j2D
 
466
 
 
467
    i2D = from_2D(1)
 
468
    j2D = from_2D(2)
 
469
    do i1D = from_1D, to_1D
 
470
       out2D(i2D, j2D) = out2D(i2D, j2D) + a1D(i1D)
 
471
       i2D = i2D + 1
 
472
       if ( i2D > to_2D(1) ) then
 
473
          i2D = from_2D(1)
 
474
          j2D = j2D + 1
 
475
       end if
 
476
    end do
 
477
 
 
478
    if ( i2D /= from_2D(1) ) &
 
479
         call die("integer: 1D+>2D failed (i)")
 
480
    if ( j2D <= to_2D(2) ) &
 
481
         call die("integer: 1D+>2D failed (j)")
 
482
 
 
483
  end subroutine aa_1d_2d_ip
 
484
  subroutine aa_1d_2d_sp(from_1D, to_1D, a1D, from_2D, to_2D, out2D)
 
485
    integer, intent(in) :: from_1D, to_1D
 
486
    real(sp), intent(in) :: a1D(:)
 
487
    integer, intent(in) :: from_2D(2), to_2D(2)
 
488
    real(sp), intent(inout) :: out2D(:,:)
 
489
 
 
490
    ! Local variables for copying
 
491
    integer :: i1D, i2D, j2D
 
492
 
 
493
    i2D = from_2D(1)
 
494
    j2D = from_2D(2)
 
495
    do i1D = from_1D, to_1D
 
496
       out2D(i2D, j2D) = out2D(i2D, j2D) + a1D(i1D)
 
497
       i2D = i2D + 1
 
498
       if ( i2D > to_2D(1) ) then
 
499
          i2D = from_2D(1)
 
500
          j2D = j2D + 1
 
501
       end if
 
502
    end do
 
503
 
 
504
    if ( i2D /= from_2D(1) ) &
 
505
         call die("real: 1D+>2D failed (i)")
 
506
    if ( j2D <= to_2D(2) ) &
 
507
         call die("real: 1D+>2D failed (j)")
 
508
 
 
509
  end subroutine aa_1d_2d_sp
 
510
  subroutine aa_1d_2d_dp(from_1D, to_1D, a1D, from_2D, to_2D, out2D)
 
511
    integer, intent(in) :: from_1D, to_1D
 
512
    real(dp), intent(in) :: a1D(:)
 
513
    integer, intent(in) :: from_2D(2), to_2D(2)
 
514
    real(dp), intent(inout) :: out2D(:,:)
 
515
 
 
516
    ! Local variables for copying
 
517
    integer :: i1D, i2D, j2D
 
518
 
 
519
    i2D = from_2D(1)
 
520
    j2D = from_2D(2)
 
521
    do i1D = from_1D, to_1D
 
522
       out2D(i2D, j2D) = out2D(i2D, j2D) + a1D(i1D)
 
523
       i2D = i2D + 1
 
524
       if ( i2D > to_2D(1) ) then
 
525
          i2D = from_2D(1)
 
526
          j2D = j2D + 1
 
527
       end if
 
528
    end do
 
529
 
 
530
    if ( i2D /= from_2D(1) ) &
 
531
         call die("double: 1D+>2D failed (i)")
 
532
    if ( j2D <= to_2D(2) ) &
 
533
         call die("double: 1D+>2D failed (j)")
 
534
 
 
535
  end subroutine aa_1d_2d_dp
 
536
 
 
537
  ! Adds a 1D array to a 3D array
 
538
  subroutine aa_1d_3d_ip(from_1D, to_1D, a1D, from_3D, to_3D, out3D)
 
539
    integer, intent(in) :: from_1D, to_1D
 
540
    integer, intent(in) :: a1D(:)
 
541
    integer, intent(in) :: from_3D(3), to_3D(3)
 
542
    integer, intent(inout) :: out3D(:,:,:)
 
543
 
 
544
    ! Local variables for copying
 
545
    integer :: i1D, i3D, j3D, k3D
 
546
 
 
547
    i3D = from_3D(1)
 
548
    j3D = from_3D(2)
 
549
    k3D = from_3D(3)
 
550
    do i1D = from_1D, to_1D
 
551
       out3D(i3D, j3D, k3D) = out3D(i3D, j3D, k3D) + a1D(i1D)
 
552
       i3D = i3D + 1
 
553
       if ( i3D > to_3D(1) ) then
 
554
          i3D = from_3D(1)
 
555
          j3D = j3D + 1
 
556
       end if
 
557
       if ( j3D > to_3D(2) ) then
 
558
          j3D = from_3D(2)
 
559
          k3D = k3D + 1
 
560
       end if
 
561
    end do
 
562
 
 
563
    if ( i3D /= from_3D(1) ) &
 
564
         call die("integer: 1D+>3D failed (i)")
 
565
    if ( j3D /= from_3D(2) ) &
 
566
         call die("integer: 1D+>3D failed (j)")
 
567
    if ( k3D <= to_3D(3) ) &
 
568
         call die("integer: 1D+>3D failed (k)")
 
569
 
 
570
  end subroutine aa_1d_3d_ip
 
571
  subroutine aa_1d_3d_sp(from_1D, to_1D, a1D, from_3D, to_3D, out3D)
 
572
    integer, intent(in) :: from_1D, to_1D
 
573
    real(sp), intent(in) :: a1D(:)
 
574
    integer, intent(in) :: from_3D(3), to_3D(3)
 
575
    real(sp), intent(inout) :: out3D(:,:,:)
 
576
 
 
577
    ! Local variables for copying
 
578
    integer :: i1D, i3D, j3D, k3D
 
579
 
 
580
    i3D = from_3D(1)
 
581
    j3D = from_3D(2)
 
582
    k3D = from_3D(3)
 
583
    do i1D = from_1D, to_1D
 
584
       out3D(i3D, j3D, k3D) = out3D(i3D, j3D, k3D) + a1D(i1D)
 
585
       i3D = i3D + 1
 
586
       if ( i3D > to_3D(1) ) then
 
587
          i3D = from_3D(1)
 
588
          j3D = j3D + 1
 
589
       end if
 
590
       if ( j3D > to_3D(2) ) then
 
591
          j3D = from_3D(2)
 
592
          k3D = k3D + 1
 
593
       end if
 
594
    end do
 
595
 
 
596
    if ( i3D /= from_3D(1) ) &
 
597
         call die("real: 1D+>3D failed (i)")
 
598
    if ( j3D /= from_3D(2) ) &
 
599
         call die("real: 1D+>3D failed (j)")
 
600
    if ( k3D <= to_3D(3) ) &
 
601
         call die("real: 1D+>3D failed (k)")
 
602
 
 
603
  end subroutine aa_1d_3d_sp
 
604
  subroutine aa_1d_3d_dp(from_1D, to_1D, a1D, from_3D, to_3D, out3D)
 
605
    integer, intent(in) :: from_1D, to_1D
 
606
    real(dp), intent(in) :: a1D(:)
 
607
    integer, intent(in) :: from_3D(3), to_3D(3)
 
608
    real(dp), intent(inout) :: out3D(:,:,:)
 
609
 
 
610
    ! Local variables for copying
 
611
    integer :: i1D, i3D, j3D, k3D
 
612
 
 
613
    i3D = from_3D(1)
 
614
    j3D = from_3D(2)
 
615
    k3D = from_3D(3)
 
616
    do i1D = from_1D, to_1D
 
617
       out3D(i3D, j3D, k3D) = out3D(i3D, j3D, k3D) + a1D(i1D)
 
618
       i3D = i3D + 1
 
619
       if ( i3D > to_3D(1) ) then
 
620
          i3D = from_3D(1)
 
621
          j3D = j3D + 1
 
622
       end if
 
623
       if ( j3D > to_3D(2) ) then
 
624
          j3D = from_3D(2)
 
625
          k3D = k3D + 1
 
626
       end if
 
627
    end do
 
628
 
 
629
    if ( i3D /= from_3D(1) ) &
 
630
         call die("double: 1D+>3D failed (i)")
 
631
    if ( j3D /= from_3D(2) ) &
 
632
         call die("double: 1D+>3D failed (j)")
 
633
    if ( k3D <= to_3D(3) ) &
 
634
         call die("double: 1D+>3D failed (k)")
 
635
 
 
636
  end subroutine aa_1d_3d_dp
 
637
 
 
638
  ! Adds a 1D array to a 4D array
 
639
  subroutine aa_1d_4d_ip(from_1D, to_1D, a1D, from_4D, to_4D, out4D)
 
640
    integer, intent(in) :: from_1D, to_1D
 
641
    integer, intent(in) :: a1D(:)
 
642
    integer, intent(in) :: from_4D(4), to_4D(4)
 
643
    integer, intent(inout) :: out4D(:,:,:,:)
 
644
 
 
645
    ! Local variables for copying
 
646
    integer :: i1D, i4D, j4D, k4D, m4D
 
647
 
 
648
    i4D = from_4D(1)
 
649
    j4D = from_4D(2)
 
650
    k4D = from_4D(3)
 
651
    m4D = from_4D(4)
 
652
    do i1D = from_1D, to_1D
 
653
       out4D(i4D, j4D, k4D, m4D) = out4D(i4D, j4D, k4D, m4D) + a1D(i1D)
 
654
       i4D = i4D + 1
 
655
       if ( i4D > to_4D(1) ) then
 
656
          i4D = from_4D(1)
 
657
          j4D = j4D + 1
 
658
       end if
 
659
       if ( j4D > to_4D(2) ) then
 
660
          j4D = from_4D(2)
 
661
          k4D = k4D + 1
 
662
       end if
 
663
       if ( k4D > to_4D(3) ) then
 
664
          k4D = from_4D(3)
 
665
          m4D = m4D + 1
 
666
       end if
 
667
    end do
 
668
 
 
669
    if ( i4D /= from_4D(1) ) &
 
670
         call die("integer: 1D+>4D failed (i)")
 
671
    if ( j4D /= from_4D(2) ) &
 
672
         call die("integer: 1D+>4D failed (j)")
 
673
    if ( k4D /= from_4D(3) ) &
 
674
         call die("integer: 1D+>4D failed (k)")
 
675
    if ( m4D <= to_4D(4) ) &
 
676
         call die("integer: 1D+>4D failed (m)")
 
677
 
 
678
  end subroutine aa_1d_4d_ip
 
679
  subroutine aa_1d_4d_sp(from_1D, to_1D, a1D, from_4D, to_4D, out4D)
 
680
    integer, intent(in) :: from_1D, to_1D
 
681
    real(sp), intent(in) :: a1D(:)
 
682
    integer, intent(in) :: from_4D(4), to_4D(4)
 
683
    real(sp), intent(inout) :: out4D(:,:,:,:)
 
684
 
 
685
    ! Local variables for copying
 
686
    integer :: i1D, i4D, j4D, k4D, m4D
 
687
 
 
688
    i4D = from_4D(1)
 
689
    j4D = from_4D(2)
 
690
    k4D = from_4D(3)
 
691
    m4D = from_4D(4)
 
692
    do i1D = from_1D, to_1D
 
693
       out4D(i4D, j4D, k4D, m4D) = out4D(i4D, j4D, k4D, m4D) + a1D(i1D)
 
694
       i4D = i4D + 1
 
695
       if ( i4D > to_4D(1) ) then
 
696
          i4D = from_4D(1)
 
697
          j4D = j4D + 1
 
698
       end if
 
699
       if ( j4D > to_4D(2) ) then
 
700
          j4D = from_4D(2)
 
701
          k4D = k4D + 1
 
702
       end if
 
703
       if ( k4D > to_4D(3) ) then
 
704
          k4D = from_4D(3)
 
705
          m4D = m4D + 1
 
706
       end if
 
707
    end do
 
708
 
 
709
    if ( i4D /= from_4D(1) ) &
 
710
         call die("real: 1D+>4D failed (i)")
 
711
    if ( j4D /= from_4D(2) ) &
 
712
         call die("real: 1D+>4D failed (j)")
 
713
    if ( k4D /= from_4D(3) ) &
 
714
         call die("real: 1D+>4D failed (k)")
 
715
    if ( m4D <= to_4D(4) ) &
 
716
         call die("real: 1D+>4D failed (m)")
 
717
 
 
718
  end subroutine aa_1d_4d_sp
 
719
  subroutine aa_1d_4d_dp(from_1D, to_1D, a1D, from_4D, to_4D, out4D)
 
720
    integer, intent(in) :: from_1D, to_1D
 
721
    real(dp), intent(in) :: a1D(:)
 
722
    integer, intent(in) :: from_4D(4), to_4D(4)
 
723
    real(dp), intent(inout) :: out4D(:,:,:,:)
 
724
 
 
725
    ! Local variables for copying
 
726
    integer :: i1D, i4D, j4D, k4D, m4D
 
727
 
 
728
    i4D = from_4D(1)
 
729
    j4D = from_4D(2)
 
730
    k4D = from_4D(3)
 
731
    m4D = from_4D(4)
 
732
    do i1D = from_1D, to_1D
 
733
       out4D(i4D, j4D, k4D, m4D) = out4D(i4D, j4D, k4D, m4D) + a1D(i1D)
 
734
       i4D = i4D + 1
 
735
       if ( i4D > to_4D(1) ) then
 
736
          i4D = from_4D(1)
 
737
          j4D = j4D + 1
 
738
       end if
 
739
       if ( j4D > to_4D(2) ) then
 
740
          j4D = from_4D(2)
 
741
          k4D = k4D + 1
 
742
       end if
 
743
       if ( k4D > to_4D(3) ) then
 
744
          k4D = from_4D(3)
 
745
          m4D = m4D + 1
 
746
       end if
 
747
    end do
 
748
 
 
749
    if ( i4D /= from_4D(1) ) &
 
750
         call die("double: 1D+>4D failed (i)")
 
751
    if ( j4D /= from_4D(2) ) &
 
752
         call die("double: 1D+>4D failed (j)")
 
753
    if ( k4D /= from_4D(3) ) &
 
754
         call die("double: 1D+>4D failed (k)")
 
755
    if ( m4D <= to_4D(4) ) &
 
756
         call die("double: 1D+>4D failed (m)")
 
757
 
 
758
  end subroutine aa_1d_4d_dp
 
759
  
 
760
  ! Adds a 2D array to a 1D array
 
761
  subroutine aa_2d_1d_ip(from_2D, to_2D, in2D, from_1D, to_1D, out1D)
 
762
    integer, intent(in) :: from_2D(2), to_2D(2)
 
763
    integer, intent(in) :: in2D(:,:)
 
764
    integer, intent(in) :: from_1D, to_1D
 
765
    integer, intent(inout) :: out1D(:)
 
766
 
 
767
    ! Local variables for copying
 
768
    integer :: i2D, j2D, i1D
 
769
 
 
770
    i1D = from_1D
 
771
    do j2D = from_2D(2), to_2D(2)
 
772
       do i2D = from_2D(1), to_2D(1)
 
773
          out1D(i1D) = out1D(i1D) + in2D(i2D, j2D)
 
774
          i1D = i1D + 1
 
775
       end do
 
776
    end do
 
777
 
 
778
    if ( i1D <= to_1D ) &
 
779
         call die("integer: 2D+>1D failed")
 
780
 
 
781
  end subroutine aa_2d_1d_ip
 
782
  subroutine aa_2d_1d_sp(from_2D, to_2D, in2D, from_1D, to_1D, out1D)
 
783
    integer, intent(in) :: from_2D(2), to_2D(2)
 
784
    real(sp), intent(in) :: in2D(:,:)
 
785
    integer, intent(in) :: from_1D, to_1D
 
786
    real(sp), intent(inout) :: out1D(:)
 
787
 
 
788
    ! Local variables for copying
 
789
    integer :: i2D, j2D, i1D
 
790
 
 
791
    i1D = from_1D
 
792
    do j2D = from_2D(2), to_2D(2)
 
793
       do i2D = from_2D(1), to_2D(1)
 
794
          out1D(i1D) = out1D(i1D) + in2D(i2D, j2D)
 
795
          i1D = i1D + 1
 
796
       end do
 
797
    end do
 
798
 
 
799
    if ( i1D <= to_1D ) &
 
800
         call die("real: 2D+>1D failed")
 
801
 
 
802
  end subroutine aa_2d_1d_sp
 
803
  subroutine aa_2d_1d_dp(from_2D, to_2D, in2D, from_1D, to_1D, out1D)
 
804
    integer, intent(in) :: from_2D(2), to_2D(2)
 
805
    real(dp), intent(in) :: in2D(:,:)
 
806
    integer, intent(in) :: from_1D, to_1D
 
807
    real(dp), intent(inout) :: out1D(:)
 
808
 
 
809
    ! Local variables for copying
 
810
    integer :: i2D, j2D, i1D
 
811
 
 
812
    i1D = from_1D
 
813
    do j2D = from_2D(2), to_2D(2)
 
814
       do i2D = from_2D(1), to_2D(1)
 
815
          out1D(i1D) = out1D(i1D) + in2D(i2D, j2D)
 
816
          i1D = i1D + 1
 
817
       end do
 
818
    end do
 
819
 
 
820
    if ( i1D <= to_1D ) &
 
821
         call die("double: 2D+>1D failed")
 
822
 
 
823
  end subroutine aa_2d_1d_dp
 
824
  
 
825
  ! Adds a 3D array to a 1D array
 
826
  subroutine aa_3d_1d_ip(from_3D, to_3D, in3D, from_1D, to_1D, out1D)
 
827
    integer, intent(in) :: from_3D(3), to_3D(3)
 
828
    integer, intent(in) :: in3D(:,:,:)
 
829
    integer, intent(in) :: from_1D, to_1D
 
830
    integer, intent(inout) :: out1D(:)
 
831
 
 
832
    ! Local variables for copying
 
833
    integer :: i3D, j3D, k3D, i1D
 
834
 
 
835
    i1D = from_1D
 
836
    do k3D = from_3D(3), to_3D(3)
 
837
       do j3D = from_3D(2), to_3D(2)
 
838
          do i3D = from_3D(1), to_3D(1)
 
839
             out1D(i1D) = out1D(i1D) + in3D(i3D, j3D, k3D)
 
840
             i1D = i1D + 1
 
841
          end do
 
842
       end do
 
843
    end do
 
844
 
 
845
    if ( i1D <= to_1D ) &
 
846
         call die("integer: 3D+>1D failed")
 
847
 
 
848
  end subroutine aa_3d_1d_ip
 
849
  subroutine aa_3d_1d_sp(from_3D, to_3D, in3D, from_1D, to_1D, out1D)
 
850
    integer, intent(in) :: from_3D(3), to_3D(3)
 
851
    real(sp), intent(in) :: in3D(:,:,:)
 
852
    integer, intent(in) :: from_1D, to_1D
 
853
    real(sp), intent(inout) :: out1D(:)
 
854
 
 
855
    ! Local variables for copying
 
856
    integer :: i3D, j3D, k3D, i1D
 
857
 
 
858
    i1D = from_1D
 
859
    do k3D = from_3D(3), to_3D(3)
 
860
       do j3D = from_3D(2), to_3D(2)
 
861
          do i3D = from_3D(1), to_3D(1)
 
862
             out1D(i1D) = out1D(i1D) + in3D(i3D, j3D, k3D)
 
863
             i1D = i1D + 1
 
864
          end do
 
865
       end do
 
866
    end do
 
867
 
 
868
    if ( i1D <= to_1D ) &
 
869
         call die("real: 3D+>1D failed")
 
870
 
 
871
  end subroutine aa_3d_1d_sp
 
872
  subroutine aa_3d_1d_dp(from_3D, to_3D, in3D, from_1D, to_1D, out1D)
 
873
    integer, intent(in) :: from_3D(3), to_3D(3)
 
874
    real(dp), intent(in) :: in3D(:,:,:)
 
875
    integer, intent(in) :: from_1D, to_1D
 
876
    real(dp), intent(inout) :: out1D(:)
 
877
 
 
878
    ! Local variables for copying
 
879
    integer :: i3D, j3D, k3D, i1D
 
880
 
 
881
    i1D = from_1D
 
882
    do k3D = from_3D(3), to_3D(3)
 
883
       do j3D = from_3D(2), to_3D(2)
 
884
          do i3D = from_3D(1), to_3D(1)
 
885
             out1D(i1D) = out1D(i1D) + in3D(i3D, j3D, k3D)
 
886
             i1D = i1D + 1
 
887
          end do
 
888
       end do
 
889
    end do
 
890
 
 
891
    if ( i1D <= to_1D ) &
 
892
         call die("double: 3D+>1D failed")
 
893
 
 
894
  end subroutine aa_3d_1d_dp
 
895
 
 
896
  ! Adds a 4D array to a 1D array
 
897
  subroutine aa_4d_1d_ip(from_4D, to_4D, in4D, from_1D, to_1D, out1D)
 
898
    integer, intent(in) :: from_4D(4), to_4D(4)
 
899
    integer, intent(in) :: in4D(:,:,:,:)
 
900
    integer, intent(in) :: from_1D, to_1D
 
901
    integer, intent(inout) :: out1D(:)
 
902
 
 
903
    ! Local variables for copying
 
904
    integer :: i4D, j4D, k4D, m4D, i1D
 
905
 
 
906
    i1D = from_1D
 
907
    do m4D = from_4D(4), to_4D(4)
 
908
       do k4D = from_4D(3), to_4D(3)
 
909
          do j4D = from_4D(2), to_4D(2)
 
910
             do i4D = from_4D(1), to_4D(1)
 
911
                out1D(i1D) = out1D(i1D) + in4D(i4D, j4D, k4D, m4D)
 
912
                i1D = i1D + 1
 
913
             end do
 
914
          end do
 
915
       end do
 
916
    end do
 
917
    
 
918
    if ( i1D <= to_1D ) &
 
919
         call die("integer: 4D+>1D failed")
 
920
 
 
921
  end subroutine aa_4d_1d_ip
 
922
  subroutine aa_4d_1d_sp(from_4D, to_4D, in4D, from_1D, to_1D, out1D)
 
923
    integer, intent(in) :: from_4D(4), to_4D(4)
 
924
    real(sp), intent(in) :: in4D(:,:,:,:)
 
925
    integer, intent(in) :: from_1D, to_1D
 
926
    real(sp), intent(inout) :: out1D(:)
 
927
 
 
928
    ! Local variables for copying
 
929
    integer :: i4D, j4D, k4D, m4D, i1D
 
930
 
 
931
    i1D = from_1D
 
932
    do m4D = from_4D(4), to_4D(4)
 
933
       do k4D = from_4D(3), to_4D(3)
 
934
          do j4D = from_4D(2), to_4D(2)
 
935
             do i4D = from_4D(1), to_4D(1)
 
936
                out1D(i1D) = out1D(i1D) + in4D(i4D, j4D, k4D, m4D)
 
937
                i1D = i1D + 1
 
938
             end do
 
939
          end do
 
940
       end do
 
941
    end do
 
942
    
 
943
    if ( i1D <= to_1D ) &
 
944
         call die("real: 4D+>1D failed")
 
945
 
 
946
  end subroutine aa_4d_1d_sp
 
947
  subroutine aa_4d_1d_dp(from_4D, to_4D, in4D, from_1D, to_1D, out1D)
 
948
    integer, intent(in) :: from_4D(4), to_4D(4)
 
949
    real(dp), intent(in) :: in4D(:,:,:,:)
 
950
    integer, intent(in) :: from_1D, to_1D
 
951
    real(dp), intent(inout) :: out1D(:)
 
952
 
 
953
    ! Local variables for copying
 
954
    integer :: i4D, j4D, k4D, m4D, i1D
 
955
 
 
956
    i1D = from_1D
 
957
    do m4D = from_4D(4), to_4D(4)
 
958
       do k4D = from_4D(3), to_4D(3)
 
959
          do j4D = from_4D(2), to_4D(2)
 
960
             do i4D = from_4D(1), to_4D(1)
 
961
                out1D(i1D) = out1D(i1D) + in4D(i4D, j4D, k4D, m4D)
 
962
                i1D = i1D + 1
 
963
             end do
 
964
          end do
 
965
       end do
 
966
    end do
 
967
 
 
968
    if ( i1D <= to_1D ) &
 
969
         call die("double: 4D+>1D failed")
 
970
 
 
971
  end subroutine aa_4d_1d_dp
 
972
 
 
973
end module m_array
 
974