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) {
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;
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) {
293 ipkey = VAT(ipc, 10);
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;
371 ipkey = VAT(ipc, 10);
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)
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) {
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;
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) {
649 ipkey = VAT(ipc, 10);
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;
726 ipkey = VAT(ipc, 10);
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);
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) {
945 MAT3(xout, *nxf, *nyf, *nzf);
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,
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.