summaryrefslogtreecommitdiff
path: root/libm/double/eigens.c
diff options
context:
space:
mode:
authorEric Andersen <andersen@codepoet.org>2001-05-10 00:40:28 +0000
committerEric Andersen <andersen@codepoet.org>2001-05-10 00:40:28 +0000
commit1077fa4d772832f77a677ce7fb7c2d513b959e3f (patch)
tree579bee13fb0b58d2800206366ec2caecbb15f3fc /libm/double/eigens.c
parent22358dd7ce7bb49792204b698f01a6f69b9c8e08 (diff)
uClibc now has a math library. muahahahaha!
-Erik
Diffstat (limited to 'libm/double/eigens.c')
-rw-r--r--libm/double/eigens.c181
1 files changed, 181 insertions, 0 deletions
diff --git a/libm/double/eigens.c b/libm/double/eigens.c
new file mode 100644
index 000000000..4035e76a1
--- /dev/null
+++ b/libm/double/eigens.c
@@ -0,0 +1,181 @@
+/* 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];
+ }
+}