24 void sortd(
int length,
double * a,
int *ind )
27 int i, ii, ij, j, m, m1, n2;
30 for ( i = 0 ; i < length; i++ ) ind[ i ] = i;
34 while ( m <= n2 ) m = 2 * m;
38 for ( j = m1-1 ; j < length; j++ ) {
55 if ( m > 0 )
goto three;
68 void eig_jacobi(
int n,
double a[3][3],
double v[3][3],
double *d,
int* ind,
int* maxsweep,
int* maxannil,
double* epsilon)
88 int i, j, k, l, p, q, sweep, annil;
89 double alpha, beta, c, s, mu1, mu2, mu3;
90 double pp, qq, pq, offdiag;
96 printf(
"In jacobi(), size of matrix = %d\n", n );
101 for( k = 0; k < n; k++ ) {
102 for( l = k+1; l < n; l++ ) {
103 offdiag = offdiag + a[k][l] * a[k][l];
107 mu1 = sqrt( offdiag ) / n ;
110 for( i = 0; i < n; i++ ) {
111 d[ i ] = a[ i ][ i ];
112 for( j = 0; j < n; j++ ) {
123 if( mu1 <= ( *epsilon * mu1 ) )
goto done;
125 for( sweep = 1; sweep < *maxsweep; sweep++ ) {
127 for( p = 0; p < n; p++ ) {
128 for( q = p+1; q < n; q++ ) {
129 if( fabs( a[ p ][ q ] ) > mu2 ) {
132 alpha = 0.5 * ( d[ p ] - d[ q ] );
133 beta = sqrt( a[ p ][ q ] * a[ p ][ q ] + alpha * alpha );
134 c = sqrt( 0.5 + fabs( alpha ) / ( 2.0 * beta ) );
138 s = - ( alpha * a[ p ][ q ] ) / ( 2.0 * beta * fabs( alpha ) * c );
144 d[ p ] = c * c * pp + s * s * qq - 2.0 * s * c * pq;
145 d[ q ] = s * s * pp + c * c * qq + 2.0 * s * c * pq;
146 a[ p ][ q ] = ( c * c - s * s ) * pq + s * c * ( pp - qq );
148 for( k = 0; k < p; k++ ) {
151 a[ k ][ p ] = c * t1 - s * t2;
152 a[ k ][ q ] = c * t2 + s * t1;
155 for( k = p+1; k < q; k++ ) {
158 a[ p ][ k ] = c * t1 - s * t2;
159 a[ k ][ q ] = c * t2 + s * t1;
162 for( k = q+1; k < n; k++ ) {
165 a[ p ][ k ] = c * t1 - s * t2;
166 a[ q ][ k ] = c * t2 + s * t1;
169 for( k = 0; k < n; k++ ) {
172 v[ p ][ k ] = c * t1 - s * t2;
173 v[ q ][ k ] = s * t1 + c * t2;
180 for( k = 0; k < n; k++ ) {
181 for( l = k+1; l < n; l++ ) {
182 offdiag = offdiag + a[k][l] * a[k][l];
185 mu3 = sqrt( offdiag ) / n;
190 printf(
"Offdiagonal sum is increasing muold= %f munow= %f\n", mu2, mu3 );
196 if( mu2 <= ( *epsilon * mu1 ) )
goto done;
void sortd(int length, double *a, int *ind)
lu decomposition and linear LMS solver
void eig_jacobi(int n, double a[3][3], double v[3][3], double *d, int *ind, int *maxsweep, int *maxannil, double *epsilon)
Implements jacobi eigenvalue/vector algorithm on a symmetric matrix stored as a 2 dimensional matrix ...