APBS  1.5
powerd.c
1 
50 #include "powerd.h"
51 
52 VPUBLIC void Vpower(int *nx, int *ny, int *nz,
53  int *iz, int *ilev,
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,
57  int *itmax, int *iters, int *iinfo) {
58 
59  int lev, level;
60  double denom, fac, rho, oldrho, error, relerr;
61 
63  double pi = 4.0 * atan( 1.0 );
64 
65  // Utility variables
66  int skipIters = 0;
67  double alpha;
68 
69  MAT2(iz, 50, 1);
70 
71  WARN_UNTESTED;
72 
73  // Recover level information
74  level = 1;
75  lev = (*ilev - 1) + level;
76 
77  // Seed vector: random to contain all components
78 
79  Vaxrand(nx, ny, nz, w1);
80 
81  Vazeros(nx, ny, nz, w2);
82  Vazeros(nx, ny, nz, w3);
83  Vazeros(nx, ny, nz, w4);
84 
85  // Compute raleigh quotient with the seed vector
86  denom = Vxnrm2(nx, ny, nz, w1);
87  fac = 1.0 / denom;
88  Vxscal(nx, ny, nz, &fac, w1);
89  Vmatvec(nx, ny, nz,
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);
93 
94  // I/O
95  if (oldrho == 0.0) {
96  if (*iinfo > 3) {
97  Vnm_print(2, "POWER: iter: estimate = %d %g\n", *iters, oldrho);
98  }
99  rho = oldrho;
100  } else {
101 
102  // Main iteration
103  *iters = 0;
104  while(1) {
105  (*iters)++;
106 
107  // Apply the matrix A
108  Vmatvec(nx, ny, nz,
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);
111 
112  Vxcopy(nx, ny, nz, w2, w1);
113 
114  // Normalize the new vector
115  denom = Vxnrm2(nx, ny, nz, w1);
116  fac = 1.0 / denom;
117  Vxscal(nx, ny, nz, &fac, w1);
118 
119  // Compute the new raleigh quotient
120  Vmatvec(nx, ny, nz,
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);
124 
125  // Stopping test ***
126  // w2=A*x, w1=x, stop = 2-norm(A*x-lamda*x)
127 
128  Vxcopy(nx, ny, nz, w1, w3);
129  Vxcopy(nx, ny, nz, w2, w4);
130  Vxscal(nx, ny, nz, &rho, w3);
131  alpha = -1.0;
132  Vxaxpy(nx, ny, nz, &alpha, w3, w4);
133  error = Vxnrm2(nx, ny, nz, w4);
134  relerr = VABS(rho - oldrho ) / VABS( rho );
135 
136  // I/O
137  if (*iinfo > 3) {
138 
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);
143  }
144 
145  if( relerr < *tol || *iters == *itmax)
146  break;
147 
148  oldrho = rho;
149  }
150  }
151 
152  // Return some stuff ***
153  *eigmax = rho;
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)));
157 }
158 
159 
160 VPUBLIC void Vipower(int *nx,int *ny,int *nz,
161  double *u, int *iz,
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,
167  int *nu1, int *nu2, int *mgsmoo,
168  int *ipc, double *rpc,
169  double *pc, double *ac, double *cc, double *tru) {
170 
171  int level, lev;
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;
176 
178  double pi = 4.0 * atan( 1.0 );
179 
180  // Utility variables
181  double alpha;
182 
183  MAT2(iz, 50, 1);
184 
185  WARN_UNTESTED;
186 
187  // Recover level information
188  level = 1;
189  lev = (*ilev - 1) + level;
190 
191  // Seed vector: random to contain all components
192  Vaxrand(nx, ny, nz, w1);
193  Vazeros(nx, ny, nz, w2);
194  Vazeros(nx, ny, nz, w3);
195  Vazeros(nx, ny, nz, w4);
196  Vazeros(nx, ny, nz, RAT(w0, VAT2(iz, 1, lev)));
197  Vazeros(nx, ny, nz, RAT( u, VAT2(iz, 1, lev)));
198 
199  // Compute raleigh quotient with the seed vector ***
200  denom = Vxnrm2(nx, ny, nz, w1);
201  fac = 1.0 / denom;
202  Vxscal(nx, ny, nz, &fac, w1);
203  Vmatvec(nx, ny, nz,
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);
207 
208  // I/O
209  if (oldrho == 0.0) {
210  if (*iinfo > 3) {
211  Vnm_print(2, "Vipower: iters=%d\n", *iters);
212  Vnm_print(2, " estimate=%f\n", oldrho);
213  }
214  rho = oldrho;
215  } else {
216 
217  //main iteration
218  *iters = 0;
219  while (1) {
220  (*iters)++;
221 
222  // Apply the matrix A^{-1} (using MG solver)
223  itmax_s = 100;
224  iters_s = 0;
225  ierror_s = 0;
226  iok_s = 0;
227  iinfo_s = 0;
228  istop_s = 0;
229  mgsmoo_s = 1;
230  nu1_s = 1;
231  nu2_s = 1;
232  errtol_s = *epsiln;
233 
234  Vxcopy(nx, ny, nz, w1, RAT(w0, VAT2(iz, 1,lev)));
235  Vmvcs(nx, ny, nz, u, iz,
236  w1, w2, w3, w4,
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);
243 
244  // Normalize the new vector
245  denom = Vxnrm2(nx, ny, nz, w1);
246  fac = 1.0 / denom;
247  Vxscal(nx, ny, nz, &fac, w1);
248 
249  // Compute the new raleigh quotient
250  Vmatvec(nx, ny, nz,
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);
254 
255  // Stopping test
256  // w2=A*x, w1=x, stop = 2-norm(A*x-lamda*x) ***
257  Vxcopy(nx, ny, nz, w1, w3);
258  Vxcopy(nx, ny, nz, w2, w4);
259  Vxscal(nx, ny, nz, &rho, w3);
260  alpha = -1.0;
261  Vxaxpy(nx, ny, nz, &alpha, w3, w4);
262  error = Vxnrm2(nx, ny, nz, w4);
263  relerr = VABS(rho - oldrho ) / VABS( rho );
264 
265  // I/O
266  if (*iinfo > 3) {
267 
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);
272  }
273 
274  if (relerr < *tol || *iters == *itmax)
275  break;
276 
277  oldrho = rho;
278  }
279  }
280 
281  // Return some stuff
282  *eigmin = rho;
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)));
287 }
288 
289 VEXTERNC void Vmpower(int *nx, int *ny, int *nz,
290  double *u, int *iz,
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,
294  int *mgsolv, int *iok, int *iinfo,
295  double *epsiln, double *errtol, double *omega,
296  int *nu1, int *nu2, int *mgsmoo, int *ipc, double *rpc,
297  double *pc, double *ac, double *cc, double *fc, double *tru) {
298 
299  // Local variables
300  int lev, level;
301  double denom, fac, rho, oldrho, error;
302  double relerr;
303  int itmax_s, iters_s, ierror_s, iok_s, iinfo_s, istop_s;
304  double alpha;
305 
306  MAT2(iz, 50, 1);
307 
308  // Recover level information
309  level = 1;
310  lev = (*ilev - 1) + level;
311 
312  // Seed vector: random to contain all components
313  Vaxrand(nx, ny, nz, w1);
314  Vazeros(nx, ny, nz, w2);
315  Vazeros(nx, ny, nz, w3);
316  Vazeros(nx, ny, nz, w4);
317  Vazeros(nx, ny, nz, RAT(u, VAT2(iz, 1, lev)));
318 
319  // NOTE: we destroy "fc" on this level due to lack of vectors... ***
320  Vazeros(nx,ny,nz,RAT(fc, VAT2(iz, 1, lev)));
321 
322  // Normalize the seed vector
323  denom = Vxnrm2(nx, ny, nz, w1);
324  fac = 1.0 / denom;
325  Vxscal(nx, ny, nz, &fac, w1);
326 
327  // Compute raleigh quotient with the seed vector
328  Vxcopy(nx, ny, nz, w1, RAT(u, VAT2(iz, 1, lev)));
329  itmax_s = 1;
330  iters_s = 0;
331  ierror_s = 0;
332  iok_s = 0;
333  iinfo_s = 0;
334  istop_s = 1;
335  Vmvcs(nx, ny, nz,
336  u, iz, w0, w2, w3, w4,
337  &istop_s, &itmax_s, &iters_s, &ierror_s,
338  nlev, ilev, nlev_real, mgsolv,
339  &iok_s, &iinfo_s,
340  epsiln, errtol, omega, nu1, nu2, mgsmoo,
341  ipc, rpc,
342  pc, ac, cc, fc, tru);
343  oldrho = Vxdot(nx, ny, nz, w1, RAT(u, VAT2(iz, 1, lev)));
344 
345  // I/O
346  if (oldrho == 0.0) {
347  if (*iinfo > 3) {
348  Vnm_print(2, "Vmp0ower: iter=%d, estimate=%f", *iters, oldrho);
349  }
350  rho = oldrho;
351 
352  } else {
353 
354  // Main iteration
355  *iters = 0;
356  while (1) {
357  (*iters)++;
358 
359  // Apply the matrix M
360  Vxcopy(nx, ny, nz, w1, RAT(u, VAT2(iz, 1, lev)));
361  itmax_s = 1;
362  iters_s = 0;
363  ierror_s = 0;
364  iok_s = 0;
365  iinfo_s = 0;
366  istop_s = 1;
367  Vmvcs(nx, ny, nz,
368  u, iz, w1, w2, w3, w4,
369  &istop_s, &itmax_s, &iters_s, &ierror_s,
370  nlev, ilev, nlev_real, mgsolv,
371  &iok_s, &iinfo_s,
372  epsiln, errtol, omega, nu1, nu2, mgsmoo,
373  ipc, rpc,
374  pc, ac, cc, fc, tru);
375  Vxcopy(nx, ny, nz, RAT(u, VAT2(iz, 1, lev)), w1);
376 
377  // Normalize the new vector
378  denom = Vxnrm2(nx, ny, nz, w1);
379  fac = 1.0 / denom;
380  Vxscal(nx, ny, nz, &fac, w1);
381 
382  // Compute the new raleigh quotient
383  Vxcopy(nx, ny, nz, w1, RAT(u, VAT2(iz, 1, lev)));
384  itmax_s = 1;
385  iters_s = 0;
386  ierror_s = 0;
387  iok_s = 0;
388  iinfo_s = 0;
389  istop_s = 1;
390  Vmvcs(nx, ny, nz,
391  u, iz, w0, w2, w3, w4,
392  &istop_s, &itmax_s, &iters_s, &ierror_s,
393  nlev, ilev, nlev_real, mgsolv,
394  &iok_s, &iinfo_s,
395  epsiln, errtol, omega, nu1, nu2, mgsmoo,
396  ipc, rpc,
397  pc, ac, cc, fc, tru);
398  Vxcopy(nx, ny, nz, RAT(u, VAT2(iz, 1, lev)), w2);
399  rho = Vxdot(nx, ny, nz, w1, w2);
400 
401  // Stopping test
402  // w2=A*x, w1=x, stop = 2-norm(A*x-lamda*x)
403  alpha = -1.0;
404  Vxcopy(nx, ny, nz, w1, w3);
405  Vxcopy(nx, ny, nz, w2, w4);
406  Vxscal(nx, ny, nz, &rho, w3);
407  Vxaxpy(nx, ny, nz, &alpha, w3, w4);
408  error = Vxnrm2(nx, ny, nz, w4);
409  relerr = VABS( rho - oldrho ) / VABS( rho );
410 
411  // I/O
412  if (*iinfo > 3) {
413  Vnm_print(2, "Vmpower: iter=%d; error=%f; relerr=%f; estimate=%f",
414  *iters, error, relerr, rho);
415  }
416 
417  if ((relerr < *tol) || (*iters == *itmax)) {
418  break;
419  }
420 
421  oldrho = rho;
422  }
423  }
424 
425  *eigmax = rho;
426 }
int mgsolv
Definition: vpmgp.h:174
int iinfo
Definition: vpmgp.h:130
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:52
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:313
VPUBLIC void Vazeros(int *nx, int *ny, int *nz, double *x)
Zero out operation for a grid function, including boundary values.
Definition: mikpckd.c:190
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:52
int mgsmoo
Definition: vpmgp.h:160
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:52
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:168
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:160
int nu1
Definition: vpmgp.h:158
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:52
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:107
VPUBLIC void Vaxrand(int *nx, int *ny, int *nz, double *x)
Fill grid function with random values, including boundary values.
Definition: mikpckd.c:286
int ny
Definition: vpmgp.h:84
int nx
Definition: vpmgp.h:83
VPUBLIC double Vxnrm2(int *nx, int *ny, int *nz, double *x)
Norm operation for a grid function with boundary values.
Definition: mikpckd.c:147
int nu2
Definition: vpmgp.h:159
int itmax
Definition: vpmgp.h:122
int nlev
Definition: vpmgp.h:86
double errtol
Definition: vpmgp.h:121
int nz
Definition: vpmgp.h:85