File:	/home/sbrandt/cactus/Cactus/arrangements/Carpet/CarpetLib/src/prolongate_3d_dgfe_rf2.cc
1	:#include <cctk.h>
2	:#include <cctk_Parameters.h>
3	:
4	:#include <algorithm>
5	:#include <cassert>
6	:#include <cmath>
7	:#include <cstdlib>
8	:
9	:#include <hrscc.hh>
10	:
11	:#include "operator_prototypes_3d.hh"
12	:#include "typeprops.hh"
13	:
14	:using namespace std;
15	:#ifdef HRSCC_GLL_ELEMENT_HH
16	:using namespace hrscc;
17	:#endif
18	:
19	:
20	:
21	:namespace CarpetLib {
22	:
23	:#define SRCIND3(i,j,k)                          \
24	:  index3 (i, j, k,                              \
25	:          srcipadext, srcjpadext, srckpadext,   \
26	:          srciext, srcjext, srckext)
27	:#define DSTIND3(i,j,k)                          \
28	:  index3 (i, j, k,                              \
29	:          dstipadext, dstjpadext, dstkpadext,   \
30	:          dstiext, dstjext, dstkext)
31	:#define SRCOFF3(i,j,k)                          \
32	:  offset3 (i, j, k,                             \
33	:           srciext, srcjext, srckext)
34	:#define DSTOFF3(i,j,k)                          \
35	:  offset3 (i, j, k,                             \
36	:           dstiext, dstjext, dstkext)
37	:
38	:
39	:
40	:  template<typename T, int ORDER>
41	:  void
42	:  prolongate_3d_dgfe_rf2(T const *restrict const src,
43	:                         ivect3 const& restrict srcpadext,
44	:                         ivect3 const& restrict srcext,
45	:                         T *restrict const dst,
46	:                         ivect3 const& restrict dstpadext,
47	:                         ivect3 const& restrict dstext,
48	:                         ibbox3 const& restrict srcbbox,
49	:                         ibbox3 const& restrict dstbbox,
50	:                         ibbox3 const& restrict,
51	:                         ibbox3 const& restrict regbbox,
52	:                         void* const extraargs)
53	:  {
54	:#ifdef HRSCC_GLL_ELEMENT_HH
55	:    assert(not extraargs);
56	:
57	:    static_assert(ORDER>=0, "ORDER must be non-negative");
58	:
59	:
60	:
61	:    if (any(srcbbox.stride() <= regbbox.stride() or
62	:            dstbbox.stride() != regbbox.stride()))
63	:    {
64	:      CCTK_WARN(0, "Internal error: strides disagree");
65	:    }
66	:
67	:    if (any(srcbbox.stride() != reffact2 * dstbbox.stride())) {
68	:      CCTK_WARN(0, "Internal error: source strides are not twice the destination strides");
69	:    }
70	:
71	:    if (any(dstbbox.stride() % 2 != 0)) {
72	:      CCTK_WARN(0, "Internal error: destination strides are not even");
73	:    }
74	:
75	:    // This could be handled, but is likely to point to an error
76	:    // elsewhere
77	:    if (regbbox.empty()) {
78	:      CCTK_WARN(0, "Internal error: region extent is empty");
79	:    }
80	:
81	:
82	:
83	:    bvect3 const is_upper_face =
84	:      (regbbox.lower() - srcbbox.lower() + regbbox.stride() / 2) %
85	:      srcbbox.stride() != 0;
86	:
87	:    ivect3 const regext = regbbox.shape() / regbbox.stride();
88	:    assert(all((regbbox.lower() - srcbbox.lower() +
89	:                either(is_upper_face, -1, +1) * regbbox.stride() / 2) %
90	:               srcbbox.stride() == 0));
91	:    ivect3 const srcoff =
92	:      (regbbox.lower()- srcbbox.lower() +
93	:       either(is_upper_face, -1, +1) * regbbox.stride() / 2) /
94	:      srcbbox.stride();
95	:    assert(all((regbbox.lower() - dstbbox.lower()) % regbbox.stride() == 0));
96	:    ivect3 const dstoff =
97	:      (regbbox.lower() - dstbbox.lower()) / regbbox.stride();
98	:
99	:
100	:
101	:    if (not regbbox.is_contained_in(srcbbox) or
102	:        not regbbox.is_contained_in(dstbbox))
103	:    {
104	:      cerr << "ORDER=" << ORDER << "\n"
105	:           << "regbbox=" << regbbox << "\n"
106	:           << "dstbbox=" << dstbbox << "\n"
107	:           << "regbbox=" << regbbox << "\n"
108	:           << "srcbbox=" << srcbbox << "\n";
109	:      CCTK_WARN(0, "Internal error: region extent is not contained in array extent");
110	:    }
111	:
112	:
113	:
114	:    ptrdiff_t const srcipadext = srcpadext[0];
115	:    ptrdiff_t const srcjpadext = srcpadext[1];
116	:    ptrdiff_t const srckpadext = srcpadext[2];
117	:
118	:    ptrdiff_t const dstipadext = dstpadext[0];
119	:    ptrdiff_t const dstjpadext = dstpadext[1];
120	:    ptrdiff_t const dstkpadext = dstpadext[2];
121	:
122	:    ptrdiff_t const srciext = srcext[0];
123	:    ptrdiff_t const srcjext = srcext[1];
124	:    ptrdiff_t const srckext = srcext[2];
125	:
126	:    ptrdiff_t const dstiext = dstext[0];
127	:    ptrdiff_t const dstjext = dstext[1];
128	:    ptrdiff_t const dstkext = dstext[2];
129	:
130	:    ptrdiff_t const regiext = regext[0];
131	:    ptrdiff_t const regjext = regext[1];
132	:    ptrdiff_t const regkext = regext[2];
133	:
134	:    ptrdiff_t const srcioff = srcoff[0];
135	:    ptrdiff_t const srcjoff = srcoff[1];
136	:    ptrdiff_t const srckoff = srcoff[2];
137	:
138	:    ptrdiff_t const dstioff = dstoff[0];
139	:    ptrdiff_t const dstjoff = dstoff[1];
140	:    ptrdiff_t const dstkoff = dstoff[2];
141	:
142	:
143	:
144	:    if (regext[0] == 1) {
145	:      // 2D prolongation on x face
146	:
147	:      assert(not is_upper_face[1]);
148	:      assert(not is_upper_face[2]);
149	:
150	:      // Ensure we traverse an even integer number of elements
151	:      assert(regext[1] % (2*(ORDER+1)) == 0);
152	:      assert(regext[2] % (2*(ORDER+1)) == 0);
153	:
154	:      int const srcdi = 1;      // 2d face
155	:      int const srcdj = SRCOFF3(0,1,0) - SRCOFF3(0,0,0);
156	:      int const srcdk = SRCOFF3(0,0,1) - SRCOFF3(0,0,0);
157	:      int const dstdi = 1;      // 2d face
158	:      int const dstdj = DSTOFF3(0,1,0) - DSTOFF3(0,0,0);
159	:      int const dstdk = DSTOFF3(0,0,1) - DSTOFF3(0,0,0);
160	:      int const srcstr2d[2] = {srcdj, srcdk};
161	:      int const dststr2d[2] = {dstdj, dstdk};
162	:
163	:      // Loop over fine region
164	:      ptrdiff_t const i=0;
165	:#pragma omp parallel for collapse(2)
166	:      // Zwicky's Intel compiler 11.1 ices on ptrdiff_t
167	:      for (/*ptrdiff_t*/int k=0; k<regkext; k+=2*(ORDER+1)) {
168	:        for (/*ptrdiff_t*/int j=0; j<regjext; j+=2*(ORDER+1)) {
169	:          GLLElement<ORDER>::prolongate_2D
170	:            (&src[SRCIND3(srcioff+i, srcjoff+j, srckoff+k)], srcstr2d,
171	:             &dst[DSTIND3(dstioff+2*i, dstjoff+2*j, dstkoff+2*k)], dststr2d);
172	:        }
173	:      }
174	:
175	:    } else if (regext[1] == 1) {
176	:      // 2D prolongation on y face
177	:
178	:      assert(not is_upper_face[0]);
179	:      assert(not is_upper_face[2]);
180	:
181	:      // Ensure we traverse an even integer number of elements
182	:      assert(regext[0] % (2*(ORDER+1)) == 0);
183	:      assert(regext[2] % (2*(ORDER+1)) == 0);
184	:
185	:      // int const srcdi = SRCOFF3(1,0,0) - SRCOFF3(0,0,0);
186	:      int const srcdi = 1;
187	:      assert(srcdi == SRCOFF3(1,0,0) - SRCOFF3(0,0,0));
188	:      int const srcdj = 1;      // 2d face
189	:      int const srcdk = SRCOFF3(0,0,1) - SRCOFF3(0,0,0);
190	:      // int const dstdi = DSTOFF3(1,0,0) - DSTOFF3(0,0,0);
191	:      int const dstdi = 1;
192	:      assert(dstdi == DSTOFF3(1,0,0) - DSTOFF3(0,0,0));
193	:      int const dstdj = 1;      // 2d face
194	:      int const dstdk = DSTOFF3(0,0,1) - DSTOFF3(0,0,0);
195	:      int const srcstr2d[2]= {srcdi, srcdk};
196	:      int const dststr2d[2]= {dstdi, dstdk};
197	:
198	:      // Loop over fine region
199	:      ptrdiff_t const j=0;
200	:#pragma omp parallel for collapse(2)
201	:      // Zwicky's Intel compiler 11.1 ices on ptrdiff_t
202	:      for (/*ptrdiff_t*/int k=0; k<regkext; k+=2*(ORDER+1)) {
203	:        for (/*ptrdiff_t*/int i=0; i<regiext; i+=2*(ORDER+1)) {
204	:          GLLElement<ORDER>::prolongate_2D
205	:            (&src[SRCIND3(srcioff+i, srcjoff+j, srckoff+k)], srcstr2d,
206	:             &dst[DSTIND3(dstioff+2*i, dstjoff+2*j, dstkoff+2*k)], dststr2d);
207	:        }
208	:      }
209	:
210	:    } else if (regext[2] == 1) {
211	:      // 2D prolongation on z face
212	:
213	:      assert(not is_upper_face[0]);
214	:      assert(not is_upper_face[1]);
215	:
216	:      // Ensure we traverse an even integer number of elements
217	:      assert(regext[0] % (2*(ORDER+1)) == 0);
218	:      assert(regext[1] % (2*(ORDER+1)) == 0);
219	:
220	:      // int const srcdi = SRCOFF3(1,0,0) - SRCOFF3(0,0,0);
221	:      int const srcdi = 1;
222	:      assert(srcdi == SRCOFF3(1,0,0) - SRCOFF3(0,0,0));
223	:      int const srcdj = SRCOFF3(0,1,0) - SRCOFF3(0,0,0);
224	:      int const srcdk = 1;      // 2d face
225	:      // int const dstdi = DSTOFF3(1,0,0) - DSTOFF3(0,0,0);
226	:      int const dstdi = 1;
227	:      assert(dstdi == DSTOFF3(1,0,0) - DSTOFF3(0,0,0));
228	:      int const dstdj = DSTOFF3(0,1,0) - DSTOFF3(0,0,0);
229	:      int const dstdk = 1;      // 2d face
230	:      int const srcstr2d[2]= {srcdi, srcdj};
231	:      int const dststr2d[2]= {dstdi, dstdj};
232	:
233	:      // Loop over fine region
234	:      ptrdiff_t const k=0;
235	:#pragma omp parallel for collapse(2)
236	:      // Zwicky's Intel compiler 11.1 ices on ptrdiff_t
237	:      for (/*ptrdiff_t*/int j=0; j<regjext; j+=2*(ORDER+1)) {
238	:        for (/*ptrdiff_t*/int i=0; i<regiext; i+=2*(ORDER+1)) {
239	:          GLLElement<ORDER>::prolongate_2D
240	:            (&src[SRCIND3(srcioff+i, srcjoff+j, srckoff+k)], srcstr2d,
241	:             &dst[DSTIND3(dstioff+2*i, dstjoff+2*j, dstkoff+2*k)], dststr2d);
242	:        }
243	:      }
244	:
245	:    } else {
246	:      // 3D prolongation
247	:
248	:      assert(not any(is_upper_face));
249	:
250	:      // Ensure we traverse an even integer number of elements
251	:      assert(all(regext % (2*(ORDER+1)) == 0));
252	:
253	:      // int const srcdi = SRCOFF3(1,0,0) - SRCOFF3(0,0,0);
254	:      int const srcdi = 1;
255	:      assert(srcdi == SRCOFF3(1,0,0) - SRCOFF3(0,0,0));
256	:      int const srcdj = SRCOFF3(0,1,0) - SRCOFF3(0,0,0);
257	:      int const srcdk = SRCOFF3(0,0,1) - SRCOFF3(0,0,0);
258	:      // int const dstdi = DSTOFF3(1,0,0) - DSTOFF3(0,0,0);
259	:      int const dstdi = 1;
260	:      assert(dstdi == DSTOFF3(1,0,0) - DSTOFF3(0,0,0));
261	:      int const dstdj = DSTOFF3(0,1,0) - DSTOFF3(0,0,0);
262	:      int const dstdk = DSTOFF3(0,0,1) - DSTOFF3(0,0,0);
263	:      int const srcstr[3] = {srcdi, srcdj, srcdk};
264	:      int const dststr[3] = {dstdi, dstdj, dstdk};
265	:
266	:      // Loop over fine region
267	:#pragma omp parallel for collapse(3)
268	:      // Zwicky's Intel compiler 11.1 ices on ptrdiff_t
269	:      for (/*ptrdiff_t*/int k=0; k<regkext; k+=2*(ORDER+1)) {
270	:        for (/*ptrdiff_t*/int j=0; j<regjext; j+=2*(ORDER+1)) {
271	:          for (/*ptrdiff_t*/int i=0; i<regiext; i+=2*(ORDER+1)) {
272	:            GLLElement<ORDER>::prolongate_full
273	:              (&src[SRCIND3(srcioff+i, srcjoff+j, srckoff+k)], srcstr,
274	:               &dst[DSTIND3(dstioff+2*i, dstjoff+2*j, dstkoff+2*k)], dststr);
275	:          }
276	:        }
277	:      }
278	:
279	:    }
280	:
281	:#else
282	:    // HRSCCore is not available
283	:    assert(0);
284	:#endif
285	:  }
286	:
287	:
288	:
289	:#define TYPECASE(N,T)                                           \
290	:                                                                \
291	:  template                                                      \
292	:  void                                                          \
293	:  prolongate_3d_dgfe_rf2<T,5>(T const *restrict const src,      \
294	:                              ivect3 const& restrict srcpadext, \
295	:                              ivect3 const& restrict srcext,    \
296	:                              T *restrict const dst,            \
297	:                              ivect3 const& restrict dstpadext, \
298	:                              ivect3 const& restrict dstext,    \
299	:                              ibbox3 const& restrict srcbbox,   \
300	:                              ibbox3 const& restrict dstbbox,   \
301	:                              ibbox3 const& restrict,           \
302	:                              ibbox3 const& restrict regbbox,   \
303	:                              void *extraargs);
304	:
305	:#define CARPET_NO_INT
306	:#define CARPET_NO_COMPLEX
307	:#include "typecase.hh"
308	:#undef TYPECASE
309	:
310	:
311	:
312	:  template<>
313	:  void
314	:  prolongate_3d_dgfe_rf2<CCTK_COMPLEX,5>(CCTK_COMPLEX const *restrict const src,
315	:                                         ivect3 const& restrict srcpadext,
316	:                                         ivect3 const& restrict srcext,
317	:                                         CCTK_COMPLEX *restrict const dst,
318	:                                         ivect3 const& restrict dstpadext,
319	:                                         ivect3 const& restrict dstext,
320	:                                         ibbox3 const& restrict srcbbox,
321	:                                         ibbox3 const& restrict dstbbox,
322	:                                         ibbox3 const& restrict,
323	:                                         ibbox3 const& restrict regbbox,
324	:                                         void *extraargs)
325	:  {
326	:    assert(0);
327	:  }
328	:
329	:} // namespace CarpetLib