掃き出し

別に掃除するってわけじゃあないよ

ガウスの消去法のプログラム。

ピボット選択をとれば、単純消去法にチェンジ出来る！（する必要はないけど）

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

#define N 3

void print_vector(double b[], int n);
void print_matrix(double a[], int n);
int gauss(double a[], double b[], double x[], int n);

int main(void)
{
    double a[N*N] = {
        2, 3, -1,
        1, 1, 1,
        3, -1, -2
    };
    double b[N] = {-2, 7, 0};
    double x[N];
    int i;

    printf("A = \n");
    print_matrix(a, N);
    printf("b = \n");
    print_vector(b, N);
    if (gauss(a, b, x, N))
        for (i = 0;i < N;i++)
            printf("x%d = %6.2f\n", i+1, x[i]);
    else
        printf("計算できませんでした.\n");

    return 0;
}

void print_vector(double b[], int n)
{
    int i;
    for (i = 0;i < n;i++)
        printf("%6.2f ", b[i]);
    printf("\n");
}

void print_matrix(double a[], int n)
{
    int i;
    for (i = 0;i < n;i++)
        print_vector(&a[i*n], n);
}

int gauss(double a[], double b[], double x[], int n)
{
    int i, j, k, idx;
    double pivot, tmp;
    double eps = 0.00001;

    for (i = 0;i < n;i++) {
        //ピボットの選択

        pivot = a[i*n+i];
        idx = i;
        for (j = i+1;j < n;j++) {
            if (pivot < a[j*n+i]) {
                pivot = a[j*n+i];
                idx = j;
            }
        }
        //ピボットが0しかなかったら計算できないのでリターン

        if (fabs(pivot) < eps)
            return 0;
        //行の交換

        if (idx != i) {
            tmp = b[i];
            b[i] = b[idx];
            b[idx] = tmp;
            for (j = i;j < n;j++) {
                tmp = a[idx*n+j];
                a[idx*n+j] = a[i*n+j];
                a[i*n+j] = tmp;
            }
        }
        //前進消去

        b[i] /= a[i*n+i];
        for (j = n-1;j >= i;j--)
            a[i*n+j] /= a[i*n+i];
        for (j = i+1;j < n;j++) {
            b[j] -= b[i] *  a[j*n+i] / a[i*n+i];
            for (k = n-1;k >= i;k--)
                a[j*n+k] -= a[i*n+k] *  a[j*n+i] / a[i*n+i];
        }
    }

    //後退代入

    for (i = n-1;i >= 0;i--) {
        x[i] = b[i];
        for (j = n-1;j > i;j--)
            x[i] -= x[j] * a[i*n+j];
    }
    return 1;
}

なんか微妙にインデントがずれてるんだが・・・・・

今度はガウス・ザイデルとSORをやると思う