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