diff options
Diffstat (limited to 'Carpet')
-rw-r--r-- | Carpet/CarpetLib/src/bbox.cc | 7 | ||||
-rw-r--r-- | Carpet/CarpetLib/src/bboxset.cc | 6 | ||||
-rw-r--r-- | Carpet/CarpetLib/src/copy_3d_real8.F77 | 29 | ||||
-rw-r--r-- | Carpet/CarpetLib/src/dh.cc | 271 | ||||
-rw-r--r-- | Carpet/CarpetLib/src/gdata.cc | 7 | ||||
-rw-r--r-- | Carpet/CarpetLib/src/ggf.cc | 4 | ||||
-rw-r--r-- | Carpet/CarpetLib/src/prolongate_3d_real8.F77 | 105 | ||||
-rw-r--r-- | Carpet/CarpetLib/src/prolongate_3d_real8_2tl.F77 | 29 | ||||
-rw-r--r-- | Carpet/CarpetLib/src/prolongate_3d_real8_2tl_o3.F77 | 43 | ||||
-rw-r--r-- | Carpet/CarpetLib/src/prolongate_3d_real8_3tl.F77 | 29 | ||||
-rw-r--r-- | Carpet/CarpetLib/src/prolongate_3d_real8_3tl_o3.F77 | 43 | ||||
-rw-r--r-- | Carpet/CarpetLib/src/prolongate_3d_real8_o3.F77 | 43 | ||||
-rw-r--r-- | Carpet/CarpetLib/src/restrict_3d_real8.F77 | 29 |
13 files changed, 484 insertions, 161 deletions
diff --git a/Carpet/CarpetLib/src/bbox.cc b/Carpet/CarpetLib/src/bbox.cc index 0369aec9e..dcae93c02 100644 --- a/Carpet/CarpetLib/src/bbox.cc +++ b/Carpet/CarpetLib/src/bbox.cc @@ -1,4 +1,4 @@ -// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/bbox.cc,v 1.12 2003/02/24 17:11:29 schnetter Exp $ +// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/bbox.cc,v 1.13 2003/02/25 22:57:00 schnetter Exp $ #include <assert.h> @@ -117,6 +117,7 @@ bbox<T,D> bbox<T,D>::operator& (const bbox& b) const { // Containment template<class T, int D> bool bbox<T,D>::is_contained_in (const bbox& b) const { + if (empty()) return true; // no alignment check return all(lower()>=b.lower() && upper()<=b.upper()); } @@ -140,7 +141,7 @@ bbox<T,D> bbox<T,D>::expand (const vect<T,D>& lo, const vect<T,D>& hi) const { // Find the smallest b-compatible box around *this template<class T, int D> bbox<T,D> bbox<T,D>::expanded_for (const bbox& b) const { - assert (! empty()); + if (empty()) return bbox(b.lower(), b.lower()-b.stride(), b.stride()); const vect<T,D> str = b.stride(); const vect<T,D> loff = ((lower() - b.lower()) % str + str) % str; const vect<T,D> uoff = ((upper() - b.lower()) % str + str) % str; @@ -152,7 +153,7 @@ bbox<T,D> bbox<T,D>::expanded_for (const bbox& b) const { // Find the largest b-compatible box inside *this template<class T, int D> bbox<T,D> bbox<T,D>::contracted_for (const bbox& b) const { - assert (! empty()); + if (empty()) return bbox(b.lower(), b.lower()-b.stride(), b.stride()); const vect<T,D> str = b.stride(); const vect<T,D> loff = ((lower() - b.lower()) % str + str) % str; const vect<T,D> uoff = ((upper() - b.lower()) % str + str) % str; diff --git a/Carpet/CarpetLib/src/bboxset.cc b/Carpet/CarpetLib/src/bboxset.cc index 8de16a015..b440bb59d 100644 --- a/Carpet/CarpetLib/src/bboxset.cc +++ b/Carpet/CarpetLib/src/bboxset.cc @@ -1,4 +1,4 @@ -// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/bboxset.cc,v 1.11 2003/01/03 15:49:36 schnetter Exp $ +// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/bboxset.cc,v 1.12 2003/02/25 22:57:00 schnetter Exp $ #include <assert.h> @@ -318,7 +318,9 @@ bool bboxset<T,D>::operator!= (const bboxset<T,D>& s) const { // Output template<class T,int D> void bboxset<T,D>::output (ostream& os) const { - os << "bboxset<T," << D << ">:size=" << size() << "," << "set=" << bs; + T Tdummy; + os << "bboxset<" << typestring(Tdummy) << "," << D << ">:" + << "size=" << size() << "," << "set=" << bs; } diff --git a/Carpet/CarpetLib/src/copy_3d_real8.F77 b/Carpet/CarpetLib/src/copy_3d_real8.F77 index 230814c04..d115418b8 100644 --- a/Carpet/CarpetLib/src/copy_3d_real8.F77 +++ b/Carpet/CarpetLib/src/copy_3d_real8.F77 @@ -1,7 +1,20 @@ c -*-Fortran-*- -c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/copy_3d_real8.F77,v 1.5 2003/02/24 17:43:10 schnetter Exp $ +c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/copy_3d_real8.F77,v 1.6 2003/02/25 22:57:00 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 copy_3d_real8 ( $ src, srciext, srcjext, srckext, @@ -27,6 +40,8 @@ c bbox(:,3) is stride integer i, j, k integer d + character msg*1000 + do d=1,3 @@ -43,16 +58,16 @@ c bbox(:,3) is stride $ .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 (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 - 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 end do if (srciext.ne.(srcbbox(1,2)-srcbbox(1,1))/srcbbox(1,3)+1 @@ -85,6 +100,10 @@ c Loop over region do j = 0, regjext-1 do i = 0, regiext-1 + CHKIDX (srcioff+i+1, srcjoff+j+1, srckoff+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+i+1, srcjoff+j+1, srckoff+k+1) diff --git a/Carpet/CarpetLib/src/dh.cc b/Carpet/CarpetLib/src/dh.cc index 68f4c18c4..7b856565a 100644 --- a/Carpet/CarpetLib/src/dh.cc +++ b/Carpet/CarpetLib/src/dh.cc @@ -1,4 +1,4 @@ -// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dh.cc,v 1.25 2003/01/10 23:07:16 schnetter Exp $ +// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dh.cc,v 1.26 2003/02/25 22:57:00 schnetter Exp $ #include <assert.h> @@ -76,6 +76,7 @@ void dh<D>::recompose () { if (h.outer_boundaries[rl][c][d][0]) ldist[d] = 0; if (h.outer_boundaries[rl][c][d][1]) udist[d] = 0; } + assert (! intr.empty()); boxes[rl][c][ml].exterior = intr.expand(ldist, udist); // Boundaries (ghost zones only) @@ -116,10 +117,7 @@ void dh<D>::recompose () { for (int rl=0; rl<h.reflevels(); ++rl) { for (int c=0; c<h.components(rl); ++c) { for (int ml=0; ml<h.mglevels(rl,c); ++ml) { - const ibbox intr = boxes[rl][c][ml].interior; - const ibbox extr = boxes[rl][c][ml].exterior; - - const ibset bnds = boxes[rl][c][ml].boundaries; + const ibset& bnds = boxes[rl][c][ml].boundaries; // Sync boxes for (int cc=0; cc<h.components(rl); ++cc) { @@ -133,6 +131,16 @@ void dh<D>::recompose () { } } // for cc + } // for ml + } // for c + } // for rl + + for (int rl=0; rl<h.reflevels(); ++rl) { + for (int c=0; c<h.components(rl); ++c) { + for (int ml=0; ml<h.mglevels(rl,c); ++ml) { + const ibbox& intr = boxes[rl][c][ml].interior; + const ibbox& extr = boxes[rl][c][ml].exterior; + // Multigrid boxes if (ml>0) { const ibbox intrf = boxes[rl][c][ml-1].interior; @@ -142,7 +150,9 @@ void dh<D>::recompose () { // (the restriction must fill all of the interior of the // coarse grid, and may use the exterior of the fine grid) const ibbox recv = intr; + assert (! recv.empty()); const ibbox send = recv.expanded_for(extrf); + assert (! send.empty()); assert (send.is_contained_in(extrf)); boxes[rl][c][ml-1].send_mg_coarse.push_back(send); boxes[rl][c][ml ].recv_mg_fine .push_back(recv); @@ -153,30 +163,57 @@ void dh<D>::recompose () { // grid, and may fill only the interior of the fine grid, // and the bbox must be as large as possible) const ibbox recv = extr.contracted_for(intrf) & intrf; + assert (! recv.empty()); const ibbox send = recv.expanded_for(extr); + assert (! send.empty()); boxes[rl][c][ml-1].recv_mg_coarse.push_back(recv); boxes[rl][c][ml ].send_mg_fine .push_back(send); } } // if not finest multigrid level - + + } // for ml + } // for c + } // for rl + + for (int rl=0; rl<h.reflevels(); ++rl) { + for (int c=0; c<h.components(rl); ++c) { + for (int ml=0; ml<h.mglevels(rl,c); ++ml) { + const ibbox& intr = boxes[rl][c][ml].interior; + const ibbox& extr = boxes[rl][c][ml].exterior; + // Refinement boxes if (rl<h.reflevels()-1) { for (int cc=0; cc<h.components(rl+1); ++cc) { const ibbox intrf = boxes[rl+1][cc][ml].interior; -// const ibbox extrf = boxes[rl+1][cc][ml].exterior; // Prolongation (interior) { // (the prolongation may use the exterior of the coarse // grid, and must fill all of the interior of the fine // grid) - const ibbox recv = extr.contracted_for(intrf) & intrf; - const ibbox send = recv.expanded_for(extr); - assert (send.is_contained_in(extr)); - boxes[rl+1][cc][ml].recv_ref_coarse[c ].push_back(recv); - boxes[rl ][c ][ml].send_ref_fine [cc].push_back(send); + const int pss = prolongation_stencil_size(); + ibset recvs = (extr.expand(-pss,-pss).contracted_for(intrf) + & intrf); + const iblistvect& rrc = boxes[rl+1][cc][ml].recv_ref_coarse; + for (typename iblistvect::const_iterator lvi=rrc.begin(); + lvi!=rrc.end(); ++lvi) { + for (typename iblist::const_iterator li=lvi->begin(); + li!=lvi->end(); ++li) { + recvs -= *li; + } + } + assert (recvs.setsize() <= 1); + if (recvs.setsize() == 1) { + const ibbox recv = *recvs.begin(); + const ibbox send = recv.expanded_for(extr); + assert (! send.empty()); + assert (send.is_contained_in(extr)); + boxes[rl+1][cc][ml].recv_ref_coarse[c ].push_back(recv); + boxes[rl ][c ][ml].send_ref_fine [cc].push_back(send); + } } // Prolongation (boundaries) { + const int pss = prolongation_stencil_size(); ibset bndsf = boxes[rl+1][cc][ml].boundaries; // coarsify boundaries of fine component for (typename ibset::const_iterator bi=bndsf.begin(); @@ -185,13 +222,34 @@ void dh<D>::recompose () { // (the prolongation may use the exterior of the // coarse grid, and must fill all of the boundary of // the fine grid) - const int pss = prolongation_stencil_size(); - const ibbox recv = (extr.expand(-pss,-pss).contracted_for(bndf) - & bndf); - const ibbox send = recv.expanded_for(extr); - assert (send.is_contained_in(extr)); - boxes[rl+1][cc][ml].recv_ref_bnd_coarse[c ].push_back(recv); - boxes[rl ][c ][ml].send_ref_bnd_fine [cc].push_back(send); + ibset recvs = (extr.expand(-pss,-pss).contracted_for(bndf) + & bndf); + const iblistvect& rrbc + = boxes[rl+1][cc][ml].recv_ref_bnd_coarse; + for (typename iblistvect::const_iterator lvi=rrbc.begin(); + lvi!=rrbc.end(); ++lvi) { + for (typename iblist::const_iterator li=lvi->begin(); + li!=lvi->end(); ++li) { + recvs -= *li; + } + } + const iblistvect& rs = boxes[rl+1][cc][ml].recv_sync; + for (typename iblistvect::const_iterator lvi=rs.begin(); + lvi!=rs.end(); ++lvi) { + for (typename iblist::const_iterator li=lvi->begin(); + li!=lvi->end(); ++li) { + recvs -= *li; + } + } + assert (recvs.setsize() <= 1); + if (recvs.setsize() == 1) { + const ibbox recv = *recvs.begin(); + const ibbox send = recv.expanded_for(extr); + assert (! send.empty()); + assert (send.is_contained_in(extr)); + boxes[rl+1][cc][ml].recv_ref_bnd_coarse[c ].push_back(recv); + boxes[rl ][c ][ml].send_ref_bnd_fine [cc].push_back(send); + } } } // Restriction (interior) @@ -201,23 +259,25 @@ void dh<D>::recompose () { // grid, and the bbox must be as large as possible) // (the restriction must not fill points that are used // to prolongate the boundaries) - const ibbox recv = intrf.contracted_for(intr) & intr; - - const iblist& sendlist = boxes[rl][c][ml].send_ref_bnd_fine[cc]; - ibset recv2s (recv); - for (typename iblist::const_iterator sli = sendlist.begin(); - sli != sendlist.end(); - ++sli) { - recv2s -= *sli; - } - assert (recv2s.setsize() == 1); - const ibbox recv2 = *recv2s.begin(); -// const ibbox send = recv.expanded_for(intrf); - const ibbox send = recv2.expanded_for(intrf); - boxes[rl+1][cc][ml].send_ref_coarse[c ].push_back(send); - boxes[rl ][c ][ml].recv_ref_fine [cc].push_back(recv); + ibset recvs = intrf.contracted_for(intr) & intr; + const iblist& sendlist + = boxes[rl][c][ml].send_ref_bnd_fine[cc]; + for (typename iblist::const_iterator sli = sendlist.begin(); + sli != sendlist.end(); + ++sli) { + recvs -= *sli; + } + assert (recvs.setsize() <= 1); + if (recvs.setsize() == 1) { + const ibbox recv = *recvs.begin(); + assert (! recv.empty()); + const ibbox send = recv.expanded_for(intrf); + assert (! send.empty()); + boxes[rl+1][cc][ml].send_ref_coarse[c ].push_back(send); + boxes[rl ][c ][ml].recv_ref_fine [cc].push_back(recv); + } } - + } // for cc } // if not finest refinement level @@ -231,44 +291,41 @@ void dh<D>::recompose () { // Boundaries that are not synced, or are neither synced nor // prolonged to from coarser grids (outer boundaries) - ibset& sync_not = boxes[rl][c][ml].sync_not; - ibset& recv_not = boxes[rl][c][ml].recv_not; - - // The whole boundary - sync_not = boxes[rl][c][ml].boundaries; - recv_not = boxes[rl][c][ml].boundaries; - - // Subtract boxes received during synchronisation - const iblistvect& recv_sync = boxes[rl][c][ml].recv_sync; - for (typename iblistvect::const_iterator lvi=recv_sync.begin(); - lvi!=recv_sync.end(); ++lvi) { - for (typename iblist::const_iterator li=lvi->begin(); - li!=lvi->end(); ++li) { - sync_not -= *li; - recv_not -= *li; - } - } - - // Subtract all boxes received - const iblistvect& recv_ref_bnd_coarse - = boxes[rl][c][ml].recv_ref_bnd_coarse; - for (typename iblistvect::const_iterator - lvi=recv_ref_bnd_coarse.begin(); - lvi!=recv_ref_bnd_coarse.end(); ++lvi) { - for (typename iblist::const_iterator li=lvi->begin(); - li!=lvi->end(); ++li) { - recv_not -= *li; - } - } - - // Check that no boundaries are left over - if (rl==0) assert (sync_not.empty()); - assert (recv_not.empty()); + ibset& sync_not = boxes[rl][c][ml].sync_not; + ibset& recv_not = boxes[rl][c][ml].recv_not; + + // The whole boundary + sync_not = boxes[rl][c][ml].boundaries; + recv_not = boxes[rl][c][ml].boundaries; + + // Subtract boxes received during synchronisation + const iblistvect& recv_sync = boxes[rl][c][ml].recv_sync; + for (typename iblistvect::const_iterator lvi=recv_sync.begin(); + lvi!=recv_sync.end(); ++lvi) { + for (typename iblist::const_iterator li=lvi->begin(); + li!=lvi->end(); ++li) { + sync_not -= *li; + recv_not -= *li; + } + } + + // Subtract boxes received during prolongation + const iblistvect& recv_ref_bnd_coarse + = boxes[rl][c][ml].recv_ref_bnd_coarse; + for (typename iblistvect::const_iterator + lvi=recv_ref_bnd_coarse.begin(); + lvi!=recv_ref_bnd_coarse.end(); ++lvi) { + for (typename iblist::const_iterator li=lvi->begin(); + li!=lvi->end(); ++li) { + recv_not -= *li; + } + } } // for ml } // for c } // for rl + // Calculate bases bases.resize(h.reflevels()); for (int rl=0; rl<h.reflevels(); ++rl) { if (h.components(rl)==0) { @@ -334,6 +391,84 @@ void dh<D>::recompose () { } } // if output_bboxes + for (int rl=0; rl<h.reflevels(); ++rl) { + for (int c=0; c<h.components(rl); ++c) { + for (int ml=0; ml<h.mglevels(rl,c); ++ml) { + + // Assert that all boundaries are synced or received + { + const ibset& sync_not = boxes[rl][c][ml].sync_not; + const ibset& recv_not = boxes[rl][c][ml].recv_not; + + // Check that no boundaries are left over + if (rl==0) assert (sync_not.empty()); + assert (recv_not.empty()); + } + + // Assert that the interior is received exactly once during + // prolongation, and that nothing else is received + { + if (rl==0) { + const iblistvect& recv_ref_coarse + = boxes[rl][c][ml].recv_ref_coarse; + assert (recv_ref_coarse.empty()); + } else { // rl!=0 + const iblistvect& recv_ref_coarse + = boxes[rl][c][ml].recv_ref_coarse; + ibset intr = boxes[rl][c][ml].interior; + for (typename iblistvect::const_iterator + lvi=recv_ref_coarse.begin(); + lvi!=recv_ref_coarse.end(); ++lvi) { + for (typename iblist::const_iterator li=lvi->begin(); + li!=lvi->end(); ++li) { + const int old_sz = intr.size(); + const int this_sz = li->size(); + intr -= *li; + const int new_sz = intr.size(); + assert (new_sz + this_sz == old_sz); + } + } + assert (intr.empty()); + } + } + + // Assert that the boundaries are received at most once during + // prolongation and synchronisation, and that nothing else is + // received + { + const iblistvect& recv_sync = boxes[rl][c][ml].recv_sync; + const iblistvect& recv_ref_bnd_coarse + = boxes[rl][c][ml].recv_ref_bnd_coarse; + ibset bnds = boxes[rl][c][ml].boundaries; + for (typename iblistvect::const_iterator lvi=recv_sync.begin(); + lvi!=recv_sync.end(); ++lvi) { + for (typename iblist::const_iterator li=lvi->begin(); + li!=lvi->end(); ++li) { + const int old_sz = bnds.size(); + const int this_sz = li->size(); + bnds -= *li; + const int new_sz = bnds.size(); + assert (new_sz + this_sz == old_sz); + } + } + for (typename iblistvect::const_iterator + lvi=recv_ref_bnd_coarse.begin(); + lvi!=recv_ref_bnd_coarse.end(); ++lvi) { + for (typename iblist::const_iterator li=lvi->begin(); + li!=lvi->end(); ++li) { + const int old_sz = bnds.size(); + const int this_sz = li->size(); + bnds -= *li; + const int new_sz = bnds.size(); + assert (new_sz + this_sz == old_sz); + } + } + } + + } // for ml + } // for c + } // for rl + for (typename list<ggf<D>*>::iterator f=gfs.begin(); f!=gfs.end(); ++f) { (*f)->recompose(); diff --git a/Carpet/CarpetLib/src/gdata.cc b/Carpet/CarpetLib/src/gdata.cc index 86f5f89c3..f17534838 100644 --- a/Carpet/CarpetLib/src/gdata.cc +++ b/Carpet/CarpetLib/src/gdata.cc @@ -1,4 +1,4 @@ -// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gdata.cc,v 1.21 2003/01/03 15:49:36 schnetter Exp $ +// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gdata.cc,v 1.22 2003/02/25 22:57:00 schnetter Exp $ #include <assert.h> @@ -43,6 +43,8 @@ void gdata<D>::copy_from (const gdata* src, const ibbox& box) assert (all((box.lower()-extent().lower())%box.stride() == 0 && (box.lower()-src->extent().lower())%box.stride() == 0)); + if (box.empty()) return; + if (proc() == src->proc()) { // copy on same processor @@ -87,6 +89,9 @@ void gdata<D> assert (all(box.upper()<=srcs[t]->extent().upper())); } + assert (! box.empty()); + if (box.empty()) return; + if (proc() == srcs[0]->proc()) { // interpolate on same processor diff --git a/Carpet/CarpetLib/src/ggf.cc b/Carpet/CarpetLib/src/ggf.cc index fcc9e546b..79a7b105f 100644 --- a/Carpet/CarpetLib/src/ggf.cc +++ b/Carpet/CarpetLib/src/ggf.cc @@ -1,4 +1,4 @@ -// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/ggf.cc,v 1.20 2003/01/03 15:49:36 schnetter Exp $ +// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/ggf.cc,v 1.21 2003/02/25 22:57:00 schnetter Exp $ #include <assert.h> #include <stdlib.h> @@ -32,8 +32,6 @@ ggf<D>::ggf (const string name, th<D>& t, dh<D>& d, assert (t.h == &d.h); d.add(this); - -// recompose(); } // Destructors diff --git a/Carpet/CarpetLib/src/prolongate_3d_real8.F77 b/Carpet/CarpetLib/src/prolongate_3d_real8.F77 index 5ac6d8ce9..2dd5768e6 100644 --- a/Carpet/CarpetLib/src/prolongate_3d_real8.F77 +++ b/Carpet/CarpetLib/src/prolongate_3d_real8.F77 @@ -1,8 +1,21 @@ c -*-Fortran-*- -c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8.F77,v 1.9 2003/02/24 17:43:10 schnetter Exp $ +c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8.F77,v 1.10 2003/02/25 22:57:00 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 prolongate_3d_real8 ( $ src, srciext, srcjext, srckext, $ dst, dstiext, dstjext, dstkext, @@ -34,11 +47,13 @@ c bbox(:,3) is stride integer i0, j0, k0 integer fi, fj, fk integer ifac(2), jfac(2), kfac(2) -c$$$ integer ii, jj, kk + integer ii, jj, kk integer fac CCTK_REAL8 res integer d + character msg*1000 + do d=1,3 @@ -58,16 +73,16 @@ c$$$ integer ii, jj, kk $ .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 (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 - 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 end do if (srciext.ne.(srcbbox(1,2)-srcbbox(1,1))/srcbbox(1,3)+1 @@ -122,44 +137,48 @@ c Loop over fine region res = 0 -c$$$ do kk=1,2 -c$$$ do jj=1,2 -c$$$ do ii=1,2 -c$$$ -c$$$ fac = ifac(ii) * jfac(jj) * kfac(kk) -c$$$ -c$$$ if (fac.ne.0) then -c$$$ res = res + fac * src(i0+ii, j0+jj, k0+kk) -c$$$ end if -c$$$ -c$$$ end do -c$$$ end do -c$$$ end do - - fac = ifac(1) * jfac(1) * kfac(1) - if (fac.ne.0) res = res + fac * src(i0+1, j0+1, k0+1) - - fac = ifac(2) * jfac(1) * kfac(1) - if (fac.ne.0) res = res + fac * src(i0+2, j0+1, k0+1) - - fac = ifac(1) * jfac(2) * kfac(1) - if (fac.ne.0) res = res + fac * src(i0+1, j0+2, k0+1) - - fac = ifac(2) * jfac(2) * kfac(1) - if (fac.ne.0) res = res + fac * src(i0+2, j0+2, k0+1) - - fac = ifac(1) * jfac(1) * kfac(2) - if (fac.ne.0) res = res + fac * src(i0+1, j0+1, k0+2) - - fac = ifac(2) * jfac(1) * kfac(2) - if (fac.ne.0) res = res + fac * src(i0+2, j0+1, k0+2) - - fac = ifac(1) * jfac(2) * kfac(2) - if (fac.ne.0) res = res + fac * src(i0+1, j0+2, k0+2) + do kk=1,2 + do jj=1,2 + do ii=1,2 + + fac = ifac(ii) * jfac(jj) * kfac(kk) + + if (fac.ne.0) then + CHKIDX (i0+ii, j0+jj, k0+kk, \ + srciext,srcjext,srckext, "source") + res = res + fac * src(i0+ii, j0+jj, k0+kk) + end if + + end do + end do + end do - fac = ifac(2) * jfac(2) * kfac(2) - if (fac.ne.0) res = res + fac * src(i0+2, j0+2, k0+2) +c$$$ fac = ifac(1) * jfac(1) * kfac(1) +c$$$ if (fac.ne.0) res = res + fac * src(i0+1, j0+1, k0+1) +c$$$ +c$$$ fac = ifac(2) * jfac(1) * kfac(1) +c$$$ if (fac.ne.0) res = res + fac * src(i0+2, j0+1, k0+1) +c$$$ +c$$$ fac = ifac(1) * jfac(2) * kfac(1) +c$$$ if (fac.ne.0) res = res + fac * src(i0+1, j0+2, k0+1) +c$$$ +c$$$ fac = ifac(2) * jfac(2) * kfac(1) +c$$$ if (fac.ne.0) res = res + fac * src(i0+2, j0+2, k0+1) +c$$$ +c$$$ fac = ifac(1) * jfac(1) * kfac(2) +c$$$ if (fac.ne.0) res = res + fac * src(i0+1, j0+1, k0+2) +c$$$ +c$$$ fac = ifac(2) * jfac(1) * kfac(2) +c$$$ if (fac.ne.0) res = res + fac * src(i0+2, j0+1, k0+2) +c$$$ +c$$$ fac = ifac(1) * jfac(2) * kfac(2) +c$$$ if (fac.ne.0) res = res + fac * src(i0+1, j0+2, k0+2) +c$$$ +c$$$ fac = ifac(2) * jfac(2) * kfac(2) +c$$$ if (fac.ne.0) res = res + fac * src(i0+2, j0+2, k0+2) + CHKIDX (dstioff+i+1, dstjoff+j+1, dstkoff+k+1, \ + dstiext,dstjext,dstkext, "destination") dst (dstioff+i+1, dstjoff+j+1, dstkoff+k+1) = dstdiv * res end do diff --git a/Carpet/CarpetLib/src/prolongate_3d_real8_2tl.F77 b/Carpet/CarpetLib/src/prolongate_3d_real8_2tl.F77 index 8fab444ae..a89ccfbba 100644 --- a/Carpet/CarpetLib/src/prolongate_3d_real8_2tl.F77 +++ b/Carpet/CarpetLib/src/prolongate_3d_real8_2tl.F77 @@ -1,7 +1,20 @@ c -*-Fortran-*- -c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_2tl.F77,v 1.8 2003/02/24 17:43:10 schnetter Exp $ +c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_2tl.F77,v 1.9 2003/02/25 22:57:00 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 prolongate_3d_real8_2tl ( $ src1, t1, src2, t2, srciext, srcjext, srckext, @@ -48,6 +61,8 @@ c bbox(:,3) is stride CCTK_REAL8 res integer d + character msg*1000 + do d=1,3 @@ -67,16 +82,16 @@ c bbox(:,3) is stride $ .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 (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 - 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 end do if (srciext.ne.(srcbbox(1,2)-srcbbox(1,1))/srcbbox(1,3)+1 @@ -151,6 +166,8 @@ c Loop over fine region fac = ifac(ii) * jfac(jj) * kfac(kk) if (fac.ne.0) then + CHKIDX (i0+ii, j0+jj, k0+kk, \ + srciext,srcjext,srckext, "source") res = res $ + fac * s1fac * src1(i0+ii, j0+jj, k0+kk) $ + fac * s2fac * src2(i0+ii, j0+jj, k0+kk) @@ -160,6 +177,8 @@ c Loop over fine region end do end do + CHKIDX (dstioff+i+1, dstjoff+j+1, dstkoff+k+1, \ + dstiext,dstjext,dstkext, "destination") dst (dstioff+i+1, dstjoff+j+1, dstkoff+k+1) = dstdiv * res end do diff --git a/Carpet/CarpetLib/src/prolongate_3d_real8_2tl_o3.F77 b/Carpet/CarpetLib/src/prolongate_3d_real8_2tl_o3.F77 index f48345562..4d9f1f899 100644 --- a/Carpet/CarpetLib/src/prolongate_3d_real8_2tl_o3.F77 +++ b/Carpet/CarpetLib/src/prolongate_3d_real8_2tl_o3.F77 @@ -1,7 +1,20 @@ c -*-Fortran-*- -c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_2tl_o3.F77,v 1.10 2003/02/24 17:43:10 schnetter Exp $ +c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_2tl_o3.F77,v 1.11 2003/02/25 22:57:00 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 prolongate_3d_real8_2tl_o3 ( $ src1, t1, src2, t2, srciext, srcjext, srckext, @@ -50,6 +63,8 @@ c bbox(:,3) is stride CCTK_REAL8 res integer d + character msg*1000 + do d=1,3 @@ -69,23 +84,33 @@ c bbox(:,3) is stride $ .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 regkext = (regbbox(d,2) - regbbox(d,1)) / regbbox(d,3) + 1 dstkfac = srcbbox(d,3) / dstbbox(d,3) srckoff = (regbbox(d,1) - srcbbox(d,1)) / dstbbox(d,3) offsetlo = regbbox(d,3) - if (mod(srckoff + 0, dstkfac).eq.0) offsetlo = 0 + if (mod(srckoff + 0, dstkfac).eq.0) then + offsetlo = 0 + if (regkext.gt.1) then + offsetlo = regbbox(d,3) + end if + end if offsethi = regbbox(d,3) - if (mod(srckoff + regkext-1, dstkfac).eq.0) offsethi = 0 + if (mod(srckoff + regkext-1, dstkfac).eq.0) then + offsethi = 0 + if (regkext.gt.1) then + offsethi = regbbox(d,3) + end if + end if if (regbbox(d,1)-offsetlo.lt.srcbbox(d,1) $ .or. regbbox(d,2)+offsethi.gt.srcbbox(d,2) $ .or. regbbox(d,1).lt.dstbbox(d,1) $ .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 - 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 end do if (srciext.ne.(srcbbox(1,2)-srcbbox(1,1))/srcbbox(1,3)+1 @@ -166,6 +191,8 @@ c Loop over fine region fac = ifac(ii) * jfac(jj) * kfac(kk) if (fac.ne.0) then + CHKIDX (i0+ii-1, j0+jj-1, k0+kk-1, \ + srciext,srcjext,srckext, "source") res = res $ + fac * s1fac * src1(i0+ii-1, j0+jj-1, k0+kk-1) $ + fac * s2fac * src2(i0+ii-1, j0+jj-1, k0+kk-1) @@ -175,6 +202,8 @@ c Loop over fine region end do end do + CHKIDX (dstioff+i+1, dstjoff+j+1, dstkoff+k+1, \ + dstiext,dstjext,dstkext, "destination") dst (dstioff+i+1, dstjoff+j+1, dstkoff+k+1) = dstdiv * res end do diff --git a/Carpet/CarpetLib/src/prolongate_3d_real8_3tl.F77 b/Carpet/CarpetLib/src/prolongate_3d_real8_3tl.F77 index 2275e4b7f..820dbe89c 100644 --- a/Carpet/CarpetLib/src/prolongate_3d_real8_3tl.F77 +++ b/Carpet/CarpetLib/src/prolongate_3d_real8_3tl.F77 @@ -1,7 +1,20 @@ c -*-Fortran-*- -c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_3tl.F77,v 1.6 2003/02/24 17:43:10 schnetter Exp $ +c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_3tl.F77,v 1.7 2003/02/25 22:57:00 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 prolongate_3d_real8_3tl ( $ src1, t1, src2, t2, src3, t3, srciext, srcjext, srckext, @@ -50,6 +63,8 @@ c bbox(:,3) is stride CCTK_REAL8 res integer d + character msg*1000 + do d=1,3 @@ -69,16 +84,16 @@ c bbox(:,3) is stride $ .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 (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 - 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 end do if (srciext.ne.(srcbbox(1,2)-srcbbox(1,1))/srcbbox(1,3)+1 @@ -154,6 +169,8 @@ c Loop over fine region fac = ifac(ii) * jfac(jj) * kfac(kk) if (fac.ne.0) then + CHKIDX (i0+ii, j0+jj, k0+kk, \ + srciext,srcjext,srckext, "source") res = res $ + fac * s1fac * src1(i0+ii, j0+jj, k0+kk) $ + fac * s2fac * src2(i0+ii, j0+jj, k0+kk) @@ -164,6 +181,8 @@ c Loop over fine region end do end do + CHKIDX (dstioff+i+1, dstjoff+j+1, dstkoff+k+1, \ + dstiext,dstjext,dstkext, "destination") dst (dstioff+i+1, dstjoff+j+1, dstkoff+k+1) = dstdiv * res end do diff --git a/Carpet/CarpetLib/src/prolongate_3d_real8_3tl_o3.F77 b/Carpet/CarpetLib/src/prolongate_3d_real8_3tl_o3.F77 index b5bb94384..053328536 100644 --- a/Carpet/CarpetLib/src/prolongate_3d_real8_3tl_o3.F77 +++ b/Carpet/CarpetLib/src/prolongate_3d_real8_3tl_o3.F77 @@ -1,7 +1,20 @@ c -*-Fortran-*- -c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_3tl_o3.F77,v 1.10 2003/02/24 17:43:10 schnetter Exp $ +c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_3tl_o3.F77,v 1.11 2003/02/25 22:57:00 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 prolongate_3d_real8_3tl_o3 ( $ src1, t1, src2, t2, src3, t3, srciext, srcjext, srckext, @@ -52,6 +65,8 @@ c bbox(:,3) is stride CCTK_REAL8 res integer d + character msg*1000 + do d=1,3 @@ -71,23 +86,33 @@ c bbox(:,3) is stride $ .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 regkext = (regbbox(d,2) - regbbox(d,1)) / regbbox(d,3) + 1 dstkfac = srcbbox(d,3) / dstbbox(d,3) srckoff = (regbbox(d,1) - srcbbox(d,1)) / dstbbox(d,3) offsetlo = regbbox(d,3) - if (mod(srckoff + 0, dstkfac).eq.0) offsetlo = 0 + if (mod(srckoff + 0, dstkfac).eq.0) then + offsetlo = 0 + if (regkext.gt.1) then + offsetlo = regbbox(d,3) + end if + end if offsethi = regbbox(d,3) - if (mod(srckoff + regkext-1, dstkfac).eq.0) offsethi = 0 + if (mod(srckoff + regkext-1, dstkfac).eq.0) then + offsethi = 0 + if (regkext.gt.1) then + offsethi = regbbox(d,3) + end if + end if if (regbbox(d,1)-offsetlo.lt.srcbbox(d,1) $ .or. regbbox(d,2)+offsethi.gt.srcbbox(d,2) $ .or. regbbox(d,1).lt.dstbbox(d,1) $ .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 - 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 end do if (srciext.ne.(srcbbox(1,2)-srcbbox(1,1))/srcbbox(1,3)+1 @@ -169,6 +194,8 @@ c Loop over fine region fac = ifac(ii) * jfac(jj) * kfac(kk) if (fac.ne.0) then + CHKIDX (i0+ii-1, j0+jj-1, k0+kk-1, \ + srciext,srcjext,srckext, "source") res = res $ + fac * s1fac * src1(i0+ii-1, j0+jj-1, k0+kk-1) $ + fac * s2fac * src2(i0+ii-1, j0+jj-1, k0+kk-1) @@ -179,6 +206,8 @@ c Loop over fine region end do end do + CHKIDX (dstioff+i+1, dstjoff+j+1, dstkoff+k+1, \ + dstiext,dstjext,dstkext, "destination") dst (dstioff+i+1, dstjoff+j+1, dstkoff+k+1) = dstdiv * res end do diff --git a/Carpet/CarpetLib/src/prolongate_3d_real8_o3.F77 b/Carpet/CarpetLib/src/prolongate_3d_real8_o3.F77 index b728980b7..65d6f1f75 100644 --- a/Carpet/CarpetLib/src/prolongate_3d_real8_o3.F77 +++ b/Carpet/CarpetLib/src/prolongate_3d_real8_o3.F77 @@ -1,7 +1,20 @@ c -*-Fortran-*- -c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_o3.F77,v 1.7 2003/02/24 17:43:10 schnetter Exp $ +c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_o3.F77,v 1.8 2003/02/25 22:57:00 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 prolongate_3d_real8_o3 ( $ src, srciext, srcjext, srckext, @@ -41,6 +54,8 @@ c bbox(:,3) is stride CCTK_REAL8 res integer d + character msg*1000 + do d=1,3 @@ -60,23 +75,33 @@ c bbox(:,3) is stride $ .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 regkext = (regbbox(d,2) - regbbox(d,1)) / regbbox(d,3) + 1 dstkfac = srcbbox(d,3) / dstbbox(d,3) srckoff = (regbbox(d,1) - srcbbox(d,1)) / dstbbox(d,3) offsetlo = regbbox(d,3) - if (mod(srckoff + 0, dstkfac).eq.0) offsetlo = 0 + if (mod(srckoff + 0, dstkfac).eq.0) then + offsetlo = 0 + if (regkext.gt.1) then + offsetlo = regbbox(d,3) + end if + end if offsethi = regbbox(d,3) - if (mod(srckoff + regkext-1, dstkfac).eq.0) offsethi = 0 + if (mod(srckoff + regkext-1, dstkfac).eq.0) then + offsethi = 0 + if (regkext.gt.1) then + offsethi = regbbox(d,3) + end if + end if if (regbbox(d,1)-offsetlo.lt.srcbbox(d,1) $ .or. regbbox(d,2)+offsethi.gt.srcbbox(d,2) $ .or. regbbox(d,1).lt.dstbbox(d,1) $ .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 - 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 end do if (srciext.ne.(srcbbox(1,2)-srcbbox(1,1))/srcbbox(1,3)+1 @@ -144,6 +169,8 @@ c Loop over fine region fac = ifac(ii) * jfac(jj) * kfac(kk) if (fac.ne.0) then + CHKIDX (i0+ii-1, j0+jj-1, k0+kk-1, \ + srciext,srcjext,srckext, "source") res = res + fac * src(i0+ii-1, j0+jj-1, k0+kk-1) end if @@ -151,6 +178,8 @@ c Loop over fine region end do end do + CHKIDX (dstioff+i+1, dstjoff+j+1, dstkoff+k+1, \ + dstiext,dstjext,dstkext, "destination") dst (dstioff+i+1, dstjoff+j+1, dstkoff+k+1) = dstdiv * res end do diff --git a/Carpet/CarpetLib/src/restrict_3d_real8.F77 b/Carpet/CarpetLib/src/restrict_3d_real8.F77 index b846b18a6..782594a84 100644 --- a/Carpet/CarpetLib/src/restrict_3d_real8.F77 +++ b/Carpet/CarpetLib/src/restrict_3d_real8.F77 @@ -1,7 +1,20 @@ c -*-Fortran-*- -c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/restrict_3d_real8.F77,v 1.5 2003/02/24 17:43:10 schnetter Exp $ +c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/restrict_3d_real8.F77,v 1.6 2003/02/25 22:57:00 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, @@ -29,6 +42,8 @@ c bbox(:,3) is stride integer i, j, k integer d + character msg*1000 + do d=1,3 @@ -48,16 +63,16 @@ c bbox(:,3) is stride $ .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 (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 - 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 end do if (srciext.ne.(srcbbox(1,2)-srcbbox(1,1))/srcbbox(1,3)+1 @@ -94,6 +109,10 @@ c Loop over coarse region 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) |