fork download
copy 
  1. #define _CRT_SECURE_NO_WARNINGS // Visual Studio のときは入れてください
  2. #include <math.h>
  3. #include <stdio.h>
  4.  
  5. #ifndef SIZE
  6. #define SIZE (16)
  7. #endif
  8.  
  9. double matA[SIZE][SIZE], matE[SIZE][SIZE], vecX[SIZE];
  10.  
  11. void dump(int n) {
  12.     int j, k;
  13.     for (k = 0; k < n; k++) {
  14.         printf("\n");
  15.         for (j = 0; j < n; j++)
  16.             printf(" %10.3le", matA[k][j]);
  17.         printf(" %10.3le ", vecX[k]);
  18.         for (j = 0; j < n; j++)
  19.             printf(" %10.3le", matE[k][j]);
  20.     }
  21.     printf("\n");
  22. }
  23.  
  24. int main(int argc, char **argv) {
  25.     double C, D, ratio, big;
  26.     int n, i, j, k, r, pivot_row;
  27.     char copyright[] = "Written by ********: ****/**/**\n";
  28.     int ch, verbose;
  29.  
  30. // getopt は Visual Studio に含まれていない。自作するか探すこと
  31.     while ((ch = getopt(argc, argv, "vh")) != -1) {
  32.         if (ch == 'v')
  33.             verbose = 1;
  34.         if (ch == 'h') {
  35.             printf("Options:\n");
  36.             printf("\t-v: dumps intermediate results.\n");
  37.             printf("\t-h: shows help message, as you see now.\n");
  38.             printf("\n");
  39.             return 1;
  40.         }
  41.     }
  42. // ターミナルから指定できないので
  43.     verbose = 1;
  44.  
  45.     /* データ入力 */
  46.     scanf("%d", &n);
  47.     if (n > SIZE) {
  48.         printf("N should be less than or equal to %d.\n", SIZE);
  49.         return 1;
  50.     }
  51.  
  52.     for (i = 0; i < n; i++) {
  53.         for (j = 0; j < n; j++) {
  54.             scanf("%lf", &(matA[i][j]));
  55.         }
  56.         scanf("%lf", &(vecX[i]));
  57.     }
  58.  
  59.     printf("Inputs:\n");
  60.  
  61.     for (i = 0; i < n; i++) {
  62.         for (j = 0; j < n; j++) {
  63.             printf(" %12.5le", matA[i][j]);
  64.         }
  65.         printf(" | %12.5le\n", vecX[i]);
  66.     }
  67.  
  68.     /* 単位行列の設定 */
  69.     for (i = 0; i < n; i++) {
  70.         for (j = 0; j < n; j++) {
  71.             matE[i][j] = 0.0;
  72.         }
  73.         matE[i][i] = 1.0;
  74.     }
  75.  
  76.     /* Solve */
  77.     if (verbose)
  78.         dump(n);
  79.  
  80.     for (i = 0; i < n; i++) {
  81.  
  82.         printf("############ i=%d", i);
  83.         if (verbose)
  84.             dump(n);
  85.         printf("-------\n");
  86.  
  87. // C はここで表示すべきでない。ピボット選択による並べ替えが終わった後。
  88. //      C = matA[i][i];
  89. //      printf("C=%.2f\n", C);
  90.  
  91.         /* ピボット選択による並べ替え(始) */
  92.  
  93.         big = 0.0;
  94.  
  95.         /* i列で最大値を検索 */
  96.  
  97. //      if (C < 1.0e-04 && -1.0e-04 < C) { // 削除。ここで判定するとピボット選択が始まらない
  98.             for (k = i; k < n; k++) { // k = 0 を k = i に訂正。既にピボット選択が終わった行は残す
  99.                 if (fabs(matA[k][i]) > big) {
  100.                     big = fabs(matA[k][i]);
  101.                     pivot_row = k;
  102.                 }
  103.             }
  104.             /* 検索した最大値が0または小さい値であったら計算できないと返す */
  105.  
  106.             if (big < 1.0e-04 && -1.0e-04 < big) {
  107.                 printf("Invalid matrix");
  108.                 return 0;
  109.             }
  110.  
  111.             /*最大値を含む行と入れ替え*/
  112.  
  113.             if (i != pivot_row) {
  114.                 for (k = 0; k < n; k++) {
  115.                     D = matA[i][k];
  116.                     matA[i][k] = matA[pivot_row][k];
  117.                     matA[pivot_row][k] = D;
  118.  
  119.                     D = matE[i][k]; // matE の入れ替えが抜けている
  120.                     matE[i][k] = matE[pivot_row][k];
  121.                     matE[pivot_row][k] = D;
  122.                 }
  123.                 D = vecX[i]; // [k] を [i] に訂正
  124.                 vecX[i] = vecX[pivot_row];
  125.                 vecX[pivot_row] = D;
  126.             }
  127. //      } // 削除
  128.         /*ピボット選択による並べ替え(終)*/
  129.  
  130.         C = matA[i][i];        // ここに移動
  131.         printf("C=%.2f\n", C); // ここに移動
  132.  
  133.         /********************************************/
  134.         /** i 列目前半の処理 **/
  135.         /********************************************/
  136.  
  137.         for (k = 0; k < n; k++) {
  138.             matA[i][k] = matA[i][k] / C;
  139.             matE[i][k] = matE[i][k] / C;
  140.         }
  141.         vecX[i] = vecX[i] / C;
  142.  
  143.         if (verbose)
  144.             dump(n);
  145.  
  146.         for (k = 0; k < n; k++) { /* 行のループ */
  147.             if (i != k) {
  148.                 D = matA[k][i];
  149.  
  150.                 /********************************************/
  151.                 /** i 列目後半の処理 **/
  152.                 /********************************************/
  153.                 for (j = 0; j < n; j++) { /* 列のループ */
  154.                     matA[k][j] = matA[k][j] - D * matA[i][j];
  155.                     matE[k][j] = matE[k][j] - D * matE[i][j];
  156.                 }
  157.                 vecX[k] = vecX[k] - D * vecX[i];
  158.             }
  159.         }
  160.  
  161.         if (verbose)
  162.             dump(n);
  163.     }
  164.  
  165.     /* 計算結果出力 */
  166.     printf("\nResults:\n");
  167.     for (i = 0; i < n; i++)
  168.         printf(" %12.5le", vecX[i]);
  169.     printf("\nInverse:\n");
  170.     for (i = 0; i < n; i++) {
  171.         for (j = 0; j < n; j++)
  172.             printf(" %12.5le", matE[i][j]);
  173.         printf("\n");
  174.     }
  175.     printf("\n");
  176. }
  177.  
Success #stdin #stdout 0s 4408KB
comments (?)
stdin
copy 
2

1 1 5
2 4 16
stdout
copy 
Inputs:
  1.00000e+00  1.00000e+00 |  5.00000e+00
  2.00000e+00  4.00000e+00 |  1.60000e+01

  1.000e+00  1.000e+00  5.000e+00   1.000e+00  0.000e+00
  2.000e+00  4.000e+00  1.600e+01   0.000e+00  1.000e+00
############ i=0
  1.000e+00  1.000e+00  5.000e+00   1.000e+00  0.000e+00
  2.000e+00  4.000e+00  1.600e+01   0.000e+00  1.000e+00
-------
C=2.00

  1.000e+00  2.000e+00  8.000e+00   0.000e+00  5.000e-01
  1.000e+00  1.000e+00  5.000e+00   1.000e+00  0.000e+00

  1.000e+00  2.000e+00  8.000e+00   0.000e+00  5.000e-01
  0.000e+00 -1.000e+00 -3.000e+00   1.000e+00 -5.000e-01
############ i=1
  1.000e+00  2.000e+00  8.000e+00   0.000e+00  5.000e-01
  0.000e+00 -1.000e+00 -3.000e+00   1.000e+00 -5.000e-01
-------
C=-1.00

  1.000e+00  2.000e+00  8.000e+00   0.000e+00  5.000e-01
 -0.000e+00  1.000e+00  3.000e+00  -1.000e+00  5.000e-01

  1.000e+00  0.000e+00  2.000e+00   2.000e+00 -5.000e-01
 -0.000e+00  1.000e+00  3.000e+00  -1.000e+00  5.000e-01

Results:
  2.00000e+00  3.00000e+00
Inverse:
  2.00000e+00 -5.00000e-01
 -1.00000e+00  5.00000e-01
https://ideone.com/xB1hyr
language:
C (gcc 8.3)
created:
4 years 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!