diff options
Diffstat (limited to 'libm/double/eigens.c')
| -rw-r--r-- | libm/double/eigens.c | 181 |
1 files changed, 0 insertions, 181 deletions
diff --git a/libm/double/eigens.c b/libm/double/eigens.c deleted file mode 100644 index 4035e76a1..000000000 --- a/libm/double/eigens.c +++ /dev/null @@ -1,181 +0,0 @@ -/* eigens.c - * - * Eigenvalues and eigenvectors of a real symmetric matrix - * - * - * - * SYNOPSIS: - * - * int n; - * double A[n*(n+1)/2], EV[n*n], E[n]; - * void eigens( A, EV, E, n ); - * - * - * - * DESCRIPTION: - * - * The algorithm is due to J. vonNeumann. - * - * A[] is a symmetric matrix stored in lower triangular form. - * That is, A[ row, column ] = A[ (row*row+row)/2 + column ] - * or equivalently with row and column interchanged. The - * indices row and column run from 0 through n-1. - * - * EV[] is the output matrix of eigenvectors stored columnwise. - * That is, the elements of each eigenvector appear in sequential - * memory order. The jth element of the ith eigenvector is - * EV[ n*i+j ] = EV[i][j]. - * - * E[] is the output matrix of eigenvalues. The ith element - * of E corresponds to the ith eigenvector (the ith row of EV). - * - * On output, the matrix A will have been diagonalized and its - * orginal contents are destroyed. - * - * ACCURACY: - * - * The error is controlled by an internal parameter called RANGE - * which is set to 1e-10. After diagonalization, the - * off-diagonal elements of A will have been reduced by - * this factor. - * - * ERROR MESSAGES: - * - * None. - * - */ - -#include <math.h> -#ifdef ANSIPROT -extern double sqrt ( double ); -extern double fabs ( double ); -#else -double sqrt(), fabs(); -#endif - -void eigens( A, RR, E, N ) -double A[], RR[], E[]; -int N; -{ -int IND, L, LL, LM, M, MM, MQ, I, J, IA, LQ; -int IQ, IM, IL, NLI, NMI; -double ANORM, ANORMX, AIA, THR, ALM, ALL, AMM, X, Y; -double SINX, SINX2, COSX, COSX2, SINCS, AIL, AIM; -double RLI, RMI; -static double RANGE = 1.0e-10; /*3.0517578e-5;*/ - - -/* Initialize identity matrix in RR[] */ -for( J=0; J<N*N; J++ ) - RR[J] = 0.0; -MM = 0; -for( J=0; J<N; J++ ) - { - RR[MM + J] = 1.0; - MM += N; - } - -ANORM=0.0; -for( I=0; I<N; I++ ) - { - for( J=0; J<N; J++ ) - { - if( I != J ) - { - IA = I + (J*J+J)/2; - AIA = A[IA]; - ANORM += AIA * AIA; - } - } - } -if( ANORM <= 0.0 ) - goto done; -ANORM = sqrt( ANORM + ANORM ); -ANORMX = ANORM * RANGE / N; -THR = ANORM; - -while( THR > ANORMX ) -{ -THR=THR/N; - -do -{ /* while IND != 0 */ -IND = 0; - -for( L=0; L<N-1; L++ ) - { - -for( M=L+1; M<N; M++ ) - { - MQ=(M*M+M)/2; - LM=L+MQ; - ALM=A[LM]; - if( fabs(ALM) < THR ) - continue; - - IND=1; - LQ=(L*L+L)/2; - LL=L+LQ; - MM=M+MQ; - ALL=A[LL]; - AMM=A[MM]; - X=(ALL-AMM)/2.0; - Y=-ALM/sqrt(ALM*ALM+X*X); - if(X < 0.0) - Y=-Y; - SINX = Y / sqrt( 2.0 * (1.0 + sqrt( 1.0-Y*Y)) ); - SINX2=SINX*SINX; - COSX=sqrt(1.0-SINX2); - COSX2=COSX*COSX; - SINCS=SINX*COSX; - -/* ROTATE L AND M COLUMNS */ -for( I=0; I<N; I++ ) - { - IQ=(I*I+I)/2; - if( (I != M) && (I != L) ) - { - if(I > M) - IM=M+IQ; - else - IM=I+MQ; - if(I >= L) - IL=L+IQ; - else - IL=I+LQ; - AIL=A[IL]; - AIM=A[IM]; - X=AIL*COSX-AIM*SINX; - A[IM]=AIL*SINX+AIM*COSX; - A[IL]=X; - } - NLI = N*L + I; - NMI = N*M + I; - RLI = RR[ NLI ]; - RMI = RR[ NMI ]; - RR[NLI]=RLI*COSX-RMI*SINX; - RR[NMI]=RLI*SINX+RMI*COSX; - } - - X=2.0*ALM*SINCS; - A[LL]=ALL*COSX2+AMM*SINX2-X; - A[MM]=ALL*SINX2+AMM*COSX2+X; - A[LM]=(ALL-AMM)*SINCS+ALM*(COSX2-SINX2); - } /* for M=L+1 to N-1 */ - } /* for L=0 to N-2 */ - - } -while( IND != 0 ); - -} /* while THR > ANORMX */ - -done: ; - -/* Extract eigenvalues from the reduced matrix */ -L=0; -for( J=1; J<=N; J++ ) - { - L=L+J; - E[J-1]=A[L-1]; - } -} |