#include <cmath>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <cassert>

int gauss(double *a, int m, int n, double eps = 1e-12);
bool readMatrix(std::istream &s, int &m, int &n, double* &a);
bool writeMatrix(std::ostream &s, int m, int n, const double *a);

using namespace std;

int main() {
    int m, n;
    double *a;
    ifstream input("input.txt");
    if (!input.is_open()) {
        cerr << "Cannot open an input file" << endl;
        return (-1);
    }
    if (!readMatrix(input, m, n, a)) {
        cerr << "Illegal format of matrix" << endl;
        return (-1);
    }
    cout << "Initial matrix:" << endl;
    writeMatrix(cout, m, n, a);
    
    int r = gauss(a, m, n);
    ofstream output("output.txt");
    if (!output.is_open()) {
        cerr << "Cannot open an output file" << endl;
        return (-1);
    }
    writeMatrix(output, m, n, a);
    cout << "Rank of matrix: " << r << endl;
    cout << "Echelonized matrix:" << endl;
    writeMatrix(cout, m, n, a);
    
    delete[] a;
    return 0;
}

bool readMatrix(istream &s, int &m, int &n, double* &a) {
    a = nullptr;
    if (!(s >> m >> n))
        return false;
    a = new double[m*n];
    for (int i = 0; i < m; ++i) {
        for (int j = 0; j < n; ++j) {
            if (!(s >> a[i*n + j]))
                return false;
        }
    }        
    return true;
}

bool writeMatrix(ostream &s, int m, int n, const double *a) {
    if (!(s << m << " " << n << endl))
        return false;
    s << fixed << setprecision(6);
    for (int i = 0; i < m; ++i) {
        for (int j = 0; j < n; ++j) {
            if (j > 0)
                if (!(s << " "))
                    return false;
            if (!(s  << setw(12) << a[i*n + j]))
                return false;
        }
        if (!(s << endl))
            return false;
    }        
    return true;    
}

int gauss(double *a, int m, int n, double eps) {
    int i = 0;
    int j = 0;
    while (i < m && j < n) {
        // Find the max. element in j-th column, 
        // starting from i-th row
        int imax = i;
        double amax = fabs(a[i*n + j]);
        for (int k = i+1; k < m; ++k) {
            if (fabs(a[k*n + j]) > amax) {
                amax = fabs(a[k*n + j]);
                imax = k;
            }
        }
        if (amax <= eps) {
            ++j; continue;
        }        
        
        assert(fabs(a[imax*n + j]) > eps);
        if (i != imax) {
            // Swap rows i and imax and change the sign 
            // of one row (to preserve the determinant)
            for (int l = j; l < n; ++l) {
                double tmp = a[i*n + l];
                a[i*n + l] = a[imax*n + l];
                a[imax*n + l] = (-tmp); // change the sign
            }
        }
        
        assert(fabs(a[i*n + j]) > eps);
        double r = a[i*n + j];  // resolving element
        
        // Annihilate the j-th column, starting 
        // from the (i+1)-th row
        for (int k = i+1; k < m; ++k) {
            double c = a[k*n + j]/r;
            a[k*n + j] = 0.;
            for (int l = j+1; l < n; ++l) {
                a[k*n + l] -= c*a[i*n + l];
            }
        }
        
        ++i; ++j;
    }
    return i;    
}