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. // Standard LU decomposition as shown in the lecture
  17. void lu_decomposition(double B[N][N], double L[N][N], double U[N][N]) {
  18.     for (int i = 0; i < N; i++) {
  19.         for (int j = 0; j < N; j++) {
  20.             L[i][j] = (i == j) ? 1.0 : 0.0;
  21.             U[i][j] = 0.0;
  22.         }
  23.     }
  24.     for (int i = 0; i < N; i++) {
  25.         for (int j = i; j < N; j++) {
  26.             double sum = 0.0;
  27.             for (int k = 0; k < i; k++) sum += L[i][k] * U[k][j];
  28.             U[i][j] = B[i][j] - sum;
  29.         }
  30.         for (int j = i + 1; j < N; j++) {
  31.             double sum = 0.0;
  32.             for (int k = 0; k < i; k++) sum += L[j][k] * U[k][i];
  33.             L[j][i] = (B[j][i] - sum) / U[i][i];
  34.         }
  35.     }
  36. }
  37.  
  38. void forward_substitution(double L[N][N], double x[N], double y[N]) {
  39.     for (int i = 0; i < N; i++) {
  40.         double sum = 0.0;
  41.         for (int j = 0; j < i; j++) sum += L[i][j] * y[j];
  42.         y[i] = x[i] - sum;
  43.     }
  44. }
  45.  
  46. void backward_substitution(double U[N][N], double y[N], double x[N]) {
  47.     for (int i = N - 1; i >= 0; i--) {
  48.         double sum = 0.0;
  49.         for (int j = i + 1; j < N; j++) sum += U[i][j] * x[j];
  50.         x[i] = (y[i] - sum) / U[i][i];
  51.     }
  52. }
  53.  
  54. void find_remaining_eigen(double lambda_hat) {
  55.     double B[N][N], L[N][N], U[N][N];
  56.     double y[N], x_next[N];
  57.     double x[N] = {1.0, 1.0, 1.0, 1.0}; // Initial vector x^(0)
  58.  
  59.     for (int i = 0; i < N; i++) {
  60.         for (int j = 0; j < N; j++) {
  61.             B[i][j] = A[i][j] - (i == j ? lambda_hat : 0.0);
  62.         }
  63.     }
  64.  
  65.     lu_decomposition(B, L, U);
  66.  
  67.     printf("\n==================================================\n");
  68.     printf(" Shift (lambda_hat) = %.3f\n", lambda_hat);
  69.     printf("==================================================\n");
  70.  
  71.     int iter = 0;
  72.     double diff;
  73.     double lambda = 0.0;
  74.  
  75.     do {
  76.         forward_substitution(L, x, y);      // Ly = x
  77.         backward_substitution(U, y, x_next); // Ux_next = y
  78.  
  79.         // FIX: Scale by the first element x_next[0] exactly like Slide 57
  80.         double scale_factor = x_next[0]; 
  81.  
  82.         for (int i = 0; i < N; i++) {
  83.             x_next[i] /= scale_factor;
  84.         }
  85.  
  86.         // Convergence criteria calculation
  87.         diff = 0.0;
  88.         for (int i = 0; i < N; i++) {
  89.             diff += fabs(x_next[i] - x[i]);
  90.         }
  91.  
  92.         for (int i = 0; i < N; i++) {
  93.             x[i] = x_next[i];
  94.         }
  95.  
  96.         // Slide 10 Formula: lambda = lambda_hat + 1 / scale_factor
  97.         lambda = lambda_hat + (1.0 / scale_factor);
  98.         iter++;
  99.  
  100.         printf("Iter %2d: x = [%.4f, %.4f, %.4f, %.4f], lambda = %.6f\n", 
  101.                iter, x[0], x[1], x[2], x[3], lambda);
  102.  
  103.     } while (diff > EPS && iter < MAX_ITER);
  104.  
  105.     printf("\n Solution Found:\n");
  106.     printf(" Eigenvalue (lambda) = %.6f\n", lambda);
  107.     printf(" Eigenvector (x)    = [%.4f, %.4f, %.4f, %.4f]\n", x[0], x[1], x[2], x[3]);
  108. }
  109.  
  110. int main() {
  111.     // Test a negative shift to pull out the remaining real eigenvalue cleanly
  112.     find_remaining_eigen(-6.0);
  113.     return 0;
  114. }
  115.  
Success #stdin #stdout 0s 5324KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
==================================================
 Shift (lambda_hat) = -6.000
==================================================
Iter  1: x = [1.0000, 1.8667, -2.8000, -0.2667], lambda = -5.000000
Iter  2: x = [1.0000, 1.7750, -2.9619, -0.2426], lambda = -5.548495
Iter  3: x = [1.0000, 1.7585, -2.9594, -0.2276], lambda = -5.499430
Iter  4: x = [1.0000, 1.7568, -2.9590, -0.2254], lambda = -5.490355
Iter  5: x = [1.0000, 1.7567, -2.9590, -0.2251], lambda = -5.489339
Iter  6: x = [1.0000, 1.7567, -2.9590, -0.2251], lambda = -5.489256
Iter  7: x = [1.0000, 1.7567, -2.9590, -0.2251], lambda = -5.489253
Iter  8: x = [1.0000, 1.7567, -2.9590, -0.2251], lambda = -5.489254

 Solution Found:
 Eigenvalue (lambda) = -5.489254
 Eigenvector (x)    = [1.0000, 1.7567, -2.9590, -0.2251]
https://ideone.com/ZDmCeP
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!