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);
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) {
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);
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;
182 #pragma omp parallel for private(i, j, k, tmpO, tmpU, tmpD) 183 for (k=2; k<=*
nz-1; k++) {
184 for (j=2; j<=*
ny-1; j++) {
185 for(i=2; i<=*
nx-1; i++) {
187 - VAT3( oN, i, j, k) * VAT3(x, i, j+1, k)
188 - VAT3( oN, i, j-1, k) * VAT3(x, i, j-1, k)
189 - VAT3( oE, i, j, k) * VAT3(x, i+1, j, k)
190 - VAT3( oE, i-1, j, k) * VAT3(x, i-1, j, k)
191 - VAT3( oNE, i, j, k) * VAT3(x, i+1, j+1, k)
192 - VAT3( oNW, i, j, k) * VAT3(x, i-1, j+1, k)
193 - VAT3( oNW, i+1, j-1, k) * VAT3(x, i+1, j-1, k)
194 - VAT3( oNE, i-1, j-1, k) * VAT3(x, i-1, j-1, k);
197 - VAT3( uC, i, j, k) * VAT3(x, i, j, k+1)
198 - VAT3( uN, i, j, k) * VAT3(x, i, j+1, k+1)
199 - VAT3( uS, i, j, k) * VAT3(x, i, j-1, k+1)
200 - VAT3( uE, i, j, k) * VAT3(x, i+1, j, k+1)
201 - VAT3( uW, i, j, k) * VAT3(x, i-1, j, k+1)
202 - VAT3( uNE, i, j, k) * VAT3(x, i+1, j+1, k+1)
203 - VAT3( uNW, i, j, k) * VAT3(x, i-1, j+1, k+1)
204 - VAT3( uSE, i, j, k) * VAT3(x, i+1, j-1, k+1)
205 - VAT3( uSW, i, j, k) * VAT3(x, i-1, j-1, k+1);
208 - VAT3( uC, i, j, k-1) * VAT3(x, i, j, k-1)
209 - VAT3( uS, i, j+1, k-1) * VAT3(x, i, j+1, k-1)
210 - VAT3( uN, i, j-1, k-1) * VAT3(x, i, j-1, k-1)
211 - VAT3( uW, i+1, j, k-1) * VAT3(x, i+1, j, k-1)
212 - VAT3( uE, i-1, j, k-1) * VAT3(x, i-1, j, k-1)
213 - VAT3( uSW, i+1, j+1, k-1) * VAT3(x, i+1, j+1, k-1)
214 - VAT3( uSE, i-1, j+1, k-1) * VAT3(x, i-1, j+1, k-1)
215 - VAT3( uNW, i+1, j-1, k-1) * VAT3(x, i+1, j-1, k-1)
216 - VAT3( uNE, i-1, j-1, k-1) * VAT3(x, i-1, j-1, k-1);
218 VAT3(y, i, j, k) = tmpO + tmpU + tmpD
219 + (VAT3(oC, i, j, k) + VAT3(cc, i, j, k)) * VAT3(x, i, j, k);
228 int *ipc,
double *rpc,
229 double *ac,
double *cc,
double *x,
double *y,
double *w1) {
234 numdia = VAT(ipc, 11);
241 }
else if (numdia == 27) {
247 Vnm_print(2,
"MATVEC: invalid stencil type given...");
253 VPUBLIC
void Vnmatvec7(
int *
nx,
int *
ny,
int *
nz,
254 int *ipc,
double *rpc,
255 double *ac,
double *cc,
256 double *x,
double *y,
double *w1) {
258 MAT2(ac, *
nx * *
ny * *
nz, 1);
265 RAT2(ac, 1, 2), RAT2(ac, 1, 3), RAT2(ac, 1, 4),
271 VPUBLIC
void Vnmatvecd7_1s(
int *
nx,
int *
ny,
int *
nz,
272 int *ipc,
double *rpc,
273 double *oC,
double *cc,
274 double *oE,
double *oN,
double *uC,
275 double *x,
double *y,
double *w1) {
292 ipkey = VAT(ipc, 10);
296 #pragma omp parallel for private(i, j, k) 297 for (k=2; k<=*
nz-1; k++)
298 for (j=2; j<=*
ny-1; j++)
299 for(i=2; i<=*
nx-1; i++)
301 - VAT3(oN, i, j, k) * VAT3(x, i, j+1, k)
302 - VAT3(oN, i, j-1, k) * VAT3(x, i, j-1, k)
303 - VAT3(oE, i, j, k) * VAT3(x, i+1, j, k)
304 - VAT3(oE, i-1, j, k) * VAT3(x, i-1, j, k)
305 - VAT3(uC, i, j, k-1) * VAT3(x, i, j, k-1)
306 - VAT3(uC, i, j, k) * VAT3(x, i, j, k+1)
307 + VAT3(oC, i, j, k) * VAT3(x, i, j, k)
312 VPUBLIC
void Vnmatvec27(
int *
nx,
int *
ny,
int *
nz,
313 int *ipc,
double *rpc,
314 double *ac,
double *cc,
315 double *x,
double *y,
double *w1) {
317 MAT2(ac, *
nx * *
ny * *
nz, 1);
322 Vnmatvecd27_1s(
nx,
ny,
nz,
325 RAT2(ac, 1, 2), RAT2(ac, 1, 3), RAT2(ac, 1, 4),
326 RAT2(ac, 1, 5), RAT2(ac, 1, 6),
327 RAT2(ac, 1, 7), RAT2(ac, 1, 8), RAT2(ac, 1, 9), RAT2(ac, 1,10),
328 RAT2(ac, 1,11), RAT2(ac, 1,12), RAT2(ac, 1,13), RAT2(ac, 1,14),
334 VPUBLIC
void Vnmatvecd27_1s(
int *
nx,
int *
ny,
int *
nz,
335 int *ipc,
double *rpc,
336 double *oC,
double *cc,
337 double *oE,
double *oN,
double *uC,
338 double *oNE,
double *oNW,
339 double *uE,
double *uW,
double *uN,
double *uS,
340 double *uNE,
double *uNW,
double *uSE,
double *uSW,
341 double *x,
double *y,
double *w1) {
346 double tmpO, tmpU, tmpD;
370 ipkey = VAT(ipc, 10);
374 #pragma omp parallel for private(i, j, k, tmpO, tmpU, tmpD) 375 for (k=2; k<=*
nz-1; k++) {
376 for (j=2; j<=*
ny-1; j++) {
377 for(i=2; i<=*
nx-1; i++) {
380 - VAT3( oN, i, j, k) * VAT3(x, i, j+1, k)
381 - VAT3( oN, i, j-1, k) * VAT3(x, i, j-1, k)
382 - VAT3( oE, i, j, k) * VAT3(x, i+1, j, k)
383 - VAT3( oE, i-1, j, k) * VAT3(x, i-1, j, k)
384 - VAT3(oNE, i, j, k) * VAT3(x, i+1, j+1, k)
385 - VAT3(oNW, i, j, k) * VAT3(x, i-1, j+1, k)
386 - VAT3(oNW, i+1, j-1, k) * VAT3(x, i+1, j-1, k)
387 - VAT3(oNE, i-1, j-1, k) * VAT3(x, i-1, j-1, k);
390 - VAT3( uC, i, j, k) * VAT3(x, i, j, k+1)
391 - VAT3( uN, i, j, k) * VAT3(x, i, j+1, k+1)
392 - VAT3( uS, i, j, k) * VAT3(x, i, j-1, k+1)
393 - VAT3( uE, i, j, k) * VAT3(x, i+1, j, k+1)
394 - VAT3( uW, i, j, k) * VAT3(x, i-1, j, k+1)
395 - VAT3(uNE, i, j, k) * VAT3(x, i+1, j+1, k+1)
396 - VAT3(uNW, i, j, k) * VAT3(x, i-1, j+1, k+1)
397 - VAT3(uSE, i, j, k) * VAT3(x, i+1, j-1, k+1)
398 - VAT3(uSW, i, j, k) * VAT3(x, i-1, j-1, k+1);
401 - VAT3( uC, i, j, k-1) * VAT3(x, i, j, k-1)
402 - VAT3( uS, i, j+1, k-1) * VAT3(x, i, j+1, k-1)
403 - VAT3( uN, i, j-1, k-1) * VAT3(x, i, j-1, k-1)
404 - VAT3( uW, i+1, j, k-1) * VAT3(x, i+1, j, k-1)
405 - VAT3( uE, i-1, j, k-1) * VAT3(x, i-1, j, k-1)
406 - VAT3(uSW, i+1, j+1, k-1) * VAT3(x, i+1, j+1, k-1)
407 - VAT3(uSE, i-1, j+1, k-1) * VAT3(x, i-1, j+1, k-1)
408 - VAT3(uNW, i+1, j-1, k-1) * VAT3(x, i+1, j-1, k-1)
409 - VAT3(uNE, i-1, j-1, k-1) * VAT3(x, i-1, j-1, k-1);
411 VAT3(y, i, j, k) = tmpO + tmpU + tmpD
412 + VAT3(oC, i, j, k) * VAT3(x, i, j, k)
422 int *ipc,
double *rpc,
423 double *ac,
double *cc,
double *fc,
424 double *x,
double *r) {
429 numdia = VAT(ipc, 11);
431 Vmresid7(
nx,
ny,
nz, ipc, rpc, ac, cc, fc, x, r);
432 }
else if (numdia == 27) {
433 Vmresid27(
nx,
ny,
nz, ipc, rpc, ac, cc, fc, x, r);
435 Vnm_print(2,
"Vmresid: invalid stencil type given...\n");
441 VPUBLIC
void Vmresid7(
int *
nx,
int *
ny,
int *
nz,
442 int *ipc,
double *rpc,
443 double *ac,
double *cc,
double *fc,
444 double *x,
double *r) {
446 MAT2(ac, *
nx * *
ny * *
nz, 1);
451 RAT2(ac, 1,1), cc, fc,
452 RAT2(ac, 1,2), RAT2(ac, 1,3), RAT2(ac, 1,4),
456 VPUBLIC
void Vmresid7_1s(
int *
nx,
int *
ny,
int *
nz,
457 int *ipc,
double *rpc,
458 double *oC,
double *cc,
double *fc,
459 double *oE,
double *oN,
double *uC,
460 double *x,
double *r) {
474 #pragma omp parallel for private(i, j, k) 475 for (k=2; k<=*
nz-1; k++) {
476 for (j=2; j<=*
ny-1; j++) {
477 for(i=2; i<=*
nx-1; i++) {
478 VAT3(r, i,j,k) = VAT3(fc, i, j, k)
479 + VAT3( oN, i, j, k) * VAT3(x, i, j+1, k)
480 + VAT3( oN, i, j-1, k) * VAT3(x, i, j-1, k)
481 + VAT3( oE, i, j, k) * VAT3(x, i+1, j, k)
482 + VAT3( oE, i-1, j, k) * VAT3(x, i-1, j, k)
483 + VAT3( uC, i, j, k-1) * VAT3(x, i, j, k-1)
484 + VAT3( uC, i, j, k) * VAT3(x, i, j, k+1)
485 - (VAT3(oC, i, j, k) + VAT3(cc, i, j, k)) * VAT3(x, i, j, k);
493 VPUBLIC
void Vmresid27(
int *
nx,
int *
ny,
int *
nz,
494 int *ipc,
double *rpc,
495 double *ac,
double *cc,
double *fc,
496 double *x,
double *r) {
498 MAT2(ac, *
nx * *
ny * *
nz, 1);
503 RAT2(ac, 1, 1), cc, fc,
504 RAT2(ac, 1, 2), RAT2(ac, 1, 3), RAT2(ac, 1, 4),
505 RAT2(ac, 1, 5), RAT2(ac, 1, 6),
506 RAT2(ac, 1, 7), RAT2(ac, 1, 8), RAT2(ac, 1, 9), RAT2(ac, 1,10),
507 RAT2(ac, 1,11), RAT2(ac, 1,12), RAT2(ac, 1,13), RAT2(ac, 1,14),
513 VPUBLIC
void Vmresid27_1s(
int *
nx,
int *
ny,
int *
nz,
514 int *ipc,
double *rpc,
515 double *oC,
double *cc,
double *fc,
516 double *oE,
double *oN,
double *uC,
517 double *oNE,
double *oNW,
518 double *uE,
double *uW,
double *uN,
double *uS,
519 double *uNE,
double *uNW,
double *uSE,
double *uSW,
520 double *x,
double *r) {
524 double tmpO, tmpU, tmpD;
547 #pragma omp parallel for private(i, j, k, tmpO, tmpU, tmpD) 548 for (k=2; k<=*
nz-1; k++) {
549 for (j=2; j<=*
ny-1; j++) {
550 for(i=2; i<=*
nx-1; i++) {
553 + VAT3( oN, i, j, k) * VAT3(x, i, j+1, k)
554 + VAT3( oN, i, j-1, k) * VAT3(x, i, j-1, k)
555 + VAT3( oE, i, j, k) * VAT3(x, i+1, j, k)
556 + VAT3( oE, i-1, j, k) * VAT3(x, i-1, j, k)
557 + VAT3( oNE, i, j, k) * VAT3(x, i+1, j+1, k)
558 + VAT3( oNW, i, j, k) * VAT3(x, i-1, j+1, k)
559 + VAT3( oNW, i+1, j-1, k) * VAT3(x, i+1, j-1, k)
560 + VAT3( oNE, i-1, j-1, k) * VAT3(x, i-1, j-1, k);
563 + VAT3( uC, i, j, k) * VAT3(x, i, j, k+1)
564 + VAT3( uN, i, j, k) * VAT3(x, i, j+1, k+1)
565 + VAT3( uS, i, j, k) * VAT3(x, i, j-1, k+1)
566 + VAT3( uE, i, j, k) * VAT3(x, i+1, j, k+1)
567 + VAT3( uW, i, j, k) * VAT3(x, i-1, j, k+1)
568 + VAT3( uNE, i, j, k) * VAT3(x, i+1, j+1, k+1)
569 + VAT3( uNW, i, j, k) * VAT3(x, i-1, j+1, k+1)
570 + VAT3( uSE, i, j, k) * VAT3(x, i+1, j-1, k+1)
571 + VAT3( uSW, i, j, k) * VAT3(x, i-1, j-1, k+1);
574 + VAT3( uC, i, j, k-1) * VAT3(x, i, j, k-1)
575 + VAT3( uS, i, j+1, k-1) * VAT3(x, i, j+1, k-1)
576 + VAT3( uN, i, j-1, k-1) * VAT3(x, i, j-1, k-1)
577 + VAT3( uW, i+1, j, k-1) * VAT3(x, i+1, j, k-1)
578 + VAT3( uE, i-1, j, k-1) * VAT3(x, i-1, j, k-1)
579 + VAT3( uSW, i+1, j+1, k-1) * VAT3(x, i+1, j+1, k-1)
580 + VAT3( uSE, i-1, j+1, k-1) * VAT3(x, i-1, j+1, k-1)
581 + VAT3( uNW, i+1, j-1, k-1) * VAT3(x, i+1, j-1, k-1)
582 + VAT3( uNE, i-1, j-1, k-1) * VAT3(x, i-1, j-1, k-1);
584 VAT3(r, i, j, k) = VAT3(fc, i, j, k) + tmpO + tmpU + tmpD
585 - (VAT3(oC, i, j, k) + VAT3(cc, i, j, k)) * VAT3(x, i, j, k);
594 int *ipc,
double *rpc,
595 double *ac,
double *cc,
double *fc,
596 double *x,
double *r,
double *w1) {
601 numdia = VAT(ipc, 11);
603 Vnmresid7(
nx,
ny,
nz, ipc, rpc, ac, cc, fc, x, r, w1);
604 }
else if (numdia == 27) {
605 Vnmresid27(
nx,
ny,
nz, ipc, rpc, ac, cc, fc, x, r, w1);
607 Vnm_print(2,
"Vnmresid: invalid stencil type given...\n");
613 VPUBLIC
void Vnmresid7(
int *
nx,
int *
ny,
int *
nz,
614 int *ipc,
double *rpc,
615 double *ac,
double *cc,
double *fc,
616 double *x,
double *r,
double *w1) {
618 MAT2(ac, *
nx * *
ny * *
nz, 1);
623 RAT2(ac, 1, 1), cc, fc,
624 RAT2(ac, 1, 2), RAT2(ac, 1, 3), RAT2(ac, 1, 4),
628 VPUBLIC
void Vnmresid7_1s(
int *
nx,
int *
ny,
int *
nz,
629 int *ipc,
double *rpc,
630 double *oC,
double *cc,
double *fc,
631 double *oE,
double *oN,
double *uC,
632 double *x,
double *r,
double *w1) {
648 ipkey = VAT(ipc, 10);
652 for (k=2; k<=*
nz-1; k++) {
653 for (j=2; j<=*
ny-1; j++) {
654 for (i=2; i<=*
nx-1; i++) {
655 VAT3(r, i, j, k) = VAT3(fc, i, j, k)
656 + VAT3(oN, i, j, k) * VAT3(x, i, j+1, k)
657 + VAT3(oN, i, j-1, k) * VAT3(x, i, j-1, k)
658 + VAT3(oE, i, j, k) * VAT3(x, i+1, j, k)
659 + VAT3(oE, i-1, j, k) * VAT3(x, i-1, j, k)
660 + VAT3(uC, i, j, k-1) * VAT3(x, i, j, k-1)
661 + VAT3(uC, i, j, k) * VAT3(x, i, j, k+1)
662 - VAT3(oC, i, j, k) * VAT3(x, i, j, k)
671 VPUBLIC
void Vnmresid27(
int *
nx,
int *
ny,
int *
nz,
672 int *ipc,
double *rpc,
673 double *ac,
double *cc,
double *fc,
674 double *x,
double *r,
double *w1) {
676 MAT2(ac, *
nx * *
ny * *
nz, 1);
681 RAT2(ac, 1, 1), cc, fc,
682 RAT2(ac, 1, 2), RAT2(ac, 1, 3), RAT2(ac, 1, 4),
683 RAT2(ac, 1, 5), RAT2(ac, 1, 6),
684 RAT2(ac, 1, 7), RAT2(ac, 1, 8), RAT2(ac, 1, 9), RAT2(ac, 1, 10),
685 RAT2(ac, 1, 11), RAT2(ac, 1, 12), RAT2(ac, 1, 13), RAT2(ac, 1, 14),
691 VPUBLIC
void Vnmresid27_1s(
int *
nx,
int *
ny,
int *
nz,
692 int *ipc,
double *rpc,
693 double *oC,
double *cc,
double *fc,
694 double *oE,
double *oN,
double *uC,
695 double *oNE,
double *oNW,
696 double *uE,
double *uW,
double *uN,
double *uS,
697 double *uNE,
double *uNW,
double *uSE,
double *uSW,
698 double *x,
double *r,
double *w1) {
702 double tmpO, tmpU, tmpD;
725 ipkey = VAT(ipc, 10);
729 for (k=2; k<=*
nz-1; k++) {
730 for (j=2; j<=*
ny-1; j++) {
731 for (i=2; i<=*
nx-1; i++) {
734 + VAT3( oN, i, j, k) * VAT3(x, i, j+1, k)
735 + VAT3( oN, i, j-1, k) * VAT3(x, i, j-1, k)
736 + VAT3( oE, i, j, k) * VAT3(x, i+1, j, k)
737 + VAT3( oE, i-1, j, k) * VAT3(x, i-1, j, k)
738 + VAT3(oNE, i, j, k) * VAT3(x, i+1, j+1, k)
739 + VAT3(oNW, i, j, k) * VAT3(x, i-1, j+1, k)
740 + VAT3(oNW, i+1, j-1, k) * VAT3(x, i+1, j-1, k)
741 + VAT3(oNE, i-1, j-1, k) * VAT3(x, i-1, j-1, k);
744 + VAT3( uC, i, j, k) * VAT3(x, i, j, k+1)
745 + VAT3( uN, i, j, k) * VAT3(x, i, j+1, k+1)
746 + VAT3( uS, i, j, k) * VAT3(x, i, j-1, k+1)
747 + VAT3( uE, i, j, k) * VAT3(x, i+1, j, k+1)
748 + VAT3( uW, i, j, k) * VAT3(x, i-1, j, k+1)
749 + VAT3(uNE, i, j, k) * VAT3(x, i+1, j+1, k+1)
750 + VAT3(uNW, i, j, k) * VAT3(x, i-1, j+1, k+1)
751 + VAT3(uSE, i, j, k) * VAT3(x, i+1, j-1, k+1)
752 + VAT3(uSW, i, j, k) * VAT3(x, i-1, j-1, k+1);
755 + VAT3( uC, i, j, k-1) * VAT3(x, i, j, k-1)
756 + VAT3( uS, i, j+1, k-1) * VAT3(x, i, j+1, k-1)
757 + VAT3( uN, i, j-1, k-1) * VAT3(x, i, j-1, k-1)
758 + VAT3( uW, i+1, j, k-1) * VAT3(x, i+1, j, k-1)
759 + VAT3( uE, i-1, j, k-1) * VAT3(x, i-1, j, k-1)
760 + VAT3(uSW, i+1, j+1, k-1) * VAT3(x, i+1, j+1, k-1)
761 + VAT3(uSE, i-1, j+1, k-1) * VAT3(x, i-1, j+1, k-1)
762 + VAT3(uNW, i+1, j-1, k-1) * VAT3(x, i+1, j-1, k-1)
763 + VAT3(uNE, i-1, j-1, k-1) * VAT3(x, i-1, j-1, k-1);
768 - VAT3(oC, i, j, k) * VAT3(x, i, j, k)
777 VPUBLIC
void Vrestrc(
int *nxf,
int *nyf,
int *nzf,
779 double *xin,
double *xout,
double *pc) {
783 Vrestrc2(nxf, nyf, nzf,
786 RAT2(pc, 1, 1), RAT2(pc, 1, 2), RAT2(pc, 1, 3), RAT2(pc, 1, 4), RAT2(pc, 1, 5),
787 RAT2(pc, 1, 6), RAT2(pc, 1, 7), RAT2(pc, 1, 8), RAT2(pc, 1, 9),
788 RAT2(pc, 1,10), RAT2(pc, 1,11), RAT2(pc, 1,12), RAT2(pc, 1,13), RAT2(pc, 1,14),
789 RAT2(pc, 1,15), RAT2(pc, 1,16), RAT2(pc, 1,17), RAT2(pc, 1,18),
790 RAT2(pc, 1,19), RAT2(pc, 1,20), RAT2(pc, 1,21), RAT2(pc, 1,22), RAT2(pc, 1,23),
791 RAT2(pc, 1,24), RAT2(pc, 1,25), RAT2(pc, 1,26), RAT2(pc, 1,27));
796 VEXTERNC
void Vrestrc2(
int *nxf,
int *nyf,
int *nzf,
798 double *xin,
double *xout,
799 double *oPC,
double *oPN,
double *oPS,
double *oPE,
double *oPW,
800 double *oPNE,
double *oPNW,
double *oPSE,
double *oPSW,
801 double *uPC,
double *uPN,
double *uPS,
double *uPE,
double *uPW,
802 double *uPNE,
double *uPNW,
double *uPSE,
double *uPSW,
803 double *dPC,
double *dPN,
double *dPS,
double *dPE,
double *dPW,
804 double *dPNE,
double *dPNW,
double *dPSE,
double *dPSW) {
810 double tmpO, tmpU, tmpD;
813 MAT3(xin, *nxf, *nyf, *nzf);
852 dimfac = VPOW(2.0, idimenshun);
855 #pragma omp parallel for private(k, kk, j, jj, i, ii, tmpO, tmpU, tmpD) 856 for (k=2; k<=*
nzc-1; k++) {
857 kk = (k - 1) * 2 + 1;
859 for (j=2; j<=*
nyc-1; j++) {
860 jj = (j - 1) * 2 + 1;
862 for (i=2; i<=*
nxc-1; i++) {
863 ii = (i - 1) * 2 + 1;
867 + VAT3( oPC, i, j, k) * VAT3(xin, ii, jj, kk)
868 + VAT3( oPN, i, j, k) * VAT3(xin, ii, jj+1, kk)
869 + VAT3( oPS, i, j, k) * VAT3(xin, ii, jj-1, kk)
870 + VAT3( oPE, i, j, k) * VAT3(xin, ii+1, jj, kk)
871 + VAT3( oPW, i, j, k) * VAT3(xin, ii-1, jj, kk)
872 + VAT3(oPNE, i, j, k) * VAT3(xin, ii+1, jj+1, kk)
873 + VAT3(oPNW, i, j, k) * VAT3(xin, ii-1, jj+1, kk)
874 + VAT3(oPSE, i, j, k) * VAT3(xin, ii+1, jj-1, kk)
875 + VAT3(oPSW, i, j, k) * VAT3(xin, ii-1, jj-1, kk);
878 + VAT3( uPC, i, j, k) * VAT3(xin, ii, jj, kk+1)
879 + VAT3( uPN, i, j, k) * VAT3(xin, ii, jj+1, kk+1)
880 + VAT3( uPS, i, j, k) * VAT3(xin, ii, jj-1, kk+1)
881 + VAT3( uPE, i, j, k) * VAT3(xin, ii+1, jj, kk+1)
882 + VAT3( uPW, i, j, k) * VAT3(xin, ii-1, jj, kk+1)
883 + VAT3(uPNE, i, j, k) * VAT3(xin, ii+1, jj+1, kk+1)
884 + VAT3(uPNW, i, j, k) * VAT3(xin, ii-1, jj+1, kk+1)
885 + VAT3(uPSE, i, j, k) * VAT3(xin, ii+1, jj-1, kk+1)
886 + VAT3(uPSW, i, j, k) * VAT3(xin, ii-1, jj-1, kk+1);
889 + VAT3( dPC, i, j, k) * VAT3(xin, ii, jj, kk-1)
890 + VAT3( dPN, i, j, k) * VAT3(xin, ii, jj+1, kk-1)
891 + VAT3( dPS, i, j, k) * VAT3(xin, ii, jj-1, kk-1)
892 + VAT3( dPE, i, j, k) * VAT3(xin, ii+1, jj, kk-1)
893 + VAT3( dPW, i, j, k) * VAT3(xin, ii-1, jj, kk-1)
894 + VAT3(dPNE, i, j, k) * VAT3(xin, ii+1, jj+1, kk-1)
895 + VAT3(dPNW, i, j, k) * VAT3(xin, ii-1, jj+1, kk-1)
896 + VAT3(dPSE, i, j, k) * VAT3(xin, ii+1, jj-1, kk-1)
897 + VAT3(dPSW, i, j, k) * VAT3(xin, ii-1, jj-1, kk-1);
899 VAT3(xout, i, j, k) = tmpO + tmpU + tmpD;
911 int *nxf,
int *nyf,
int *nzf,
912 double *xin,
double *xout,
920 RAT2(pc, 1, 1), RAT2(pc, 1, 2), RAT2(pc, 1, 3), RAT2(pc, 1, 4), RAT2(pc, 1, 5),
921 RAT2(pc, 1, 6), RAT2(pc, 1, 7), RAT2(pc, 1, 8), RAT2(pc, 1, 9),
922 RAT2(pc, 1,10), RAT2(pc, 1,11), RAT2(pc, 1,12), RAT2(pc, 1,13), RAT2(pc, 1,14),
923 RAT2(pc, 1,15), RAT2(pc, 1,16), RAT2(pc, 1,17), RAT2(pc, 1,18),
924 RAT2(pc, 1,19), RAT2(pc, 1,20), RAT2(pc, 1,21), RAT2(pc, 1,22), RAT2(pc, 1,23),
925 RAT2(pc, 1,24), RAT2(pc, 1,25), RAT2(pc, 1,26), RAT2(pc, 1,27));
930 VPUBLIC
void VinterpPMG2(
int *
nxc,
int *
nyc,
int *
nzc,
931 int *nxf,
int *nyf,
int *nzf,
932 double *xin,
double *xout,
933 double *oPC,
double *oPN,
double *oPS,
double *oPE,
double *oPW,
934 double *oPNE,
double *oPNW,
double *oPSE,
double *oPSW,
935 double *uPC,
double *uPN,
double *uPS,
double *uPE,
double *uPW,
936 double *uPNE,
double *uPNW,
double *uPSE,
double *uPSW,
937 double *dPC,
double *dPN,
double *dPS,
double *dPE,
double *dPW,
938 double *dPNE,
double *dPNW,
double *dPSE,
double *dPSW) {
944 MAT3(xout, *nxf, *nyf, *nzf);
987 for (k=1; k<=*nzf-2; k+=2) {
988 kk = (k - 1) / 2 + 1;
990 for (j=1; j<=*nyf-2; j+=2) {
991 jj = (j - 1) / 2 + 1;
993 for (i=1; i<=*nxf-2; i+=2) {
994 ii = (i - 1) / 2 + 1;
1001 VAT3(xout, i, j, k) = VAT3(xin, ii, jj, kk);
1009 VAT3(xout, i+1, j, k) = VAT3(oPE, ii, jj, kk) * VAT3(xin, ii, jj, kk)
1010 + VAT3(oPW, ii+1, jj, kk) * VAT3(xin, ii+1, jj, kk);
1014 VAT3(xout, i, j+1, k) = VAT3(oPN, ii, jj, kk) * VAT3(xin, ii, jj, kk)
1015 + VAT3(oPS, ii, jj+1, kk) * VAT3(xin, ii, jj+1, kk);
1019 VAT3(xout, i, j, k+1) = VAT3(uPC, ii, jj, kk) * VAT3(xin, ii, jj, kk)
1020 + VAT3(dPC, ii, jj, kk+1) * VAT3(xin, ii, jj, kk+1);
1029 VAT3(xout, i+1, j+1, k) = VAT3(oPNE, ii, jj, kk) * VAT3(xin, ii, jj, kk)
1030 + VAT3(oPNW, ii+1, jj, kk) * VAT3(xin, ii+1, jj, kk)
1031 + VAT3(oPSE, ii, jj+1, kk) * VAT3(xin, ii, jj+1, kk)
1032 + VAT3(oPSW, ii+1, jj+1, kk) * VAT3(xin, ii+1, jj+1, kk);
1036 VAT3(xout, i+1, j, k+1) = VAT3(uPE, ii, jj, kk) * VAT3(xin, ii, jj, kk)
1037 + VAT3(uPW, ii+1, jj, kk) * VAT3(xin, ii+1, jj, kk)
1038 + VAT3(dPE, ii, jj, kk+1) * VAT3(xin, ii, jj, kk+1)
1039 + VAT3(dPW, ii+1, jj, kk+1) * VAT3(xin, ii+1, jj, kk+1);
1043 VAT3(xout, i, j+1, k+1) = VAT3(uPN, ii, jj, kk) * VAT3(xin, ii, jj, kk)
1044 + VAT3(uPS, ii, jj+1, kk) * VAT3(xin, ii, jj+1, kk)
1045 + VAT3(dPN, ii, jj,kk+1) * VAT3(xin, ii, jj, kk+1)
1046 + VAT3(dPS, ii, jj+1,kk+1) * VAT3(xin, ii, jj+1, kk+1);
1054 VAT3(xout, i+1,j+1,k+1) =
1055 + VAT3(uPNE, ii, jj, kk) * VAT3(xin, ii, jj, kk)
1056 + VAT3(uPNW, ii+1, jj, kk) * VAT3(xin, ii+1, jj, kk)
1057 + VAT3(uPSE, ii, jj+1, kk) * VAT3(xin, ii, jj+1, kk)
1058 + VAT3(uPSW, ii+1, jj+1, kk) * VAT3(xin, ii+1, jj+1, kk)
1059 + VAT3(dPNE, ii, jj, kk+1) * VAT3(xin, ii, jj, kk+1)
1060 + VAT3(dPNW, ii+1, jj, kk+1) * VAT3(xin, ii+1, jj, kk+1)
1061 + VAT3(dPSE, ii, jj+1, kk+1) * VAT3(xin, ii, jj+1, kk+1)
1062 + VAT3(dPSW, ii+1, jj+1, kk+1) * VAT3(xin, ii+1, jj+1, kk+1);
1073 VPUBLIC
void Vextrac(
int *nxf,
int *nyf,
int *nzf,
1075 double *xin,
double *xout) {
1080 MAT3( xin, *nxf, *nyf, *nzf);
1087 for (k=2; k<=*
nzc-1; k++) {
1088 kk = (k - 1) * 2 + 1;
1090 for (j=2; j<=*
nyc-1; j++) {
1091 jj = (j - 1) * 2 + 1;
1093 for (i=2; i<=*
nxc-1; i++) {
1094 ii = (i - 1) * 2 + 1;
1097 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.