55 double *w0,
double *w1,
double *w2,
double *w3,
double *w4,
56 int *
istop,
int *
itmax,
int *iters,
int *ierror,
57 int *
nlev,
int *ilev,
int *nlev_real,
59 double *epsiln,
double *
errtol,
double *omega,
61 int *ipc,
double *rpc,
62 double *pc,
double *ac,
double *cc,
double *fc,
double *tru) {
65 int level, itmxd, nlevd, iterd, iokd;
82 Vmkcors(&numlev, &nxf, &nyf, &nzf, &
nxc, &
nyc, &
nzc);
87 Vnm_print(2,
"Vfmvfas: starting: (%d,%d,%d) (%d,%d,%d)\n",
91 for (level = *nlev_real; level >= *ilev + 1; level--) {
96 nlevd = *nlev_real - level + 1;
108 &istpd, &itmxd, &iterd, ierror,
109 &nlevd, &level, nlev_real,
114 pc, ac, cc, fc, tru);
118 Vmkfine(&numlev, &
nxc, &
nyc, &
nzc, &nxf, &nyf, &nzf);
123 RAT( x, VAT2(iz, 1, level)),
124 RAT( x, VAT2(iz, 1, level-1)),
125 RAT(pc, VAT2(iz, 11, level-1)));
142 ierror,
nlev, &level, nlev_real,
147 pc, ac, cc, fc, tru);
155 double *w0,
double *w1,
double *w2,
double *w3,
double *w4,
156 int *
istop,
int *
itmax,
int *iters,
int *ierror,
157 int *
nlev,
int *ilev,
int *nlev_real,
159 double *epsiln,
double *
errtol,
double *omega,
161 int *ipc,
double *rpc,
162 double *pc,
double *ac,
double *cc,
double *fc,
double *tru) {
166 int itmax_s, iters_s, nuuu, ivariv, mgsmoo_s, iresid;
171 double rsden, rsnrm, orsnrm;
183 lev = (*ilev - 1) + level;
191 Vmkcors(&numlev, &nxf, &nyf, &nzf, &
nxc, &
nyc, &
nzc);
195 Vnm_print(2,
"Vmvfas: starting: (%d, %d, %d) (%d, %d, %d)\n",
201 Vprtstp(*iok, -1, 0.0, 0.0, 0.0);
218 }
else if (*
istop == 1) {
226 RAT(ipc, VAT2(iz, 5, lev)), RAT(rpc, VAT2(iz, 6, lev)),
227 RAT( ac, VAT2(iz, 7, lev)), RAT( cc, VAT2(iz, 1, lev)),
228 RAT( fc, VAT2(iz, 1, lev)),
231 rsden =
Vxnrm1(&nxf, &nyf, &nzf, w2);
233 }
else if (*
istop == 2) {
234 rsden = VSQRT(nxf * nyf * nzf);
235 }
else if (*
istop == 3) {
236 rsden =
Vxnrm2(&nxf, &nyf, &nzf, RAT(tru, VAT2(iz, 1, lev)));
237 }
else if (*
istop == 4) {
238 rsden =
Vxnrm2(&nxf, &nyf, &nzf, RAT(tru, VAT2(iz, 1, lev)));
239 }
else if (*
istop == 5) {
241 RAT(ipc, VAT2(iz, 5, lev)), RAT(rpc, VAT2(iz, 6, lev)),
242 RAT( ac, VAT2(iz, 7, lev)), RAT( cc, VAT2(iz, 1, lev)),
243 RAT(tru, VAT2(iz, 1, lev)),
245 rsden = VSQRT(
Vxdot(&nxf, &nyf, &nzf,RAT(tru, VAT2(iz, 1, lev)), w1));
247 Vnm_print(2,
"Vmvfas: bad istop value: %d\n", *
istop);
251 Vnm_print(2,
"Vmfas: rhs is zero on finest level\n");
256 Vprtstp(*iok, 0, rsnrm, rsden, orsnrm);
275 Vazeros(&nxf, &nyf, &nzf,RAT(x, VAT2(iz, 1, lev)));
277 RAT(ipc, VAT2(iz, 5, lev)), RAT(rpc, VAT2(iz, 6, lev)),
278 RAT( ac, VAT2(iz, 7, lev)), RAT( cc, VAT2(iz, 1, lev)),
279 RAT( fc, VAT2(iz, 1, lev)),
280 RAT( x, VAT2(iz, 1, lev)),
282 &itmax_s, &iters_s, &errtol_s, omega,
283 &iresid, &iadjoint, &mgsmoo_s);
291 RAT(ipc, VAT2(iz, 5, lev)),
292 RAT(rpc, VAT2(iz, 6, lev)),
293 RAT( ac, VAT2(iz, 7, lev)),
294 RAT( cc, VAT2(iz, 1, lev)),
295 RAT( fc, VAT2(iz, 1, lev)),
296 RAT( x, VAT2(iz, 1, lev)),
298 rsnrm =
Vxnrm1(&nxf, &nyf, &nzf, w1);
299 }
else if (*
istop == 1) {
301 RAT(ipc, VAT2(iz, 5, lev)),
302 RAT(rpc, VAT2(iz, 6, lev)),
303 RAT( ac, VAT2(iz, 7, lev)),
304 RAT( cc, VAT2(iz, 1, lev)),
305 RAT( fc, VAT2(iz, 1, lev)),
306 RAT( x, VAT2(iz, 1, lev)),
308 rsnrm =
Vxnrm1(&nxf, &nyf, &nzf, w1);
309 }
else if (*
istop == 2) {
311 RAT(tru, VAT2(iz, 1, lev)),
315 &alpha, RAT(x, VAT2(iz, 1, lev)), w1);
316 rsnrm =
Vxnrm1(&nxf, &nyf, &nzf, w1);
318 RAT(x, VAT2(iz, 1, lev)),
319 RAT(tru, VAT2(iz, 1, lev)));
320 }
else if (*
istop == 3) {
322 RAT(tru, VAT2(iz, 1, lev)),
326 &alpha, RAT(x, VAT2(iz, 1, lev)), w1);
327 rsnrm =
Vxnrm2(&nxf, &nyf, &nzf, w1);
328 }
else if (*
istop == 4) {
330 RAT(tru, VAT2(iz, 1, lev)),
334 &alpha, RAT(x, VAT2(iz, 1, lev)), w1);
335 rsnrm =
Vxnrm2(&nxf, &nyf, &nzf, w1);
336 }
else if (*
istop == 5) {
338 RAT(tru, VAT2(iz, 1, lev)),
342 &alpha, RAT(x, VAT2(iz, 1, lev)), w1);
344 RAT(ipc, VAT2(iz, 5, lev)),
345 RAT(rpc, VAT2(iz, 6, lev)),
346 RAT( ac, VAT2(iz, 7, lev)),
347 RAT( cc, VAT2(iz, 1, lev)),
349 rsnrm = VSQRT(
Vxdot(&nxf, &nyf, &nzf, w1, w2));
351 Vnm_print(2,
"Vmvcs: bad istop value: %d\n", *
istop);
354 Vprtstp(*iok, *iters, rsnrm, rsden, orsnrm);
369 lev = (*ilev - 1) + level;
376 nuuu = Vivariv(
nu1, &lev);
378 RAT(ipc, VAT2(iz, 5, lev)), RAT(rpc, VAT2(iz, 6, lev)),
379 RAT( ac, VAT2(iz, 7, lev)), RAT( cc, VAT2(iz, 1, lev)),
380 RAT( fc, VAT2(iz, 1, lev)), RAT( x, VAT2(iz, 1, lev)),
382 &nuuu, &iters_s, &errtol_s, omega,
383 &iresid, &iadjoint,
mgsmoo);
384 Vxcopy(&nxf, &nyf, &nzf, w1, RAT(w0, VAT2(iz, 1, lev)));
391 for (level = 2; level <= *
nlev; level++) {
392 lev = (*ilev - 1) + level;
396 Vmkcors(&numlev, &nxf, &nyf, &nzf, &
nxc, &
nyc, &
nzc);
400 w1, RAT(w0, VAT2(iz, 1, lev)),
401 RAT(pc, VAT2(iz, 11, lev-1)));
405 RAT(x, VAT2(iz, 1, lev-1)), RAT(w4, VAT2(iz, 1, lev)));
414 RAT(ipc, VAT2(iz, 5, lev)), RAT(rpc, VAT2(iz, 6, lev)),
415 RAT( ac, VAT2(iz, 7, lev)), RAT( cc, VAT2(iz, 1, lev)),
416 RAT( w4, VAT2(iz, 1, lev)), RAT( fc, VAT2(iz, 1, lev)),
421 Vxaxpy(&nxf, &nyf, &nzf, &alpha,
422 RAT(w0, VAT2(iz, 1, lev)), RAT(fc, VAT2(iz, 1, lev)));
425 if (level != *
nlev) {
429 RAT(w4, VAT2(iz, 1, lev)), RAT(x, VAT2(iz, 1, lev)));
434 nuuu = Vivariv (
nu1,&lev);
436 RAT(ipc, VAT2(iz, 5, lev)), RAT(rpc, VAT2(iz, 6, lev)),
437 RAT( ac, VAT2(iz, 7, lev)), RAT( cc, VAT2(iz, 1, lev)),
438 RAT( fc, VAT2(iz, 1, lev)),
439 RAT( x, VAT2(iz, 1, lev)),
441 &nuuu, &iters_s, &errtol_s, omega,
442 &iresid, &iadjoint,
mgsmoo);
456 lev = (*ilev - 1) + level;
465 Vazeros(&nxf, &nyf, &nzf, RAT(x, VAT2(iz, 1, lev)));
467 RAT(ipc, VAT2(iz, 5, lev)), RAT(rpc, VAT2(iz, 6, lev)),
468 RAT( ac, VAT2(iz, 7, lev)), RAT( cc, VAT2(iz, 1, lev)),
469 RAT( fc, VAT2(iz, 1, lev)),
470 RAT( x, VAT2(iz, 1, lev)),
472 &itmax_s, &iters_s, &errtol_s, omega,
473 &iresid, &iadjoint, &mgsmoo_s);
480 for (level = *
nlev - 1; level >= 1; level--) {
481 lev = (*ilev - 1) + level;
485 Vmkfine(&numlev, &nxf, &nyf, &nzf, &
nxc, &
nyc, &
nzc);
489 Vxaxpy(&nxf, &nyf, &nzf, &alpha,
490 RAT(w4, VAT2(iz, 1, lev + 1)), RAT(x, VAT2(iz, 1, lev + 1)));
493 Vlinesearch(&nxf, &nyf, &nzf, &xdamp,
494 RAT(ipc, VAT2(iz, 5, lev + 1)),
495 RAT(rpc, VAT2(iz, 6, lev + 1)),
496 RAT( ac, VAT2(iz, 7, lev + 1)),
497 RAT( cc, VAT2(iz, 1, lev + 1)),
498 RAT( fc, VAT2(iz, 1, lev + 1)),
499 RAT( x, VAT2(iz, 1, lev + 1)),
500 RAT( w4, VAT2(iz, 1, lev + 1)),
501 RAT( w0, VAT2(iz, 1, lev + 1)),
506 RAT(x, VAT2(iz, 1, lev + 1)),
508 RAT(pc, VAT2(iz, 11, lev)));
516 Vxaxpy(&nxf, &nyf, &nzf, &xdamp, w1, RAT(x, VAT2(iz, 1, lev)));
523 nuuu = Vivariv(
nu2, &lev);
525 RAT(ipc, VAT2(iz, 5, lev)), RAT(rpc, VAT2(iz, 6, lev)),
526 RAT( ac, VAT2(iz, 7, lev)), RAT( cc, VAT2(iz, 1, lev)),
527 RAT( fc, VAT2(iz, 1, lev)),
528 RAT( x, VAT2(iz, 1, lev)),
530 &nuuu, &iters_s, &errtol_s, omega,
531 &iresid, &iadjoint,
mgsmoo);
547 RAT(ipc, VAT2(iz, 5, lev)), RAT(rpc, VAT2(iz, 6, lev)),
548 RAT( ac, VAT2(iz, 7, lev)), RAT( cc, VAT2(iz, 1, lev)),
549 RAT( fc, VAT2(iz, 1, lev)), RAT( x, VAT2(iz, 1, lev)),
551 rsnrm =
Vxnrm1(&nxf, &nyf, &nzf, w1);
552 }
else if (*
istop == 1) {
554 RAT(ipc, VAT2(iz, 5, lev)), RAT(rpc, VAT2(iz, 6, lev)),
555 RAT( ac, VAT2(iz, 7, lev)), RAT( cc, VAT2(iz, 1, lev)),
556 RAT( fc, VAT2(iz, 1, lev)), RAT( x, VAT2(iz, 1, lev)),
558 rsnrm =
Vxnrm1(&nxf, &nyf, &nzf,w1);
559 }
else if (*
istop == 2) {
560 Vxcopy(&nxf, &nyf, &nzf, RAT(tru, VAT2(iz, 1, lev)), w1);
562 Vxaxpy(&nxf, &nyf, &nzf, &alpha, RAT(x, VAT2(iz, 1, lev)), w1);
563 rsnrm =
Vxnrm1(&nxf, &nyf, &nzf, w1);
565 RAT( x, VAT2(iz, 1, lev)),
566 RAT(tru, VAT2(iz, 1, lev)));
567 }
else if (*
istop == 3) {
568 Vxcopy(&nxf, &nyf, &nzf, RAT(tru, VAT2(iz, 1, lev)), w1);
570 Vxaxpy(&nxf, &nyf, &nzf, &alpha, RAT(x, VAT2(iz, 1, lev)), w1);
571 rsnrm =
Vxnrm2(&nxf, &nyf, &nzf,w1);
572 }
else if (*
istop == 4) {
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 == 5) {
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);
582 RAT(ipc, VAT2(iz, 5, lev)), RAT(rpc, VAT2(iz, 6, lev)),
583 RAT( ac, VAT2(iz, 7, lev)), RAT( cc, VAT2(iz, 1, lev)),
585 rsnrm = VSQRT(
Vxdot(&nxf, &nyf, &nzf,w1,w2));
587 VABORT_MSG1(
"Bad istop value: %d", *
istop);
590 Vprtstp(*iok, *iters, rsnrm, rsden, orsnrm);
592 if ((rsnrm / rsden) <= *
errtol)
596 if (*iters >= *
itmax) {
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 Vfmvfas(int *nx, int *ny, int *nz, double *x, int *iz, double *w0, double *w1, double *w2, double *w3, double *w4, int *istop, int *itmax, int *iters, int *ierror, int *nlev, int *ilev, int *nlev_real, int *mgsolv, int *iok, int *iinfo, double *epsiln, double *errtol, double *omega, int *nu1, int *nu2, int *mgsmoo, int *ipc, double *rpc, double *pc, double *ac, double *cc, double *fc, double *tru)
Multigrid nonlinear solve iteration routine.
VPUBLIC void Vazeros(int *nx, int *ny, int *nz, double *x)
Zero out operation for a grid function, including boundary values.
VPUBLIC void Vxcopy(int *nx, int *ny, int *nz, double *x, double *y)
A collection of useful low-level routines (timing, etc).
VPUBLIC double Vxdot(int *nx, int *ny, int *nz, double *x, double *y)
Inner product operation for a grid function with boundary values.
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 Vmvfas(int *nx, int *ny, int *nz, double *x, int *iz, double *w0, double *w1, double *w2, double *w3, double *w4, int *istop, int *itmax, int *iters, int *ierror, int *nlev, int *ilev, int *nlev_real, int *mgsolv, int *iok, int *iinfo, double *epsiln, double *errtol, double *omega, int *nu1, int *nu2, int *mgsmoo, int *ipc, double *rpc, double *pc, double *ac, double *cc, double *fc, double *tru)
Nonlinear multilevel method.
VPUBLIC void Vxaxpy(int *nx, int *ny, int *nz, double *alpha, double *x, double *y)
saxpy operation for a grid function with boundary values.
VPUBLIC double Vxnrm2(int *nx, int *ny, int *nz, double *x)
Norm operation for a grid function with boundary values.
VPUBLIC double Vxnrm1(int *nx, int *ny, int *nz, double *x)
Norm operation for a grid function with boundary values.
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.
VEXTERNC void Vnsmooth(int *nx, int *ny, int *nz, int *ipc, double *rpc, double *ac, double *cc, double *fc, double *x, double *w1, double *w2, double *r, int *itmax, int *iters, double *errtol, double *omega, int *iresid, int *iadjoint, int *meth)
call the appropriate non-linear smoothing 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.