58 int *nx,
int *ny,
int *nz,
59 int *nlev,
int *ipkey,
int *iinfo,
61 int *mgprol,
int *mgcoar,
int *mgsolv,
int *mgdisc,
int *ipc,
62 double *rpc,
double *pc,
double *ac,
double *cc,
double *fc,
63 double *xf,
double *yf,
double *zf,
64 double *gxcf,
double *gycf,
double *gzcf,
65 double *a1cf,
double *a2cf,
double *a3cf,
66 double *ccf,
double *fcf,
double *tcf
91 if (*ido == 0 || *ido == 2) {
97 VMESSAGE3(
"Fine: (%03d, %03d, %03d)", nxx, nyy, nzz);
101 ipkey, mgdisc, &numdia,
102 RAT(ipc, VAT2(iz, 5,lev)), RAT(rpc, VAT2(iz, 6,lev)),
103 RAT(ac, VAT2(iz, 7,lev)), RAT(cc, VAT2(iz, 1,lev)), RAT(fc, VAT2(iz, 1,lev)),
104 RAT(xf, VAT2(iz, 8,lev)), RAT(yf, VAT2(iz, 9,lev)), RAT(zf, VAT2(iz, 10,lev)),
105 RAT(gxcf, VAT2(iz, 2,lev)), RAT(gycf, VAT2(iz, 3,lev)), RAT(gzcf, VAT2(iz, 4,lev)),
106 RAT(a1cf, VAT2(iz, 1,lev)), RAT(a2cf, VAT2(iz, 1,lev)), RAT(a3cf, VAT2(iz, 1,lev)),
107 RAT(ccf, VAT2(iz, 1,lev)), RAT(fcf, VAT2(iz, 1,lev)));
109 VMESSAGE2(
"Operator stencil (lev, numdia) = (%d, %d)", lev, numdia);
112 VAT2(iz, 7, lev+1) = VAT2(iz, 7, lev) + numdia * nxx * nyy * nzz;
117 RAT(ipc, VAT2(iz, 5,lev)), RAT(rpc, VAT2(iz, 6,lev)), RAT(ac, VAT2(iz, 7,lev)));
122 if (*ido == 1 || *ido == 2 || *ido == 3) {
124 for (lev=2; lev<=*nlev; lev++) {
131 Vmkcors(&i, &nxold, &nyold, &nzold, &nxx, &nyy, &nzz);
135 VbuildP(&nxold, &nyold, &nzold,
138 RAT(ipc, VAT2(iz, 5,lev-1)), RAT(rpc, VAT2(iz, 6,lev-1)),
139 RAT(pc, VAT2(iz, 11,lev-1)), RAT(ac, VAT2(iz, 7,lev-1)),
140 RAT(xf, VAT2(iz, 8,lev-1)), RAT(yf, VAT2(iz, 9,lev-1)), RAT(zf, VAT2(iz, 10,lev-1)));
147 VMESSAGE3(
"Stand: (%03d, %03d, %03d)", nxx, nyy, nzz);
151 Vbuildcopy0(&nxx, &nyy, &nzz,
152 &nxold, &nyold, &nzold,
153 RAT(xf, VAT2(iz, 8,lev )), RAT(yf, VAT2(iz, 9,lev )), RAT(zf, VAT2(iz, 10,lev )),
154 RAT(gxcf, VAT2(iz, 2,lev )), RAT(gycf, VAT2(iz, 3,lev )), RAT(gzcf, VAT2(iz, 4,lev )),
155 RAT(a1cf, VAT2(iz, 1,lev )), RAT(a2cf, VAT2(iz, 1,lev )), RAT(a3cf, VAT2(iz, 1,lev )),
156 RAT(ccf, VAT2(iz, 1,lev )), RAT(fcf, VAT2(iz, 1,lev )), RAT(tcf, VAT2(iz, 1,lev )),
157 RAT(xf, VAT2(iz, 8,lev-1)), RAT(yf, VAT2(iz, 9,lev-1)), RAT(zf, VAT2(iz, 10,lev-1)),
158 RAT(gxcf, VAT2(iz, 2,lev-1)), RAT(gycf, VAT2(iz, 3,lev-1)), RAT(gzcf, VAT2(iz, 4,lev-1)),
159 RAT(a1cf, VAT2(iz, 1,lev-1)), RAT(a2cf, VAT2(iz, 1,lev-1)), RAT(a3cf, VAT2(iz, 1,lev-1)),
160 RAT(ccf, VAT2(iz, 1,lev-1)), RAT(fcf, VAT2(iz, 1,lev-1)), RAT(tcf, VAT2(iz, 1,lev-1)));
163 ipkey, mgdisc, &numdia,
164 RAT(ipc, VAT2(iz, 5,lev)), RAT(rpc, VAT2(iz, 6,lev)),
165 RAT(ac, VAT2(iz, 7,lev)), RAT(cc, VAT2(iz, 1,lev)), RAT(fc, VAT2(iz, 1,lev)),
166 RAT(xf, VAT2(iz, 8,lev)), RAT(yf, VAT2(iz, 9,lev)), RAT(zf, VAT2(iz, 10,lev)),
167 RAT(gxcf, VAT2(iz, 2,lev)), RAT(gycf, VAT2(iz, 3,lev)), RAT(gzcf, VAT2(iz, 4,lev)),
168 RAT(a1cf, VAT2(iz, 1,lev)), RAT(a2cf, VAT2(iz, 1,lev)), RAT(a3cf, VAT2(iz, 1,lev)),
169 RAT(ccf, VAT2(iz, 1,lev)), RAT(fcf, VAT2(iz, 1,lev)));
173 else if (*mgcoar == 1) {
177 VMESSAGE3(
"Harm0: (%03d, %03d, %03d)", nxx, nyy, nzz);
179 Vbuildharm0(&nxx, &nyy, &nzz, &nxold, &nyold, &nzold,
180 RAT(xf, VAT2(iz, 8, lev )), RAT(yf, VAT2(iz, 9, lev )), RAT(zf, VAT2(iz, 10, lev )),
181 RAT(gxcf, VAT2(iz, 2, lev )), RAT(gycf, VAT2(iz, 3, lev )), RAT(gzcf, VAT2(iz, 4, lev )),
182 RAT(a1cf, VAT2(iz, 1, lev )), RAT(a2cf, VAT2(iz, 1, lev )), RAT(a3cf, VAT2(iz, 1, lev )),
183 RAT(ccf, VAT2(iz, 1, lev )), RAT(fcf, VAT2(iz, 1, lev )), RAT(tcf, VAT2(iz, 1, lev )),
184 RAT(xf, VAT2(iz, 8, lev-1)), RAT(yf, VAT2(iz, 9, lev-1)), RAT(zf, VAT2(iz, 10, lev-1)),
185 RAT(gxcf, VAT2(iz, 2, lev-1)), RAT(gycf, VAT2(iz, 3, lev-1)), RAT(gzcf, VAT2(iz, 4, lev-1)),
186 RAT(a1cf, VAT2(iz, 1, lev-1)), RAT(a2cf, VAT2(iz, 1, lev-1)), RAT(a3cf, VAT2(iz, 1, lev-1)),
187 RAT(ccf, VAT2(iz, 1, lev-1)), RAT(fcf, VAT2(iz, 1, lev-1)), RAT(tcf, VAT2(iz, 1, lev-1)));
190 ipkey, mgdisc, &numdia,
191 RAT(ipc, VAT2(iz, 5,lev)), RAT(rpc, VAT2(iz, 6,lev)),
192 RAT(ac, VAT2(iz, 7,lev)), RAT(cc, VAT2(iz, 1,lev)), RAT(fc, VAT2(iz, 1,lev)),
193 RAT(xf, VAT2(iz, 8,lev)), RAT(yf, VAT2(iz, 9,lev)), RAT(zf, VAT2(iz, 10,lev)),
194 RAT(gxcf, VAT2(iz, 2,lev)), RAT(gycf, VAT2(iz, 3,lev)), RAT(gzcf, VAT2(iz, 4,lev)),
195 RAT(a1cf, VAT2(iz, 1,lev)), RAT(a2cf, VAT2(iz, 1,lev)), RAT(a3cf, VAT2(iz, 1,lev)),
196 RAT(ccf, VAT2(iz, 1,lev)), RAT(fcf, VAT2(iz, 1,lev)));
200 else if (*mgcoar == 2) {
204 VMESSAGE3(
"Galer: (%03d, %03d, %03d)", nxx, nyy, nzz);
209 RAT(pc, VAT2(iz, 11,lev-1)),
210 RAT(ipc, VAT2(iz, 5,lev-1)), RAT(rpc, VAT2(iz, 6,lev-1)),
211 RAT(ac, VAT2(iz, 7,lev-1)), RAT(cc, VAT2(iz, 1,lev-1)), RAT(fc, VAT2(iz, 1,lev-1)),
212 RAT(ipc, VAT2(iz, 5,lev )), RAT(rpc, VAT2(iz, 6,lev )),
213 RAT(ac, VAT2(iz, 7,lev )), RAT(cc, VAT2(iz, 1,lev )), RAT(fc, VAT2(iz, 1,lev )));
217 Vextrac(&nxold, &nyold, &nzold,
219 RAT(tcf, VAT2(iz, 1,lev-1)), RAT(tcf, VAT2(iz, 1,lev)));
222 VABORT_MSG1(
"Bad mgcoar value given: %d", *mgcoar);
228 VAT2(iz, 7, lev+1) = VAT2(iz, 7,lev) + numdia * nxx * nyy * nzz;
234 RAT(ipc, VAT2(iz, 5,lev)), RAT(rpc, VAT2(iz, 6,lev)), RAT(ac, VAT2(iz, 7,lev)));
244 RAT(ipc, VAT2(iz, 5,lev )), RAT(rpc, VAT2(iz, 6,lev )), RAT(ac, VAT2(iz, 7,lev )),
245 RAT(ipc, VAT2(iz, 5,lev+1)), RAT(rpc, VAT2(iz, 6,lev+1)), RAT(ac, VAT2(iz, 7,lev+1)));
248 VERRMSG0(
"Changing your mgsolv to iterative");
257VPUBLIC
void Vbuildstr(
int *nx,
int *ny,
int *nz,
int *nlev,
int *iz) {
259 int nxold, nyold, nzold;
260 int nxnew, nynew, nznew;
272 n = nxnew * nynew * nznew;
278 VAT2(iz, 1, lev) = 1;
279 VAT2(iz, 2, lev) = 1;
280 VAT2(iz, 3, lev) = 1;
281 VAT2(iz, 4, lev) = 1;
282 VAT2(iz, 5, lev) = 1;
283 VAT2(iz, 6, lev) = 1;
284 VAT2(iz, 7, lev) = 1;
285 VAT2(iz, 8, lev) = 1;
286 VAT2(iz, 9, lev) = 1;
287 VAT2(iz, 10, lev) = 1;
288 VAT2(iz, 11, lev) = 1;
291 VAT2(iz, 1, lev + 1) = VAT2(iz, 1, lev) + n;
292 VAT2(iz, 2, lev + 1) = VAT2(iz, 2, lev) + 4 * nynew * nznew;
293 VAT2(iz, 3, lev + 1) = VAT2(iz, 3, lev) + 4 * nxnew * nznew;
294 VAT2(iz, 4, lev + 1) = VAT2(iz, 4, lev) + 4 * nxnew * nynew;
295 VAT2(iz, 5, lev + 1) = VAT2(iz, 5, lev) + 100;
296 VAT2(iz, 6, lev + 1) = VAT2(iz, 6, lev) + 100;
297 VAT2(iz, 8, lev + 1) = VAT2(iz, 8, lev) + nxnew;
298 VAT2(iz, 9, lev + 1) = VAT2(iz, 9, lev) + nynew;
299 VAT2(iz, 10, lev + 1) = VAT2(iz, 10, lev) + nznew;
310 for (lev=2; lev<=*nlev; lev++) {
315 Vmkcors(&numlev, &nxold, &nyold, &nzold, &nxnew, &nynew, &nznew);
316 n = nxnew * nynew * nznew;
319 VAT2(iz, 1, lev + 1) = VAT2(iz, 1, lev) + n;
320 VAT2(iz, 2, lev + 1) = VAT2(iz, 2, lev) + 4 * nynew * nznew;
321 VAT2(iz, 3, lev + 1) = VAT2(iz, 3, lev) + 4 * nxnew * nznew;
322 VAT2(iz, 4, lev + 1) = VAT2(iz, 4, lev) + 4 * nxnew * nynew;
323 VAT2(iz, 5, lev + 1) = VAT2(iz, 5, lev) + 100;
324 VAT2(iz, 6, lev + 1) = VAT2(iz, 6, lev) + 100;
325 VAT2(iz, 7, lev + 1) = VAT2(iz, 7, lev) + 4 * n;
326 VAT2(iz, 8, lev + 1) = VAT2(iz, 8, lev) + nxnew;
327 VAT2(iz, 9, lev + 1) = VAT2(iz, 9, lev) + nynew;
328 VAT2(iz, 10, lev + 1) = VAT2(iz, 10, lev) + nznew;
331 VAT2(iz, 11, lev) = VAT2(iz, 11, lev - 1) + 27 * n;
342 int *nxc,
int *nyc,
int *nzc,
343 int *ipkey,
int *numdia,
344 double *pcFF,
int *ipcFF,
double *rpcFF,
345 double *acFF,
double *ccFF,
double *fcFF,
346 int *ipc,
double *rpc,
347 double *ac,
double *cc,
double *fc) {
352 numdia_loc = VAT(ipcFF, 11);
363 VAT(ipc, 10) = *ipkey;
377VPUBLIC
void Vmkcors(
int *numlev,
378 int *nxold,
int *nyold,
int *nzold,
379 int *nxnew,
int *nynew,
int *nznew) {
380 int nxtmp, nytmp, nztmp;
388 for (i=1; i<=*numlev; i++) {
392 Vcorsr(&nxtmp,nxnew);
393 Vcorsr(&nytmp,nynew);
394 Vcorsr(&nztmp,nznew);
400VPUBLIC
void Vcorsr(
int *nold,
int *nnew) {
403 *nnew = (*nold - 1) / 2 + 1;
406 if ((*nnew - 1) * 2 != *nold - 1) {
407 Vnm_print(2,
"Vcorsr: WARNING! The grid dimensions are not\n");
408 Vnm_print(2,
"Vcorsr: consistent with nlev! Your\n");
409 Vnm_print(2,
"Vcorsr: calculation will only work if you\n");
410 Vnm_print(2,
"Vcorsr: are performing a mg-dummy run.\n");
414 Vnm_print(2,
"Vcorsr: ERROR! The grid dimensions are not\n");
415 Vnm_print(2,
"Vcorsr: consistent with nlev!\n");
416 Vnm_print(2,
"Vcorsr: Grid coarsened below zero.\n");
422VPUBLIC
void Vmkfine(
int *numlev,
423 int *nxold,
int *nyold,
int *nzold,
424 int *nxnew,
int *nynew,
int *nznew) {
426 int nxtmp, nytmp, nztmp, i;
433 for (i=1; i<=*numlev; i++) {
439 Vfiner(&nxtmp, nxnew);
440 Vfiner(&nytmp, nynew);
441 Vfiner(&nztmp, nznew);
448VPUBLIC
void Vfiner(
int *nold,
int *nnew) {
450 *nnew = (*nold - 1) * 2 + 1;
455VPUBLIC
int Vivariv(
int *nu,
int *level) {
471VPUBLIC
int Vmaxlev(
int n1,
int n2,
int n3) {
485 iden = (int)VPOW(2, lev - 1);
486 n1c = (n1 - 1) / iden + 1;
487 n2c = (n2 - 1) / iden + 1;
488 n3c = (n3 - 1) / iden + 1;
489 if (((n1c - 1) * iden != (n1 - 1)) || (n1c <= 2) )
491 if (((n2c - 1) * iden != (n2 - 1)) || (n2c <= 2) )
493 if (((n3c - 1) * iden != (n3 - 1)) || (n3c <= 2) )
496 }
while (idone != 1);
503VPUBLIC
void Vprtstp(
int iok,
int iters,
504 double rsnrm,
double rsden,
double orsnrm) {
507 double contrac = 0.0;
517 else if (iters == -1) {
518 Vnm_tstop(40,
"MG iteration");
531 VERRMSG0(
"Vprtstp: avoided division by zero\n");
533 relres = rsnrm / rsden;
539 VERRMSG0(
"avoided division by zero\n");
541 contrac = rsnrm / orsnrm;
545 if (iok == 1 || iok == 2) {
546 VMESSAGE1(
"iteration = %d", iters);
547 VMESSAGE1(
"relative residual = %e", relres);
548 VMESSAGE1(
"contraction number = %e", contrac);
555VPUBLIC
void Vpackmg(
int *iparm,
double *rparm,
size_t *nrwk,
int *niwk,
556 int *nx,
int *ny,
int *nz,
int *nlev,
int *nu1,
int *nu2,
int *mgkey,
557 int *itmax,
int *istop,
int *ipcon,
int *nonlin,
int *mgsmoo,
int *mgprol,
558 int *mgcoar,
int *mgsolv,
int *mgdisc,
int *iinfo,
double *errtol,
559 int *ipkey,
double *omegal,
double *omegan,
int *irite,
int *iperf) {
564 VAT(iparm, 1) = *nrwk;
565 VAT(iparm, 2) = *niwk;
569 VAT(iparm, 6) = *nlev;
570 VAT(iparm, 7) = *nu1;
571 VAT(iparm, 8) = *nu2;
572 VAT(iparm, 9) = *mgkey;
573 VAT(iparm, 10) = *itmax;
574 VAT(iparm, 11) = *istop;
575 VAT(iparm, 12) = *iinfo;
576 VAT(iparm, 13) = *irite;
577 VAT(iparm, 14) = *ipkey;
578 VAT(iparm, 15) = *ipcon;
579 VAT(iparm, 16) = *nonlin;
580 VAT(iparm, 17) = *mgprol;
581 VAT(iparm, 18) = *mgcoar;
582 VAT(iparm, 19) = *mgdisc;
583 VAT(iparm, 20) = *mgsmoo;
584 VAT(iparm, 21) = *mgsolv;
585 VAT(iparm, 22) = *iperf;
588 VAT(rparm, 1) = *errtol;
589 VAT(rparm, 9) = *omegal;
590 VAT(rparm, 10) = *omegan;
595VEXTERNC
void Vbuildcopy0(
int *nx,
int *ny,
int *nz,
596 int *nxf,
int *nyf,
int *nzf,
597 double *xc,
double *yc,
double *zc,
598 double *gxc,
double *gyc,
double *gzc,
599 double *a1c,
double *a2c,
double *a3c,
600 double *cc,
double *fc,
double *tc,
601 double *xf,
double *yf,
double *zf,
602 double *gxcf,
double *gycf,
double *gzcf,
603 double *a1cf,
double *a2cf,
double *a3cf,
604 double *ccf,
double *fcf,
double *tcf) {
609 int iadd, jadd, kadd;
611 MAT3( gxc, *ny, *nz, 4);
612 MAT3( gyc, *nx, *nz, 4);
613 MAT3( gzc, *nx, *ny, 4);
614 MAT3( a1c, *nx, *ny, *nz);
615 MAT3( a2c, *nx, *ny, *nz);
616 MAT3( a3c, *nx, *ny, *nz);
617 MAT3( cc, *nx, *ny, *nz);
618 MAT3( fc, *nx, *ny, *nz);
619 MAT3( tc, *nx, *ny, *nz);
620 MAT3(gxcf, *nyf, *nzf, 4);
621 MAT3(gycf, *nxf, *nzf, 4);
622 MAT3(gzcf, *nxf, *nyf, 4);
623 MAT3(a1cf, *nxf, *nyf, *nzf);
624 MAT3(a2cf, *nxf, *nyf, *nzf);
625 MAT3(a3cf, *nxf, *nyf, *nzf);
626 MAT3( tcf, *nxf, *nyf, *nzf);
627 MAT3( ccf, *nxf, *nyf, *nzf);
628 MAT3( fcf, *nxf, *nyf, *nzf);
633 iadd = (*nxf - 1) / (*nx - 1);
634 jadd = (*nyf - 1) / (*ny - 1);
635 kadd = (*nzf - 1) / (*nz - 1);
637 if (iadd != 2 || jadd != 2 || kadd != 2) {
638 Vnm_print(2,
"Vbuildcopy0: Problem with grid dimensions...\n");
642 for (k=1; k<=*nz; k++) {
644 VAT(zc, k) = VAT(zf, kk);
646 for (j=1; j<=*ny; j++) {
648 VAT(yc, j) = VAT(yf, jj);
650 for (i=1; i<=*nx; i++){
652 VAT(xc, i) = VAT(xf, ii);
655 VAT3( tc, i, j, k) = VAT3( tcf, ii, jj, kk);
658 VAT3( cc, i, j, k) = VAT3( ccf, ii, jj, kk);
661 VAT3( fc, i, j, k) = VAT3( fcf, ii, jj, kk);
664 VAT3(a1c, i, j, k) = VAT3(a1cf, ii, jj, kk);
667 VAT3(a2c, i, j, k) = VAT3(a2cf, ii, jj, kk);
670 VAT3(a3c, i, j, k) = VAT3(a3cf, ii, jj, kk);
676 for (k=1; k<=*nz; k++) {
679 for (j=1; j<=*ny; j++) {
682 VAT3(gxc, j, k, 1) = VAT3(gxcf, jj, kk, 1);
683 VAT3(gxc, j, k, 2) = VAT3(gxcf, jj, kk, 2);
684 VAT3(gxc, j, k, 3) = VAT3(gxcf, jj, kk, 3);
685 VAT3(gxc, j, k, 4) = VAT3(gxcf, jj, kk, 4);
690 for (k=1; k<=*nz; k++) {
693 for (i=1; i<=*nx; i++) {
696 VAT3(gyc, i, k, 1) = VAT3(gycf, ii, kk, 1);
697 VAT3(gyc, i, k, 2) = VAT3(gycf, ii, kk, 2);
698 VAT3(gyc, i, k, 3) = VAT3(gycf, ii, kk, 3);
699 VAT3(gyc, i, k, 4) = VAT3(gycf, ii, kk, 4);
704 for (j=1; j<=*ny; j++) {
707 for (i=1; i<=*nx; i++) {
710 VAT3(gzc, i, j, 1) = VAT3(gzcf, ii, jj, 1);
711 VAT3(gzc, i, j, 2) = VAT3(gzcf, ii, jj, 2);
712 VAT3(gzc, i, j, 3) = VAT3(gzcf, ii, jj, 3);
713 VAT3(gzc, i, j, 4) = VAT3(gzcf, ii, jj, 4);
718VPUBLIC
void Vbuildharm0(
int *nx,
int *ny,
int *nz,
719 int *nxf,
int *nyf,
int *nzf,
720 double *xc,
double *yc,
double *zc,
721 double *gxc,
double *gyc,
double *gzc,
722 double *a1c,
double *a2c,
double *a3c,
723 double *cc,
double *fc,
double *tc,
724 double *xf,
double *yf,
double *zf,
725 double *gxcf,
double *gycf,
double *gzcf,
726 double *a1cf,
double *a2cf,
double *a3cf,
727 double *ccf,
double *fcf,
double *tcf) {
729 Vnm_print(2,
"WARNING: FUNCTION IS NOT FULLY IMPLEMENTED YET!!!");
733 int iadd, jadd, kadd;
735 MAT3( gxc, *ny, *nz, 4);
736 MAT3( gyc, *nx, *nz, 4);
737 MAT3( gzc, *nx, *ny, 4);
739 MAT3( a1c, *nx, *ny, *nz);
740 MAT3( a2c, *nx, *ny, *nz);
741 MAT3( a3c, *nx, *ny, *nz);
743 MAT3( cc, *nx, *ny, *nz);
744 MAT3( fc, *nx, *ny, *nz);
745 MAT3( tc, *nx, *ny, *nz);
747 MAT3(gxcf, *nyf, *nzf, 4);
748 MAT3(gycf, *nxf, *nzf, 4);
749 MAT3(gzcf, *nxf, *nyf, 4);
750 MAT3(a1cf, *nxf, *nyf, *nzf);
751 MAT3(a2cf, *nxf, *nyf, *nzf);
752 MAT3(a3cf, *nxf, *nyf, *nzf);
753 MAT3( tcf, *nxf, *nyf, *nzf);
754 MAT3( ccf, *nxf, *nyf, *nzf);
755 MAT3( fcf, *nxf, *nyf, *nzf);
759 double a, b, c, d, e, f, g, h;
762 iadd = (*nxf - 1) / (*nx - 1);
763 jadd = (*nyf - 1) / (*ny - 1);
764 kadd = (*nzf - 1) / (*nz - 1);
765 if (iadd !=2 || jadd != 2 || kadd != 2) {
766 Vnm_print(2,
"BUILDHARM0: problem with grid dimensions...\n");
770 for (k=1; k<=*nz; k++) {
772 VAT(zc, k) = VAT(zf, kk);
774 for (j=1; j<=*ny; j++) {
776 VAT(yc, j) = VAT(yf,jj);
778 for (i=1; i<=*nx; i++) {
780 VAT(xc, i) = VAT(xf, ii);
783 VAT3(tc, i, j, k) = VAT3(tcf, ii, jj, kk);
786 VAT3(cc, i, j, k) = VAT3(ccf, ii, jj, kk);
801 VAT3(fc, i, j, k) = VAT3(fcf, ii, jj, kk);
815 VAT3(a1c, i, j, k) = (
816 +0.500 *
HARMO2(VAT3(a1cf, ii, jj, kk),
817 VAT3(a1cf, VMIN2(*nxf, ii+1), jj, kk))
818 +0.125 *
HARMO2(VAT3(a1cf, ii, jj, VMAX2(1, kk-1)),
819 VAT3(a1cf, VMIN2(*nxf, ii+1), jj, VMAX2(1, kk-1)))
820 +0.125 *
HARMO2(VAT3(a1cf, ii, jj, VMIN2(*nzf, kk+1)),
821 VAT3(a1cf, VMIN2(*nxf, ii+1), jj, VMIN2(*nzf, kk+1)))
822 +0.125 *
HARMO2(VAT3(a1cf, ii, VMAX2(1, jj-1), kk),
823 VAT3(a1cf, VMIN2(*nxf, ii+1), VMAX2(1, jj-1), kk))
824 +0.125 *
HARMO2(VAT3(a1cf, ii, VMIN2(*nyf, jj+1), kk),
825 VAT3(a1cf, VMIN2(*nxf, ii+1), VMIN2(*nyf, jj+1), kk))
829 VAT3(a2c, i, j, k) = (
830 +0.500 *
HARMO2(VAT3(a2cf, ii, jj, kk),
831 VAT3(a2cf, ii, VMIN2(*nyf, jj+1), kk))
832 +0.125 *
HARMO2(VAT3(a2cf, ii, jj, VMAX2(1, kk-1)),
833 VAT3(a2cf, ii, VMIN2(*nyf, jj+1), VMAX2(1, kk-1)))
834 +0.125 *
HARMO2(VAT3(a2cf, ii, jj, VMIN2(*nzf, kk+1)),
835 VAT3(a2cf, ii, VMIN2(*nyf, jj+1), VMIN2(*nzf, kk+1)))
836 +0.125 *
HARMO2(VAT3(a2cf, VMAX2(1, ii-1), jj, kk),
837 VAT3(a2cf, VMAX2(1, ii-1), VMIN2(nyf, jj+1), kk))
838 +0.125 *
HARMO2(VAT3(a2cf, VMIN2(*nxf, ii+1), jj, kk),
839 VAT3(a2cf, VMIN2(*nxf, ii+1), VMIN2(*nyf, jj+1), kk))
843 VAT3(a3c, i, j, k) = (
844 +0.500 *
HARMO2(VAT3(a3cf, ii, jj, kk),
845 VAT3(a3cf, ii, jj, VMIN2(*nzf, kk+1)))
846 +0.125 *
HARMO2(VAT3(a3cf, ii, VMAX2(1, jj-1), kk),
847 VAT3(a3cf, ii, VMAX2(1, jj-1), VMIN2(*nzf, kk+1)))
848 +0.125 *
HARMO2(VAT3(a3cf, ii, VMIN2(*nyf, jj+1), kk),
849 VAT3(a3cf, ii, VMIN2(*nyf, jj+1), VMIN2(*nzf, kk+1)))
850 +0.125 *
HARMO2(VAT3(a3cf, VMAX2(1, ii-1), jj, kk),
851 VAT3(a3cf, VMAX2(1, ii-1), jj, VMIN2(*nzf, kk+1)))
852 +0.125 *
HARMO2(VAT3(a3cf, VMIN2(*nxf, ii+1), jj, kk),
853 VAT3(a3cf, VMIN2(*nxf, ii+1), jj, VMIN2(*nzf, kk+1)))
860 for (k=1; k<=*nz; k++) {
863 for (j=1; j<=*ny; j++) {
866 VAT3(gxc, j, k, 1) = VAT3(gxcf, jj, kk, 1);
867 VAT3(gxc, j, k, 2) = VAT3(gxcf, jj, kk, 2);
868 VAT3(gxc, j, k, 3) = VAT3(gxcf, jj, kk, 3);
869 VAT3(gxc, j, k, 4) = VAT3(gxcf, jj, kk, 4);
874 for (k=1; k<=*nz; k++) {
877 for (i=1; i<=*nx; i++) {
879 VAT3(gyc, i, k, 1) = VAT3(gycf, ii, kk, 1);
880 VAT3(gyc, i, k, 2) = VAT3(gycf, ii, kk, 2);
881 VAT3(gyc, i, k, 3) = VAT3(gycf, ii, kk, 3);
882 VAT3(gyc, i, k, 4) = VAT3(gycf, ii, kk, 4);
887 for (j=1; j<=*ny; j++) {
890 for (i=1; i<=*nx; i++) {
893 VAT3(gzc, i, j, 1) = VAT3(gzcf, ii, jj, 1);
894 VAT3(gzc, i, j, 2) = VAT3(gzcf, ii, jj, 2);
895 VAT3(gzc, i, j, 3) = VAT3(gzcf, ii, jj, 3);
896 VAT3(gzc, i, j, 4) = VAT3(gzcf, ii, jj, 4);
904VPUBLIC
void Vbuildalg(
int *nx,
int *ny,
int *nz,
905 int *mode,
int *nlev,
int *iz,
906 int *ipc,
double *rpc,
907 double *ac,
double *cc,
double *fc,
908 double *x,
double *y,
double *tmp) {
911 int nxold, nyold, nzold;
923 if ((*mode == 1) || (*mode == 2)) {
925 RAT(ipc, VAT2(iz, 5, lev)), RAT(rpc, VAT2(iz, 6, lev)),
926 RAT( ac, VAT2(iz, 7, lev)), RAT( cc, VAT2(iz, 1, lev)),
927 RAT( x, VAT2(iz, 1, lev)), RAT( y, VAT2(iz, 1, lev)),
931 RAT(ipc, VAT2(iz, 5, lev)), RAT(rpc, VAT2(iz, 6, lev)),
932 RAT( ac, VAT2(iz, 7, lev)), RAT( cc, VAT2(iz, 1, lev)),
933 RAT( x, VAT2(iz, 1, lev)), RAT( y, VAT2(iz, 1, lev)));
937 for (lev=2; lev <= *nlev; lev++) {
943 Vmkcors(&numlev, &nxold, &nyold, &nzold, &nxx, &nyy, &nzz);
946 if ((*mode == 1) || (*mode == 2)) {
948 RAT(ipc, VAT2(iz, 5, lev)), RAT(rpc, VAT2(iz, 6, lev)),
949 RAT( ac, VAT2(iz, 7, lev)), RAT( cc, VAT2(iz, 1, lev)),
950 RAT( x, VAT2(iz, 1, lev)), RAT( y, VAT2(iz, 1, lev)),
954 RAT(ipc, VAT2(iz, 5, lev)), RAT(rpc, VAT2(iz, 6, lev)),
955 RAT( ac, VAT2(iz, 7, lev)), RAT( cc, VAT2(iz, 1, lev)),
956 RAT( x, VAT2(iz, 1, lev)), RAT( y, VAT2(iz, 1, lev)));
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...
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 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 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 Vprtmatd(int *nx, int *ny, int *nz, int *ipc, double *rpc, double *ac)
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 Vbuildstr(int *nx, int *ny, int *nz, int *nlev, int *iz)
Build the nexted operator framework in the array iz.
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 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 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.
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.