/**
 * This is homework5, which deals with A equation using
 * 1. Gauss-Jordan Elimination
 * 2. LU Decomposition
 * 3. Singular Value Decomposition
 **/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

// compile with gcc -g -o main *.c -lm

// 1. Ax=b 풀기
// 여기에서 Ai와 bi는 course homepage에서 확인 가능
// 2. iterative improvement method
// 3. A의 inverse와 determinant 구하기

int gaussj(float **a, int n, float **b, int m);
void ludcmp(float **a, int n, int *indx, float *d);
void lubksb(float **a, int n, int *indx, float *d);
void svdcmp(float **a, int m, int n, float w[], float **v);
void svbksb(float **u, float w[], float **v, int m, int n, float b[], float x[]);
void mprove(float **a, float **alud, int n, int indx[], float b[], float x[]);

int main(void){
    int three = 1;
    for(; three <= 3; three++){
        char filename[20];
        sprintf(filename, "lineq%d.dat", three);
        printf("\n*********************\nfile %s:\n*********************", filename);
        FILE* file = fopen(filename, "r");
        int row, col;
        fscanf(file, "%d %d", &row, &col);
        float **A = (float**)malloc((row+1)*sizeof(float*));
        for(int i = 0; i <= row; i++)
            A[i] = (float*)malloc((col+1)*sizeof(float));
        float *b = (float*)malloc((col+1)*sizeof(float));

        for(int i = 1; i <= row; i++)
            for(int j = 1; j <= col; j++)
                fscanf(file, "%f", &A[i][j]);

        for(int i = 1; i <= col; i++)
            fscanf(file, "%f\t", &b[i]);

        // Gauss-Jordan elimination
        // need to do memcpy, since the function changes the value passed by pointer
        float **tmp_A = (float**)malloc((row+1)*sizeof(float*));
        for(int i = 0; i <= row; i++)
            tmp_A[i] = (float*)malloc((col+1)*sizeof(float));
        for(int i = 1; i <= row; i++)
            for(int j = 1; j <= col; j++)
                tmp_A[i][j] = A[i][j];
        float **B = (float**)malloc((row+1)*sizeof(float*));
        for(int i = 0; i <= row; i++)
            B[i] = (float*)malloc((2)*sizeof(float));
        for(int i = 1; i <= row; i++)
            B[i][1] = b[i];
        // do the Gauss-Jordan elimination
        if(gaussj(tmp_A, row, B, 1) < 0)
            printf("\nGauss-Joran Elimination: Singular Matrix");
        else{
            printf("\nInverse Matrix:\n");
            for(int i = 1; i <= row; i++){
                for(int j = 1; j <= col; j++){
                    printf("%f\t", tmp_A[i][j]);
                }
                printf("\n");
            }
            printf("\nGauss-Joran elimination:\n");
            for(int i = 1; i <= row; i++)
                printf("%f\t", B[i][1]);
        }

        // LU Decomposition
        // ludcmp function uses Crout's method
        printf("\nLU Decomposition:\n");
        for(int i = 1; i <= row; i++)
            for(int j = 1; j <= col; j++)
                tmp_A[i][j] = A[i][j];
        int *idx = (int*)malloc((row+1)*sizeof(int));
        float *tmp_b = (float*)malloc((row+1)*sizeof(float));
        for(int i = 1; i <= row; i++)
            tmp_b[i] = b[i];
        // do LU decomposition
        // ludcmp returns L and U A, and Crout's method set
        // U's diagonal A into 1
        float tmp_d;
        ludcmp(tmp_A, row, idx, &tmp_d);
        float det = 1.0;
        for(int i = 1; i <= row; i++)
            det *= tmp_A[i][i];
        lubksb(tmp_A, row, idx, tmp_b);
        for(int i = 1; i <= row; i++)
            printf("%f\t", tmp_b[i]);

        // Singular Value Decomposition
        printf("\nSingular Value Decomposition:\n");
        for(int i = 1; i <= row; i++)
            for(int j = 1; j <= col; j++)
                tmp_A[i][j] = A[i][j];
        float *w = (float*)malloc((row+1)*sizeof(float));
        float **V = (float**)malloc((row+1)*sizeof(float*));
        for(int i = 1; i <= row; i++)
            tmp_b[i] = b[i];
        float *x = (float*)malloc((row+1)*sizeof(float));
        for(int i = 0; i <= row; i++)
            V[i] = (float*)malloc((col+1)*sizeof(float));
        svdcmp(tmp_A, row, col, w, V);
        svbksb(tmp_A, w, V, row, col, tmp_b, x);
        for(int i = 1; i <= row; i++)
            printf("%f\t", x[i]);
        printf("\n");

        printf("determinant: %f\n", det);

        // mprove
        printf("\nmprove with LU Decomposition:\n");
        for(int i = 1; i <= row; i++)
            for(int j = 1; j <= col; j++)
                tmp_A[i][j] = A[i][j];
        for(int i = 1; i <= row; i++)
            tmp_b[i] = b[i];
        ludcmp(tmp_A, row, idx, &tmp_d);
        lubksb(tmp_A, row, idx, tmp_b);
        mprove(A, tmp_A, row, idx, b, tmp_b);
        for(int i = 1; i <= row; i++)
            printf("%f\t", tmp_b[i]);
        printf("\n");

        // for(int i = 0; i <= row; i++){
        //     free(tmp_A[i]);
        //     free(V[i]);
        //     free(A[i]);
        // }
        // free(tmp_A);
        // free(tmp_b);
        // free(V);
        // free(w);
        // free(x);
        // free(A);
        // free(b);
    }
    return 0;
}