57 double *rpc,
double *pc,
double *ac,
double *cc,
double *fc,
58 double *xf,
double *yf,
double *zf,
59 double *gxcf,
double *gycf,
double *gzcf,
60 double *a1cf,
double *a2cf,
double *a3cf,
61 double *ccf,
double *fcf,
double *tcf
86 if (*ido == 0 || *ido == 2) {
92 VMESSAGE3(
"Fine: (%03d, %03d, %03d)", nxx, nyy, nzz);
97 RAT(ipc, VAT2(iz, 5,lev)), RAT(rpc, VAT2(iz, 6,lev)),
98 RAT(ac, VAT2(iz, 7,lev)), RAT(cc, VAT2(iz, 1,lev)), RAT(fc, VAT2(iz, 1,lev)),
99 RAT(xf, VAT2(iz, 8,lev)), RAT(yf, VAT2(iz, 9,lev)), RAT(zf, VAT2(iz, 10,lev)),
100 RAT(gxcf, VAT2(iz, 2,lev)), RAT(gycf, VAT2(iz, 3,lev)), RAT(gzcf, VAT2(iz, 4,lev)),
101 RAT(a1cf, VAT2(iz, 1,lev)), RAT(a2cf, VAT2(iz, 1,lev)), RAT(a3cf, VAT2(iz, 1,lev)),
102 RAT(ccf, VAT2(iz, 1,lev)), RAT(fcf, VAT2(iz, 1,lev)));
104 VMESSAGE2(
"Operator stencil (lev, numdia) = (%d, %d)", lev, numdia);
107 VAT2(iz, 7, lev+1) = VAT2(iz, 7, lev) + numdia * nxx * nyy * nzz;
112 RAT(ipc, VAT2(iz, 5,lev)), RAT(rpc, VAT2(iz, 6,lev)), RAT(ac, VAT2(iz, 7,lev)));
117 if (*ido == 1 || *ido == 2 || *ido == 3) {
119 for (lev=2; lev<=*
nlev; lev++) {
126 Vmkcors(&i, &nxold, &nyold, &nzold, &nxx, &nyy, &nzz);
130 VbuildP(&nxold, &nyold, &nzold,
133 RAT(ipc, VAT2(iz, 5,lev-1)), RAT(rpc, VAT2(iz, 6,lev-1)),
134 RAT(pc, VAT2(iz, 11,lev-1)), RAT(ac, VAT2(iz, 7,lev-1)),
135 RAT(xf, VAT2(iz, 8,lev-1)), RAT(yf, VAT2(iz, 9,lev-1)), RAT(zf, VAT2(iz, 10,lev-1)));
142 VMESSAGE3(
"Stand: (%03d, %03d, %03d)", nxx, nyy, nzz);
146 Vbuildcopy0(&nxx, &nyy, &nzz,
147 &nxold, &nyold, &nzold,
148 RAT(xf, VAT2(iz, 8,lev )), RAT(yf, VAT2(iz, 9,lev )), RAT(zf, VAT2(iz, 10,lev )),
149 RAT(gxcf, VAT2(iz, 2,lev )), RAT(gycf, VAT2(iz, 3,lev )), RAT(gzcf, VAT2(iz, 4,lev )),
150 RAT(a1cf, VAT2(iz, 1,lev )), RAT(a2cf, VAT2(iz, 1,lev )), RAT(a3cf, VAT2(iz, 1,lev )),
151 RAT(ccf, VAT2(iz, 1,lev )), RAT(fcf, VAT2(iz, 1,lev )), RAT(tcf, VAT2(iz, 1,lev )),
152 RAT(xf, VAT2(iz, 8,lev-1)), RAT(yf, VAT2(iz, 9,lev-1)), RAT(zf, VAT2(iz, 10,lev-1)),
153 RAT(gxcf, VAT2(iz, 2,lev-1)), RAT(gycf, VAT2(iz, 3,lev-1)), RAT(gzcf, VAT2(iz, 4,lev-1)),
154 RAT(a1cf, VAT2(iz, 1,lev-1)), RAT(a2cf, VAT2(iz, 1,lev-1)), RAT(a3cf, VAT2(iz, 1,lev-1)),
155 RAT(ccf, VAT2(iz, 1,lev-1)), RAT(fcf, VAT2(iz, 1,lev-1)), RAT(tcf, VAT2(iz, 1,lev-1)));
159 RAT(ipc, VAT2(iz, 5,lev)), RAT(rpc, VAT2(iz, 6,lev)),
160 RAT(ac, VAT2(iz, 7,lev)), RAT(cc, VAT2(iz, 1,lev)), RAT(fc, VAT2(iz, 1,lev)),
161 RAT(xf, VAT2(iz, 8,lev)), RAT(yf, VAT2(iz, 9,lev)), RAT(zf, VAT2(iz, 10,lev)),
162 RAT(gxcf, VAT2(iz, 2,lev)), RAT(gycf, VAT2(iz, 3,lev)), RAT(gzcf, VAT2(iz, 4,lev)),
163 RAT(a1cf, VAT2(iz, 1,lev)), RAT(a2cf, VAT2(iz, 1,lev)), RAT(a3cf, VAT2(iz, 1,lev)),
164 RAT(ccf, VAT2(iz, 1,lev)), RAT(fcf, VAT2(iz, 1,lev)));
172 VMESSAGE3(
"Harm0: (%03d, %03d, %03d)", nxx, nyy, nzz);
174 Vbuildharm0(&nxx, &nyy, &nzz, &nxold, &nyold, &nzold,
175 RAT(xf, VAT2(iz, 8, lev )), RAT(yf, VAT2(iz, 9, lev )), RAT(zf, VAT2(iz, 10, lev )),
176 RAT(gxcf, VAT2(iz, 2, lev )), RAT(gycf, VAT2(iz, 3, lev )), RAT(gzcf, VAT2(iz, 4, lev )),
177 RAT(a1cf, VAT2(iz, 1, lev )), RAT(a2cf, VAT2(iz, 1, lev )), RAT(a3cf, VAT2(iz, 1, lev )),
178 RAT(ccf, VAT2(iz, 1, lev )), RAT(fcf, VAT2(iz, 1, lev )), RAT(tcf, VAT2(iz, 1, lev )),
179 RAT(xf, VAT2(iz, 8, lev-1)), RAT(yf, VAT2(iz, 9, lev-1)), RAT(zf, VAT2(iz, 10, lev-1)),
180 RAT(gxcf, VAT2(iz, 2, lev-1)), RAT(gycf, VAT2(iz, 3, lev-1)), RAT(gzcf, VAT2(iz, 4, lev-1)),
181 RAT(a1cf, VAT2(iz, 1, lev-1)), RAT(a2cf, VAT2(iz, 1, lev-1)), RAT(a3cf, VAT2(iz, 1, lev-1)),
182 RAT(ccf, VAT2(iz, 1, lev-1)), RAT(fcf, VAT2(iz, 1, lev-1)), RAT(tcf, VAT2(iz, 1, lev-1)));
186 RAT(ipc, VAT2(iz, 5,lev)), RAT(rpc, VAT2(iz, 6,lev)),
187 RAT(ac, VAT2(iz, 7,lev)), RAT(cc, VAT2(iz, 1,lev)), RAT(fc, VAT2(iz, 1,lev)),
188 RAT(xf, VAT2(iz, 8,lev)), RAT(yf, VAT2(iz, 9,lev)), RAT(zf, VAT2(iz, 10,lev)),
189 RAT(gxcf, VAT2(iz, 2,lev)), RAT(gycf, VAT2(iz, 3,lev)), RAT(gzcf, VAT2(iz, 4,lev)),
190 RAT(a1cf, VAT2(iz, 1,lev)), RAT(a2cf, VAT2(iz, 1,lev)), RAT(a3cf, VAT2(iz, 1,lev)),
191 RAT(ccf, VAT2(iz, 1,lev)), RAT(fcf, VAT2(iz, 1,lev)));
199 VMESSAGE3(
"Galer: (%03d, %03d, %03d)", nxx, nyy, nzz);
204 RAT(pc, VAT2(iz, 11,lev-1)),
205 RAT(ipc, VAT2(iz, 5,lev-1)), RAT(rpc, VAT2(iz, 6,lev-1)),
206 RAT(ac, VAT2(iz, 7,lev-1)), RAT(cc, VAT2(iz, 1,lev-1)), RAT(fc, VAT2(iz, 1,lev-1)),
207 RAT(ipc, VAT2(iz, 5,lev )), RAT(rpc, VAT2(iz, 6,lev )),
208 RAT(ac, VAT2(iz, 7,lev )), RAT(cc, VAT2(iz, 1,lev )), RAT(fc, VAT2(iz, 1,lev )));
212 Vextrac(&nxold, &nyold, &nzold,
214 RAT(tcf, VAT2(iz, 1,lev-1)), RAT(tcf, VAT2(iz, 1,lev)));
217 VABORT_MSG1(
"Bad mgcoar value given: %d", *
mgcoar);
223 VAT2(iz, 7, lev+1) = VAT2(iz, 7,lev) + numdia * nxx * nyy * nzz;
229 RAT(ipc, VAT2(iz, 5,lev)), RAT(rpc, VAT2(iz, 6,lev)), RAT(ac, VAT2(iz, 7,lev)));
239 RAT(ipc, VAT2(iz, 5,lev )), RAT(rpc, VAT2(iz, 6,lev )), RAT(ac, VAT2(iz, 7,lev )),
240 RAT(ipc, VAT2(iz, 5,lev+1)), RAT(rpc, VAT2(iz, 6,lev+1)), RAT(ac, VAT2(iz, 7,lev+1)));
243 VERRMSG0(
"Changing your mgsolv to iterative");
254 int nxold, nyold, nzold;
255 int nxnew, nynew, nznew;
267 n = nxnew * nynew * nznew;
273 VAT2(iz, 1, lev) = 1;
274 VAT2(iz, 2, lev) = 1;
275 VAT2(iz, 3, lev) = 1;
276 VAT2(iz, 4, lev) = 1;
277 VAT2(iz, 5, lev) = 1;
278 VAT2(iz, 6, lev) = 1;
279 VAT2(iz, 7, lev) = 1;
280 VAT2(iz, 8, lev) = 1;
281 VAT2(iz, 9, lev) = 1;
282 VAT2(iz, 10, lev) = 1;
283 VAT2(iz, 11, lev) = 1;
286 VAT2(iz, 1, lev + 1) = VAT2(iz, 1, lev) + n;
287 VAT2(iz, 2, lev + 1) = VAT2(iz, 2, lev) + 4 * nynew * nznew;
288 VAT2(iz, 3, lev + 1) = VAT2(iz, 3, lev) + 4 * nxnew * nznew;
289 VAT2(iz, 4, lev + 1) = VAT2(iz, 4, lev) + 4 * nxnew * nynew;
290 VAT2(iz, 5, lev + 1) = VAT2(iz, 5, lev) + 100;
291 VAT2(iz, 6, lev + 1) = VAT2(iz, 6, lev) + 100;
292 VAT2(iz, 8, lev + 1) = VAT2(iz, 8, lev) + nxnew;
293 VAT2(iz, 9, lev + 1) = VAT2(iz, 9, lev) + nynew;
294 VAT2(iz, 10, lev + 1) = VAT2(iz, 10, lev) + nznew;
305 for (lev=2; lev<=*
nlev; lev++) {
310 Vmkcors(&numlev, &nxold, &nyold, &nzold, &nxnew, &nynew, &nznew);
311 n = nxnew * nynew * nznew;
314 VAT2(iz, 1, lev + 1) = VAT2(iz, 1, lev) + n;
315 VAT2(iz, 2, lev + 1) = VAT2(iz, 2, lev) + 4 * nynew * nznew;
316 VAT2(iz, 3, lev + 1) = VAT2(iz, 3, lev) + 4 * nxnew * nznew;
317 VAT2(iz, 4, lev + 1) = VAT2(iz, 4, lev) + 4 * nxnew * nynew;
318 VAT2(iz, 5, lev + 1) = VAT2(iz, 5, lev) + 100;
319 VAT2(iz, 6, lev + 1) = VAT2(iz, 6, lev) + 100;
320 VAT2(iz, 7, lev + 1) = VAT2(iz, 7, lev) + 4 * n;
321 VAT2(iz, 8, lev + 1) = VAT2(iz, 8, lev) + nxnew;
322 VAT2(iz, 9, lev + 1) = VAT2(iz, 9, lev) + nynew;
323 VAT2(iz, 10, lev + 1) = VAT2(iz, 10, lev) + nznew;
326 VAT2(iz, 11, lev) = VAT2(iz, 11, lev - 1) + 27 * n;
338 int *
ipkey,
int *numdia,
339 double *pcFF,
int *ipcFF,
double *rpcFF,
340 double *acFF,
double *ccFF,
double *fcFF,
341 int *ipc,
double *rpc,
342 double *ac,
double *cc,
double *fc) {
347 numdia_loc = VAT(ipcFF, 11);
358 VAT(ipc, 10) = *
ipkey;
372 VPUBLIC
void Vmkcors(
int *numlev,
373 int *nxold,
int *nyold,
int *nzold,
374 int *nxnew,
int *nynew,
int *nznew) {
375 int nxtmp, nytmp, nztmp;
383 for (i=1; i<=*numlev; i++) {
387 Vcorsr(&nxtmp,nxnew);
388 Vcorsr(&nytmp,nynew);
389 Vcorsr(&nztmp,nznew);
395 VPUBLIC
void Vcorsr(
int *nold,
int *nnew) {
398 *nnew = (*nold - 1) / 2 + 1;
401 if ((*nnew - 1) * 2 != *nold - 1) {
402 Vnm_print(2,
"Vcorsr: WARNING! The grid dimensions are not\n");
403 Vnm_print(2,
"Vcorsr: consistent with nlev! Your\n");
404 Vnm_print(2,
"Vcorsr: calculation will only work if you\n");
405 Vnm_print(2,
"Vcorsr: are performing a mg-dummy run.\n");
409 Vnm_print(2,
"Vcorsr: ERROR! The grid dimensions are not\n");
410 Vnm_print(2,
"Vcorsr: consistent with nlev!\n");
411 Vnm_print(2,
"Vcorsr: Grid coarsened below zero.\n");
417 VPUBLIC
void Vmkfine(
int *numlev,
418 int *nxold,
int *nyold,
int *nzold,
419 int *nxnew,
int *nynew,
int *nznew) {
421 int nxtmp, nytmp, nztmp, i;
428 for (i=1; i<=*numlev; i++) {
434 Vfiner(&nxtmp, nxnew);
435 Vfiner(&nytmp, nynew);
436 Vfiner(&nztmp, nznew);
443 VPUBLIC
void Vfiner(
int *nold,
int *nnew) {
445 *nnew = (*nold - 1) * 2 + 1;
450 VPUBLIC
int Vivariv(
int *nu,
int *level) {
466 VPUBLIC
int Vmaxlev(
int n1,
int n2,
int n3) {
480 iden = (int)VPOW(2, lev - 1);
481 n1c = (n1 - 1) / iden + 1;
482 n2c = (n2 - 1) / iden + 1;
483 n3c = (n3 - 1) / iden + 1;
484 if (((n1c - 1) * iden != (n1 - 1)) || (n1c <= 2) )
486 if (((n2c - 1) * iden != (n2 - 1)) || (n2c <= 2) )
488 if (((n3c - 1) * iden != (n3 - 1)) || (n3c <= 2) )
491 }
while (idone != 1);
498 VPUBLIC
void Vprtstp(
int iok,
int iters,
499 double rsnrm,
double rsden,
double orsnrm) {
502 double contrac = 0.0;
512 else if (iters == -1) {
513 Vnm_tstop(40,
"MG iteration");
526 VERRMSG0(
"Vprtstp: avoided division by zero\n");
528 relres = rsnrm / rsden;
534 VERRMSG0(
"avoided division by zero\n");
536 contrac = rsnrm / orsnrm;
540 if (iok == 1 || iok == 2) {
541 VMESSAGE1(
"iteration = %d", iters);
542 VMESSAGE1(
"relative residual = %e", relres);
543 VMESSAGE1(
"contraction number = %e", contrac);
559 VAT(iparm, 1) = *
nrwk;
560 VAT(iparm, 2) = *
niwk;
564 VAT(iparm, 6) = *
nlev;
565 VAT(iparm, 7) = *
nu1;
566 VAT(iparm, 8) = *
nu2;
567 VAT(iparm, 9) = *
mgkey;
568 VAT(iparm, 10) = *
itmax;
569 VAT(iparm, 11) = *
istop;
570 VAT(iparm, 12) = *
iinfo;
571 VAT(iparm, 13) = *
irite;
572 VAT(iparm, 14) = *
ipkey;
573 VAT(iparm, 15) = *
ipcon;
580 VAT(iparm, 22) = *
iperf;
590 VEXTERNC
void Vbuildcopy0(
int *
nx,
int *
ny,
int *
nz,
591 int *nxf,
int *nyf,
int *nzf,
592 double *xc,
double *yc,
double *zc,
593 double *gxc,
double *gyc,
double *gzc,
594 double *a1c,
double *a2c,
double *a3c,
595 double *cc,
double *fc,
double *tc,
596 double *xf,
double *yf,
double *zf,
597 double *gxcf,
double *gycf,
double *gzcf,
598 double *a1cf,
double *a2cf,
double *a3cf,
599 double *ccf,
double *fcf,
double *tcf) {
604 int iadd, jadd, kadd;
606 MAT3( gxc, *
ny, *
nz, 4);
607 MAT3( gyc, *
nx, *
nz, 4);
608 MAT3( gzc, *
nx, *
ny, 4);
609 MAT3( a1c, *
nx, *
ny, *
nz);
610 MAT3( a2c, *
nx, *
ny, *
nz);
611 MAT3( a3c, *
nx, *
ny, *
nz);
615 MAT3(gxcf, *nyf, *nzf, 4);
616 MAT3(gycf, *nxf, *nzf, 4);
617 MAT3(gzcf, *nxf, *nyf, 4);
618 MAT3(a1cf, *nxf, *nyf, *nzf);
619 MAT3(a2cf, *nxf, *nyf, *nzf);
620 MAT3(a3cf, *nxf, *nyf, *nzf);
621 MAT3( tcf, *nxf, *nyf, *nzf);
622 MAT3( ccf, *nxf, *nyf, *nzf);
623 MAT3( fcf, *nxf, *nyf, *nzf);
628 iadd = (*nxf - 1) / (*
nx - 1);
629 jadd = (*nyf - 1) / (*
ny - 1);
630 kadd = (*nzf - 1) / (*
nz - 1);
632 if (iadd != 2 || jadd != 2 || kadd != 2) {
633 Vnm_print(2,
"Vbuildcopy0: Problem with grid dimensions...\n");
637 for (k=1; k<=*
nz; k++) {
639 VAT(zc, k) = VAT(zf, kk);
641 for (j=1; j<=*
ny; j++) {
643 VAT(yc, j) = VAT(yf, jj);
645 for (i=1; i<=*
nx; i++){
647 VAT(xc, i) = VAT(xf, ii);
650 VAT3( tc, i, j, k) = VAT3( tcf, ii, jj, kk);
653 VAT3( cc, i, j, k) = VAT3( ccf, ii, jj, kk);
656 VAT3( fc, i, j, k) = VAT3( fcf, ii, jj, kk);
659 VAT3(a1c, i, j, k) = VAT3(a1cf, ii, jj, kk);
662 VAT3(a2c, i, j, k) = VAT3(a2cf, ii, jj, kk);
665 VAT3(a3c, i, j, k) = VAT3(a3cf, ii, jj, kk);
671 for (k=1; k<=*
nz; k++) {
674 for (j=1; j<=*
ny; j++) {
677 VAT3(gxc, j, k, 1) = VAT3(gxcf, jj, kk, 1);
678 VAT3(gxc, j, k, 2) = VAT3(gxcf, jj, kk, 2);
679 VAT3(gxc, j, k, 3) = VAT3(gxcf, jj, kk, 3);
680 VAT3(gxc, j, k, 4) = VAT3(gxcf, jj, kk, 4);
685 for (k=1; k<=*
nz; k++) {
688 for (i=1; i<=*
nx; i++) {
691 VAT3(gyc, i, k, 1) = VAT3(gycf, ii, kk, 1);
692 VAT3(gyc, i, k, 2) = VAT3(gycf, ii, kk, 2);
693 VAT3(gyc, i, k, 3) = VAT3(gycf, ii, kk, 3);
694 VAT3(gyc, i, k, 4) = VAT3(gycf, ii, kk, 4);
699 for (j=1; j<=*
ny; j++) {
702 for (i=1; i<=*
nx; i++) {
705 VAT3(gzc, i, j, 1) = VAT3(gzcf, ii, jj, 1);
706 VAT3(gzc, i, j, 2) = VAT3(gzcf, ii, jj, 2);
707 VAT3(gzc, i, j, 3) = VAT3(gzcf, ii, jj, 3);
708 VAT3(gzc, i, j, 4) = VAT3(gzcf, ii, jj, 4);
713 VPUBLIC
void Vbuildharm0(
int *
nx,
int *
ny,
int *
nz,
714 int *nxf,
int *nyf,
int *nzf,
715 double *xc,
double *yc,
double *zc,
716 double *gxc,
double *gyc,
double *gzc,
717 double *a1c,
double *a2c,
double *a3c,
718 double *cc,
double *fc,
double *tc,
719 double *xf,
double *yf,
double *zf,
720 double *gxcf,
double *gycf,
double *gzcf,
721 double *a1cf,
double *a2cf,
double *a3cf,
722 double *ccf,
double *fcf,
double *tcf) {
724 Vnm_print(2,
"WARNING: FUNCTION IS NOT FULLY IMPLEMENTED YET!!!");
728 int iadd, jadd, kadd;
730 MAT3( gxc, *
ny, *
nz, 4);
731 MAT3( gyc, *
nx, *
nz, 4);
732 MAT3( gzc, *
nx, *
ny, 4);
734 MAT3( a1c, *
nx, *
ny, *
nz);
735 MAT3( a2c, *
nx, *
ny, *
nz);
736 MAT3( a3c, *
nx, *
ny, *
nz);
742 MAT3(gxcf, *nyf, *nzf, 4);
743 MAT3(gycf, *nxf, *nzf, 4);
744 MAT3(gzcf, *nxf, *nyf, 4);
745 MAT3(a1cf, *nxf, *nyf, *nzf);
746 MAT3(a2cf, *nxf, *nyf, *nzf);
747 MAT3(a3cf, *nxf, *nyf, *nzf);
748 MAT3( tcf, *nxf, *nyf, *nzf);
749 MAT3( ccf, *nxf, *nyf, *nzf);
750 MAT3( fcf, *nxf, *nyf, *nzf);
754 double a, b, c, d, e, f, g, h;
757 iadd = (*nxf - 1) / (*
nx - 1);
758 jadd = (*nyf - 1) / (*
ny - 1);
759 kadd = (*nzf - 1) / (*
nz - 1);
760 if (iadd !=2 || jadd != 2 || kadd != 2) {
761 Vnm_print(2,
"BUILDHARM0: problem with grid dimensions...\n");
765 for (k=1; k<=*
nz; k++) {
767 VAT(zc, k) = VAT(zf, kk);
769 for (j=1; j<=*
ny; j++) {
771 VAT(yc, j) = VAT(yf,jj);
773 for (i=1; i<=*
nx; i++) {
775 VAT(xc, i) = VAT(xf, ii);
778 VAT3(tc, i, j, k) = VAT3(tcf, ii, jj, kk);
781 VAT3(cc, i, j, k) = VAT3(ccf, ii, jj, kk);
796 VAT3(fc, i, j, k) = VAT3(fcf, ii, jj, kk);
810 VAT3(a1c, i, j, k) = (
811 +0.500 *
HARMO2(VAT3(a1cf, ii, jj, kk),
812 VAT3(a1cf, VMIN2(*nxf, ii+1), jj, kk))
813 +0.125 *
HARMO2(VAT3(a1cf, ii, jj, VMAX2(1, kk-1)),
814 VAT3(a1cf, VMIN2(*nxf, ii+1), jj, VMAX2(1, kk-1)))
815 +0.125 *
HARMO2(VAT3(a1cf, ii, jj, VMIN2(*nzf, kk+1)),
816 VAT3(a1cf, VMIN2(*nxf, ii+1), jj, VMIN2(*nzf, kk+1)))
817 +0.125 *
HARMO2(VAT3(a1cf, ii, VMAX2(1, jj-1), kk),
818 VAT3(a1cf, VMIN2(*nxf, ii+1), VMAX2(1, jj-1), kk))
819 +0.125 *
HARMO2(VAT3(a1cf, ii, VMIN2(*nyf, jj+1), kk),
820 VAT3(a1cf, VMIN2(*nxf, ii+1), VMIN2(*nyf, jj+1), kk))
824 VAT3(a2c, i, j, k) = (
825 +0.500 *
HARMO2(VAT3(a2cf, ii, jj, kk),
826 VAT3(a2cf, ii, VMIN2(*nyf, jj+1), kk))
827 +0.125 *
HARMO2(VAT3(a2cf, ii, jj, VMAX2(1, kk-1)),
828 VAT3(a2cf, ii, VMIN2(*nyf, jj+1), VMAX2(1, kk-1)))
829 +0.125 *
HARMO2(VAT3(a2cf, ii, jj, VMIN2(*nzf, kk+1)),
830 VAT3(a2cf, ii, VMIN2(*nyf, jj+1), VMIN2(*nzf, kk+1)))
831 +0.125 *
HARMO2(VAT3(a2cf, VMAX2(1, ii-1), jj, kk),
832 VAT3(a2cf, VMAX2(1, ii-1), VMIN2(nyf, jj+1), kk))
833 +0.125 *
HARMO2(VAT3(a2cf, VMIN2(*nxf, ii+1), jj, kk),
834 VAT3(a2cf, VMIN2(*nxf, ii+1), VMIN2(*nyf, jj+1), kk))
838 VAT3(a3c, i, j, k) = (
839 +0.500 *
HARMO2(VAT3(a3cf, ii, jj, kk),
840 VAT3(a3cf, ii, jj, VMIN2(*nzf, kk+1)))
841 +0.125 *
HARMO2(VAT3(a3cf, ii, VMAX2(1, jj-1), kk),
842 VAT3(a3cf, ii, VMAX2(1, jj-1), VMIN2(*nzf, kk+1)))
843 +0.125 *
HARMO2(VAT3(a3cf, ii, VMIN2(*nyf, jj+1), kk),
844 VAT3(a3cf, ii, VMIN2(*nyf, jj+1), VMIN2(*nzf, kk+1)))
845 +0.125 *
HARMO2(VAT3(a3cf, VMAX2(1, ii-1), jj, kk),
846 VAT3(a3cf, VMAX2(1, ii-1), jj, VMIN2(*nzf, kk+1)))
847 +0.125 *
HARMO2(VAT3(a3cf, VMIN2(*nxf, ii+1), jj, kk),
848 VAT3(a3cf, VMIN2(*nxf, ii+1), jj, VMIN2(*nzf, kk+1)))
855 for (k=1; k<=*
nz; k++) {
858 for (j=1; j<=*
ny; j++) {
861 VAT3(gxc, j, k, 1) = VAT3(gxcf, jj, kk, 1);
862 VAT3(gxc, j, k, 2) = VAT3(gxcf, jj, kk, 2);
863 VAT3(gxc, j, k, 3) = VAT3(gxcf, jj, kk, 3);
864 VAT3(gxc, j, k, 4) = VAT3(gxcf, jj, kk, 4);
869 for (k=1; k<=*
nz; k++) {
872 for (i=1; i<=*
nx; i++) {
874 VAT3(gyc, i, k, 1) = VAT3(gycf, ii, kk, 1);
875 VAT3(gyc, i, k, 2) = VAT3(gycf, ii, kk, 2);
876 VAT3(gyc, i, k, 3) = VAT3(gycf, ii, kk, 3);
877 VAT3(gyc, i, k, 4) = VAT3(gycf, ii, kk, 4);
882 for (j=1; j<=*
ny; j++) {
885 for (i=1; i<=*
nx; i++) {
888 VAT3(gzc, i, j, 1) = VAT3(gzcf, ii, jj, 1);
889 VAT3(gzc, i, j, 2) = VAT3(gzcf, ii, jj, 2);
890 VAT3(gzc, i, j, 3) = VAT3(gzcf, ii, jj, 3);
891 VAT3(gzc, i, j, 4) = VAT3(gzcf, ii, jj, 4);
899 VPUBLIC
void Vbuildalg(
int *
nx,
int *
ny,
int *
nz,
900 int *mode,
int *
nlev,
int *iz,
901 int *ipc,
double *rpc,
902 double *ac,
double *cc,
double *fc,
903 double *x,
double *y,
double *tmp) {
906 int nxold, nyold, nzold;
918 if ((*mode == 1) || (*mode == 2)) {
920 RAT(ipc, VAT2(iz, 5, lev)), RAT(rpc, VAT2(iz, 6, lev)),
921 RAT( ac, VAT2(iz, 7, lev)), RAT( cc, VAT2(iz, 1, lev)),
922 RAT( x, VAT2(iz, 1, lev)), RAT( y, VAT2(iz, 1, lev)),
926 RAT(ipc, VAT2(iz, 5, lev)), RAT(rpc, VAT2(iz, 6, lev)),
927 RAT( ac, VAT2(iz, 7, lev)), RAT( cc, VAT2(iz, 1, lev)),
928 RAT( x, VAT2(iz, 1, lev)), RAT( y, VAT2(iz, 1, lev)));
932 for (lev=2; lev <= *
nlev; lev++) {
938 Vmkcors(&numlev, &nxold, &nyold, &nzold, &nxx, &nyy, &nzz);
941 if ((*mode == 1) || (*mode == 2)) {
943 RAT(ipc, VAT2(iz, 5, lev)), RAT(rpc, VAT2(iz, 6, lev)),
944 RAT( ac, VAT2(iz, 7, lev)), RAT( cc, VAT2(iz, 1, lev)),
945 RAT( x, VAT2(iz, 1, lev)), RAT( y, VAT2(iz, 1, lev)),
949 RAT(ipc, VAT2(iz, 5, lev)), RAT(rpc, VAT2(iz, 6, lev)),
950 RAT( ac, VAT2(iz, 7, lev)), RAT( cc, VAT2(iz, 1, lev)),
951 RAT( x, VAT2(iz, 1, lev)), RAT( y, VAT2(iz, 1, lev)));
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 Vbuildband(int *key, int *nx, int *ny, int *nz, int *ipc, double *rpc, double *ac, int *ipcB, double *rpcB, double *acB)
Banded matrix builder.
#define HARMO2(a, b)
Multigrid subroutines.
VPUBLIC void VbuildA(int *nx, int *ny, int *nz, int *ipkey, int *mgdisc, int *numdia, int *ipc, double *rpc, double *ac, double *cc, double *fc, double *xf, double *yf, double *zf, double *gxcf, double *gycf, double *gzcf, double *a1cf, double *a2cf, double *a3cf, double *ccf, double *fcf)
Build the Laplacian.
VPUBLIC void Vbuildops(int *nx, int *ny, int *nz, int *nlev, int *ipkey, int *iinfo, int *ido, int *iz, int *mgprol, int *mgcoar, int *mgsolv, int *mgdisc, int *ipc, double *rpc, double *pc, double *ac, double *cc, double *fc, double *xf, double *yf, double *zf, double *gxcf, double *gycf, double *gzcf, double *a1cf, double *a2cf, double *a3cf, double *ccf, double *fcf, double *tcf)
Build operators, boundary arrays, modify affine vectors ido==0: do only fine level ido==1: do only co...
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 VbuildP(int *nxf, int *nyf, int *nzf, int *nxc, int *nyc, int *nzc, int *mgprol, int *ipc, double *rpc, double *pc, double *ac, double *xf, double *yf, double *zf)
Builds prolongation matrix.
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 Vpackmg(int *iparm, double *rparm, size_t *nrwk, int *niwk, int *nx, int *ny, int *nz, int *nlev, int *nu1, int *nu2, int *mgkey, int *itmax, int *istop, int *ipcon, int *nonlin, int *mgsmoo, int *mgprol, int *mgcoar, int *mgsolv, int *mgdisc, int *iinfo, double *errtol, int *ipkey, double *omegal, double *omegan, int *irite, int *iperf)
Print out a column-compressed sparse matrix in Harwell-Boeing format.
VPUBLIC void VbuildG(int *nxf, int *nyf, int *nzf, int *nxc, int *nyc, int *nzc, int *numdia, double *pcFF, double *acFF, double *ac)
Build Galerkin matrix structures.
VPUBLIC void Vbuildgaler0(int *nxf, int *nyf, int *nzf, int *nxc, int *nyc, int *nzc, int *ipkey, int *numdia, double *pcFF, int *ipcFF, double *rpcFF, double *acFF, double *ccFF, double *fcFF, int *ipc, double *rpc, double *ac, double *cc, double *fc)
Form the Galerkin coarse grid system.
VPUBLIC void Vbuildstr(int *nx, int *ny, int *nz, int *nlev, int *iz)
Build the nexted operator framework in the array iz.
VPUBLIC void Vprtmatd(int *nx, int *ny, int *nz, int *ipc, double *rpc, double *ac)
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.