fork(1) 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分解 */
  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. /* 前進代入 */
  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. /* 後退代入 */
  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(" 行列1 に対する逆反復法を実行中 (近似固有値 lambda_hat = %.3f %+.3fi)\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.         /* 絶対値が最大の要素を探索 */
  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.         // 収束状況の出力レイアウト
  152.         printf("反復 %2d: x = [(%.4f%+.4fi), (%.4f%+.4fi), (%.4f%+.4fi), (%.4f%+.4fi)], 固有値 = %.6f%+.6fi\n", 
  153.                iter + 1, 
  154.                creal(x[0]), cimag(x[0]),
  155.                creal(x[1]), cimag(x[1]),
  156.                creal(x[2]), cimag(x[2]),
  157.                creal(x[3]), cimag(x[3]),
  158.                creal(lambda), cimag(lambda));
  159.  
  160.         iter++;
  161.  
  162.     }while(diff>EPS && iter<MAX_ITER);
  163.  
  164.     // 最終結果出力ブロック
  165.     printf("\n 計算結果:\n");
  166.     printf(" 収束固有値 (lambda)   = %10.6f %+.6fi\n", creal(lambda), cimag(lambda));
  167.     printf(" 収束固有ベクトル (x)  = [(%.4f%+.4fi), (%.4f%+.4fi), (%.4f%+.4fi), (%.4f%+.4fi)]\n", 
  168.            creal(x[0]), cimag(x[0]),
  169.            creal(x[1]), cimag(x[1]),
  170.            creal(x[2]), cimag(x[2]),
  171.            creal(x[3]), cimag(x[3]));
  172.     printf("==================================================================\n");
  173.  
  174.     // 固有値・固有ベクトルであることの確認（残差チェック）
  175.     printf("\n 固有対の検証チェック (A*x - lambda*x):\n");
  176.     for (int i = 0; i < N; i++) {
  177.         double complex ax_i = 0.0 + 0.0*I;
  178.         for (int j = 0; j < N; j++) {
  179.             ax_i += A[i][j] * x[j];
  180.         }
  181.         double complex lx_i = lambda * x[i];
  182.         double complex residue = ax_i - lx_i;
  183.         printf(" 第%d行: A*x = (%9.5f%+.5fi), lambda*x = (%9.5f%+.5fi), 残差絶対値 = %e\n", 
  184.                i + 1, 
  185.                creal(ax_i), cimag(ax_i), 
  186.                creal(lx_i), cimag(lx_i), 
  187.                cabs(residue));
  188.     }
  189.     printf("==================================================================\n");
  190. }
  191.  
  192. int main()
  193. {
  194.     printf("=== 行列1: 複素シフトによる残りの固有値スペックトルの探索 ===\n");
  195.  
  196.     // 近似固有値の設定（必要に応じて自由に変更可能）
  197.     inverse_iteration(-5.4 + 0.0*I);
  198.     inverse_iteration(-2.1 + 1.7*I);
  199.     inverse_iteration(-2.1 - 1.7*I);
  200.  
  201.     return 0;
  202. }
  203.  
Success #stdin #stdout 0.01s 5324KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
=== 行列1: 複素シフトによる残りの固有値スペックトルの探索 ===

==================================================================
 行列1 に対する逆反復法を実行中 (近似固有値 lambda_hat = -5.400 +0.000i)
==================================================================
反復  1: x = [(-0.3354+0.0000i), (-0.5824-0.0000i), (1.0000+0.0000i), (0.0721+0.0000i)], 固有値 = -5.350763+0.000000i
反復  2: x = [(-0.3379-0.0000i), (-0.5939-0.0000i), (1.0000-0.0000i), (0.0763-0.0000i)], 固有値 = -5.492612+0.000000i
反復  3: x = [(-0.3379+0.0000i), (-0.5937-0.0000i), (1.0000-0.0000i), (0.0761-0.0000i)], 固有値 = -5.489178+0.000000i
反復  4: x = [(-0.3379-0.0000i), (-0.5937-0.0000i), (1.0000-0.0000i), (0.0761-0.0000i)], 固有値 = -5.489255+0.000000i
反復  5: x = [(-0.3379+0.0000i), (-0.5937-0.0000i), (1.0000-0.0000i), (0.0761-0.0000i)], 固有値 = -5.489254+0.000000i
反復  6: x = [(-0.3379-0.0000i), (-0.5937-0.0000i), (1.0000-0.0000i), (0.0761-0.0000i)], 固有値 = -5.489254+0.000000i

 計算結果:
 収束固有値 (lambda)   =  -5.489254 +0.000000i
 収束固有ベクトル (x)  = [(-0.3379-0.0000i), (-0.5937-0.0000i), (1.0000-0.0000i), (0.0761-0.0000i)]
==================================================================

 固有対の検証チェック (A*x - lambda*x):
 第1行: A*x = (  1.85508+0.00000i), lambda*x = (  1.85508+0.00000i), 残差絶対値 = 8.030443e-11
 第2行: A*x = (  3.25883+0.00000i), lambda*x = (  3.25883+0.00000i), 残差絶対値 = 2.493117e-12
 第3行: A*x = ( -5.48925+0.00000i), lambda*x = ( -5.48925+0.00000i), 残差絶対値 = 0.000000e+00
 第4行: A*x = ( -0.41759+0.00000i), lambda*x = ( -0.41759+0.00000i), 残差絶対値 = 1.728266e-10
==================================================================

==================================================================
 行列1 に対する逆反復法を実行中 (近似固有値 lambda_hat = -2.100 +1.700i)
==================================================================
反復  1: x = [(-0.4181+0.0652i), (-0.3703+0.0217i), (1.0000-0.0000i), (-0.2179-0.1941i)], 固有値 = -2.115342+1.685657i
反復  2: x = [(-0.4141+0.0673i), (-0.3707+0.0224i), (1.0000+0.0000i), (-0.2172-0.1894i)], 固有値 = -2.103743+1.674409i
反復  3: x = [(-0.4141+0.0673i), (-0.3707+0.0224i), (1.0000+0.0000i), (-0.2172-0.1894i)], 固有値 = -2.103484+1.674116i
反復  4: x = [(-0.4141+0.0673i), (-0.3707+0.0224i), (1.0000+0.0000i), (-0.2172-0.1894i)], 固有値 = -2.103481+1.674115i
反復  5: x = [(-0.4141+0.0673i), (-0.3707+0.0224i), (1.0000+0.0000i), (-0.2172-0.1894i)], 固有値 = -2.103481+1.674115i

 計算結果:
 収束固有値 (lambda)   =  -2.103481 +1.674115i
 収束固有ベクトル (x)  = [(-0.4141+0.0673i), (-0.3707+0.0224i), (1.0000+0.0000i), (-0.2172-0.1894i)]
==================================================================

 固有対の検証チェック (A*x - lambda*x):
 第1行: A*x = (  0.75849-0.83477i), lambda*x = (  0.75849-0.83477i), 残差絶対値 = 4.659298e-12
 第2行: A*x = (  0.74222-0.66768i), lambda*x = (  0.74222-0.66768i), 残差絶対値 = 2.149513e-11
 第3行: A*x = ( -2.10348+1.67411i), lambda*x = ( -2.10348+1.67411i), 残差絶対値 = 8.005932e-16
 第4行: A*x = (  0.77391+0.03470i), lambda*x = (  0.77391+0.03470i), 残差絶対値 = 5.726044e-11
==================================================================

==================================================================
 行列1 に対する逆反復法を実行中 (近似固有値 lambda_hat = -2.100 -1.700i)
==================================================================
反復  1: x = [(-0.4181-0.0652i), (-0.3703-0.0217i), (1.0000-0.0000i), (-0.2179+0.1941i)], 固有値 = -2.115342-1.685657i
反復  2: x = [(-0.4141-0.0673i), (-0.3707-0.0224i), (1.0000-0.0000i), (-0.2172+0.1894i)], 固有値 = -2.103743-1.674409i
反復  3: x = [(-0.4141-0.0673i), (-0.3707-0.0224i), (1.0000-0.0000i), (-0.2172+0.1894i)], 固有値 = -2.103484-1.674116i
反復  4: x = [(-0.4141-0.0673i), (-0.3707-0.0224i), (1.0000-0.0000i), (-0.2172+0.1894i)], 固有値 = -2.103481-1.674115i
反復  5: x = [(-0.4141-0.0673i), (-0.3707-0.0224i), (1.0000-0.0000i), (-0.2172+0.1894i)], 固有値 = -2.103481-1.674115i

 計算結果:
 収束固有値 (lambda)   =  -2.103481 -1.674115i
 収束固有ベクトル (x)  = [(-0.4141-0.0673i), (-0.3707-0.0224i), (1.0000-0.0000i), (-0.2172+0.1894i)]
==================================================================

 固有対の検証チェック (A*x - lambda*x):
 第1行: A*x = (  0.75849+0.83477i), lambda*x = (  0.75849+0.83477i), 残差絶対値 = 4.659298e-12
 第2行: A*x = (  0.74222+0.66768i), lambda*x = (  0.74222+0.66768i), 残差絶対値 = 2.149513e-11
 第3行: A*x = ( -2.10348-1.67411i), lambda*x = ( -2.10348-1.67411i), 残差絶対値 = 8.005932e-16
 第4行: A*x = (  0.77391-0.03470i), lambda*x = (  0.77391-0.03470i), 残差絶対値 = 5.726044e-11
==================================================================
https://ideone.com/Tfwzr4
language:
C (gcc 8.3)
created:
11 hours 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!