EVSL  1.1.0
EigenValues Slicing Library
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Macros
ratfilter.c
Go to the documentation of this file.
1 #include <stdio.h>
2 #include <stdlib.h>
3 #include <string.h>
4 #include <complex.h>
5 #include <math.h>
6 #include "def.h"
7 #include "blaslapack.h"
8 #include "struct.h"
9 #include "internal_proto.h"
10 
23 void contQuad(int method, int n, complex double* zk) {
24  int i, m, INFO;
25  double *beta, *D, *Z, *WORK;
26  char JOBZ = 'V';
27  complex double tmp2;
28  if (method == 0) {
29  m = n-1;
30  Malloc(D, n, double);
31  Malloc(Z, n*n, double);
32  Malloc(WORK, 2*n-2, double);
33  for (i=0; i<n; i++) {
34  D[i] = 0.0;
35  }
36  Malloc(beta, m, double);
37  for (i=0; i<m; i++) {
38  beta[i] = 0.5/(sqrt(1-pow(2*(i+1),-2)));
39  }
40  DSTEV(&JOBZ, &n, D, beta, Z, &n, WORK, &INFO);
41  for (i=0; i<n; i++) {
42  tmp2 = I*M_PI/2.0*(1.0-D[i]);
43  zk[i] = cexp(tmp2);
44  }
45  free(beta);
46  free(D);
47  free(Z);
48  free(WORK);
49  } else if (method == 1) {
50  for (i=0; i<n; i++) {
51  tmp2 = M_PI*I*(2*(i+1)-1)/(2.0*n);
52  zk[i] = cexp(tmp2);
53  }
54  }
55 }
56 
69 void ratf2p2(int n, int *mulp, complex double *zk, complex double* alp, int m,
70  double *z, double *xx) {
71  complex double y, x, t;
72  int i, j, k, k1, k2;
73  for (k2=0; k2<m; k2++) {
74  k = 0;
75  y = 0.0 + 0.0*I;
76  for (j=0; j<n; j++) {
77  x = 1.0 / (z[k2]-zk[j]);
78  k1 = k + mulp[j];
79  t = 0.0+0.0*I;
80  for (i=k1-1;i>=k;i--) {
81  t = x*(alp[i]+t);
82  }
83  y = y+t;
84  k = k1;
85  }
86  xx[k2] = 2.0*creal(y);
87  }
88 }
89 
90 
91 
95 void pfe2(complex double s1, complex double s2, int k1, int k2,
96  complex double* alp, complex double* bet) {
97  int i;
98  complex double d, xp;
99  if (cabs(s1-s2) < 1.0e-12 * (cabs(s1)+cabs(s2))) {
100  for (i=0; i<k1+k2; i++) {
101  alp[i] = 0.0;
102  }
103  alp[k1+k2-1] = 1;
104  } else if ((k1 == 1) && (k2 == 1)) {
105  d = s1-s2;
106  alp[k1-1] = 1.0 / d;
107  bet[k1-1] = -alp[k1-1];
108  } else {
109  d = 1.0 + 0.0*I;
110  xp = 1.0 / cpow((s1-s2),k2);
111  for (i=0; i<k1; i++) {
112  alp[k1-i-1] = d * xp;
113  xp = xp / (s1-s2);
114  d = -d * (k2+i) / (i+1.0);
115  }
116  d = 1.0 + 0.0*I;
117  xp = 1.0 / cpow((s2-s1),k1);
118  for (i=0; i<k2; i++) {
119  bet[k2-i-1] = d * xp;
120  xp = xp / (s2-s1);
121  d = -d * (k1+i) / (i+1.0);
122  }
123  }
124 }
125 
129 complex double integg2(complex double s1, complex double s2,
130  complex double* alp, int k1, complex double* bet,
131  int k2, double a, double b) {
132  complex double t, t1, t0, scal;
133  int k;
134  t = 0.0 + 0.0*I;
135  t1 = 0.0 + 0.0*I;
136  t0 = 0.0 +0.0*I;
137  for (k=0; k<k1; k++) {
138  scal = alp[k];
139  if (k==0) {
140  t0 = scal*clog((b-s1)/(a-s1));
141  } else {
142  t = t - (scal*1.0/k) * (1.0/cpow((b-s1),k)-1.0/cpow((a-s1),k));
143  }
144  }
145  for (k=0; k<k2; k++) {
146  scal = bet[k];
147  if (k==0) {
148  t1 = scal*clog((b-s2)/(a-s2));
149  } else {
150  t = t - (scal*1.0/k)*(1.0/cpow((b-s2),k)-1.0/cpow((a-s2),k));
151  }
152  }
153  t = t + (t0+t1);
154  return t;
155 }
156 
157 
169 void weights(int n, complex double* zk, int* mulp, double lambda,
170  complex double* omega) {
171  int INFO;
172  int nrhs = 1;
173  int *ipiv;
174  int m;
175  double mu = 10.0;
176  int i, j, ii, jj, ki, kj, n1, n2, nf=0, k1, k2;
177  complex double s1, s2, s3, t;
178  complex double *rhs, *A, *B, *mat, *alp, *bet;
179  double scaling;
180  for (i=0; i<n; i++) {
181  nf += mulp[i];
182  }
183  m = 2*nf;
184  Malloc(ipiv, m, int);
185  Malloc(rhs, m, complex double);
186  Malloc(A, nf*nf, complex double);
187  Malloc(B, nf*nf, complex double);
188  Malloc(mat, 4*nf*nf, complex double);
189  for(i=0; i<nf; i++) {
190  for(j=0; j<nf; j++) {
191  A[i*nf+j] = 0.0 + 0.0*I;
192  B[i*nf+j] = 0.0 + 0.0*I;
193  }
194  }
195  if (fabs(lambda) < 1.0e-12) {
196  lambda = 1.0e-5;
197  }
198  ki = 0;
199  for (ii=0; ii<n; ii++) {
200  s1 = zk[ii];
201  n1 = mulp[ii];
202  kj = 0;
203  s3 = conj(s1);
204  for (i=0; i<n1; i++) {
205  if (i==0) {
206  rhs[ki+i] = lambda*clog((s3-1.0)/(s3+1.0));
207  } else {
208  rhs[ki+i] = -lambda*(1.0/(i))*(1/cpow((1.0-s3),i)-1.0/cpow((-1.0-s3),i));
209  }
210  }
211  for (jj=0; jj<n; jj++) {
212  s2 = zk[jj];
213  n2 = mulp[jj];
214  for (i=0; i<n1; i++) {
215  for (j=0; j<n2; j++) {
216  s3 = conj(s2);
217  if (cabs(s1-s3) < 1.0e-12*(cabs(s1)+cabs(s3))) {
218  Malloc(alp, i+j+2, complex double);
219  Malloc(bet, 1, complex double);
220  k1 = i+1+j+1;
221  k2 = 0;
222  } else {
223  Malloc(alp, i+1, complex double);
224  Malloc(bet, j+1, complex double);
225  k1 = i+1;
226  k2 = j+1;
227  }
228  pfe2(s1, s3, k1, k2, alp, bet);
229  t = integg2(s1, s3, alp, k1, bet, k2, -mu, mu);
230  t += (lambda-1)*integg2(s1, s3, alp, k1, bet, k2, -1.0, 1.0);
231  A[(ki+i)*nf+kj+j] = t;
232  free(bet);
233  free(alp);
234  if (cabs(s1-s2) < 1.0e-12*(cabs(s1)+cabs(s2))) {
235  Malloc(alp, i+j+2, complex double);
236  Malloc(bet, 1, complex double);
237  k1 = i+1+j+1;
238  k2 = 0;
239  } else {
240  Malloc(alp, i+1, complex double);
241  Malloc(bet, j+1, complex double);
242  k1 = i+1;
243  k2 = j+1;
244  }
245  pfe2(s1, s2, k1, k2, alp, bet);
246  t = integg2(s1, s2, alp, k1, bet, k2, -mu, mu);
247  t += (lambda-1)*integg2(s1, s2, alp, k1, bet, k2, -1.0, 1.0);
248  B[(ki+i)*nf+kj+j] = t;
249  free(alp);
250  free(bet);
251  }
252  }
253  kj = kj+n2;
254  }
255  ki = ki+n1;
256  }
257  for (i=nf; i<2*nf; i++) {
258  rhs[i] = conj(rhs[i-nf]);
259  }
260  /*---form mat = [A,B;conj(B),conj(A)]---*/
261  /* Note that mat is stored column-wise for lapack routine */
262  for (i=0; i<nf; i++) {
263  for(j=0; j<nf; j++) {
264  mat[i+j*m] = conj(A[i*nf+j]);
265  }
266  }
267  for (i=0; i<nf; i++) {
268  for (j=nf; j<m; j++) {
269  mat[i+j*m] = conj(B[i*nf+j-nf]);
270  }
271  }
272  for (i=nf; i<m; i++) {
273  for (j=0; j<nf; j++) {
274  mat[i+j*m] = B[(i-nf)*nf+j];
275  }
276  }
277  for (i=nf; i<m; i++) {
278  for (j=nf; j<m; j++) {
279  mat[i+j*m] = A[(i-nf)*nf+j-nf];
280  }
281  }
282  ZGESV(&m, &nrhs, mat, &m, ipiv, rhs, &m, &INFO);
283  for(i=0;i<nf;i++) {
284  omega[i] = rhs[i];
285  }
286 
287  /* Scale coefs to let the filter pass through [-1, 0.5] */
288  double aa = 1.0;
289  ratf2p2(n, mulp, zk, omega, 1, &aa, &scaling);
290  scaling = 0.5 / scaling;
291  for (i=0; i<nf; i++) {
292  omega[i] *= scaling;
293  }
294 
295  free(A);
296  free(B);
297  free(rhs);
298  free(mat);
299  free(ipiv);
300 }
301 
302 
313 int scaleweigthts(int n, double a, double b, complex double *zk, int* mulp,
314  complex double* omegaM) {
315  int i, j, k, nf=0;
316  double c, h;
317  c = 0.5 * (a + b);
318  h = 0.5 * (b - a);
319  for (i=0; i<n; i++) {
320  nf += mulp[i];
321  zk[i] = h*zk[i]+c;
322  }
323  /* Transform the coefs for multiple poles */
324  double tmp;
325  k = -1;
326  for (i=0; i<n; i++) {
327  for (j=0; j<mulp[i]; j++) {
328  k = k+1;
329  omegaM[k] = omegaM[k]*cpow(h,j+1);
330  }
331  }
332  /* Scale ration function to let it pass through [a, 1/2] */
333  ratf2p2(n, mulp, zk, omegaM, 1, &a, &tmp);
334  tmp = 0.5 / tmp;
335  for (i=0; i<nf; i++) {
336  omegaM[i] = omegaM[i] * tmp;
337  }
338  return 0;
339 }
340 
341 
346  // -------------------- this sets default values for ratparams struct.
347  rat->num = 1; // number of the poles
348  rat->pw = 2; // default multplicity of each pole
349  rat->method = 1; // using poles from mid-point rule
350  rat->beta = 0.01; // beta in LS approximation
351  rat->bar = 0.5; // this is fixed for rational filter
352  rat->aa = -1.0; // left endpoint of interval
353  rat->bb = 1.0; // right endpoint of interval
354  //rat->cc = 0.0; // center of interval
355  //rat->dd = 1.0; // width of interval
356 }
357 
375 int find_ratf(double *intv, ratparams *rat) {
376  complex double *omega; // weights of the poles
377  complex double *zk; // location of the poles
378  int *mulp; // multiplicity of the each pole
379  int n = rat->num, i, pow = 0, pw = rat->pw, method = rat->method;
380  double beta = rat->beta;
381  /*-------------------- A few parameters to be set or reset */
382  Malloc(mulp, n, int);
383  Malloc(zk, n, complex double);
384  for (i=0; i<n; i++) { // set the multiplicity of each pole
385  mulp[i] = pw;
386  pow += mulp[i];
387  }
388  rat->zk = zk;
389  rat->mulp = mulp;
390  rat->pow = pow; // total multiplicity of the poles
391  Malloc(omega, pow, complex double);
392  rat->omega = omega;
393  //-------------------- intervals related
394  if (check_intv(intv, stdout) < 0) {
395  return -1;
396  }
397  double aa, bb;
398  aa = max(intv[0], intv[2]); bb = min(intv[1], intv[3]);
399  if (intv[0] < intv[2] || intv[1] > intv[3]) {
400  fprintf(stdout, " warning [%s (%d)]: interval (%e, %e) is adjusted to (%e, %e)\n",
401  __FILE__, __LINE__, intv[0], intv[1], aa, bb);
402  }
403  //double lmin = intv[2], lmax = intv[3];
404  /*-------------------- */
405  rat->aa = aa;
406  rat->bb = bb;
407  /*-------------------- cc, rr: center and half-width of [aa, bb] */
408  //double cc = 0.5 * (aa + bb);
409  //double dd = 0.5 * (bb - aa);
410  //rat->cc = cc;
411  //rat->dd = dd;
412  /*------------ compute the location of the poles */
413  contQuad(method, n, zk);
414  /*------------ compute expansion coefficients of rational filter on [-1, 1] */
415  weights(n, zk, mulp, beta, omega);
416  /*-------------------- compute expansion coefficients on [aa, bb]*/
417  scaleweigthts(n, aa, bb, zk, mulp, omega);
418 
419  rat->ASIGBsol = NULL;
420 
421  return 0;
422 }
423 
424 void free_rat(ratparams *rat) {
425  free(rat->mulp);
426  free(rat->omega);
427  free(rat->zk);
428  free(rat->ASIGBsol);
429 }
430 
448 void RatFiltApply(int n, ratparams *rat, double *b, double *x, double *w6) {
449  double tt = evsl_timer();
450  const int ifGenEv = evsldata.ifGenEv;
451  int jj, kk, k=0, kf;
452  int *mulp = rat->mulp;
453  int num = rat->num;
454  complex double *omega = rat->omega;
455  double dtwo = 2.0;
456  double done = 1.0;
457  int one = 1;
458 
459  double *xr, *xz, *bz, *br, *yr=NULL, *yz=NULL;
460  double zkr, zkc;
461  xr = w6;
462  xz = xr + n;
463  bz = xz + n;
464  br = bz + n;
465  if (ifGenEv) {
466  yr = br + n;
467  yz = yr + n;
468  }
469  /*------------------ loop through each pole */
470  for (kk=0; kk<num; kk++) {
471  /*---------------- solver for A-s[kk]*B */
472  EVSLASIGMABSol *sol = &rat->ASIGBsol[kk];
473  kf = k + mulp[kk];
474  /*------------------ power loop */
475  for (jj=kf-1; jj>=k; jj--) {
476  /*---------------- weight */
477  zkr = creal(omega[jj]);
478  zkc = cimag(omega[jj]);
479  /*---------------- initilize the right hand side */
480  memcpy(br, b, n*sizeof(double));
481  memcpy(bz, b, n*sizeof(double));
482  DSCAL(&n, &zkr, br, &one);
483  DSCAL(&n, &zkc, bz, &one);
484  if (jj != kf-1) {
485  DAXPY(&n, &done, xr, &one, br, &one);
486  DAXPY(&n, &done, xz, &one, bz, &one);
487  }
488  /*---------------- solve shifted system */
489  if (ifGenEv) {
490  if (jj > k) {
491  //(sol->func)(n, br, bz, yr, yz, sol->data);
492  solve_ASigB(sol, n, br, bz, yr, yz);
493  matvec_B(yr, xr);
494  matvec_B(yz, xz);
495  } else {
496  /*------------- jj == k */
497  //(sol->func)(n, br, bz, xr, xz, sol->data);
498  solve_ASigB(sol, n, br, bz, xr, xz);
499  }
500  } else {
501  //(sol->func)(n, br, bz, xr, xz, sol->data);
502  solve_ASigB(sol, n, br, bz, xr, xz);
503  }
504  }
505  /*------------------ solution (real part) */
506  if (kk) {
507  DAXPY(&n, &dtwo, xr, &one, x, &one);
508  } else {
509  memcpy(x, xr, n*sizeof(double));
510  DSCAL(&n, &dtwo, x, &one);
511  }
512  k = kf;
513  }
514 
515  evslstat.n_ratAv ++;
516  evslstat.t_ratAv += evsl_timer() - tt;
517 }
518 
double t_ratAv
Definition: struct.h:212
complex double * zk
Definition: struct.h:126
int num
Definition: struct.h:114
void free_rat(ratparams *rat)
Definition: ratfilter.c:424
double bar
Definition: struct.h:120
void DSTEV(char *jobz, int *n, double *diagonal, double *subdiagonal, double *V, int *ldz, double *work, int *info)
defs in EVSL
int ifGenEv
Definition: struct.h:166
int pow
Definition: struct.h:123
void DSCAL(int *n, double *a, double *x, int *incx)
int pw
Definition: struct.h:115
int method
Definition: struct.h:116
This file contains function prototypes and constant definitions internally used in EVSL...
#define min(a, b)
Definition: def.h:62
EVSLASIGMABSol * ASIGBsol
Definition: struct.h:127
void weights(int n, complex double *zk, int *mulp, double lambda, complex double *omega)
Compute the LS weight for each multiple pole.
Definition: ratfilter.c:169
#define Malloc(base, nmem, type)
Definition: def.h:22
evslStat evslstat
global statistics of EVSL
Definition: evsl.c:21
structs used in evsl
user-provided function and data for solving (A - SIGMA*B) x = b
Definition: struct.h:102
size_t n_ratAv
Definition: struct.h:220
double aa
Definition: struct.h:118
void DAXPY(int *n, double *alpha, double *x, int *incx, double *y, int *incy)
#define max(a, b)
Definition: def.h:56
int scaleweigthts(int n, double a, double b, complex double *zk, int *mulp, complex double *omegaM)
Compute the weights and pole locations on [a, b].
Definition: ratfilter.c:313
Defs for blaslapack routines.
complex double integg2(complex double s1, complex double s2, complex double *alp, int k1, complex double *bet, int k2, double a, double b)
Integration of 1/[(z-s1)^k1 (z-s2)^k2] from a to b.
Definition: ratfilter.c:129
void pfe2(complex double s1, complex double s2, int k1, int k2, complex double *alp, complex double *bet)
Get the fraction expansion of 1/[(z-s1)^k1 (z-s2)^k2].
Definition: ratfilter.c:95
evslData evsldata
global variable of EVSL
Definition: evsl.c:15
double beta
Definition: struct.h:117
void set_ratf_def(ratparams *rat)
Sets default values for ratparams struct.
Definition: ratfilter.c:345
#define M_PI
Definition: def.h:13
int * mulp
Definition: struct.h:122
int find_ratf(double *intv, ratparams *rat)
Definition: ratfilter.c:375
void contQuad(int method, int n, complex double *zk)
Compute the locations of the poles.
Definition: ratfilter.c:23
void ratf2p2(int n, int *mulp, complex double *zk, complex double *alp, int m, double *z, double *xx)
Compute the function value of the multiple pole rational filter at real locations.
Definition: ratfilter.c:69
void RatFiltApply(int n, ratparams *rat, double *b, double *x, double *w6)
Apply rational filter R to a vetor b.
Definition: ratfilter.c:448
complex double * omega
Definition: struct.h:125
double evsl_timer()
evsl timer for mac
Definition: mactime.c:14
double bb
Definition: struct.h:119
parameters for rational filter
Definition: struct.h:112
void ZGESV(int *n, int *nrow, complex double *A, int *m, int *ipiv, complex double *rhs, int *k, int *INFO)