fork download
copy 
  1. #include <stdio.h>
  2. #include <math.h>
  3.  
  4. #define N 4         
  5. #define EPS 1e-6    
  6. #define MAX_ITER 100
  7.  
  8. void print_matrix(double mat[N][N]) {
  9.     for (int i = 0; i < N; i++) {
  10.         for (int j = 0; j < N; j++) {
  11.             printf("%7.4f ", mat[i][j]);
  12.         }
  13.         printf("\n");
  14.     }
  15. }
  16.  
  17. int main() {
  18.     double A[N][N] = {
  19.         {5.0, 4.0, 1.0, 1.0},
  20.         {4.0, 5.0, 1.0, 1.0},
  21.         {1.0, 1.0, 4.0, 2.0},
  22.         {1.0, 1.0, 2.0, 4.0}
  23.     };
  24.  
  25.     double A_orig[N][N];
  26.     for (int i = 0; i < N; i++)
  27.         for (int j = 0; j < N; j++)
  28.             A_orig[i][j] = A[i][j];
  29.  
  30.     double P[N][N] = {0};
  31.     for (int i = 0; i < N; i++) P[i][i] = 1.0;
  32.  
  33.     int iter = 0;
  34.  
  35.     // Clean Tabled Header for Convergence History [cite: 13]
  36.     printf("=== CONVERGENCE HISTORY ===\n");
  37.     printf("%-5s   %-16s   %-8s\n", "Iter", "Max Off-Diagonal", "Position");
  38.     printf("---------------------------------------\n");
  39.  
  40.     while (iter < MAX_ITER) {
  41.         int p = 0, q = 1;
  42.         double max_val = fabs(A[0][1]);
  43.         for (int i = 0; i < N; i++) {
  44.             for (int j = i + 1; j < N; j++) {
  45.                 if (fabs(A[i][j]) > max_val) {
  46.                     max_val = fabs(A[i][j]);
  47.                     p = i;
  48.                     q = j;
  49.                 }
  50.             }
  51.         }
  52.  
  53.         // Tidy row output alignment
  54.         printf("%-5d   %-16.6f   A[%d][%d]\n", iter, max_val, p, q);
  55.  
  56.         if (max_val < EPS) {
  57.             printf("---------------------------------------\n");
  58.             printf("Status: Successfully converged.\n\n");
  59.             break;
  60.         }
  61.  
  62.         double phi, cos_t, sin_t;
  63.         if (fabs(A[p][p] - A[q][q]) < 1e-12) {
  64.             phi = acos(-1.0) / 4.0; 
  65.         } else {
  66.             phi = 0.5 * atan2(2.0 * A[p][q], A[p][p] - A[q][q]);
  67.         }
  68.         cos_t = cos(phi);
  69.         sin_t = sin(phi);
  70.  
  71.         double Ap_old = A[p][p];
  72.         double Aq_old = A[q][q];
  73.  
  74.         A[p][p] = Ap_old * cos_t * cos_t + Aq_old * sin_t * sin_t + 2.0 * A[p][q] * sin_t * cos_t;
  75.         A[q][q] = Ap_old * sin_t * sin_t + Aq_old * cos_t * cos_t - 2.0 * A[p][q] * sin_t * cos_t;
  76.         A[p][q] = A[q][p] = 0.0; 
  77.  
  78.         for (int i = 0; i < N; i++) {
  79.             if (i != p && i != q) {
  80.                 double a_ip = A[i][p];
  81.                 double a_iq = A[i][q];
  82.                 A[i][p] = A[p][i] = a_ip * cos_t + a_iq * sin_t;
  83.                 A[i][q] = A[q][i] = -a_ip * sin_t + a_iq * cos_t;
  84.             }
  85.         }
  86.  
  87.         for (int i = 0; i < N; i++) {
  88.             double p_ip = P[i][p];
  89.             double p_iq = P[i][q];
  90.             P[i][p] = p_ip * cos_t + p_iq * sin_t;
  91.             P[i][q] = -p_ip * sin_t + p_iq * cos_t;
  92.         }
  93.  
  94.         iter++;
  95.     }
  96.  
  97.     // === Output Final Results ===
  98.     printf("=== EIGENVALUES & EIGENVECTORS ===\n");
  99.     for (int j = 0; j < N; j++) {
  100.         printf("Eigenvalue %d: %.4f\n", j + 1, A[j][j]);
  101.         printf("Eigenvector %d:\n", j + 1);
  102.         for (int i = 0; i < N; i++) {
  103.             printf("  %7.4f\n", P[i][j]); 
  104.         }
  105.         printf("\n");
  106.     }
  107.  
  108.     // === Automated Verification Step ===
  109.     printf("=== CODE-LEVEL VERIFICATION (Ax = lambda * x) ===\n");
  110.  
  111.     int all_passed = 1;
  112.     double tolerance = 1e-5; // Threshold for acceptable math verification error
  113.  
  114.     for (int j = 0; j < N; j++) {
  115.         double lambda = A[j][j];
  116.  
  117.         for (int i = 0; i < N; i++) {
  118.             double Ax_i = 0.0;
  119.             for (int k = 0; k < N; k++) {
  120.                 Ax_i += A_orig[i][k] * P[k][j];
  121.             }
  122.             double lambda_x_i = lambda * P[i][j];
  123.  
  124.             // Programmatically verify the absolute difference
  125.             if (fabs(Ax_i - lambda_x_i) > tolerance) {
  126.                 all_passed = 0;
  127.             }
  128.         }
  129.     }
  130.  
  131.     if (all_passed) {
  132.         printf("Verification Status: PASSED (All equations hold within error tolerances)\n");
  133.     } else {
  134.         printf("Verification Status: FAILED (Discrepancy detected)\n");
  135.     }
  136.  
  137.     return 0;
  138. }
  139.  
Success #stdin #stdout 0s 5320KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
=== CONVERGENCE HISTORY ===
Iter    Max Off-Diagonal   Position
---------------------------------------
0       4.000000           A[0][1]
1       2.000000           A[2][3]
2       2.000000           A[0][2]
3       0.000000           A[0][3]
---------------------------------------
Status: Successfully converged.

=== EIGENVALUES & EIGENVECTORS ===
Eigenvalue 1: 10.0000
Eigenvector 1:
   0.6325
   0.6325
   0.3162
   0.3162

Eigenvalue 2: 1.0000
Eigenvector 2:
  -0.7071
   0.7071
   0.0000
   0.0000

Eigenvalue 3: 5.0000
Eigenvector 3:
  -0.3162
  -0.3162
   0.6325
   0.6325

Eigenvalue 4: 2.0000
Eigenvector 4:
   0.0000
   0.0000
  -0.7071
   0.7071

=== CODE-LEVEL VERIFICATION (Ax = lambda * x) ===
Verification Status: PASSED (All equations hold within error tolerances)
https://ideone.com/GdK9Db
language:
C (gcc 8.3)
created:
4 days ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!