fork download
copy 
  1. #include <stdio.h>
  2. #include <math.h>
  3.  
  4. #define N 4
  5. #define EPSILON 1e-6
  6. #define MAX_ITER 1000
  7.  
  8. void matrixVectorMultiply(double mat[N][N], double vec[N], double result[N]) {
  9.     for (int i = 0; i < N; i++) {
  10.         result[i] = 0;
  11.         for (int j = 0; j < N; j++) {
  12.             result[i] += mat[i][j] * vec[j];
  13.         }
  14.     }
  15. }
  16.  
  17. double vectorNorm(double vec[N]) {
  18.     double sum = 0;
  19.     for (int i = 0; i < N; i++) {
  20.         sum += vec[i] * vec[i];
  21.     }
  22.     return sqrt(sum);
  23. }
  24.  
  25. void normalizeVector(double vec[N]) {
  26.     double norm = vectorNorm(vec);
  27.     for (int i = 0; i < N; i++) {
  28.         vec[i] /= norm;
  29.     }
  30. }
  31.  
  32. double powerMethod(double mat[N][N], double eigenvector[N]) {
  33.     double vec[N] = {1, 1, 1, 1}; // 初期ベクトル
  34.     double lambda = 0;
  35.     double lambda_old;
  36.  
  37.     for (int iter = 0; iter < MAX_ITER; iter++) {
  38.         double temp[N];
  39.         matrixVectorMultiply(mat, vec, temp);
  40.  
  41.         lambda_old = lambda;
  42.         lambda = vectorNorm(temp);
  43.  
  44.         for (int i = 0; i < N; i++) {
  45.             vec[i] = temp[i];
  46.         }
  47.         normalizeVector(vec);
  48.  
  49.         if (fabs(lambda - lambda_old) < EPSILON) {
  50.             break;
  51.         }
  52.     }
  53.  
  54.     for (int i = 0; i < N; i++) {
  55.         eigenvector[i] = vec[i];
  56.     }
  57.  
  58.     return lambda;
  59. }
  60.  
  61. int main() {
  62.     double mat[N][N] = {
  63.         {16, -1, 1, 2},
  64.         {2, 12, 1, -1},
  65.         {1, 3, 24, 2},
  66.         {4, -2, 1, 20}
  67.     };
  68.  
  69.     double eigenvector[N];
  70.     double eigenvalue = powerMethod(mat, eigenvector);
  71.  
  72.     printf("最大固有値: %lf\n", eigenvalue);
  73.     printf("対応する固有ベクトル: [");
  74.     for (int i = 0; i < N; i++) {
  75.         printf("%lf", eigenvector[i]);
  76.         if (i < N - 1) {
  77.             printf(", ");
  78.         }
  79.     }
  80.     printf("]\n");
  81.  
  82.     return 0;
  83. }
  84.  
  85.  
Success #stdin #stdout 0s 5292KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
最大固有値: 25.013627
対応する固有ベクトル: [0.159495, 0.074987, 0.942110, 0.285249]
https://ideone.com/y6wXX3
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!