/**
 * 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 구하기

void 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 = 2;
    for(; three <= 3; three++){
        char filename[20];
        sprintf(filename, "lineq%d.dat", three);
        printf("\nfile %s:\n", filename);
        FILE* file = fopen(filename, "r");
        int row, col;
        fscanf(file, "%d %d", &row, &col);
        row++, col++;
        float **A = (float**)malloc(row*sizeof(float*));
        for(int i = 0; i < row; i++)
            A[i] = (float*)malloc(col*sizeof(float));
        float *b = (float*)malloc(col*sizeof(float));
        row--, col--;
        printf("\nGiven matrix A: \n");
        for(int i = 1; i <= row; i++){
            for(int j = 1; j <= col; j++){
                fscanf(file, "%f", &A[i][j]);
                printf("%.2f\t", A[i][j]);
            }
            printf("\n");
        }
        printf("Given matrix b:\n");
        for(int i = 1; i <= col; i++){
            fscanf(file, "%f\t", &b[i]);
            printf("%.2f ", b[i]);
        }
        printf("\n");
        // Gauss-Jordan elimination
        printf("\nGauss-Joran elimination:\n");
        // need to do memcpy, since the function changes the value passed by pointer
        row++, col++;
        float **tmp_A = (float**)malloc(row*sizeof(float*));
        for(int i = 0; i < row; i++){
            tmp_A[i] = (float*)malloc(col*sizeof(float));
            memcpy(tmp_A[i], A[i], sizeof(float)*col);
        }
        float **B = (float**)malloc(row*sizeof(float*));
        for(int i = 1; i <= row; i++)
            B[i] = (float*)malloc((col+1)*sizeof(float));
        memcpy(B, A, sizeof(float)*row*col);
        for(int i = 1; i <= col; i++)
            B[row - 1][i] = b[i];
        row--, col--;
        // do the Gauss-Jordan elimination
        gaussj(tmp_A, row, B, col);
        for(int i = 1; i <= row; i++)
            printf("%.2f\t", B[row][i]);
        printf("\n");


        // LU Decomposition
        // ludcmp function uses Crout's method
        printf("\nLU Decomposition:\n");
        row++, col++;
        for(int i = 0; i < row; i++)
            memcpy(tmp_A[i], A[i], sizeof(float)*col);
        int *idx = (int*)malloc(row*sizeof(int));
        float *tmp_b = (float*)malloc(col*sizeof(float));
        memcpy(tmp_b, b, sizeof(b));
        row--, col--;
        // 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);
        lubksb(tmp_A, row, idx, tmp_b);
        for(int i = 1; i <= row; i++)
            printf("%.2f\t", tmp_b[i]);
        printf("\n");

        // Singular Value Decomposition
        printf("\nSingular Value Decomposition:\n");
        row++, col++;
        for(int i = 0; i < row; i++)
            memcpy(tmp_A[i], A[i], sizeof(float)*col);
        float *w = (float*)malloc(row * sizeof(float));
        float **V = (float**)malloc(row*sizeof(float*));
        memcpy(tmp_b, b, sizeof(b));
        float *x = (float*)malloc(row * sizeof(float));
        for(int i = 0; i < row; i++)
            V[i] = (float*)malloc(col*sizeof(float));
        row--, col--;
        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("%.2f\t", x[i]);
        printf("\n");

        row++, col++;
        for(int i = 0; i < row; i++){
            free(tmp_A[i]);
            free(A[i]);
        }
        free(tmp_A);
        free(A);
    }
    return 0;
}