162VPUBLIC
void Vnewton(
int *nx,
int *ny,
int *nz,
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,
167 int *mgsolv,
int *iok,
int *iinfo,
168 double *epsiln,
double *errtol,
double *omega,
169 int *nu1,
int *nu2,
int *mgsmoo,
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");