APBS 3.0.0
Loading...
Searching...
No Matches
powerd.c
1
55#include "powerd.h"
56
57VPUBLIC void Vpower(int *nx, int *ny, int *nz,
58 int *iz, int *ilev,
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) {
63
64 int lev, level;
65 double denom, fac, rho, oldrho, error, relerr;
66
68 double pi = 4.0 * atan( 1.0 );
69
70 // Utility variables
71 int skipIters = 0;
72 double alpha;
73
74 MAT2(iz, 50, 1);
75
76 WARN_UNTESTED;
77
78 // Recover level information
79 level = 1;
80 lev = (*ilev - 1) + level;
81
82 // Seed vector: random to contain all components
83
84 Vaxrand(nx, ny, nz, w1);
85
86 Vazeros(nx, ny, nz, w2);
87 Vazeros(nx, ny, nz, w3);
88 Vazeros(nx, ny, nz, w4);
89
90 // Compute raleigh quotient with the seed vector
91 denom = Vxnrm2(nx, ny, nz, w1);
92 fac = 1.0 / denom;
93 Vxscal(nx, ny, nz, &fac, w1);
94 Vmatvec(nx, ny, nz,
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);
98
99 // I/O
100 if (oldrho == 0.0) {
101 if (*iinfo > 3) {
102 Vnm_print(2, "POWER: iter: estimate = %d %g\n", *iters, oldrho);
103 }
104 rho = oldrho;
105 } else {
106
107 // Main iteration
108 *iters = 0;
109 while(1) {
110 (*iters)++;
111
112 // Apply the matrix A
113 Vmatvec(nx, ny, nz,
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);
116
117 Vxcopy(nx, ny, nz, w2, w1);
118
119 // Normalize the new vector
120 denom = Vxnrm2(nx, ny, nz, w1);
121 fac = 1.0 / denom;
122 Vxscal(nx, ny, nz, &fac, w1);
123
124 // Compute the new raleigh quotient
125 Vmatvec(nx, ny, nz,
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);
129
130 // Stopping test ***
131 // w2=A*x, w1=x, stop = 2-norm(A*x-lamda*x)
132
133 Vxcopy(nx, ny, nz, w1, w3);
134 Vxcopy(nx, ny, nz, w2, w4);
135 Vxscal(nx, ny, nz, &rho, w3);
136 alpha = -1.0;
137 Vxaxpy(nx, ny, nz, &alpha, w3, w4);
138 error = Vxnrm2(nx, ny, nz, w4);
139 relerr = VABS(rho - oldrho ) / VABS( rho );
140
141 // I/O
142 if (*iinfo > 3) {
143
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);
148 }
149
150 if( relerr < *tol || *iters == *itmax)
151 break;
152
153 oldrho = rho;
154 }
155 }
156
157 // Return some stuff ***
158 *eigmax = rho;
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)));
162}
163
164
165VPUBLIC void Vipower(int *nx,int *ny,int *nz,
166 double *u, int *iz,
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) {
175
176 int level, lev;
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;
181
183 double pi = 4.0 * atan( 1.0 );
184
185 // Utility variables
186 double alpha;
187
188 MAT2(iz, 50, 1);
189
190 WARN_UNTESTED;
191
192 // Recover level information
193 level = 1;
194 lev = (*ilev - 1) + level;
195
196 // Seed vector: random to contain all components
197 Vaxrand(nx, ny, nz, w1);
198 Vazeros(nx, ny, nz, w2);
199 Vazeros(nx, ny, nz, w3);
200 Vazeros(nx, ny, nz, w4);
201 Vazeros(nx, ny, nz, RAT(w0, VAT2(iz, 1, lev)));
202 Vazeros(nx, ny, nz, RAT( u, VAT2(iz, 1, lev)));
203
204 // Compute raleigh quotient with the seed vector ***
205 denom = Vxnrm2(nx, ny, nz, w1);
206 fac = 1.0 / denom;
207 Vxscal(nx, ny, nz, &fac, w1);
208 Vmatvec(nx, ny, nz,
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);
212
213 // I/O
214 if (oldrho == 0.0) {
215 if (*iinfo > 3) {
216 Vnm_print(2, "Vipower: iters=%d\n", *iters);
217 Vnm_print(2, " estimate=%f\n", oldrho);
218 }
219 rho = oldrho;
220 } else {
221
222 //main iteration
223 *iters = 0;
224 while (1) {
225 (*iters)++;
226
227 // Apply the matrix A^{-1} (using MG solver)
228 itmax_s = 100;
229 iters_s = 0;
230 ierror_s = 0;
231 iok_s = 0;
232 iinfo_s = 0;
233 istop_s = 0;
234 mgsmoo_s = 1;
235 nu1_s = 1;
236 nu2_s = 1;
237 errtol_s = *epsiln;
238
239 Vxcopy(nx, ny, nz, w1, RAT(w0, VAT2(iz, 1,lev)));
240 Vmvcs(nx, ny, nz, u, iz,
241 w1, w2, w3, w4,
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);
248
249 // Normalize the new vector
250 denom = Vxnrm2(nx, ny, nz, w1);
251 fac = 1.0 / denom;
252 Vxscal(nx, ny, nz, &fac, w1);
253
254 // Compute the new raleigh quotient
255 Vmatvec(nx, ny, nz,
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);
259
260 // Stopping test
261 // w2=A*x, w1=x, stop = 2-norm(A*x-lamda*x) ***
262 Vxcopy(nx, ny, nz, w1, w3);
263 Vxcopy(nx, ny, nz, w2, w4);
264 Vxscal(nx, ny, nz, &rho, w3);
265 alpha = -1.0;
266 Vxaxpy(nx, ny, nz, &alpha, w3, w4);
267 error = Vxnrm2(nx, ny, nz, w4);
268 relerr = VABS(rho - oldrho ) / VABS( rho );
269
270 // I/O
271 if (*iinfo > 3) {
272
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);
277 }
278
279 if (relerr < *tol || *iters == *itmax)
280 break;
281
282 oldrho = rho;
283 }
284 }
285
286 // Return some stuff
287 *eigmin = rho;
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)));
292}
293
294VEXTERNC void Vmpower(int *nx, int *ny, int *nz,
295 double *u, int *iz,
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) {
303
304 // Local variables
305 int lev, level;
306 double denom, fac, rho, oldrho, error;
307 double relerr;
308 int itmax_s, iters_s, ierror_s, iok_s, iinfo_s, istop_s;
309 double alpha;
310
311 MAT2(iz, 50, 1);
312
313 // Recover level information
314 level = 1;
315 lev = (*ilev - 1) + level;
316
317 // Seed vector: random to contain all components
318 Vaxrand(nx, ny, nz, w1);
319 Vazeros(nx, ny, nz, w2);
320 Vazeros(nx, ny, nz, w3);
321 Vazeros(nx, ny, nz, w4);
322 Vazeros(nx, ny, nz, RAT(u, VAT2(iz, 1, lev)));
323
324 // NOTE: we destroy "fc" on this level due to lack of vectors... ***
325 Vazeros(nx,ny,nz,RAT(fc, VAT2(iz, 1, lev)));
326
327 // Normalize the seed vector
328 denom = Vxnrm2(nx, ny, nz, w1);
329 fac = 1.0 / denom;
330 Vxscal(nx, ny, nz, &fac, w1);
331
332 // Compute raleigh quotient with the seed vector
333 Vxcopy(nx, ny, nz, w1, RAT(u, VAT2(iz, 1, lev)));
334 itmax_s = 1;
335 iters_s = 0;
336 ierror_s = 0;
337 iok_s = 0;
338 iinfo_s = 0;
339 istop_s = 1;
340 Vmvcs(nx, ny, nz,
341 u, iz, w0, w2, w3, w4,
342 &istop_s, &itmax_s, &iters_s, &ierror_s,
343 nlev, ilev, nlev_real, mgsolv,
344 &iok_s, &iinfo_s,
345 epsiln, errtol, omega, nu1, nu2, mgsmoo,
346 ipc, rpc,
347 pc, ac, cc, fc, tru);
348 oldrho = Vxdot(nx, ny, nz, w1, RAT(u, VAT2(iz, 1, lev)));
349
350 // I/O
351 if (oldrho == 0.0) {
352 if (*iinfo > 3) {
353 Vnm_print(2, "Vmp0ower: iter=%d, estimate=%f", *iters, oldrho);
354 }
355 rho = oldrho;
356
357 } else {
358
359 // Main iteration
360 *iters = 0;
361 while (1) {
362 (*iters)++;
363
364 // Apply the matrix M
365 Vxcopy(nx, ny, nz, w1, RAT(u, VAT2(iz, 1, lev)));
366 itmax_s = 1;
367 iters_s = 0;
368 ierror_s = 0;
369 iok_s = 0;
370 iinfo_s = 0;
371 istop_s = 1;
372 Vmvcs(nx, ny, nz,
373 u, iz, w1, w2, w3, w4,
374 &istop_s, &itmax_s, &iters_s, &ierror_s,
375 nlev, ilev, nlev_real, mgsolv,
376 &iok_s, &iinfo_s,
377 epsiln, errtol, omega, nu1, nu2, mgsmoo,
378 ipc, rpc,
379 pc, ac, cc, fc, tru);
380 Vxcopy(nx, ny, nz, RAT(u, VAT2(iz, 1, lev)), w1);
381
382 // Normalize the new vector
383 denom = Vxnrm2(nx, ny, nz, w1);
384 fac = 1.0 / denom;
385 Vxscal(nx, ny, nz, &fac, w1);
386
387 // Compute the new raleigh quotient
388 Vxcopy(nx, ny, nz, w1, RAT(u, VAT2(iz, 1, lev)));
389 itmax_s = 1;
390 iters_s = 0;
391 ierror_s = 0;
392 iok_s = 0;
393 iinfo_s = 0;
394 istop_s = 1;
395 Vmvcs(nx, ny, nz,
396 u, iz, w0, w2, w3, w4,
397 &istop_s, &itmax_s, &iters_s, &ierror_s,
398 nlev, ilev, nlev_real, mgsolv,
399 &iok_s, &iinfo_s,
400 epsiln, errtol, omega, nu1, nu2, mgsmoo,
401 ipc, rpc,
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);
405
406 // Stopping test
407 // w2=A*x, w1=x, stop = 2-norm(A*x-lamda*x)
408 alpha = -1.0;
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 );
415
416 // I/O
417 if (*iinfo > 3) {
418 Vnm_print(2, "Vmpower: iter=%d; error=%f; relerr=%f; estimate=%f",
419 *iters, error, relerr, rho);
420 }
421
422 if ((relerr < *tol) || (*iters == *itmax)) {
423 break;
424 }
425
426 oldrho = rho;
427 }
428 }
429
430 *eigmax = rho;
431}
VPUBLIC void Vaxrand(int *nx, int *ny, int *nz, double *x)
Fill grid function with random values, including boundary values.
Definition mikpckd.c:291
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.
Definition powerd.c:165
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.
Definition matvecd.c:57
VPUBLIC double Vxdot(int *nx, int *ny, int *nz, double *x, double *y)
Inner product operation for a grid function with boundary values.
Definition mikpckd.c:173
VPUBLIC void Vxscal(int *nx, int *ny, int *nz, double *fac, double *x)
Scale operation for a grid function with boundary values.
Definition mikpckd.c:318
VPUBLIC void Vxcopy(int *nx, int *ny, int *nz, double *x, double *y)
A collection of useful low-level routines (timing, etc).
Definition mikpckd.c:57
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.
Definition mgcsd.c:57
VPUBLIC void Vxaxpy(int *nx, int *ny, int *nz, double *alpha, double *x, double *y)
saxpy operation for a grid function with boundary values.
Definition mikpckd.c:112
VPUBLIC double Vxnrm2(int *nx, int *ny, int *nz, double *x)
Norm operation for a grid function with boundary values.
Definition mikpckd.c:152
VPUBLIC void Vazeros(int *nx, int *ny, int *nz, double *x)
Zero out operation for a grid function, including boundary values.
Definition mikpckd.c:195
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.
Definition powerd.c:57