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,
59 double *epsiln,
double *
errtol,
double *omega,
62 int *ipc,
double *rpc,
63 double *pc,
double *ac,
double *cc,
double *fc,
double *tru) {
98 lev = (*ilev - 1) + level;
105 Vmkcors(&numlev, &nxf, &nyf, &nzf, &
nxc, &
nyc, &
nzc);
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);
115 Vprtstp(*iok, -1, 0.0, 0.0, 0.0);
134 else if (*
istop == 1) {
135 rsden =
Vxnrm1(&nxf, &nyf, &nzf, RAT(fc, VAT2(iz, 1,lev)));
137 else if (*
istop == 2) {
138 rsden = VSQRT(nxf * nyf * nzf);
140 else if (*
istop == 3) {
141 rsden =
Vxnrm2(&nxf, &nyf, &nzf, RAT(tru, VAT2(iz, 1,lev)));
143 else if (*
istop == 4) {
144 rsden =
Vxnrm2(&nxf, &nyf, &nzf, RAT(tru, VAT2(iz, 1,lev)));
146 else if (*
istop == 5) {
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));
154 VABORT_MSG1(
"Bad istop value: %d", *
istop);
159 VERRMSG0(
"rhs is zero on finest level");
165 Vprtstp(*iok, 0, rsnrm, rsden, orsnrm);
188 Vazeros(&nxf, &nyf, &nzf, RAT(x, VAT2(iz, 1,lev)));
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);
198 VWARN_MSG2(iters_s <= itmax_s,
199 "Exceeded maximum iterations: iters_s=%d, itmax_s=%d",
202 }
else if (*
mgsolv == 1) {
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);
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)));
221 VABORT_MSG1(
"Invalid coarse solver requested: %d", *
mgsolv);
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);
239 rsnrm =
Vxnrm1(&nxf, &nyf, &nzf, w1);
242 else if (*
istop == 1) {
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)),
249 rsnrm =
Vxnrm1(&nxf, &nyf, &nzf, w1);
252 else if (*
istop == 2) {
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)));
263 else if (*
istop == 3) {
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);
272 else if (*
istop == 4) {
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);
281 else if (*
istop == 5) {
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);
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)),
292 rsnrm = VSQRT(
Vxdot(&nxf, &nyf, &nzf, w1, w2));
296 VABORT_MSG1(
"Bad istop value: %d\n", *
istop);
298 Vprtstp(*iok, *iters, rsnrm, rsden, orsnrm);
314 lev = (*ilev - 1) + level;
321 nuuu = Vivariv(
nu1, &lev);
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,
329 &iresid, &iadjoint,
mgsmoo);
331 Vxcopy(&nxf, &nyf, &nzf, w1, RAT(w0, VAT2(iz, 1,lev)));
340 for (level=2; level<=*
nlev; level++) {
342 lev = (*ilev - 1) + level;
346 Vmkcors(&numlev, &nxf, &nyf, &nzf, &
nxc, &
nyc, &
nzc);
351 w1, RAT(w0, VAT2(iz, 1,lev)), RAT(pc, VAT2(iz, 11,lev-1)));
359 if (level != *
nlev) {
363 Vazeros(&nxf, &nyf, &nzf, RAT(x, VAT2(iz, 1,lev)));
368 nuuu = Vivariv(
nu1, &lev);
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,
375 &iresid, &iadjoint,
mgsmoo);
388 lev = (*ilev - 1) + level;
400 Vazeros(&nxf, &nyf, &nzf, RAT(x, VAT2(iz, 1,lev)));
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,
407 &iresid, &iadjoint, &mgsmoo_s);
410 VWARN_MSG2(iters_s <= itmax_s,
411 "Exceeded maximum iterations: iters_s=%d, itmax_s=%d",
413 }
else if (*
mgsolv == 1) {
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);
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)));
432 VABORT_MSG1(
"Invalid coarse solver requested: %d", *
mgsolv);
442 for (level=*
nlev-1; level>=1; level--) {
444 lev = (*ilev - 1) + level;
455 RAT(x, VAT2(iz, 1,lev+1)), w1, RAT(pc, VAT2(iz, 11,lev)));
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);
465 xnum =
Vxdot(&nxf, &nyf, &nzf,
466 RAT(x, VAT2(iz, 1,lev+1)), RAT(w0, VAT2(iz, 1,lev+1)));
468 xden =
Vxdot(&nxf, &nyf, &nzf,
469 RAT(x, VAT2(iz, 1,lev+1)), w2);
480 &xdamp, w1, RAT(x, VAT2(iz, 1,lev)));
487 nuuu = Vivariv(
nu2, &lev);
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);
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);
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) {
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);
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);
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);
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);
546 Vxaxpy(&nxf, &nyf, &nzf, &alpha, RAT(x, VAT2(iz, 1,lev)), w1);
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)),
551 rsnrm = VSQRT(
Vxdot(&nxf, &nyf, &nzf, w1, w2));
553 VABORT_MSG1(
"Bad istop value: %d", *
istop);
555 Vprtstp(*iok, *iters, rsnrm, rsden, orsnrm);
557 }
while (*iters<*
itmax && (rsnrm/rsden) > *
errtol);
559 *ierror = *iters < *
itmax ? 0 : 1;
VPUBLIC void VfboundPMG00(int *nx, int *ny, int *nz, double *x)
Initialize a grid function to have a zero boundary value.
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.
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.
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.
VPUBLIC void Vazeros(int *nx, int *ny, int *nz, double *x)
Zero out operation for a grid function, including boundary values.
VPUBLIC void Vxcopy(int *nx, int *ny, int *nz, double *x, double *y)
A collection of useful low-level routines (timing, etc).
VPUBLIC double Vxdot(int *nx, int *ny, int *nz, double *x, double *y)
Inner product operation for a grid function with boundary values.
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.
VPUBLIC void Vxaxpy(int *nx, int *ny, int *nz, double *alpha, double *x, double *y)
saxpy operation for a grid function with boundary values.
VPUBLIC void Vdpbsl(double *abd, int *lda, int *n, int *m, double *b)
LINPACK interface.
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...
VPUBLIC double Vxnrm2(int *nx, int *ny, int *nz, double *x)
Norm operation for a grid function with boundary values.
VPUBLIC double Vxnrm1(int *nx, int *ny, int *nz, double *x)
Norm operation for a grid function with boundary values.
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.
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...
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.