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