fork download
copy 
  1. #include <stdio.h>
  2. #include <math.h>
  3.  
  4. #define N 4
  5. #define EPS 1e-8
  6. #define MAX_ITER 100
  7.  
  8. double A[N][N] = {
  9.     {1.0, 2.0, 3.0, 5.0},
  10.     {3.0, -5.0, 1.0, 4.0},
  11.     {5.0, 9.0, 2.0, -6.0},
  12.     {1.0, 7.0, 4.0, 1.0}
  13. };
  14.  
  15. /* LU decomposition */
  16. void lu_decomposition(double B[N][N],
  17.                       double L[N][N],
  18.                       double U[N][N])
  19. {
  20.     for(int i=0;i<N;i++){
  21.         for(int j=0;j<N;j++){
  22.             L[i][j]=(i==j)?1.0:0.0;
  23.             U[i][j]=0.0;
  24.         }
  25.     }
  26.  
  27.     for(int i=0;i<N;i++){
  28.  
  29.         for(int j=i;j<N;j++){
  30.             double sum = 0.0;
  31.  
  32.             for(int k=0;k<i;k++)
  33.                 sum += L[i][k]*U[k][j];
  34.  
  35.             U[i][j] = B[i][j] - sum;
  36.         }
  37.  
  38.         for(int j=i+1;j<N;j++){
  39.             double sum = 0.0;
  40.  
  41.             for(int k=0;k<i;k++)
  42.                 sum += L[j][k]*U[k][i];
  43.  
  44.             L[j][i] = (B[j][i]-sum)/U[i][i];
  45.         }
  46.     }
  47. }
  48.  
  49. /* Forward substitution */
  50. void forward_substitution(double L[N][N],
  51.                           double x[N],
  52.                           double y[N])
  53. {
  54.     for(int i=0;i<N;i++){
  55.  
  56.         double sum = 0.0;
  57.  
  58.         for(int j=0;j<i;j++)
  59.             sum += L[i][j]*y[j];
  60.  
  61.         y[i] = x[i] - sum;
  62.     }
  63. }
  64.  
  65. /* Backward substitution */
  66. void backward_substitution(double U[N][N],
  67.                            double y[N],
  68.                            double x[N])
  69. {
  70.     for(int i=N-1;i>=0;i--){
  71.  
  72.         double sum = 0.0;
  73.  
  74.         for(int j=i+1;j<N;j++)
  75.             sum += U[i][j]*x[j];
  76.  
  77.         x[i] = (y[i]-sum)/U[i][i];
  78.     }
  79. }
  80.  
  81. void inverse_iteration(double lambda_hat)
  82. {
  83.     double B[N][N];
  84.     double L[N][N];
  85.     double U[N][N];
  86.  
  87.     double x[N] = {1.0,1.0,1.0,1.0};
  88.  
  89.     double y[N];
  90.     double x_next[N];
  91.  
  92.     for(int i=0;i<N;i++){
  93.         for(int j=0;j<N;j++){
  94.  
  95.             B[i][j] = A[i][j];
  96.  
  97.             if(i==j)
  98.                 B[i][j] -= lambda_hat;
  99.         }
  100.     }
  101.  
  102.     lu_decomposition(B,L,U);
  103.  
  104.     printf("\n=====================================\n");
  105.     printf("shift = %.6f\n", lambda_hat);
  106.     printf("=====================================\n");
  107.  
  108.     int iter=0;
  109.     double diff;
  110.     double lambda=0.0;
  111.  
  112.     do{
  113.  
  114.         forward_substitution(L,x,y);
  115.         backward_substitution(U,y,x_next);
  116.  
  117.         double max_mag = 0.0;
  118.         int idx = 0;
  119.  
  120.         for(int i=0;i<N;i++){
  121.             if(fabs(x_next[i]) > max_mag){
  122.                 max_mag = fabs(x_next[i]);
  123.                 idx = i;
  124.             }
  125.         }
  126.  
  127.         double scale = x_next[idx];
  128.  
  129.         for(int i=0;i<N;i++)
  130.             x_next[i] /= scale;
  131.  
  132.         diff = 0.0;
  133.  
  134.         for(int i=0;i<N;i++)
  135.             diff += fabs(x_next[i]-x[i]);
  136.  
  137.         for(int i=0;i<N;i++)
  138.             x[i]=x_next[i];
  139.  
  140.         lambda = lambda_hat + 1.0/scale;
  141.  
  142.         printf("Iter %2d : lambda = %.10f\n",
  143.                iter+1, lambda);
  144.  
  145.         iter++;
  146.  
  147.     }while(diff>EPS && iter<MAX_ITER);
  148.  
  149.     printf("\nFinal eigenvalue:\n");
  150.     printf("%.10f\n", lambda);
  151.  
  152.     printf("\nEigenvector:\n");
  153.     for(int i=0;i<N;i++)
  154.         printf("%.10f\n", x[i]);
  155. }
  156.  
  157. int main()
  158. {
  159.     inverse_iteration(-5.0);
  160.     inverse_iteration(-2.0);
  161.     inverse_iteration(-3.0);
  162.  
  163.     return 0;
  164. }
Success #stdin #stdout 0s 5288KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
=====================================
shift = -5.000000
=====================================
Iter  1 : lambda = -4.7713414634
Iter  2 : lambda = -5.6168969581
Iter  3 : lambda = -5.4734716409
Iter  4 : lambda = -5.4907285135
Iter  5 : lambda = -5.4892147143
Iter  6 : lambda = -5.4892322604
Iter  7 : lambda = -5.4892600355
Iter  8 : lambda = -5.4892526322
Iter  9 : lambda = -5.4892539128
Iter 10 : lambda = -5.4892537468
Iter 11 : lambda = -5.4892537615

Final eigenvalue:
-5.4892537615

Eigenvector:
-0.3379471337
-0.5936744248
1.0000000000
0.0760747116

=====================================
shift = -2.000000
=====================================
Iter  1 : lambda = -1.3181818182
Iter  2 : lambda = -3.6547277937
Iter  3 : lambda = 2.0495068652
Iter  4 : lambda = -2.7988121282
Iter  5 : lambda = 1.6517282729
Iter  6 : lambda = -2.6960303207
Iter  7 : lambda = 4.2205861066
Iter  8 : lambda = -2.4422094713
Iter  9 : lambda = -9.5197434891
Iter 10 : lambda = -1.6333113017
Iter 11 : lambda = -10.6555698873
Iter 12 : lambda = -1.6755024438
Iter 13 : lambda = -12.1277283710
Iter 14 : lambda = -1.7185880746
Iter 15 : lambda = -9.2870369763
Iter 16 : lambda = -1.6026239193
Iter 17 : lambda = -6.6554821684
Iter 18 : lambda = -1.3675767497
Iter 19 : lambda = -5.3516814496
Iter 20 : lambda = -1.1053660617
Iter 21 : lambda = -4.5539078404
Iter 22 : lambda = -0.8012649705
Iter 23 : lambda = -4.0014053883
Iter 24 : lambda = -0.4321782092
Iter 25 : lambda = -3.5851903318
Iter 26 : lambda = 0.0412921022
Iter 27 : lambda = -3.2513574975
Iter 28 : lambda = 0.6937746655
Iter 29 : lambda = -2.9698810116
Iter 30 : lambda = 1.6876328469
Iter 31 : lambda = -2.7223778830
Iter 32 : lambda = 3.4584340828
Iter 33 : lambda = -2.4965882240
Iter 34 : lambda = -9.3163573848
Iter 35 : lambda = -1.6234905829
Iter 36 : lambda = -10.3787965319
Iter 37 : lambda = -1.6661096581
Iter 38 : lambda = -11.7585350324
Iter 39 : lambda = -1.7089014901
Iter 40 : lambda = -10.2990744710
Iter 41 : lambda = -1.6523320419
Iter 42 : lambda = -7.0725227797
Iter 43 : lambda = -1.4217792894
Iter 44 : lambda = -5.5830790951
Iter 45 : lambda = -1.1666853658
Iter 46 : lambda = -4.7044450227
Iter 47 : lambda = -0.8735188617
Iter 48 : lambda = -4.1098540686
Iter 49 : lambda = -0.5215308592
Iter 50 : lambda = -3.6692287012
Iter 51 : lambda = -0.0760234789
Iter 52 : lambda = -3.3202491260
Iter 53 : lambda = 0.5270805638
Iter 54 : lambda = -3.0290145059
Iter 55 : lambda = 1.4223671019
Iter 56 : lambda = -2.7751761140
Iter 57 : lambda = 2.9512393434
Iter 58 : lambda = -2.5454166744
Iter 59 : lambda = -9.1096955282
Iter 60 : lambda = -1.6139042740
Iter 61 : lambda = -10.1205439462
Iter 62 : lambda = -1.6566627495
Iter 63 : lambda = -11.4148854083
Iter 64 : lambda = -1.6993029478
Iter 65 : lambda = -11.6252653995
Iter 66 : lambda = -1.7012871487
Iter 67 : lambda = -7.5636018501
Iter 68 : lambda = -1.4747886477
Iter 69 : lambda = -5.8424943778
Iter 70 : lambda = -1.2261469454
Iter 71 : lambda = -4.8684012342
Iter 72 : lambda = -0.9429149670
Iter 73 : lambda = -4.2256850112
Iter 74 : lambda = -0.6063629025
Iter 75 : lambda = -3.7576981757
Iter 76 : lambda = -0.1857862455
Iter 77 : lambda = -3.3919471270
Iter 78 : lambda = 0.3741786555
Iter 79 : lambda = -3.0899714074
Iter 80 : lambda = 1.1861106366
Iter 81 : lambda = -2.8291509141
Iter 82 : lambda = 2.5217206535
Iter 83 : lambda = -2.5949583158
Iter 84 : lambda = -8.9116147259
Iter 85 : lambda = -1.6042463777
Iter 86 : lambda = -9.8759765479
Iter 87 : lambda = -1.6472100477
Iter 88 : lambda = -11.0944261778
Iter 89 : lambda = -1.6897636843
Iter 90 : lambda = -12.7401135922
Iter 91 : lambda = -1.7332252150
Iter 92 : lambda = -8.1513133049
Iter 93 : lambda = -1.5267157068
Iter 94 : lambda = -6.1358141626
Iter 95 : lambda = -1.2839211629
Iter 96 : lambda = -5.0479380345
Iter 97 : lambda = -1.0097179256
Iter 98 : lambda = -4.3498722663
Iter 99 : lambda = -0.6871277837
Iter 100 : lambda = -3.8511033425

Final eigenvalue:
-3.8511033425

Eigenvector:
-0.4843350499
-0.3940645422
1.0000000000
-0.0195254644

=====================================
shift = -3.000000
=====================================
Iter  1 : lambda = -2.5000000000
Iter  2 : lambda = -13.0000000000
Iter  3 : lambda = -1.9393939394
Iter  4 : lambda = -2.1309255079
Iter  5 : lambda = -13.3537268817
Iter  6 : lambda = -2.4695420024
Iter  7 : lambda = -1.5167648190
Iter  8 : lambda = -13.1347169456
Iter  9 : lambda = -2.6074939788
Iter 10 : lambda = -1.0420120549
Iter 11 : lambda = -10.4247871184
Iter 12 : lambda = -2.5634436370
Iter 13 : lambda = -0.6340273136
Iter 14 : lambda = -7.9239470043
Iter 15 : lambda = -2.4366313254
Iter 16 : lambda = -0.2281217831
Iter 17 : lambda = -6.4075895649
Iter 18 : lambda = -2.2917211965
Iter 19 : lambda = 0.2323351010
Iter 20 : lambda = -5.4375708325
Iter 21 : lambda = -2.1395352557
Iter 22 : lambda = 0.8115209985
Iter 23 : lambda = -4.7684325556
Iter 24 : lambda = -1.9831552119
Iter 25 : lambda = 1.6103438625
Iter 26 : lambda = -4.2750213949
Iter 27 : lambda = -1.8223292390
Iter 28 : lambda = 2.8348290160
Iter 29 : lambda = -3.8908771906
Iter 30 : lambda = -1.6551479411
Iter 31 : lambda = -10.8399211563
Iter 32 : lambda = -2.4080121352
Iter 33 : lambda = -1.4786165085
Iter 34 : lambda = -11.2062409158
Iter 35 : lambda = -2.4918726629
Iter 36 : lambda = -1.2886980248
Iter 37 : lambda = -11.6393985925
Iter 38 : lambda = -2.5648689616
Iter 39 : lambda = -1.0800170628
Iter 40 : lambda = -12.1706959838
Iter 41 : lambda = -2.6300889373
Iter 42 : lambda = -0.8452200362
Iter 43 : lambda = -12.8520747750
Iter 44 : lambda = -2.6897427248
Iter 45 : lambda = -0.5738264791
Iter 46 : lambda = -8.6963646775
Iter 47 : lambda = -2.5184694241
Iter 48 : lambda = -0.2501575590
Iter 49 : lambda = -6.7692781081
Iter 50 : lambda = -2.3516381593
Iter 51 : lambda = 0.1506041682
Iter 52 : lambda = -5.6565392712
Iter 53 : lambda = -2.1894958517
Iter 54 : lambda = 0.6705226722
Iter 55 : lambda = -4.9208763133
Iter 56 : lambda = -2.0289479913
Iter 57 : lambda = 1.3874325251
Iter 58 : lambda = -4.3900781900
Iter 59 : lambda = -1.8670205451
Iter 60 : lambda = 2.4637654530
Iter 61 : lambda = -3.9824780981
Iter 62 : lambda = -1.7006359871
Iter 63 : lambda = 4.3052307814
Iter 64 : lambda = -3.6542602138
Iter 65 : lambda = -1.5263725341
Iter 66 : lambda = -11.1056829311
Iter 67 : lambda = -2.4710815475
Iter 68 : lambda = -1.3401650044
Iter 69 : lambda = -11.5185163328
Iter 70 : lambda = -2.5466127662
Iter 71 : lambda = -1.1368879866
Iter 72 : lambda = -12.0200318979
Iter 73 : lambda = -2.6136384449
Iter 74 : lambda = -0.9097226883
Iter 75 : lambda = -12.6552732627
Iter 76 : lambda = -2.6745749152
Iter 77 : lambda = -0.6491263012
Iter 78 : lambda = -9.4650097272
Iter 79 : lambda = -2.5632857591
Iter 80 : lambda = -0.3410441208
Iter 81 : lambda = -7.1648420209
Iter 82 : lambda = -2.3946828609
Iter 83 : lambda = 0.0364076405
Iter 84 : lambda = -5.9005131871
Iter 85 : lambda = -2.2316251389
Iter 86 : lambda = 0.5195998826
Iter 87 : lambda = -5.0887808000
Iter 88 : lambda = -2.0709492683
Iter 89 : lambda = 1.1741398365
Iter 90 : lambda = -4.5145966204
Iter 91 : lambda = -1.9096720573
Iter 92 : lambda = 2.1321349054
Iter 93 : lambda = -4.0800550286
Iter 94 : lambda = -1.7447652176
Iter 95 : lambda = 3.7058010795
Iter 96 : lambda = -3.7340895610
Iter 97 : lambda = -1.5729227905
Iter 98 : lambda = -11.0086550562
Iter 99 : lambda = -2.4495282760
Iter 100 : lambda = -1.3902820753

Final eigenvalue:
-1.3902820753

Eigenvector:
-0.3854644092
-0.3611375032
1.0000000000
-0.2978795833
https://ideone.com/WbJtoE
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!