home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD2.mdf / c / library / dos / math / cephes / eigen / eigens.c < prev    next >
Encoding:
Text File  |  1992-11-17  |  3.2 KB  |  179 lines
  1. /*                            eigens.c
  2.  *
  3.  *    Eigenvalues and eigenvectors of a real symmetric matrix
  4.  *
  5.  *
  6.  *
  7.  * SYNOPSIS:
  8.  *
  9.  * int n;
  10.  * double A[n*(n+1)/2], EV[n*n], E[n];
  11.  * void eigens( A, EV, E, n );
  12.  *
  13.  *
  14.  *
  15.  * DESCRIPTION:
  16.  *
  17.  * The algorithm is due to J. vonNeumann.
  18.  *                   -     -
  19.  * A[] is a symmetric matrix stored in lower triangular form.
  20.  * That is, A[ row, column ] = A[ (row*row+row)/2 + column ]
  21.  * or equivalently with row and column interchanged.  The
  22.  * indices row and column run from 0 through n-1.
  23.  *
  24.  * EV[] is the output matrix of eigenvectors stored columnwise.
  25.  * That is, the elements of each eigenvector appear in sequential
  26.  * memory order.  The jth element of the ith eigenvector is
  27.  * EV[ n*i+j ] = EV[i][j].
  28.  *
  29.  * E[] is the output matrix of eigenvalues.  The ith element
  30.  * of E corresponds to the ith eigenvector (the ith row of EV).
  31.  *
  32.  * On output, the matrix A will have been diagonalized and its
  33.  * orginal contents are destroyed.
  34.  *
  35.  * ACCURACY:
  36.  *
  37.  * The error is controlled by an internal parameter called RANGE
  38.  * which is set to 1e-10.  After diagonalization, the
  39.  * off-diagonal elements of A will have been reduced by
  40.  * this factor.
  41.  *
  42.  * ERROR MESSAGES:
  43.  *
  44.  * None.
  45.  *
  46.  */
  47. /*
  48. Copyright 1973, 1991 by Stephen L. Moshier
  49. Copyleft version.
  50. */
  51.  
  52. void eigens( A, RR, E, N )
  53. double A[], RR[], E[];
  54. int N;
  55. {
  56. int IND, L, LL, LM, M, MM, MQ, I, J, K, IA, LQ;
  57. int IQ, IM, IL, NLI, NMI;
  58. double ANORM, ANORMX, AIA, THR, ALM, QI, ALL, AMM, X, Y;
  59. double SINX, SINX2, COSX, COSX2, SINCS, AIL, AIM;
  60. double RLI, RMI, Q, V;
  61. double sqrt(), fabs();
  62. static double RANGE = 1.0e-10; /*3.0517578e-5;*/
  63.  
  64.  
  65. /* Initialize identity matrix in RR[] */
  66. for( J=0; J<N*N; J++ )
  67.     RR[J] = 0.0;
  68. MM = 0;
  69. for( J=0; J<N; J++ )
  70.     {
  71.     RR[MM + J] = 1.0;
  72.     MM += N;
  73.     }
  74.  
  75. ANORM=0.0;
  76. for( I=0; I<N; I++ )
  77.     {
  78.     for( J=0; J<N; J++ )
  79.         {
  80.         if( I != J )
  81.             {
  82.             IA = I + (J*J+J)/2;
  83.             AIA = A[IA];
  84.             ANORM += AIA * AIA;
  85.             }
  86.         }
  87.     }
  88. if( ANORM <= 0.0 )
  89.     goto done;
  90. ANORM = sqrt( ANORM + ANORM );
  91. ANORMX = ANORM * RANGE / N;
  92. THR = ANORM;
  93.  
  94. while( THR > ANORMX )
  95. {
  96. THR=THR/N;
  97.  
  98. do
  99. { /* while IND != 0 */
  100. IND = 0;
  101.  
  102. for( L=0; L<N-1; L++ )
  103.     {
  104.  
  105. for( M=L+1; M<N; M++ )
  106.     {
  107.     MQ=(M*M+M)/2;
  108.     LM=L+MQ;
  109.     ALM=A[LM];
  110.     if( fabs(ALM) < THR )
  111.         continue;
  112.  
  113.     IND=1;
  114.     LQ=(L*L+L)/2;
  115.     LL=L+LQ;
  116.     MM=M+MQ;
  117.     ALL=A[LL];
  118.     AMM=A[MM];
  119.     X=(ALL-AMM)/2.0;
  120.     Y=-ALM/sqrt(ALM*ALM+X*X);
  121.     if(X < 0.0)
  122.         Y=-Y;
  123.     SINX = Y / sqrt( 2.0 * (1.0 + sqrt( 1.0-Y*Y)) );
  124.     SINX2=SINX*SINX;
  125.     COSX=sqrt(1.0-SINX2);
  126.     COSX2=COSX*COSX;
  127.     SINCS=SINX*COSX;
  128.  
  129. /*       ROTATE L AND M COLUMNS */
  130. for( I=0; I<N; I++ )
  131.     {
  132.     IQ=(I*I+I)/2;
  133.     if( (I != M) && (I != L) )
  134.         {
  135.         if(I > M)
  136.             IM=M+IQ;
  137.         else
  138.             IM=I+MQ;
  139.         if(I >= L)
  140.             IL=L+IQ;
  141.         else
  142.             IL=I+LQ;
  143.         AIL=A[IL];
  144.         AIM=A[IM];
  145.         X=AIL*COSX-AIM*SINX;
  146.         A[IM]=AIL*SINX+AIM*COSX;
  147.         A[IL]=X;
  148.         }
  149.     NLI = N*L + I;
  150.     NMI = N*M + I;
  151.     RLI = RR[ NLI ];
  152.     RMI = RR[ NMI ];
  153.     RR[NLI]=RLI*COSX-RMI*SINX;
  154.     RR[NMI]=RLI*SINX+RMI*COSX;
  155.     }
  156.  
  157.     X=2.0*ALM*SINCS;
  158.     A[LL]=ALL*COSX2+AMM*SINX2-X;
  159.     A[MM]=ALL*SINX2+AMM*COSX2+X;
  160.     A[LM]=(ALL-AMM)*SINCS+ALM*(COSX2-SINX2);
  161.     } /* for M=L+1 to N-1 */
  162.     } /* for L=0 to N-2 */
  163.  
  164.     }
  165. while( IND != 0 );
  166.  
  167. } /* while THR > ANORMX */
  168.  
  169. done:    ;
  170.  
  171. /* Extract eigenvalues from the reduced matrix */
  172. L=0;
  173. for( J=1; J<=N; J++ )
  174.     {
  175.     L=L+J;
  176.     E[J-1]=A[L-1];
  177.     }
  178. }
  179.