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. // Function to perform LU decomposition exactly like your working code
  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. // Forward Substitution
  39. void forward_substitution(double L[N][N], double x[N], double y[N]) {
  40.     for (int i = 0; i < N; i++) {
  41.         double sum = 0.0;
  42.         for (int j = 0; j < i; j++) sum += L[i][j] * y[j];
  43.         y[i] = x[i] - sum;
  44.     }
  45. }
  46.  
  47. // Backward Substitution
  48. void backward_substitution(double U[N][N], double y[N], double x[N]) {
  49.     for (int i = N - 1; i >= 0; i--) {
  50.         double sum = 0.0;
  51.         for (int j = i + 1; j < N; j++) sum += U[i][j] * x[j];
  52.         x[i] = (y[i] - sum) / U[i][i];
  53.     }
  54. }
  55.  
  56. // Inverse Iteration Execution matching Matrix 2 structure exactly
  57. void find_remaining_eigen(double lambda_hat) {
  58.     double B[N][N], L[N][N], U[N][N];
  59.     double y[N], x_next[N];
  60.     double x[N] = {1.0, 1.0, 1.0, 1.0};
  61.  
  62.     for (int i = 0; i < N; i++) {
  63.         for (int j = 0; j < N; j++) {
  64.             B[i][j] = A[i][j] - (i == j ? lambda_hat : 0.0);
  65.         }
  66.     }
  67.  
  68.     lu_decomposition(B, L, U);
  69.  
  70.     printf("\n==================================================================\n");
  71.     printf(" Running Inverse Iteration for Matrix 1 (lambda_hat = %.3f)\n", lambda_hat);
  72.     printf("==================================================================\n");
  73.  
  74.     int iter = 0;
  75.     double diff;
  76.     double lambda = 0.0;
  77.  
  78.     do {
  79.         forward_substitution(L, x, y);
  80.         backward_substitution(U, y, x_next);
  81.  
  82.         // Find scaling factor (max absolute value)
  83.         double max_val = 0.0;
  84.         for (int i = 0; i < N; i++) {
  85.             if (fabs(x_next[i]) > fabs(max_val)) {
  86.                 max_val = x_next[i];
  87.             }
  88.         }
  89.  
  90.         // Normalize
  91.         for (int i = 0; i < N; i++) {
  92.             x_next[i] /= max_val;
  93.         }
  94.  
  95.         diff = 0.0;
  96.         for (int i = 0; i < N; i++) {
  97.             diff += fabs(x_next[i] - x[i]);
  98.         }
  99.  
  100.         for (int i = 0; i < N; i++) {
  101.             x[i] = x_next[i];
  102.         }
  103.  
  104.         lambda = lambda_hat + (1.0 / max_val);
  105.  
  106.         printf("Iter %2d: x = [%.4f, %.4f, %.4f, %.4f], lambda = %.6f\n", 
  107.                iter + 1, x[0], x[1], x[2], x[3], lambda);
  108.  
  109.         iter++;
  110.     } while (diff > EPS && iter < MAX_ITER);
  111.  
  112.     printf("\n Discovered Answer:\n");
  113.     printf(" Eigenvalue (lambda)   = %10.6f\n", lambda);
  114.     printf(" Eigenvector (x)        = [%.4f, %.4f, %.4f, %.4f]\n", x[0], x[1], x[2], x[3]);
  115.     printf("==================================================================\n");
  116. }
  117.  
  118. int main() {
  119.     printf("=== Matrix 1: Finding Missing Spectrum Elements via Real-Shifts ===\n");
  120.  
  121.     // We pass three distinct real shift values to uncover the remaining spots
  122.     find_remaining_eigen(-6.0); 
  123.     find_remaining_eigen(-2.0); 
  124.     find_remaining_eigen(1.0);  
  125.  
  126.     return 0;
  127. }
  128.  
Success #stdin #stdout 0s 5312KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
=== Matrix 1: Finding Missing Spectrum Elements via Real-Shifts ===

==================================================================
 Running Inverse Iteration for Matrix 1 (lambda_hat = -6.000)
==================================================================
Iter  1: x = [-0.3571, -0.6667, 1.0000, 0.0952], lambda = -6.357143
Iter  2: x = [-0.3376, -0.5993, 1.0000, 0.0819], lambda = -5.573171
Iter  3: x = [-0.3379, -0.5942, 1.0000, 0.0769], lambda = -5.499016
Iter  4: x = [-0.3380, -0.5937, 1.0000, 0.0762], lambda = -5.490283
Iter  5: x = [-0.3379, -0.5937, 1.0000, 0.0761], lambda = -5.489342
Iter  6: x = [-0.3379, -0.5937, 1.0000, 0.0761], lambda = -5.489259
Iter  7: x = [-0.3379, -0.5937, 1.0000, 0.0761], lambda = -5.489254

 Discovered Answer:
 Eigenvalue (lambda)   =  -5.489254
 Eigenvector (x)        = [-0.3379, -0.5937, 1.0000, 0.0761]
==================================================================

==================================================================
 Running Inverse Iteration for Matrix 1 (lambda_hat = -2.000)
==================================================================
Iter  1: x = [-0.2273, -0.3636, 1.0000, -0.1818], lambda = -1.318182
Iter  2: x = [-0.4549, -0.4291, 1.0000, -0.0802], lambda = -3.654728
Iter  3: x = [-0.2920, -0.2851, 1.0000, -0.6792], lambda = 2.049507
Iter  4: x = [-0.4508, -0.3721, 1.0000, -0.1340], lambda = -2.798812
Iter  5: x = [-0.2536, -0.3284, 1.0000, -0.6458], lambda = 1.651728
Iter  6: x = [-0.4360, -0.3802, 1.0000, -0.1510], lambda = -2.696030
Iter  7: x = [-0.1634, -0.2832, 1.0000, -0.9311], lambda = 4.220586
Iter  8: x = [-0.4282, -0.3749, 1.0000, -0.1787], lambda = -2.442209
Iter  9: x = [-0.0362, 0.1406, -0.6510, 1.0000], lambda = -9.519743
Iter 10: x = [-0.4195, -0.3726, 1.0000, -0.2020], lambda = -1.633311
Iter 11: x = [-0.3099, -0.0821, -0.0919, 1.0000], lambda = -10.655570
Iter 12: x = [-0.4112, -0.3697, 1.0000, -0.2255], lambda = -1.675502
Iter 13: x = [-0.6678, -0.3704, 0.6361, 1.0000], lambda = -12.127728
Iter 14: x = [-0.4028, -0.3669, 1.0000, -0.2492], lambda = -1.718588
Iter 15: x = [-0.7027, -0.4668, 1.0000, 0.5954], lambda = -9.287037
Iter 16: x = [-0.3940, -0.3640, 1.0000, -0.2739], lambda = -1.602624
Iter 17: x = [-0.5970, -0.4316, 1.0000, 0.2977], lambda = -6.655482
Iter 18: x = [-0.3846, -0.3608, 1.0000, -0.3004], lambda = -1.367577
Iter 19: x = [-0.5446, -0.4141, 1.0000, 0.1502], lambda = -5.351681
Iter 20: x = [-0.3740, -0.3573, 1.0000, -0.3301], lambda = -1.105366
Iter 21: x = [-0.5126, -0.4035, 1.0000, 0.0600], lambda = -4.553908
Iter 22: x = [-0.3618, -0.3533, 1.0000, -0.3645], lambda = -0.801265
Iter 23: x = [-0.4904, -0.3961, 1.0000, -0.0025], lambda = -4.001405
Iter 24: x = [-0.3470, -0.3483, 1.0000, -0.4063], lambda = -0.432178
Iter 25: x = [-0.4737, -0.3905, 1.0000, -0.0496], lambda = -3.585190
Iter 26: x = [-0.3279, -0.3420, 1.0000, -0.4598], lambda = 0.041292
Iter 27: x = [-0.4602, -0.3860, 1.0000, -0.0874], lambda = -3.251357
Iter 28: x = [-0.3017, -0.3333, 1.0000, -0.5336], lambda = 0.693775
Iter 29: x = [-0.4489, -0.3823, 1.0000, -0.1192], lambda = -2.969881
Iter 30: x = [-0.2618, -0.3200, 1.0000, -0.6460], lambda = 1.687633
Iter 31: x = [-0.4390, -0.3790, 1.0000, -0.1472], lambda = -2.722378
Iter 32: x = [-0.1907, -0.2963, 1.0000, -0.8463], lambda = 3.458434
Iter 33: x = [-0.4299, -0.3759, 1.0000, -0.1727], lambda = -2.496588
Iter 34: x = [0.0148, 0.1803, -0.7532, 1.0000], lambda = -9.316357
Iter 35: x = [-0.4214, -0.3731, 1.0000, -0.1968], lambda = -1.623491
Iter 36: x = [-0.2431, -0.0278, -0.2282, 1.0000], lambda = -10.378797
Iter 37: x = [-0.4130, -0.3703, 1.0000, -0.2203], lambda = -1.666110
Iter 38: x = [-0.5780, -0.2981, 0.4536, 1.0000], lambda = -11.758535
Iter 39: x = [-0.4047, -0.3675, 1.0000, -0.2438], lambda = -1.708901
Iter 40: x = [-0.7434, -0.4803, 1.0000, 0.7098], lambda = -10.299074
Iter 41: x = [-0.3960, -0.3646, 1.0000, -0.2682], lambda = -1.652332
Iter 42: x = [-0.6138, -0.4372, 1.0000, 0.3449], lambda = -7.072523
Iter 43: x = [-0.3867, -0.3616, 1.0000, -0.2943], lambda = -1.421779
Iter 44: x = [-0.5539, -0.4172, 1.0000, 0.1764], lambda = -5.583079
Iter 45: x = [-0.3765, -0.3581, 1.0000, -0.3232], lambda = -1.166685
Iter 46: x = [-0.5186, -0.4055, 1.0000, 0.0770], lambda = -4.704445
Iter 47: x = [-0.3647, -0.3542, 1.0000, -0.3563], lambda = -0.873519
Iter 48: x = [-0.4947, -0.3975, 1.0000, 0.0097], lambda = -4.109854
Iter 49: x = [-0.3506, -0.3495, 1.0000, -0.3961], lambda = -0.521531
Iter 50: x = [-0.4770, -0.3916, 1.0000, -0.0401], lambda = -3.669229
Iter 51: x = [-0.3327, -0.3436, 1.0000, -0.4465], lambda = -0.076023
Iter 52: x = [-0.4630, -0.3870, 1.0000, -0.0796], lambda = -3.320249
Iter 53: x = [-0.3084, -0.3355, 1.0000, -0.5148], lambda = 0.527081
Iter 54: x = [-0.4513, -0.3831, 1.0000, -0.1125], lambda = -3.029015
Iter 55: x = [-0.2725, -0.3235, 1.0000, -0.6160], lambda = 1.422367
Iter 56: x = [-0.4411, -0.3797, 1.0000, -0.1412], lambda = -2.775176
Iter 57: x = [-0.2110, -0.3030, 1.0000, -0.7890], lambda = 2.951239
Iter 58: x = [-0.4319, -0.3766, 1.0000, -0.1672], lambda = -2.545417
Iter 59: x = [0.0650, 0.2207, -0.8553, 1.0000], lambda = -9.109696
Iter 60: x = [-0.4232, -0.3737, 1.0000, -0.1916], lambda = -1.613904
Iter 61: x = [-0.1804, 0.0227, -0.3558, 1.0000], lambda = -10.120544
Iter 62: x = [-0.4149, -0.3709, 1.0000, -0.2151], lambda = -1.656663
Iter 63: x = [-0.4946, -0.2308, 0.2838, 1.0000], lambda = -11.414885
Iter 64: x = [-0.4065, -0.3682, 1.0000, -0.2386], lambda = -1.699303
Iter 65: x = [-0.7967, -0.4981, 1.0000, 0.8598], lambda = -11.625265
Iter 66: x = [-0.3980, -0.3653, 1.0000, -0.2627], lambda = -1.701287
Iter 67: x = [-0.6335, -0.4437, 1.0000, 0.4004], lambda = -7.563602
Iter 68: x = [-0.3889, -0.3623, 1.0000, -0.2883], lambda = -1.474789
Iter 69: x = [-0.5643, -0.4207, 1.0000, 0.2057], lambda = -5.842494
Iter 70: x = [-0.3789, -0.3589, 1.0000, -0.3164], lambda = -1.226147
Iter 71: x = [-0.5252, -0.4077, 1.0000, 0.0955], lambda = -4.868401
Iter 72: x = [-0.3675, -0.3552, 1.0000, -0.3485], lambda = -0.942915
Iter 73: x = [-0.4994, -0.3991, 1.0000, 0.0228], lambda = -4.225685
Iter 74: x = [-0.3540, -0.3506, 1.0000, -0.3866], lambda = -0.606363
Iter 75: x = [-0.4806, -0.3928, 1.0000, -0.0301], lambda = -3.757698
Iter 76: x = [-0.3371, -0.3450, 1.0000, -0.4341], lambda = -0.185786
Iter 77: x = [-0.4659, -0.3879, 1.0000, -0.0715], lambda = -3.391947
Iter 78: x = [-0.3146, -0.3375, 1.0000, -0.4975], lambda = 0.374179
Iter 79: x = [-0.4538, -0.3839, 1.0000, -0.1056], lambda = -3.089971
Iter 80: x = [-0.2820, -0.3267, 1.0000, -0.5893], lambda = 1.186111
Iter 81: x = [-0.4433, -0.3804, 1.0000, -0.1351], lambda = -2.829151
Iter 82: x = [-0.2283, -0.3088, 1.0000, -0.7404], lambda = 2.521721
Iter 83: x = [-0.4339, -0.3773, 1.0000, -0.1616], lambda = -2.594958
Iter 84: x = [0.1131, 0.2595, -0.9532, 1.0000], lambda = -8.911615
Iter 85: x = [-0.4251, -0.3743, 1.0000, -0.1862], lambda = -1.604246
Iter 86: x = [-0.1210, 0.0706, -0.4767, 1.0000], lambda = -9.875977
Iter 87: x = [-0.4167, -0.3715, 1.0000, -0.2099], lambda = -1.647210
Iter 88: x = [-0.4168, -0.1680, 0.1254, 1.0000], lambda = -11.094426
Iter 89: x = [-0.4084, -0.3688, 1.0000, -0.2333], lambda = -1.689764
Iter 90: x = [-0.8163, -0.4903, 0.9386, 1.0000], lambda = -12.740114
Iter 91: x = [-0.3999, -0.3659, 1.0000, -0.2572], lambda = -1.733225
Iter 92: x = [-0.6571, -0.4516, 1.0000, 0.4669], lambda = -8.151313
Iter 93: x = [-0.3909, -0.3630, 1.0000, -0.2824], lambda = -1.526716
Iter 94: x = [-0.5761, -0.4246, 1.0000, 0.2389], lambda = -6.135814
Iter 95: x = [-0.3812, -0.3597, 1.0000, -0.3099], lambda = -1.283921
Iter 96: x = [-0.5324, -0.4101, 1.0000, 0.1159], lambda = -5.047938
Iter 97: x = [-0.3702, -0.3560, 1.0000, -0.3409], lambda = -1.009718
Iter 98: x = [-0.5044, -0.4007, 1.0000, 0.0369], lambda = -4.349872
Iter 99: x = [-0.3572, -0.3517, 1.0000, -0.3774], lambda = -0.687128
Iter 100: x = [-0.4843, -0.3941, 1.0000, -0.0195], lambda = -3.851103

 Discovered Answer:
 Eigenvalue (lambda)   =  -3.851103
 Eigenvector (x)        = [-0.4843, -0.3941, 1.0000, -0.0195]
==================================================================

==================================================================
 Running Inverse Iteration for Matrix 1 (lambda_hat = 1.000)
==================================================================
Iter  1: x = [-nan, -nan, -nan, -nan], lambda = inf

 Discovered Answer:
 Eigenvalue (lambda)   =        inf
 Eigenvector (x)        = [-nan, -nan, -nan, -nan]
==================================================================
https://ideone.com/tUscum
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!