fork download
copy 
  1. #include <stdio.h>
  2. #include <math.h>
  3.  
  4. #define N 4         
  5. #define EPS 1e-6    // Convergence threshold
  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.     // Initial symmetric matrix A [cite: 13]
  19.     double A[N][N] = {
  20.         {5.0, 4.0, 1.0, 1.0},
  21.         {4.0, 5.0, 1.0, 1.0},
  22.         {1.0, 1.0, 4.0, 2.0},
  23.         {1.0, 1.0, 2.0, 4.0}
  24.     };
  25.  
  26.     // Backup the original matrix for verification (Ax = lambda * x)
  27.     double A_orig[N][N];
  28.     for (int i = 0; i < N; i++)
  29.         for (int j = 0; j < N; j++)
  30.             A_orig[i][j] = A[i][j];
  31.  
  32.     // Initialize P as an Identity Matrix. Columns will store eigenvectors[cite: 39].
  33.     double P[N][N] = {0};
  34.     for (int i = 0; i < N; i++) P[i][i] = 1.0;
  35.  
  36.     int iter = 0;
  37.  
  38.     // Tidy, aligned table header for convergence history
  39.     printf("=== CONVERGENCE HISTORY ===\n");
  40.     printf("%-5s   %-16s   %-8s\n", "Iter", "Max Off-Diagonal", "Position");
  41.     printf("---------------------------------------\n");
  42.  
  43.     while (iter < MAX_ITER) {
  44.         // 1. Find the maximum off-diagonal element A[p][q] [cite: 76]
  45.         int p = 0, q = 1;
  46.         double max_val = fabs(A[0][1]);
  47.         for (int i = 0; i < N; i++) {
  48.             for (int j = i + 1; j < N; j++) {
  49.                 if (fabs(A[i][j]) > max_val) {
  50.                     max_val = fabs(A[i][j]);
  51.                     p = i;
  52.                     q = j;
  53.                 }
  54.             }
  55.         }
  56.  
  57.         // Tabled log row output
  58.         printf("%-5d   %-16.6f   A[%d][%d]\n", iter, max_val, p, q);
  59.  
  60.         // Check convergence status [cite: 76]
  61.         if (max_val < EPS) {
  62.             printf("---------------------------------------\n");
  63.             printf("Status: Successfully converged.\n\n");
  64.             break;
  65.         }
  66.  
  67.         // 2. Calculate rotation angle using acos(-1.0) dynamically [cite: 54]
  68.         double phi, cos_t, sin_t;
  69.         if (fabs(A[p][p] - A[q][q]) < 1e-12) {
  70.             phi = acos(-1.0) / 4.0; 
  71.         } else {
  72.             phi = 0.5 * atan2(2.0 * A[p][q], A[p][p] - A[q][q]);
  73.         }
  74.         cos_t = cos(phi);
  75.         sin_t = sin(phi);
  76.  
  77.         // 3. Update Matrix A [cite: 51]
  78.         double Ap_old = A[p][p];
  79.         double Aq_old = A[q][q];
  80.  
  81.         A[p][p] = Ap_old * cos_t * cos_t + Aq_old * sin_t * sin_t + 2.0 * A[p][q] * sin_t * cos_t;
  82.         A[q][q] = Ap_old * sin_t * sin_t + Aq_old * cos_t * cos_t - 2.0 * A[p][q] * sin_t * cos_t;
  83.         A[p][q] = A[q][p] = 0.0; 
  84.  
  85.         for (int i = 0; i < N; i++) {
  86.             if (i != p && i != q) {
  87.                 double a_ip = A[i][p];
  88.                 double a_iq = A[i][q];
  89.                 A[i][p] = A[p][i] = a_ip * cos_t + a_iq * sin_t;
  90.                 A[i][q] = A[q][i] = -a_ip * sin_t + a_iq * cos_t;
  91.             }
  92.         }
  93.  
  94.         // 4. Update Cumulative Eigenvector Matrix P [cite: 77]
  95.         for (int i = 0; i < N; i++) {
  96.             double p_ip = P[i][p];
  97.             double p_iq = P[i][q];
  98.             P[i][p] = p_ip * cos_t + p_iq * sin_t;
  99.             P[i][q] = -p_ip * sin_t + p_iq * cos_t;
  100.         }
  101.  
  102.         iter++;
  103.     }
  104.  
  105.     // === Output Final Results ===
  106.     printf("=== EIGENVALUES & EIGENVECTORS ===\n");
  107.     for (int j = 0; j < N; j++) {
  108.         printf("Eigenvalue %d (固有値): %f\n", j + 1, A[j][j]); // [cite: 24]
  109.         printf("Eigenvector %d (固有ベクトル):\n[ ", j + 1);
  110.         for (int i = 0; i < N; i++) {
  111.             printf("%f ", P[i][j]); // Prints column vectors horizontally [cite: 39]
  112.         }
  113.         printf("]\n\n");
  114.     }
  115.  
  116.     // === Automated Verification Step (Ax = lambda * x) ===
  117.     printf("=== 固有対の検証チェック (A*x - lambda*x) ===\n");
  118.  
  119.     for (int j = 0; j < N; j++) {
  120.         double lambda = A[j][j];
  121.         printf("Eigenvalue %d (%f):\n", j + 1, lambda);
  122.  
  123.         for (int i = 0; i < N; i++) {
  124.             // Compute elements of Ax
  125.             double Ax_i = 0.0;
  126.             for (int k = 0; k < N; k++) {
  127.                 Ax_i += A_orig[i][k] * P[k][j];
  128.             }
  129.             // Compute elements of lambda * x
  130.             double lambda_x_i = lambda * P[i][j];
  131.  
  132.             // Calculate absolute residual error
  133.             double residual = fabs(Ax_i - lambda_x_i);
  134.  
  135.             // Print out the numerical matching comparison with scientific notation residuals
  136.             printf("  第 %d 行: A*x = %9.6f, lambda*x = %9.6f, 残差絶対値 = %e\n", 
  137.                    i + 1, Ax_i, lambda_x_i, residual);
  138.         }
  139.         printf("\n");
  140.     }
  141.  
  142.     return 0;
  143. }
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.000000
Eigenvector 1 (固有ベクトル):
[ 0.632456 0.632456 0.316228 0.316228 ]

Eigenvalue 2 (固有値): 1.000000
Eigenvector 2 (固有ベクトル):
[ -0.707107 0.707107 0.000000 0.000000 ]

Eigenvalue 3 (固有値): 5.000000
Eigenvector 3 (固有ベクトル):
[ -0.316228 -0.316228 0.632456 0.632456 ]

Eigenvalue 4 (固有値): 2.000000
Eigenvector 4 (固有ベクトル):
[ 0.000000 0.000000 -0.707107 0.707107 ]

=== 固有対の検証チェック (A*x - lambda*x) ===
Eigenvalue 1 (10.000000):
  第 1 行: A*x =  6.324555, lambda*x =  6.324555, 残差絶対値 = 8.881784e-16
  第 2 行: A*x =  6.324555, lambda*x =  6.324555, 残差絶対値 = 1.776357e-15
  第 3 行: A*x =  3.162278, lambda*x =  3.162278, 残差絶対値 = 4.440892e-16
  第 4 行: A*x =  3.162278, lambda*x =  3.162278, 残差絶対値 = 8.881784e-16

Eigenvalue 2 (1.000000):
  第 1 行: A*x = -0.707107, lambda*x = -0.707107, 残差絶対値 = 4.440892e-16
  第 2 行: A*x =  0.707107, lambda*x =  0.707107, 残差絶対値 = 3.330669e-16
  第 3 行: A*x =  0.000000, lambda*x =  0.000000, 残差絶対値 = 1.110223e-16
  第 4 行: A*x =  0.000000, lambda*x =  0.000000, 残差絶対値 = 1.110223e-16

Eigenvalue 3 (5.000000):
  第 1 行: A*x = -1.581139, lambda*x = -1.581139, 残差絶対値 = 4.440892e-16
  第 2 行: A*x = -1.581139, lambda*x = -1.581139, 残差絶対値 = 2.220446e-16
  第 3 行: A*x =  3.162278, lambda*x =  3.162278, 残差絶対値 = 4.440892e-16
  第 4 行: A*x =  3.162278, lambda*x =  3.162278, 残差絶対値 = 0.000000e+00

Eigenvalue 4 (2.000000):
  第 1 行: A*x =  0.000000, lambda*x =  0.000000, 残差絶対値 = 1.110223e-16
  第 2 行: A*x =  0.000000, lambda*x =  0.000000, 残差絶対値 = 1.110223e-16
  第 3 行: A*x = -1.414214, lambda*x = -1.414214, 残差絶対値 = 2.220446e-16
  第 4 行: A*x =  1.414214, lambda*x =  1.414214, 残差絶対値 = 2.220446e-16
https://ideone.com/2lojcK
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!