APBS  1.5
newtond.c
1 
56 #include "newtond.h"
57 
58 VPUBLIC void Vfnewton(int *nx,int *ny,int *nz,
59  double *x, int *iz,
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,
63  int *mgsolv, int *iok, int *iinfo,
64  double *epsiln, double *errtol, double *omega,
65  int *nu1, int *nu2, int *mgsmoo,
66  double *cprime, double *rhs, double *xtmp,
67  int *ipc, double *rpc,
68  double *pc, double *ac, double *cc, double *fc, double *tru) {
69 
70  int level, itmxd, nlevd, iterd, iokd;
71  int nxf, nyf, nzf;
72  int nxc, nyc, nzc;
73  int istpd;
74  int numlev;
75  double errd;
76 
77  MAT2( iz, 50, 1);
78 
79  // Recover gridsizes ***
80  nxf = *nx;
81  nyf = *ny;
82  nzf = *nz;
83 
84  numlev = *nlev - 1;
85  Vmkcors(&numlev, &nxf, &nyf, &nzf, &nxc, &nyc, &nzc);
86 
87  // Move up grids: interpolate solution to finer, do newton
88  if (*iinfo > 1) {
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);
92  }
93 
94  for (level=*nlev_real; level<=*ilev+1; level--) {
95 
96  // Call mv cycle
97  errd = *errtol;
98  itmxd = 1000;
99  nlevd = *nlev_real - level + 1;
100  iterd = 0;
101  iokd = *iok;
102  istpd = *istop;
103 
104  Vnewton(&nxc, &nyc, &nzc,
105  x, iz,
106  w0, w1, w2, w3,
107  &istpd, &itmxd, &iterd, ierror,
108  &nlevd, &level, nlev_real,
109  mgsolv, &iokd, iinfo,
110  epsiln, &errd, omega,
111  nu1, nu2, mgsmoo,
112  cprime, rhs, xtmp,
113  ipc, rpc,
114  pc, ac, cc, fc, tru);
115 
116 
117  // Find new grid size ***
118  numlev = 1;
119  Vmkfine(&numlev, &nxc, &nyc, &nzc, &nxf, &nyf, &nzf);
120 
121  // Interpolate to next finer grid (use correct bc's)
122  VinterpPMG(&nxc, &nyc, &nzc,
123  &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)));
127 
128  /*
129  Commented out fortran code. May need to implement later
130  call ninterpPMG(nxc,nyc,nzc,nxf,nyf,nzf,
131  x(iz(1,level)),x(iz(1,level-1)),pc(iz(11,level-1)),
132  ipc(iz(5,level-1)),rpc(iz(6,level-1)),
133  ac(iz(7,level-1)),cc(iz(1,level-1)),fc(iz(1,level-1)))
134  */
135 
136  // New grid size
137  nxc = nxf;
138  nyc = nyf;
139  nzc = nzf;
140  }
141 
142 
143 
144  // Call mv cycle
145  level = *ilev;
146 
147  Vnewton(nx, ny, nz,
148  x, iz,
149  w0, w1, w2, w3,
150  istop, itmax, iters, ierror,
151  nlev, &level, nlev_real,
152  mgsolv, iok, iinfo,
153  epsiln, errtol, omega,
154  nu1, nu2, mgsmoo,
155  cprime, rhs, xtmp,
156  ipc, rpc,
157  pc, ac, cc, fc, tru);
158 }
159 
160 
161 
162 VPUBLIC void Vnewton(int *nx, int *ny, int *nz,
163  double *x, int *iz,
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) {
173 
174  int level, lev;
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;
178 
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;
183 
184  // Utility and temproary parameters
185  double alpha;
186 
187  MAT2(iz, 50, 1);
188 
189  // Recover level information
190  level = 1;
191  lev = (*ilev - 1) + level;
192 
193  // Do some i/o if requested
194  if (*iinfo > 1) {
195  VMESSAGE3("Starting: (%d, %d, %d)", *nx, *ny, *nz);
196  }
197 
198  if (*iok != 0) {
199  Vprtstp(*iok, -1, 0.0, 0.0, 0.0);
200  }
201 
202  /**************************************************************
203  *** note: if (iok!=0) then: use a stopping test. ***
204  *** else: use just the itmax to stop iteration. ***
205  **************************************************************
206  *** istop=0 most efficient (whatever it is) ***
207  *** istop=1 relative residual ***
208  *** istop=2 rms difference of successive iterates ***
209  *** istop=3 relative true error (provided for testing) ***
210  **************************************************************/
211 
212  // Compute denominator for stopping criterion
213  if (*istop == 0) {
214  rsden = 1.0;
215  } else if (*istop == 1) {
216 
217  // Compute initial residual with zero initial guess
218  // this is analogous to the linear case where one can
219  // simply take norm of rhs for a zero initial guess
220 
221  Vazeros(nx, ny, nz, w1);
222 
223  Vnmresid(nx, ny, nz,
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)),
227  w1, w2, w3);
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) {
236  Vnmatvec(nx, ny, nz,
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)),
240  w1, w2);
241  rsden = VSQRT(Vxdot(nx, ny, nz, RAT(tru, VAT2(iz, 1, lev)), w1));
242  } else {
243  VABORT_MSG1("Bad istop value: %d\n", *istop);
244  }
245 
246  if (rsden == 0.0) {
247  rsden = 1.0;
248  VWARN_MSG0(rsden != 0, "rhs is zero");
249  }
250  rsnrm = rsden;
251  orsnrm = rsnrm;
252 
253  if (*iok != 0) {
254  Vprtstp(*iok, 0, rsnrm, rsden, orsnrm);
255  }
256 
257  /*********************************************************************
258  *** begin newton iteration
259  *********************************************************************/
260 
261  // Now compute residual with the initial guess
262 
263  Vnmresid(nx, ny, nz,
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)),
267  w0, w2);
268  xnorm_old = Vxnrm1(nx, ny, nz, w0);
269  if (*iok != 0) {
270  xnorm_den = rsden;
271  } else {
272  xnorm_den = xnorm_old;
273  }
274 
275 
276 
277  /*********************************************************************
278  *** begin the loop
279  *********************************************************************/
280 
281  // Setup for the looping
282  VMESSAGE0("Damping enabled");
283  idamp = 1;
284  *iters = 0;
285 
286  //30
287  while(1) {
288 
289  (*iters)++;
290 
291  // Save iterate if stop test will use it on next iter
292  if (*istop == 2) {
293  Vxcopy(nx, ny, nz,
294  RAT(x, VAT2(iz, 1, lev)), RAT(tru, VAT2(iz, 1, lev)));
295  }
296 
297  // Compute the current jacobian system and rhs
298  ipkey = VAT(ipc, 10);
299  Vgetjac(nx, ny, nz, nlev_real, iz, ilev, &ipkey,
300  x, w0, cprime, rhs, cc, pc);
301 
302  // Determine number of correct digits in current residual
303  // Algorithm 5.3 in the thesis, test version (1')
304  // Global-superlinear convergence
305  bigc = 1.0;
306  ord = 2.0;
307 
308  /* NAB 06-18-01: If complex problems are not converging, set this to
309  * machine epsilon. This makes it use the exact jacobian rather than
310  * the appropriate form (as here)
311  */
312  errtol_s = VMIN2((0.9 * xnorm_old), (bigc * VPOW(xnorm_old, ord)));
313  VMESSAGE1("Using errtol_s: %f", errtol_s);
314 
315  // Do a linear multigrid solve of the newton equations
316  Vazeros(nx, ny, nz, RAT(xtmp, VAT2(iz, 1, lev)));
317 
318  itmax_s = 1000;
319  istop_s = 0;
320  iters_s = 0;
321  ierror_s = 0;
322 
323  // NAB 06-18-01 -- What this used to be:
324  iok_s = 0;
325  iinfo_s = 0;
326  if ((*iinfo >= 2) && (*ilev == 1))
327  iok_s = 2;
328 
329  // What it's changed to:
330  if (*iinfo >= 2)
331  iinfo_s = *iinfo;
332  iok_s = 2;
333 
334  // End of NAB hack.
335 
336  Vmvcs(nx, ny, nz,
337  xtmp, iz,
338  w0, w1, w2, w3,
339  &istop_s, &itmax_s, &iters_s, &ierror_s,
340  nlev, ilev, nlev_real, mgsolv,
341  &iok_s, &iinfo_s,
342  epsiln, &errtol_s, omega,
343  nu1, nu2, mgsmoo,
344  ipc, rpc, pc, ac, cprime, rhs, tru);
345 
346  /**************************************************************
347  *** note: rhs and cprime are now available as temp vectors ***
348  **************************************************************/
349 
350  // If damping is still enabled -- doit
351  if (idamp == 1) {
352 
353  // Try the correction
354  Vxcopy(nx, ny, nz,
355  RAT(x, VAT2(iz, 1, lev)), w1);
356  damp = 1.0;
357  Vxaxpy(nx, ny, nz, &damp, RAT(xtmp, VAT2(iz, 1, lev)), w1);
358 
359  Vnmresid(nx, ny, nz,
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)),
363  w1, w0,
364  RAT(rhs, VAT2(iz, 1, lev)));
365  xnorm_new = Vxnrm1(nx, ny, nz, w0);
366 
367  // Damping is still enabled -- doit
368  damp = 1.0;
369  iter_d = 0;
370  itmax_d = 10;
371  mode = 0;
372 
373  VMESSAGE1("Attempting damping, relres = %f", xnorm_new / xnorm_den);
374 
375  while(iter_d < itmax_d) {
376  if (mode == 0) {
377  if (xnorm_new < xnorm_old) {
378  mode = 1;
379  }
380  } else if (xnorm_new > xnorm_med) {
381  break;
382  }
383 
384  // Keep old soln and residual around, and its norm
385  Vxcopy(nx, ny, nz, w1, w2);
386  Vxcopy(nx, ny, nz, w0, w3);
387  xnorm_med = xnorm_new;
388 
389  // New damped correction, residual, and its norm
390  Vxcopy(nx, ny, nz,
391  RAT(x, VAT2(iz, 1, lev)), w1);
392  damp = damp / 2.0;
393  Vxaxpy(nx, ny, nz, &damp, RAT(xtmp, VAT2(iz, 1, lev)), w1);
394 
395  Vnmresid(nx, ny, nz,
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)),
399  w1, w0,
400  RAT(rhs, VAT2(iz, 1, lev)));
401  xnorm_new = Vxnrm1(nx, ny, nz, w0);
402 
403  // Next iter...
404  iter_d = iter_d + 1;
405  VMESSAGE1("Attempting damping, relres = %f",
406  xnorm_new / xnorm_den);
407 
408  }
409 
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;
414 
415  VMESSAGE1("Damping accepted, relres = %f", xnorm_new / xnorm_den);
416 
417  // Determine whether or not to disable damping
418  if ((iter_d - 1) == 0) {
419  VMESSAGE0("Damping disabled");
420  idamp = 0;
421  }
422  } else {
423 
424  // Damping is disabled -- accept the newton step
425  damp = 1.0;
426 
427  Vxaxpy(nx, ny, nz, &damp,
428  RAT(xtmp, VAT2(iz, 1, lev)), RAT(x, VAT2(iz, 1, lev)));
429 
430  Vnmresid(nx, ny, nz,
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)),
434  w0,
435  RAT(rhs, VAT2(iz, 1, lev)));
436 
437  xnorm_new = Vxnrm1(nx, ny, nz, w0);
438  xnorm_old = xnorm_new;
439  }
440 
441  // Compute/check the current stopping test ***
442  if (iok != 0) {
443 
444  orsnrm = rsnrm;
445 
446  if (*istop == 0) {
447  rsnrm = xnorm_new;
448  } else if (*istop == 1) {
449  rsnrm = xnorm_new;
450  } else if (*istop == 2) {
451  Vxcopy(nx, ny, nz, RAT(tru, VAT2(iz, 1, lev)), w1);
452  alpha = -1.0;
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);
457  alpha = -1.0;
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);
462  alpha = -1.0;
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);
467  alpha = -1.0;
468  Vxaxpy(nx, ny, nz, &alpha, RAT(x, VAT2(iz, 1, lev)), w1);
469  Vnmatvec(nx, ny, nz,
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)),
472  w1, w2, w3);
473  rsnrm = VSQRT(Vxdot(nx, ny, nz, w1, w2));
474  } else {
475  VABORT_MSG1("Bad istop value: %d", *istop);
476  }
477 
478  Vprtstp(*iok, *iters, rsnrm, rsden, orsnrm);
479 
480  if ((rsnrm/rsden) <= *errtol)
481  break;
482  }
483 
484  // Check iteration count ***
485  if (*iters >= *itmax)
486  break;
487  }
488 
489  // Condition estimate of final jacobian
490  if (*iinfo > 2) {
491 
492  Vnm_print(2, "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n");
493  Vnm_print(2, "% Vnewton: JACOBIAN ANALYSIS ==> (%d, %d, %d)\n",
494  *nx, *ny, *nz );
495 
496  // Largest eigenvalue of the jacobian matrix
497  Vnm_print(2, "% Vnewton: Power calculating rho(JAC)\n");
498  itmax_p = 1000;
499  errtol_p = 1.0e-4;
500  iters_p = 0;
501  iinfo_p = *iinfo;
502 
503  Vpower(nx, ny, nz,
504  iz, ilev,
505  ipc, rpc, ac, cprime,
506  w0, w1, w2, w3,
507  &rho_max, &rho_max_mod,
508  &errtol_p, &itmax_p, &iters_p, &iinfo_p);
509 
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);
513 
514  // Smallest eigenvalue of the system matrix A ***
515  Vnm_print(2, "% Vnewton: ipower calculating lambda_min(JAC)...\n");
516  itmax_p = 1000;
517  errtol_p = 1.0e-4;
518  iters_p = 0;
519  iinfo_p = *iinfo;
520 
521  Vazeros(nx, ny, nz, xtmp);
522 
523  Vipower(nx, ny, nz,
524  xtmp, iz,
525  w0, w1, w2, w3,
526  rhs, &rho_min, &rho_min_mod,
527  &errtol_p, &itmax_p, &iters_p,
528  nlev, ilev, nlev_real, mgsolv,
529  &iok_p, &iinfo_p,
530  epsiln, errtol, omega,
531  nu1, nu2, mgsmoo,
532  ipc, rpc,
533  pc, ac, cprime, tru);
534 
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);
538 
539  // Condition number estimate
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");
545  }
546 }
547 
548 
549 
550 VPUBLIC void Vgetjac(int *nx, int *ny, int *nz,
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) {
555 
556  int nxx, nyy, nzz;
557  int nxold, nyold, nzold;
558  int level, numlev;
559 
560  MAT2(iz, 50, 1);
561 
562  // Setup
563  nxx = *nx;
564  nyy = *ny;
565  nzz = *nz;
566 
567  // Form the rhs of the newton system -- just current residual
568  Vxcopy(nx, ny, nz, r, RAT(rhs, VAT2(iz, 1,*lev)));
569 
570  // Get nonlinear part of the jacobian operator
571  Vdc_vec(RAT(cc, VAT2(iz, 1,*lev)), RAT(x, VAT2(iz, 1,*lev)),
572  RAT(cprime, VAT2(iz, 1,*lev)),
573  nx, ny, nz, ipkey);
574 
575  // Build the (nlev-1) level operators
576  for (level=*lev+1; level<=*nlev_real; level++) {
577  nxold = nxx;
578  nyold = nyy;
579  nzold = nzz;
580 
581  numlev = 1;
582  Vmkcors(&numlev, &nxold, &nyold, &nzold, &nxx, &nyy, &nzz);
583 
584  // Make the coarse grid rhs functions
585  Vrestrc(&nxold, &nyold, &nzold,
586  &nxx, &nyy, &nzz,
587  RAT(rhs, VAT2(iz, 1,level-1)),
588  RAT(rhs, VAT2(iz, 1,level )),
589  RAT( pc, VAT2(iz, 11,level-1)));
590 
591  // Make the coarse grid helmholtz terms
592  Vrestrc(&nxold, &nyold, &nzold,
593  &nxx, &nyy, &nzz,
594  RAT(cprime, VAT2(iz, 1, level-1)),
595  RAT(cprime, VAT2(iz, 1, level )),
596  RAT( pc, VAT2(iz, 11, level-1)));
597  }
598 }
599 
int mgsolv
Definition: vpmgp.h:174
int iinfo
Definition: vpmgp.h:130
int nzc
Definition: vpmgp.h:98
int nxc
Definition: vpmgp.h:96
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.
Definition: matvecd.c:910
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.
Definition: mgcsd.c:52
int ipkey
Definition: vpmgp.h:109
VPUBLIC void Vazeros(int *nx, int *ny, int *nz, double *x)
Zero out operation for a grid function, including boundary values.
Definition: mikpckd.c:190
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)
Definition: mypdec.c:336
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.
Definition: powerd.c:52
int mgsmoo
Definition: vpmgp.h:160
VPUBLIC void Vxcopy(int *nx, int *ny, int *nz, double *x, double *y)
A collection of useful low-level routines (timing, etc).
Definition: mikpckd.c:52
VPUBLIC double Vxdot(int *nx, int *ny, int *nz, double *x, double *y)
Inner product operation for a grid function with boundary values.
Definition: mikpckd.c:168
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.
Definition: powerd.c:160
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.
Definition: matvecd.c:227
int nu1
Definition: vpmgp.h:158
VPUBLIC void Vxaxpy(int *nx, int *ny, int *nz, double *alpha, double *x, double *y)
saxpy operation for a grid function with boundary values.
Definition: mikpckd.c:107
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.
Definition: newtond.c:162
int nyc
Definition: vpmgp.h:97
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.
Definition: newtond.c:550
int istop
Definition: vpmgp.h:123
int ny
Definition: vpmgp.h:84
int nx
Definition: vpmgp.h:83
VPUBLIC double Vxnrm2(int *nx, int *ny, int *nz, double *x)
Norm operation for a grid function with boundary values.
Definition: mikpckd.c:147
VPUBLIC double Vxnrm1(int *nx, int *ny, int *nz, double *x)
Norm operation for a grid function with boundary values.
Definition: mikpckd.c:126
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.
Definition: matvecd.c:593
int nu2
Definition: vpmgp.h:159
int itmax
Definition: vpmgp.h:122
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.
Definition: newtond.c:58
int nlev
Definition: vpmgp.h:86
double errtol
Definition: vpmgp.h:121
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.
Definition: matvecd.c:777
int nz
Definition: vpmgp.h:85