京大理学部 物理科学課題演習 C1 （その５）

#ピボット選択

　[前回](https://qiita.com/sasakitakanori/items/a755414cfb221033bc9d)の課題で見たとおり、対角要素が $0$ になった場合そのままでは計算を続けることができません。そこで、方程式の順番を入れ替えることになります。これは拡大行列の行を入れ替えることを意味します。ところで、このような入れ替えは必ずしも対角要素が $0$ のときしか行ってはいけないというものでもありません。一般には $k$ 行目の方程式を軸に消去を行う場合、$a_{ik}\ (i = k, \cdots, n)$ の中から絶対値が最大のものを選び、それを含む行と $k$ 行を入れ替えてから消去を行います。これを**部分ピボット選択**とよびます。

　一方、$a_{ik}\ (i = k, \cdots, n)$ の中からだけでなく、$a_{ij}\ (i = k, \cdots, n,\ \ j = k, \cdots, n)$ の中から絶対値が最大のものを選び（それを $a_{imax\ jmax}$ としましょう）、$k$ 行と $imax$ 行、$k$ 列と $jmax$ 列を入れ替える方法を**完全ピボット選択**とよびます。

　部分ピボット選択の方は、単なる方程式の入れ替えですから単純です。しかし完全ピボット選択の方は、列の入れ替えを含むためプログラムがかなり複雑になります。というのは、列の入れ替えは変数の順番の入れ替えに相当するため、最後に得られた解が最初の順番と変わってしまうからです。したがって、何列目と何列目を入れ替えたのかを記憶しておいて、最後に解の順番を並べ替えなければなりません。

　ここでは簡単なピボット選択のプログラムを作ってみましょう。

```c
#include <stdio.h>
#include <math.h>
#define N 3

void simeq(double *apt, double *bpt, int n);

int main(void)
{
    static double a[N][N] = {0.0, 1.0, 1.0, 2.0, 3.0, 1.0, 1.0, 1.0, 3.0};
    static double b[N] = {2.0, 6.0, 5.0};

    int i, j, k;

    simeq(&a[0][0], b, N);

    for (k = 0; k < N; k++) {
        printf("x%d = %4.1lf\n", k, b[k]);
    }
    for (i = 0; i < N; i++) {
        for (j = 0; j < N; j++) {
            printf("%4.1lf", a[i][j]);
        }
        printf("\n");
    }

    return 0;
}

/* 連立一次方程式を解くための関数（ガウス・ジョルダンの消去法、部分ピボット版） */
void simeq(double *apt, double *bpt, int n)
{
    int i, j, k, ipivot, ip;
    double amax, atmp, btmp;

    for (k = 0; k < n; k++) {
        ipivot = k;
        amax = fabs(*(apt + k * n + k));
        for (ip = k + 1; ip < n; ip++) {
            if (amax < fabs(*(apt + ip * n + k))) {
                ipivot = ip;
                amax = fabs(*(apt + ip * n + k));
            }
        }
        if (ipivot != k) {
            for (j = 0; j < n; j++) {
                atmp = *(apt + k * n + j);
                *(apt + k * n + j) = *(apt + ipivot * n + j);
                *(apt + ipivot * n + j) = atmp;
            }
            btmp = *(bpt + k);
            *(bpt + k) = *(bpt + ipivot);
            *(bpt + ipivot) = btmp;
        }
        for (j = k + 1; j < n; j++) {
            *(apt + k * n + j) = *(apt + k * n + j) / *(apt + k * n +k);
        }
        *(bpt + k) = *(bpt + k) / *(apt + k * n + k);
        for (i = 0; i < n; i++) {
            if (i != k) {
                for (j = k + 1; j < n; j++) {
                    *(apt + i * n + j) = *(apt + i * n + j) - *(apt + i * n + k) * *(apt + k * n + j);
                }
                *(bpt + i) = *(bpt + i) - *(apt + i * n + k) * *(bpt + k);
            }
        }
    }
}
```

　このプログラムをコンパイルする際には、`math.h` を `include` しているので

```console
$ gcc ipivot.c -lm -o ipivot
```

のように `-lm` オプションを忘れないようにしましょう。

　ところで、このプログラムは係数行列が正則な場合にしか対応していません。係数行列の行列式が $0$ の場合、異常終了のコンディションコード（ `icond = 1` ）を返すようにプログラムを変更しましょう。正常終了の場合は `icond = 0` とします。まず引数の並びは

```c
メイン関数側： 　simeq(&a[0][0], b, N, icond);
simeq 関数側：　void simeq(double *apt, double *bpt, int n, int icond);
```

のように変わります。さて、係数行列が非正則であるかどうかをどのように判断すればよいのでしょうか。実際に行列式を計算するのは面倒ですから、ここでは選択されたピボット $a_{ipivot\ k}$ を用いて判断しましょう。すなわち $a_{ipivot\ k}$ が $0$ かあるいは $0$ に近い値の場合、係数行列は非正則だと判断するわけです。ピボットを選択する `for` ループの後で

```c
if (fabs(*(apt + ipivot * n + j)) < 1.0e-5) {
    icond = 1;
    return;
}
```

といった手続きを入れればよいわけです。

　一般には部分ピボット選択で十分です。完全ピボット選択の場合は最大値の探索範囲が広がるため計算時間も増大します。したがって、完全ピボット選択についてはここでは解説しません。興味のある人は適当な教科書等を参照して勉強してみてください。

##課題

　係数行列が非正則の場合にも対応したプログラムを完成させなさい。