57VEXTERNC
void Vmvcs(
int *nx,
int *ny,
int *nz,
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,
67 int *ipc,
double *rpc,
68 double *pc,
double *ac,
double *cc,
double *fc,
double *tru) {
103 lev = (*ilev - 1) + level;
110 Vmkcors(&numlev, &nxf, &nyf, &nzf, &nxc, &nyc, &nzc);
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);
120 Vprtstp(*iok, -1, 0.0, 0.0, 0.0);
139 else if (*istop == 1) {
140 rsden =
Vxnrm1(&nxf, &nyf, &nzf, RAT(fc, VAT2(iz, 1,lev)));
142 else if (*istop == 2) {
143 rsden = VSQRT(nxf * nyf * nzf);
145 else if (*istop == 3) {
146 rsden =
Vxnrm2(&nxf, &nyf, &nzf, RAT(tru, VAT2(iz, 1,lev)));
148 else if (*istop == 4) {
149 rsden =
Vxnrm2(&nxf, &nyf, &nzf, RAT(tru, VAT2(iz, 1,lev)));
151 else if (*istop == 5) {
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));
159 VABORT_MSG1(
"Bad istop value: %d", *istop);
164 VERRMSG0(
"rhs is zero on finest level");
170 Vprtstp(*iok, 0, rsnrm, rsden, orsnrm);
193 Vazeros(&nxf, &nyf, &nzf, RAT(x, VAT2(iz, 1,lev)));
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);
203 VWARN_MSG2(iters_s <= itmax_s,
204 "Exceeded maximum iterations: iters_s=%d, itmax_s=%d",
207 }
else if (*mgsolv == 1) {
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);
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)));
226 VABORT_MSG1(
"Invalid coarse solver requested: %d", *mgsolv);
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);
244 rsnrm =
Vxnrm1(&nxf, &nyf, &nzf, w1);
247 else if (*istop == 1) {
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)),
254 rsnrm =
Vxnrm1(&nxf, &nyf, &nzf, w1);
257 else if (*istop == 2) {
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)));
268 else if (*istop == 3) {
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);
277 else if (*istop == 4) {
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);
286 else if (*istop == 5) {
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);
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)),
297 rsnrm = VSQRT(
Vxdot(&nxf, &nyf, &nzf, w1, w2));
301 VABORT_MSG1(
"Bad istop value: %d\n", *istop);
303 Vprtstp(*iok, *iters, rsnrm, rsden, orsnrm);
319 lev = (*ilev - 1) + level;
326 nuuu = Vivariv(nu1, &lev);
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,
334 &iresid, &iadjoint, mgsmoo);
336 Vxcopy(&nxf, &nyf, &nzf, w1, RAT(w0, VAT2(iz, 1,lev)));
345 for (level=2; level<=*nlev; level++) {
347 lev = (*ilev - 1) + level;
351 Vmkcors(&numlev, &nxf, &nyf, &nzf, &nxc, &nyc, &nzc);
356 w1, RAT(w0, VAT2(iz, 1,lev)), RAT(pc, VAT2(iz, 11,lev-1)));
364 if (level != *nlev) {
368 Vazeros(&nxf, &nyf, &nzf, RAT(x, VAT2(iz, 1,lev)));
373 nuuu = Vivariv(nu1, &lev);
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,
380 &iresid, &iadjoint, mgsmoo);
393 lev = (*ilev - 1) + level;
405 Vazeros(&nxf, &nyf, &nzf, RAT(x, VAT2(iz, 1,lev)));
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,
412 &iresid, &iadjoint, &mgsmoo_s);
415 VWARN_MSG2(iters_s <= itmax_s,
416 "Exceeded maximum iterations: iters_s=%d, itmax_s=%d",
418 }
else if (*mgsolv == 1) {
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);
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)));
437 VABORT_MSG1(
"Invalid coarse solver requested: %d", *mgsolv);
447 for (level=*nlev-1; level>=1; level--) {
449 lev = (*ilev - 1) + level;
460 RAT(x, VAT2(iz, 1,lev+1)), w1, RAT(pc, VAT2(iz, 11,lev)));
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);
470 xnum =
Vxdot(&nxf, &nyf, &nzf,
471 RAT(x, VAT2(iz, 1,lev+1)), RAT(w0, VAT2(iz, 1,lev+1)));
473 xden =
Vxdot(&nxf, &nyf, &nzf,
474 RAT(x, VAT2(iz, 1,lev+1)), w2);
485 &xdamp, w1, RAT(x, VAT2(iz, 1,lev)));
492 nuuu = Vivariv(nu2, &lev);
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);
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);
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) {
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);
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);
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);
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);
551 Vxaxpy(&nxf, &nyf, &nzf, &alpha, RAT(x, VAT2(iz, 1,lev)), w1);
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)),
556 rsnrm = VSQRT(
Vxdot(&nxf, &nyf, &nzf, w1, w2));
558 VABORT_MSG1(
"Bad istop value: %d", *istop);
560 Vprtstp(*iok, *iters, rsnrm, rsden, orsnrm);
562 }
while (*iters<*itmax && (rsnrm/rsden) > *errtol);
564 *ierror = *iters < *itmax ? 0 : 1;