fork download
copy 
  1. #include <stdio.h>
  2. #include <stdlib.h>
  3. #include <math.h>
  4. #include <time.h>
  5.  
  6. #define SIZE 4        // 行列のサイズ
  7. #define MAX_ITER 1000 // 最大反復回数
  8. #define TOL 0.0000001      // 収束判定の許容誤差
  9.  
  10.  
  11. // QR分解を行うためのヘルパー関数（グラムシュミット直交化を使用）
  12. void qr_decompose(double A[SIZE][SIZE], double Q[SIZE][SIZE], double R[SIZE][SIZE]) {
  13.     int i, j, k;
  14.     double norm;
  15.  
  16.     // 初期化
  17.     for (i = 0; i < SIZE; i++) {
  18.         for (j = 0; j < SIZE; j++) {
  19.             Q[i][j] = 0.0;
  20.             R[i][j] = 0.0;
  21.         }
  22.     }
  23.  
  24.     // グラムシュミット直交化
  25.     for (k = 0; k < SIZE; k++) {
  26.         for (i = 0; i < SIZE; i++) {
  27.             Q[i][k] = A[i][k];
  28.         }
  29.         for (j = 0; j < k; j++) {
  30.             R[j][k] = 0.0;
  31.             for (i = 0; i < SIZE; i++) {
  32.                 R[j][k] += Q[i][j] * A[i][k];
  33.             }
  34.             for (i = 0; i < SIZE; i++) {
  35.                 Q[i][k] -= R[j][k] * Q[i][j];
  36.             }
  37.         }
  38.         norm = 0.0;
  39.         for (i = 0; i < SIZE; i++) {
  40.             norm += Q[i][k] * Q[i][k];
  41.         }
  42.         R[k][k] = sqrt(norm);
  43.         for (i = 0; i < SIZE; i++) {
  44.             Q[i][k] /= R[k][k];
  45.         }
  46.     }
  47. }
  48.  
  49. // シフトなしQR法で固有値を計算
  50. void qr_algorithm(double A[SIZE][SIZE], double eigenvalues[SIZE]) {
  51.     double Q[SIZE][SIZE], R[SIZE][SIZE], A_next[SIZE][SIZE];
  52.     int i, j, k, l;
  53.  
  54.     for (k = 0; k < MAX_ITER; k++) {
  55.         qr_decompose(A, Q, R);
  56.  
  57.         // 新しい行列 A = R * Q を計算
  58.         for (i = 0; i < SIZE; i++) {
  59.             for (j = 0; j < SIZE; j++) {
  60.                 A_next[i][j] = 0.0;
  61.                 for (l = 0; l < SIZE; l++) {
  62.                     A_next[i][j] += R[i][l] * Q[l][j];
  63.                 }
  64.             }
  65.         }
  66.  
  67.         // 収束判定（非対角成分が小さいか確認）
  68.         double diff_norm = 0.0;
  69. for (i = 0; i < SIZE; i++) {
  70.     for (j = 0; j < SIZE; j++) {
  71.         if (i != j) {
  72.             diff_norm += A_next[i][j] * A_next[i][j];
  73.         }
  74.     }
  75. }
  76. if (sqrt(diff_norm) < TOL) break;
  77.  
  78.         // 更新
  79.         for (i = 0; i < SIZE; i++) {
  80.             for (j = 0; j < SIZE; j++) {
  81.                 A[i][j] = A_next[i][j];
  82.             }
  83.         }
  84.     }
  85.  
  86.     // 固有値を対角成分として取得
  87.     for (i = 0; i < SIZE; i++) {
  88.         eigenvalues[i] = A[i][i];
  89.     }
  90. }
  91.  
  92.    void solve_linear_system(double A[SIZE][SIZE], double b[SIZE], double x[SIZE]) {
  93.     double augmented[SIZE][SIZE + 1];
  94.     int i, j, k;
  95.  
  96.     // 拡大係数行列の作成
  97.     for (i = 0; i < SIZE; i++) {
  98.         for (j = 0; j < SIZE; j++) {
  99.             augmented[i][j] = A[i][j];
  100.         }
  101.         augmented[i][SIZE] = b[i]; // 最右列に b を追加
  102.     }
  103.  
  104.     // 前進消去（部分ピボット選択付き）
  105.     for (k = 0; k < SIZE; k++) {
  106.         // ピボット選択
  107.         int max_row = k;
  108.         for (i = k + 1; i < SIZE; i++) {
  109.             if (fabs(augmented[i][k]) > fabs(augmented[max_row][k])) {
  110.                 max_row = i;
  111.             }
  112.         }
  113.  
  114.         // 行交換（必要な場合）
  115.         if (max_row != k) {
  116.             for (j = 0; j <= SIZE; j++) {
  117.                 double temp = augmented[k][j];
  118.                 augmented[k][j] = augmented[max_row][j];
  119.                 augmented[max_row][j] = temp;
  120.             }
  121.         }
  122.  
  123.         // ピボット列を基準に前進消去
  124.         for (i = k + 1; i < SIZE; i++) {
  125.             double factor = augmented[i][k] / augmented[k][k];
  126.             for (j = k; j <= SIZE; j++) {
  127.                 augmented[i][j] -= factor * augmented[k][j];
  128.             }
  129.         }
  130.     }
  131.  
  132.     // 後退代入
  133.     for (i = SIZE - 1; i >= 0; i--) {
  134.         x[i] = augmented[i][SIZE]; // 拡大係数行列の最後の列が定数項
  135.         for (j = i + 1; j < SIZE; j++) {
  136.             x[i] -= augmented[i][j] * x[j];
  137.         }
  138.         x[i] /= augmented[i][i];
  139.     }
  140. }
  141.  
  142.  
  143. void inverse_iteration(double A[SIZE][SIZE], double mu, double eigenvector[SIZE]) {
  144.     double I[SIZE][SIZE], shifted_A[SIZE][SIZE], y[SIZE], x[SIZE];
  145.     int i, j, k;
  146.  
  147.     // 単位行列 I の作成
  148.     for (i = 0; i < SIZE; i++) {
  149.         for (j = 0; j < SIZE; j++) {
  150.             I[i][j] = (i == j) ? 1.0 : 0.0;
  151.         }
  152.     }
  153.  
  154.     // 初期ベクトルをランダムに設定
  155.  
  156.     for (i = 0; i < SIZE; i++) {
  157.           srand((unsigned int)time(NULL)); 
  158.         x[i] = (double)rand() / RAND_MAX; // 0.0～1.0の範囲で初期化
  159.     }
  160.  
  161.     // (A - mu * I) の計算
  162.     for (i = 0; i < SIZE; i++) {
  163.         for (j = 0; j < SIZE; j++) {
  164.             shifted_A[i][j] = A[i][j] - mu * I[i][j];
  165.         }
  166.     }
  167.  
  168.     // 反復計算
  169.     for (k = 0; k < MAX_ITER; k++) {
  170.         // (A - mu * I)y = x を解く
  171.         solve_linear_system(shifted_A, x, y);
  172.  
  173.         // 正規化
  174.         double norm = 0.0;
  175.         for (i = 0; i < SIZE; i++) {
  176.             norm += y[i] * y[i];
  177.         }
  178.         norm = sqrt(norm);
  179.         for (i = 0; i < SIZE; i++) {
  180.             y[i] /= norm; // 正規化
  181.         }
  182.  
  183.         // 収束判定
  184.         double diff = 0.0;
  185.         for (i = 0; i < SIZE; i++) {
  186.             diff += fabs(y[i] - x[i]);
  187.         }
  188.         if (diff < TOL) {
  189.             break;
  190.         }
  191.  
  192.         // x を更新
  193.         for (i = 0; i < SIZE; i++) {
  194.             x[i] = y[i];
  195.         }
  196.     }
  197.  
  198.     // 固有ベクトルを返す
  199.     for (i = 0; i < SIZE; i++) {
  200.         eigenvector[i] = y[i];
  201.     }
  202. }
  203.  
  204. // メイン関数
  205. int main() {
  206.     double A[SIZE][SIZE] = {
  207.         {16, -1, 1, 2},
  208.         {2, 12, 1, -1},
  209.         {1, 3, 24, 2},
  210.         {4, -2, 1, 20}
  211.     };
  212.     double eigenvalues[SIZE], eigenvector[SIZE];
  213.  
  214.     printf("初期行列:\n");
  215.     for (int i = 0; i < SIZE; i++) {
  216.         for (int j = 0; j < SIZE; j++) {
  217.             printf("%8.4f ", A[i][j]);
  218.         }
  219.         printf("\n");
  220.     }
  221.  
  222.     // シフトなしQR法で固有値を計算
  223.     qr_algorithm(A, eigenvalues);
  224.  
  225.     printf("\n固有値:\n");
  226.     for (int i = 0; i < SIZE; i++) {
  227.         printf("lambda[%d] = %8.4f\n", i, eigenvalues[i]);
  228.     }
  229.  
  230.  
  231.     return 0;
  232. }
  233.  
Success #stdin #stdout 0s 5288KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
初期行列:
 16.0000  -1.0000   1.0000   2.0000 
  2.0000  12.0000   1.0000  -1.0000 
  1.0000   3.0000  24.0000   2.0000 
  4.0000  -2.0000   1.0000  20.0000 

固有値:
lambda[0] =  25.0136
lambda[1] =  20.8385
lambda[2] =  14.1274
lambda[3] =  12.0204
https://ideone.com/WLszkr
language:
C (gcc 8.3)
created:
4 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!