APBS  1.4.1
mgcsd.c
1 
50 #include "mgcsd.h"
51 
52 VEXTERNC void Vmvcs(int *nx, int *ny, int *nz,
53  double *x,
54  int *iz,
55  double *w0, double *w1, double *w2, double *w3,
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,
61  int *mgsmoo,
62  int *ipc, double *rpc,
63  double *pc, double *ac, double *cc, double *fc, double *tru) {
64 
65  int level; // @todo: doc
66  int lev; // @todo: doc
67  int itmax_s; // @todo: doc
68  int iters_s; // @todo: doc
69  int nuuu; // @todo: doc
70  int mgsmoo_s; // @todo: doc
71  int iresid; // @todo: doc
72  int nxf; // @todo: doc
73  int nyf; // @todo: doc
74  int nzf; // @todo: doc
75  int nxc; // @todo: doc
76  int nyc; // @todo: doc
77  int nzc; // @todo: doc
78  int lpv; // @todo: doc
79  int n; // @todo: doc
80  int m; // @todo: doc
81  int iadjoint; // @todo: doc
82  double errtol_s; // @todo: doc
83  double rsden; // @todo: doc
84  double rsnrm; // @todo: doc
85  double orsnrm; // @todo: doc
86  double xnum; // @todo: doc
87  double xden; // @todo: doc
88  double xdamp; // @todo: doc
89  int lda; // @todo: doc
90 
91  double alpha; // A utility variable used to pass a parameter to xaxpy
92  int numlev; // A utility variable used to pass a parameter to mkcors
93 
94  MAT2(iz, 50, 1);
95 
96  // Recover level information
97  level = 1;
98  lev = (*ilev - 1) + level;
99 
100  // Recover grid sizes
101  nxf = *nx;
102  nyf = *ny;
103  nzf = *nz;
104  numlev = *nlev - 1;
105  Vmkcors(&numlev, &nxf, &nyf, &nzf, &nxc, &nyc, &nzc);
106 
107  // Do some i/o if requested
108  if (*iinfo > 1) {
109  VMESSAGE0("Starting mvcs operation");
110  VMESSAGE3("Fine Grid Size: (%d, %d, %d)", nxf, nyf, nzf);
111  VMESSAGE3("Coarse Grid Size: (%d, %d, %d)", nxc, nyc, nzc);
112  }
113 
114  if (*iok != 0) {
115  Vprtstp(*iok, -1, 0.0, 0.0, 0.0);
116 
117  }
118 
119  /* **************************************************************
120  * *** Note: if (iok != 0) then: use a stopping test. ***
121  * *** else: use just the itmax to stop iteration. ***
122  * **************************************************************
123  * *** istop=0 most efficient (whatever it is) ***
124  * *** istop=1 relative residual ***
125  * *** istop=2 rms difference of successive iterates ***
126  * *** istop=3 relative true error (provided for testing) ***
127  * **************************************************************/
128 
129  // Compute denominator for stopping criterion
130  if (*iok != 0) {
131  if (*istop == 0) {
132  rsden = 1.0;
133  }
134  else if (*istop == 1) {
135  rsden = Vxnrm1(&nxf, &nyf, &nzf, RAT(fc, VAT2(iz, 1,lev)));
136  }
137  else if (*istop == 2) {
138  rsden = VSQRT(nxf * nyf * nzf);
139  }
140  else if (*istop == 3) {
141  rsden = Vxnrm2(&nxf, &nyf, &nzf, RAT(tru, VAT2(iz, 1,lev)));
142  }
143  else if (*istop == 4) {
144  rsden = Vxnrm2(&nxf, &nyf, &nzf, RAT(tru, VAT2(iz, 1,lev)));
145  }
146  else if (*istop == 5) {
147  Vmatvec(&nxf, &nyf, &nzf,
148  RAT(ipc, VAT2(iz, 5,lev)), RAT(rpc, VAT2(iz, 6,lev)),
149  RAT(ac, VAT2(iz, 7,lev)), RAT(cc, VAT2(iz, 1,lev)),
150  RAT(tru, VAT2(iz, 1,lev)), w1);
151  rsden = VSQRT(Vxdot(&nxf, &nyf, &nzf, RAT(tru, VAT2(iz, 1,lev)), w1));
152  }
153  else {
154  VABORT_MSG1("Bad istop value: %d", *istop);
155  }
156 
157  if (rsden == 0.0) {
158  rsden = 1.0;
159  VERRMSG0("rhs is zero on finest level");
160  }
161  rsnrm = rsden;
162  orsnrm = rsnrm;
163  iters_s = 0;
164 
165  Vprtstp(*iok, 0, rsnrm, rsden, orsnrm);
166  }
167 
168 
169 
170  /* *********************************************************************
171  * *** solve directly if nlev = 1
172  * *********************************************************************/
173 
174  // Solve directly if on the coarse grid
175  if (*nlev == 1) {
176 
177  // Use iterative method?
178  if (*mgsolv == 0) {
179 
180  // solve on coarsest grid with cghs, mgsmoo_s=4 (no residual)
181  iresid = 0;
182  iadjoint = 0;
183  itmax_s = 100;
184  iters_s = 0;
185  errtol_s = *epsiln;
186  mgsmoo_s = 4;
187 
188  Vazeros(&nxf, &nyf, &nzf, RAT(x, VAT2(iz, 1,lev)));
189 
190  Vsmooth(&nxf, &nyf, &nzf,
191  RAT(ipc, VAT2(iz, 5,lev)), RAT(rpc, VAT2(iz, 6,lev)),
192  RAT(ac, VAT2(iz, 7,lev)), RAT(cc, VAT2(iz, 1,lev)), RAT(fc, VAT2(iz, 1,lev)),
193  RAT(x, VAT2(iz, 1,lev)), w1, w2, w3,
194  &itmax_s, &iters_s, &errtol_s, omega,
195  &iresid, &iadjoint, &mgsmoo_s);
196 
197  // Check for trouble on the coarse grid
198  VWARN_MSG2(iters_s <= itmax_s,
199  "Exceeded maximum iterations: iters_s=%d, itmax_s=%d",
200  iters_s, itmax_s);
201 
202  } else if (*mgsolv == 1) {
203 
204  // Use direct method?
205 
206  // Setup lpv to access the factored/banded operator
207  lpv = lev + 1;
208 
209  // setup for banded format
210  n = *RAT(ipc, (VAT2(iz, 5,lpv) - 1) + 1);
211  m = *RAT(ipc, (VAT2(iz, 5,lpv) - 1) + 2);
212  lda = *RAT(ipc, (VAT2(iz, 5,lpv) - 1) + 3);
213 
214  // Call dpbsl to solve
215  Vxcopy_small(&nxf, &nyf, &nzf, RAT(fc, VAT2(iz, 1,lev)), w1);
216  Vdpbsl(RAT(ac, VAT2(iz, 7,lpv)), &lda, &n, &m, w1);
217  Vxcopy_large(&nxf, &nyf, &nzf, w1, RAT(x, VAT2(iz, 1,lev)));
218  VfboundPMG00(&nxf, &nyf, &nzf, RAT(x, VAT2(iz, 1,lev)));
219 
220  } else {
221  VABORT_MSG1("Invalid coarse solver requested: %d", *mgsolv);
222  }
223 
224 
225  // Compute the stopping test
226  *iters = 1;
227  if (*iok != 0) {
228 
229  orsnrm = rsnrm;
230 
231  if (*istop == 0) {
232 
233  Vmresid(&nxf, &nyf, &nzf,
234  RAT(ipc, VAT2(iz, 5, lev)), RAT(rpc, VAT2(iz, 6, lev)),
235  RAT( ac, VAT2(iz, 7, lev)), RAT(cc , VAT2(iz, 1, lev)),
236  RAT( fc, VAT2(iz, 1,lev)),
237  RAT( x, VAT2(iz, 1, lev)), w1);
238 
239  rsnrm = Vxnrm1(&nxf, &nyf, &nzf, w1);
240  }
241 
242  else if (*istop == 1) {
243 
244  Vmresid(&nxf, &nyf, &nzf,
245  RAT(ipc, VAT2(iz, 5, lev)), RAT(rpc, VAT2(iz, 6, lev)),
246  RAT( ac, VAT2(iz, 7, lev)), RAT( cc, VAT2(iz, 1, lev)),
247  RAT( fc, VAT2(iz, 1, lev)), RAT( x, VAT2(iz, 1, lev)),
248  w1);
249  rsnrm = Vxnrm1(&nxf, &nyf, &nzf, w1);
250  }
251 
252  else if (*istop == 2) {
253 
254  alpha = -1.0;
255 
256  Vxcopy(&nxf, &nyf, &nzf, RAT(tru, VAT2(iz, 1,lev)), w1);
257  Vxaxpy(&nxf, &nyf, &nzf, &alpha,
258  RAT(x, VAT2(iz, 1,lev)), w1);
259  rsnrm = Vxnrm1(&nxf, &nyf, &nzf, w1);
260  Vxcopy(&nxf, &nyf, &nzf, RAT(x, VAT2(iz, 1,lev)), RAT(tru, VAT2(iz, 1,lev)));
261  }
262 
263  else if (*istop == 3) {
264 
265  alpha = -1.0;
266 
267  Vxcopy(&nxf, &nyf, &nzf, RAT(tru, VAT2(iz, 1,lev)), w1);
268  Vxaxpy(&nxf, &nyf, &nzf, &alpha, RAT(x, VAT2(iz, 1,lev)), w1);
269  rsnrm = Vxnrm2(&nxf, &nyf, &nzf, w1);
270  }
271 
272  else if (*istop == 4) {
273 
274  alpha = -1.0;
275 
276  Vxcopy(&nxf, &nyf, &nzf, RAT(tru, VAT2(iz, 1,lev)), w1);
277  Vxaxpy(&nxf, &nyf, &nzf, &alpha, RAT(x, VAT2(iz, 1,lev)), w1);
278  rsnrm = Vxnrm2(&nxf, &nyf, &nzf, w1);
279  }
280 
281  else if (*istop == 5) {
282 
283  alpha = -1.0;
284 
285  Vxcopy(&nxf, &nyf, &nzf, RAT(tru, VAT2(iz, 1,lev)), w1);
286  Vxaxpy(&nxf, &nyf, &nzf, &alpha, RAT(x, VAT2(iz, 1,lev)), w1);
287 
288  Vmatvec(&nxf, &nyf, &nzf,
289  RAT(ipc, VAT2(iz, 5,lev)), RAT(rpc, VAT2(iz, 6,lev)),
290  RAT(ac, VAT2(iz, 7,lev)), RAT(cc, VAT2(iz, 1,lev)),
291  w1, w2);
292  rsnrm = VSQRT(Vxdot(&nxf, &nyf, &nzf, w1, w2));
293  }
294 
295  else {
296  VABORT_MSG1("Bad istop value: %d\n", *istop);
297  }
298  Vprtstp(*iok, *iters, rsnrm, rsden, orsnrm);
299  }
300  return;
301  }
302 
303 
304  /* *********************************************************************
305  * *** begin mg iteration (note nxf,nyf,nzf changes during loop)
306  * *********************************************************************/
307 
308  // Setup for the v-cycle looping
309  *iters = 0;
310  do {
311 
312  // Finest level initialization
313  level = 1;
314  lev = (*ilev - 1) + level;
315 
316  // nu1 pre-smoothings on fine grid (with residual)
317  iresid = 1;
318  iadjoint = 0;
319  iters_s = 0;
320  errtol_s = 0.0;
321  nuuu = Vivariv(nu1, &lev);
322 
323  Vsmooth(&nxf, &nyf, &nzf,
324  RAT(ipc, VAT2(iz, 5,lev)), RAT(rpc, VAT2(iz, 6,lev)),
325  RAT(ac, VAT2(iz, 7,lev)), RAT(cc, VAT2(iz, 1,lev)), RAT(fc, VAT2(iz, 1,lev)),
326  RAT(x, VAT2(iz, 1,lev)), w2, w3, w1,
327  &nuuu, &iters_s,
328  &errtol_s, omega,
329  &iresid, &iadjoint, mgsmoo);
330 
331  Vxcopy(&nxf, &nyf, &nzf, w1, RAT(w0, VAT2(iz, 1,lev)));
332 
333 
334 
335  /* *********************************************************************
336  * begin cycling down to coarse grid
337  * *********************************************************************/
338 
339  // Go down grids: restrict resid to coarser and smooth
340  for (level=2; level<=*nlev; level++) {
341 
342  lev = (*ilev - 1) + level;
343 
344  // Find new grid size
345  numlev = 1;
346  Vmkcors(&numlev, &nxf, &nyf, &nzf, &nxc, &nyc, &nzc);
347 
348  // Restrict residual to coarser grid ***
349  Vrestrc(&nxf, &nyf, &nzf,
350  &nxc, &nyc, &nzc,
351  w1, RAT(w0, VAT2(iz, 1,lev)), RAT(pc, VAT2(iz, 11,lev-1)));
352 
354  nxf = nxc;
355  nyf = nyc;
356  nzf = nzc;
357 
358  // if not on coarsest level yet...
359  if (level != *nlev) {
360 
361  // nu1 pre-smoothings on this level (with residual)
362  // (w1 has residual...)
363  Vazeros(&nxf, &nyf, &nzf, RAT(x, VAT2(iz, 1,lev)));
364  iresid = 1;
365  iadjoint = 0;
366  iters_s = 0;
367  errtol_s = 0.0;
368  nuuu = Vivariv(nu1, &lev);
369  Vsmooth(&nxf, &nyf, &nzf,
370  RAT(ipc, VAT2(iz, 5,lev)), RAT(rpc, VAT2(iz, 6,lev)),
371  RAT(ac, VAT2(iz, 7,lev)), RAT(cc, VAT2(iz, 1,lev)), RAT(w0, VAT2(iz, 1,lev)),
372  RAT(x, VAT2(iz, 1,lev)), w2, w3, w1,
373  &nuuu, &iters_s,
374  &errtol_s, omega ,
375  &iresid, &iadjoint, mgsmoo);
376  }
377  // End of cycling down to coarse grid loop
378  }
379 
380 
381 
382  /* *********************************************************************
383  * begin coarse grid
384  * *********************************************************************/
385 
386  // Coarsest level
387  level = *nlev;
388  lev = (*ilev - 1) + level;
389 
390  // Use iterative method?
391  if (*mgsolv == 0) {
392 
393  // solve on coarsest grid with cghs, mgsmoo_s=4 (no residual)
394  iresid = 0;
395  iadjoint = 0;
396  itmax_s = 100;
397  iters_s = 0;
398  errtol_s = *epsiln;
399  mgsmoo_s = 4;
400  Vazeros(&nxf, &nyf, &nzf, RAT(x, VAT2(iz, 1,lev)));
401  Vsmooth(&nxf, &nyf, &nzf,
402  RAT(ipc, VAT2(iz, 5,lev)), RAT(rpc, VAT2(iz, 6,lev)),
403  RAT(ac, VAT2(iz, 7,lev)), RAT(cc, VAT2(iz, 1,lev)), RAT(w0, VAT2(iz, 1,lev)),
404  RAT(x, VAT2(iz, 1,lev)), w1, w2, w3,
405  &itmax_s, &iters_s,
406  &errtol_s, omega,
407  &iresid, &iadjoint, &mgsmoo_s);
408 
409  // Check for trouble on the coarse grid
410  VWARN_MSG2(iters_s <= itmax_s,
411  "Exceeded maximum iterations: iters_s=%d, itmax_s=%d",
412  iters_s, itmax_s);
413  } else if (*mgsolv == 1) {
414 
415  // use direct method?
416 
417  // Setup lpv to access the factored/banded operator
418  lpv = lev + 1;
419 
420  // Setup for banded format
421  n = VAT(ipc, (VAT2(iz, 5, lpv) - 1) + 1);
422  m = VAT(ipc, (VAT2(iz, 5, lpv) - 1) + 2);
423  lda = VAT(ipc, (VAT2(iz, 5, lpv) - 1) + 3);
424 
425  // Call dpbsl to solve
426  Vxcopy_small(&nxf, &nyf, &nzf, RAT(w0, VAT2(iz, 1,lev)), w1);
427  Vdpbsl(RAT(ac, VAT2(iz, 7,lpv)), &lda, &n, &m, w1);
428  Vxcopy_large(&nxf, &nyf, &nzf, w1, RAT(x, VAT2(iz, 1,lev)));
429  VfboundPMG00(&nxf, &nyf, &nzf, RAT(x, VAT2(iz, 1,lev)));
430 
431  } else {
432  VABORT_MSG1("Invalid coarse solver requested: %d", *mgsolv);
433  }
434 
435 
436  /* *********************************************************************
437  * begin cycling back to fine grid
438  * *********************************************************************/
439 
440 
441  // Move up grids: interpolate resid to finer and smooth
442  for (level=*nlev-1; level>=1; level--) {
443 
444  lev = (*ilev - 1) + level;
445 
446  // Find new grid size
447  numlev = 1;
448  Vmkfine(&numlev,
449  &nxf, &nyf, &nzf,
450  &nxc, &nyc, &nzc);
451 
452  // Interpolate to next finer grid
453  VinterpPMG(&nxf, &nyf, &nzf,
454  &nxc, &nyc, &nzc,
455  RAT(x, VAT2(iz, 1,lev+1)), w1, RAT(pc, VAT2(iz, 11,lev)));
456 
457  /* Compute the hackbusch/reusken damping parameter
458  * which is equivalent to the standard linear cg steplength
459  */
460  Vmatvec(&nxf, &nyf, &nzf,
461  RAT(ipc, VAT2(iz, 5,lev+1)), RAT(rpc, VAT2(iz, 6,lev+1)),
462  RAT(ac, VAT2(iz, 7,lev+1)), RAT(cc, VAT2(iz, 1,lev+1)),
463  RAT(x, VAT2(iz, 1,lev+1)), w2);
464 
465  xnum = Vxdot(&nxf, &nyf, &nzf,
466  RAT(x, VAT2(iz, 1,lev+1)), RAT(w0, VAT2(iz, 1,lev+1)));
467 
468  xden = Vxdot(&nxf, &nyf, &nzf,
469  RAT(x, VAT2(iz, 1,lev+1)), w2);
470  xdamp = xnum / xden;
471 
472  // New grid size
473  nxf = nxc;
474  nyf = nyc;
475  nzf = nzc;
476 
477  // perform the coarse grid correction
478  // xdamp = 1.0d0
479  Vxaxpy(&nxf, &nyf, &nzf,
480  &xdamp, w1, RAT(x, VAT2(iz, 1,lev)));
481 
482  // nu2 post-smoothings for correction (no residual)
483  iresid = 0;
484  iadjoint = 1;
485  iters_s = 0;
486  errtol_s = 0.0;
487  nuuu = Vivariv(nu2, &lev);
488  if (level == 1) {
489  Vsmooth(&nxf, &nyf, &nzf,
490  RAT(ipc, VAT2(iz, 5,lev)), RAT(rpc, VAT2(iz, 6,lev)),
491  RAT(ac, VAT2(iz, 7,lev)), RAT(cc, VAT2(iz, 1,lev)), RAT(fc, VAT2(iz, 1,lev)),
492  RAT(x, VAT2(iz, 1,lev)),w1,w2,w3,
493  &nuuu, &iters_s, &errtol_s, omega,
494  &iresid, &iadjoint, mgsmoo);
495  } else {
496  Vsmooth(&nxf, &nyf, &nzf,
497  RAT(ipc, VAT2(iz, 5,lev)), RAT(rpc, VAT2(iz, 6,lev)),
498  RAT(ac, VAT2(iz, 7,lev)), RAT(cc, VAT2(iz, 1,lev)), RAT(w0, VAT2(iz, 1,lev)),
499  RAT(x, VAT2(iz, 1,lev)), w1, w2, w3,
500  &nuuu, &iters_s, &errtol_s, omega,
501  &iresid, &iadjoint, mgsmoo);
502  }
503  }
504 
505  /* *********************************************************************
506  * iteration complete: do some i/o
507  * *********************************************************************/
508 
509  // Increment the iteration counter
510  (*iters)++;
511 
512  // Compute/check the current stopping test
513  if (iok != 0) {
514  orsnrm = rsnrm;
515  if (*istop == 0) {
516  Vmresid(&nxf, &nyf, &nzf,
517  RAT(ipc, VAT2(iz, 5,lev)), RAT(rpc, VAT2(iz, 6,lev)),
518  RAT(ac, VAT2(iz, 7,lev)), RAT(cc, VAT2(iz, 1,lev)), RAT(fc, VAT2(iz, 1,lev)),
519  RAT(x, VAT2(iz, 1,lev)), w1);
520  rsnrm = Vxnrm1(&nxf, &nyf, &nzf, w1);
521  } else if(*istop == 1) {
522  Vmresid(&nxf, &nyf, &nzf,
523  RAT(ipc, VAT2(iz, 5,lev)), RAT(rpc, VAT2(iz, 6,lev)),
524  RAT(ac, VAT2(iz, 7,lev)), RAT(cc, VAT2(iz, 1,lev)), RAT(fc, VAT2(iz, 1,lev)),
525  RAT(x, VAT2(iz, 1,lev)), w1);
526  rsnrm = Vxnrm1(&nxf, &nyf, &nzf, w1);
527  } else if (*istop == 2) {
528  Vxcopy(&nxf, &nyf, &nzf, RAT(tru, VAT2(iz, 1,lev)), w1);
529  alpha = -1.0;
530  Vxaxpy(&nxf, &nyf, &nzf, &alpha, RAT(x, VAT2(iz, 1,lev)), w1);
531  rsnrm = Vxnrm1(&nxf, &nyf, &nzf, w1);
532  Vxcopy(&nxf, &nyf, &nzf, RAT(x, VAT2(iz, 1,lev)), RAT(tru, VAT2(iz, 1,lev)));
533  } else if (*istop == 3) {
534  Vxcopy(&nxf, &nyf, &nzf, RAT(tru, VAT2(iz, 1,lev)), w1);
535  alpha = -1.0;
536  Vxaxpy(&nxf, &nyf, &nzf, &alpha, RAT(x, VAT2(iz, 1,lev)), w1);
537  rsnrm = Vxnrm2(&nxf, &nyf, &nzf, w1);
538  } else if (*istop == 4) {
539  Vxcopy(&nxf, &nyf, &nzf, RAT(tru, VAT2(iz, 1,lev)), w1);
540  alpha = -1.0;
541  Vxaxpy(&nxf, &nyf, &nzf, &alpha, RAT(x, VAT2(iz, 1,lev)), w1);
542  rsnrm = Vxnrm2(&nxf, &nyf, &nzf, w1);
543  } else if (*istop == 5) {
544  Vxcopy(&nxf, &nyf, &nzf, RAT(tru, VAT2(iz, 1,lev)), w1);
545  alpha = -1.0;
546  Vxaxpy(&nxf, &nyf, &nzf, &alpha, RAT(x, VAT2(iz, 1,lev)), w1);
547  Vmatvec(&nxf, &nyf, &nzf,
548  RAT(ipc, VAT2(iz, 5,lev)), RAT(rpc, VAT2(iz, 6,lev)),
549  RAT(ac, VAT2(iz, 7,lev)), RAT(cc, VAT2(iz, 1,lev)),
550  w1, w2);
551  rsnrm = VSQRT(Vxdot(&nxf, &nyf, &nzf, w1, w2));
552  } else {
553  VABORT_MSG1("Bad istop value: %d", *istop);
554  }
555  Vprtstp(*iok, *iters, rsnrm, rsden, orsnrm);
556  }
557  } while (*iters<*itmax && (rsnrm/rsden) > *errtol);
558 
559  *ierror = *iters < *itmax ? 0 : 1;
560 }
VPUBLIC void VfboundPMG00(int *nx, int *ny, int *nz, double *x)
Initialize a grid function to have a zero boundary value.
Definition: mikpckd.c:253
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:911
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
VEXTERNC void Vsmooth(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)
Multigrid smoothing functions.
Definition: smoothd.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
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
int nu1
Definition: vpmgp.h:158
VPUBLIC void Vmatvec(int *nx, int *ny, int *nz, int *ipc, double *rpc, double *ac, double *cc, double *x, double *y)
Matrix-vector multiplication routines.
Definition: matvecd.c:52
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 Vdpbsl(double *abd, int *lda, int *n, int *m, double *b)
LINPACK interface.
Definition: mlinpckd.c:52
int nyc
Definition: vpmgp.h:97
VPUBLIC void Vxcopy_small(int *nx, int *ny, int *nz, double *x, double *y)
Copy operation for a grid function with boundary values. Quite simply copies one 3d matrix to another...
Definition: mikpckd.c:70
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
int nu2
Definition: vpmgp.h:159
int itmax
Definition: vpmgp.h:122
VPUBLIC void Vmresid(int *nx, int *ny, int *nz, int *ipc, double *rpc, double *ac, double *cc, double *fc, double *x, double *r)
Break the matrix data-structure into diagonals and then call the residual routine.
Definition: matvecd.c:422
VPUBLIC void Vxcopy_large(int *nx, int *ny, int *nz, double *x, double *y)
Copy operation for a grid function with boundary values. Quite simply copies one 3d matrix to another...
Definition: mikpckd.c:86
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