Output of file : EIGENS.C contained in archive : CEPHES22.ZIP
/* 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.
*
*/
/*
Copyright 1973, 1991 by Stephen L. Moshier
Copyleft version.
*/

void eigens( A, RR, E, N )
double A[], RR[], E[];
int N;
{
int IND, L, LL, LM, M, MM, MQ, I, J, K, IA, LQ;
int IQ, IM, IL, NLI, NMI;
double ANORM, ANORMX, AIA, THR, ALM, QI, ALL, AMM, X, Y;
double SINX, SINX2, COSX, COSX2, SINCS, AIL, AIM;
double RLI, RMI, Q, V;
double sqrt(), fabs();
static double RANGE = 1.0e-10; /*3.0517578e-5;*/

/* Initialize identity matrix in RR[] */
for( J=0; J RR[J] = 0.0;
MM = 0;
for( J=0; J {
RR[MM + J] = 1.0;
MM += N;
}

ANORM=0.0;
for( I=0; I {
for( J=0; 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 {

for( M=L+1; 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 {
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];
}
}

3 Responses to “Category : C Source CodeArchive   : CEPHES22.ZIPFilename : EIGENS.C”

1. Very nice! Thank you for this wonderful archive. I wonder why I found it only now. Long live the BBS file archives!

2. This is so awesome! 😀 I’d be cool if you could download an entire archive of this at once, though.

3. But one thing that puzzles me is the “mtswslnkmcjklsdlsbdmMICROSOFT” string. There is an article about it here. It is definitely worth a read: http://www.os2museum.com/wp/mtswslnk/