diff options
Diffstat (limited to 'Carpet/CarpetLib/src/restrict_3d_real8.F77')
| -rw-r--r-- | Carpet/CarpetLib/src/restrict_3d_real8.F77 | 128 |
1 files changed, 128 insertions, 0 deletions
diff --git a/Carpet/CarpetLib/src/restrict_3d_real8.F77 b/Carpet/CarpetLib/src/restrict_3d_real8.F77 new file mode 100644 index 000000000..81f4cfd0a --- /dev/null +++ b/Carpet/CarpetLib/src/restrict_3d_real8.F77 @@ -0,0 +1,128 @@ +c -*-Fortran-*- +c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/restrict_3d_real8.F77,v 1.7 2004/03/11 12:03:09 schnetter Exp $ + +#include "cctk.h" + + + +#define CHKIDX(i,j,k, imax,jmax,kmax, where) \ + if ((i).lt.1 .or. (i).gt.(imax) \ + .or. (j).lt.1 .or. (j).gt.(jmax) \ + .or. (k).lt.1 .or. (k).gt.(kmax)) then &&\ + write (msg, '(a, " array index out of bounds: shape is (",i4,",",i4,",",i4,"), index is (",i4,",",i4,",",i4,")")') \ + (where), (imax), (jmax), (kmax), (i), (j), (k) &&\ + call CCTK_WARN (0, msg(1:len_trim(msg))) &&\ + end if + + + + subroutine restrict_3d_real8 ( + $ src, srciext, srcjext, srckext, + $ dst, dstiext, dstjext, dstkext, + $ srcbbox, dstbbox, regbbox) + + implicit none + + 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 srcifac, srcjfac, srckfac + + integer srcioff, srcjoff, srckoff + integer dstioff, dstjoff, dstkoff + + integer i, j, k + integer d + + character msg*1000 + + + + 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 (mod(dstbbox(d,3), srcbbox(d,3)).ne.0) then + call CCTK_WARN (0, "Internal error: source strides are not integer multiples of the destination 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 + + srcifac = dstbbox(1,3) / srcbbox(1,3) + srcjfac = dstbbox(2,3) / srcbbox(2,3) + srckfac = dstbbox(3,3) / srcbbox(3,3) + + 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 + + CHKIDX (srcioff+srcifac*i+1, srcjoff+srcjfac*j+1, srckoff+srckfac*k+1, \ + srciext,srcjext,srckext, "source") + CHKIDX (dstioff+i+1, dstjoff+j+1, dstkoff+k+1, \ + dstiext,dstjext,dstkext, "destination") + dst (dstioff+i+1, dstjoff+j+1, dstkoff+k+1) + $ = src (srcioff+srcifac*i+1, srcjoff+srcjfac*j+1, srckoff+srckfac*k+1) + + end do + end do + end do + + end |