A row reduction algorithm implement in C++

matrix.cpp

#include <iostream>
#include <vector>
#include <fstream>

#define DATA_TYPE float

using namespace std;

class Position
{
  public:
    int row;
    int col;
    Position(int row, int col)
    {
        this->row = row;
        this->col = col;
    }
    int GetValue(vector<vector<DATA_TYPE> > matrix)
    {
        return matrix[row][col];
    }
};
void dump_matrix(vector<vector<DATA_TYPE> > matrix)
{
    for (vector<vector<DATA_TYPE> >::iterator i = matrix.begin(); i < matrix.end(); i++)
    {
        for (vector<DATA_TYPE>::iterator j = i->begin(); j < i->end(); j++)
        {
            cout << *j << " ";
        }
        cout << endl;
    }
}
void row_scale(vector<vector<DATA_TYPE> >::iterator row, DATA_TYPE ratio)
{
    for (vector<DATA_TYPE>::iterator i = row->begin(); i < row->end(); i++)
    {
        *i *= ratio;
    }
}
void row_replace(vector<vector<DATA_TYPE> >::iterator base,
                 vector<vector<DATA_TYPE> >::iterator op, DATA_TYPE ratio)
{
    for (vector<DATA_TYPE>::iterator i = base->begin(); i < base->end(); i++)
    {
        *i += *(op->begin() + (i - base->begin())) * ratio;
    }
}
void row_exchange(vector<vector<DATA_TYPE> >::iterator a,
                  vector<vector<DATA_TYPE> >::iterator b)
{
    for (vector<DATA_TYPE>::iterator i = a->begin(); i < a->end(); i++)
    {

        float temp = *(b->begin() + (i - a->begin()));
        *(b->begin() + (i - a->begin())) = *i;
        *i = temp;
    }
}
int is_nonzero_column(vector<vector<DATA_TYPE> > matrix, int column_id, int rows, int next_row_id)
{
    for (int row = next_row_id; row < rows; row++)
    {
        if (matrix[row][column_id] != 0)
            return row;
    }
    return -1;
}

int main()
{
    vector<vector<DATA_TYPE> > matrix;
    int m, n;
    cin >> m >> n;
    for (int i = 0; i < m; i++)
    {
        vector<DATA_TYPE> row;
        for (int j = 0; j < n; j++)
        {
            DATA_TYPE a;
            cin >> a;
            row.push_back(a);
        }
        matrix.push_back(row);
    }
    int next_row_id = 0;
    vector<Position> pivot_positions;
    for (int col = 0; col < n; col++)
    {
        // check if this is a nonzero column
        int nonzero_row_id = is_nonzero_column(matrix, col, m, next_row_id);
        if (nonzero_row_id >= 0)
        {
            cout << "discovered leftmost nonzero column " << col
                 << ", and topmost nonzero row " << nonzero_row_id << endl;
            if (nonzero_row_id != next_row_id)
            {
                cout << "exchanged row " << next_row_id << " and row " << nonzero_row_id << endl;
                row_exchange(matrix.begin() + next_row_id, matrix.begin() + nonzero_row_id);
                nonzero_row_id = next_row_id;
            }
            cout << "selected pivot at position (row=" << nonzero_row_id << ",col=" << col << ") with value "
                 << matrix[nonzero_row_id][col] << endl;
            pivot_positions.push_back(Position(nonzero_row_id, col));
            for (int row = next_row_id; row < m; row++)
            {
                if (matrix[row][col] == 0)
                    continue;
                if (row == nonzero_row_id)
                    continue;
                cout << "adding " << -matrix[row][col] / matrix[nonzero_row_id][col]
                     << " times of row " << nonzero_row_id << " onto row " << row << endl;
                row_replace(matrix.begin() + row,
                            matrix.begin() + nonzero_row_id,
                            -matrix[row][col] / matrix[nonzero_row_id][col]);
            }
            next_row_id++;
        }
        else
        {
            continue;
        }
    }
    // now it comes to echelon form
    // and we start to make it into reduced echelon form
    for (vector<Position>::iterator pos = pivot_positions.end() - 1;
         pos >= pivot_positions.begin(); pos--)
    {
        if (pos->GetValue(matrix) != 1)
        {
            cout << "scaling row " << pos->row << " by "
                 << 1 / matrix[pos->row][pos->col] << endl;
            row_scale(matrix.begin() + pos->row, 1 / matrix[pos->row][pos->col]);
        }
        for (int row = 0; row < m; row++)
        {
            if (matrix[row][pos->col] != 0 && row != pos->row)
            {
                cout << "adding " << -matrix[row][pos->col] / matrix[pos->row][pos->col]
                     << " times of row " << pos->row << " onto row " << row << endl;
                row_replace(matrix.begin() + row, matrix.begin() + pos->row,
                            -matrix[row][pos->col] / matrix[pos->row][pos->col]);
            }
        }
    }
    dump_matrix(matrix);
    return 0;
}