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