# Linear System Solver using Gauss-Jordan Elimination

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;
}
}

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; }

}
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);
}