57VPUBLIC
void Vpower(
int *nx,
int *ny,
int *nz,
59 int *ipc,
double *rpc,
double *ac,
double *cc,
60 double *w1,
double *w2,
double *w3,
double *w4,
61 double *eigmax,
double *eigmax_model,
double *tol,
62 int *itmax,
int *iters,
int *iinfo) {
65 double denom, fac, rho, oldrho, error, relerr;
68 double pi = 4.0 * atan( 1.0 );
80 lev = (*ilev - 1) + level;
91 denom =
Vxnrm2(nx, ny, nz, w1);
93 Vxscal(nx, ny, nz, &fac, w1);
95 RAT(ipc, VAT2(iz, 5, lev)), RAT(rpc, VAT2(iz, 6,lev)),
96 RAT(ac, VAT2(iz, 7, lev)), RAT(cc, VAT2(iz, 1,lev)), w1, w2);
97 oldrho =
Vxdot(nx, ny, nz, w1, w2);
102 Vnm_print(2,
"POWER: iter: estimate = %d %g\n", *iters, oldrho);
114 RAT(ipc, VAT2(iz, 5, lev)), RAT(rpc, VAT2(iz, 6, lev)),
115 RAT(ac, VAT2(iz, 7, lev)), RAT(cc, VAT2(iz, 1, lev)), w1, w2);
117 Vxcopy(nx, ny, nz, w2, w1);
120 denom =
Vxnrm2(nx, ny, nz, w1);
122 Vxscal(nx, ny, nz, &fac, w1);
126 RAT(ipc, VAT2(iz, 5, lev)), RAT(rpc, VAT2(iz, 6, lev)),
127 RAT( ac, VAT2(iz, 7, lev)), RAT( cc, VAT2(iz, 1, lev)), w1, w2);
128 rho =
Vxdot(nx, ny, nz, w1, w2);
133 Vxcopy(nx, ny, nz, w1, w3);
134 Vxcopy(nx, ny, nz, w2, w4);
135 Vxscal(nx, ny, nz, &rho, w3);
137 Vxaxpy(nx, ny, nz, &alpha, w3, w4);
138 error =
Vxnrm2(nx, ny, nz, w4);
139 relerr = VABS(rho - oldrho ) / VABS( rho );
144 Vnm_print(2,
"POWER: iters =%d\n", *iters);
145 Vnm_print(2,
" error =%g\n", error);
146 Vnm_print(2,
" relerr =%g\n", relerr);
147 Vnm_print(2,
" rho =%g\n", rho);
150 if( relerr < *tol || *iters == *itmax)
159 fac = VPOW(2.0, *ilev - 1);
160 *eigmax_model = fac * (6.0 - 2.0 * VCOS((*nx - 2) * pi / (*nx - 1))
161 - 2.0 * VCOS((*ny - 2) * pi / (*ny - 1)));
167 double *w0,
double *w1,
double *w2,
double *w3,
double *w4,
168 double *eigmin,
double *eigmin_model,
double *tol,
169 int *itmax,
int *iters,
170 int *nlev,
int *ilev,
int *nlev_real,
int *mgsolv,
171 int *iok,
int *iinfo,
double *epsiln,
double *errtol,
double *omega,
172 int *nu1,
int *nu2,
int *mgsmoo,
173 int *ipc,
double *rpc,
174 double *pc,
double *ac,
double *cc,
double *tru) {
177 double denom, fac, rho, oldrho;
178 double error, relerr, errtol_s;
179 int itmax_s, iters_s, ierror_s, iok_s, iinfo_s, istop_s;
180 int nu1_s, nu2_s, mgsmoo_s;
183 double pi = 4.0 * atan( 1.0 );
194 lev = (*ilev - 1) + level;
201 Vazeros(nx, ny, nz, RAT(w0, VAT2(iz, 1, lev)));
202 Vazeros(nx, ny, nz, RAT( u, VAT2(iz, 1, lev)));
205 denom =
Vxnrm2(nx, ny, nz, w1);
207 Vxscal(nx, ny, nz, &fac, w1);
209 RAT(ipc, VAT2(iz, 5, lev)), RAT(rpc, VAT2(iz, 6, lev)),
210 RAT( ac, VAT2(iz, 7, lev)), RAT( cc, VAT2(iz, 1, lev)), w1, w2);
211 oldrho =
Vxdot(nx, ny, nz, w1, w2);
216 Vnm_print(2,
"Vipower: iters=%d\n", *iters);
217 Vnm_print(2,
" estimate=%f\n", oldrho);
239 Vxcopy(nx, ny, nz, w1, RAT(w0, VAT2(iz, 1,lev)));
240 Vmvcs(nx, ny, nz, u, iz,
242 &istop_s, &itmax_s, &iters_s, &ierror_s,
243 nlev, ilev, nlev_real, mgsolv,
244 &iok_s, &iinfo_s, epsiln,
245 &errtol_s, omega, &nu1_s, &nu2_s, &mgsmoo_s,
246 ipc, rpc, pc, ac, cc, w0, tru);
247 Vxcopy(nx, ny, nz, RAT(u, VAT2(iz, 1, lev)), w1);
250 denom =
Vxnrm2(nx, ny, nz, w1);
252 Vxscal(nx, ny, nz, &fac, w1);
256 RAT(ipc, VAT2(iz, 5, lev)), RAT(rpc, VAT2(iz, 6, lev)),
257 RAT(ac, VAT2(iz, 7,lev)), RAT(cc, VAT2(iz, 1, lev)), w1, w2);
258 rho =
Vxdot(nx, ny, nz, w1, w2);
262 Vxcopy(nx, ny, nz, w1, w3);
263 Vxcopy(nx, ny, nz, w2, w4);
264 Vxscal(nx, ny, nz, &rho, w3);
266 Vxaxpy(nx, ny, nz, &alpha, w3, w4);
267 error =
Vxnrm2(nx, ny, nz, w4);
268 relerr = VABS(rho - oldrho ) / VABS( rho );
273 Vnm_print(2,
"POWER: iters =%d\n", *iters);
274 Vnm_print(2,
" error =%g\n", error);
275 Vnm_print(2,
" relerr =%g\n", relerr);
276 Vnm_print(2,
" rho =%g\n", rho);
279 if (relerr < *tol || *iters == *itmax)
288 fac = VPOW(2.0, *ilev - 1);
289 *eigmin_model = fac * (6.0 - 2.0 * VCOS(pi / (*nx - 1))
290 - 2.0 * VCOS(pi / (*ny - 1))
291 - 2.0 * VCOS(pi / (*nz - 1)));
294VEXTERNC
void Vmpower(
int *nx,
int *ny,
int *nz,
296 double *w0,
double *w1,
double *w2,
double *w3,
double *w4,
297 double *eigmax,
double *tol,
298 int *itmax,
int *iters,
int *nlev,
int *ilev,
int *nlev_real,
299 int *mgsolv,
int *iok,
int *iinfo,
300 double *epsiln,
double *errtol,
double *omega,
301 int *nu1,
int *nu2,
int *mgsmoo,
int *ipc,
double *rpc,
302 double *pc,
double *ac,
double *cc,
double *fc,
double *tru) {
306 double denom, fac, rho, oldrho, error;
308 int itmax_s, iters_s, ierror_s, iok_s, iinfo_s, istop_s;
315 lev = (*ilev - 1) + level;
322 Vazeros(nx, ny, nz, RAT(u, VAT2(iz, 1, lev)));
325 Vazeros(nx,ny,nz,RAT(fc, VAT2(iz, 1, lev)));
328 denom =
Vxnrm2(nx, ny, nz, w1);
330 Vxscal(nx, ny, nz, &fac, w1);
333 Vxcopy(nx, ny, nz, w1, RAT(u, VAT2(iz, 1, lev)));
341 u, iz, w0, w2, w3, w4,
342 &istop_s, &itmax_s, &iters_s, &ierror_s,
343 nlev, ilev, nlev_real, mgsolv,
345 epsiln, errtol, omega, nu1, nu2, mgsmoo,
347 pc, ac, cc, fc, tru);
348 oldrho =
Vxdot(nx, ny, nz, w1, RAT(u, VAT2(iz, 1, lev)));
353 Vnm_print(2,
"Vmp0ower: iter=%d, estimate=%f", *iters, oldrho);
365 Vxcopy(nx, ny, nz, w1, RAT(u, VAT2(iz, 1, lev)));
373 u, iz, w1, w2, w3, w4,
374 &istop_s, &itmax_s, &iters_s, &ierror_s,
375 nlev, ilev, nlev_real, mgsolv,
377 epsiln, errtol, omega, nu1, nu2, mgsmoo,
379 pc, ac, cc, fc, tru);
380 Vxcopy(nx, ny, nz, RAT(u, VAT2(iz, 1, lev)), w1);
383 denom =
Vxnrm2(nx, ny, nz, w1);
385 Vxscal(nx, ny, nz, &fac, w1);
388 Vxcopy(nx, ny, nz, w1, RAT(u, VAT2(iz, 1, lev)));
396 u, iz, w0, w2, w3, w4,
397 &istop_s, &itmax_s, &iters_s, &ierror_s,
398 nlev, ilev, nlev_real, mgsolv,
400 epsiln, errtol, omega, nu1, nu2, mgsmoo,
402 pc, ac, cc, fc, tru);
403 Vxcopy(nx, ny, nz, RAT(u, VAT2(iz, 1, lev)), w2);
404 rho =
Vxdot(nx, ny, nz, w1, w2);
409 Vxcopy(nx, ny, nz, w1, w3);
410 Vxcopy(nx, ny, nz, w2, w4);
411 Vxscal(nx, ny, nz, &rho, w3);
412 Vxaxpy(nx, ny, nz, &alpha, w3, w4);
413 error =
Vxnrm2(nx, ny, nz, w4);
414 relerr = VABS( rho - oldrho ) / VABS( rho );
418 Vnm_print(2,
"Vmpower: iter=%d; error=%f; relerr=%f; estimate=%f",
419 *iters, error, relerr, rho);
422 if ((relerr < *tol) || (*iters == *itmax)) {
VPUBLIC void Vaxrand(int *nx, int *ny, int *nz, double *x)
Fill grid function with random values, including 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 double Vxdot(int *nx, int *ny, int *nz, double *x, double *y)
Inner product operation for a grid function with boundary values.
VPUBLIC void Vxscal(int *nx, int *ny, int *nz, double *fac, double *x)
Scale operation for a grid function with boundary values.
VPUBLIC void Vxcopy(int *nx, int *ny, int *nz, double *x, double *y)
A collection of useful low-level routines (timing, etc).
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 Vxaxpy(int *nx, int *ny, int *nz, double *alpha, double *x, double *y)
saxpy operation for a grid function with boundary values.
VPUBLIC double Vxnrm2(int *nx, int *ny, int *nz, double *x)
Norm 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.