home *** CD-ROM | disk | FTP | other *** search
/ Stars of Shareware: Programmierung / SOURCE.mdf / programm / msdos / c / cephes22 / eigen / eigtst.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-11-17  |  1.7 KB  |  116 lines
  1. /* Test program for eigens.c
  2.  */
  3.  
  4. #define N 11
  5. #define DEBUG 0
  6. static double A[N*N], A0[N*N], X[N*N], D[N*N];
  7. static double RR[N*N];
  8. static double B[N], Y[N];
  9. static int IPS[N];
  10. double maxoffd(), sqrt();
  11.  
  12. main()
  13. {
  14. long ntrials;
  15. int i, j, k;
  16. double x, e, merr;
  17.  
  18. merr = 0.0;
  19. ntrials = 0;
  20.  
  21. loop:
  22.  
  23. /* Fill A[] with random numbers */
  24. k = N*(N+1)/2;
  25. for( i=0; i<k; i++ )
  26.     {
  27.     drand( &x );
  28.     x -= 1.5;
  29.     x = x + x;
  30.     A[i] = x;
  31.     }
  32. /* Unpack A into a symmetric matrix A0 */
  33. tritosquare( N, A, A0 );
  34.  
  35. eigens( A, RR, B, N );
  36. mtransp( N, RR, RR );
  37.  
  38. #if DEBUG
  39. /* Print the computed eigenvalues */
  40. for( i=0; i<N; i++ )
  41.     printf( "%6.2f", B[i] );
  42. /* Print the reduced triangular matrix */
  43. printf( "\n" );
  44. k = 0;
  45. for( i=0; i<N; i++ )
  46.     {
  47.     for( j=0; j<=i; j++ )
  48.         {
  49.         printf( "%6.2f", A[k] );
  50.         k += 1;
  51.         }
  52.     printf( "\n" );
  53.     }
  54.  
  55. /* Print the original symmetric matrix */
  56. printf( "\n" );
  57. for( i=0; i<N; i++ )
  58.     {
  59.     for( j=0; j<N; j++ )
  60.         printf( "%6.2f", A0[N*i+j] );
  61.     printf( "\n" );
  62.     }
  63.  
  64. /* Print the matrix of eigenvectors */
  65. printf( "\n" );
  66. for( i=0; i<N; i++ )
  67.     {
  68.     for( j=0; j<N; j++ )
  69.         printf( "%6.2f", RR[N*i+j] );
  70.     printf( "\n" );
  71.     }
  72.  
  73. #endif
  74.  
  75. /* Multiply the original matrix by the eigenvectors */
  76. mmmpy( N, N, A0, RR, D );
  77. /* Divide each result vector by the corresponding eigenvalue. */
  78. k = 0;
  79. for( i=0; i<N; i++ )
  80.     {
  81.     x = B[i];
  82.     for( j=0; j<N; j++ )
  83.         {
  84.         D[k] /= x;
  85.         k += 1;
  86.         }
  87.     }
  88.  
  89. #if DEBUG
  90. printf( "\n" );
  91. for( i=0; i<N; i++ )
  92.     {
  93.     for( j=0; j<N; j++ )
  94.         printf( "%6.2f", D[N*i+j] );
  95.     printf( "\n" );
  96.     }
  97. #endif
  98.  
  99. /* Compare the result with the computed eigenvector */
  100. mtransp( N, D, D );
  101. e = 0.0;
  102. for( i=0; i<N*N; i++ )
  103.     {
  104.     x = D[i] - RR[i];
  105.     if( x < 0 )
  106.         x = -x;
  107.     if( x > e )
  108.         e = x;
  109.     }
  110. if( e > merr )
  111.     merr = e;
  112. ntrials += 1;
  113. printf( "%4ld %.3e\n", ntrials, merr );
  114. goto loop;
  115. }
  116.