source: branches/profile/GDE/PHYML/eigen.c

Last change on this file was 4073, checked in by westram, 13 years ago
  • phyml 2.4.5
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 22.7 KB
Line 
1/***********************************************************
2*  This eigen() routine works for eigenvalue/vector analysis
3*         for real general square matrix A
4*         A will be destroyed
5*         rr,ri are vectors containing eigenvalues
6*         vr,vi are matrices containing (right) eigenvectors
7*
8*              A*[vr+vi*i] = [vr+vi*i] * diag{rr+ri*i}
9*
10*  Algorithm: Handbook for Automatic Computation, vol 2
11*             by Wilkinson and Reinsch, 1971
12*             most of source codes were taken from a public domain
13*             solftware called MATCALC.
14*  Credits:   to the authors of MATCALC
15*
16*  return     -1 not converged
17*              0 no complex eigenvalues/vectors
18*              1 complex eigenvalues/vectors
19*  Tianlin Wang at University of Illinois
20*  Thu May  6 15:22:31 CDT 1993
21***************************************************************/
22
23#include "utilities.h"
24#include "eigen.h"
25
26#define BASE        2    /* base of floating point arithmetic */
27#define DIGITS     40    /* no. of digits to the base BASE in the fraction */
28/*
29#define DIGITS     53
30*/
31#define MAXITER    30    /* max2. no. of iterations to converge */
32
33#define pos(i,j,n)      ((i)*(n)+(j))
34
35
36/* rr/vr : real parts of eigen values/vectors */
37/* ri/vi : imaginary part s of eigen values/vectors */
38
39
40int eigen(int job, double *A, int n, double *rr, double *ri, 
41          double *vr, double *vi, double *work)
42{   
43/* job=0: eigen values only
44       1: both eigen values and eigen vectors
45   double w[n*2]: work space
46*/
47    int low,hi,i,j,k, it, istate=0;
48    double tiny=sqrt(pow((double)BASE,(double)(1-DIGITS))), t; 
49
50/*     printf("EIGEN\n"); */
51
52    balance(A,n,&low,&hi,work);
53    elemhess(job,A,n,low,hi,vr,vi, (int*)(work+n));
54    if (-1 == realeig(job,A,n,low,hi,rr,ri,vr,vi)) return (-1);
55    if (job) unbalance(n,vr,vi,low,hi,work);
56
57/* sort, added by Z. Yang */
58   for (i=0; i<n; i++) {
59       for (j=i+1,it=i,t=rr[i]; j<n; j++)
60           if (t<rr[j]) { t=rr[j]; it=j; }
61       rr[it]=rr[i];   rr[i]=t;
62       t=ri[it];       ri[it]=ri[i];  ri[i]=t;
63       for (k=0; k<n; k++) {
64          t=vr[k*n+it];  vr[k*n+it]=vr[k*n+i];  vr[k*n+i]=t;
65          t=vi[k*n+it];  vi[k*n+it]=vi[k*n+i];  vi[k*n+i]=t;
66       }
67       if (fabs(ri[i])>tiny) istate=1;
68   }
69
70    return (istate) ;
71}
72
73/* complex funcctions
74*/
75
76complex compl (double re,double im)
77{
78    complex r;
79
80    r.re = re;
81    r.im = im;
82    return(r);
83}
84
85complex _conj (complex a)
86{
87    a.im = -a.im;
88    return(a);
89}
90
91#define csize(a) (fabs(a.re)+fabs(a.im))
92
93complex cplus (complex a, complex b)
94{
95   complex c;
96   c.re = a.re+b.re; 
97   c.im = a.im+b.im;   
98   return (c);
99}
100
101complex cminus (complex a, complex b)
102{
103   complex c;
104   c.re = a.re-b.re; 
105   c.im = a.im-b.im;   
106   return (c);
107}
108
109complex cby (complex a, complex b)
110{
111   complex c;
112   c.re = a.re*b.re-a.im*b.im ;
113   c.im = a.re*b.im+a.im*b.re ;
114   return (c);
115}
116
117complex cdiv (complex a,complex b)
118{
119    double ratio, den;
120    complex c;
121
122    if (fabs(b.re) <= fabs(b.im)) {
123        ratio = b.re / b.im;
124        den = b.im * (1 + ratio * ratio);
125        c.re = (a.re * ratio + a.im) / den;
126        c.im = (a.im * ratio - a.re) / den;
127    }
128    else {
129        ratio = b.im / b.re;
130        den = b.re * (1 + ratio * ratio);
131        c.re = (a.re + a.im * ratio) / den;
132        c.im = (a.im - a.re * ratio) / den;
133    }
134    return(c);
135}
136
137/* complex local_cexp (complex a) */
138/* { */
139/*    complex c; */
140/*    c.re = exp(a.re); */
141/*    if (fabs(a.im)==0) c.im = 0;  */
142/*    else  { c.im = c.re*sin(a.im); c.re*=cos(a.im); } */
143/*    return (c); */
144/* } */
145
146complex cfactor (complex x, double a)
147{
148   complex c;
149   c.re = a*x.re; 
150   c.im = a*x.im;
151   return (c);
152}
153
154int cxtoy (complex *x, complex *y, int n)
155{
156   int i;
157   For (i,n) y[i]=x[i];
158   return (0);
159}
160
161int cmatby (complex *a, complex *b, complex *c, int n,int m,int k)
162/* a[n*m], b[m*k], c[n*k]  ......  c = a*b
163*/
164{
165   int i,j,i1;
166   complex t;
167
168   For (i,n)  For(j,k) {
169       for (i1=0,t=compl(0,0); i1<m; i1++) 
170           t = cplus (t, cby(a[i*m+i1],b[i1*k+j]));
171       c[i*k+j] = t;
172   }
173   return (0);
174}
175
176int cmatinv( complex *x, int n, int m, double *space)
177{
178/* x[n*m]  ... m>=n
179*/
180   int i,j,k, *irow=(int*) space;
181   double xmaxsize, ee=1e-20;
182   complex xmax, t,t1;
183
184   For(i,n)  {
185       xmaxsize = 0.;
186       for (j=i; j<n; j++) {
187          if ( xmaxsize < csize (x[j*m+i]))  {
188               xmaxsize = csize (x[j*m+i]);
189               xmax = x[j*m+i];
190               irow[i] = j;
191          }
192       }
193       if (xmaxsize < ee)   {
194           printf("\nDet goes to zero at %8d!\t\n", i+1);
195           return(-1);
196       }
197       if (irow[i] != i) {
198           For(j,m) {
199                t = x[i*m+j];
200                x[i*m+j] = x[irow[i]*m+j];
201                x[ irow[i]*m+j] = t;
202           }
203       }
204       t = cdiv (compl(1,0), x[i*m+i]);
205       For(j,n) {
206           if (j == i) continue;
207           t1 = cby (t,x[j*m+i]);
208           For(k,m)  x[j*m+k] = cminus (x[j*m+k], cby(t1,x[i*m+k]));
209           x[j*m+i] = cfactor (t1, -1);
210       }
211       For(j,m)   x[i*m+j] = cby (x[i*m+j], t);
212       x[i*m+i] = t;
213   }                         
214   for (i=n-1; i>=0; i--) {
215        if (irow[i] == i) continue;
216        For(j,n)  {
217           t = x[j*m+i];
218           x[j*m+i] = x[j*m+irow[i]];
219           x[ j*m+irow[i]] = t;
220        }
221   }
222   return (0);
223}
224
225
226void balance(double *mat, int n,int *low, int *hi, double *scale)
227{
228/* Balance a matrix for calculation of eigenvalues and eigenvectors
229*/
230    double c,f,g,r,s;
231    int i,j,k,l,done;
232        /* search for rows isolating an eigenvalue and push them down */
233    for (k = n - 1; k >= 0; k--) {
234        for (j = k; j >= 0; j--) {
235            for (i = 0; i <= k; i++) {
236                if (i != j && fabs(mat[pos(j,i,n)]) != 0) break;
237            }
238
239            if (i > k) {
240                scale[k] = j;
241
242                if (j != k) {
243                    for (i = 0; i <= k; i++) {
244                       c = mat[pos(i,j,n)];
245                       mat[pos(i,j,n)] = mat[pos(i,k,n)];
246                       mat[pos(i,k,n)] = c;
247                    }
248
249                    for (i = 0; i < n; i++) {
250                       c = mat[pos(j,i,n)];
251                       mat[pos(j,i,n)] = mat[pos(k,i,n)];
252                       mat[pos(k,i,n)] = c;
253                    }
254                }
255                break;
256            }
257        }
258        if (j < 0) break;
259    }
260
261    /* search for columns isolating an eigenvalue and push them left */
262
263    for (l = 0; l <= k; l++) {
264        for (j = l; j <= k; j++) {
265            for (i = l; i <= k; i++) {
266                if (i != j && fabs(mat[pos(i,j,n)]) != 0) break;
267            }
268            if (i > k) {
269                scale[l] = j;
270                if (j != l) {
271                    for (i = 0; i <= k; i++) {
272                       c = mat[pos(i,j,n)];
273                       mat[pos(i,j,n)] = mat[pos(i,l,n)];
274                       mat[pos(i,l,n)] = c;
275                    }
276
277                    for (i = l; i < n; i++) {
278                       c = mat[pos(j,i,n)];
279                       mat[pos(j,i,n)] = mat[pos(l,i,n)];
280                       mat[pos(l,i,n)] = c;
281                    }
282                }
283
284                break;
285            }
286        }
287
288        if (j > k) break;
289    }
290
291    *hi = k;
292    *low = l;
293
294    /* balance the submatrix in rows l through k */
295
296    for (i = l; i <= k; i++) {
297        scale[i] = 1;
298    }
299
300    do {
301        for (done = 1,i = l; i <= k; i++) {
302            for (c = 0,r = 0,j = l; j <= k; j++) {
303                if (j != i) {
304                    c += fabs(mat[pos(j,i,n)]);
305                    r += fabs(mat[pos(i,j,n)]);
306                }
307            }
308
309            if (c != 0 && r != 0) {
310                g = r / BASE;
311                f = 1;
312                s = c + r;
313
314                while (c < g) {
315                    f *= BASE;
316                    c *= BASE * BASE;
317                }
318
319                g = r * BASE;
320
321                while (c >= g) {
322                    f /= BASE;
323                    c /= BASE * BASE;
324                }
325
326                if ((c + r) / f < 0.95 * s) {
327                    done = 0;
328                    g = 1 / f;
329                    scale[i] *= f;
330
331                    for (j = l; j < n; j++) {
332                        mat[pos(i,j,n)] *= g;
333                    }
334
335                    for (j = 0; j <= k; j++) {
336                        mat[pos(j,i,n)] *= f;
337                    }
338                }
339            }
340        }
341    } while (!done);
342}
343
344
345/*
346 * Transform back eigenvectors of a balanced matrix
347 * into the eigenvectors of the original matrix
348 */
349void unbalance(int n,double *vr,double *vi, int low, int hi, double *scale)
350{
351    int i,j,k;
352    double tmp;
353
354    for (i = low; i <= hi; i++) {
355        for (j = 0; j < n; j++) {
356            vr[pos(i,j,n)] *= scale[i];
357            vi[pos(i,j,n)] *= scale[i];
358        }
359    }
360
361    for (i = low - 1; i >= 0; i--) {
362        if ((k = (int)scale[i]) != i) {
363            for (j = 0; j < n; j++) {
364                tmp = vr[pos(i,j,n)];
365                vr[pos(i,j,n)] = vr[pos(k,j,n)];
366                vr[pos(k,j,n)] = tmp;
367
368                tmp = vi[pos(i,j,n)];
369                vi[pos(i,j,n)] = vi[pos(k,j,n)];
370                vi[pos(k,j,n)] = tmp;       
371            }
372        }
373    }
374
375    for (i = hi + 1; i < n; i++) {
376        if ((k = (int)scale[i]) != i) {
377            for (j = 0; j < n; j++) {
378                tmp = vr[pos(i,j,n)];
379                vr[pos(i,j,n)] = vr[pos(k,j,n)];
380                vr[pos(k,j,n)] = tmp;
381
382                tmp = vi[pos(i,j,n)];
383                vi[pos(i,j,n)] = vi[pos(k,j,n)];
384                vi[pos(k,j,n)] = tmp;       
385            }
386        }
387    }
388}
389
390/*
391 * Reduce the submatrix in rows and columns low through hi of real matrix mat to
392 * Hessenberg form by elementary similarity transformations
393 */
394void elemhess(int job,double *mat,int n,int low,int hi, double *vr,
395              double *vi, int *work)
396{
397/* work[n] */
398    int i,j,m;
399    double x,y;
400
401    for (m = low + 1; m < hi; m++) {
402        for (x = 0,i = m,j = m; j <= hi; j++) {
403            if (fabs(mat[pos(j,m-1,n)]) > fabs(x)) {
404                x = mat[pos(j,m-1,n)];
405                i = j;
406            }
407        }
408
409        if ((work[m] = i) != m) {
410            for (j = m - 1; j < n; j++) {
411               y = mat[pos(i,j,n)];
412               mat[pos(i,j,n)] = mat[pos(m,j,n)];
413               mat[pos(m,j,n)] = y;
414            }
415
416            for (j = 0; j <= hi; j++) {
417               y = mat[pos(j,i,n)];
418               mat[pos(j,i,n)] = mat[pos(j,m,n)];
419               mat[pos(j,m,n)] = y;
420            }
421        }
422
423        if (x != 0) {
424            for (i = m + 1; i <= hi; i++) {
425                if ((y = mat[pos(i,m-1,n)]) != 0) {
426                    y = mat[pos(i,m-1,n)] = y / x;
427
428                    for (j = m; j < n; j++) {
429                        mat[pos(i,j,n)] -= y * mat[pos(m,j,n)];
430                    }
431
432                    for (j = 0; j <= hi; j++) {
433                        mat[pos(j,m,n)] += y * mat[pos(j,i,n)];
434                    }
435                }
436            }
437        }
438    }
439    if (job) {
440       for (i=0; i<n; i++) {
441          for (j=0; j<n; j++) {
442             vr[pos(i,j,n)] = 0.0; vi[pos(i,j,n)] = 0.0;
443          }
444          vr[pos(i,i,n)] = 1.0;
445       }
446
447       for (m = hi - 1; m > low; m--) {
448          for (i = m + 1; i <= hi; i++) {
449             vr[pos(i,m,n)] = mat[pos(i,m-1,n)];
450          }
451
452         if ((i = work[m]) != m) {
453            for (j = m; j <= hi; j++) {
454               vr[pos(m,j,n)] = vr[pos(i,j,n)];
455               vr[pos(i,j,n)] = 0.0;
456            }
457            vr[pos(i,m,n)] = 1.0;
458         }
459      }
460   }
461}
462
463/*
464 * Calculate eigenvalues and eigenvectors of a real upper Hessenberg matrix
465 * Return 1 if converges successfully and 0 otherwise
466 */
467 
468int realeig(int job,double *mat,int n,int low, int hi, double *valr,
469      double *vali, double *vr,double *vi)
470{
471   complex v;
472   double p=0,q=0,r=0,s=0,t,w,x,y,z=0,ra,sa,norm,eps;
473   int niter,en,i,j,k,l,m;
474   double precision  = pow((double)BASE,(double)(1-DIGITS));
475
476   eps = precision;
477   for (i=0; i<n; i++) {
478      valr[i]=0.0;
479      vali[i]=0.0;
480   }
481      /* store isolated roots and calculate norm */
482   for (norm = 0,i = 0; i < n; i++) {
483      for (j = MAX(0,i-1); j < n; j++) {
484         norm += fabs(mat[pos(i,j,n)]);
485      }
486      if (i < low || i > hi) valr[i] = mat[pos(i,i,n)];
487   }
488   t = 0;
489   en = hi;
490
491   while (en >= low) {
492      niter = 0;
493      for (;;) {
494
495       /* look for single small subdiagonal element */
496
497         for (l = en; l > low; l--) {
498            s = fabs(mat[pos(l-1,l-1,n)]) + fabs(mat[pos(l,l,n)]);
499            if (s == 0) s = norm;
500            if (fabs(mat[pos(l,l-1,n)]) <= eps * s) break;
501         }
502
503         /* form shift */
504
505         x = mat[pos(en,en,n)];
506
507         if (l == en) {             /* one root found */
508            valr[en] = x + t;
509            if (job) mat[pos(en,en,n)] = x + t;
510            en--;
511            break;
512         }
513
514         y = mat[pos(en-1,en-1,n)];
515         w = mat[pos(en,en-1,n)] * mat[pos(en-1,en,n)];
516
517         if (l == en - 1) {                /* two roots found */
518            p = (y - x) / 2;
519            q = p * p + w;
520            z = sqrt(fabs(q));
521            x += t;
522            if (job) {
523               mat[pos(en,en,n)] = x;
524               mat[pos(en-1,en-1,n)] = y + t;
525            }
526            if (q < 0) {                /* complex pair */
527               valr[en-1] = x+p;
528               vali[en-1] = z;
529               valr[en] = x+p;
530               vali[en] = -z;
531            }
532            else {                      /* real pair */
533               z = (p < 0) ? p - z : p + z;
534               valr[en-1] = x + z;
535               valr[en] = (z == 0) ? x + z : x - w / z;
536               if (job) {
537                  x = mat[pos(en,en-1,n)];
538                  s = fabs(x) + fabs(z);
539                  p = x / s;
540                  q = z / s;
541                  r = sqrt(p*p+q*q);
542                  p /= r;
543                  q /= r;
544                  for (j = en - 1; j < n; j++) {
545                     z = mat[pos(en-1,j,n)];
546                     mat[pos(en-1,j,n)] = q * z + p *
547                     mat[pos(en,j,n)];
548                     mat[pos(en,j,n)] = q * mat[pos(en,j,n)] - p*z;
549                  }
550                  for (i = 0; i <= en; i++) {
551                     z = mat[pos(i,en-1,n)];
552                     mat[pos(i,en-1,n)] = q * z + p * mat[pos(i,en,n)];
553                     mat[pos(i,en,n)] = q * mat[pos(i,en,n)] - p*z;
554                  }
555                  for (i = low; i <= hi; i++) {
556                     z = vr[pos(i,en-1,n)];
557                     vr[pos(i,en-1,n)] = q*z + p*vr[pos(i,en,n)];
558                     vr[pos(i,en,n)] = q*vr[pos(i,en,n)] - p*z;
559                  }
560               }
561            }
562            en -= 2;
563            break;
564         }
565         if (niter == MAXITER) return(-1);
566         if (niter != 0 && niter % 10 == 0) {
567            t += x;
568            for (i = low; i <= en; i++) mat[pos(i,i,n)] -= x;
569            s = fabs(mat[pos(en,en-1,n)]) + fabs(mat[pos(en-1,en-2,n)]);
570            x = y = 0.75 * s;
571            w = -0.4375 * s * s;
572         }
573         niter++;
574           /* look for two consecutive small subdiagonal elements */
575         for (m = en - 2; m >= l; m--) {
576            z = mat[pos(m,m,n)];
577            r = x - z;
578            s = y - z;
579            p = (r * s - w) / mat[pos(m+1,m,n)] + mat[pos(m,m+1,n)];
580            q = mat[pos(m+1,m+1,n)] - z - r - s;
581            r = mat[pos(m+2,m+1,n)];
582            s = fabs(p) + fabs(q) + fabs(r);
583            p /= s;
584            q /= s;
585            r /= s;
586            if (m == l || fabs(mat[pos(m,m-1,n)]) * (fabs(q)+fabs(r)) <=
587                eps * (fabs(mat[pos(m-1,m-1,n)]) + fabs(z) +
588                fabs(mat[pos(m+1,m+1,n)])) * fabs(p)) break;
589         }
590         for (i = m + 2; i <= en; i++) mat[pos(i,i-2,n)] = 0;
591         for (i = m + 3; i <= en; i++) mat[pos(i,i-3,n)] = 0;
592             /* double QR step involving rows l to en and columns m to en */
593         for (k = m; k < en; k++) {
594            if (k != m) {
595               p = mat[pos(k,k-1,n)];
596               q = mat[pos(k+1,k-1,n)];
597               r = (k == en - 1) ? 0 : mat[pos(k+2,k-1,n)];
598               if ((x = fabs(p) + fabs(q) + fabs(r)) == 0) continue;
599               p /= x;
600               q /= x;
601               r /= x;
602            }
603            s = sqrt(p*p+q*q+r*r);
604            if (p < 0) s = -s;
605            if (k != m) {
606               mat[pos(k,k-1,n)] = -s * x;
607            }
608            else if (l != m) {
609               mat[pos(k,k-1,n)] = -mat[pos(k,k-1,n)];
610            }
611            p += s;
612            x = p / s;
613            y = q / s;
614            z = r / s;
615            q /= p;
616            r /= p;
617                /* row modification */
618            for (j = k; j <= (!job ? en : n-1); j++){
619               p = mat[pos(k,j,n)] + q * mat[pos(k+1,j,n)];
620               if (k != en - 1) {
621                  p += r * mat[pos(k+2,j,n)];
622                  mat[pos(k+2,j,n)] -= p * z;
623               }
624               mat[pos(k+1,j,n)] -= p * y;
625               mat[pos(k,j,n)] -= p * x;
626            }
627            j = MIN(en,k+3);
628              /* column modification */
629            for (i = (!job ? l : 0); i <= j; i++) {
630               p = x * mat[pos(i,k,n)] + y * mat[pos(i,k+1,n)];
631               if (k != en - 1) {
632                  p += z * mat[pos(i,k+2,n)];
633                  mat[pos(i,k+2,n)] -= p*r;
634               }
635               mat[pos(i,k+1,n)] -= p*q;
636               mat[pos(i,k,n)] -= p;
637            }
638            if (job) {             /* accumulate transformations */
639               for (i = low; i <= hi; i++) {
640                  p = x * vr[pos(i,k,n)] + y * vr[pos(i,k+1,n)];
641                  if (k != en - 1) {
642                     p += z * vr[pos(i,k+2,n)];
643                     vr[pos(i,k+2,n)] -= p*r;
644                  }
645                  vr[pos(i,k+1,n)] -= p*q;
646                  vr[pos(i,k,n)] -= p;
647               }
648            }
649         }
650      }
651   }
652
653   if (!job) return(0);
654   if (norm != 0) {
655       /* back substitute to find vectors of upper triangular form */
656      for (en = n-1; en >= 0; en--) {
657         p = valr[en];
658         if ((q = vali[en]) < 0) {            /* complex vector */
659            m = en - 1;
660            if (fabs(mat[pos(en,en-1,n)]) > fabs(mat[pos(en-1,en,n)])) {
661               mat[pos(en-1,en-1,n)] = q / mat[pos(en,en-1,n)];
662               mat[pos(en-1,en,n)] = (p - mat[pos(en,en,n)]) /
663                     mat[pos(en,en-1,n)];
664            }
665            else {
666               v = cdiv(compl(0.0,-mat[pos(en-1,en,n)]),
667                    compl(mat[pos(en-1,en-1,n)]-p,q));
668               mat[pos(en-1,en-1,n)] = v.re;
669               mat[pos(en-1,en,n)] = v.im;
670            }
671            mat[pos(en,en-1,n)] = 0;
672            mat[pos(en,en,n)] = 1;
673            for (i = en - 2; i >= 0; i--) {
674               w = mat[pos(i,i,n)] - p;
675               ra = 0;
676               sa = mat[pos(i,en,n)];
677               for (j = m; j < en; j++) {
678                  ra += mat[pos(i,j,n)] * mat[pos(j,en-1,n)];
679                  sa += mat[pos(i,j,n)] * mat[pos(j,en,n)];
680               }
681               if (vali[i] < 0) {
682                  z = w;
683                  r = ra;
684                  s = sa;
685               }
686               else {
687                  m = i;
688                  if (vali[i] == 0) {
689                     v = cdiv(compl(-ra,-sa),compl(w,q));
690                     mat[pos(i,en-1,n)] = v.re;
691                     mat[pos(i,en,n)] = v.im;
692                  }
693                  else {                      /* solve complex equations */
694                     x = mat[pos(i,i+1,n)];
695                     y = mat[pos(i+1,i,n)];
696                     v.re = (valr[i]- p)*(valr[i]-p) + vali[i]*vali[i] - q*q;
697                     v.im = (valr[i] - p)*2*q;
698                     if ((fabs(v.re) + fabs(v.im)) == 0) {
699                        v.re = eps * norm * (fabs(w) +
700                                fabs(q) + fabs(x) + fabs(y) + fabs(z));
701                     }
702                     v = cdiv(compl(x*r-z*ra+q*sa,x*s-z*sa-q*ra),v);
703                     mat[pos(i,en-1,n)] = v.re;
704                     mat[pos(i,en,n)] = v.im;
705                     if (fabs(x) > fabs(z) + fabs(q)) {
706                        mat[pos(i+1,en-1,n)] = 
707                             (-ra - w * mat[pos(i,en-1,n)] +
708                             q * mat[pos(i,en,n)]) / x;
709                        mat[pos(i+1,en,n)] = (-sa - w * mat[pos(i,en,n)] -
710                             q * mat[pos(i,en-1,n)]) / x;
711                     }
712                     else {
713                        v = cdiv(compl(-r-y*mat[pos(i,en-1,n)],
714                             -s-y*mat[pos(i,en,n)]),compl(z,q));
715                        mat[pos(i+1,en-1,n)] = v.re;
716                        mat[pos(i+1,en,n)] = v.im;
717                     }
718                  }
719               }
720            }
721         }
722         else if (q == 0) {                             /* real vector */
723            m = en;
724            mat[pos(en,en,n)] = 1;
725            for (i = en - 1; i >= 0; i--) {
726               w = mat[pos(i,i,n)] - p;
727               r = mat[pos(i,en,n)];
728               for (j = m; j < en; j++) {
729                  r += mat[pos(i,j,n)] * mat[pos(j,en,n)];
730               }
731               if (vali[i] < 0) {
732                  z = w;
733                  s = r;
734               }
735               else {
736                  m = i;
737                  if (vali[i] == 0) {
738                     if ((t = w) == 0) t = eps * norm;
739                     mat[pos(i,en,n)] = -r / t;
740                  }
741                  else {            /* solve real equations */
742                     x = mat[pos(i,i+1,n)];
743                     y = mat[pos(i+1,i,n)];
744                     q = (valr[i] - p) * (valr[i] - p) + vali[i]*vali[i];
745                     t = (x * s - z * r) / q;
746                     mat[pos(i,en,n)] = t;
747                     if (fabs(x) <= fabs(z)) {
748                        mat[pos(i+1,en,n)] = (-s - y * t) / z;
749                     }
750                     else {
751                        mat[pos(i+1,en,n)] = (-r - w * t) / x;
752                     }
753                  }
754               }
755            }
756         }
757      }
758             /* vectors of isolated roots */
759      for (i = 0; i < n; i++) {
760         if (i < low || i > hi) {
761            for (j = i; j < n; j++) {
762               vr[pos(i,j,n)] = mat[pos(i,j,n)];
763            }
764         }
765      }
766       /* multiply by transformation matrix */
767
768      for (j = n-1; j >= low; j--) {
769         m = MIN(j,hi);
770         for (i = low; i <= hi; i++) {
771            for (z = 0,k = low; k <= m; k++) {
772               z += vr[pos(i,k,n)] * mat[pos(k,j,n)];
773            }
774            vr[pos(i,j,n)] = z;
775         }
776      }
777   }
778    /* rearrange complex eigenvectors */
779   for (j = 0; j < n; j++) {
780      if (vali[j] != 0) {
781         for (i = 0; i < n; i++) {
782            vi[pos(i,j,n)] = vr[pos(i,j+1,n)];
783            vr[pos(i,j+1,n)] = vr[pos(i,j,n)];
784            vi[pos(i,j+1,n)] = -vi[pos(i,j,n)];
785         }
786         j++;
787      }
788   }
789   return(0);
790}
Note: See TracBrowser for help on using the repository browser.