ガウス

自分用メモ

http://home.f01.itscom.net/toge/programingreport/articlerecord.html

http://www9.plala.or.jp/sgwr-t/c_sub/ConvHTML.html

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

typedef struct {
    int h;
    int w;
    double **m;
} Matrix;

void error_message(char *m)
{
    printf("エラー : %s\n", m);
    exit(1);
}

Matrix *create_matrix(int height, int width)
{
    Matrix *mat;
    int i;
    
    mat = (Matrix*)malloc(sizeof(mat));
    if (mat == NULL)
        error_message("Matrix メモリ確保失敗");
    mat->m = (double**)malloc(sizeof(double*)*height);
    if (mat->m == NULL)
        error_message("m メモリ確保失敗");
    for (i = 0;i < width;i++) {
        mat->m[i] = (double*)malloc(sizeof(double)*width);
        if (mat->m[i] == NULL)
            error_message("各行 メモリ確保失敗");
    }
    mat->h = height;
    mat->w = width;
    
    return mat;
}

void free_matrix(Matrix *mat)
{
    int i;
    
    for (i = 0;i < mat->h;i++)
        free(mat->m[i]);
    free(mat->m);
    free(mat);
}

int gauss(Matrix *temp)
{
    int i, j, k, d;
    double alpha;
    double copy;
    double max;
    int size;

    if(temp->h > temp->w)
        size = temp->w;
    else
        size = temp->h;
    
    for(i=0; i<size; i++){
        d = i;
        max = fabs(temp->m[i][i]);

        for(j=i+1; j<temp->h; j++){
            if(max < fabs(temp->m[j][i])){
                d = j;
                max = fabs(temp->m[j][i]);
            }
        }
        for(j=i; j<temp->w; j++){
            copy = temp->m[i][j];
            temp->m[i][j] = temp->m[d][j];
            temp->m[d][j] = copy;
        }

        for(j=i+1; j<temp->h; j++){
            if(!temp->m[i][i]) {
                return 1;
            }
            alpha = temp->m[j][i] / temp->m[i][i];

            for(k=i; k<temp->w; k++){
                temp->m[j][k] -= alpha*temp->m[i][k];
            }
        }
    }

    return 0;
}

double *backward(Matrix *matrix)
{
    int i, j;
    double b;
    double *temp;
    int width, height, size;

    width = matrix->w;
    size = height = matrix->h;

    temp = (double*)malloc(sizeof(double) * height);
    temp[size-1] = matrix->m[height-1][width-1] / matrix->m[height-1][width-2];

    for(i = size-2; i >= 0; i--){
        b = matrix->m[i][width-1];
        for(j=i+1; j<size; j++){
            b -= matrix->m[i][j]*temp[j];
        }
        if(!matrix->m[i][i]) {
            return NULL;
        }
        temp[i] = b/matrix->m[i][i];
    }

    return temp;
}

int main(void)
{
    Matrix *mat;
    int height, width;
    int i, k;
    double *d;

    scanf("%d %d", &height, &width);
    mat = create_matrix(height, width);
    srand(time(NULL));
    for (i = 0;i < height;i++) {
        for (k = 0;k < width;k++) {
            scanf("%lf", &mat->m[i][k]);
            printf("%6.2f ", mat->m[i][k]);
        }
        printf("\n");
    }
    gauss(mat);
    for (i = 0;i < height;i++) {
        for (k = 0;k < width;k++) {
            printf("%6.2f ", mat->m[i][k]);
        }
        printf("\n");
    }
    d = backward(mat);
    for (i = 0;i < mat->h;i++)
        printf("x%2d = %lf\n",i+1,d[i]);
    
    return 0;
}