home *** CD-ROM | disk | FTP | other *** search
/ RISC DISC 1 / RISC_DISC_1.iso / pd_share / code / meschach / !Meschach / c / symmeig < prev    next >
Encoding:
Text File  |  1994-02-15  |  5.8 KB  |  212 lines
  1.  
  2. /**************************************************************************
  3. **
  4. ** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  5. **
  6. **                 Meschach Library
  7. ** 
  8. ** This Meschach Library is provided "as is" without any express 
  9. ** or implied warranty of any kind with respect to this software. 
  10. ** In particular the authors shall not be liable for any direct, 
  11. ** indirect, special, incidental or consequential damages arising 
  12. ** in any way from use of the software.
  13. ** 
  14. ** Everyone is granted permission to copy, modify and redistribute this
  15. ** Meschach Library, provided:
  16. **  1.  All copies contain this copyright notice.
  17. **  2.  All modified copies shall carry a notice stating who
  18. **      made the last modification and the date of such modification.
  19. **  3.  No charge is made for this software or works derived from it.  
  20. **      This clause shall not be construed as constraining other software
  21. **      distributed on the same medium as this software, nor is a
  22. **      distribution fee considered a charge.
  23. **
  24. ***************************************************************************/
  25.  
  26.  
  27. /*
  28.     File containing routines for symmetric eigenvalue problems
  29. */
  30.  
  31. #include    <stdio.h>
  32. #include    <math.h>
  33. #include    "matrix.h"
  34. #include        "matrix2.h"
  35.  
  36.  
  37. static char rcsid[] = "$Id: symmeig.c,v 1.5 1994/02/16 03:23:39 des Exp $";
  38.  
  39.  
  40.  
  41. #define    SQRT2    1.4142135623730949
  42. #define    sgn(x)    ( (x) >= 0 ? 1 : -1 )
  43.  
  44. /* trieig -- finds eigenvalues of symmetric tridiagonal matrices
  45.     -- matrix represented by a pair of vectors a (diag entries)
  46.         and b (sub- & super-diag entries)
  47.     -- eigenvalues in a on return */
  48. VEC    *trieig(a,b,Q)
  49. VEC    *a, *b;
  50. MAT    *Q;
  51. {
  52.     int    i, i_min, i_max, n, split;
  53.     Real    *a_ve, *b_ve;
  54.     Real    b_sqr, bk, ak1, bk1, ak2, bk2, z;
  55.     Real    c, c2, cs, s, s2, d, mu;
  56.  
  57.     if ( ! a || ! b )
  58.         error(E_NULL,"trieig");
  59.     if ( a->dim != b->dim + 1 || ( Q && Q->m != a->dim ) )
  60.         error(E_SIZES,"trieig");
  61.     if ( Q && Q->m != Q->n )
  62.         error(E_SQUARE,"trieig");
  63.  
  64.     n = a->dim;
  65.     a_ve = a->ve;        b_ve = b->ve;
  66.  
  67.     i_min = 0;
  68.     while ( i_min < n )        /* outer while loop */
  69.     {
  70.         /* find i_max to suit;
  71.             submatrix i_min..i_max should be irreducible */
  72.         i_max = n-1;
  73.         for ( i = i_min; i < n-1; i++ )
  74.             if ( b_ve[i] == 0.0 )
  75.             {    i_max = i;    break;    }
  76.         if ( i_max <= i_min )
  77.         {
  78.             /* printf("# i_min = %d, i_max = %d\n",i_min,i_max); */
  79.             i_min = i_max + 1;
  80.             continue;    /* outer while loop */
  81.         }
  82.  
  83.         /* printf("# i_min = %d, i_max = %d\n",i_min,i_max); */
  84.  
  85.         /* repeatedly perform QR method until matrix splits */
  86.         split = FALSE;
  87.         while ( ! split )        /* inner while loop */
  88.         {
  89.  
  90.             /* find Wilkinson shift */
  91.             d = (a_ve[i_max-1] - a_ve[i_max])/2;
  92.             b_sqr = b_ve[i_max-1]*b_ve[i_max-1];
  93.             mu = a_ve[i_max] - b_sqr/(d + sgn(d)*sqrt(d*d+b_sqr));
  94.             /* printf("# Wilkinson shift = %g\n",mu); */
  95.  
  96.             /* initial Givens' rotation */
  97.             givens(a_ve[i_min]-mu,b_ve[i_min],&c,&s);
  98.             s = -s;
  99.             /* printf("# c = %g, s = %g\n",c,s); */
  100.             if ( fabs(c) < SQRT2 )
  101.             {    c2 = c*c;    s2 = 1-c2;    }
  102.             else
  103.             {    s2 = s*s;    c2 = 1-s2;    }
  104.             cs = c*s;
  105.             ak1 = c2*a_ve[i_min]+s2*a_ve[i_min+1]-2*cs*b_ve[i_min];
  106.             bk1 = cs*(a_ve[i_min]-a_ve[i_min+1]) +
  107.                         (c2-s2)*b_ve[i_min];
  108.             ak2 = s2*a_ve[i_min]+c2*a_ve[i_min+1]+2*cs*b_ve[i_min];
  109.             bk2 = ( i_min < i_max-1 ) ? c*b_ve[i_min+1] : 0.0;
  110.             z  = ( i_min < i_max-1 ) ? -s*b_ve[i_min+1] : 0.0;
  111.             a_ve[i_min] = ak1;
  112.             a_ve[i_min+1] = ak2;
  113.             b_ve[i_min] = bk1;
  114.             if ( i_min < i_max-1 )
  115.             b_ve[i_min+1] = bk2;
  116.             if ( Q )
  117.             rot_cols(Q,i_min,i_min+1,c,-s,Q);
  118.             /* printf("# z = %g\n",z); */
  119.             /* printf("# a [temp1] =\n");    v_output(a); */
  120.             /* printf("# b [temp1] =\n");    v_output(b); */
  121.  
  122.             for ( i = i_min+1; i < i_max; i++ )
  123.             {
  124.             /* get Givens' rotation for sub-block -- k == i-1 */
  125.             givens(b_ve[i-1],z,&c,&s);
  126.             s = -s;
  127.             /* printf("# c = %g, s = %g\n",c,s); */
  128.  
  129.             /* perform Givens' rotation on sub-block */
  130.                 if ( fabs(c) < SQRT2 )
  131.                 {    c2 = c*c;    s2 = 1-c2;    }
  132.                 else
  133.                 {    s2 = s*s;    c2 = 1-s2;    }
  134.                 cs = c*s;
  135.             bk  = c*b_ve[i-1] - s*z;
  136.             ak1 = c2*a_ve[i]+s2*a_ve[i+1]-2*cs*b_ve[i];
  137.             bk1 = cs*(a_ve[i]-a_ve[i+1]) +
  138.                         (c2-s2)*b_ve[i];
  139.             ak2 = s2*a_ve[i]+c2*a_ve[i+1]+2*cs*b_ve[i];
  140.             bk2 = ( i+1 < i_max ) ? c*b_ve[i+1] : 0.0;
  141.             z  = ( i+1 < i_max ) ? -s*b_ve[i+1] : 0.0;
  142.             a_ve[i] = ak1;    a_ve[i+1] = ak2;
  143.             b_ve[i] = bk1;
  144.             if ( i < i_max-1 )
  145.                 b_ve[i+1] = bk2;
  146.             if ( i > i_min )
  147.                 b_ve[i-1] = bk;
  148.             if ( Q )
  149.                 rot_cols(Q,i,i+1,c,-s,Q);
  150.                 /* printf("# a [temp2] =\n");    v_output(a); */
  151.                 /* printf("# b [temp2] =\n");    v_output(b); */
  152.             }
  153.  
  154.             /* test to see if matrix should be split */
  155.             for ( i = i_min; i < i_max; i++ )
  156.             if ( fabs(b_ve[i]) < MACHEPS*
  157.                     (fabs(a_ve[i])+fabs(a_ve[i+1])) )
  158.             {   b_ve[i] = 0.0;    split = TRUE;    }
  159.  
  160.             /* printf("# a =\n");    v_output(a); */
  161.             /* printf("# b =\n");    v_output(b); */
  162.         }
  163.     }
  164.  
  165.     return a;
  166. }
  167.  
  168. /* symmeig -- computes eigenvalues of a dense symmetric matrix
  169.     -- A **must** be symmetric on entry
  170.     -- eigenvalues stored in out
  171.     -- Q contains orthogonal matrix of eigenvectors
  172.     -- returns vector of eigenvalues */
  173. VEC    *symmeig(A,Q,out)
  174. MAT    *A, *Q;
  175. VEC    *out;
  176. {
  177.     int    i;
  178.     static MAT    *tmp = MNULL;
  179.     static VEC    *b   = VNULL, *diag = VNULL, *beta = VNULL;
  180.  
  181.     if ( ! A )
  182.         error(E_NULL,"symmeig");
  183.     if ( A->m != A->n )
  184.         error(E_SQUARE,"symmeig");
  185.     if ( ! out || out->dim != A->m )
  186.         out = v_resize(out,A->m);
  187.  
  188.     tmp  = m_copy(A,tmp);
  189.     b    = v_resize(b,A->m - 1);
  190.     diag = v_resize(diag,(u_int)A->m);
  191.     beta = v_resize(beta,(u_int)A->m);
  192.     MEM_STAT_REG(tmp,TYPE_MAT);
  193.     MEM_STAT_REG(b,TYPE_VEC);
  194.     MEM_STAT_REG(diag,TYPE_VEC);
  195.     MEM_STAT_REG(beta,TYPE_VEC);
  196.  
  197.     Hfactor(tmp,diag,beta);
  198.     if ( Q )
  199.         makeHQ(tmp,diag,beta,Q);
  200.  
  201.     for ( i = 0; i < A->m - 1; i++ )
  202.     {
  203.         out->ve[i] = tmp->me[i][i];
  204.         b->ve[i] = tmp->me[i][i+1];
  205.     }
  206.     out->ve[i] = tmp->me[i][i];
  207.     trieig(out,b,Q);
  208.  
  209.     return out;
  210. }
  211.  
  212.