C言語ガウスの消去法

秋学期の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

行動規範の内容に同意します

回答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;`

追記

C
1    /* 前進消去（ピボット選択）*/
2    for(k = 0; k < N-1; k++){  /* 第ｋステップ */
3        m = k;
4        mmax = fabs(*(ak+k));
5		ai = ak;                          // ★
6        for(i = k+1; i < N; i++){  /* ピボット選択 */
7		   ai += N;                       // ★
8           if(fabs(*(ai+k)) > mmax){      // ★
9              m = i;
10              mmax = fabs(*(ai+k));       // ★
11           }
12        }
13        printf("column %d: row %d is max\n", k, m);
14        if(m != k){  /* 第ｋ行と第m行を交換　*/
15		   double *am = A + m*N;          // ★
16           for(i = k; i < N; i++){/* 係数行列　*/
17               temp = *(ak+i);
18               *(ak+i) = *(am+i);         // ★
19               *(am+i) = temp;            // ★
20           }
21           /* 既知ベクトル */
22           temp = b[k];
23           b[k] = b[m];
24           b[m] = temp;
25        }
26        /* 第k列消去 */
27		ai = ak;                          // ★
28        for (i = k + 1; i < N; i++){
29			ai += N;                      // ★
30            r = *(ai+k) / *(ak+k);        // ★
31            for (j = k; j < N; j++){
32                *(ai+j) -= *(ak+j) * r;   // ★
33            }
34            b[i] -= b[k] * r;
35        }
36        PrintAb(A, b, N);
37        ak+=N;
38    }