157VPUBLIC
void Vmvfas(
int *nx,
int *ny,
int *nz,
160 double *w0,
double *w1,
double *w2,
double *w3,
double *w4,
161 int *istop,
int *itmax,
int *iters,
int *ierror,
162 int *nlev,
int *ilev,
int *nlev_real,
163 int *mgsolv,
int *iok,
int *iinfo,
164 double *epsiln,
double *errtol,
double *omega,
165 int *nu1,
int *nu2,
int *mgsmoo,
166 int *ipc,
double *rpc,
167 double *pc,
double *ac,
double *cc,
double *fc,
double *tru) {
171 int itmax_s, iters_s, nuuu, ivariv, mgsmoo_s, iresid;
176 double rsden, rsnrm, orsnrm;
188 lev = (*ilev - 1) + level;
196 Vmkcors(&numlev, &nxf, &nyf, &nzf, &nxc, &nyc, &nzc);
200 Vnm_print(2,
"Vmvfas: starting: (%d, %d, %d) (%d, %d, %d)\n",
201 nxf, nyf, nzf, nxc, nyc, nzc);
206 Vprtstp(*iok, -1, 0.0, 0.0, 0.0);
223 }
else if (*istop == 1) {
231 RAT(ipc, VAT2(iz, 5, lev)), RAT(rpc, VAT2(iz, 6, lev)),
232 RAT( ac, VAT2(iz, 7, lev)), RAT( cc, VAT2(iz, 1, lev)),
233 RAT( fc, VAT2(iz, 1, lev)),
236 rsden =
Vxnrm1(&nxf, &nyf, &nzf, w2);
238 }
else if (*istop == 2) {
239 rsden = VSQRT(nxf * nyf * nzf);
240 }
else if (*istop == 3) {
241 rsden =
Vxnrm2(&nxf, &nyf, &nzf, RAT(tru, VAT2(iz, 1, lev)));
242 }
else if (*istop == 4) {
243 rsden =
Vxnrm2(&nxf, &nyf, &nzf, RAT(tru, VAT2(iz, 1, lev)));
244 }
else if (*istop == 5) {
246 RAT(ipc, VAT2(iz, 5, lev)), RAT(rpc, VAT2(iz, 6, lev)),
247 RAT( ac, VAT2(iz, 7, lev)), RAT( cc, VAT2(iz, 1, lev)),
248 RAT(tru, VAT2(iz, 1, lev)),
250 rsden = VSQRT(
Vxdot(&nxf, &nyf, &nzf,RAT(tru, VAT2(iz, 1, lev)), w1));
252 Vnm_print(2,
"Vmvfas: bad istop value: %d\n", *istop);
256 Vnm_print(2,
"Vmfas: rhs is zero on finest level\n");
261 Vprtstp(*iok, 0, rsnrm, rsden, orsnrm);
280 Vazeros(&nxf, &nyf, &nzf,RAT(x, VAT2(iz, 1, lev)));
282 RAT(ipc, VAT2(iz, 5, lev)), RAT(rpc, VAT2(iz, 6, lev)),
283 RAT( ac, VAT2(iz, 7, lev)), RAT( cc, VAT2(iz, 1, lev)),
284 RAT( fc, VAT2(iz, 1, lev)),
285 RAT( x, VAT2(iz, 1, lev)),
287 &itmax_s, &iters_s, &errtol_s, omega,
288 &iresid, &iadjoint, &mgsmoo_s);
296 RAT(ipc, VAT2(iz, 5, lev)),
297 RAT(rpc, VAT2(iz, 6, lev)),
298 RAT( ac, VAT2(iz, 7, lev)),
299 RAT( cc, VAT2(iz, 1, lev)),
300 RAT( fc, VAT2(iz, 1, lev)),
301 RAT( x, VAT2(iz, 1, lev)),
303 rsnrm =
Vxnrm1(&nxf, &nyf, &nzf, w1);
304 }
else if (*istop == 1) {
306 RAT(ipc, VAT2(iz, 5, lev)),
307 RAT(rpc, VAT2(iz, 6, lev)),
308 RAT( ac, VAT2(iz, 7, lev)),
309 RAT( cc, VAT2(iz, 1, lev)),
310 RAT( fc, VAT2(iz, 1, lev)),
311 RAT( x, VAT2(iz, 1, lev)),
313 rsnrm =
Vxnrm1(&nxf, &nyf, &nzf, w1);
314 }
else if (*istop == 2) {
316 RAT(tru, VAT2(iz, 1, lev)),
320 &alpha, RAT(x, VAT2(iz, 1, lev)), w1);
321 rsnrm =
Vxnrm1(&nxf, &nyf, &nzf, w1);
323 RAT(x, VAT2(iz, 1, lev)),
324 RAT(tru, VAT2(iz, 1, lev)));
325 }
else if (*istop == 3) {
327 RAT(tru, VAT2(iz, 1, lev)),
331 &alpha, RAT(x, VAT2(iz, 1, lev)), w1);
332 rsnrm =
Vxnrm2(&nxf, &nyf, &nzf, w1);
333 }
else if (*istop == 4) {
335 RAT(tru, VAT2(iz, 1, lev)),
339 &alpha, RAT(x, VAT2(iz, 1, lev)), w1);
340 rsnrm =
Vxnrm2(&nxf, &nyf, &nzf, w1);
341 }
else if (*istop == 5) {
343 RAT(tru, VAT2(iz, 1, lev)),
347 &alpha, RAT(x, VAT2(iz, 1, lev)), w1);
349 RAT(ipc, VAT2(iz, 5, lev)),
350 RAT(rpc, VAT2(iz, 6, lev)),
351 RAT( ac, VAT2(iz, 7, lev)),
352 RAT( cc, VAT2(iz, 1, lev)),
354 rsnrm = VSQRT(
Vxdot(&nxf, &nyf, &nzf, w1, w2));
356 Vnm_print(2,
"Vmvcs: bad istop value: %d\n", *istop);
359 Vprtstp(*iok, *iters, rsnrm, rsden, orsnrm);
374 lev = (*ilev - 1) + level;
381 nuuu = Vivariv(nu1, &lev);
383 RAT(ipc, VAT2(iz, 5, lev)), RAT(rpc, VAT2(iz, 6, lev)),
384 RAT( ac, VAT2(iz, 7, lev)), RAT( cc, VAT2(iz, 1, lev)),
385 RAT( fc, VAT2(iz, 1, lev)), RAT( x, VAT2(iz, 1, lev)),
387 &nuuu, &iters_s, &errtol_s, omega,
388 &iresid, &iadjoint, mgsmoo);
389 Vxcopy(&nxf, &nyf, &nzf, w1, RAT(w0, VAT2(iz, 1, lev)));
396 for (level = 2; level <= *nlev; level++) {
397 lev = (*ilev - 1) + level;
401 Vmkcors(&numlev, &nxf, &nyf, &nzf, &nxc, &nyc, &nzc);
404 Vrestrc(&nxf, &nyf, &nzf, &nxc, &nyc, &nzc,
405 w1, RAT(w0, VAT2(iz, 1, lev)),
406 RAT(pc, VAT2(iz, 11, lev-1)));
409 Vextrac(&nxf, &nyf, &nzf, &nxc, &nyc, &nzc,
410 RAT(x, VAT2(iz, 1, lev-1)), RAT(w4, VAT2(iz, 1, lev)));
419 RAT(ipc, VAT2(iz, 5, lev)), RAT(rpc, VAT2(iz, 6, lev)),
420 RAT( ac, VAT2(iz, 7, lev)), RAT( cc, VAT2(iz, 1, lev)),
421 RAT( w4, VAT2(iz, 1, lev)), RAT( fc, VAT2(iz, 1, lev)),
426 Vxaxpy(&nxf, &nyf, &nzf, &alpha,
427 RAT(w0, VAT2(iz, 1, lev)), RAT(fc, VAT2(iz, 1, lev)));
430 if (level != *nlev) {
434 RAT(w4, VAT2(iz, 1, lev)), RAT(x, VAT2(iz, 1, lev)));
439 nuuu = Vivariv (nu1,&lev);
441 RAT(ipc, VAT2(iz, 5, lev)), RAT(rpc, VAT2(iz, 6, lev)),
442 RAT( ac, VAT2(iz, 7, lev)), RAT( cc, VAT2(iz, 1, lev)),
443 RAT( fc, VAT2(iz, 1, lev)),
444 RAT( x, VAT2(iz, 1, lev)),
446 &nuuu, &iters_s, &errtol_s, omega,
447 &iresid, &iadjoint, mgsmoo);
461 lev = (*ilev - 1) + level;
470 Vazeros(&nxf, &nyf, &nzf, RAT(x, VAT2(iz, 1, lev)));
472 RAT(ipc, VAT2(iz, 5, lev)), RAT(rpc, VAT2(iz, 6, lev)),
473 RAT( ac, VAT2(iz, 7, lev)), RAT( cc, VAT2(iz, 1, lev)),
474 RAT( fc, VAT2(iz, 1, lev)),
475 RAT( x, VAT2(iz, 1, lev)),
477 &itmax_s, &iters_s, &errtol_s, omega,
478 &iresid, &iadjoint, &mgsmoo_s);
485 for (level = *nlev - 1; level >= 1; level--) {
486 lev = (*ilev - 1) + level;
490 Vmkfine(&numlev, &nxf, &nyf, &nzf, &nxc, &nyc, &nzc);
494 Vxaxpy(&nxf, &nyf, &nzf, &alpha,
495 RAT(w4, VAT2(iz, 1, lev + 1)), RAT(x, VAT2(iz, 1, lev + 1)));
498 Vlinesearch(&nxf, &nyf, &nzf, &xdamp,
499 RAT(ipc, VAT2(iz, 5, lev + 1)),
500 RAT(rpc, VAT2(iz, 6, lev + 1)),
501 RAT( ac, VAT2(iz, 7, lev + 1)),
502 RAT( cc, VAT2(iz, 1, lev + 1)),
503 RAT( fc, VAT2(iz, 1, lev + 1)),
504 RAT( x, VAT2(iz, 1, lev + 1)),
505 RAT( w4, VAT2(iz, 1, lev + 1)),
506 RAT( w0, VAT2(iz, 1, lev + 1)),
510 VinterpPMG(&nxf, &nyf, &nzf, &nxc, &nyc, &nzc,
511 RAT(x, VAT2(iz, 1, lev + 1)),
513 RAT(pc, VAT2(iz, 11, lev)));
521 Vxaxpy(&nxf, &nyf, &nzf, &xdamp, w1, RAT(x, VAT2(iz, 1, lev)));
528 nuuu = Vivariv(nu2, &lev);
530 RAT(ipc, VAT2(iz, 5, lev)), RAT(rpc, VAT2(iz, 6, lev)),
531 RAT( ac, VAT2(iz, 7, lev)), RAT( cc, VAT2(iz, 1, lev)),
532 RAT( fc, VAT2(iz, 1, lev)),
533 RAT( x, VAT2(iz, 1, lev)),
535 &nuuu, &iters_s, &errtol_s, omega,
536 &iresid, &iadjoint, mgsmoo);
552 RAT(ipc, VAT2(iz, 5, lev)), RAT(rpc, VAT2(iz, 6, lev)),
553 RAT( ac, VAT2(iz, 7, lev)), RAT( cc, VAT2(iz, 1, lev)),
554 RAT( fc, VAT2(iz, 1, lev)), RAT( x, VAT2(iz, 1, lev)),
556 rsnrm =
Vxnrm1(&nxf, &nyf, &nzf, w1);
557 }
else if (*istop == 1) {
559 RAT(ipc, VAT2(iz, 5, lev)), RAT(rpc, VAT2(iz, 6, lev)),
560 RAT( ac, VAT2(iz, 7, lev)), RAT( cc, VAT2(iz, 1, lev)),
561 RAT( fc, VAT2(iz, 1, lev)), RAT( x, VAT2(iz, 1, lev)),
563 rsnrm =
Vxnrm1(&nxf, &nyf, &nzf,w1);
564 }
else if (*istop == 2) {
565 Vxcopy(&nxf, &nyf, &nzf, RAT(tru, VAT2(iz, 1, lev)), w1);
567 Vxaxpy(&nxf, &nyf, &nzf, &alpha, RAT(x, VAT2(iz, 1, lev)), w1);
568 rsnrm =
Vxnrm1(&nxf, &nyf, &nzf, w1);
570 RAT( x, VAT2(iz, 1, lev)),
571 RAT(tru, VAT2(iz, 1, lev)));
572 }
else if (*istop == 3) {
573 Vxcopy(&nxf, &nyf, &nzf, RAT(tru, VAT2(iz, 1, lev)), w1);
575 Vxaxpy(&nxf, &nyf, &nzf, &alpha, RAT(x, VAT2(iz, 1, lev)), w1);
576 rsnrm =
Vxnrm2(&nxf, &nyf, &nzf,w1);
577 }
else if (*istop == 4) {
578 Vxcopy(&nxf, &nyf, &nzf, RAT(tru, VAT2(iz, 1, lev)), w1);
580 Vxaxpy(&nxf, &nyf, &nzf, &alpha, RAT(x, VAT2(iz, 1, lev)), w1);
581 rsnrm =
Vxnrm2(&nxf, &nyf, &nzf,w1);
582 }
else if (*istop == 5) {
583 Vxcopy(&nxf, &nyf, &nzf, RAT(tru, VAT2(iz, 1, lev)), w1);
585 Vxaxpy(&nxf, &nyf, &nzf, &alpha, RAT(x, VAT2(iz, 1, lev)), w1);
587 RAT(ipc, VAT2(iz, 5, lev)), RAT(rpc, VAT2(iz, 6, lev)),
588 RAT( ac, VAT2(iz, 7, lev)), RAT( cc, VAT2(iz, 1, lev)),
590 rsnrm = VSQRT(
Vxdot(&nxf, &nyf, &nzf,w1,w2));
592 VABORT_MSG1(
"Bad istop value: %d", *istop);
595 Vprtstp(*iok, *iters, rsnrm, rsden, orsnrm);
597 if ((rsnrm / rsden) <= *errtol)
601 if (*iters >= *itmax) {