/* Test program for eigens.c
 */

#define N 11
#define DEBUG 0
static double A[N*N], A0[N*N], X[N*N], D[N*N];
static double RR[N*N];
static double B[N], Y[N];
static int IPS[N];
double maxoffd(), sqrt();

main()
{
long ntrials;
int i, j, k;
double x, e, merr;

merr = 0.0;
ntrials = 0;

loop:

/* Fill A[] with random numbers */
k = N*(N+1)/2;
for( i=0; i<k; i++ )
 {
 drand( &x );
 x -= 1.5;
 x = x + x;
 A[i] = x;
 }
/* Unpack A into a symmetric matrix A0 */
tritosquare( N, A, A0 );

eigens( A, RR, B, N );
mtransp( N, RR, RR );

#if DEBUG
/* Print the computed eigenvalues */
for( i=0; i<N; i++ )
 printf( "%6.2f", B[i] );
/* Print the reduced triangular matrix */
printf( "\n" );
k = 0;
for( i=0; i<N; i++ )
 {
 for( j=0; j<=i; j++ )
  {
  printf( "%6.2f", A[k] );
  k += 1;
  }
 printf( "\n" );
 }

/* Print the original symmetric matrix */
printf( "\n" );
for( i=0; i<N; i++ )
 {
 for( j=0; j<N; j++ )
  printf( "%6.2f", A0[N*i+j] );
 printf( "\n" );
 }

/* Print the matrix of eigenvectors */
printf( "\n" );
for( i=0; i<N; i++ )
 {
 for( j=0; j<N; j++ )
  printf( "%6.2f", RR[N*i+j] );
 printf( "\n" );
 }

#endif

/* Multiply the original matrix by the eigenvectors */
mmmpy( N, N, A0, RR, D );
/* Divide each result vector by the corresponding eigenvalue. */
k = 0;
for( i=0; i<N; i++ )
 {
 x = B[i];
 for( j=0; j<N; j++ )
  {
  D[k] /= x;
  k += 1;
  }
 }

#if DEBUG
printf( "\n" );
for( i=0; i<N; i++ )
 {
 for( j=0; j<N; j++ )
  printf( "%6.2f", D[N*i+j] );
 printf( "\n" );
 }
#endif

/* Compare the result with the computed eigenvector */
mtransp( N, D, D );
e = 0.0;
for( i=0; i<N*N; i++ )
 {
 x = D[i] - RR[i];
 if( x < 0 )
  x = -x;
 if( x > e )
  e = x;
 }
if( e > merr )
 merr = e;
ntrials += 1;
printf( "%4ld %.3e\n", ntrials, merr );
goto loop;
}