// class Matrix that uses the standar library STL
#ifndef MATRIX_H
#define MATRIX_H 1

#include <iostream>
#include <stdexcept>
#include <cassert>
#include <vector>

class Matrix {
private:
    int m;  // Number of rows
    int n;  // Number of columns
    std::vector<double> a;
public:
    Matrix(
        int rows = 1, int cols = 1, 
        const double *elems = nullptr
    ):
        m(rows),
        n(cols),
        a(m*n)
    {
        if (elems != nullptr) {
            for (int i = 0; i < m*n; ++i) {
                a[i] = elems[i];
            }
        }
    }
    
    Matrix(
        int rows, int cols, 
        const std::vector<double>& elems
    ):
        m(rows),
        n(cols),
        a(m*n)
    {
        for (int i = 0; i < m*n; ++i) {
            a[i] = elems[i];
        }
    }
    
    int rows() const {
        return m;
    }
    
    int cols() const {
        return n;
    }
    
    void setDimensions(
        int newRows, int newCols
    ) {
        m = newRows;
        n = newCols;
        a.resize(m*n);
    }
    
    // x = a[i][j];
    const double* operator[](int i) const {
        return &(a[i*n]);
    }
      
    // a[i][j] = y;
    double* operator[](int i) {
        return &(a[i*n]);
    }
    
    const double& at(int i, int j) const {
        if (
            i < 0 || i >= m ||
            j < 0 || j >= n
        )
            throw std::out_of_range(
                "Incorrect matrix indices"
            ); 
        return a[i*n + j];
    }
    
    double& at(int i, int j) {
        if (
            i < 0 || i >= m ||
            j < 0 || j >= n
        )
            throw std::out_of_range(
                "Incorrect matrix indices"
            ); 
        return a[i*n + j];
    }
    
    // Transform the matrix to the row echelon form
    // using elementary transforms of matrix rows.
    // Return the rank of matrix.
    int gaussElimination(double eps = 1e-12);
    double det() const;
    
private:
    // Transform the row echelon form of matrix
    // to the reduced form
    void reversePass(int r, double eps = 1e-12);
public:
    // Compute reduced row echelon form of matrix
    int rref(Matrix& res, double eps = 1e-12) const; 
    
    // Solve the linear system
    // A*x = b, where matrix A = *this
    // Return false for singular matrix.
    bool solveLinearSystem(
        const double* b,
        double* x,
        double eps = 1e-12
    ) const;
    
    bool inverseMatrix(
        Matrix& res,
        double eps = 1e-12
    ) const; 

    // Elementary transfors of matrix rows
    // Swap two matrix rows and multiply one of them
    // by -1 so that determinant be invariant.
    void swapRows(int i, int k);
    
    // Add k-th row of matrix multiplied by the coeff. c
    // to i-th row
    void addRows(int i, int k, double c);
    
    // Multiply the i-th row of matrix by the coeff. c
    void multiplyRow(int i, double c);
};

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

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

#endif