Linear System Solver using Gauss-Jordan Elimination

Question

I'm learning about Gauss-Jordan Elimination, and decided to write a program to automate the process to help solidify my understanding. Along with the algorithm itself, I wrote a Matrix class that has methods for the 2 required row operations (Adding a row to a row, and multiplying a row by a constant). My main concerns that I'd like thoughts on:

The Matrix:

• The method I'm using to set up the matrix is ugly. Each row requires a call to setRow_To, which makes for (probably) longer-than-necessary code.

• The naming scheme for some of the functions is ugly. My goal was to make them "fill-in-the-blanks" style to make them easier to comprehend, but I don't think it helped.

  • Ex: add_xRow_toRow takes a multiplier as its first argument (first _), a row number as the second argument (second _), then a second row number for its last argument.

• 2D vectors were a pain to deal with last time I tried, so I went with the "pseudo-2D" approach that simulates a 2D vector using the equation getIndexFor to get the index for a given row/col. It works, but leads to some ugly code elsewhere when getIndexFor isn't sufficient.

• In add_xRow_toRow, I put the calls to getRowStartIndexFor and getRowEndIndexFor on the same line or else it spreads out over 4 lines. Is it really neater?

The Solver:

• I'm unsure if this should even be its own class. I went in thinking that it would be more complicated, and that it would require aux functions to assist it. Turns out I was able to fit it all in a single function (its constructor).

• It came out as a long-ish procedure that isn't the prettiest to read. Is there anything I can do to help readability?

• Once a column has been "eliminated", and contains only 0s and 1s, it's possible for a 0 to be multiplied by a negative number, which yields a -0 (I didn't even know -0 was meaningful). From my testing, this doesn't seem to affect the results, but it doesn't leave you with a perfect identity matrix, which is the end goal. Should I be concerned about this? This should be an easy fix (n == -0 ? n * -1 : n), but I don't want to put it in unless it's necessary, or at least beneficial.

Matrix.h

#ifndef MATRIX_H
#define MATRIX_H

#include <vector>

typedef unsigned int Index;
typedef double CellData;
typedef std::vector<CellData> Row;

class Matrix {
    Index width;
    Index height;

    std::vector<CellData> arr;

    Index getIndexFor(Index row, Index col) const;

    Index getRowStartIndexFor(Index rowN) const;
    Index getRowEndIndexFor(Index rowN) const;

public:

    Matrix(Index rows, Index cols);

    Index getWidth() const;
    Index getHeight() const;

    void setRow_To(Index rowN, Row newRow);

    CellData getAt(Index x, Index y) const;
    void updateAt(Index x, Index y, CellData newVal);

    void divRow_By(Index rowN, double divisor);
    void add_xRow_toRow(double mult, Index sourceRowN, Index targRowN);

    void display(unsigned int spacing = 5, unsigned int decPlaces = 3) const;
};

#endif

Matrix.cpp

#include "Matrix.h"

#include <iostream>
#include <iomanip>
#include <stdexcept>

Index Matrix::getIndexFor(Index row, Index col) const {
    return row * width + col;
}

Index Matrix::getRowStartIndexFor(Index rowN) const {
    return getIndexFor(rowN, 0);
}

Index Matrix::getRowEndIndexFor(Index rowN) const {
    return getIndexFor(rowN, width - 1);
}


Matrix::Matrix(Index rows, Index cols) :
    width(cols),
    height(rows) {
    arr.resize(rows * cols);
}

Index Matrix::getWidth() const {
    return width;
}

Index Matrix::getHeight() const {
    return height;
}

void Matrix::setRow_To(Index row, std::vector<double> newRow) {
    if (newRow.size() != width) {
        throw std::invalid_argument(
            "The new row vector must be as long as the matrix is wide."
        );
    }

    for (Index col = 0; col < width; col++) {
        updateAt(row, col, newRow[col]);
    }
}

double Matrix::getAt(Index row, Index col) const {
    return arr[ getIndexFor(row, col) ];
}

void Matrix::updateAt(Index row, Index col, double newVal) {
    arr[ getIndexFor(row, col) ] = newVal;
}

// Row Operation 2
void Matrix::divRow_By(Index rowN, double divisor) {
    Index startI = getRowStartIndexFor(rowN);
    Index endI = getRowEndIndexFor(rowN);

    for (Index i = startI; i <= endI; i++) {
        arr[i] /= divisor;
    }
}

// Modified Row Operations 2 + 3 combined
void Matrix::add_xRow_toRow(double mult, Index sourceRowN, Index targRowN) {
    Index sStartI = getRowStartIndexFor(sourceRowN), sEndI = getRowEndIndexFor(sourceRowN);
    Index tStartI = getRowStartIndexFor(targRowN), tEndI = getRowEndIndexFor(targRowN);

    for (Index sRow = sStartI, tRow = tStartI; sRow <= sEndI; sRow++, tRow++) {
        double toAdd = arr[sRow] * mult;
        arr[tRow] += toAdd;
    }
}

void Matrix::display(unsigned int spacing, unsigned int decPlaces) const {
    using namespace std;
    for (int row = 0; row < height; row++) {
        for (int col = 0; col < width; col++) {
            cout << left << setw(spacing) << setprecision(decPlaces) << getAt(row, col) << " ";
        }
        cout << std::endl;
    }
}

GJESolver.h:

#ifndef GJESOLVER_H
#define GJESOLVER_H

#include "Matrix.h"

#define PRINT_STEPS false

class GJESolver {

    Matrix m;

public:

    GJESolver(Matrix);

    Matrix getSolvedMatrix() const;
};

#endif

GJESolver.cpp:

#include "GJESolver.h"

#include <iostream>


GJESolver::GJESolver(Matrix mToSolve) :
    m(mToSolve) {

    //width - 1 because we don't need to process the last column
    for (Index curCol = 0; curCol < m.getWidth() - 1; curCol++) {

        //Put a 1 along the main diagonal
        if (PRINT_STEPS) { std::cout << "Dividing row " << curCol << " by " << m.getAt(curCol, curCol) << std::endl; }

        m.divRow_By(curCol, m.getAt(curCol, curCol));

        for (Index curRow = 0; curRow < m.getHeight(); curRow++) {
            if (curRow == curCol) continue;

            CellData multiplier = -m.getAt(curRow,curCol);

            if (PRINT_STEPS) { std::cout << "Adding " << multiplier << "x row " << curCol << " to row " << curRow << std::endl << std::endl; }

            m.add_xRow_toRow(multiplier, curCol, curRow);
        }
        if (PRINT_STEPS) {
            std::cout << "\n\n";
            m.display(6, 3);
            std::cout << "\n\n";
        }

    }
}

Matrix GJESolver::getSolvedMatrix() const {
    return m;
}

main.cpp:

#include <iostream>

#include "Matrix.h"
#include "GJESolver.h"

int main(int argc, char* argv[]) {
    using namespace std;
    Matrix m(3, 4);

    m.setRow_To(0, Row { 2, 4, 10, 7 });
    m.setRow_To(1, Row { 8, 7, 9, 3 });
    m.setRow_To(2, Row { 13, 2, 6, 11 });

    m.display(6,3);
    cout << "\n\n";

    GJESolver gjesolver(m);

    Matrix solvedMatrix = gjesolver.getSolvedMatrix();

    solvedMatrix.display(10, 6);
}

Answer 1

General Observations

Based on how long ago this was asked and how active you have been on the Code Review Community you probably know most of this already.

You might want to read up on negative zero, this Stack Overflow post predates this question.

2D vectors were a pain to deal with last time I tried, so I went with the "pseudo-2D" approach that simulates a 2D vector using the equation getIndexFor to get the index for a given row/col. It works, but leads to some ugly code elsewhere when getIndexFor isn't sufficient.

I would have stuck with using real 2D vectors, especially since you created a type for a row of the matrix. The Index functions would not have been necessary, and it added complexity to the Matrix class.

    std::vector<Row> localMatrix;

I'm unsure if this should even be its own class. I went in thinking that it would be more complicated, and that it would require aux functions to assist it. Turns out I was able to fit it all in a single function (its constructor).

This is a design issue. It is possible you need to work on Object Oriented design more. Based on the complexity of the Matrix class it is quite possible you could have added a solver function to the Matrix class.

Never solve a problem in a constructor, there is only one thing to do in a constructor and that is set the value of all necessary values.

The code would perform better if it was using iterators rather than idexes in the loops.

Prefer Using Statements Over C Style typedef

Since C++11 the using statement has replaced the use of the C programming language typedef. In Matrix.h

typedef unsigned int Index;
typedef double CellData;
typedef std::vector<CellData> Row;

in modern C++ this would be written as:

using Index = unsigned int;
using CellData = double;
using Row = std::vector<CellData>;

Inconsistent Use of Index Type

Generally the Index type is used properly, but there is one function where it isn't used and should be:

void Matrix::display(unsigned int spacing, unsigned int decPlaces) const {
    using namespace std;
    for (int row = 0; row < height; row++) {
        for (int col = 0; col < width; col++) {
            cout << left << setw(spacing) << setprecision(decPlaces) << getAt(row, col) << " ";
        }
        cout << std::endl;
    }
}

Both loops in the display function should be using Index instead of int as the loop control variables.

The Index Type
Rather than using unsigned int it would be better to use std::size_t. In addition to being unsigned it can store the maximum size of a theoretically possible object of any type (including array).

Rather than definining the Index type I would just use std::size_t.

Avoid using namespace std;

If you are coding professionally you probably should get out of the habit of using the using namespace std; statement. The code will more clearly define where cout and other identifiers are coming from (std::cin, std::cout). As you start using namespaces in your code it is better to identify where each function comes from because there may be function name collisions from different namespaces. The identifiercout you may override within your own classes, and you may override the operator << in your own classes as well. This stack overflow question discusses this in more detail.