57VPUBLIC
void Vmgdriv(
int* iparm,
double* rparm,
58 int* iwork,
double* rwork,
double* u,
59 double* xf,
double* yf,
double* zf,
60 double* gxcf,
double* gycf,
double* gzcf,
61 double* a1cf,
double* a2cf,
double* a3cf,
62 double* ccf,
double* fcf,
double* tcf) {
108 nrwk = VAT(iparm, 1);
109 niwk = VAT(iparm, 2);
113 nlev = VAT(iparm, 6);
116 VASSERT_MSG1(nlev > 0,
"nlev must be positive: %d", nlev);
117 VASSERT_MSG1( nx > 0,
"nx must be positive: %d", nx);
118 VASSERT_MSG1( ny > 0,
"nv must be positive: %d", ny);
119 VASSERT_MSG1( nz > 0,
"nz must be positive: %d", nz);
121 mxlv = Vmaxlev(nx, ny, nz);
124 "number of levels exceeds maximum: %d > %d",
130 mgcoar = VAT(iparm, 18);
131 mgdisc = VAT(iparm, 19);
132 mgsolv = VAT(iparm, 21);
134 Vmgsz(&mgcoar, &mgdisc, &mgsolv,
140 &n_rpc, &n_iz, &n_ipc,
146 "real workspace exceeds maximum size: %d > %d",
151 "integer workspace exceeds maximum size: %d > %d",
161 k_cc = k_rpc + n_rpc;
164 k_ac = k_pc + 27 * narrc;
169 iz = RAT(iwork, k_iz);
170 ipc = RAT(iwork, k_ipc);
172 rpc = RAT(rwork, k_rpc);
173 pc = RAT(rwork, k_pc);
174 ac = RAT(rwork, k_ac);
175 cc = RAT(rwork, k_cc);
176 fc = RAT(rwork, k_fc);
191 int *nx,
int *ny,
int *nz,
193 int *iz,
int *ipc,
double *rpc,
194 double *pc,
double *ac,
double *cc,
double *fc,
195 double *xf,
double *yf,
double *zf,
196 double *gxcf,
double *gycf,
double *gzcf,
197 double *a1cf,
double *a2cf,
double *a3cf,
198 double *ccf,
double *fcf,
double *tcf) {
232 double tsetupf = 0.0;
233 double tsetupc = 0.0;
244 double errtol_p = 0.0;
246 double rho_min = 0.0;
247 double rho_max = 0.0;
248 double rho_min_mod = 0.0;
249 double rho_max_mod = 0.0;
266 int nlev = VAT(iparm, 6);
272 mgkey = VAT(iparm, 9);
273 itmax = VAT(iparm, 10);
274 istop = VAT(iparm, 11);
275 iinfo = VAT(iparm, 12);
276 ipkey = VAT(iparm, 14);
277 mode = VAT(iparm, 16);
278 mgprol = VAT(iparm, 17);
279 mgcoar = VAT(iparm, 18);
280 mgdisc = VAT(iparm, 19);
281 mgsmoo = VAT(iparm, 20);
282 mgsolv = VAT(iparm, 21);
283 iperf = VAT(iparm, 22);
286 errtol = VAT(rparm, 1);
287 omegal = VAT(rparm, 9);
288 omegan = VAT(rparm, 10);
291 Vprtstp(0, -99, 0.0, 0.0, 0.0);
297 Vnm_tstart(30,
"Vmgdrv2: fine problem setup");
302 &nlev, &ipkey, &iinfo, &ido, iz,
303 &mgprol, &mgcoar, &mgsolv, &mgdisc,
304 ipc, rpc, pc, ac, cc, fc,
311 Vnm_tstop(30,
"Vmgdrv2: fine problem setup");
314 Vnm_tstart(30,
"Vmgdrv2: coarse problem setup");
319 &nlev, &ipkey, &iinfo, &ido, iz,
320 &mgprol, &mgcoar, &mgsolv, &mgdisc,
321 ipc, rpc, pc, ac, cc, fc,
328 Vnm_tstop(30,
"Vmgdrv2: coarse problem setup");
331 epsiln = Vnm_epsmac();
350 for (level=1; level <= nlev_real; level++) {
351 nlevd = nlev_real - level + 1;
358 Vmkcors(&numlev, &nxf, &nyf, &nzf, &nxc, &nyc, &nzc);
367 VMESSAGE3(
"Analysis ==> (%3d, %3d, %3d)", nxf, nyf, nzf);
372 if (iperf == 1 || iperf == 3) {
375 VMESSAGE0(
"Power calculating rho(A)");
385 a1cf, a2cf, a3cf, ccf,
386 &rho_max, &rho_max_mod, &errtol_p,
387 &itmax_p, &iters_p, &iinfo_p);
390 VMESSAGE1(
"Power iters = %d", iters_p);
391 VMESSAGE1(
"Power eigmax = %f", rho_max);
392 VMESSAGE1(
"Power (MODEL) = %f", rho_max_mod);
397 VMESSAGE0(
"Ipower calculating lambda_min(A)...");
406 Vipower(&nxf, &nyf, &nzf, u, iz,
407 a1cf, a2cf, a3cf, ccf, fcf,
408 &rho_min, &rho_min_mod, &errtol_p, &itmax_p, &iters_p,
409 &nlevd, &level, &nlev_real, &mgsolv,
410 &iok_p, &iinfo_p, &epsiln, &errtol, &omegal,
412 ipc, rpc, pc, ac, cc, tcf);
415 VMESSAGE1(
"Ipower iters = %d", iters_p);
416 VMESSAGE1(
"Ipower eigmin = %f", rho_min);
417 VMESSAGE1(
"Ipower (MODEL) = %f", rho_min_mod);
420 VMESSAGE1(
"Condition number = %f", rho_max / rho_min);
421 VMESSAGE1(
"Condition (MODEL) = %f", rho_max_mod / rho_min_mod);
427 if (iperf == 2 || iperf == 3) {
430 VMESSAGE0(
"Mpower calculating rho(M)");
437 Vazeros(&nxf, &nyf, &nzf, RAT(u, VAT2(iz, 1, level)));
440 Vmpower(&nxf, &nyf, &nzf, u, iz,
441 a1cf, a2cf, a3cf, ccf, fcf,
442 &rho_p, &errtol_p, &itmax_p, &iters_p,
443 &nlevd, &level, &nlev_real, &mgsolv,
444 &iok_p, &iinfo_p, &epsiln,
445 &errtol, &omegal, &nu1, &nu2, &mgsmoo,
446 ipc, rpc, pc, ac, cc, fc, tcf);
449 VMESSAGE1(
"Mpower iters = %d", iters_p);
450 VMESSAGE1(
"Mpower rho(M) = %f", rho_p);
456 Vazeros(&nxf, &nyf, &nzf, RAT(u, VAT2(iz, 1, level)));
468 if (istop == 4 || istop == 5 || iperf != 0 ) {
471 VMESSAGE0(
"Generating algebraic RHS from your soln...");
476 Vbuildalg(nx, ny, nz, &mode, &nlev, iz,
477 ipc, rpc, ac, cc, ccf, tcf, fc, fcf);
486 Vnm_tstart(30,
"Vmgdrv2: solve");
489 if (mode == 0 || mode == 2) {
497 u, iz, a1cf, a2cf, a3cf, ccf,
498 &istop, &itmax, &iters, &ierror, &nlev,
499 &ilev, &nlev_real, &mgsolv,
500 &iok, &iinfo, &epsiln, &errtol, &omegal,
502 ipc, rpc, pc, ac, cc, fc, tcf);
504 }
else if (mgkey == 1) {
507 u, iz, a1cf, a2cf, a3cf, ccf,
508 &istop, &itmax, &iters, &ierror, &nlev,
509 &ilev, &nlev_real, &mgsolv,
510 &iok, &iinfo, &epsiln, &errtol, &omegal,
512 ipc, rpc, pc, ac, cc, fc, tcf);
515 VABORT_MSG1(
"Bad mgkey given: %d", mgkey);
519 if (mode == 1 || mode == 2) {
528 u, iz, a1cf, a2cf, a3cf, ccf, fcf,
529 &istop, &itmax, &iters, &ierror, &nlev,
530 &ilev, &nlev_real, &mgsolv,
531 &iok, &iinfo, &epsiln, &errtol, &omegan,
533 ipc, rpc, pc, ac, cc, fc, tcf);
535 }
else if (mgkey == 1) {
539 a1cf, a2cf, a3cf, ccf, fcf,
540 &istop, &itmax, &iters, &ierror, &nlev,
541 &ilev, &nlev_real, &mgsolv,
542 &iok, &iinfo, &epsiln, &errtol, &omegan,
544 ipc, rpc, pc, ac, cc, fc, tcf);
547 VABORT_MSG1(
"Bad mgkey given: %d", mgkey);
552 Vnm_tstop(30,
"Vmgdrv2: solve");
557 VfboundPMG(&ibound, nx, ny, nz, u, gxcf, gycf, gzcf);
562VPUBLIC
void Vmgsz(
int *mgcoar,
int *mgdisc,
int *mgsolv,
563 int *nx,
int *ny,
int *nz,
565 int *nxc,
int *nyc,
int *nzc,
567 int *narr,
int *narrc,
568 int *n_rpc,
int *n_iz,
int *n_ipc,
569 int *iretot,
int *iintot) {
577 int nc_band, num_band, n_band;
580 int num_nf_oper, num_narrc_oper;
586 *nf = *nx * *ny * *nz;
598 for (level=2; level<=*nlev; level++) {
603 Vmkcors(&numlev, &nxf, &nyf, &nzf, nxc, nyc, nzc);
611 *narr += nxf * nyf * nzf;
613 *nc = *nxc * *nyc * *nzc;
614 *narrc = *narr - *nf;
619 }
else if (*mgdisc == 1) {
622 Vnm_print(2,
"Vmgsz: invalid mgdisc parameter: %d\n", *mgdisc);
626 if ((*mgcoar == 0 || *mgcoar == 1) && *mgdisc == 0) {
628 }
else if (*mgcoar == 2) {
631 Vnm_print(2,
"Vmgsz: invalid mgcoar parameter: %d\n", *mgcoar);
637 }
else if (*mgsolv == 1) {
638 if ((*mgcoar == 0 || *mgcoar == 1) && *mgdisc == 0) {
639 num_band = 1 + (*nxc - 2) * (*nyc - 2);
641 num_band = 1 + (*nxc - 2) * (*nyc - 2) + (*nxc - 2) + 1;
643 nc_band = (*nxc - 2) * (*nyc - 2) * (*nzc - 2);
644 n_band = nc_band * num_band;
646 Vnm_print(2,
"Vmgsz: invalid mgsolv parameter: %d\n", *mgsolv);
650 *n_rpc = 100 * (*nlev + 1);
653 *iretot = num_narr * *narr
654 + (num_nf + num_nf_oper) * *nf
655 + (num_narrc + num_narrc_oper) * *narrc
660 *n_iz = 50 * (*nlev + 1);
661 *n_ipc = 100 * (*nlev + 1);
664 *iintot = *n_iz + *n_ipc;
VPUBLIC void Vfmvfas(int *nx, int *ny, int *nz, double *x, int *iz, double *w0, double *w1, double *w2, double *w3, double *w4, 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)
Multigrid nonlinear solve iteration routine.
VPUBLIC void VfboundPMG(int *ibound, int *nx, int *ny, int *nz, double *x, double *gxc, double *gyc, double *gzc)
Initialize a grid function to have a certain boundary value,.
VPUBLIC void Vbuildops(int *nx, int *ny, int *nz, int *nlev, int *ipkey, int *iinfo, int *ido, int *iz, int *mgprol, int *mgcoar, int *mgsolv, int *mgdisc, int *ipc, double *rpc, double *pc, double *ac, double *cc, double *fc, double *xf, double *yf, double *zf, double *gxcf, double *gycf, double *gzcf, double *a1cf, double *a2cf, double *a3cf, double *ccf, double *fcf, double *tcf)
Build operators, boundary arrays, modify affine vectors ido==0: do only fine level ido==1: do only co...
VPUBLIC void Vipower(int *nx, int *ny, int *nz, double *u, int *iz, double *w0, double *w1, double *w2, double *w3, double *w4, double *eigmin, double *eigmin_model, double *tol, int *itmax, int *iters, 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 *tru)
Standard inverse power method for minimum eigenvalue estimation.
VPUBLIC void Vmvfas(int *nx, int *ny, int *nz, double *x, int *iz, double *w0, double *w1, double *w2, double *w3, double *w4, 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)
Nonlinear multilevel method.
VPUBLIC void Vbuildstr(int *nx, int *ny, int *nz, int *nlev, int *iz)
Build the nexted operator framework in the array iz.
VPUBLIC void Vmgsz(int *mgcoar, int *mgdisc, int *mgsolv, int *nx, int *ny, int *nz, int *nlev, int *nxc, int *nyc, int *nzc, int *nf, int *nc, int *narr, int *narrc, int *n_rpc, int *n_iz, int *n_ipc, int *iretot, int *iintot)
This routine computes the required sizes of the real and integer work arrays for the multigrid code....
VPUBLIC void VfboundPMG00(int *nx, int *ny, int *nz, double *x)
Initialize a grid function to have a zero boundary value.
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.
VPUBLIC void Vazeros(int *nx, int *ny, int *nz, double *x)
Zero out operation for a grid function, including boundary values.
VPUBLIC void Vpower(int *nx, int *ny, int *nz, int *iz, int *ilev, int *ipc, double *rpc, double *ac, double *cc, double *w1, double *w2, double *w3, double *w4, double *eigmax, double *eigmax_model, double *tol, int *itmax, int *iters, int *iinfo)
Power methods for eigenvalue estimation.
VPUBLIC void Vmgdriv2(int *iparm, double *rparm, int *nx, int *ny, int *nz, double *u, int *iz, int *ipc, double *rpc, double *pc, double *ac, double *cc, double *fc, double *xf, double *yf, double *zf, double *gxcf, double *gycf, double *gzcf, double *a1cf, double *a2cf, double *a3cf, double *ccf, double *fcf, double *tcf)
Solves the pde using the multi-grid method.
VPUBLIC void Vmgdriv(int *iparm, double *rparm, int *iwork, double *rwork, double *u, double *xf, double *yf, double *zf, double *gxcf, double *gycf, double *gzcf, double *a1cf, double *a2cf, double *a3cf, double *ccf, double *fcf, double *tcf)
Multilevel solver driver.