54 int *ipc,
double *rpc,
double *ac,
double *cc,
55 double *w1,
double *w2,
double *w3,
double *w4,
56 double *eigmax,
double *eigmax_model,
double *tol,
60 double denom, fac, rho, oldrho, error, relerr;
63 double pi = 4.0 * atan( 1.0 );
75 lev = (*ilev - 1) + level;
86 denom =
Vxnrm2(nx, ny, nz, w1);
88 Vxscal(nx, ny, nz, &fac, w1);
90 RAT(ipc, VAT2(iz, 5, lev)), RAT(rpc, VAT2(iz, 6,lev)),
91 RAT(ac, VAT2(iz, 7, lev)), RAT(cc, VAT2(iz, 1,lev)), w1, w2);
92 oldrho =
Vxdot(nx, ny, nz, w1, w2);
97 Vnm_print(2,
"POWER: iter: estimate = %d %g\n", *iters, oldrho);
109 RAT(ipc, VAT2(iz, 5, lev)), RAT(rpc, VAT2(iz, 6, lev)),
110 RAT(ac, VAT2(iz, 7, lev)), RAT(cc, VAT2(iz, 1, lev)), w1, w2);
112 Vxcopy(nx, ny, nz, w2, w1);
115 denom =
Vxnrm2(nx, ny, nz, w1);
117 Vxscal(nx, ny, nz, &fac, w1);
121 RAT(ipc, VAT2(iz, 5, lev)), RAT(rpc, VAT2(iz, 6, lev)),
122 RAT( ac, VAT2(iz, 7, lev)), RAT( cc, VAT2(iz, 1, lev)), w1, w2);
123 rho =
Vxdot(nx, ny, nz, w1, w2);
128 Vxcopy(nx, ny, nz, w1, w3);
129 Vxcopy(nx, ny, nz, w2, w4);
130 Vxscal(nx, ny, nz, &rho, w3);
132 Vxaxpy(nx, ny, nz, &alpha, w3, w4);
133 error =
Vxnrm2(nx, ny, nz, w4);
134 relerr = VABS(rho - oldrho ) / VABS( rho );
139 Vnm_print(2,
"POWER: iters =%d\n", *iters);
140 Vnm_print(2,
" error =%g\n", error);
141 Vnm_print(2,
" relerr =%g\n", relerr);
142 Vnm_print(2,
" rho =%g\n", rho);
145 if( relerr < *tol || *iters == *itmax)
154 fac = VPOW(2.0, *ilev - 1);
155 *eigmax_model = fac * (6.0 - 2.0 * VCOS((*nx - 2) * pi / (*nx - 1))
156 - 2.0 * VCOS((*ny - 2) * pi / (*ny - 1)));
162 double *w0,
double *w1,
double *w2,
double *w3,
double *w4,
163 double *eigmin,
double *eigmin_model,
double *tol,
164 int *
itmax,
int *iters,
165 int *
nlev,
int *ilev,
int *nlev_real,
int *
mgsolv,
166 int *iok,
int *
iinfo,
double *epsiln,
double *
errtol,
double *omega,
168 int *ipc,
double *rpc,
169 double *pc,
double *ac,
double *cc,
double *tru) {
172 double denom, fac, rho, oldrho;
173 double error, relerr, errtol_s;
174 int itmax_s, iters_s, ierror_s, iok_s, iinfo_s, istop_s;
175 int nu1_s, nu2_s, mgsmoo_s;
178 double pi = 4.0 * atan( 1.0 );
189 lev = (*ilev - 1) + level;
196 Vazeros(nx, ny, nz, RAT(w0, VAT2(iz, 1, lev)));
197 Vazeros(nx, ny, nz, RAT( u, VAT2(iz, 1, lev)));
200 denom =
Vxnrm2(nx, ny, nz, w1);
202 Vxscal(nx, ny, nz, &fac, w1);
204 RAT(ipc, VAT2(iz, 5, lev)), RAT(rpc, VAT2(iz, 6, lev)),
205 RAT( ac, VAT2(iz, 7, lev)), RAT( cc, VAT2(iz, 1, lev)), w1, w2);
206 oldrho =
Vxdot(nx, ny, nz, w1, w2);
211 Vnm_print(2,
"Vipower: iters=%d\n", *iters);
212 Vnm_print(2,
" estimate=%f\n", oldrho);
234 Vxcopy(nx, ny, nz, w1, RAT(w0, VAT2(iz, 1,lev)));
235 Vmvcs(nx, ny, nz, u, iz,
237 &istop_s, &itmax_s, &iters_s, &ierror_s,
238 nlev, ilev, nlev_real, mgsolv,
239 &iok_s, &iinfo_s, epsiln,
240 &errtol_s, omega, &nu1_s, &nu2_s, &mgsmoo_s,
241 ipc, rpc, pc, ac, cc, w0, tru);
242 Vxcopy(nx, ny, nz, RAT(u, VAT2(iz, 1, lev)), w1);
245 denom =
Vxnrm2(nx, ny, nz, w1);
247 Vxscal(nx, ny, nz, &fac, w1);
251 RAT(ipc, VAT2(iz, 5, lev)), RAT(rpc, VAT2(iz, 6, lev)),
252 RAT(ac, VAT2(iz, 7,lev)), RAT(cc, VAT2(iz, 1, lev)), w1, w2);
253 rho =
Vxdot(nx, ny, nz, w1, w2);
257 Vxcopy(nx, ny, nz, w1, w3);
258 Vxcopy(nx, ny, nz, w2, w4);
259 Vxscal(nx, ny, nz, &rho, w3);
261 Vxaxpy(nx, ny, nz, &alpha, w3, w4);
262 error =
Vxnrm2(nx, ny, nz, w4);
263 relerr = VABS(rho - oldrho ) / VABS( rho );
268 Vnm_print(2,
"POWER: iters =%d\n", *iters);
269 Vnm_print(2,
" error =%g\n", error);
270 Vnm_print(2,
" relerr =%g\n", relerr);
271 Vnm_print(2,
" rho =%g\n", rho);
274 if (relerr < *tol || *iters == *itmax)
283 fac = VPOW(2.0, *ilev - 1);
284 *eigmin_model = fac * (6.0 - 2.0 * VCOS(pi / (*nx - 1))
285 - 2.0 * VCOS(pi / (*ny - 1))
286 - 2.0 * VCOS(pi / (*nz - 1)));
289 VEXTERNC
void Vmpower(
int *
nx,
int *
ny,
int *
nz,
291 double *w0,
double *w1,
double *w2,
double *w3,
double *w4,
292 double *eigmax,
double *tol,
293 int *
itmax,
int *iters,
int *
nlev,
int *ilev,
int *nlev_real,
295 double *epsiln,
double *
errtol,
double *omega,
297 double *pc,
double *ac,
double *cc,
double *fc,
double *tru) {
301 double denom, fac, rho, oldrho, error;
303 int itmax_s, iters_s, ierror_s, iok_s, iinfo_s, istop_s;
310 lev = (*ilev - 1) + level;
336 u, iz, w0, w2, w3, w4,
337 &istop_s, &itmax_s, &iters_s, &ierror_s,
342 pc, ac, cc, fc, tru);
343 oldrho =
Vxdot(
nx,
ny,
nz, w1, RAT(u, VAT2(iz, 1, lev)));
348 Vnm_print(2,
"Vmp0ower: iter=%d, estimate=%f", *iters, oldrho);
368 u, iz, w1, w2, w3, w4,
369 &istop_s, &itmax_s, &iters_s, &ierror_s,
374 pc, ac, cc, fc, tru);
391 u, iz, w0, w2, w3, w4,
392 &istop_s, &itmax_s, &iters_s, &ierror_s,
397 pc, ac, cc, fc, tru);
409 relerr = VABS( rho - oldrho ) / VABS( rho );
413 Vnm_print(2,
"Vmpower: iter=%d; error=%f; relerr=%f; estimate=%f",
414 *iters, error, relerr, rho);
417 if ((relerr < *tol) || (*iters == *
itmax)) {
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 Vxscal(int *nx, int *ny, int *nz, double *fac, double *x)
Scale operation for a grid function with boundary values.
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 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 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 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 Vaxrand(int *nx, int *ny, int *nz, double *x)
Fill grid function with random values, including boundary values.
VPUBLIC double Vxnrm2(int *nx, int *ny, int *nz, double *x)
Norm operation for a grid function with boundary values.