C++ matrix calculator

Question

I wrote a program that is able to perform calculations with matrices:

main.cpp (Just for test purposes)

#include <iostream> 
#include <vector>
#include "matrix.hpp"

int main() { 
    std::cout << "Width (matrix 1):\n";
    int width = getDimension();

    std::cout << "Height (matrix 1):\n";
    int height = getDimension();

    std::vector<std::vector<double>> matrix(height, std::vector<double> (width));

    //Now, the user has to enter the matrix line by line, seperated by commas
    for(int i = 1; i <= height; i++) {
        getUserInput(matrix, i, width);
    }

    //Output
    printMatrix(matrix);
    std::cout << "\n";
    std::cout << "Determinant is " << getDeterminant(matrix) << "\n";
    std::cout << "\nMatrix * Matrix =\n"; 
    printMatrix(getMultiplication(matrix, matrix));
    std::cout << "\nTransposed Matrix:\n"; 
    printMatrix(getTranspose(matrix));
    std::cout << "\nCofactor-Matrix:\n"; 
    printMatrix(getCofactor(matrix));
    std::cout << "\nInverse-Matrix:\n"; 
    printMatrix(getInverse(matrix));
    return 0;
}

matrix.cpp (responsible for the actual calculations)

#include <iostream>
#include <vector>
#include <math.h>
#include <iomanip>
#include <stdexcept>
#include "matrix.hpp"

//Bareiss-Algorithm with O(n^3); not Laplace-expansion with O(n!)
double getDeterminant(std::vector<std::vector<double>> vect) {
    for(int i = 0; i < vect.size(); i++) {
        for(int j = i + 1; j < vect.size(); j++) {
            for(int k = i + 1; k < vect.size(); k++) {
                vect[j][k] = vect[j][k] * vect[i][i] - vect[j][i] * vect[i][k];
                if(i != 0) {
                    vect[j][k] /= vect[i - 1][i - 1];
                }
            }
        }
    }
    return vect[vect.size() - 1][vect.size() - 1];
}

//O(n^3) - much faster isn't possible
std::vector<std::vector<double>> getMultiplication(const std::vector<std::vector<double>> matrix1, const std::vector<std::vector<double>> matrix2) {

    if(matrix1[0].size() != matrix2.size()) {
        throw std::runtime_error("Matrices cannot be multiplicated");
    }
    int height1 = matrix1.size();

    int width2 = matrix2[0].size();
    int height2 = matrix2.size();

    std::vector<std::vector<double>> solution(height1, std::vector<double> (width2));

    //Formula for matrix-multiplication
    for(int i = 0; i < height1; i++) {
        for(int k = 0; k < width2; k++) {
            for(int j = 0; j < height2; j++) {
                solution[i][k] = matrix1[i][j] * matrix2[j][k]; 
            }
        }
    }

    return solution;
}

//O(n^2) - faster is not possible; 
std::vector<std::vector<double>> getTranspose(const std::vector<std::vector<double>> matrix1) {

    //Transpose-matrix: height = width(matrix), width = height(matrix)
    std::vector<std::vector<double>> solution(matrix1[0].size(), std::vector<double> (matrix1.size()));

    //Filling solution-matrix
    for(size_t i = 0; i < matrix1.size(); i++) {
        for(size_t j = 0; j < matrix1[0].size(); j++) {
            solution[j][i] = matrix1[i][j];
        }
    }
    return solution;
}

std::vector<std::vector<double>> getCofactor(const std::vector<std::vector<double>> vect) {
    if(vect.size() != vect[0].size()) {
        throw std::runtime_error("Matrix is not quadratic");
    } 

    std::vector<std::vector<double>> solution(vect.size(), std::vector<double> (vect.size()));
    std::vector<std::vector<double>> subVect(vect.size() - 1, std::vector<double> (vect.size() - 1));

    for(std::size_t i = 0; i < vect.size(); i++) {
        for(std::size_t j = 0; j < vect[0].size(); j++) {

            int p = 0;
            for(size_t x = 0; x < vect.size(); x++) {
                if(x == i) {
                    continue;
                }
                int q = 0;

                for(size_t y = 0; y < vect.size(); y++) {
                    if(y == j) {
                        continue;
                    }

                    subVect[p][q] = vect[x][y];
                    q++;
                }
                p++;
            }
            solution[i][j] = pow(-1, i + j) * getDeterminant(subVect);
        }
    }
    return solution;
}

std::vector<std::vector<double>> getInverse(const std::vector<std::vector<double>> vect) {
    double det = getDeterminant(vect);
    if(det == 0) {
        throw std::runtime_error("Determinant is 0");
    } 
    std::vector<std::vector<double>> solution(vect.size(), std::vector<double> (vect.size()));

    solution = getTranspose(getCofactor(vect));

    for(size_t i = 0; i < vect.size(); i++) {
        for(size_t j = 0; j < vect.size(); j++) {
            solution[i][j] *= (1/det);
        }
    }

    return solution;
}

int getDimension() {
    int dimension;
    std::cout << "Please enter dimension of Matrix: ";
    std::cin >> dimension;
    std::cout << "\n";

    if(dimension < 0 || std::cin.fail()) {
        std::cin.clear();
        std::cin.ignore();
        std::cout << "ERROR: Dimension cannot be < 0.\n";
        return getDimension();
    }

    return dimension;
}

void printMatrix(const std::vector<std::vector<double>> vect) {
    for(std::size_t i = 0; i < vect.size(); i++) {
        for(std::size_t j = 0; j < vect[0].size(); j++) {
            std::cout << std::setw(8) << vect[i][j] << " ";
        }
        std::cout << "\n";
    }
}

void getUserInput(std::vector<std::vector<double>>& vect, int i, int dimension) {
    std::string str = "";
    std::cout << "Enter line " << i << " only seperated by commas: ";
    std::cin >> str;
    std::cout << "\n";
    str = str + ',';

    std::string number = "";
    int count = 0;

    for(std::size_t k = 0; k < str.length(); k++) {
        if(str[k] != ',') {
            number = number + str[k];
        }
        else if(count < dimension) {
            if(number.find_first_not_of("0123456789.-") != std::string::npos) {
                std::cout << "ERROR: Not only numbers entered.\n";
                getUserInput(vect, i, dimension);
                break;
            }
            else if(number.find_first_not_of("0123456789-.") == std::string::npos) {
                vect[i - 1][count] = std::stod(number);
                number = "";
                count++;
            }
            else {
                std::cout << "ERROR: Not enough numbers entered.\n";
                getUserInput(vect, i, dimension);
                break;
            }       
        }
        else {
            std::cout << "ERROR: Too many numbers entered.\n";
            getUserInput(vect, i, dimension);
            break;
        }
    }    
}

matrix.hpp

#ifndef MATRIX_HPP
#define MATRIX_HPP

#include <vector>

double getDeterminant(const std::vector<std::vector<double>> vect);
std::vector<std::vector<double>> getMultiplication(const std::vector<std::vector<double>> matrix1, const std::vector<std::vector<double>> matrix2);
std::vector<std::vector<double>> getTranspose(const std::vector<std::vector<double>> matrix1);
std::vector<std::vector<double>> getCofactor(const std::vector<std::vector<double>> vect);
std::vector<std::vector<double>> getInverse(const std::vector<std::vector<double>> vect);
int getDimension();
void printMatrix(const std::vector<std::vector<double>> vect);
void getUserInput(std::vector<std::vector<double>>& vect, int i, int dimension);
#endif

My questions:

What's the runtime (O-notation) of getCofactor() and getInverse()?
Are there more efficient ways to do this tasks?
How to generally improve the code?

G. Sliepen · Accepted Answer · 2020-02-19 21:50:48Z

Create a `class Matrix`

Instead of having a vector of vectors, and have global functions that manipulate those, create a class Matrix and add member functions to it to manipulate matrices. You should probably also create overloads for arithmetic operations, so you can write things like auto matrix3 = matrix1 + matrix2;

Have a look at existing C++ matrix libraries to see what is possible.

Complexity of `getCofactor()` and `getInverse()`

In getCofactor(), you have four nested for-loops that only skip iterations if x == i or y == j, so that would make it O(N^4). However, you are calling getDeterminant() inside the second outer-most loop on a subvector of size (N-1)^2, so that would make it O(N^2 * (N-1)^3) = O(N^5).

In getInverse(), the complexity is dominated by the call to getCofactor(), so it is also O(N^5). This is really bad for large matrices, since the Gauss-Jordan method which you should learn in first year Mathematics studies is just O(N^3).

Wikipedia has a list of matrix algebra algorithm complexities that you can check to see what is possible. Note however that the best looking algorithms with weird exponents are probably hard to implement and might actually be much less efficient for small matrices than the simpler algorithms. The best algorithm for you therefore depends on what sizes of matrices your programs are going to work with.

Thanks for your reply. I just noticed that the getDeterminant-function sometimes returns "-nan". What does that mean? How can I fix it? — user214772
– user214772, Commented Feb 19, 2020 at 22:26
A nan is most likely the result of a floating point division by zero. — G. Sliepen
– G. Sliepen, Commented Feb 20, 2020 at 6:37

Stack Exchange Network

C++ matrix calculator

1 Answer 1

Create a `class Matrix`

Complexity of `getCofactor()` and `getInverse()`

You must log in to answer this question.

Hot Network Questions

C++ matrix calculator

1 Answer 1

Create a class Matrix

Complexity of getCofactor() and getInverse()

You must log in to answer this question.

Related

Hot Network Questions

Create a `class Matrix`

Complexity of `getCofactor()` and `getInverse()`