fork download
copy 
  1. #include <stdio.h>
  2. #include <math.h>
  3. #include <complex.h>
  4.  
  5. #define N 4
  6. #define EPS 1e-8
  7. #define MAX_ITER 100
  8.  
  9. double complex 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 */
  17. void lu_decomposition(double complex B[N][N],
  18.                       double complex L[N][N],
  19.                       double complex U[N][N])
  20. {
  21.     for(int i=0;i<N;i++){
  22.         for(int j=0;j<N;j++){
  23.             L[i][j]=(i==j)?1.0+0.0*I:0.0+0.0*I;
  24.             U[i][j]=0.0+0.0*I;
  25.         }
  26.     }
  27.  
  28.     for(int i=0;i<N;i++){
  29.  
  30.         for(int j=i;j<N;j++){
  31.             double complex sum = 0.0 + 0.0*I;
  32.  
  33.             for(int k=0;k<i;k++)
  34.                 sum += L[i][k]*U[k][j];
  35.  
  36.             U[i][j] = B[i][j] - sum;
  37.         }
  38.  
  39.         for(int j=i+1;j<N;j++){
  40.             double complex sum = 0.0 + 0.0*I;
  41.  
  42.             for(int k=0;k<i;k++)
  43.                 sum += L[j][k]*U[k][i];
  44.  
  45.             L[j][i] = (B[j][i]-sum)/U[i][i];
  46.         }
  47.     }
  48. }
  49.  
  50. /* Forward substitution */
  51. void forward_substitution(double complex L[N][N],
  52.                           double complex x[N],
  53.                           double complex y[N])
  54. {
  55.     for(int i=0;i<N;i++){
  56.  
  57.         double complex sum = 0.0 + 0.0*I;
  58.  
  59.         for(int j=0;j<i;j++)
  60.             sum += L[i][j]*y[j];
  61.  
  62.         y[i] = x[i] - sum;
  63.     }
  64. }
  65.  
  66. /* Backward substitution */
  67. void backward_substitution(double complex U[N][N],
  68.                            double complex y[N],
  69.                            double complex x[N])
  70. {
  71.     for(int i=N-1;i>=0;i--){
  72.  
  73.         double complex sum = 0.0 + 0.0*I;
  74.  
  75.         for(int j=i+1;j<N;j++)
  76.             sum += U[i][j]*x[j];
  77.  
  78.         x[i] = (y[i]-sum)/U[i][i];
  79.     }
  80. }
  81.  
  82. void inverse_iteration(double complex lambda_hat)
  83. {
  84.     double complex B[N][N];
  85.     double complex L[N][N];
  86.     double complex U[N][N];
  87.  
  88.     double complex x[N] =
  89.     {
  90.         1.0+0.0*I,
  91.         1.0+0.0*I,
  92.         1.0+0.0*I,
  93.         1.0+0.0*I
  94.     };
  95.  
  96.     double complex y[N];
  97.     double complex x_next[N];
  98.  
  99.     for(int i=0;i<N;i++){
  100.         for(int j=0;j<N;j++){
  101.  
  102.             B[i][j] = A[i][j];
  103.  
  104.             if(i==j)
  105.                 B[i][j] -= lambda_hat;
  106.         }
  107.     }
  108.  
  109.     lu_decomposition(B,L,U);
  110.  
  111.     printf("\n=====================================\n");
  112.     printf("shift = %.6f %+ .6fi\n",
  113.            creal(lambda_hat), cimag(lambda_hat));
  114.     printf("=====================================\n");
  115.  
  116.     int iter=0;
  117.     double diff;
  118.     double complex lambda = 0.0 + 0.0*I;
  119.  
  120.     do{
  121.  
  122.         forward_substitution(L,x,y);
  123.         backward_substitution(U,y,x_next);
  124.  
  125.         /* largest magnitude component */
  126.         double max_mag = 0.0;
  127.         int idx = 0;
  128.  
  129.         for(int i=0;i<N;i++){
  130.             if(cabs(x_next[i]) > max_mag){
  131.                 max_mag = cabs(x_next[i]);
  132.                 idx = i;
  133.             }
  134.         }
  135.  
  136.         double complex scale = x_next[idx];
  137.  
  138.         for(int i=0;i<N;i++)
  139.             x_next[i] /= scale;
  140.  
  141.         diff = 0.0;
  142.  
  143.         for(int i=0;i<N;i++)
  144.             diff += cabs(x_next[i]-x[i]);
  145.  
  146.         for(int i=0;i<N;i++)
  147.             x[i]=x_next[i];
  148.  
  149.         lambda = lambda_hat + 1.0/scale;
  150.  
  151.         printf("Iter %2d : lambda = %.8f %+ .8fi\n",
  152.                iter+1,
  153.                creal(lambda),
  154.                cimag(lambda));
  155.  
  156.         iter++;
  157.  
  158.     }while(diff>EPS && iter<MAX_ITER);
  159.  
  160.     printf("\nFinal eigenvalue:\n");
  161.     printf("%.10f %+ .10fi\n",
  162.            creal(lambda),
  163.            cimag(lambda));
  164.  
  165.     printf("\nEigenvector:\n");
  166.  
  167.     for(int i=0;i<N;i++){
  168.         printf("(%.6f %+ .6fi)\n",
  169.                creal(x[i]),
  170.                cimag(x[i]));
  171.     }
  172. }
  173.  
  174. int main()
  175. {
  176.     inverse_iteration(-5.4 + 0.0*I);
  177.  
  178.     inverse_iteration(-2.1 + 1.7*I);
  179.  
  180.     inverse_iteration(-2.1 - 1.7*I);
  181.  
  182.     return 0;
  183. }
  184.  
Success #stdin #stdout 0s 5316KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
=====================================
shift = -5.400000 +0.000000i
=====================================
Iter  1 : lambda = -5.35076268 +0.00000000i
Iter  2 : lambda = -5.49261247 +0.00000000i
Iter  3 : lambda = -5.48917817 +0.00000000i
Iter  4 : lambda = -5.48925515 +0.00000000i
Iter  5 : lambda = -5.48925374 +0.00000000i
Iter  6 : lambda = -5.48925376 +0.00000000i

Final eigenvalue:
-5.4892537610 +0.0000000000i

Eigenvector:
(-0.337947 -0.000000i)
(-0.593674 -0.000000i)
(1.000000 -0.000000i)
(0.076075 -0.000000i)

=====================================
shift = -2.100000 +1.700000i
=====================================
Iter  1 : lambda = -2.11534171 +1.68565706i
Iter  2 : lambda = -2.10374252 +1.67440919i
Iter  3 : lambda = -2.10348353 +1.67411596i
Iter  4 : lambda = -2.10348057 +1.67411478i
Iter  5 : lambda = -2.10348055 +1.67411478i

Final eigenvalue:
-2.1034805502 +1.6741147833i

Eigenvector:
(-0.414119 +0.067262i)
(-0.370680 +0.022401i)
(1.000000 +0.000000i)
(-0.217207 -0.189366i)

=====================================
shift = -2.100000 -1.700000i
=====================================
Iter  1 : lambda = -2.11534171 -1.68565706i
Iter  2 : lambda = -2.10374252 -1.67440919i
Iter  3 : lambda = -2.10348353 -1.67411596i
Iter  4 : lambda = -2.10348057 -1.67411478i
Iter  5 : lambda = -2.10348055 -1.67411478i

Final eigenvalue:
-2.1034805502 -1.6741147833i

Eigenvector:
(-0.414119 -0.067262i)
(-0.370680 -0.022401i)
(1.000000 -0.000000i)
(-0.217207 +0.189366i)
https://ideone.com/gB3Cmr
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!