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 (terminates when max off-diagonal element falls below this)
  6. #define MAX_ITER 100
  7.  
  8. // Function to print a matrix
  9. void print_matrix(double mat[N][N]) {
  10.     for (int i = 0; i < N; i++) {
  11.         for (int j = 0; j < N; j++) {
  12.             printf("%7.4f ", mat[i][j]);
  13.         }
  14.         printf("\n");
  15.     }
  16. }
  17.  
  18. int main() {
  19.     // Initial Matrix A Definition
  20.     double A[N][N] = {
  21.         {5.0, 4.0, 1.0, 1.0},
  22.         {4.0, 5.0, 1.0, 1.0},
  23.         {1.0, 1.0, 4.0, 2.0},
  24.         {1.0, 1.0, 2.0, 4.0}
  25.     };
  26.  
  27.     // Backup the original matrix A for verification later
  28.     double A_orig[N][N];
  29.     for (int i = 0; i < N; i++)
  30.         for (int j = 0; j < N; j++)
  31.             A_orig[i][j] = A[i][j];
  32.  
  33.     // Eigenvector matrix P (initialized as an Identity Matrix)
  34.     double P[N][N] = {0};
  35.     for (int i = 0; i < N; i++) P[i][i] = 1.0;
  36.  
  37.     int iter = 0;
  38.     printf("=== Jacobi Method Convergence Process ===\n\n");
  39.  
  40.     while (iter < MAX_ITER) {
  41.         // 1. Find the maximum off-diagonal element A[p][q]
  42.         int p = 0, q = 1;
  43.         double max_val = fabs(A[0][1]);
  44.         for (int i = 0; i < N; i++) {
  45.             for (int j = i + 1; j < N; j++) {
  46.                 if (fabs(A[i][j]) > max_val) {
  47.                     max_val = fabs(A[i][j]);
  48.                     p = i;
  49.                     q = j;
  50.                 }
  51.             }
  52.         }
  53.  
  54.         // Output convergence status to track reduction of off-diagonal elements
  55.         printf("Iteration %d: Max off-diagonal = %f (at A[%d][%d])\n", iter, max_val, p, q);
  56.  
  57.         // Convergence check
  58.         if (max_val < EPS) {
  59.             printf("\n--> Successfully converged! (Total iterations: %d)\n\n", iter);
  60.             break;
  61.         }
  62.  
  63.         // 2. Calculate the rotation angle theta (phi)
  64.         double phi, cos_t, sin_t;
  65.         if (fabs(A[p][p] - A[q][q]) < 1e-12) {
  66.             phi = acos(-1.0) / 4.0; // 45 degrees if diagonal elements are equal
  67.         } else {
  68.             phi = 0.5 * atan2(2.0 * A[p][q], A[p][p] - A[q][q]);
  69.         }
  70.         cos_t = cos(phi);
  71.         sin_t = sin(phi);
  72.  
  73.         // 3. Update Matrix A (P^T * A * P)
  74.         double Ap_old = A[p][p];
  75.         double Aq_old = A[q][q];
  76.  
  77.         A[p][p] = Ap_old * cos_t * cos_t + Aq_old * sin_t * sin_t + 2.0 * A[p][q] * sin_t * cos_t;
  78.         A[q][q] = Ap_old * sin_t * sin_t + Aq_old * cos_t * cos_t - 2.0 * A[p][q] * sin_t * cos_t;
  79.         A[p][q] = A[q][p] = 0.0; // Zero out the target element
  80.  
  81.         for (int i = 0; i < N; i++) {
  82.             if (i != p && i != q) {
  83.                 double a_ip = A[i][p];
  84.                 double a_iq = A[i][q];
  85.                 A[i][p] = A[p][i] = a_ip * cos_t + a_iq * sin_t;
  86.                 A[i][q] = A[q][i] = -a_ip * sin_t + a_iq * cos_t;
  87.             }
  88.         }
  89.  
  90.         // 4. Update Eigenvector Matrix P (P = P * P_k)
  91.         for (int i = 0; i < N; i++) {
  92.             double p_ip = P[i][p];
  93.             double p_iq = P[i][q];
  94.             P[i][p] = p_ip * cos_t + p_iq * sin_t;
  95.             P[i][q] = -p_ip * sin_t + p_iq * cos_t;
  96.         }
  97.  
  98.         iter++;
  99.     }
  100.  
  101.     // === Display Results ===
  102.     printf("=== Final Results ===\n");
  103.     printf("Diagonalized Matrix A (Diagonal Elements = Eigenvalues):\n");
  104.     print_matrix(A);
  105.     printf("\n");
  106.  
  107.     // === Verification of Ax = lambda * x ===
  108.     printf("=== Verification (Ax - lambda * x) ===\n");
  109.     printf("Note: If the values are extremely close to 0, the math holds true!\n\n");
  110.  
  111.     for (int j = 0; j < N; j++) {
  112.         double lambda = A[j][j];
  113.         printf("Eigenvalue %d: %7.4f Verification Error:\n", j + 1, lambda);
  114.  
  115.         for (int i = 0; i < N; i++) {
  116.             // Compute Ax
  117.             double Ax = 0.0;
  118.             for (int k = 0; k < N; k++) {
  119.                 Ax += A_orig[i][k] * P[k][j];
  120.             }
  121.             // Compute lambda * x
  122.             double rx = lambda * P[i][j];
  123.  
  124.             // Output the difference (Residual Error)
  125.             printf("  Row %d: %10.3e\n", i + 1, Ax - rx);
  126.         }
  127.         printf("\n");
  128.     }
  129.  
  130.     return 0;
  131. }
  132.  
Success #stdin #stdout 0s 5324KB
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! (Total iterations: 3)

=== Final Results ===
Diagonalized Matrix A (Diagonal Elements = Eigenvalues):
10.0000  0.0000  0.0000  0.0000 
 0.0000  1.0000  0.0000  0.0000 
 0.0000  0.0000  5.0000 -0.0000 
 0.0000  0.0000 -0.0000  2.0000 

=== Verification (Ax - lambda * x) ===
Note: If the values are extremely close to 0, the math holds true!

Eigenvalue 1: 10.0000 Verification Error:
  Row 1:  8.882e-16
  Row 2:  1.776e-15
  Row 3:  4.441e-16
  Row 4:  8.882e-16

Eigenvalue 2:  1.0000 Verification Error:
  Row 1:  4.441e-16
  Row 2:  3.331e-16
  Row 3:  1.110e-16
  Row 4:  1.110e-16

Eigenvalue 3:  5.0000 Verification Error:
  Row 1:  4.441e-16
  Row 2:  2.220e-16
  Row 3: -4.441e-16
  Row 4:  0.000e+00

Eigenvalue 4:  2.0000 Verification Error:
  Row 1:  1.110e-16
  Row 2:  1.110e-16
  Row 3:  2.220e-16
  Row 4:  2.220e-16
https://ideone.com/uj6gh2
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!