57VPUBLIC
void Vmatvec(
int *nx,
int *ny,
int *nz,
58 int *ipc,
double *rpc,
59 double *ac,
double *cc,
60 double *x,
double *y) {
65 numdia = VAT(ipc, 11);
72 }
else if (numdia == 27) {
78 Vnm_print(2,
"MATVEC: invalid stencil type given...");
82VPUBLIC
void Vmatvec7(
int *nx,
int *ny,
int *nz,
83 int *ipc,
double *rpc,
84 double *ac,
double *cc,
85 double *x,
double *y) {
87 MAT2(ac, *nx * *ny * *nz, 1);
89 Vmatvec7_1s(nx, ny, nz,
92 RAT2(ac, 1, 2), RAT2(ac, 1, 3), RAT2(ac, 1, 4),
98VEXTERNC
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) {
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);
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++) {
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);
134VPUBLIC
void Vmatvec27(
int *nx,
int *ny,
int *nz,
135 int *ipc,
double *rpc,
136 double *ac,
double *cc,
137 double *x,
double *y) {
139 MAT2(ac, *nx * *ny * *nz, 1);
141 Vmatvec27_1s(nx, ny, nz,
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),
153VPUBLIC
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) {
164 double tmpO, tmpU, tmpD;
166 MAT3(cc, *nx, *ny, *nz);
167 MAT3(x, *nx, *ny, *nz);
168 MAT3(y, *nx, *ny, *nz);
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);
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);
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++) {
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);
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);
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);
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);
233 int *ipc,
double *rpc,
234 double *ac,
double *cc,
double *x,
double *y,
double *w1) {
239 numdia = VAT(ipc, 11);
242 Vnmatvec7(nx, ny, nz,
246 }
else if (numdia == 27) {
247 Vnmatvec27(nx, ny, nz,
252 Vnm_print(2,
"MATVEC: invalid stencil type given...");
258VPUBLIC
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) {
263 MAT2(ac, *nx * *ny * *nz, 1);
267 Vnmatvecd7_1s(nx, ny, nz,
270 RAT2(ac, 1, 2), RAT2(ac, 1, 3), RAT2(ac, 1, 4),
276VPUBLIC
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) {
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);
297 ipkey = VAT(ipc, 10);
298 Vc_vec(cc, x, w1, nx, ny, nz, &ipkey);
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++)
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)
317VPUBLIC
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) {
322 MAT2(ac, *nx * *ny * *nz, 1);
327 Vnmatvecd27_1s(nx, ny, nz,
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),
339VPUBLIC
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) {
351 double tmpO, tmpU, tmpD;
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);
375 ipkey = VAT(ipc, 10);
376 Vc_vec(cc, x, w1, nx, ny, nz, &ipkey);
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++) {
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);
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);
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);
416 VAT3(y, i, j, k) = tmpO + tmpU + tmpD
417 + VAT3(oC, i, j, k) * VAT3(x, i, j, k)
426VPUBLIC
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) {
434 numdia = VAT(ipc, 11);
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);
440 Vnm_print(2,
"Vmresid: invalid stencil type given...\n");
446VPUBLIC
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) {
451 MAT2(ac, *nx * *ny * *nz, 1);
454 Vmresid7_1s(nx, ny, nz,
456 RAT2(ac, 1,1), cc, fc,
457 RAT2(ac, 1,2), RAT2(ac, 1,3), RAT2(ac, 1,4),
461VPUBLIC
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) {
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);
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);
498VPUBLIC
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) {
503 MAT2(ac, *nx * *ny * *nz, 1);
506 Vmresid27_1s(nx,ny,nz,
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),
518VPUBLIC
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) {
529 double tmpO, tmpU, tmpD;
531 MAT3(cc, *nx, *ny, *nz);
532 MAT3(fc, *nx, *ny, *nz);
533 MAT3(x, *nx, *ny, *nz);
534 MAT3(r, *nx, *ny, *nz);
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);
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);
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++) {
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);
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);
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);
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);
599 int *ipc,
double *rpc,
600 double *ac,
double *cc,
double *fc,
601 double *x,
double *r,
double *w1) {
606 numdia = VAT(ipc, 11);
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);
612 Vnm_print(2,
"Vnmresid: invalid stencil type given...\n");
618VPUBLIC
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) {
623 MAT2(ac, *nx * *ny * *nz, 1);
626 Vnmresid7_1s(nx, ny, nz,
628 RAT2(ac, 1, 1), cc, fc,
629 RAT2(ac, 1, 2), RAT2(ac, 1, 3), RAT2(ac, 1, 4),
633VPUBLIC
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) {
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);
653 ipkey = VAT(ipc, 10);
654 Vc_vec(cc, x, w1, nx, ny, nz, &ipkey);
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)
676VPUBLIC
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) {
681 MAT2(ac, *nx * *ny * *nz, 1);
684 Vnmresid27_1s(nx, ny, nz,
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),
696VPUBLIC
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) {
707 double tmpO, tmpU, tmpD;
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);
730 ipkey = VAT(ipc, 10);
731 Vc_vec(cc, x, w1, nx, ny, nz, &ipkey);
734 for (k=2; k<=*nz-1; k++) {
735 for (j=2; j<=*ny-1; j++) {
736 for (i=2; i<=*nx-1; i++) {
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);
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);
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);
773 - VAT3(oC, i, j, k) * VAT3(x, i, j, k)
782VPUBLIC
void Vrestrc(
int *nxf,
int *nyf,
int *nzf,
783 int *nxc,
int *nyc,
int *nzc,
784 double *xin,
double *xout,
double *pc) {
786 MAT2(pc, *nxc * *nyc * *nzc, 1 );
788 Vrestrc2(nxf, nyf, nzf,
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));
801VEXTERNC
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) {
815 double tmpO, tmpU, tmpD;
818 MAT3(xin, *nxf, *nyf, *nzf);
819 MAT3(xout, *nxc, *nyc, *nzc);
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);
827 MAT3(oPNE, *nxc, *nyc, *nzc);
828 MAT3(oPNW, *nxc, *nyc, *nzc);
829 MAT3(oPSE, *nxc, *nyc, *nzc);
830 MAT3(oPSW, *nxc, *nyc, *nzc);
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);
838 MAT3(uPNE, *nxc, *nyc, *nzc);
839 MAT3(uPNW, *nxc, *nyc, *nzc);
840 MAT3(uPSE, *nxc, *nyc, *nzc);
841 MAT3(uPSW, *nxc, *nyc, *nzc);
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);
849 MAT3(dPNE, *nxc, *nyc, *nzc);
850 MAT3(dPNW, *nxc, *nyc, *nzc);
851 MAT3(dPSE, *nxc, *nyc, *nzc);
852 MAT3(dPSW, *nxc, *nyc, *nzc);
857 dimfac = VPOW(2.0, idimenshun);
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;
864 for (j=2; j<=*nyc-1; j++) {
865 jj = (j - 1) * 2 + 1;
867 for (i=2; i<=*nxc-1; i++) {
868 ii = (i - 1) * 2 + 1;
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);
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);
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);
904 VAT3(xout, i, j, k) = tmpO + tmpU + tmpD;
916 int *nxf,
int *nyf,
int *nzf,
917 double *xin,
double *xout,
920 MAT2(pc, *nxc * *nyc * *nzc, 1);
922 VinterpPMG2(nxc, nyc, nzc,
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));
935VPUBLIC
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) {
948 MAT3( xin, *nxc, *nyc, *nzc);
949 MAT3(xout, *nxf, *nyf, *nzf);
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);
957 MAT3(oPNE, *nxc, *nyc, *nzc);
958 MAT3(oPNW, *nxc, *nyc, *nzc);
959 MAT3(oPSE, *nxc, *nyc, *nzc);
960 MAT3(oPSW, *nxc, *nyc, *nzc);
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);
968 MAT3(uPNE, *nxc, *nyc, *nzc);
969 MAT3(uPNW, *nxc, *nyc, *nzc);
970 MAT3(uPSE, *nxc, *nyc, *nzc);
971 MAT3(uPSW, *nxc, *nyc, *nzc);
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);
979 MAT3(dPNE, *nxc, *nyc, *nzc);
980 MAT3(dPNW, *nxc, *nyc, *nzc);
981 MAT3(dPSE, *nxc, *nyc, *nzc);
982 MAT3(dPSW, *nxc, *nyc, *nzc);
992 for (k=1; k<=*nzf-2; k+=2) {
993 kk = (k - 1) / 2 + 1;
995 for (j=1; j<=*nyf-2; j+=2) {
996 jj = (j - 1) / 2 + 1;
998 for (i=1; i<=*nxf-2; i+=2) {
999 ii = (i - 1) / 2 + 1;
1006 VAT3(xout, i, j, k) = VAT3(xin, ii, jj, kk);
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);
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);
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);
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);
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);
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);
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);
1079 int *nxc,
int *nyc,
int *nzc,
1080 double *xin,
double *xout) {
1085 MAT3( xin, *nxf, *nyf, *nzf);
1086 MAT3(xout, *nxc, *nyc, *nzc);
1092 for (k=2; k<=*nzc-1; k++) {
1093 kk = (k - 1) * 2 + 1;
1095 for (j=2; j<=*nyc-1; j++) {
1096 jj = (j - 1) * 2 + 1;
1098 for (i=2; i<=*nxc-1; i++) {
1099 ii = (i - 1) * 2 + 1;
1102 VAT3(xout, i, j, k) = VAT3(xin, ii, jj, kk);
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 Vc_vec(double *coef, double *uin, double *uout, int *nx, int *ny, int *nz, int *ipkey)
Define the nonlinearity (vector version)
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 VfboundPMG00(int *nx, int *ny, int *nz, double *x)
Initialize a grid function to have a zero boundary value.
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 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 Vrestrc(int *nxf, int *nyf, int *nzf, int *nxc, int *nyc, int *nzc, double *xin, double *xout, double *pc)
Apply the restriction operator.
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.
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.