// Work with matrices
// reading from file, writing to file/stdout,
// Gauss method, inverse matrix, solution of linear system
// Matrix is represented by a linear array a. Let
//     m -- number of matrix rows,
//     n -- number of matrix columns;
// then an element with indices (i, j) is an element of array a
// with index i*n + j:
//     a_ij == a[i*n + j].
//
#include <stdio.h>
#include <math.h>
#include <assert.h>

bool readMatrix(FILE* f, double **a, int* m, int* n); 
bool writeMatrix(FILE* f, const double *a, int m, int n);

int gaussMethod(double *a, int m, int n); // Return a matrix rank

void swapMatrixRows(
    double *a, int m, int n,
    int i, int k
);
void addMatrixRows(         // row[k] += row[i]*coeff 
    double *a, int m, int n,
    int k, int i, double coeff
);
void multiplyMatrixRow(
    double *a, int m, int n,
    int i, double coeff
);

double det(double *a, int m, int n);
bool solveLinearSystem(
    const double *a, int n, // matrix n*n
    const double* b,        // column of free terms
    double* x               // out: solution of the system a*x = b
);
bool inverseMatrix(
    const double *a, int n, // matrix n*n
    double* inverse         // out: inverse matrix
);

const double EPS = 1e-8;

int main() {
    int m, n;
    double *a = 0;
    
    FILE* f = fopen("input.txt", "rt");
    if (f == NULL) {
        perror("Cannot open an input file");
        return (-1);
    }
    
    if (!readMatrix(f, &a, &m, &n)) {
        printf("Incorrect format of matrix.\n");
        delete[] a;
        fclose(f);
        return (-1);
    }
    fclose(f);
        
    printf("Matrix of size %d*%d:\n", m, n);
    writeMatrix(stdout, a, m, n);
    
    double* b = new double[m*n];
    for (int i = 0; i < m*n; ++i) {
        b[i] = a[i];
    }
    
    int rank = gaussMethod(a, m, n);
    
    printf("Row echelon form of matrix:\n");
    writeMatrix(stdout, a, m, n);

    printf("Rank of matrix: %d\n", rank);
    if (m == n) {
        printf("Determinant of matrix: %f\n", det(a, m, n));
    }
    
    // Writing the resulting matrix to the file "output.txt"
    FILE* g = fopen("output.txt", "wt");
    if (g == NULL) {
        perror("Cannot open an output file");
    } else {
        if (!writeMatrix(g, a, m, n))
            perror("Could not write a matrix to the output file");
        fclose(g);
    }

    delete[] a;
    
    if (m == n && rank == m) {
        // Compute an inverse matrix
        double* inverse = new double[n*n];
        bool res = inverseMatrix(b, n, inverse);
        assert(res);
        
        printf("Inverse matrix:\n");
        writeMatrix(stdout, inverse, n, n);
        
        // Output an inverse matrix into the file "inverse.txt"
        FILE* h = fopen("inverse.txt", "wt");
        if (h == NULL) {
            perror("Cannot open an output file");
        } else {
            if (!writeMatrix(h, inverse, n, n))
                perror("Could not write a matrix to the output file");
            fclose(h);
        }
        
        double* x = new double[n];
        double* r = new double[n];
        
        printf("Enter right column of linear system:\n");
        for (int i = 0; i < n; ++i) {
            if (scanf("%lf", &(r[i])) < 1) {
                printf("Incorrect input.\n");
                // exit(-1);
                goto LExit;
            }
        }
        res = solveLinearSystem(
            b, n, r, x
        );
        assert(res);
        
        printf("Solution of the system:\n");
        for (int i = 0; i < n; ++i) {
            printf("%g ", x[i]); 
        }
        printf("\n");

        LExit: ;
        delete[] r;
        delete[] x;
        delete[] b;
        delete[] inverse;
    }

    return 0;
}

bool readMatrix(FILE* f, double **a, int* m, int* n) {
    // 2 3     
    // 11 12 13
    // 21 22 23
    int rows, cols;
    if (fscanf(f, "%d%d", &rows, &cols) < 2) {
        return false;
    }
    *m = rows; *n = cols;
    double *matrix = new double[rows*cols];
    *a = matrix;
    for (int i = 0; i < rows; ++i) {
        // Read i-th row of matrix
        for (int j = 0; j < cols; ++j) {
            if (fscanf(f, "%lf", matrix + i*cols + j) < 1) {
                // Read error
                return false;
            }
        }
    }
    return true;
}
    
bool writeMatrix(FILE* f, const double *a, int m, int n) {
    if (fprintf(f, "%d %d\n", m, n) <= 0) {
        return false;
    } 
    
    for (int i = 0; i < m; ++i) {
        for (int j = 0; j < n; ++j) {
            if (fprintf(f, "%12.4f", a[i*n + j]) <= 0)
                return false;
        }
        fprintf(f, "\n");
    }    
    return true;
}

int gaussMethod(double *a, int m, int n) { // Return a matrix rank 
    int i = 0;
    int j = 0;
    while (i < m && j < n) {
        // 1. Find the maximal element in j-th column
        int k = i;
        double maxElem = fabs(a[k*n + j]);
        for (int l = i+1; l < m; ++l) {
            if (fabs(a[l*n + j]) > maxElem) {
                maxElem = fabs(a[l*n + j]);
                k = l;
            }
        }
        
        if (maxElem <= EPS) {
            // Zero colunm
            for (int l = i; l < m; ++l) {
                a[l*n + j] = 0.;
            }
            ++j;
            continue;
        }
        
        if (k != i) {
            // Swap rows i, k
            /*
            for (int l = j; l < n; ++l) {
                double tmp = a[i*n + l];
                a[i*n + l] = a[k*n + l];
                a[k*n + l] = (-tmp);
            }
            */
            swapMatrixRows(a, m, n, i, k);
        }            
            
        for (int l = i+1; l < m; ++l) {
            // Annihilate the element a_lj
            /*
            double r = a[l*n + j]/a[i*n + j];            
            a[l*n + j] = 0.;    
            for (int s = j+1; s < n; ++s) {
                a[l*n + s] -= r*a[i*n + s];
            }
            */
            double r = (-a[l*n + j]/a[i*n + j]);
            addMatrixRows(
                a, m, n,
                l, i, r
            );
            a[l*n + j] = 0.;
        }
        ++i; ++j;
    }
    return i;
}
            
double det(double *a, int m, int n) {
    double d = 1.;
    int k = n;
    if (m < n)
        k = m;
    for (int i = 0; i < k; ++i) {
        d *= a[i*n + i];
    }
    return d;
}

void swapMatrixRows(
    double *a, int m, int n,
    int i, int k
) {
    for (int j = 0; j < n; ++j) {
        double tmp = a[i*n + j];
        a[i*n + j] = a[k*n + j];
        a[k*n + j] = (-tmp);
    }
}

void addMatrixRows(         // row[k] += row[i]*coeff 
    double *a, int m, int n,
    int k, int i, double coeff
) {
    for (int j = 0; j < n; ++j) {
        a[k*n + j] += a[i*n + j]*coeff;
    }
}

void multiplyMatrixRow(
    double *a, int m, int n,
    int i, double coeff
) {
    for (int j = 0; j < n; ++j) {
        a[i*n + j] *= coeff;
    }
}

bool inverseMatrix(
    const double *a, int n, // matrix n*n
    double* inverse         // out: inverse matrix
) {
    int m = n;
    int n2 = 2*n;
    double* b = new double[m*n2];
    for (int i = 0; i < m; ++i) {
        for (int j = 0; j < n; ++j) {
            b[i*n2 + j] = a[i*n + j];
        }
        for (int j = 0; j < n; ++j) {
            if (j == i) {
                b[i*n2 + (n+j)] = 1.;
            } else {
                b[i*n2 + (n+j)] = 0.;
            }
        }
    }
    
    printf("Extended matrix:\n");
    writeMatrix(stdout, b, m, n2);
    
    // Do this in parallel!
    gaussMethod(b, m, n2);
    // End of paralel computations
    
    if (fabs(b[(n-1)*n2 + (n-1)]) <= EPS) {
        printf("Singular matrix cannot be inverted.\n");
        delete[] b;
        return false;
    }
    
    printf("Extended matrix after Gauss:\n");
    writeMatrix(stdout, b, m, n2);
    
    // Reverse pass in Gauss method
    for (int i = m - 1; i >= 0; --i) {
        assert(fabs(b[i*n2 + i]) > EPS);
        double coeff = 1./b[i*n2 + i];
        multiplyMatrixRow(b, m, n2, i, coeff);
        assert(fabs(b[i*n2 + i] - 1.) <= EPS); // Unit diagonal element 

        // Do this in parallel!
        for (int k = 0; k < i; ++k) {
            double coeff = (-b[k*n2 + i]);
            addMatrixRows(
                b, m, n2,
                k, i, coeff
            );
        }
        // End of parallel computing
    }
    
    printf("Extended matrix after reverse pass:\n");
    writeMatrix(stdout, b, m, n2);
    
    // Output an inverse matrix
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < n; ++j) {
            inverse[i*n + j] = b[i*n2 + (n+j)];
        }
    }
    return true;
}    
        
bool solveLinearSystem(
    const double *a, int n, // matrix n*n
    const double* b,        // column of the system right side
    double* x               // Solution of the system a*x = b
) {    
    int m = n;
    int n1 = n + 1;
    double* c = new double[m*n1];
    for (int i = 0; i < m; ++i) {
        for (int j = 0; j < n; ++j) {
            c[i*n1 + j] = a[i*n + j];
        }
        c[i*n1 + n] = b[i];
    }
    
    printf("Extended matrix:\n");
    writeMatrix(stdout, c, m, n1);
    
    // Do this in parallel!
    gaussMethod(c, m, n1);
    // End of paralel computations
    
    if (fabs(c[(n-1)*n1 + (n-1)]) <= EPS) {
        printf("The system with singular matrix cannot be solved.\n");
        delete[] b;
        return false;
    }
    
    printf("Extended matrix after Gauss:\n");
    writeMatrix(stdout, c, m, n1);
    
    // Reverse pass in Gauss method
    for (int i = m - 1; i >= 0; --i) {
        assert(fabs(c[i*n1 + i]) > EPS);
        double coeff = 1./c[i*n1 + i];
        multiplyMatrixRow(c, m, n1, i, coeff);
        assert(fabs(c[i*n1 + i] - 1.) <= EPS); // Unit diagonal element 

        // Do this in parallel!
        for (int k = 0; k < i; ++k) {
            double coeff = (-c[k*n1 + i]);
            addMatrixRows(
                c, m, n1,
                k, i, coeff
            );
        }
        // End of parallel computing
    }
    
    printf("Extended matrix after reverse pass:\n");
    writeMatrix(stdout, c, m, n1);
    
    // Output the solution of the system
    for (int i = 0; i < n; ++i) {
        x[i] = c[i*n1 + n];
    }
    return true;
}