APBS 3.0.0
Loading...
Searching...
No Matches
mgcsd.c
1
55#include "mgcsd.h"
56
57VEXTERNC 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