###### tags: `code`
# c14086028-homework-6-2
###### tags: `use'C' to code the 'L-AL:Gaussian Elimination' to find inverse`

#include<stdio.h>
#include<stdlib.h>
# define N 20
void initial(double [N][N] ,double [N][N], double* );
void clear(double[N][N],double [N][N],double [N]);
void clear2(double[N][N],double [N][N], double[N]);
double matrixsolver(double[N][N], double [N][N],double[N]);
int main() {
    double u = 0;
    double A[N][N] = { {} };
    double B[N][N] = { {} };
    double F[N] = {};
    initial(A,B,F);
    /*print initial*/
    printf("initial linear problem is \n");
    printf("-----------\n");
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
            printf("%1.3lf\t", A[i][j]);
        }
        printf("\n");
    }
    printf("-----------\n");
    for (int i = 0; i < N; i++) {
        printf("%1.3lf\n", F[i]);
    }
    printf("\n");
    printf("-----------\n");

    /*solve and print it */

    u=matrixsolver(A,B,F);
    if (u == 0) {
        printf("solution is\n");
        printf("-----------\n");
        for (int i = 0; i < N; i++) {
            for (int j = 0; j < N; j++) {
                printf("%1.3lf\t", A[i][j]);
            }
            printf("\n");
        }
        printf("inverse of A is\n");
        printf("-----------\n");
        for (int i = 0; i < N; i++) {
            for (int j = 0; j < N; j++) {
                printf("%1.3lf\t", B[i][j]);
            }
            printf("\n");
        }
        printf("-----------\n");
        for (int i = 0; i < N; i++) {
            printf("%1.3lf\n", F[i]);
        }
    }
    else {
        printf("it has a lot of solution\n");
    }
    


    system("pause");
    return 0;
}

void initial(double A[N][N],double B[N][N],double F[N]) {
    for (int i = 0; i < N; i++) {
        A[i][i] = -2;
        if ((i + 1) < N) {
            A[i][i + 1] = 1;
            A[i + 1][i] = 1;
        }
        F[i] = (i == 19);
    }
    for (int i = 0; i < N; i++) {
        B[i][i] = 1;
    }
}

void clear(double A[N][N],double B[N][N],double F[N]) {
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
            if (A[i][j] != 0) {
                double t = 0;
                t = A[i][j];
                for (int k = 0; k < N; k++) {
                    A[i][k] = A[i][k] / t;
                }
                for (int k = 0; k < N; k++) {
                    B[i][k] = B[i][k] / t;
                }
                F[i] = F[i] / t;
                break;
            }
            else {
                continue;
            }
        }
   }
    
}

void clear2(double A[N][N],double B[N][N] ,double F[N]) {
    for (int i = N - 1; i >= 0; i--) {
        for (int j = N - 1; j >= 0; j--) {
            if (A[i][j] != 0) {
                double t = 0;
                t = A[i][j];
                for (int k = N-1; k >= 0; k--) {
                    A[i][k] = A[i][k] / t;
                }
                for (int k = N-1; k >= 0; k--) {
                    B[i][k] = B[i][k] / t;
                }
                F[i] = F[i] / t;
                break;
            }
            else {
                continue;
            }
        }
    }
  
}

double matrixsolver(double A[N][N], double B[N][N],double F[N]) {
    clear(A,B,F);
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
            if (A[i][j] != 0) {
                for (int k = i; k < N; k++) {
                    if (A[k][j] != 0 && k!=i) {

                        for (int l = 0; l < N; l++) {
                            
                            A[k][l] = A[k][l] - A[i][l];

                        }
                        for (int l = 0; l < N; l++) {

                            B[k][l] = B[k][l] - B[i][l];

                        }
                            F[k] = F[k] - F[i];
                    }
                    else {
                        continue;
                    }
                }
                break;
            }
            else {
                continue;
            }
        }
        clear(A,B,F);
    }

    clear2(A,B,F);
    for (int i = N-1; i >=0; i--) {
        for (int j = N-1; j >= 0; j--) {
            if (A[i][j] != 0) {
                for (int k = i; k >=0 ; k--) {
                    if (A[k][j] != 0 && k != i) {

                        for (int l = 0; l < N; l++) {

                            A[k][l] = A[k][l] - A[i][l];
                        }
                        for (int l = 0; l < N; l++) {

                            B[k][l] = B[k][l] - B[i][l];
                        }
                        F[k] = F[k] - F[i];
                    }
                    else {
                        continue;
                    }
                }
                break;
            }
            else {
                continue;
            }
        }
        clear2(A,B,F);
    }
   
    for (int e = 0; e < N; e++) {
        if (A[e][e] == 0) {
            return 1;
        }
    }
    return 0;
}