APBS 3.0.0
Loading...
Searching...
No Matches
newtond.c
1
56#include "newtond.h"
57
58VPUBLIC 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
162VPUBLIC 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
550VPUBLIC 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}
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:165
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
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:173
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
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:341
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:57
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:915
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:57
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:112
VPUBLIC double Vxnrm2(int *nx, int *ny, int *nz, double *x)
Norm operation for a grid function with boundary values.
Definition mikpckd.c:152
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:598
VPUBLIC void Vazeros(int *nx, int *ny, int *nz, double *x)
Zero out operation for a grid function, including boundary values.
Definition mikpckd.c:195
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:57
VPUBLIC double Vxnrm1(int *nx, int *ny, int *nz, double *x)
Norm operation for a grid function with boundary values.
Definition mikpckd.c:131
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:782
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
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:232