行列×逆行列＝良い！

行列色々

行列の計算する時はMatrix構造体作ってやってるけど無駄かもね～
そういえば行列式が0の時を考慮してなかったな・・・

多分InfInfって出ると思うけど、まぁいっか

行列式の関数は一応あるけど、使ってないよ～

相変わらず色々なサイトから引っ張ってきて、都合のいいように改変してますごめんなさい

行列の積を求める方法ってのは、基本に忠実にやればO（n³）なんですが、ストラッセンの行列積アルゴリズムを使えば大体O（n^2.81）になるそうです

あのアルゴリズム実装するの難しそう。。。。っていうかMATLAB使っちゃえばいいんだけどね！

MATLAB最高！畳み込み最高！

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

typedef struct {
    int n;
    double **m;
} Matrix;

Matrix *create_matrix(int n);
Matrix *copy_matrix(Matrix *m);
void free_matrix(Matrix *m);
Matrix *inverse_matrix(Matrix *m);
Matrix *multi_matrix(Matrix *m1, Matrix *m2);
double det_matrix(Matrix *m);
void print_matrix(Matrix *mat);

int main(void)
{
    Matrix *mat, *inv_mat, *mul_mat;
    int i, k, n;

    printf("n -> ");
    scanf("%d", &n);

    mat = create_matrix(n);

    srand(time(NULL));
    for (i = 0;i < mat->n;i++) {
        for (k = 0;k < mat->n;k++) {
            mat->m[i][k] = (double)(rand()%150)/(rand()%50+1);
            //printf("%d行%d列目 : ", i+1, k+1);
            //scanf("%lf", &mat->m[i][k]);
        }
    }
    
    printf("行列A\n");
    print_matrix(mat);
    
    inv_mat = inverse_matrix(mat);
    
    printf("逆行列A-1\n");
    print_matrix(inv_mat);
    
    mul_mat = multi_matrix(mat, inv_mat);
    
    printf("AA-1\n");
    print_matrix(mul_mat);
    
    free_matrix(mat);
    free_matrix(inv_mat);
    free_matrix(mul_mat);

    return 0;
}


Matrix *create_matrix(int n)
{
    int i;
    Matrix *mat;
    
    mat = (Matrix*)malloc(sizeof(Matrix));
    if (mat == NULL) {
        printf("Error : Allocation error\n");
        exit(1);    
    }
    mat->m = (double**)calloc(n ,sizeof(double*));
    if (mat->m == NULL) {
        printf("Error : Allocation error\n");
        exit(1);        
    }
    for (i = 0;i < n;i++) {
        mat->m[i] = (double*)calloc(n, sizeof(double));
        if (mat->m[i] == NULL) {
            printf("Error : Allocation error\n");
            exit(1);
        }
    }
    mat->n = n;
    
    return mat;
}

Matrix *copy_matrix(Matrix *mat)
{
    Matrix *cpy_mat;
    int i, k;
    
    cpy_mat = create_matrix(mat->n);
    for (i = 0;i < mat->n;i++)
        for (k = 0;k < mat->n;k++)
            cpy_mat->m[i][k] = mat->m[i][k];
    cpy_mat->n = mat->n;
    
    return cpy_mat;
}

Matrix *inverse_matrix(Matrix *m)
{
    Matrix *inv_mat, *mat;
    double tmp;
    int i, j, k, n = m->n;
    
    mat = copy_matrix(m);
    inv_mat = create_matrix(n);
    for (i = 0;i < n;i++)
        inv_mat->m[i][i] = 1.0;

    for (i = 0;i < n;i++) {
        tmp = 1 / mat->m[i][i];
        for (j = 0;j < n;j++) {
            mat->m[i][j] *= tmp;
            inv_mat->m[i][j] *= tmp;
        }
        for (j = 0;j < n;j++) {
            if (i != j) {
                tmp = mat->m[j][i];
                for (k = 0;k < n;k++) {
                    mat->m[j][k] -= mat->m[i][k]*tmp;
                    inv_mat->m[j][k] -= inv_mat->m[i][k]*tmp;
                }
            }
        }
    }
    free_matrix(mat);
    
    return inv_mat;
}

double det_matrix(Matrix *m)
{
    int i, j, k, n = m->n;
    double det = 1.0, tmp;
    Matrix *mat;
    
    mat = copy_matrix(m);
    
    for (i = 0;i < n;i++) {
        for (j = 0;j < n;j++) {
            if (i < j) {
                tmp = mat->m[j][i] / mat->m[i][i];
                for (k = 0;k < n;k++)
                    mat->m[j][k] -= mat->m[i][k] * tmp;
            }
        }
    }
    for (i = 0;i < n;i++)
        det *= mat->m[i][i];

    free_matrix(mat);

    return det;
}

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

void print_matrix(Matrix *mat)
{
    int i, k;

    for (i = 0;i < mat->n;i++) {
        for (k = 0;k < mat->n;k++)
            printf("%8.3f", mat->m[i][k]);
        printf("\n");
    }
}

Matrix *multi_matrix(Matrix *m1, Matrix *m2)
{
    int i, j, k, n = m1->n;
    Matrix *result;

    result = create_matrix(n);
    
    for (i = 0;i < n;i++) {
        for (j = 0;j < n;j++) {
            for (k = 0;k < n;k++) {
                result->m[i][j] += m1->m[i][k] * m2->m[k][j];
            }
        }
    }
    result->n = n;
    
    return result;
}