fork download
copy 
  1. #include <stdio.h>
  2. #include <math.h>
  3.  
  4. #define N 4         // Matrix size
  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.     // Original matrix A
  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.     // Keep a backup of the original matrix A 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 of P will become our eigenvectors.
  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.     printf("=== Jacobi Method Convergence Process ===\n\n");
  38.  
  39.     while (iter < MAX_ITER) {
  40.         // 1. Find the maximum off-diagonal element A[p][q]
  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.         // Print track of convergence
  54.         printf("Iteration %d: Max off-diagonal = %f (at A[%d][%d])\n", iter, max_val, p, q);
  55.  
  56.         if (max_val < EPS) {
  57.             printf("\n--> Successfully converged after %d iterations.\n\n", iter);
  58.             break;
  59.         }
  60.  
  61.         // 2. Calculate the rotation angle theta (using acos(-1.0) for Pi)
  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.         // 3. Update Matrix A (P_k^T * A * P_k)
  72.         double Ap_old = A[p][p];
  73.         double Aq_old = A[q][q];
  74.  
  75.         A[p][p] = Ap_old * cos_t * cos_t + Aq_old * sin_t * sin_t + 2.0 * A[p][q] * sin_t * cos_t;
  76.         A[q][q] = Ap_old * sin_t * sin_t + Aq_old * cos_t * cos_t - 2.0 * A[p][q] * sin_t * cos_t;
  77.         A[p][q] = A[q][p] = 0.0; 
  78.  
  79.         for (int i = 0; i < N; i++) {
  80.             if (i != p && i != q) {
  81.                 double a_ip = A[i][p];
  82.                 double a_iq = A[i][q];
  83.                 A[i][p] = A[p][i] = a_ip * cos_t + a_iq * sin_t;
  84.                 A[i][q] = A[q][i] = -a_ip * sin_t + a_iq * cos_t;
  85.             }
  86.         }
  87.  
  88.         // 4. Update Eigenvector Matrix P (P = P * P_k)
  89.         // We modify the columns p and q of matrix P
  90.         for (int i = 0; i < N; i++) {
  91.             double p_ip = P[i][p];
  92.             double p_iq = P[i][q];
  93.             P[i][p] = p_ip * cos_t + p_iq * sin_t;
  94.             P[i][q] = -p_ip * sin_t + p_iq * cos_t;
  95.         }
  96.  
  97.         iter++;
  98.     }
  99.  
  100.     // === Output Final Eigenvalues & Eigenvectors ===
  101.     printf("=== Final Results ===\n");
  102.     for (int j = 0; j < N; j++) {
  103.         double lambda = A[j][j]; // Diagonal elements are eigenvalues
  104.         printf("Eigenvalue %d: %7.4f\n", j + 1, lambda);
  105.         printf("Eigenvector %d (Column vector):\n", j + 1);
  106.         for (int i = 0; i < N; i++) {
  107.             printf("  [%7.4f]\n", P[i][j]); // P[i][j] extracts the j-th column
  108.         }
  109.         printf("\n");
  110.     }
  111.  
  112.     // === Verification (Ax = lambda * x) ===
  113.     printf("=== Code-Level Verification (Ax - lambda * x) ===\n");
  114.     printf("If the values below are essentially 0, Ax = lambda * x holds true perfectly.\n\n");
  115.  
  116.     for (int j = 0; j < N; j++) {
  117.         double lambda = A[j][j];
  118.         printf("Verifying Eigenvalue %d (%7.4f):\n", j + 1, lambda);
  119.  
  120.         for (int i = 0; i < N; i++) {
  121.             // Compute i-th row of (A_orig * j-th column of P)
  122.             double Ax_i = 0.0;
  123.             for (int k = 0; k < N; k++) {
  124.                 Ax_i += A_orig[i][k] * P[k][j];
  125.             }
  126.             // Compute i-th row of (lambda * j-th column of P)
  127.             double lambda_x_i = lambda * P[i][j];
  128.  
  129.             // Difference should approach 0
  130.             printf("  Row %d difference: %10.3e\n", i + 1, Ax_i - lambda_x_i);
  131.         }
  132.         printf("\n");
  133.     }
  134.  
  135.     return 0;
  136. }
Success #stdin #stdout 0s 5320KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
=== Jacobi Method Convergence Process ===

Iteration 0: Max off-diagonal = 4.000000 (at A[0][1])
Iteration 1: Max off-diagonal = 2.000000 (at A[2][3])
Iteration 2: Max off-diagonal = 2.000000 (at A[0][2])
Iteration 3: Max off-diagonal = 0.000000 (at A[0][3])

--> Successfully converged after 3 iterations.

=== Final Results ===
Eigenvalue 1: 10.0000
Eigenvector 1 (Column vector):
  [ 0.6325]
  [ 0.6325]
  [ 0.3162]
  [ 0.3162]

Eigenvalue 2:  1.0000
Eigenvector 2 (Column vector):
  [-0.7071]
  [ 0.7071]
  [ 0.0000]
  [ 0.0000]

Eigenvalue 3:  5.0000
Eigenvector 3 (Column vector):
  [-0.3162]
  [-0.3162]
  [ 0.6325]
  [ 0.6325]

Eigenvalue 4:  2.0000
Eigenvector 4 (Column vector):
  [ 0.0000]
  [ 0.0000]
  [-0.7071]
  [ 0.7071]

=== Code-Level Verification (Ax - lambda * x) ===
If the values below are essentially 0, Ax = lambda * x holds true perfectly.

Verifying Eigenvalue 1 (10.0000):
  Row 1 difference:  8.882e-16
  Row 2 difference:  1.776e-15
  Row 3 difference:  4.441e-16
  Row 4 difference:  8.882e-16

Verifying Eigenvalue 2 ( 1.0000):
  Row 1 difference:  4.441e-16
  Row 2 difference:  3.331e-16
  Row 3 difference:  1.110e-16
  Row 4 difference:  1.110e-16

Verifying Eigenvalue 3 ( 5.0000):
  Row 1 difference:  4.441e-16
  Row 2 difference:  2.220e-16
  Row 3 difference: -4.441e-16
  Row 4 difference:  0.000e+00

Verifying Eigenvalue 4 ( 2.0000):
  Row 1 difference:  1.110e-16
  Row 2 difference:  1.110e-16
  Row 3 difference:  2.220e-16
  Row 4 difference:  2.220e-16
https://ideone.com/CRt3oc
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!