APBS  1.5
matvecd.c
1 
50 #include "matvecd.h"
51 
52 VPUBLIC void Vmatvec(int *nx, int *ny, int *nz,
53  int *ipc, double *rpc,
54  double *ac, double *cc,
55  double *x, double *y) {
56 
57  int numdia;
58 
59  // Do in one step
60  numdia = VAT(ipc, 11);
61 
62  if (numdia == 7) {
63  Vmatvec7(nx, ny, nz,
64  ipc, rpc,
65  ac, cc,
66  x, y);
67  } else if (numdia == 27) {
68  Vmatvec27(nx, ny, nz,
69  ipc, rpc,
70  ac, cc,
71  x, y);
72  } else {
73  Vnm_print(2, "MATVEC: invalid stencil type given...");
74  }
75 }
76 
77 VPUBLIC void Vmatvec7(int *nx, int *ny, int *nz,
78  int *ipc, double *rpc,
79  double *ac, double *cc,
80  double *x, double *y) {
81 
82  MAT2(ac, *nx * *ny * *nz, 1);
83 
84  Vmatvec7_1s(nx, ny, nz,
85  ipc, rpc,
86  RAT2(ac, 1, 1), cc,
87  RAT2(ac, 1, 2), RAT2(ac, 1, 3), RAT2(ac, 1, 4),
88  x, y);
89 }
90 
91 
92 
93 VEXTERNC void Vmatvec7_1s(int *nx, int *ny, int *nz,
94  int *ipc, double *rpc,
95  double *oC, double *cc,
96  double *oE, double *oN, double *uC,
97  double *x, double *y) {
98 
99  int i, j, k;
100 
101  MAT3(oE, *nx, *ny, *nz);
102  MAT3(oN, *nx, *ny, *nz);
103  MAT3(uC, *nx, *ny, *nz);
104  MAT3(cc, *nx, *ny, *nz);
105  MAT3(oC, *nx, *ny, *nz);
106  MAT3(x, *nx, *ny, *nz);
107  MAT3(y, *nx, *ny, *nz);
108 
109  // Do it
110  #pragma omp parallel for private(i, j, k)
111  for (k=2; k<=*nz-1; k++) {
112  for (j=2; j<=*ny-1; j++) {
113  for(i=2; i<=*nx-1; i++) {
114  VAT3(y, i, j, k) =
115  - VAT3( oN, i, j, k) * VAT3(x, i, j+1, k)
116  - VAT3( oN, i, j-1, k) * VAT3(x, i, j-1, k)
117  - VAT3( oE, i, j, k) * VAT3(x, i+1, j, k)
118  - VAT3( oE, i-1, j, k) * VAT3(x, i-1, j, k)
119  - VAT3( uC, i, j, k-1) * VAT3(x, i, j,k-1)
120  - VAT3( uC, i, j, k) * VAT3(x, i, j,k+1)
121  + (VAT3(oC, i, j, k) + VAT3(cc, i, j, k)) * VAT3(x, i, j, k);
122  }
123  }
124  }
125 }
126 
127 
128 
129 VPUBLIC void Vmatvec27(int *nx, int *ny, int *nz,
130  int *ipc, double *rpc,
131  double *ac, double *cc,
132  double *x, double *y) {
133 
134  MAT2(ac, *nx * *ny * *nz, 1);
135 
136  Vmatvec27_1s(nx, ny, nz,
137  ipc, rpc,
138  RAT2(ac, 1, 1), cc,
139  RAT2(ac, 1, 2), RAT2(ac, 1, 3), RAT2(ac, 1, 4),
140  RAT2(ac, 1, 5), RAT2(ac, 1, 6),
141  RAT2(ac, 1, 7), RAT2(ac, 1, 8), RAT2(ac, 1, 9), RAT2(ac, 1,10),
142  RAT2(ac, 1,11), RAT2(ac, 1,12), RAT2(ac, 1,13), RAT2(ac, 1,14),
143  x, y);
144 }
145 
146 
147 
148 VPUBLIC void Vmatvec27_1s(int *nx, int *ny, int *nz,
149  int *ipc, double *rpc,
150  double *oC, double *cc,
151  double *oE, double *oN, double *uC,
152  double *oNE, double *oNW,
153  double *uE, double *uW, double *uN, double *uS,
154  double *uNE, double *uNW, double *uSE, double *uSW,
155  double *x, double *y) {
156 
157  int i, j, k;
158 
159  double tmpO, tmpU, tmpD;
160 
161  MAT3(cc, *nx, *ny, *nz);
162  MAT3(x, *nx, *ny, *nz);
163  MAT3(y, *nx, *ny, *nz);
164 
165  MAT3(oC, *nx, *ny, *nz);
166  MAT3(oE, *nx, *ny, *nz);
167  MAT3(oN, *nx, *ny, *nz);
168  MAT3(oNE, *nx, *ny, *nz);
169  MAT3(oNW, *nx, *ny, *nz);
170 
171  MAT3(uC, *nx, *ny, *nz);
172  MAT3(uE, *nx, *ny, *nz);
173  MAT3(uW, *nx, *ny, *nz);
174  MAT3(uN, *nx, *ny, *nz);
175  MAT3(uS, *nx, *ny, *nz);
176  MAT3(uNE, *nx, *ny, *nz);
177  MAT3(uNW, *nx, *ny, *nz);
178  MAT3(uSE, *nx, *ny, *nz);
179  MAT3(uSW, *nx, *ny, *nz);
180 
181  // Do it
182  #pragma omp parallel for private(i, j, k, tmpO, tmpU, tmpD)
183  for (k=2; k<=*nz-1; k++) {
184  for (j=2; j<=*ny-1; j++) {
185  for(i=2; i<=*nx-1; i++) {
186  tmpO =
187  - VAT3( oN, i, j, k) * VAT3(x, i, j+1, k)
188  - VAT3( oN, i, j-1, k) * VAT3(x, i, j-1, k)
189  - VAT3( oE, i, j, k) * VAT3(x, i+1, j, k)
190  - VAT3( oE, i-1, j, k) * VAT3(x, i-1, j, k)
191  - VAT3( oNE, i, j, k) * VAT3(x, i+1, j+1, k)
192  - VAT3( oNW, i, j, k) * VAT3(x, i-1, j+1, k)
193  - VAT3( oNW, i+1, j-1, k) * VAT3(x, i+1, j-1, k)
194  - VAT3( oNE, i-1, j-1, k) * VAT3(x, i-1, j-1, k);
195 
196  tmpU =
197  - VAT3( uC, i, j, k) * VAT3(x, i, j, k+1)
198  - VAT3( uN, i, j, k) * VAT3(x, i, j+1, k+1)
199  - VAT3( uS, i, j, k) * VAT3(x, i, j-1, k+1)
200  - VAT3( uE, i, j, k) * VAT3(x, i+1, j, k+1)
201  - VAT3( uW, i, j, k) * VAT3(x, i-1, j, k+1)
202  - VAT3( uNE, i, j, k) * VAT3(x, i+1, j+1, k+1)
203  - VAT3( uNW, i, j, k) * VAT3(x, i-1, j+1, k+1)
204  - VAT3( uSE, i, j, k) * VAT3(x, i+1, j-1, k+1)
205  - VAT3( uSW, i, j, k) * VAT3(x, i-1, j-1, k+1);
206 
207  tmpD =
208  - VAT3( uC, i, j, k-1) * VAT3(x, i, j, k-1)
209  - VAT3( uS, i, j+1, k-1) * VAT3(x, i, j+1, k-1)
210  - VAT3( uN, i, j-1, k-1) * VAT3(x, i, j-1, k-1)
211  - VAT3( uW, i+1, j, k-1) * VAT3(x, i+1, j, k-1)
212  - VAT3( uE, i-1, j, k-1) * VAT3(x, i-1, j, k-1)
213  - VAT3( uSW, i+1, j+1, k-1) * VAT3(x, i+1, j+1, k-1)
214  - VAT3( uSE, i-1, j+1, k-1) * VAT3(x, i-1, j+1, k-1)
215  - VAT3( uNW, i+1, j-1, k-1) * VAT3(x, i+1, j-1, k-1)
216  - VAT3( uNE, i-1, j-1, k-1) * VAT3(x, i-1, j-1, k-1);
217 
218  VAT3(y, i, j, k) = tmpO + tmpU + tmpD
219  + (VAT3(oC, i, j, k) + VAT3(cc, i, j, k)) * VAT3(x, i, j, k);
220  }
221  }
222  }
223 }
224 
225 
226 
227 VEXTERNC void Vnmatvec(int *nx, int *ny, int *nz,
228  int *ipc, double *rpc,
229  double *ac, double *cc, double *x, double *y, double *w1) {
230 
231  int numdia;
232 
233  // Do in one step
234  numdia = VAT(ipc, 11);
235 
236  if (numdia == 7) {
237  Vnmatvec7(nx, ny, nz,
238  ipc, rpc,
239  ac, cc,
240  x, y, w1);
241  } else if (numdia == 27) {
242  Vnmatvec27(nx, ny, nz,
243  ipc, rpc,
244  ac, cc,
245  x, y, w1);
246  } else {
247  Vnm_print(2, "MATVEC: invalid stencil type given...");
248  }
249 }
250 
251 
252 
253 VPUBLIC void Vnmatvec7(int *nx, int *ny, int *nz,
254  int *ipc, double *rpc,
255  double *ac, double *cc,
256  double *x, double *y, double *w1) {
257 
258  MAT2(ac, *nx * *ny * *nz, 1);
259 
260  WARN_UNTESTED;
261 
262  Vnmatvecd7_1s(nx, ny, nz,
263  ipc, rpc,
264  RAT2(ac, 1, 1), cc,
265  RAT2(ac, 1, 2), RAT2(ac, 1, 3), RAT2(ac, 1, 4),
266  x, y, w1);
267 }
268 
269 
270 
271 VPUBLIC void Vnmatvecd7_1s(int *nx, int *ny, int *nz,
272  int *ipc, double *rpc,
273  double *oC, double *cc,
274  double *oE, double *oN, double *uC,
275  double *x, double *y, double *w1) {
276 
277  int i, j, k;
278  int ipkey;
279 
280  MAT3(oE, *nx, *ny, *nz);
281  MAT3(oN, *nx, *ny, *nz);
282  MAT3(uC, *nx, *ny, *nz);
283  MAT3(cc, *nx, *ny, *nz);
284  MAT3(oC, *nx, *ny, *nz);
285  MAT3( x, *nx, *ny, *nz);
286  MAT3( y, *nx, *ny, *nz);
287  MAT3(w1, *nx, *ny, *nz);
288 
289  WARN_UNTESTED;
290 
291  // first get vector nonlinear term to avoid subroutine calls
292  ipkey = VAT(ipc, 10);
293  Vc_vec(cc, x, w1, nx, ny, nz, &ipkey);
294 
295  // The operator
296  #pragma omp parallel for private(i, j, k)
297  for (k=2; k<=*nz-1; k++)
298  for (j=2; j<=*ny-1; j++)
299  for(i=2; i<=*nx-1; i++)
300  VAT3(y, i, j, k) =
301  - VAT3(oN, i, j, k) * VAT3(x, i, j+1, k)
302  - VAT3(oN, i, j-1, k) * VAT3(x, i, j-1, k)
303  - VAT3(oE, i, j, k) * VAT3(x, i+1, j, k)
304  - VAT3(oE, i-1, j, k) * VAT3(x, i-1, j, k)
305  - VAT3(uC, i, j, k-1) * VAT3(x, i, j, k-1)
306  - VAT3(uC, i, j, k) * VAT3(x, i, j, k+1)
307  + VAT3(oC, i, j, k) * VAT3(x, i, j, k)
308  + VAT3(w1, i, j, k);
309 }
310 
311 
312 VPUBLIC void Vnmatvec27(int *nx, int *ny, int *nz,
313  int *ipc, double *rpc,
314  double *ac, double *cc,
315  double *x, double *y, double *w1) {
316 
317  MAT2(ac, *nx * *ny * *nz, 1);
318 
319  WARN_UNTESTED;
320 
321  // Do in one step
322  Vnmatvecd27_1s(nx, ny, nz,
323  ipc, rpc,
324  RAT2(ac, 1, 1), cc,
325  RAT2(ac, 1, 2), RAT2(ac, 1, 3), RAT2(ac, 1, 4),
326  RAT2(ac, 1, 5), RAT2(ac, 1, 6),
327  RAT2(ac, 1, 7), RAT2(ac, 1, 8), RAT2(ac, 1, 9), RAT2(ac, 1,10),
328  RAT2(ac, 1,11), RAT2(ac, 1,12), RAT2(ac, 1,13), RAT2(ac, 1,14),
329  x, y, w1);
330 }
331 
332 
333 
334 VPUBLIC void Vnmatvecd27_1s(int *nx, int *ny, int *nz,
335  int *ipc, double *rpc,
336  double *oC, double *cc,
337  double *oE, double *oN, double *uC,
338  double *oNE, double *oNW,
339  double *uE, double *uW, double *uN, double *uS,
340  double *uNE, double *uNW, double *uSE, double *uSW,
341  double *x, double *y, double *w1) {
342 
343  int i, j, k;
344  int ipkey;
345 
346  double tmpO, tmpU, tmpD;
347 
348  MAT3( oE, *nx, *ny, *nz);
349  MAT3( oN, *nx, *ny, *nz);
350  MAT3( uC, *nx, *ny, *nz);
351  MAT3(oNE, *nx, *ny, *nz);
352  MAT3(oNW, *nx, *ny, *nz);
353  MAT3( uE, *nx, *ny, *nz);
354  MAT3( uW, *nx, *ny, *nz);
355  MAT3( uN, *nx, *ny, *nz);
356  MAT3( uS, *nx, *ny, *nz);
357  MAT3(uNE, *nx, *ny, *nz);
358  MAT3(uNW, *nx, *ny, *nz);
359  MAT3(uSE, *nx, *ny, *nz);
360  MAT3(uSW, *nx, *ny, *nz);
361  MAT3( cc, *nx, *ny, *nz);
362  MAT3( oC, *nx, *ny, *nz);
363  MAT3( x, *nx, *ny, *nz);
364  MAT3( y, *nx, *ny, *nz);
365  MAT3( w1, *nx, *ny, *nz);
366 
367  WARN_UNTESTED;
368 
369  // First get vector noNlinear term to avoid subroutine calls
370  ipkey = VAT(ipc, 10);
371  Vc_vec(cc, x, w1, nx, ny, nz, &ipkey);
372 
373  // The operator
374  #pragma omp parallel for private(i, j, k, tmpO, tmpU, tmpD)
375  for (k=2; k<=*nz-1; k++) {
376  for (j=2; j<=*ny-1; j++) {
377  for(i=2; i<=*nx-1; i++) {
378 
379  tmpO =
380  - VAT3( oN, i, j, k) * VAT3(x, i, j+1, k)
381  - VAT3( oN, i, j-1, k) * VAT3(x, i, j-1, k)
382  - VAT3( oE, i, j, k) * VAT3(x, i+1, j, k)
383  - VAT3( oE, i-1, j, k) * VAT3(x, i-1, j, k)
384  - VAT3(oNE, i, j, k) * VAT3(x, i+1, j+1, k)
385  - VAT3(oNW, i, j, k) * VAT3(x, i-1, j+1, k)
386  - VAT3(oNW, i+1, j-1, k) * VAT3(x, i+1, j-1, k)
387  - VAT3(oNE, i-1, j-1, k) * VAT3(x, i-1, j-1, k);
388 
389  tmpU =
390  - VAT3( uC, i, j, k) * VAT3(x, i, j, k+1)
391  - VAT3( uN, i, j, k) * VAT3(x, i, j+1, k+1)
392  - VAT3( uS, i, j, k) * VAT3(x, i, j-1, k+1)
393  - VAT3( uE, i, j, k) * VAT3(x, i+1, j, k+1)
394  - VAT3( uW, i, j, k) * VAT3(x, i-1, j, k+1)
395  - VAT3(uNE, i, j, k) * VAT3(x, i+1, j+1, k+1)
396  - VAT3(uNW, i, j, k) * VAT3(x, i-1, j+1, k+1)
397  - VAT3(uSE, i, j, k) * VAT3(x, i+1, j-1, k+1)
398  - VAT3(uSW, i, j, k) * VAT3(x, i-1, j-1, k+1);
399 
400  tmpD =
401  - VAT3( uC, i, j, k-1) * VAT3(x, i, j, k-1)
402  - VAT3( uS, i, j+1, k-1) * VAT3(x, i, j+1, k-1)
403  - VAT3( uN, i, j-1, k-1) * VAT3(x, i, j-1, k-1)
404  - VAT3( uW, i+1, j, k-1) * VAT3(x, i+1, j, k-1)
405  - VAT3( uE, i-1, j, k-1) * VAT3(x, i-1, j, k-1)
406  - VAT3(uSW, i+1, j+1, k-1) * VAT3(x, i+1, j+1, k-1)
407  - VAT3(uSE, i-1, j+1, k-1) * VAT3(x, i-1, j+1, k-1)
408  - VAT3(uNW, i+1, j-1, k-1) * VAT3(x, i+1, j-1, k-1)
409  - VAT3(uNE, i-1, j-1, k-1) * VAT3(x, i-1, j-1, k-1);
410 
411  VAT3(y, i, j, k) = tmpO + tmpU + tmpD
412  + VAT3(oC, i, j, k) * VAT3(x, i, j, k)
413  + VAT3(w1, i, j, k);
414  }
415  }
416  }
417 }
418 
419 
420 
421 VPUBLIC void Vmresid(int *nx, int *ny, int *nz,
422  int *ipc, double *rpc,
423  double *ac, double *cc, double *fc,
424  double *x, double *r) {
425 
426  int numdia;
427 
428  // Do in one step
429  numdia = VAT(ipc, 11);
430  if (numdia == 7) {
431  Vmresid7(nx, ny, nz, ipc, rpc, ac, cc, fc, x, r);
432  } else if (numdia == 27) {
433  Vmresid27(nx, ny, nz, ipc, rpc, ac, cc, fc, x, r);
434  } else {
435  Vnm_print(2, "Vmresid: invalid stencil type given...\n");
436  }
437 }
438 
439 
440 
441 VPUBLIC void Vmresid7(int *nx, int *ny, int *nz,
442  int *ipc, double *rpc,
443  double *ac, double *cc, double *fc,
444  double *x, double *r) {
445 
446  MAT2(ac, *nx * *ny * *nz, 1);
447 
448  // Do in one step
449  Vmresid7_1s(nx, ny, nz,
450  ipc, rpc,
451  RAT2(ac, 1,1), cc, fc,
452  RAT2(ac, 1,2), RAT2(ac, 1,3), RAT2(ac, 1,4),
453  x,r);
454 }
455 
456 VPUBLIC void Vmresid7_1s(int *nx, int *ny, int *nz,
457  int *ipc, double *rpc,
458  double *oC, double *cc, double *fc,
459  double *oE, double *oN, double *uC,
460  double *x, double *r) {
461 
462  int i, j, k;
463 
464  MAT3(oE, *nx, *ny, *nz);
465  MAT3(oN, *nx, *ny, *nz);
466  MAT3(uC, *nx, *ny, *nz);
467  MAT3(cc, *nx, *ny, *nz);
468  MAT3(fc, *nx, *ny, *nz);
469  MAT3(oC, *nx, *ny, *nz);
470  MAT3(x, *nx, *ny, *nz);
471  MAT3(r, *nx, *ny, *nz);
472 
473  // Do it
474  #pragma omp parallel for private(i, j, k)
475  for (k=2; k<=*nz-1; k++) {
476  for (j=2; j<=*ny-1; j++) {
477  for(i=2; i<=*nx-1; i++) {
478  VAT3(r, i,j,k) = VAT3(fc, i, j, k)
479  + VAT3( oN, i, j, k) * VAT3(x, i, j+1, k)
480  + VAT3( oN, i, j-1, k) * VAT3(x, i, j-1, k)
481  + VAT3( oE, i, j, k) * VAT3(x, i+1, j, k)
482  + VAT3( oE, i-1, j, k) * VAT3(x, i-1, j, k)
483  + VAT3( uC, i, j, k-1) * VAT3(x, i, j, k-1)
484  + VAT3( uC, i, j, k) * VAT3(x, i, j, k+1)
485  - (VAT3(oC, i, j, k) + VAT3(cc, i, j, k)) * VAT3(x, i, j, k);
486  }
487  }
488  }
489 }
490 
491 
492 
493 VPUBLIC void Vmresid27(int *nx, int *ny, int *nz,
494  int *ipc, double *rpc,
495  double *ac, double *cc, double *fc,
496  double *x, double *r) {
497 
498  MAT2(ac, *nx * *ny * *nz, 1);
499 
500  // Do in one step
501  Vmresid27_1s(nx,ny,nz,
502  ipc, rpc,
503  RAT2(ac, 1, 1), cc, fc,
504  RAT2(ac, 1, 2), RAT2(ac, 1, 3), RAT2(ac, 1, 4),
505  RAT2(ac, 1, 5), RAT2(ac, 1, 6),
506  RAT2(ac, 1, 7), RAT2(ac, 1, 8), RAT2(ac, 1, 9), RAT2(ac, 1,10),
507  RAT2(ac, 1,11), RAT2(ac, 1,12), RAT2(ac, 1,13), RAT2(ac, 1,14),
508  x,r);
509 }
510 
511 
512 
513 VPUBLIC void Vmresid27_1s(int *nx, int *ny, int *nz,
514  int *ipc, double *rpc,
515  double *oC, double *cc, double *fc,
516  double *oE, double *oN, double *uC,
517  double *oNE, double *oNW,
518  double *uE, double *uW, double *uN, double *uS,
519  double *uNE, double *uNW, double *uSE, double *uSW,
520  double *x, double *r) {
521 
522  int i, j, k;
523 
524  double tmpO, tmpU, tmpD;
525 
526  MAT3(cc, *nx, *ny, *nz);
527  MAT3(fc, *nx, *ny, *nz);
528  MAT3(x, *nx, *ny, *nz);
529  MAT3(r, *nx, *ny, *nz);
530 
531  MAT3(oC, *nx, *ny, *nz);
532  MAT3(oE, *nx, *ny, *nz);
533  MAT3(oN, *nx, *ny, *nz);
534  MAT3(oNE, *nx, *ny, *nz);
535  MAT3(oNW, *nx, *ny, *nz);
536 
537  MAT3(uC, *nx, *ny, *nz);
538  MAT3(uE, *nx, *ny, *nz);
539  MAT3(uW, *nx, *ny, *nz);
540  MAT3(uN, *nx, *ny, *nz);
541  MAT3(uS, *nx, *ny, *nz);
542  MAT3(uNE, *nx, *ny, *nz);
543  MAT3(uNW, *nx, *ny, *nz);
544  MAT3(uSE, *nx, *ny, *nz);
545  MAT3(uSW, *nx, *ny, *nz);
546 
547  #pragma omp parallel for private(i, j, k, tmpO, tmpU, tmpD)
548  for (k=2; k<=*nz-1; k++) {
549  for (j=2; j<=*ny-1; j++) {
550  for(i=2; i<=*nx-1; i++) {
551 
552  tmpO =
553  + VAT3( oN, i, j, k) * VAT3(x, i, j+1, k)
554  + VAT3( oN, i, j-1, k) * VAT3(x, i, j-1, k)
555  + VAT3( oE, i, j, k) * VAT3(x, i+1, j, k)
556  + VAT3( oE, i-1, j, k) * VAT3(x, i-1, j, k)
557  + VAT3( oNE, i, j, k) * VAT3(x, i+1, j+1, k)
558  + VAT3( oNW, i, j, k) * VAT3(x, i-1, j+1, k)
559  + VAT3( oNW, i+1, j-1, k) * VAT3(x, i+1, j-1, k)
560  + VAT3( oNE, i-1, j-1, k) * VAT3(x, i-1, j-1, k);
561 
562  tmpU =
563  + VAT3( uC, i, j, k) * VAT3(x, i, j, k+1)
564  + VAT3( uN, i, j, k) * VAT3(x, i, j+1, k+1)
565  + VAT3( uS, i, j, k) * VAT3(x, i, j-1, k+1)
566  + VAT3( uE, i, j, k) * VAT3(x, i+1, j, k+1)
567  + VAT3( uW, i, j, k) * VAT3(x, i-1, j, k+1)
568  + VAT3( uNE, i, j, k) * VAT3(x, i+1, j+1, k+1)
569  + VAT3( uNW, i, j, k) * VAT3(x, i-1, j+1, k+1)
570  + VAT3( uSE, i, j, k) * VAT3(x, i+1, j-1, k+1)
571  + VAT3( uSW, i, j, k) * VAT3(x, i-1, j-1, k+1);
572 
573  tmpD =
574  + VAT3( uC, i, j, k-1) * VAT3(x, i, j, k-1)
575  + VAT3( uS, i, j+1, k-1) * VAT3(x, i, j+1, k-1)
576  + VAT3( uN, i, j-1, k-1) * VAT3(x, i, j-1, k-1)
577  + VAT3( uW, i+1, j, k-1) * VAT3(x, i+1, j, k-1)
578  + VAT3( uE, i-1, j, k-1) * VAT3(x, i-1, j, k-1)
579  + VAT3( uSW, i+1, j+1, k-1) * VAT3(x, i+1, j+1, k-1)
580  + VAT3( uSE, i-1, j+1, k-1) * VAT3(x, i-1, j+1, k-1)
581  + VAT3( uNW, i+1, j-1, k-1) * VAT3(x, i+1, j-1, k-1)
582  + VAT3( uNE, i-1, j-1, k-1) * VAT3(x, i-1, j-1, k-1);
583 
584  VAT3(r, i, j, k) = VAT3(fc, i, j, k) + tmpO + tmpU + tmpD
585  - (VAT3(oC, i, j, k) + VAT3(cc, i, j, k)) * VAT3(x, i, j, k);
586  }
587  }
588  }
589 }
590 
591 
592 
593 VPUBLIC void Vnmresid(int *nx, int *ny, int *nz,
594  int *ipc, double *rpc,
595  double *ac, double *cc, double *fc,
596  double *x, double *r, double *w1) {
597 
598  int numdia;
599 
600  // Do in oNe step ***
601  numdia = VAT(ipc, 11);
602  if (numdia == 7) {
603  Vnmresid7(nx, ny, nz, ipc, rpc, ac, cc, fc, x, r, w1);
604  } else if (numdia == 27) {
605  Vnmresid27(nx, ny, nz, ipc, rpc, ac, cc, fc, x, r, w1);
606  } else {
607  Vnm_print(2, "Vnmresid: invalid stencil type given...\n");
608  }
609 }
610 
611 
612 
613 VPUBLIC void Vnmresid7(int *nx, int *ny, int *nz,
614  int *ipc, double *rpc,
615  double *ac, double *cc, double *fc,
616  double *x, double *r, double *w1) {
617 
618  MAT2(ac, *nx * *ny * *nz, 1);
619 
620  // Do in oNe step
621  Vnmresid7_1s(nx, ny, nz,
622  ipc, rpc,
623  RAT2(ac, 1, 1), cc, fc,
624  RAT2(ac, 1, 2), RAT2(ac, 1, 3), RAT2(ac, 1, 4),
625  x, r, w1);
626 }
627 
628 VPUBLIC void Vnmresid7_1s(int *nx, int *ny, int *nz,
629  int *ipc, double *rpc,
630  double *oC, double *cc, double *fc,
631  double *oE, double *oN, double *uC,
632  double *x, double *r, double *w1) {
633 
634  int i, j, k;
635  int ipkey;
636 
637  MAT3(oE, *nx, *ny, *nz);
638  MAT3(oN, *nx, *ny, *nz);
639  MAT3(uC, *nx, *ny, *nz);
640  MAT3(cc, *nx, *ny, *nz);
641  MAT3(fc, *nx, *ny, *nz);
642  MAT3(oC, *nx, *ny, *nz);
643  MAT3( x, *nx, *ny, *nz);
644  MAT3( r, *nx, *ny, *nz);
645  MAT3(w1, *nx, *ny, *nz);
646 
647  // First get vector nonlinear term to avoid subroutine calls
648  ipkey = VAT(ipc, 10);
649  Vc_vec(cc, x, w1, nx, ny, nz, &ipkey);
650 
651  // The residual
652  for (k=2; k<=*nz-1; k++) {
653  for (j=2; j<=*ny-1; j++) {
654  for (i=2; i<=*nx-1; i++) {
655  VAT3(r, i, j, k) = VAT3(fc, i, j, k)
656  + VAT3(oN, i, j, k) * VAT3(x, i, j+1, k)
657  + VAT3(oN, i, j-1, k) * VAT3(x, i, j-1, k)
658  + VAT3(oE, i, j, k) * VAT3(x, i+1, j, k)
659  + VAT3(oE, i-1, j, k) * VAT3(x, i-1, j, k)
660  + VAT3(uC, i, j, k-1) * VAT3(x, i, j, k-1)
661  + VAT3(uC, i, j, k) * VAT3(x, i, j, k+1)
662  - VAT3(oC, i, j, k) * VAT3(x, i, j, k)
663  - VAT3(w1, i, j, k);
664  }
665  }
666  }
667 }
668 
669 
670 
671 VPUBLIC void Vnmresid27(int *nx, int *ny, int *nz,
672  int *ipc, double *rpc,
673  double *ac, double *cc, double *fc,
674  double *x, double *r, double *w1) {
675 
676  MAT2(ac, *nx * *ny * *nz, 1);
677 
678  // Do in oNe step
679  Vnmresid27_1s(nx, ny, nz,
680  ipc, rpc,
681  RAT2(ac, 1, 1), cc, fc,
682  RAT2(ac, 1, 2), RAT2(ac, 1, 3), RAT2(ac, 1, 4),
683  RAT2(ac, 1, 5), RAT2(ac, 1, 6),
684  RAT2(ac, 1, 7), RAT2(ac, 1, 8), RAT2(ac, 1, 9), RAT2(ac, 1, 10),
685  RAT2(ac, 1, 11), RAT2(ac, 1, 12), RAT2(ac, 1, 13), RAT2(ac, 1, 14),
686  x, r, w1);
687 }
688 
689 
690 
691 VPUBLIC void Vnmresid27_1s(int *nx, int *ny, int *nz,
692  int *ipc, double *rpc,
693  double *oC, double *cc, double *fc,
694  double *oE, double *oN, double *uC,
695  double *oNE, double *oNW,
696  double *uE, double *uW, double *uN, double *uS,
697  double *uNE, double *uNW, double *uSE, double *uSW,
698  double *x, double *r, double *w1) {
699 
700  int i, j, k;
701  int ipkey;
702  double tmpO, tmpU, tmpD;
703 
704  MAT3( oC, *nx, *ny, *nz);
705  MAT3( cc, *nx, *ny, *nz);
706  MAT3( fc, *nx, *ny, *nz);
707  MAT3( oE, *nx, *ny, *nz);
708  MAT3( oN, *nx, *ny, *nz);
709  MAT3( uC, *nx, *ny, *nz);
710  MAT3(oNE, *nx, *ny, *nz);
711  MAT3(oNW, *nx, *ny, *nz);
712  MAT3( uE, *nx, *ny, *nz);
713  MAT3( uW, *nx, *ny, *nz);
714  MAT3( uN, *nx, *ny, *nz);
715  MAT3( uS, *nx, *ny, *nz);
716  MAT3(uNE, *nx, *ny, *nz);
717  MAT3(uNW, *nx, *ny, *nz);
718  MAT3(uSE, *nx, *ny, *nz);
719  MAT3(uSW, *nx, *ny, *nz);
720  MAT3( x, *nx, *ny, *nz);
721  MAT3( r, *nx, *ny, *nz);
722  MAT3( w1, *nx, *ny, *nz);
723 
724  // First get vector noNlinear term to avoid subroutine calls
725  ipkey = VAT(ipc, 10);
726  Vc_vec(cc, x, w1, nx, ny, nz, &ipkey);
727 
728  // The residual
729  for (k=2; k<=*nz-1; k++) {
730  for (j=2; j<=*ny-1; j++) {
731  for (i=2; i<=*nx-1; i++) {
732 
733  tmpO =
734  + VAT3( oN, i, j, k) * VAT3(x, i, j+1, k)
735  + VAT3( oN, i, j-1, k) * VAT3(x, i, j-1, k)
736  + VAT3( oE, i, j, k) * VAT3(x, i+1, j, k)
737  + VAT3( oE, i-1, j, k) * VAT3(x, i-1, j, k)
738  + VAT3(oNE, i, j, k) * VAT3(x, i+1, j+1, k)
739  + VAT3(oNW, i, j, k) * VAT3(x, i-1, j+1, k)
740  + VAT3(oNW, i+1, j-1, k) * VAT3(x, i+1, j-1, k)
741  + VAT3(oNE, i-1, j-1, k) * VAT3(x, i-1, j-1, k);
742 
743  tmpU =
744  + VAT3( uC, i, j, k) * VAT3(x, i, j, k+1)
745  + VAT3( uN, i, j, k) * VAT3(x, i, j+1, k+1)
746  + VAT3( uS, i, j, k) * VAT3(x, i, j-1, k+1)
747  + VAT3( uE, i, j, k) * VAT3(x, i+1, j, k+1)
748  + VAT3( uW, i, j, k) * VAT3(x, i-1, j, k+1)
749  + VAT3(uNE, i, j, k) * VAT3(x, i+1, j+1, k+1)
750  + VAT3(uNW, i, j, k) * VAT3(x, i-1, j+1, k+1)
751  + VAT3(uSE, i, j, k) * VAT3(x, i+1, j-1, k+1)
752  + VAT3(uSW, i, j, k) * VAT3(x, i-1, j-1, k+1);
753 
754  tmpD =
755  + VAT3( uC, i, j, k-1) * VAT3(x, i, j, k-1)
756  + VAT3( uS, i, j+1, k-1) * VAT3(x, i, j+1, k-1)
757  + VAT3( uN, i, j-1, k-1) * VAT3(x, i, j-1, k-1)
758  + VAT3( uW, i+1, j, k-1) * VAT3(x, i+1, j, k-1)
759  + VAT3( uE, i-1, j, k-1) * VAT3(x, i-1, j, k-1)
760  + VAT3(uSW, i+1, j+1, k-1) * VAT3(x, i+1, j+1, k-1)
761  + VAT3(uSE, i-1, j+1, k-1) * VAT3(x, i-1, j+1, k-1)
762  + VAT3(uNW, i+1, j-1, k-1) * VAT3(x, i+1, j-1, k-1)
763  + VAT3(uNE, i-1, j-1, k-1) * VAT3(x, i-1, j-1, k-1);
764 
765  VAT3(r, i, j, k) =
766  + tmpO + tmpU + tmpD
767  + VAT3(fc, i, j, k)
768  - VAT3(oC, i, j, k) * VAT3(x, i, j, k)
769  - VAT3(w1, i, j, k);
770  }
771  }
772  }
773 }
774 
775 
776 
777 VPUBLIC void Vrestrc(int *nxf, int *nyf, int *nzf,
778  int *nxc, int *nyc, int *nzc,
779  double *xin, double *xout, double *pc) {
780 
781  MAT2(pc, *nxc * *nyc * *nzc, 1 );
782 
783  Vrestrc2(nxf, nyf, nzf,
784  nxc, nyc, nzc,
785  xin, xout,
786  RAT2(pc, 1, 1), RAT2(pc, 1, 2), RAT2(pc, 1, 3), RAT2(pc, 1, 4), RAT2(pc, 1, 5),
787  RAT2(pc, 1, 6), RAT2(pc, 1, 7), RAT2(pc, 1, 8), RAT2(pc, 1, 9),
788  RAT2(pc, 1,10), RAT2(pc, 1,11), RAT2(pc, 1,12), RAT2(pc, 1,13), RAT2(pc, 1,14),
789  RAT2(pc, 1,15), RAT2(pc, 1,16), RAT2(pc, 1,17), RAT2(pc, 1,18),
790  RAT2(pc, 1,19), RAT2(pc, 1,20), RAT2(pc, 1,21), RAT2(pc, 1,22), RAT2(pc, 1,23),
791  RAT2(pc, 1,24), RAT2(pc, 1,25), RAT2(pc, 1,26), RAT2(pc, 1,27));
792 }
793 
794 
795 
796 VEXTERNC void Vrestrc2(int *nxf, int *nyf, int *nzf,
797  int *nxc, int *nyc, int *nzc,
798  double *xin, double *xout,
799  double *oPC, double *oPN, double *oPS, double *oPE, double *oPW,
800  double *oPNE, double *oPNW, double *oPSE, double *oPSW,
801  double *uPC, double *uPN, double *uPS, double *uPE, double *uPW,
802  double *uPNE, double *uPNW, double *uPSE, double *uPSW,
803  double *dPC, double *dPN, double *dPS, double *dPE, double *dPW,
804  double *dPNE, double *dPNW, double *dPSE, double *dPSW) {
805 
806  int i, j, k;
807  int ii, jj, kk;
808  int idimenshun = 3;
809 
810  double tmpO, tmpU, tmpD;
811  double dimfac;
812 
813  MAT3(xin, *nxf, *nyf, *nzf);
814  MAT3(xout, *nxc, *nyc, *nzc);
815 
816  MAT3(oPC, *nxc, *nyc, *nzc);
817  MAT3(oPN, *nxc, *nyc, *nzc);
818  MAT3(oPS, *nxc, *nyc, *nzc);
819  MAT3(oPE, *nxc, *nyc, *nzc);
820  MAT3(oPW, *nxc, *nyc, *nzc);
821 
822  MAT3(oPNE, *nxc, *nyc, *nzc);
823  MAT3(oPNW, *nxc, *nyc, *nzc);
824  MAT3(oPSE, *nxc, *nyc, *nzc);
825  MAT3(oPSW, *nxc, *nyc, *nzc);
826 
827  MAT3(uPC, *nxc, *nyc, *nzc);
828  MAT3(uPN, *nxc, *nyc, *nzc);
829  MAT3(uPS, *nxc, *nyc, *nzc);
830  MAT3(uPE, *nxc, *nyc, *nzc);
831  MAT3(uPW, *nxc, *nyc, *nzc);
832 
833  MAT3(uPNE, *nxc, *nyc, *nzc);
834  MAT3(uPNW, *nxc, *nyc, *nzc);
835  MAT3(uPSE, *nxc, *nyc, *nzc);
836  MAT3(uPSW, *nxc, *nyc, *nzc);
837 
838  MAT3(dPC, *nxc, *nyc, *nzc);
839  MAT3(dPN, *nxc, *nyc, *nzc);
840  MAT3(dPS, *nxc, *nyc, *nzc);
841  MAT3(dPE, *nxc, *nyc, *nzc);
842  MAT3(dPW, *nxc, *nyc, *nzc);
843 
844  MAT3(dPNE, *nxc, *nyc, *nzc);
845  MAT3(dPNW, *nxc, *nyc, *nzc);
846  MAT3(dPSE, *nxc, *nyc, *nzc);
847  MAT3(dPSW, *nxc, *nyc, *nzc);
848 
849  // Verify correctness of the input boundary points
850  VfboundPMG00(nxf, nyf, nzf, xin);
851 
852  dimfac = VPOW(2.0, idimenshun);
853 
854  // Handle the interior points as average of 5 finer grid pts ***
855  #pragma omp parallel for private(k, kk, j, jj, i, ii, tmpO, tmpU, tmpD)
856  for (k=2; k<=*nzc-1; k++) {
857  kk = (k - 1) * 2 + 1;
858 
859  for (j=2; j<=*nyc-1; j++) {
860  jj = (j - 1) * 2 + 1;
861 
862  for (i=2; i<=*nxc-1; i++) {
863  ii = (i - 1) * 2 + 1;
864 
865  // Compute the restriction
866  tmpO =
867  + VAT3( oPC, i, j, k) * VAT3(xin, ii, jj, kk)
868  + VAT3( oPN, i, j, k) * VAT3(xin, ii, jj+1, kk)
869  + VAT3( oPS, i, j, k) * VAT3(xin, ii, jj-1, kk)
870  + VAT3( oPE, i, j, k) * VAT3(xin, ii+1, jj, kk)
871  + VAT3( oPW, i, j, k) * VAT3(xin, ii-1, jj, kk)
872  + VAT3(oPNE, i, j, k) * VAT3(xin, ii+1, jj+1, kk)
873  + VAT3(oPNW, i, j, k) * VAT3(xin, ii-1, jj+1, kk)
874  + VAT3(oPSE, i, j, k) * VAT3(xin, ii+1, jj-1, kk)
875  + VAT3(oPSW, i, j, k) * VAT3(xin, ii-1, jj-1, kk);
876 
877  tmpU =
878  + VAT3( uPC, i, j, k) * VAT3(xin, ii, jj, kk+1)
879  + VAT3( uPN, i, j, k) * VAT3(xin, ii, jj+1, kk+1)
880  + VAT3( uPS, i, j, k) * VAT3(xin, ii, jj-1, kk+1)
881  + VAT3( uPE, i, j, k) * VAT3(xin, ii+1, jj, kk+1)
882  + VAT3( uPW, i, j, k) * VAT3(xin, ii-1, jj, kk+1)
883  + VAT3(uPNE, i, j, k) * VAT3(xin, ii+1, jj+1, kk+1)
884  + VAT3(uPNW, i, j, k) * VAT3(xin, ii-1, jj+1, kk+1)
885  + VAT3(uPSE, i, j, k) * VAT3(xin, ii+1, jj-1, kk+1)
886  + VAT3(uPSW, i, j, k) * VAT3(xin, ii-1, jj-1, kk+1);
887 
888  tmpD =
889  + VAT3( dPC, i, j, k) * VAT3(xin, ii, jj, kk-1)
890  + VAT3( dPN, i, j, k) * VAT3(xin, ii, jj+1, kk-1)
891  + VAT3( dPS, i, j, k) * VAT3(xin, ii, jj-1, kk-1)
892  + VAT3( dPE, i, j, k) * VAT3(xin, ii+1, jj, kk-1)
893  + VAT3( dPW, i, j, k) * VAT3(xin, ii-1, jj, kk-1)
894  + VAT3(dPNE, i, j, k) * VAT3(xin, ii+1, jj+1, kk-1)
895  + VAT3(dPNW, i, j, k) * VAT3(xin, ii-1, jj+1, kk-1)
896  + VAT3(dPSE, i, j, k) * VAT3(xin, ii+1, jj-1, kk-1)
897  + VAT3(dPSW, i, j, k) * VAT3(xin, ii-1, jj-1, kk-1);
898 
899  VAT3(xout, i, j, k) = tmpO + tmpU + tmpD;
900  }
901  }
902  }
903 
904  // Verify correctness of the output boundary points
905  VfboundPMG00(nxc, nyc, nzc, xout);
906 }
907 
908 
909 
910 VPUBLIC void VinterpPMG(int *nxc, int *nyc, int *nzc,
911  int *nxf, int *nyf, int *nzf,
912  double *xin, double *xout,
913  double *pc) {
914 
915  MAT2(pc, *nxc * *nyc * *nzc, 1);
916 
917  VinterpPMG2(nxc, nyc, nzc,
918  nxf, nyf, nzf,
919  xin, xout,
920  RAT2(pc, 1, 1), RAT2(pc, 1, 2), RAT2(pc, 1, 3), RAT2(pc, 1, 4), RAT2(pc, 1, 5),
921  RAT2(pc, 1, 6), RAT2(pc, 1, 7), RAT2(pc, 1, 8), RAT2(pc, 1, 9),
922  RAT2(pc, 1,10), RAT2(pc, 1,11), RAT2(pc, 1,12), RAT2(pc, 1,13), RAT2(pc, 1,14),
923  RAT2(pc, 1,15), RAT2(pc, 1,16), RAT2(pc, 1,17), RAT2(pc, 1,18),
924  RAT2(pc, 1,19), RAT2(pc, 1,20), RAT2(pc, 1,21), RAT2(pc, 1,22), RAT2(pc, 1,23),
925  RAT2(pc, 1,24), RAT2(pc, 1,25), RAT2(pc, 1,26), RAT2(pc, 1,27));
926 }
927 
928 
929 
930 VPUBLIC void VinterpPMG2(int *nxc, int *nyc, int *nzc,
931  int *nxf, int *nyf, int *nzf,
932  double *xin, double *xout,
933  double *oPC, double *oPN, double *oPS, double *oPE, double *oPW,
934  double *oPNE, double *oPNW, double *oPSE, double *oPSW,
935  double *uPC, double *uPN, double *uPS, double *uPE, double *uPW,
936  double *uPNE, double *uPNW, double *uPSE, double *uPSW,
937  double *dPC, double *dPN, double *dPS, double *dPE, double *dPW,
938  double *dPNE, double *dPNW, double *dPSE, double *dPSW) {
939 
940  int i, j, k;
941  int ii, jj, kk;
942 
943  MAT3( xin, *nxc, *nyc, *nzc);
944  MAT3(xout, *nxf, *nyf, *nzf);
945 
946  MAT3( oPC, *nxc, *nyc, *nzc);
947  MAT3( oPN, *nxc, *nyc, *nzc);
948  MAT3( oPS, *nxc, *nyc, *nzc);
949  MAT3( oPE, *nxc, *nyc, *nzc);
950  MAT3( oPW, *nxc, *nyc, *nzc);
951 
952  MAT3(oPNE, *nxc, *nyc, *nzc);
953  MAT3(oPNW, *nxc, *nyc, *nzc);
954  MAT3(oPSE, *nxc, *nyc, *nzc);
955  MAT3(oPSW, *nxc, *nyc, *nzc);
956 
957  MAT3( uPC, *nxc, *nyc, *nzc);
958  MAT3( uPN, *nxc, *nyc, *nzc);
959  MAT3( uPS, *nxc, *nyc, *nzc);
960  MAT3( uPE, *nxc, *nyc, *nzc);
961  MAT3( uPW, *nxc, *nyc, *nzc);
962 
963  MAT3(uPNE, *nxc, *nyc, *nzc);
964  MAT3(uPNW, *nxc, *nyc, *nzc);
965  MAT3(uPSE, *nxc, *nyc, *nzc);
966  MAT3(uPSW, *nxc, *nyc, *nzc);
967 
968  MAT3( dPC, *nxc, *nyc, *nzc);
969  MAT3( dPN, *nxc, *nyc, *nzc);
970  MAT3( dPS, *nxc, *nyc, *nzc);
971  MAT3( dPE, *nxc, *nyc, *nzc);
972  MAT3( dPW, *nxc, *nyc, *nzc);
973 
974  MAT3(dPNE, *nxc, *nyc, *nzc);
975  MAT3(dPNW, *nxc, *nyc, *nzc);
976  MAT3(dPSE, *nxc, *nyc, *nzc);
977  MAT3(dPSW, *nxc, *nyc, *nzc);
978 
979  /* *********************************************************************
980  * Setup
981  * *********************************************************************/
982 
983  // Verify correctness of the input boundary points ***
984  VfboundPMG00(nxc, nyc, nzc, xin);
985 
986  // Do it
987  for (k=1; k<=*nzf-2; k+=2) {
988  kk = (k - 1) / 2 + 1;
989 
990  for (j=1; j<=*nyf-2; j+=2) {
991  jj = (j - 1) / 2 + 1;
992 
993  for (i=1; i<=*nxf-2; i+=2) {
994  ii = (i - 1) / 2 + 1;
995 
996  /* ******************************************************** *
997  * Type 1 -- Fine grid points common to a coarse grid point *
998  * ******************************************************** */
999 
1000  // Copy coinciding points from coarse grid to fine grid
1001  VAT3(xout, i, j, k) = VAT3(xin, ii, jj, kk);
1002 
1003  /* ******************************************************** *
1004  * type 2 -- fine grid points common to a coarse grid plane *
1005  * ******************************************************** */
1006 
1007  // Fine grid pts common only to y-z planes on coarse grid
1008  // (intermediate pts between 2 grid points on x-row)
1009  VAT3(xout, i+1, j, k) = VAT3(oPE, ii, jj, kk) * VAT3(xin, ii, jj, kk)
1010  + VAT3(oPW, ii+1, jj, kk) * VAT3(xin, ii+1, jj, kk);
1011 
1012  // Fine grid pts common only to x-z planes on coarse grid
1013  // (intermediate pts between 2 grid points on a y-row)
1014  VAT3(xout, i, j+1, k) = VAT3(oPN, ii, jj, kk) * VAT3(xin, ii, jj, kk)
1015  + VAT3(oPS, ii, jj+1, kk) * VAT3(xin, ii, jj+1, kk);
1016 
1017  // Fine grid pts common only to x-y planes on coarse grid
1018  // (intermediate pts between 2 grid points on a z-row)
1019  VAT3(xout, i, j, k+1) = VAT3(uPC, ii, jj, kk) * VAT3(xin, ii, jj, kk)
1020  + VAT3(dPC, ii, jj, kk+1) * VAT3(xin, ii, jj, kk+1);
1021 
1022  /* ******************************************************* *
1023  * type 3 -- fine grid points common to a coarse grid line *
1024  * ******************************************************* */
1025 
1026  // Fine grid pts common only to z planes on coarse grid
1027  // (intermediate pts between 4 grid pts on the xy-plane
1028 
1029  VAT3(xout, i+1, j+1, k) = VAT3(oPNE, ii, jj, kk) * VAT3(xin, ii, jj, kk)
1030  + VAT3(oPNW, ii+1, jj, kk) * VAT3(xin, ii+1, jj, kk)
1031  + VAT3(oPSE, ii, jj+1, kk) * VAT3(xin, ii, jj+1, kk)
1032  + VAT3(oPSW, ii+1, jj+1, kk) * VAT3(xin, ii+1, jj+1, kk);
1033 
1034  // Fine grid pts common only to y planes on coarse grid
1035  // (intermediate pts between 4 grid pts on the xz-plane
1036  VAT3(xout, i+1, j, k+1) = VAT3(uPE, ii, jj, kk) * VAT3(xin, ii, jj, kk)
1037  + VAT3(uPW, ii+1, jj, kk) * VAT3(xin, ii+1, jj, kk)
1038  + VAT3(dPE, ii, jj, kk+1) * VAT3(xin, ii, jj, kk+1)
1039  + VAT3(dPW, ii+1, jj, kk+1) * VAT3(xin, ii+1, jj, kk+1);
1040 
1041  // Fine grid pts common only to x planes on coarse grid
1042  // (intermediate pts between 4 grid pts on the yz-plane***
1043  VAT3(xout, i, j+1, k+1) = VAT3(uPN, ii, jj, kk) * VAT3(xin, ii, jj, kk)
1044  + VAT3(uPS, ii, jj+1, kk) * VAT3(xin, ii, jj+1, kk)
1045  + VAT3(dPN, ii, jj,kk+1) * VAT3(xin, ii, jj, kk+1)
1046  + VAT3(dPS, ii, jj+1,kk+1) * VAT3(xin, ii, jj+1, kk+1);
1047 
1048  /* **************************************** *
1049  * type 4 -- fine grid points not common to *
1050  * coarse grid pts/lines/planes *
1051  * **************************************** */
1052 
1053  // Completely interior points
1054  VAT3(xout, i+1,j+1,k+1) =
1055  + VAT3(uPNE, ii, jj, kk) * VAT3(xin, ii, jj, kk)
1056  + VAT3(uPNW, ii+1, jj, kk) * VAT3(xin, ii+1, jj, kk)
1057  + VAT3(uPSE, ii, jj+1, kk) * VAT3(xin, ii, jj+1, kk)
1058  + VAT3(uPSW, ii+1, jj+1, kk) * VAT3(xin, ii+1, jj+1, kk)
1059  + VAT3(dPNE, ii, jj, kk+1) * VAT3(xin, ii, jj, kk+1)
1060  + VAT3(dPNW, ii+1, jj, kk+1) * VAT3(xin, ii+1, jj, kk+1)
1061  + VAT3(dPSE, ii, jj+1, kk+1) * VAT3(xin, ii, jj+1, kk+1)
1062  + VAT3(dPSW, ii+1, jj+1, kk+1) * VAT3(xin, ii+1, jj+1, kk+1);
1063  }
1064  }
1065  }
1066 
1067  // Verify correctness of the output boundary points ***
1068  VfboundPMG00(nxf, nyf, nzf, xout);
1069 }
1070 
1071 
1072 
1073 VPUBLIC void Vextrac(int *nxf, int *nyf, int *nzf,
1074  int *nxc, int *nyc, int *nzc,
1075  double *xin, double *xout) {
1076 
1077  int i, j, k;
1078  int ii, jj, kk;
1079 
1080  MAT3( xin, *nxf, *nyf, *nzf);
1081  MAT3(xout, *nxc, *nyc, *nzc);
1082 
1083  // Verify correctness of the input boundary points
1084  VfboundPMG00(nxf, nyf, nzf, xin);
1085 
1086  // Do it
1087  for (k=2; k<=*nzc-1; k++) {
1088  kk = (k - 1) * 2 + 1;
1089 
1090  for (j=2; j<=*nyc-1; j++) {
1091  jj = (j - 1) * 2 + 1;
1092 
1093  for (i=2; i<=*nxc-1; i++) {
1094  ii = (i - 1) * 2 + 1;
1095 
1096  // Compute the restriction
1097  VAT3(xout, i, j, k) = VAT3(xin, ii, jj, kk);
1098  }
1099  }
1100  }
1101 
1102  // Verify correctness of the output boundary points
1103  VfboundPMG00(nxc, nyc, nzc, xout);
1104 }
VPUBLIC void VfboundPMG00(int *nx, int *ny, int *nz, double *x)
Initialize a grid function to have a zero boundary value.
Definition: mikpckd.c:253
VPUBLIC void Vextrac(int *nxf, int *nyf, int *nzf, int *nxc, int *nyc, int *nzc, double *xin, double *xout)
Simple injection of a fine grid function into coarse grid.
Definition: matvecd.c:1073
int nzc
Definition: vpmgp.h:98
int nxc
Definition: vpmgp.h:96
VPUBLIC void VinterpPMG(int *nxc, int *nyc, int *nzc, int *nxf, int *nyf, int *nzf, double *xin, double *xout, double *pc)
Apply the prolongation operator.
Definition: matvecd.c:910
int ipkey
Definition: vpmgp.h:109
VPUBLIC void Vc_vec(double *coef, double *uin, double *uout, int *nx, int *ny, int *nz, int *ipkey)
Define the nonlinearity (vector version)
Definition: mypdec.c:116
VEXTERNC void Vnmatvec(int *nx, int *ny, int *nz, int *ipc, double *rpc, double *ac, double *cc, double *x, double *y, double *w1)
Break the matrix data-structure into diagonals and then call the matrix-vector routine.
Definition: matvecd.c:227
VPUBLIC void Vmatvec(int *nx, int *ny, int *nz, int *ipc, double *rpc, double *ac, double *cc, double *x, double *y)
Matrix-vector multiplication routines.
Definition: matvecd.c:52
int nyc
Definition: vpmgp.h:97
int ny
Definition: vpmgp.h:84
int nx
Definition: vpmgp.h:83
VPUBLIC void Vnmresid(int *nx, int *ny, int *nz, int *ipc, double *rpc, double *ac, double *cc, double *fc, double *x, double *r, double *w1)
Break the matrix data-structure into diagonals and then call the residual routine.
Definition: matvecd.c:593
VPUBLIC void Vmresid(int *nx, int *ny, int *nz, int *ipc, double *rpc, double *ac, double *cc, double *fc, double *x, double *r)
Break the matrix data-structure into diagonals and then call the residual routine.
Definition: matvecd.c:421
VPUBLIC void Vrestrc(int *nxf, int *nyf, int *nzf, int *nxc, int *nyc, int *nzc, double *xin, double *xout, double *pc)
Apply the restriction operator.
Definition: matvecd.c:777
int nz
Definition: vpmgp.h:85