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. // MATRIX 1 definition
  9. double A[N][N] = { 
  10.     {1.0, 2.0, 3.0, 5.0}, 
  11.     {3.0, -5.0, 1.0, 4.0}, 
  12.     {5.0, 9.0, 2.0, -6.0}, 
  13.     {1.0, 7.0, 4.0, 1.0} 
  14. };
  15.  
  16. // LU Decomposition with Partial Pivoting to prevent numerical breakdown
  17. void lu_decomposition(double B[N][N], double L[N][N], double U[N][N], int P[N]) {
  18.     double A_temp[N][N];
  19.     for (int i = 0; i < N; i++) {
  20.         P[i] = i;
  21.         for (int j = 0; j < N; j++) {
  22.             A_temp[i][j] = B[i][j];
  23.             L[i][j] = (i == j) ? 1.0 : 0.0;
  24.             U[i][j] = 0.0;
  25.         }
  26.     }
  27.  
  28.     for (int i = 0; i < N; i++) {
  29.         int pivot_row = i;
  30.         double max_val = fabs(A_temp[i][i]);
  31.         for (int r = i + 1; r < N; r++) {
  32.             if (fabs(A_temp[r][i]) > max_val) {
  33.                 max_val = fabs(A_temp[r][i]);
  34.                 pivot_row = r;
  35.             }
  36.         }
  37.  
  38.         if (pivot_row != i) {
  39.             int t = P[i]; P[i] = P[pivot_row]; P[pivot_row] = t;
  40.             for (int k = 0; k < N; k++) {
  41.                 double tmp = A_temp[i][k];
  42.                 A_temp[i][k] = A_temp[pivot_row][k];
  43.                 A_temp[pivot_row][k] = tmp;
  44.             }
  45.         }
  46.  
  47.         for (int j = i; j < N; j++) {
  48.             double sum = 0.0;
  49.             for (int k = 0; k < i; k++) sum += L[i][k] * U[k][j];
  50.             U[i][j] = A_temp[i][j] - sum;
  51.         }
  52.         for (int j = i + 1; j < N; j++) {
  53.             double sum = 0.0;
  54.             for (int k = 0; k < i; k++) sum += L[j][k] * U[k][i];
  55.             L[j][i] = (A_temp[j][i] - sum) / U[i][i];
  56.         }
  57.     }
  58. }
  59.  
  60. void forward_substitution(double L[N][N], double x[N], double y[N], int P[N]) {
  61.     double px[N];
  62.     for (int i = 0; i < N; i++) px[i] = x[P[i]];
  63.     for (int i = 0; i < N; i++) {
  64.         double sum = 0.0;
  65.         for (int j = 0; j < i; j++) sum += L[i][j] * y[j];
  66.         y[i] = px[i] - sum;
  67.     }
  68. }
  69.  
  70. void backward_substitution(double U[N][N], double y[N], double x[N]) {
  71.     for (int i = N - 1; i >= 0; i--) {
  72.         double sum = 0.0;
  73.         for (int j = i + 1; j < N; j++) sum += U[i][j] * x[j];
  74.         x[i] = (y[i] - sum) / U[i][i];
  75.     }
  76. }
  77.  
  78. void inverse_iteration(double lambda_hat) {
  79.     double B[N][N], L[N][N], U[N][N];
  80.     double y[N], x_next[N];
  81.     int P[N];
  82.     double x[N] = {1.0, 1.0, 1.0, 1.0}; // Initial vector
  83.  
  84.     for (int i = 0; i < N; i++) {
  85.         for (int j = 0; j < N; j++) {
  86.             B[i][j] = A[i][j] - (i == j ? lambda_hat : 0.0);
  87.         }
  88.     }
  89.  
  90.     lu_decomposition(B, L, U, P);
  91.  
  92.     printf("\n--- Shifting value (lambda_hat) = %.2f ---\n", lambda_hat);
  93.  
  94.     int iter = 0;
  95.     double diff;
  96.     double lambda = 0.0;
  97.  
  98.     do {
  99.         forward_substitution(L, x, y, P);
  100.         backward_substitution(U, y, x_next);
  101.  
  102.         // Find the element with maximum absolute value (keeping its sign)
  103.         double max_val = x_next[0];
  104.         for (int i = 1; i < N; i++) {
  105.             if (fabs(x_next[i]) > fabs(max_val)) {
  106.                 max_val = x_next[i];
  107.             }
  108.         }
  109.  
  110.         // Normalize the vector (Slide 9/10 style)
  111.         for (int i = 0; i < N; i++) {
  112.             x_next[i] /= max_val;
  113.         }
  114.  
  115.         // Check convergence
  116.         diff = 0.0;
  117.         for (int i = 0; i < N; i++) {
  118.             diff += fabs(x_next[i] - x[i]);
  119.         }
  120.  
  121.         for (int i = 0; i < N; i++) {
  122.             x[i] = x_next[i];
  123.         }
  124.  
  125.         // Slide 10 Formula: lambda = lambda_hat + 1 / max_val
  126.         lambda = lambda_hat + (1.0 / max_val);
  127.         iter++;
  128.  
  129.         printf("Iter %2d: Lambda = %10.6f\n", iter, lambda);
  130.  
  131.     } while (diff > EPS && iter < MAX_ITER);
  132.  
  133.     printf("Result -> Eigenvalue: %10.6f\n", lambda);
  134.     printf("Result -> Eigenvector: [%.4f, %.4f, %.4f, %.4f]\n", x[0], x[1], x[2], x[3]);
  135. }
  136.  
  137. int main() {
  138.     printf("=== Matrix 1 Inverse Iteration Solutions ===\n");
  139.  
  140.     // Try running it with a couple of different real shifting guesses
  141.     inverse_iteration(8.5);  // Targets the dominant real eigenvalue
  142.     inverse_iteration(-6.0); // Targets the negative real eigenvalue
  143.  
  144.     return 0;
  145. }
Success #stdin #stdout 0s 5312KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
=== Matrix 1 Inverse Iteration Solutions ===

--- Shifting value (lambda_hat) = 8.50 ---
Iter  1: Lambda =   8.648460
Iter  2: Lambda =   8.697207
Iter  3: Lambda =   8.696197
Iter  4: Lambda =   8.696215
Iter  5: Lambda =   8.696215
Result -> Eigenvalue:   8.696215
Result -> Eigenvector: [1.0000, 0.5414, 0.6254, 0.9474]

--- Shifting value (lambda_hat) = -6.00 ---
Iter  1: Lambda =   4.984375
Iter  2: Lambda =   4.949623
Iter  3: Lambda =  -8.965511
Iter  4: Lambda =  -8.275506
Iter  5: Lambda =  -3.337073
Iter  6: Lambda =  -2.988731
Iter  7: Lambda =  -2.837259
Iter  8: Lambda =  -2.793058
Iter  9: Lambda =  -2.792802
Iter 10: Lambda =  -2.802882
Iter 11: Lambda =  -2.810887
Iter 12: Lambda =  -2.814872
Iter 13: Lambda =  -2.816180
Iter 14: Lambda =  -2.816301
Iter 15: Lambda =  -2.816101
Iter 16: Lambda =  -2.815916
Iter 17: Lambda =  -2.815816
Iter 18: Lambda =  -2.815780
Iter 19: Lambda =  -2.815775
Iter 20: Lambda =  -2.815778
Iter 21: Lambda =  -2.815782
Iter 22: Lambda =  -2.815785
Result -> Eigenvalue:  -2.815785
Result -> Eigenvector: [-0.8610, -0.1658, 1.0000, 0.1234]
https://ideone.com/KCFDwY
language:
C (gcc 8.3)
created:
3 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!