2
\$\begingroup\$

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);
}
\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.