home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 22 gnu / 22-gnu.zip / mesch12a.zip / symmeig.c < prev    next >
C/C++ Source or Header  |  1994-01-13  |  6KB  |  209 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.2 1994/01/13 05:27:47 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 = c*b_ve[i_min+1];
  110.             z  = -s*b_ve[i_min+1];
  111.             a_ve[i_min] = ak1;    a_ve[i_min+1] = ak2;
  112.             b_ve[i_min] = bk1;    b_ve[i_min+1] = bk2;
  113.             if ( Q )
  114.             rot_cols(Q,i_min,i_min+1,c,-s,Q);
  115.             /* printf("# z = %g\n",z); */
  116.             /* printf("# a [temp1] =\n");    v_output(a); */
  117.             /* printf("# b [temp1] =\n");    v_output(b); */
  118.  
  119.             for ( i = i_min+1; i < i_max; i++ )
  120.             {
  121.             /* get Givens' rotation for sub-block -- k == i-1 */
  122.             givens(b_ve[i-1],z,&c,&s);
  123.             s = -s;
  124.             /* printf("# c = %g, s = %g\n",c,s); */
  125.  
  126.             /* perform Givens' rotation on sub-block */
  127.                 if ( fabs(c) < SQRT2 )
  128.                 {    c2 = c*c;    s2 = 1-c2;    }
  129.                 else
  130.                 {    s2 = s*s;    c2 = 1-s2;    }
  131.                 cs = c*s;
  132.             bk  = c*b_ve[i-1] - s*z;
  133.             ak1 = c2*a_ve[i]+s2*a_ve[i+1]-2*cs*b_ve[i];
  134.             bk1 = cs*(a_ve[i]-a_ve[i+1]) +
  135.                         (c2-s2)*b_ve[i];
  136.             ak2 = s2*a_ve[i]+c2*a_ve[i+1]+2*cs*b_ve[i];
  137.             bk2 = c*b_ve[i+1];
  138.             z  = -s*b_ve[i+1];
  139.             a_ve[i] = ak1;    a_ve[i+1] = ak2;
  140.             b_ve[i] = bk1;
  141.             if ( i < i_max-1 )
  142.                 b_ve[i+1] = bk2;
  143.             if ( i > i_min )
  144.                 b_ve[i-1] = bk;
  145.             if ( Q )
  146.                 rot_cols(Q,i,i+1,c,-s,Q);
  147.                 /* printf("# a [temp2] =\n");    v_output(a); */
  148.                 /* printf("# b [temp2] =\n");    v_output(b); */
  149.             }
  150.  
  151.             /* test to see if matrix should be split */
  152.             for ( i = i_min; i < i_max; i++ )
  153.             if ( fabs(b_ve[i]) < MACHEPS*
  154.                     (fabs(a_ve[i])+fabs(a_ve[i+1])) )
  155.             {   b_ve[i] = 0.0;    split = TRUE;    }
  156.  
  157.             /* printf("# a =\n");    v_output(a); */
  158.             /* printf("# b =\n");    v_output(b); */
  159.         }
  160.     }
  161.  
  162.     return a;
  163. }
  164.  
  165. /* symmeig -- computes eigenvalues of a dense symmetric matrix
  166.     -- A **must** be symmetric on entry
  167.     -- eigenvalues stored in out
  168.     -- Q contains orthogonal matrix of eigenvectors
  169.     -- returns vector of eigenvalues */
  170. VEC    *symmeig(A,Q,out)
  171. MAT    *A, *Q;
  172. VEC    *out;
  173. {
  174.     int    i;
  175.     static MAT    *tmp = MNULL;
  176.     static VEC    *b   = VNULL, *diag = VNULL, *beta = VNULL;
  177.  
  178.     if ( ! A )
  179.         error(E_NULL,"symmeig");
  180.     if ( A->m != A->n )
  181.         error(E_SQUARE,"symmeig");
  182.     if ( ! out || out->dim != A->m )
  183.         out = v_resize(out,A->m);
  184.  
  185.     tmp  = m_copy(A,tmp);
  186.     b    = v_resize(b,A->m - 1);
  187.     diag = v_resize(diag,(u_int)A->m);
  188.     beta = v_resize(beta,(u_int)A->m);
  189.     MEM_STAT_REG(tmp,TYPE_MAT);
  190.     MEM_STAT_REG(b,TYPE_VEC);
  191.     MEM_STAT_REG(diag,TYPE_VEC);
  192.     MEM_STAT_REG(beta,TYPE_VEC);
  193.  
  194.     Hfactor(tmp,diag,beta);
  195.     if ( Q )
  196.         makeHQ(tmp,diag,beta,Q);
  197.  
  198.     for ( i = 0; i < A->m - 1; i++ )
  199.     {
  200.         out->ve[i] = tmp->me[i][i];
  201.         b->ve[i] = tmp->me[i][i+1];
  202.     }
  203.     out->ve[i] = tmp->me[i][i];
  204.     trieig(out,b,Q);
  205.  
  206.     return out;
  207. }
  208.  
  209.