C言語 ガウスの消去法

Question

秋学期のC言語の授業のために自主学習として勉強しているものです。

キーボードから次元 N を入力した後，N × N 行列 A と N 次元ベクトル b のデータを要素 ごとにキーボードから入力し，その行列・ベクトルを使って，連立 1 次方程式 Ax = b を 部分ピボット選択付きガウスの消去法を用いて解け.
以上のような問いに対し、以下のソースコードを書きました。
質問をしてから、何度もやり直しや修正を図り少しづつ答えには近づいてきたものの、まだ正確にはうまくいきません。よろしくお願いいたします。

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

void PrintAb(double *A, double *b, int n){
    int i,j;
    double *ai;
    ai= A;
    printf("A, b = \n");
    for( i = 0 ; i < n ; i++){
        for( j = 0 ; j < n ; j++){
            printf("%6.1f\t", *(ai+j));
        }
        printf("|%6.1f\n",*(b+i));
        ai+=n;
    }
}


int main(void) {
    int i, j, k, m, N;
    double r;
    double *ai, *ak;
    double *A, *b;
    double mmax, temp;
    printf("N= ");
    scanf("%d",&N);

    A = (double *)malloc(N*N*sizeof(double));
    if(A==NULL){
        printf("Can’t allocate memory. Y=n");
        exit(1);
    }

   b = (double *)malloc(N*sizeof(double));
    if(b==NULL){
        printf("Can’t allocate memory. Y=n");
        exit(1);
    }

    ai = A;
    ak = A;
    for(i = 0 ; i < N ; i++ ){
        for( j = 0 ; j < N ; j++ ){
            printf("A[%d][%d] = ", i , j );
            scanf("%lf", ai+j );
        }
        ai+=N;
    }

    for(i = 0 ; i < N ; i++ ){
        printf("b[%d] = ", i );
        scanf("%lf", b+i);
    }
    printf("\n");
    
    PrintAb(A, b, N);

    for(k = 0; k < N-1; k++){
        m = k;
        mmax = fabs(*(ak+k));
        for(i = k+1; i < N; i++){
           if(fabs(*(ai+i*N+k)) > mmax){
              m = i;
               mmax = fabs(*(ai+i*N+k));
           }
        }
        printf("column %d: row %d is max\n", k, m);
        if(m != k){
           for(i = k; i < N; i++){
               temp = *(ak+i);
               *(ak+i) = *(ak+m*N+i);
               *(ak+m*N+i) = temp;
           }
           temp = b[k];
           b[k] = b[m];
           b[m] = temp;
        }
        ai = A;
        ak = A;
        for (i = k + 1; i < N; i++){
            r = *(ai+i*N+k) / *(ak+k);
            for (j = k; j < N; j++){
                *(ai+i*N+j) -= *(ak+j) * r;
            }
            b[i] -= b[k] * r;
        }
        PrintAb(A, b, N);
        ak+=N;
    }

    ai = ak;
    for (i = N - 1; i >= 0; i--) {
        for (j = i + 1; j < N; j++) {
            b[i] -= *(ai+j) * b[j];
        }
        b[i] /= *(ai+i);
        PrintAb(A, b,N);
        ai-=N;
    }

    /*メモリの解放*/
    free(A);
    free(b);

    return 0;
}

しかし、私自身が書いたソースコードを実行しても正しくいきません(エラーなどは出ていません)。
(第二回目の消去までは正しいのですが、それ以降が正しくありません。)
どなたかわかる方がいらっしゃれば教えていただけると幸いです。
=正しい出力=
N= 3
A[0][0] = 1
A[0][1] = 1
A[0][2] = 2
A[1][0] = 3
A[1][1] = 2
A[1][2] = -1
A[2][0] = 2
A[2][1] = -1
A[2][2] = 1
b[0] = 4
b[1] = -1
b[2] = 5

A, b =
1.0 1.0 2.0 | 4.0
3.0 2.0 -1.0 | -1.0
2.0 -1.0 1.0 | 5.0
column 0: row 1 is max
A, b =
3.0 2.0 -1.0 | -1.0
0.0 0.3 2.3 | 4.3
0.0 -2.3 1.7 | 5.7
column 1: row 2 is max
A, b =
3.0 2.0 -1.0 | -1.0
0.0 -2.3 1.7 | 5.7
0.0 -0.0 2.6 | 5.1
A, b =
3.0 2.0 -1.0 | -1.0
0.0 -2.3 1.7 | 5.7
0.0 -0.0 2.6 | 2.0
A, b =
3.0 2.0 -1.0 | -1.0
0.0 -2.3 1.7 | -1.0
0.0 -0.0 2.6 | 2.0
A, b =
3.0 2.0 -1.0 | 1.0
0.0 -2.3 1.7 | -1.0
0.0 -0.0 2.6 | 2.0

Answer 1

m = b[k]; b[k] = b[m]; b[m] = m;

これって、明らかに変ですよね。b[k] と b[m] を交換するなら、

       `double t = b[k]; b[k] = b[m]; b[m] = t;`

追記

    /* 前進消去（ピボット選択）*/
    for(k = 0; k < N-1; k++){  /* 第ｋステップ */
        m = k;
        mmax = fabs(*(ak+k));
		ai = ak;                          // ★
        for(i = k+1; i < N; i++){  /* ピボット選択 */
		   ai += N;                       // ★
           if(fabs(*(ai+k)) > mmax){      // ★
              m = i;
              mmax = fabs(*(ai+k));       // ★
           }
        }
        printf("column %d: row %d is max\n", k, m);
        if(m != k){  /* 第ｋ行と第m行を交換　*/
		   double *am = A + m*N;          // ★
           for(i = k; i < N; i++){/* 係数行列　*/
               temp = *(ak+i);
               *(ak+i) = *(am+i);         // ★
               *(am+i) = temp;            // ★
           }
           /* 既知ベクトル */
           temp = b[k];
           b[k] = b[m];
           b[m] = temp;
        }
        /* 第k列消去 */
		ai = ak;                          // ★
        for (i = k + 1; i < N; i++){
			ai += N;                      // ★
            r = *(ai+k) / *(ak+k);        // ★
            for (j = k; j < N; j++){
                *(ai+j) -= *(ak+j) * r;   // ★
            }
            b[i] -= b[k] * r;
        }
        PrintAb(A, b, N);
        ak+=N;
    }

Comments

[pocpocpoc]

ありがとうございます。
そこの部分に関しては私自身のミスで明らかにおかしかったです。修正します。

[pocpocpoc]

追記ありがとうございます。
難しいですが、理解できるように頑張ります。