60 double *w0,
double *w1,
double *w2,
double *w3,
61 int *
istop,
int *
itmax,
int *iters,
int *ierror,
62 int *
nlev,
int *ilev,
int *nlev_real,
64 double *epsiln,
double *
errtol,
double *omega,
66 double *cprime,
double *rhs,
double *xtmp,
67 int *ipc,
double *rpc,
68 double *pc,
double *ac,
double *cc,
double *fc,
double *tru) {
70 int level, itmxd, nlevd, iterd, iokd;
85 Vmkcors(&numlev, &nxf, &nyf, &nzf, &nxc, &nyc, &nzc);
89 VMESSAGE0(
"Starting");
90 VMESSAGE3(
"Fine Grid Size: (%d, %d, %d)", nxf, nyf, nzf);
91 VMESSAGE3(
"Course Grid Size: (%d, %d, %d)", nxc, nyc, nzc);
94 for (level=*nlev_real; level<=*ilev+1; level--) {
99 nlevd = *nlev_real - level + 1;
107 &istpd, &itmxd, &iterd, ierror,
108 &nlevd, &level, nlev_real,
109 mgsolv, &iokd, iinfo,
110 epsiln, &errd, omega,
114 pc, ac, cc, fc, tru);
119 Vmkfine(&numlev, &nxc, &nyc, &nzc, &nxf, &nyf, &nzf);
124 RAT( x, VAT2(iz, 1, level )),
125 RAT( x, VAT2(iz, 1, level-1)),
126 RAT(pc, VAT2(iz, 11, level-1)));
150 istop, itmax, iters, ierror,
151 nlev, &level, nlev_real,
153 epsiln, errtol, omega,
157 pc, ac, cc, fc, tru);
164 double *w0,
double *w1,
double *w2,
double *w3,
165 int *
istop,
int *
itmax,
int *iters,
int *ierror,
166 int *
nlev,
int *ilev,
int *nlev_real,
168 double *epsiln,
double *
errtol,
double *omega,
170 double *cprime,
double *rhs,
double *xtmp,
171 int *ipc,
double *rpc,
172 double *pc,
double *ac,
double *cc,
double *fc,
double *tru) {
175 int itmax_s, iters_s, ierror_s, iok_s, iinfo_s, istop_s;
176 double errtol_s, ord, bigc;
177 double rsden, rsnrm, orsnrm;
179 double xnorm_old, xnorm_new, damp, xnorm_med, xnorm_den;
180 double rho_max, rho_min, rho_max_mod, rho_min_mod, errtol_p;
181 int iter_d, itmax_d, mode, idamp,
ipkey;
182 int itmax_p, iters_p, iok_p, iinfo_p;
191 lev = (*ilev - 1) + level;
195 VMESSAGE3(
"Starting: (%d, %d, %d)", *nx, *ny, *nz);
199 Vprtstp(*iok, -1, 0.0, 0.0, 0.0);
215 }
else if (*istop == 1) {
224 RAT(ipc, VAT2(iz, 5, lev)), RAT(rpc, VAT2(iz, 6, lev)),
225 RAT( ac, VAT2(iz, 7, lev)), RAT( cc, VAT2(iz, 1, lev)),
226 RAT( fc, VAT2(iz, 1, lev)),
228 rsden =
Vxnrm1(nx, ny, nz, w2);
229 }
else if (*istop == 2) {
230 rsden = VSQRT( *nx * *ny * *nz);
231 }
else if (*istop == 3) {
232 rsden =
Vxnrm2(nx, ny, nz, RAT(tru, VAT2(iz, 1, lev)));
233 }
else if (*istop == 4) {
234 rsden =
Vxnrm2(nx, ny, nz, RAT(tru, VAT2(iz, 1, lev)));
235 }
else if (*istop == 5) {
237 RAT(ipc, VAT2(iz, 5, lev)), RAT(rpc, VAT2(iz, 6, lev)),
238 RAT( ac, VAT2(iz, 7, lev)), RAT( cc, VAT2(iz, 1, lev)),
239 RAT(tru, VAT2(iz, 1, lev)),
241 rsden = VSQRT(
Vxdot(nx, ny, nz, RAT(tru, VAT2(iz, 1, lev)), w1));
243 VABORT_MSG1(
"Bad istop value: %d\n", *istop);
248 VWARN_MSG0(rsden != 0,
"rhs is zero");
254 Vprtstp(*iok, 0, rsnrm, rsden, orsnrm);
264 RAT(ipc, VAT2(iz, 5, lev)), RAT(rpc, VAT2(iz, 6, lev)),
265 RAT( ac, VAT2(iz, 7, lev)), RAT( cc, VAT2(iz, 1, lev)),
266 RAT( fc, VAT2(iz, 1, lev)), RAT( x, VAT2(iz, 1, lev)),
268 xnorm_old =
Vxnrm1(nx, ny, nz, w0);
272 xnorm_den = xnorm_old;
282 VMESSAGE0(
"Damping enabled");
294 RAT(x, VAT2(iz, 1, lev)), RAT(tru, VAT2(iz, 1, lev)));
298 ipkey = VAT(ipc, 10);
299 Vgetjac(nx, ny, nz, nlev_real, iz, ilev, &ipkey,
300 x, w0, cprime, rhs, cc, pc);
312 errtol_s = VMIN2((0.9 * xnorm_old), (bigc * VPOW(xnorm_old, ord)));
313 VMESSAGE1(
"Using errtol_s: %f", errtol_s);
316 Vazeros(nx, ny, nz, RAT(xtmp, VAT2(iz, 1, lev)));
326 if ((*iinfo >= 2) && (*ilev == 1))
339 &istop_s, &itmax_s, &iters_s, &ierror_s,
340 nlev, ilev, nlev_real, mgsolv,
342 epsiln, &errtol_s, omega,
344 ipc, rpc, pc, ac, cprime, rhs, tru);
355 RAT(x, VAT2(iz, 1, lev)), w1);
357 Vxaxpy(nx, ny, nz, &damp, RAT(xtmp, VAT2(iz, 1, lev)), w1);
360 RAT(ipc, VAT2(iz, 5, lev)), RAT(rpc, VAT2(iz, 6, lev)),
361 RAT( ac, VAT2(iz, 7, lev)), RAT( cc, VAT2(iz, 1, lev)),
362 RAT( fc, VAT2(iz, 1, lev)),
364 RAT(rhs, VAT2(iz, 1, lev)));
365 xnorm_new =
Vxnrm1(nx, ny, nz, w0);
373 VMESSAGE1(
"Attempting damping, relres = %f", xnorm_new / xnorm_den);
375 while(iter_d < itmax_d) {
377 if (xnorm_new < xnorm_old) {
380 }
else if (xnorm_new > xnorm_med) {
385 Vxcopy(nx, ny, nz, w1, w2);
386 Vxcopy(nx, ny, nz, w0, w3);
387 xnorm_med = xnorm_new;
391 RAT(x, VAT2(iz, 1, lev)), w1);
393 Vxaxpy(nx, ny, nz, &damp, RAT(xtmp, VAT2(iz, 1, lev)), w1);
396 RAT(ipc, VAT2(iz, 5, lev)), RAT(rpc, VAT2(iz, 6, lev)),
397 RAT( ac, VAT2(iz, 7, lev)), RAT( cc, VAT2(iz, 1, lev)),
398 RAT( fc, VAT2(iz, 1, lev)),
400 RAT(rhs, VAT2(iz, 1, lev)));
401 xnorm_new =
Vxnrm1(nx, ny, nz, w0);
405 VMESSAGE1(
"Attempting damping, relres = %f",
406 xnorm_new / xnorm_den);
410 Vxcopy(nx, ny, nz, w2, RAT(x, VAT2(iz, 1, lev)));
411 Vxcopy(nx, ny, nz, w3, w0);
412 xnorm_new = xnorm_med;
413 xnorm_old = xnorm_new;
415 VMESSAGE1(
"Damping accepted, relres = %f", xnorm_new / xnorm_den);
418 if ((iter_d - 1) == 0) {
419 VMESSAGE0(
"Damping disabled");
428 RAT(xtmp, VAT2(iz, 1, lev)), RAT(x, VAT2(iz, 1, lev)));
431 RAT(ipc, VAT2(iz, 5, lev)), RAT(rpc, VAT2(iz, 6, lev)),
432 RAT( ac, VAT2(iz, 7, lev)), RAT( cc, VAT2(iz, 1, lev)),
433 RAT( fc, VAT2(iz, 1, lev)), RAT( x, VAT2(iz, 1, lev)),
435 RAT(rhs, VAT2(iz, 1, lev)));
437 xnorm_new =
Vxnrm1(nx, ny, nz, w0);
438 xnorm_old = xnorm_new;
448 }
else if (*istop == 1) {
450 }
else if (*istop == 2) {
451 Vxcopy(nx, ny, nz, RAT(tru, VAT2(iz, 1, lev)), w1);
453 Vxaxpy(nx, ny, nz, &alpha, RAT(x, VAT2(iz, 1, lev)), w1);
454 rsnrm =
Vxnrm1(nx, ny, nz, w1);
455 }
else if (*istop == 3) {
456 Vxcopy(nx, ny, nz, RAT(tru, VAT2(iz, 1, lev)), w1);
458 Vxaxpy(nx, ny, nz, &alpha, RAT(x, VAT2(iz, 1, lev)), w1);
459 rsnrm =
Vxnrm2(nx, ny, nz, w1);
460 }
else if (*istop == 4) {
461 Vxcopy(nx, ny, nz, RAT(tru, VAT2(iz, 1, lev)), w1);
463 Vxaxpy(nx, ny, nz, &alpha, RAT(x, VAT2(iz, 1, lev)), w1);
464 rsnrm =
Vxnrm2(nx, ny, nz, w1);
465 }
else if (*istop == 5) {
466 Vxcopy(nx, ny, nz, RAT(tru, VAT2(iz, 1, lev)), w1);
468 Vxaxpy(nx, ny, nz, &alpha, RAT(x, VAT2(iz, 1, lev)), w1);
470 RAT(ipc, VAT2(iz, 5, lev)), RAT(rpc, VAT2(iz, 6, lev)),
471 RAT( ac, VAT2(iz, 7, lev)), RAT( cc, VAT2(iz, 1, lev)),
473 rsnrm = VSQRT(
Vxdot(nx, ny, nz, w1, w2));
475 VABORT_MSG1(
"Bad istop value: %d", *istop);
478 Vprtstp(*iok, *iters, rsnrm, rsden, orsnrm);
480 if ((rsnrm/rsden) <= *errtol)
485 if (*iters >= *itmax)
492 Vnm_print(2,
"%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n");
493 Vnm_print(2,
"% Vnewton: JACOBIAN ANALYSIS ==> (%d, %d, %d)\n",
497 Vnm_print(2,
"% Vnewton: Power calculating rho(JAC)\n");
505 ipc, rpc, ac, cprime,
507 &rho_max, &rho_max_mod,
508 &errtol_p, &itmax_p, &iters_p, &iinfo_p);
510 Vnm_print(2,
"% Vnewton: power iters = %d\n", iters_p);
511 Vnm_print(2,
"% Vnewton: power eigmax = %d\n", rho_max);
512 Vnm_print(2,
"% Vnewton: power (MODEL) = %d\n", rho_max_mod);
515 Vnm_print(2,
"% Vnewton: ipower calculating lambda_min(JAC)...\n");
526 rhs, &rho_min, &rho_min_mod,
527 &errtol_p, &itmax_p, &iters_p,
528 nlev, ilev, nlev_real, mgsolv,
530 epsiln, errtol, omega,
533 pc, ac, cprime, tru);
535 Vnm_print(2,
"% Vnewton: ipower iters = %d\n", iters_p);
536 Vnm_print(2,
"% Vnewton: ipower eigmin = %d\n", rho_min);
537 Vnm_print(2,
"% Vnewton: ipower (MODEL) = %d\n", rho_min_mod);
540 Vnm_print(2,
"% Vnewton: condition number = %f\n",
541 (
double)rho_max / rho_min);
542 Vnm_print(2,
"% Vnewton: condition (MODEL) = %f\n",
543 (
double)rho_max_mod / rho_min_mod);
544 Vnm_print(2,
"%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n");
551 int *nlev_real,
int *iz,
int *lev,
int *
ipkey,
552 double *x,
double *r,
553 double *cprime,
double *rhs,
554 double *cc,
double *pc) {
557 int nxold, nyold, nzold;
568 Vxcopy(nx, ny, nz, r, RAT(rhs, VAT2(iz, 1,*lev)));
571 Vdc_vec(RAT(cc, VAT2(iz, 1,*lev)), RAT(x, VAT2(iz, 1,*lev)),
572 RAT(cprime, VAT2(iz, 1,*lev)),
576 for (level=*lev+1; level<=*nlev_real; level++) {
582 Vmkcors(&numlev, &nxold, &nyold, &nzold, &nxx, &nyy, &nzz);
585 Vrestrc(&nxold, &nyold, &nzold,
587 RAT(rhs, VAT2(iz, 1,level-1)),
588 RAT(rhs, VAT2(iz, 1,level )),
589 RAT( pc, VAT2(iz, 11,level-1)));
592 Vrestrc(&nxold, &nyold, &nzold,
594 RAT(cprime, VAT2(iz, 1, level-1)),
595 RAT(cprime, VAT2(iz, 1, level )),
596 RAT( pc, VAT2(iz, 11, level-1)));
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.
VEXTERNC void Vmvcs(int *nx, int *ny, int *nz, double *x, int *iz, double *w0, double *w1, double *w2, double *w3, 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)
MG helper functions.
VPUBLIC void Vazeros(int *nx, int *ny, int *nz, double *x)
Zero out operation for a grid function, including boundary values.
VPUBLIC void Vdc_vec(double *coef, double *uin, double *uout, int *nx, int *ny, int *nz, int *ipkey)
Define the derivative of the nonlinearity (vector version)
VPUBLIC void Vpower(int *nx, int *ny, int *nz, int *iz, int *ilev, int *ipc, double *rpc, double *ac, double *cc, double *w1, double *w2, double *w3, double *w4, double *eigmax, double *eigmax_model, double *tol, int *itmax, int *iters, int *iinfo)
Power methods for eigenvalue estimation.
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.
VPUBLIC void Vipower(int *nx, int *ny, int *nz, double *u, int *iz, double *w0, double *w1, double *w2, double *w3, double *w4, double *eigmin, double *eigmin_model, double *tol, int *itmax, int *iters, 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 *tru)
Standard inverse power method for minimum eigenvalue estimation.
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 Vxaxpy(int *nx, int *ny, int *nz, double *alpha, double *x, double *y)
saxpy operation for a grid function with boundary values.
VPUBLIC void Vnewton(int *nx, int *ny, int *nz, double *x, int *iz, double *w0, double *w1, double *w2, double *w3, 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, double *cprime, double *rhs, double *xtmp, int *ipc, double *rpc, double *pc, double *ac, double *cc, double *fc, double *tru)
Inexact-newton-multilevel method.
VPUBLIC void Vgetjac(int *nx, int *ny, int *nz, int *nlev_real, int *iz, int *lev, int *ipkey, double *x, double *r, double *cprime, double *rhs, double *cc, double *pc)
Form the jacobian system.
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.
VPUBLIC void Vfnewton(int *nx, int *ny, int *nz, double *x, int *iz, double *w0, double *w1, double *w2, double *w3, 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, double *cprime, double *rhs, double *xtmp, int *ipc, double *rpc, double *pc, double *ac, double *cc, double *fc, double *tru)
Driver routines for the Newton method.
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.