APBS  3.0.0
powerd.c
1 
55 #include "powerd.h"
56 
57 VPUBLIC 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 
165 VPUBLIC 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 
294 VEXTERNC 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