ここから本文です

プログラミングが得意な方、教えてください。 以下のプログラムは、N×N行列をGaus...

アバター

ID非公開さん

2019/7/2819:00:04

プログラミングが得意な方、教えてください。
以下のプログラムは、N×N行列をGauss-Jordan法によって簡約化して、N個の未知数をもつ連立方程式きをとかプログラムです。
ここで、ピボット選択

による並び替えのプログラムの部分が正常に動作せず、困っています。なぜできないか、教えて頂けると幸いです

#include <stdio.h>
#include <math.h>

#ifndef SIZE
#define SIZE (16)
#endif


double matA[SIZE][SIZE], matE[SIZE][SIZE], vecX[SIZE];

void dump ( int n )
{
int j, k;
for ( k=0; k<n; k++ ){
printf( "\n" );
for ( j=0; j<n; j++ ) printf( " %10.3le", matA[k][j] );
printf( " %10.3le ", vecX[k] );
for ( j=0; j<n; j++ ) printf( " %10.3le", matE[k][j] );
}
printf( "\n" );
}


int main( int argc, char **argv )
{
double C, D, ratio, big;
int n, i, j, k,r, pivot_row;
char copyright[] = "Written by M.Takata: 1999/01/20\n";
int ch, verbose;

verbose = 0;
while (( ch = getopt(argc,argv,"vh")) != -1 ) {
if ( ch == 'v' ) verbose = 1;
if ( ch == 'h' ) {
printf( "Options:\n" );
printf( "\t-v: dumps intermediate results.\n" );
printf( "\t-h: shows help message, as you see now.\n" );
printf( "\n" );
return 1;
}
}


/* データ入力 */
scanf( "%d", &n );
if ( n > SIZE ) {
printf("N should be less than or equal to %d.\n", SIZE);
return 1;
}

for ( i=0; i<n; i++ ) {
for ( j=0; j<n; j++ ) {
scanf( "%lf", &(matA[i][j]) ); }
scanf( "%lf", &(vecX[i]) );
}


printf("Inputs:\n");

for ( i=0; i<n; i++ ){
for ( j=0; j<n; j++ ){
printf( " %12.5le", matA[i][j] ); }
printf( " | %12.5le\n", vecX[i] );
}


/* 単位行列の設定 */
for ( i=0; i<n; i++ ){
for ( j=0; j<n; j++ ){ matE[i][j] = 0.0; }
matE[i][i] = 1.0;
}

/* Solve */
if ( verbose ) dump( n );


for ( i=0; i<n; i++ ){

printf("############ i=%d", i);
if ( verbose ) dump( n );
printf("-------\n");

C = matA[i][i];
printf("C=%.2f\n", C);

/* ピボット選択による並べ替え(始) */

big = 0.0;

/* i列で最大値を検索 */

if ( C < 1.0e-04 && -1.0e-04 < C ) {
for(k=0; k<n; k++){
if(fabs(matA[k][i]) > big){
big = fabs(matA[k][i]);
pivot_row = k;
}
}
/* 検索した最大値が0または小さい値であったら計算できないと返す */

if(big < 1.0e-04 && -1.0e-04 < big){
printf("Invalid matrix");
return 0;
}

/*最大値を含む行と入れ替え*/

if(i != pivot_row){
for(k=0; k<n; k++){
D = matA[i][k];
matA[i][k] = matA[pivot_row][k];
matA[pivot_row][k] = D;

}
D = vecX[k];
vecX[k] = vecX[pivot_row];
vecX[pivot_row] = D;

}

}
/*ピボット選択による並べ替え(終)*/

/********************************************/
/** i 列目前半の処理 **/
/********************************************/

for ( k=0; k<n; k++ ){
matA[i][k] = matA[i][k] / C;
matE[i][k] = matE[i][k] / C;
}
vecX[i] = vecX[i] / C;

if ( verbose ) dump( n );

for ( k=0; k<n; k++ ){ /* 行のループ */
if ( i != k ){
D = matA[k][i];

/********************************************/
/** i 列目後半の処理 **/
/********************************************/
for ( j=0; j<n; j++ ){ /* 列のループ */
matA[k][j] = matA[k][j] - D * matA[i][j];
matE[k][j] = matE[k][j] - D * matE[i][j];
}
vecX[k] = vecX[k] - D * vecX[i];
}
}

if ( verbose ) dump( n );
}


/* 計算結果出力 */
printf( "\nResults:\n" );
for ( i=0; i<n; i++ ) printf( " %12.5le", vecX[i] );
printf( "\nInverse:\n" );
for ( i=0; i<n; i++ ){
for ( j=0; j<n; j++ ) printf( " %12.5le", matE[i][j] );
printf("\n");
}
printf("\n");

}

閲覧数:
40
回答数:
1
お礼:
250枚

違反報告

ベストアンサーに選ばれた回答

2019/8/321:21:56

以下のように修正します。
main に詰め込みすぎです。
https://ideone.com/xB1hyr


#define _CRT_SECURE_NO_WARNINGS // Visual Studio のときは入れてください
#include <math.h>
#include <stdio.h>

#ifndef SIZE
#define SIZE (16)
#endif

double matA[SIZE][SIZE], matE[SIZE][SIZE], vecX[SIZE];

void dump(int n) {
    int j, k;
    for (k = 0; k < n; k++) {
        printf("\n");
        for (j = 0; j < n; j++)
            printf(" %10.3le", matA[k][j]);
        printf(" %10.3le ", vecX[k]);
        for (j = 0; j < n; j++)
            printf(" %10.3le", matE[k][j]);
    }
    printf("\n");
}

int main(int argc, char **argv) {
    double C, D, ratio, big;
    int n, i, j, k, r, pivot_row;
    char copyright[] = "Written by ********: ****/**/**\n";
    int ch, verbose;

// getopt は Visual Studio に含まれていない。自作するか探すこと
    while ((ch = getopt(argc, argv, "vh")) != -1) {
        if (ch == 'v')
            verbose = 1;
        if (ch == 'h') {
            printf("Options:\n");
            printf("\t-v: dumps intermediate results.\n");
            printf("\t-h: shows help message, as you see now.\n");
            printf("\n");
            return 1;
        }
    }
// ターミナルから指定できないので
    verbose = 1;

    /* データ入力 */
    scanf("%d", &n);
    if (n > SIZE) {
        printf("N should be less than or equal to %d.\n", SIZE);
        return 1;
    }

    for (i = 0; i < n; i++) {
        for (j = 0; j < n; j++) {
            scanf("%lf", &(matA[i][j]));
        }
        scanf("%lf", &(vecX[i]));
    }

    printf("Inputs:\n");

    for (i = 0; i < n; i++) {
        for (j = 0; j < n; j++) {
            printf(" %12.5le", matA[i][j]);
        }
        printf(" | %12.5le\n", vecX[i]);
    }

    /* 単位行列の設定 */
    for (i = 0; i < n; i++) {
        for (j = 0; j < n; j++) {
            matE[i][j] = 0.0;
        }
        matE[i][i] = 1.0;
    }

    /* Solve */
    if (verbose)
        dump(n);

    for (i = 0; i < n; i++) {

        printf("############ i=%d", i);
        if (verbose)
            dump(n);
        printf("-------\n");

// C はここで表示すべきでない。ピボット選択による並べ替えが終わった後。
//      C = matA[i][i];
//      printf("C=%.2f\n", C);

        /* ピボット選択による並べ替え(始) */

        big = 0.0;

        /* i列で最大値を検索 */

//      if (C < 1.0e-04 && -1.0e-04 < C) { // 削除。ここで判定するとピボット選択が始まらない
            for (k = i; k < n; k++) { // k = 0 を k = i に訂正。既にピボット選択が終わった行は残す
                if (fabs(matA[k][i]) > big) {
                    big = fabs(matA[k][i]);
                    pivot_row = k;
                }
            }
            /* 検索した最大値が0または小さい値であったら計算できないと返す */

            if (big < 1.0e-04 && -1.0e-04 < big) {
                printf("Invalid matrix");
                return 0;
            }

            /*最大値を含む行と入れ替え*/

            if (i != pivot_row) {
                for (k = 0; k < n; k++) {
                    D = matA[i][k];
                    matA[i][k] = matA[pivot_row][k];
                    matA[pivot_row][k] = D;

                    D = matE[i][k]; // matE の入れ替えが抜けている
                    matE[i][k] = matE[pivot_row][k];
                    matE[pivot_row][k] = D;
                }
                D = vecX[i]; // [k] を [i] に訂正
                vecX[i] = vecX[pivot_row];
                vecX[pivot_row] = D;
            }
//      } // 削除
        /*ピボット選択による並べ替え(終)*/

        C = matA[i][i];        // ここに移動
        printf("C=%.2f\n", C); // ここに移動

        /********************************************/
        /** i 列目前半の処理 **/
        /********************************************/

        for (k = 0; k < n; k++) {
            matA[i][k] = matA[i][k] / C;
            matE[i][k] = matE[i][k] / C;
        }
        vecX[i] = vecX[i] / C;

        if (verbose)
            dump(n);

        for (k = 0; k < n; k++) { /* 行のループ */
            if (i != k) {
                D = matA[k][i];

                /********************************************/
                /** i 列目後半の処理 **/
                /********************************************/
                for (j = 0; j < n; j++) { /* 列のループ */
                    matA[k][j] = matA[k][j] - D * matA[i][j];
                    matE[k][j] = matE[k][j] - D * matE[i][j];
                }
                vecX[k] = vecX[k] - D * vecX[i];
            }
        }

        if (verbose)
            dump(n);
    }

    /* 計算結果出力 */
    printf("\nResults:\n");
    for (i = 0; i < n; i++)
        printf(" %12.5le", vecX[i]);
    printf("\nInverse:\n");
    for (i = 0; i < n; i++) {
        for (j = 0; j < n; j++)
            printf(" %12.5le", matE[i][j]);
        printf("\n");
    }
    printf("\n");
}

この回答は投票によってベストアンサーに選ばれました!

あわせて知りたい

この質問につけられたタグ

みんなで作る知恵袋 悩みや疑問、なんでも気軽にきいちゃおう!

Q&Aをキーワードで検索:

Yahoo! JAPANは、回答に記載された内容の信ぴょう性、正確性を保証しておりません。
お客様自身の責任と判断で、ご利用ください。
本文はここまでです このページの先頭へ

「追加する」ボタンを押してください。

閉じる

※知恵コレクションに追加された質問は選択されたID/ニックネームのMy知恵袋で確認できます。

不適切な投稿でないことを報告しました。

閉じる