53 int *ipc,
double *rpc,
54 double *ac,
double *cc,
55 double *x,
double *y) {
60 numdia = VAT(ipc, 11);
67 }
else if (numdia == 27) {
73 Vnm_print(2,
"MATVEC: invalid stencil type given...");
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) {
82 MAT2(ac, *nx * *ny * *nz, 1);
84 Vmatvec7_1s(nx, ny, nz,
87 RAT2(ac, 1, 2), RAT2(ac, 1, 3), RAT2(ac, 1, 4),
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) {
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);
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++) {
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);
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) {
134 MAT2(ac, *nx * *ny * *nz, 1);
136 Vmatvec27_1s(nx, ny, nz,
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),
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) {
159 double tmpO, tmpU, tmpD;
161 MAT3(cc, *nx, *ny, *nz);
162 MAT3(x, *nx, *ny, *nz);
163 MAT3(y, *nx, *ny, *nz);
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);
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);
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++) {
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);
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);
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);
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);
229 int *ipc,
double *rpc,
230 double *ac,
double *cc,
double *x,
double *y,
double *w1) {
235 numdia = VAT(ipc, 11);
238 Vnmatvec7(nx, ny, nz,
242 }
else if (numdia == 27) {
243 Vnmatvec27(nx, ny, nz,
248 Vnm_print(2,
"MATVEC: invalid stencil type given...");
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) {
259 MAT2(ac, *nx * *ny * *nz, 1);
263 Vnmatvecd7_1s(nx, ny, nz,
266 RAT2(ac, 1, 2), RAT2(ac, 1, 3), RAT2(ac, 1, 4),
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) {
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);
293 ipkey = VAT(ipc, 10);
294 Vc_vec(cc, x, w1, nx, ny, nz, &ipkey);
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++)
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)
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) {
318 MAT2(ac, *nx * *ny * *nz, 1);
323 Vnmatvecd27_1s(nx, ny, nz,
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),
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) {
347 double tmpO, tmpU, tmpD;
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);
371 ipkey = VAT(ipc, 10);
372 Vc_vec(cc, x, w1, nx, ny, nz, &ipkey);
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++) {
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);
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);
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);
412 VAT3(y, i, j, k) = tmpO + tmpU + tmpD
413 + VAT3(oC, i, j, k) * VAT3(x, i, j, k)
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) {
430 numdia = VAT(ipc, 11);
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);
436 Vnm_print(2,
"Vmresid: invalid stencil type given...\n");
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) {
447 MAT2(ac, *nx * *ny * *nz, 1);
450 Vmresid7_1s(nx, ny, nz,
452 RAT2(ac, 1,1), cc, fc,
453 RAT2(ac, 1,2), RAT2(ac, 1,3), RAT2(ac, 1,4),
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) {
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);
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);
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) {
499 MAT2(ac, *nx * *ny * *nz, 1);
502 Vmresid27_1s(nx,ny,nz,
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),
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) {
525 double tmpO, tmpU, tmpD;
527 MAT3(cc, *nx, *ny, *nz);
528 MAT3(fc, *nx, *ny, *nz);
529 MAT3(x, *nx, *ny, *nz);
530 MAT3(r, *nx, *ny, *nz);
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);
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);
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++) {
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);
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);
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);
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);
595 int *ipc,
double *rpc,
596 double *ac,
double *cc,
double *fc,
597 double *x,
double *r,
double *w1) {
602 numdia = VAT(ipc, 11);
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);
608 Vnm_print(2,
"Vnmresid: invalid stencil type given...\n");
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) {
619 MAT2(ac, *nx * *ny * *nz, 1);
622 Vnmresid7_1s(nx, ny, nz,
624 RAT2(ac, 1, 1), cc, fc,
625 RAT2(ac, 1, 2), RAT2(ac, 1, 3), RAT2(ac, 1, 4),
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) {
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);
649 ipkey = VAT(ipc, 10);
650 Vc_vec(cc, x, w1, nx, ny, nz, &ipkey);
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)
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) {
677 MAT2(ac, *nx * *ny * *nz, 1);
680 Vnmresid27_1s(nx, ny, nz,
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),
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) {
703 double tmpO, tmpU, tmpD;
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);
726 ipkey = VAT(ipc, 10);
727 Vc_vec(cc, x, w1, nx, ny, nz, &ipkey);
730 for (k=2; k<=*nz-1; k++) {
731 for (j=2; j<=*ny-1; j++) {
732 for (i=2; i<=*nx-1; i++) {
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);
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);
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);
769 - VAT3(oC, i, j, k) * VAT3(x, i, j, k)
778 VPUBLIC
void Vrestrc(
int *nxf,
int *nyf,
int *nzf,
780 double *xin,
double *xout,
double *pc) {
782 MAT2(pc, *nxc * *nyc * *nzc, 1 );
784 Vrestrc2(nxf, nyf, nzf,
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));
797 VEXTERNC
void Vrestrc2(
int *nxf,
int *nyf,
int *nzf,
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) {
811 double tmpO, tmpU, tmpD;
814 MAT3(xin, *nxf, *nyf, *nzf);
815 MAT3(xout, *nxc, *nyc, *nzc);
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);
823 MAT3(oPNE, *nxc, *nyc, *nzc);
824 MAT3(oPNW, *nxc, *nyc, *nzc);
825 MAT3(oPSE, *nxc, *nyc, *nzc);
826 MAT3(oPSW, *nxc, *nyc, *nzc);
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);
834 MAT3(uPNE, *nxc, *nyc, *nzc);
835 MAT3(uPNW, *nxc, *nyc, *nzc);
836 MAT3(uPSE, *nxc, *nyc, *nzc);
837 MAT3(uPSW, *nxc, *nyc, *nzc);
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);
845 MAT3(dPNE, *nxc, *nyc, *nzc);
846 MAT3(dPNW, *nxc, *nyc, *nzc);
847 MAT3(dPSE, *nxc, *nyc, *nzc);
848 MAT3(dPSW, *nxc, *nyc, *nzc);
853 dimfac = VPOW(2.0, idimenshun);
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;
860 for (j=2; j<=*nyc-1; j++) {
861 jj = (j - 1) * 2 + 1;
863 for (i=2; i<=*nxc-1; i++) {
864 ii = (i - 1) * 2 + 1;
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);
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);
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);
900 VAT3(xout, i, j, k) = tmpO + tmpU + tmpD;
912 int *nxf,
int *nyf,
int *nzf,
913 double *xin,
double *xout,
916 MAT2(pc, *nxc * *nyc * *nzc, 1);
918 VinterpPMG2(nxc, nyc, nzc,
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));
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) {
944 MAT3( xin, *nxc, *nyc, *nzc);
945 MAT3(xout, *nxf, *nyf, *nzf);
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);
953 MAT3(oPNE, *nxc, *nyc, *nzc);
954 MAT3(oPNW, *nxc, *nyc, *nzc);
955 MAT3(oPSE, *nxc, *nyc, *nzc);
956 MAT3(oPSW, *nxc, *nyc, *nzc);
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);
964 MAT3(uPNE, *nxc, *nyc, *nzc);
965 MAT3(uPNW, *nxc, *nyc, *nzc);
966 MAT3(uPSE, *nxc, *nyc, *nzc);
967 MAT3(uPSW, *nxc, *nyc, *nzc);
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);
975 MAT3(dPNE, *nxc, *nyc, *nzc);
976 MAT3(dPNW, *nxc, *nyc, *nzc);
977 MAT3(dPSE, *nxc, *nyc, *nzc);
978 MAT3(dPSW, *nxc, *nyc, *nzc);
988 for (k=1; k<=*nzf-2; k+=2) {
989 kk = (k - 1) / 2 + 1;
991 for (j=1; j<=*nyf-2; j+=2) {
992 jj = (j - 1) / 2 + 1;
994 for (i=1; i<=*nxf-2; i+=2) {
995 ii = (i - 1) / 2 + 1;
1002 VAT3(xout, i, j, k) = VAT3(xin, ii, jj, kk);
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);
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);
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);
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);
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);
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);
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);
1074 VPUBLIC
void Vextrac(
int *nxf,
int *nyf,
int *nzf,
1075 int *nxc,
int *nyc,
int *nzc,
1076 double *xin,
double *xout) {
1081 MAT3( xin, *nxf, *nyf, *nzf);
1082 MAT3(xout, *nxc, *nyc, *nzc);
1088 for (k=2; k<=*nzc-1; k++) {
1089 kk = (k - 1) * 2 + 1;
1091 for (j=2; j<=*nyc-1; j++) {
1092 jj = (j - 1) * 2 + 1;
1094 for (i=2; i<=*nxc-1; i++) {
1095 ii = (i - 1) * 2 + 1;
1098 VAT3(xout, i, j, k) = VAT3(xin, ii, jj, kk);
VPUBLIC void VfboundPMG00(int *nx, int *ny, int *nz, double *x)
Initialize a grid function to have a zero boundary value.
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.
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.
VPUBLIC void Vc_vec(double *coef, double *uin, double *uout, int *nx, int *ny, int *nz, int *ipkey)
Define the nonlinearity (vector version)
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.
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.
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.
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.
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.