APBS  1.4.1
mgfasd.c
1 
50 #include "mgfasd.h"
51 
52 VPUBLIC void Vfmvfas(int *nx, int *ny, int *nz,
53  double *x,
54  int *iz,
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,
58  int *mgsolv, int *iok, int *iinfo,
59  double *epsiln, double *errtol, double *omega,
60  int *nu1, int *nu2, int *mgsmoo,
61  int *ipc, double *rpc,
62  double *pc, double *ac, double *cc, double *fc, double *tru) {
63 
64  // Other Declarations
65  int level, itmxd, nlevd, iterd, iokd;
66  int nxf, nyf, nzf;
67  int nxc, nyc, nzc;
68  int istpd, iinfod;
69  double errd;
70 
71  // Utility variables
72  int numlev;
73 
74  MAT2(iz, 50, *nlev);
75 
76  // Recover gridsizes
77  nxf = *nx;
78  nyf = *ny;
79  nzf = *nz;
80 
81  numlev = *nlev - 1;
82  Vmkcors(&numlev, &nxf, &nyf, &nzf, &nxc, &nyc, &nzc);
83 
84  // Move up grids: interpolate solution to finer, do v cycle
85  if (*iinfo != 0) {
86 
87  Vnm_print(2, "Vfmvfas: starting: (%d,%d,%d) (%d,%d,%d)\n",
88  nxf, nyf, nzf, nxc, nyc, nzc);
89  }
90 
91  for (level = *nlev_real; level >= *ilev + 1; level--) {
92 
93  // Call mv cycle
94  errd = 1.0e-5;
95  itmxd = 1;
96  nlevd = *nlev_real - level + 1;
97  iterd = 0;
98  iokd = 2;
99  iinfod = *iinfo;
100  istpd = *istop;
101  if (*iinfo >= 2)
102  iokd = 2;
103 
104  Vmvfas(&nxc, &nyc, &nzc,
105  x,
106  iz,
107  w0, w1, w2, w3, w4,
108  &istpd, &itmxd, &iterd, ierror,
109  &nlevd, &level, nlev_real,
110  mgsolv, &iokd, &iinfod,
111  epsiln, errtol, omega,
112  nu1, nu2, mgsmoo,
113  ipc, rpc,
114  pc, ac, cc, fc, tru);
115 
116  // Find new grid size
117  numlev = 1;
118  Vmkfine(&numlev, &nxc, &nyc, &nzc, &nxf, &nyf, &nzf);
119 
120  // Interpolate to next finer grid
121  VinterpPMG(&nxc, &nyc, &nzc,
122  &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)));
126 
127  // New grid size
128  nxc = nxf;
129  nyc = nyf;
130  nzc = nzf;
131  }
132 
133 
134 
135  // Call mv cycle
136  level = *ilev;
137 
138  Vmvfas(&nxf, &nyf, &nzf,
139  x, iz,
140  w0, w1, w2, w3, w4,
141  istop, itmax, iters,
142  ierror, nlev, &level, nlev_real,
143  mgsolv, iok, iinfo,
144  epsiln, errtol, omega,
145  nu1, nu2, mgsmoo,
146  ipc, rpc,
147  pc, ac, cc, fc, tru);
148 }
149 
150 
151 
152 VPUBLIC void Vmvfas(int *nx, int *ny, int *nz,
153  double *x,
154  int *iz,
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,
158  int *mgsolv, int *iok, int *iinfo,
159  double *epsiln, double *errtol, double *omega,
160  int *nu1, int *nu2, int *mgsmoo,
161  int *ipc, double *rpc,
162  double *pc, double *ac, double *cc, double *fc, double *tru) {
163 
164  // Other declarations
165  int level, lev;
166  int itmax_s, iters_s, nuuu, ivariv, mgsmoo_s, iresid;
167  int nxf, nyf, nzf;
168  int nxc, nyc, nzc;
169  int iadjoint;
170  double errtol_s;
171  double rsden, rsnrm, orsnrm;
172  double xdamp;
173 
174  int numlev;
175  double alpha;
176 
177  MAT2(iz, 50, *nlev);
178 
179  WARN_UNTESTED;
180 
181  // Recover level information
182  level = 1;
183  lev = (*ilev - 1) + level;
184 
185  // Recover gridsizes
186  nxf = *nx;
187  nyf = *ny;
188  nzf = *nz;
189 
190  numlev = *nlev - 1;
191  Vmkcors(&numlev, &nxf, &nyf, &nzf, &nxc, &nyc, &nzc);
192 
193  // Do some i/o if requested
194  if (*iinfo != 0) {
195  Vnm_print(2, "Vmvfas: starting: (%d, %d, %d) (%d, %d, %d)\n",
196  nxf, nyf, nzf, nxc, nyc, nzc);
197  }
198 
199  // Initial wall clock
200  if (*iok != 0) {
201  Vprtstp(*iok, -1, 0.0, 0.0, 0.0);
202  }
203 
204  /**************************************************************
205  *** note: if (iok != 0) then: use a stopping test. ***
206  *** else: use just the itmax to stop iteration. ***
207  **************************************************************
208  *** istop=0 most efficient (whatever it is) ***
209  *** istop=1 relative residual ***
210  *** istop=2 rms difference of successive iterates ***
211  *** istop=3 relative true error (provided for testing) ***
212  **************************************************************/
213 
214  // Compute denominator for stopping criterion
215  if (*iok != 0) {
216  if (*istop == 0) {
217  rsden = 1.0;
218  } else if (*istop == 1) {
219 
220  // Compute initial residual with zero initial guess
221  // this is analogous to the linear case where one can
222  // simply take norm of rhs for a zero initial guess
223  Vazeros(&nxf, &nyf, &nzf, w1);
224 
225  Vnmresid(&nxf, &nyf, &nzf,
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)),
229  w1, w2, w3);
230 
231  rsden = Vxnrm1(&nxf, &nyf, &nzf, w2);
232 
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) {
240  Vnmatvec(&nxf, &nyf, &nzf,
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)),
244  w1, w2);
245  rsden = VSQRT(Vxdot(&nxf, &nyf, &nzf,RAT(tru, VAT2(iz, 1, lev)), w1));
246  } else {
247  Vnm_print(2, "Vmvfas: bad istop value: %d\n", *istop);
248  }
249  if (rsden == 0.0) {
250  rsden = 1.0;
251  Vnm_print(2, "Vmfas: rhs is zero on finest level\n");
252  }
253  rsnrm = rsden;
254  orsnrm = rsnrm;
255 
256  Vprtstp(*iok, 0, rsnrm, rsden, orsnrm);
257  }
258 
259 
260 
261  /********************************************************************
262  *** solve directly if nlev = 1
263  ********************************************************************/
264 
265  // Solve directly if on the coarse grid
266  if (*nlev == 1) {
267 
268  //solve with ncghs, mgsmoo_s=4 (no residual)
269  iresid = 0;
270  iadjoint = 0;
271  itmax_s = 100;
272  iters_s = 0;
273  errtol_s = *epsiln;
274  mgsmoo_s = *mgsmoo;
275  Vazeros(&nxf, &nyf, &nzf,RAT(x, VAT2(iz, 1, lev)));
276  Vnsmooth(&nxf, &nyf, &nzf,
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)),
281  w1, w2, w3,
282  &itmax_s, &iters_s, &errtol_s, omega,
283  &iresid, &iadjoint, &mgsmoo_s);
284 
285  // Compute the stopping test
286  *iters = 1;
287  if (*iok != 0) {
288  orsnrm = rsnrm;
289  if (*istop == 0) {
290  Vnmresid(&nxf, &nyf, &nzf,
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)),
297  w1, w2);
298  rsnrm = Vxnrm1(&nxf, &nyf, &nzf, w1);
299  } else if (*istop == 1) {
300  Vnmresid(&nxf, &nyf, &nzf,
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)),
307  w1, w2);
308  rsnrm = Vxnrm1(&nxf, &nyf, &nzf, w1);
309  } else if (*istop == 2) {
310  Vxcopy(&nxf, &nyf, &nzf,
311  RAT(tru, VAT2(iz, 1, lev)),
312  w1);
313  alpha = -1.0;
314  Vxaxpy(&nxf, &nyf, &nzf,
315  &alpha, RAT(x, VAT2(iz, 1, lev)), w1);
316  rsnrm = Vxnrm1(&nxf, &nyf, &nzf, w1);
317  Vxcopy(&nxf, &nyf, &nzf,
318  RAT(x, VAT2(iz, 1, lev)),
319  RAT(tru, VAT2(iz, 1, lev)));
320  } else if (*istop == 3) {
321  Vxcopy(&nxf, &nyf, &nzf,
322  RAT(tru, VAT2(iz, 1, lev)),
323  w1);
324  alpha = -1.0;
325  Vxaxpy(&nxf, &nyf, &nzf,
326  &alpha, RAT(x, VAT2(iz, 1, lev)), w1);
327  rsnrm = Vxnrm2(&nxf, &nyf, &nzf, w1);
328  } else if (*istop == 4) {
329  Vxcopy(&nxf, &nyf, &nzf,
330  RAT(tru, VAT2(iz, 1, lev)),
331  w1);
332  alpha = -1.0;
333  Vxaxpy(&nxf, &nyf, &nzf,
334  &alpha, RAT(x, VAT2(iz, 1, lev)), w1);
335  rsnrm = Vxnrm2(&nxf, &nyf, &nzf, w1);
336  } else if (*istop == 5) {
337  Vxcopy(&nxf, &nyf, &nzf,
338  RAT(tru, VAT2(iz, 1, lev)),
339  w1);
340  alpha = -1.0;
341  Vxaxpy(&nxf, &nyf, &nzf,
342  &alpha, RAT(x, VAT2(iz, 1, lev)), w1);
343  Vnmatvec(&nxf, &nyf, &nzf,
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)),
348  w1, w2, w3);
349  rsnrm = VSQRT(Vxdot(&nxf, &nyf, &nzf, w1, w2));
350  } else {
351  Vnm_print(2, "Vmvcs: bad istop value: %d\n", *istop);
352  }
353 
354  Vprtstp(*iok, *iters, rsnrm, rsden, orsnrm);
355  }
356  return;
357  }
358 
359  /*********************************************************************
360  *** begin mg iteration (note nxf,nyf,nzf changes during loop)
361  *********************************************************************/
362 
363  // setup for the v-cycle looping ***
364  *iters = 0;
365  while(1) {
366 
367  // Finest level initialization
368  level = 1;
369  lev = (*ilev - 1) + level;
370 
371  // Nu1 pre-smoothings on fine grid (with residual)
372  iresid = 1;
373  iadjoint = 0;
374  iters_s = 0;
375  errtol_s = 0.0;
376  nuuu = Vivariv(nu1, &lev);
377  Vnsmooth(&nxf, &nyf, &nzf,
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)),
381  w2, w3, w1,
382  &nuuu, &iters_s, &errtol_s, omega,
383  &iresid, &iadjoint, mgsmoo);
384  Vxcopy(&nxf, &nyf, &nzf, w1, RAT(w0, VAT2(iz, 1, lev)));
385 
386  /* *********************************************************************
387  * begin cycling down to coarse grid
388  * ********************************************************************/
389 
390  // Go down grids: restrict resid to coarser and smooth
391  for (level = 2; level <= *nlev; level++) {
392  lev = (*ilev - 1) + level;
393 
394  // Find new grid size
395  numlev = 1;
396  Vmkcors(&numlev, &nxf, &nyf, &nzf, &nxc, &nyc, &nzc);
397 
398  // Restrict residual to coarser grid
399  Vrestrc(&nxf, &nyf, &nzf, &nxc, &nyc, &nzc,
400  w1, RAT(w0, VAT2(iz, 1, lev)),
401  RAT(pc, VAT2(iz, 11, lev-1)));
402 
403  // Restrict (extract) solution to coarser grid
404  Vextrac(&nxf, &nyf, &nzf, &nxc, &nyc, &nzc,
405  RAT(x, VAT2(iz, 1, lev-1)), RAT(w4, VAT2(iz, 1, lev)));
406 
407  // New grid size
408  nxf = nxc;
409  nyf = nyc;
410  nzf = nzc;
411 
412  // Apply coarse grid operator to coarse grid soln
413  Vnmatvec(&nxf, &nyf, &nzf,
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)),
417  w3);
418 
419  // Build coarse grid right hand side
420  alpha = 1.0;
421  Vxaxpy(&nxf, &nyf, &nzf, &alpha,
422  RAT(w0, VAT2(iz, 1, lev)), RAT(fc, VAT2(iz, 1, lev)));
423 
424  // If not on coarsest level yet...
425  if (level != *nlev) {
426 
427  // nu1 pre-smoothings on this level (with residual)
428  Vxcopy(&nxf, &nyf, &nzf,
429  RAT(w4, VAT2(iz, 1, lev)), RAT(x, VAT2(iz, 1, lev)));
430  iresid = 1;
431  iadjoint = 0;
432  iters_s = 0;
433  errtol_s = 0.0;
434  nuuu = Vivariv (nu1,&lev);
435  Vnsmooth(&nxf, &nyf, &nzf,
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)),
440  w2, w3, w1,
441  &nuuu, &iters_s, &errtol_s, omega,
442  &iresid, &iadjoint, mgsmoo);
443  }
444 
445  // End of cycling down to coarse grid loop
446  }
447 
448 
449 
450  /* *********************************************************************
451  * begin coarse grid
452  * *********************************************************************/
453 
454  // Coarsest level
455  level = *nlev;
456  lev = (*ilev - 1) + level;
457 
458  // Solve on coarsest grid with ncghs, mgsmoo_s=4 (no residual)
459  iresid = 0;
460  iadjoint = 0;
461  itmax_s = 100;
462  iters_s = 0;
463  errtol_s = *epsiln;
464  mgsmoo_s = *mgsmoo;
465  Vazeros(&nxf, &nyf, &nzf, RAT(x, VAT2(iz, 1, lev)));
466  Vnsmooth(&nxf, &nyf, &nzf,
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)),
471  w1, w2, w3,
472  &itmax_s, &iters_s, &errtol_s, omega,
473  &iresid, &iadjoint, &mgsmoo_s);
474 
475  /* *********************************************************************
476  * begin cycling back to fine grid
477  * ********************************************************************/
478 
479  // Move up grids: interpolate resid to finer and smooth
480  for (level = *nlev - 1; level >= 1; level--) {
481  lev = (*ilev - 1) + level;
482 
483  // Find new grid size
484  numlev = 1;
485  Vmkfine(&numlev, &nxf, &nyf, &nzf, &nxc, &nyc, &nzc);
486 
487  // Form difference of new approx at the coarse grid
488  alpha = -1.0;
489  Vxaxpy(&nxf, &nyf, &nzf, &alpha,
490  RAT(w4, VAT2(iz, 1, lev + 1)), RAT(x, VAT2(iz, 1, lev + 1)));
491 
492  // Call the line search (on the coarser level)
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)),
502  w1, w2, w3);
503 
504  // Interpolate to next finer grid
505  VinterpPMG(&nxf, &nyf, &nzf, &nxc, &nyc, &nzc,
506  RAT(x, VAT2(iz, 1, lev + 1)),
507  w1,
508  RAT(pc, VAT2(iz, 11, lev)));
509 
510  // New grid size
511  nxf = nxc;
512  nyf = nyc;
513  nzf = nzc;
514 
515  // Perform the coarse grid correction
516  Vxaxpy(&nxf, &nyf, &nzf, &xdamp, w1, RAT(x, VAT2(iz, 1, lev)));
517 
518  // nu2 post-smoothings for correction (no residual) ***
519  iresid = 0;
520  iadjoint = 1;
521  iters_s = 0;
522  errtol_s = 0.0;
523  nuuu = Vivariv(nu2, &lev);
524  Vnsmooth(&nxf, &nyf, &nzf,
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)),
529  w1, w2, w3,
530  &nuuu, &iters_s, &errtol_s, omega,
531  &iresid, &iadjoint, mgsmoo);
532  }
533 
534 
535  /* *********************************************************************
536  * iteration complete: do some i/o
537  * ********************************************************************/
538 
539  // Increment the iteration counter
540  iters = iters + 1;
541 
542  // Compute/check the current stopping test
543  if (*iok != 0) {
544  orsnrm = rsnrm;
545  if (*istop == 0) {
546  Vnmresid(&nxf, &nyf, &nzf,
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)),
550  w1, w2);
551  rsnrm = Vxnrm1(&nxf, &nyf, &nzf, w1);
552  } else if (*istop == 1) {
553  Vnmresid(&nxf, &nyf, &nzf,
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)),
557  w1, w2);
558  rsnrm = Vxnrm1(&nxf, &nyf, &nzf,w1);
559  } else if (*istop == 2) {
560  Vxcopy(&nxf, &nyf, &nzf, RAT(tru, VAT2(iz, 1, lev)), w1);
561  alpha = -1.0;
562  Vxaxpy(&nxf, &nyf, &nzf, &alpha, RAT(x, VAT2(iz, 1, lev)), w1);
563  rsnrm = Vxnrm1(&nxf, &nyf, &nzf, w1);
564  Vxcopy(&nxf, &nyf, &nzf,
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);
569  alpha = -1.0;
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);
574  alpha = -1.0;
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);
579  alpha = -1.0;
580  Vxaxpy(&nxf, &nyf, &nzf, &alpha, RAT(x, VAT2(iz, 1, lev)), w1);
581  Vnmatvec(&nxf, &nyf, &nzf,
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)),
584  w1, w2, w3);
585  rsnrm = VSQRT(Vxdot(&nxf, &nyf, &nzf,w1,w2));
586  } else {
587  VABORT_MSG1("Bad istop value: %d", *istop);
588  }
589 
590  Vprtstp(*iok, *iters, rsnrm, rsden, orsnrm);
591 
592  if ((rsnrm / rsden) <= *errtol)
593  break;
594  }
595 
596  if (*iters >= *itmax) {
597  *ierror = 1;
598  break;
599  }
600  }
601 }
int mgsolv
Definition: vpmgp.h:174
int iinfo
Definition: vpmgp.h:130
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.
Definition: matvecd.c:1074
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:911
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.
Definition: mgfasd.c:52
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 Vxcopy(int *nx, int *ny, int *nz, double *x, double *y)
A collection of useful low-level routines (timing, etc).
Definition: mikpckd.c:52
int mgsmoo
Definition: vpmgp.h:160
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
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:228
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.
Definition: mgfasd.c:152
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
int nyc
Definition: vpmgp.h:97
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:594
int nu2
Definition: vpmgp.h:159
int itmax
Definition: vpmgp.h:122
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.
Definition: smoothd.c:92
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:778
int nz
Definition: vpmgp.h:85