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
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
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
107
108
109
110
111
|
c -*-Fortran-*-
c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/restrict_3d_real8_rf2.F77,v 1.2 2004/03/11 12:03:09 schnetter Exp $
#include "cctk.h"
#include "cctk_Parameters.h"
subroutine restrict_3d_real8_rf2 (
$ src, srciext, srcjext, srckext,
$ dst, dstiext, dstjext, dstkext,
$ srcbbox, dstbbox, regbbox)
implicit none
DECLARE_CCTK_PARAMETERS
integer srciext, srcjext, srckext
CCTK_REAL8 src(srciext,srcjext,srckext)
integer dstiext, dstjext, dstkext
CCTK_REAL8 dst(dstiext,dstjext,dstkext)
c bbox(:,1) is lower boundary (inclusive)
c bbox(:,2) is upper boundary (inclusive)
c bbox(:,3) is stride
integer srcbbox(3,3), dstbbox(3,3), regbbox(3,3)
integer regiext, regjext, regkext
integer srcioff, srcjoff, srckoff
integer dstioff, dstjoff, dstkoff
integer i, j, k
integer d
do d=1,3
if (srcbbox(d,3).eq.0 .or. dstbbox(d,3).eq.0
$ .or. regbbox(d,3).eq.0) then
call CCTK_WARN (0, "Internal error: stride is zero")
end if
if (srcbbox(d,3).ge.regbbox(d,3)
$ .or. dstbbox(d,3).ne.regbbox(d,3)) then
call CCTK_WARN (0, "Internal error: strides disagree")
end if
if (dstbbox(d,3).ne.srcbbox(d,3)*2) then
call CCTK_WARN (0, "Internal error: destination strides are not twice the source strides")
end if
if (mod(srcbbox(d,1), srcbbox(d,3)).ne.0
$ .or. mod(dstbbox(d,1), dstbbox(d,3)).ne.0
$ .or. mod(regbbox(d,1), regbbox(d,3)).ne.0) then
call CCTK_WARN (0, "Internal error: array origins are not integer multiples of the strides")
end if
if (regbbox(d,1).gt.regbbox(d,2)) then
c This could be handled, but is likely to point to an error elsewhere
call CCTK_WARN (0, "Internal error: region extent is empty")
end if
if (mod(srcbbox(d,2) - srcbbox(d,1), srcbbox(d,3)).ne.0
$ .or. mod(dstbbox(d,2) - dstbbox(d,1), dstbbox(d,3)).ne.0
$ .or. mod(regbbox(d,2) - regbbox(d,1), regbbox(d,3)).ne.0) then
call CCTK_WARN (0, "Internal error: array extents are not integer multiples of the strides")
end if
if (regbbox(d,1).lt.srcbbox(d,1)
$ .or. regbbox(d,1).lt.dstbbox(d,1)
$ .or. regbbox(d,2).gt.srcbbox(d,2)
$ .or. regbbox(d,2).gt.dstbbox(d,2)) then
call CCTK_WARN (0, "Internal error: region extent is not contained in array extent")
end if
end do
if (srciext.ne.(srcbbox(1,2)-srcbbox(1,1))/srcbbox(1,3)+1
$ .or. srcjext.ne.(srcbbox(2,2)-srcbbox(2,1))/srcbbox(2,3)+1
$ .or. srckext.ne.(srcbbox(3,2)-srcbbox(3,1))/srcbbox(3,3)+1
$ .or. dstiext.ne.(dstbbox(1,2)-dstbbox(1,1))/dstbbox(1,3)+1
$ .or. dstjext.ne.(dstbbox(2,2)-dstbbox(2,1))/dstbbox(2,3)+1
$ .or. dstkext.ne.(dstbbox(3,2)-dstbbox(3,1))/dstbbox(3,3)+1) then
call CCTK_WARN (0, "Internal error: array sizes don't agree with bounding boxes")
end if
regiext = (regbbox(1,2) - regbbox(1,1)) / regbbox(1,3) + 1
regjext = (regbbox(2,2) - regbbox(2,1)) / regbbox(2,3) + 1
regkext = (regbbox(3,2) - regbbox(3,1)) / regbbox(3,3) + 1
srcioff = (regbbox(1,1) - srcbbox(1,1)) / srcbbox(1,3)
srcjoff = (regbbox(2,1) - srcbbox(2,1)) / srcbbox(2,3)
srckoff = (regbbox(3,1) - srcbbox(3,1)) / srcbbox(3,3)
dstioff = (regbbox(1,1) - dstbbox(1,1)) / dstbbox(1,3)
dstjoff = (regbbox(2,1) - dstbbox(2,1)) / dstbbox(2,3)
dstkoff = (regbbox(3,1) - dstbbox(3,1)) / dstbbox(3,3)
c Loop over coarse region
do k = 0, regkext-1
do j = 0, regjext-1
do i = 0, regiext-1
if (check_array_accesses.ne.0) then
call checkindex (srcioff+2*i+1, srcjoff+2*j+1, srckoff+2*k+1, 1,1,1, srciext,srcjext,srckext, "source")
call checkindex (dstioff+i+1, dstjoff+j+1, dstkoff+k+1, 1,1,1, dstiext,dstjext,dstkext, "destination")
end if
dst(dstioff+i+1, dstjoff+j+1, dstkoff+k+1) =
$ src(srcioff+2*i+1, srcjoff+2*j+1, srckoff+2*k+1)
end do
end do
end do
end
|
