#include <stdio.h>
#include <stdlib.h>
#include <string.h>
 
typedef struct matrix_s{
    char *name;
    double **element;
    long columns;
    long rows;
}matrix_t;
 
typedef enum boolean_e{
    FALSE = 0,
    TRUE = 1
} boolean_t;
 
matrix_t *matrix_allocate(long columns, long rows, char *name)
{
    long i,j;
    matrix_t *m=calloc(1,sizeof(matrix_t));
 
    m->name=malloc(strlen(name));
    strcpy(m->name,name);
 
    m->rows = rows;
    m->columns = columns;
 
    m->element=calloc(rows,sizeof(double*));
    if(!m->name || !m->element)
    {
        fprintf(stderr,"Alloc error!");
        free(m);
        return 0;
    }
 
    for(i=0;i<rows;i++)
    {
        m->element[i]=calloc(columns,sizeof(double));
        if(!m->element[i])
        {
            fprintf(stderr,"alloc error");
            for(j=0;i<i;j++)
                free(m->element[j]);    //bisherigen wieder freigeben
            free(m->element);
            free(m->name);
            free(m);
            return 0;
        }
    }
    return m;
}
 
boolean_t matrix_insert(matrix_t *m,long columns,long rows,double val)
{
    boolean_t status=FALSE;
    if(!m) return status;
 
    if(columns<0 || rows<0 || columns>=m->columns || rows>=m->rows) return status;
 
    m->element[rows][columns] = val;
 
    return status=TRUE;
}
 
boolean_t matrix_free(matrix_t *m)
{
    boolean_t status=FALSE;
    if(!m) return status;
 
    long i;
    for(i=0;i<m->rows;i++)
        free(m->element[i]);
    free(m->element);
    if(m->name) free(m->name);
    free(m);
 
    return status=TRUE;
}
 
void matrix_print(matrix_t *m)
{
    long i,j;
 
    if(!m) return;
 
    for(i=0;i<m->rows;i++)
    {
        for(j=0;j<m->columns;j++)
            printf("%5g",m->element[i][j]);
        printf("\n");
    }
    printf("----------------\n");
}
 
void test_matrix(matrix_t *m)
{
    long i,j;
    if(!m)
    {
        fprintf(stderr, "Nullpointer\n");
        return;
    }
    for(i=0;i<m->rows;i++)
    {
        for(j=0;j<m->columns;j++)
            m->element[i][j]=i+j;
    }
    matrix_print(m);
}
 
boolean_t gauss_forward(matrix_t *m,matrix_t *rs)
{
    boolean_t status=FALSE;
    long i,j,k;
 
    if(!m || !rs)
    {
        fprintf(stderr, "Nullpointer\n");
        return status;
    }
 
    for(i=0;i<=m->columns-2;i++)
    {
        for(j=i+1;j<=m->rows-1;j++)
        {
            for(k=i+1;k<=m->columns-1;k++)
            {
                if((m->element[i][i])==0)
                {
                    fprintf(stderr,"DIV 0");
                    return 0;
                }
                m->element[j][k]=m->element[j][k]-(m->element[j][i]/m->element[i][i])*m->element[i][k];
            }
            if((m->element[i][i])==0)
            {
                    fprintf(stderr,"DIV 0");
                    return 0;
            }
            rs->element[j][0]=rs->element[j][0]-(m->element[j][i]/m->element[i][i])*rs->element[i][0];
            for(k=0;k<=i;k++)
                m->element[j][k]=0;
        }
    }
    return status=TRUE;
}
 
boolean_t gauss_backward(matrix_t *m,matrix_t *rs,matrix_t *sol)
{
    long i,j;
    double sum=0;
    boolean_t status=FALSE;
 
    if(!m || !rs || !sol)
    {
        fprintf(stderr, "Nullpointer\n");
        return status;
    }
 
    for(i=m->rows-1;i>=0;i--)
    {
        sum=0;
        for(j=i;j<=m->columns-1;j++)
        {
            sum=sum+(m->element[i][j]*sol->element[j][0]);
        }
        if((m->element[i][i])==0)
        {
                    fprintf(stderr,"DIV 0");
                    return 0;
        }
        sol->element[i][0] = (rs->element[i][0]-sum)/m->element[i][i];
    }
    return status=TRUE;
}
 
boolean_t gauss_solve(matrix_t *m, matrix_t *rs, matrix_t *sol)
{
    boolean_t status=FALSE;
 
    if(!m || !rs || !sol)
    {
        fprintf(stderr, "Nullpointer\n");
        return status;
    }
    if(m->columns != m->rows)
    {
        fprintf(stderr, "Wrong Dimension\n");
        return status;
    }
    if(!gauss_forward(m, rs))
    {
        fprintf(stderr, "Error in gauss_forward\n");
        return status;
    }
    if(!gauss_backward(m, rs, sol))
    {
        fprintf(stderr, "Error in gauss_backward\n");
        return status;
    }
    matrix_print(sol);
    return status=TRUE;
}
 
boolean_t matrix_fill(matrix_t *m,long dim)
{
    long c=0,val,i,j,k=0;
    FILE *fp;
    char x[10];
 
    printf("stdin (1) or file (2) ?\n");
    scanf("%ld",&c);
    fflush(stdin);
 
    switch(c)
    {
        case 1 :   for(i=0;i<dim;i++)
                        for(j=0;j<dim;j++)
                        {
                            printf("%ld.%ld.: ",i,j);
                            scanf("%ld",&val);
                            fflush(stdin);
                            matrix_insert(m,j,i,val);
                        }
                    break;
 
        case 2 :    if((fp=fopen("test.csv","r")))  //Trennung mit ;
                    {
                        for(i=0;i<dim;i++)
                            for(j=0;j<dim;j++)
                            {
                                while(x[k-1]!=';')
                                {
                                    x[k]=fgetc(fp);
                                    k++;
                                }
                                k=0;
                                matrix_insert(m,j,i,atof(x));
                            }
                        fclose(fp);
                    }
                    break;
    }
    return 1;
}
 
int main()
{
    long c=0, i, dim;
    double val=0;
    matrix_t *m=0, *rs=0, *sol=0;
    do
    {
        printf("1. Create + Fill\n");
        printf("2. Print\n");
        printf("3. Gauss Solve\n");
        printf("0. Exit\n");
 
        scanf("%ld", &c);
        fflush(stdin);
        switch(c)
        {
            case 1:
                printf("Dimension? ");
                scanf("%ld",&dim);
                fflush(stdin);
                m=matrix_allocate(dim,dim,"Matrix1");
                matrix_fill(m,dim);
                break;
            case 2:
                matrix_print(m);
                break;
            case 3:
                rs=matrix_allocate(1, dim, "Right Side");
                sol=matrix_allocate(1, dim, "Solve");
                for(i=0;i<dim;i++)
                {
                    printf("%d.%ld: ", 1, i+1);
                    scanf("%lf", &val);
                    rs->element[i][0]=val;
                }
                matrix_print(m);
                printf("Right Side:\n");
                matrix_print(rs);
                gauss_solve(m,rs,sol);
                matrix_free(m);
                matrix_free(rs);
                matrix_free(sol);
                break;
        }
    } while(c!=0);
    return 0;
}