#include <cmath>
#include <cassert>
#include "Matrix.h"

// Gauss elimination: transform the matrix to the
// row echelon form using elementary operations with rows
// Return the rank of matrix
int Matrix::gauss(double eps) {
    int i = 0;
    int j = 0;
    while (i < m && j < n) {
        // 1. Search for maximal element in j-th column
        int imax = i;
        double amax = fabs((*this)[i][j]);
        for (int k = i+1; k < m; ++k) {
            if (fabs((*this)[k][j]) > amax) {
                imax = k;
                amax = fabs((*this)[k][j]);
            }
        }
        if (amax <= eps) {
            ++j; continue;
        }

        if (imax != i) {
            // Swap rows i and imax
            for (int l = j; l < n; ++l) {
                double tmp = (*this)[i][l];
                (*this)[i][l] = (*this)[imax][l];
                (*this)[imax][l] = (-tmp);
            }
        }
        
        // Annihilate j-th column starting
        //     from i+1-th row
        for (int k = i+1; k < m; ++k) {
            double coeff = (*this)[k][j]/(*this)[i][j];
            for (int l = j; l < n; ++l) {
                (*this)[k][l] -= (*this)[i][l]*coeff;
            }
        }
        ++i;
        ++j;
    }
    return i;
}

// Transform the matrix to the reduced row echelon form:
// first non-zero elements of rows equal 1, and
// the elements above them equal 0. 
int Matrix::rref(double eps) {


} 

std::ostream& operator<<(std::ostream& s, const Matrix& a) {
    for (int i = 0; i < a.rows(); ++i) {
        for (int j = 0; j < a.cols(); ++j) {
            if (j > 0)
                s << " ";   // Delimiter
            s << a[i][j];
        }
        s << std::endl;
    }
    return s;
}

std::istream& operator>>(std::istream& s, Matrix& a) {
    for (int i = 0; s.good() && i < a.rows(); ++i) {
        for (int j = 0; s.good() && j < a.cols(); ++j) {
            s >> a[i][j];
        }
    }
    return s;
}