1 /******************************************************************
2
3 zero_cfrac.c generates zero points and associated residues
4 of a continued fraction expansion terminated at 2n+2 level
5 of the Ferm-Dirac function, which is derived from
6 a hypergeometric function.
7
8 This code is distributed under the constitution of GNU-GPL.
9
10 (C) Taisuke Ozaki (AIST-RICS)
11
12 Log of zero_cfrac.c:
13
14 14/July/2005 Released by T.Ozaki
15
16 ******************************************************************/
17
18 #include <stdio.h>
19 #include <stdlib.h>
20 #include <math.h>
21 #include <search.h>
22 #include <string.h>
23
24 #include "openmx_common.h"
25 #include "lapack_prototypes.h"
26 #include "mpi.h"
27
28
29 static void Eigen_DGGEVX( int n, double **a, double **s, double *eval,
30 double *resr, double *resi );
31
zero_cfrac(int N,dcomplex * zp,dcomplex * Rp)32 void zero_cfrac( int N, dcomplex *zp, dcomplex *Rp )
33 {
34
35 int i,j,M;
36 double **A,**B,*zp0,*Rp0;
37
38 /* find the number of zeros */
39
40 M = 2*N;
41
42 /* allocation of arrays */
43
44 A = (double**)malloc(sizeof(double*)*(M+2));
45 for (i=0; i<(M+2); i++){
46 A[i] = (double*)malloc(sizeof(double)*(M+2));
47 }
48
49 B = (double**)malloc(sizeof(double*)*(M+2));
50 for (i=0; i<(M+2); i++){
51 B[i] = (double*)malloc(sizeof(double)*(M+2));
52 }
53
54 zp0 = (double*)malloc(sizeof(double)*(M+2));
55 Rp0 = (double*)malloc(sizeof(double)*(M+2));
56
57 /* initialize arrays */
58
59 for (i=0; i<(M+2); i++){
60 for (j=0; j<(M+2); j++){
61 A[i][j] = 0.0;
62 B[i][j] = 0.0;
63 }
64 }
65
66 /* set matrix elements */
67
68 for (i=1; i<=M; i++){
69 B[i][i] = (2.0*(double)i - 1.0);
70 }
71
72 for (i=1; i<=(M-1); i++){
73 A[i][i+1] = -0.5;
74 A[i+1][i] = -0.5;
75 }
76
77 /* diagonalization */
78
79 {
80 int i,j;
81 char jobz = 'V';
82 char uplo ='U';
83 static INTEGER itype=1;
84 static INTEGER n,lda,ldb,lwork,info;
85 double *a,*b;
86 double *work;
87
88 n = M;
89 lda = M;
90 ldb = M;
91 lwork = 3*M;
92
93 a = (double*)malloc(sizeof(double)*n*n);
94 b = (double*)malloc(sizeof(double)*n*n);
95 work = (double*)malloc(sizeof(double)*3*n);
96
97 for (i=0; i<n; i++) {
98 for (j=0; j<n; j++) {
99 a[j*n+i] = A[i+1][j+1];
100 b[j*n+i] = B[i+1][j+1];
101 }
102 }
103
104 F77_NAME(dsygv,DSYGV)(&itype, &jobz, &uplo, &n, a, &lda, b, &ldb, zp0, work, &lwork, &info);
105
106 /*
107 printf("info=%2d\n",info);
108 */
109
110 /* shift zp0 by 1 */
111 for (i=n; i>=1; i--){
112 zp0[i]= zp0[i-1];
113 }
114
115 /* store residue */
116
117 for (i=1; i<=n; i++){
118 zp0[i] = -1.0/zp0[i];
119 }
120
121 for (i=0; i<n; i++) {
122 Rp0[i+1] = -a[i*n]*a[i*n]*zp0[i+1]*zp0[i+1]*0.250;
123 }
124
125 free(a);
126 free(b);
127 free(work);
128 }
129
130 for (i=1; i<=N; i++){
131 zp[i-1].r = 0.0;
132 zp[i-1].i = zp0[i];
133 Rp[i-1].r = Rp0[i];
134 Rp[i-1].i = 0.0;
135 }
136
137 /* print result */
138
139 /*
140 for (i=0; i<N; i++){
141 printf("i=%5d zp=%18.14e Rp=%18.14e\n",i,zp[i].i,Rp[i].r);
142 }
143
144 MPI_Finalize();
145 exit(0);
146 */
147
148 /* free of arrays */
149
150 for (i=0; i<(M+2); i++){
151 free(A[i]);
152 }
153 free(A);
154
155 for (i=0; i<(M+2); i++){
156 free(B[i]);
157 }
158 free(B);
159
160 free(zp0);
161 free(Rp0);
162
163 }
164
165
166
167
zero_cfrac2(int n,dcomplex * zp,dcomplex * Rp)168 void zero_cfrac2( int n, dcomplex *zp, dcomplex *Rp )
169 {
170 static int i,j,N;
171 static double **a,**s,*eval,*resr,*resi;
172
173 /* check input parameters */
174
175 if (n<=0){
176 printf("\ncould not find the number of zeros\n\n");
177 MPI_Finalize();
178 exit(0);
179 }
180
181 /* the total number of zeros including minus value */
182
183 N = 2*n + 1;
184
185 /* allocation of arrays */
186
187 a = (double**)malloc(sizeof(double*)*(N+2));
188 for (i=0; i<(N+2); i++){
189 a[i] = (double*)malloc(sizeof(double)*(N+2));
190 }
191
192 s = (double**)malloc(sizeof(double*)*(N+2));
193 for (i=0; i<(N+2); i++){
194 s[i] = (double*)malloc(sizeof(double)*(N+2));
195 }
196
197 eval = (double*)malloc(sizeof(double)*(n+3));
198 resr = (double*)malloc(sizeof(double)*(n+3));
199 resi = (double*)malloc(sizeof(double)*(n+3));
200
201 /* initialize arrays */
202
203 for (i=0; i<(N+2); i++){
204 for (j=0; j<(N+2); j++){
205 a[i][j] = 0.0;
206 s[i][j] = 0.0;
207 }
208 }
209
210 /* set matrix elements */
211
212 s[2][1] = 1.0;
213 s[2][2] = -0.5;
214
215 for (i=3; i<=N; i++){
216 s[i][i-1] = -0.5;
217 s[i-1][i] = 0.5;
218 }
219
220 a[1][1] = -2.0;
221 a[1][2] = 1.0;
222 a[2][2] = -1.0;
223
224 for (i=3; i<=N; i++){
225 a[i][i] = -(2.0*(double)i - 3.0);
226 }
227
228 /* diagonalization */
229
230 Eigen_DGGEVX( N, a, s, eval, resr, resi );
231
232 for (i=0; i<n; i++){
233 zp[i].r = 0.0;
234 zp[i].i = eval[i+1];
235 Rp[i].r = resr[i+1];
236 Rp[i].i = 0.0;
237 }
238
239 /* print result */
240
241 for (i=0; i<n; i++){
242 printf("i=%5d zp=%18.14e Rp=%18.14e\n",i,zp[i].i,Rp[i].r);
243 }
244
245 MPI_Finalize();
246 exit(0);
247
248
249 /* free of arrays */
250
251 for (i=0; i<(N+2); i++){
252 free(a[i]);
253 }
254 free(a);
255
256 for (i=0; i<(N+2); i++){
257 free(s[i]);
258 }
259 free(s);
260
261 free(eval);
262 free(resr);
263 free(resi);
264 }
265
266
267
268
269
270
271
272
273
274
275
276
277
Eigen_DGGEVX(int n,double ** a,double ** s,double * eval,double * resr,double * resi)278 void Eigen_DGGEVX( int n, double **a, double **s, double *eval, double *resr, double *resi )
279 {
280 static int i,j,k,l,num;
281
282 static char balanc = 'N';
283 static char jobvl = 'V';
284 static char jobvr = 'V';
285 static char sense = 'B';
286 static double *A;
287 static double *b;
288 static double *alphar;
289 static double *alphai;
290 static double *beta;
291 static double *vl;
292 static double *vr;
293 static double *lscale;
294 static double *rscale;
295 static double abnrm;
296 static double bbnrm;
297 static double *rconde;
298 static double *rcondv;
299 static double *work;
300 static double *tmpvecr,*tmpveci;
301 static INTEGER *iwork;
302 static INTEGER lda,ldb,ldvl,ldvr,ilo,ihi;
303 static INTEGER lwork,info;
304 static logical *bwork;
305 static double sumr,sumi,tmpr,tmpi;
306 static double *sortnum;
307
308 lda = n;
309 ldb = n;
310 ldvl = n;
311 ldvr = n;
312
313 A = (double*)malloc(sizeof(double)*n*n);
314 b = (double*)malloc(sizeof(double)*n*n);
315 alphar = (double*)malloc(sizeof(double)*n);
316 alphai = (double*)malloc(sizeof(double)*n);
317 beta = (double*)malloc(sizeof(double)*n);
318
319 vl = (double*)malloc(sizeof(double)*n*n);
320 vr = (double*)malloc(sizeof(double)*n*n);
321
322 lscale = (double*)malloc(sizeof(double)*n);
323 rscale = (double*)malloc(sizeof(double)*n);
324
325 rconde = (double*)malloc(sizeof(double)*n);
326 rcondv = (double*)malloc(sizeof(double)*n);
327
328 lwork = 2*n*n + 12*n + 16;
329 work = (double*)malloc(sizeof(double)*lwork);
330
331 iwork = (INTEGER*)malloc(sizeof(INTEGER)*(n+6));
332 bwork = (logical*)malloc(sizeof(logical)*n);
333
334 tmpvecr = (double*)malloc(sizeof(double)*(n+2));
335 tmpveci = (double*)malloc(sizeof(double)*(n+2));
336
337 sortnum = (double*)malloc(sizeof(double)*(n+2));
338
339 /* convert two dimensional arrays to one-dimensional arrays */
340
341 for (i=0; i<n; i++) {
342 for (j=0; j<n; j++) {
343 A[j*n+i]= a[i+1][j+1];
344 b[j*n+i]= s[i+1][j+1];
345 }
346 }
347
348 /* call dggevx_() */
349
350 F77_NAME(dggevx,DGGEVX)(
351 &balanc, &jobvl, & jobvr, &sense, &n, A, &lda, b, &ldb,
352 alphar, alphai, beta, vl, &ldvl, vr, &ldvr, &ilo, &ihi,
353 lscale, rscale, &abnrm, &bbnrm, rconde, rcondv, work,
354 &lwork, iwork, bwork, &info );
355
356 if (info!=0){
357 printf("Errors in dggevx_() info=%2d\n",info);
358 }
359
360 /*
361 for (i=0; i<n; i++){
362 printf("i=%4d %18.13f %18.13f %18.13f\n",i,alphar[i],alphai[i],beta[i]);
363 }
364 printf("\n");
365 */
366
367 num = 0;
368 for (i=0; i<n; i++){
369
370 if ( 1.0e-13<fabs(beta[i]) && 0.0<alphai[i]/beta[i] ){
371
372 /* normalize the eigenvector */
373
374 for (j=0; j<n; j++) {
375
376 sumr = 0.0;
377 sumi = 0.0;
378
379 for (k=0; k<n; k++) {
380 sumr += s[j+1][k+1]*vr[i*n +k];
381 sumi += s[j+1][k+1]*vr[(i+1)*n+k];
382 }
383
384 tmpvecr[j] = sumr;
385 tmpveci[j] = sumi;
386 }
387
388 sumr = 0.0;
389 sumi = 0.0;
390
391 for (k=0; k<n; k++) {
392 sumr += vl[i*n+k]*tmpvecr[k] + vl[(i+1)*n+k]*tmpveci[k];
393 sumi += vl[i*n+k]*tmpveci[k] - vl[(i+1)*n+k]*tmpvecr[k];
394 }
395
396 /* calculate zero point and residue */
397
398 eval[num+1] = alphai[i]/beta[i];
399 tmpr = vr[i*n]*vl[i*n] + vr[(i+1)*n]*vl[(i+1)*n];
400 tmpi = -vr[i*n]*vl[(i+1)*n] + vr[(i+1)*n]*vl[i*n];
401 resr[num+1] = tmpi/sumi;
402 resi[num+1] = -tmpr/sumi;
403
404 num++;
405 }
406 else{
407 /*
408 printf("i=%4d %18.13f %18.13f %18.13f\n",i+1,alphar[i],alphai[i],beta[i]);
409 */
410 }
411 }
412
413 /* check round-off error */
414
415 for (i=1; i<=num; i++){
416 if (1.0e-8<fabs(resi[i])){
417 printf("Could not calculate zero points and residues due to round-off error\n");
418 MPI_Finalize();
419 exit(0);
420 }
421 }
422
423 /* sorting */
424
425 qsort_double(num,eval,resr);
426
427 /* free arraies */
428
429 free(A);
430 free(b);
431 free(alphar);
432 free(alphai);
433 free(beta);
434
435 free(vl);
436 free(vr);
437
438 free(lscale);
439 free(rscale);
440
441 free(rconde);
442 free(rcondv);
443
444 free(work);
445
446 free(iwork);
447 free(bwork);
448
449 free(tmpvecr);
450 free(tmpveci);
451 free(sortnum);
452 }
453
454