user12493956
user12493956

Reputation:

C++ - Calculating the inverse of a matrix

I've tried to write a program that should be able to calculate the inverse of a matrix: Here's what I have so far:

#include <iostream>
#include <vector>
#include <math.h>
#include <iomanip>
#include <stdexcept>

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

    if(dimension == 0) {
        return 1;
    }

    if(dimension == 1) {
        return vect[0][0];
    }

    //Formula for 2x2-matrix
    if(dimension == 2) {
        return vect[0][0] * vect[1][1] - vect[0][1] * vect[1][0];
    }

    double result = 0;
    int sign = 1;
    for(int i = 0; i < dimension; i++) {

        //Submatrix
        std::vector<std::vector<double>> subVect(dimension - 1, std::vector<double> (dimension - 1));
        for(int m = 1; m < dimension; m++) {
            int z = 0;
            for(int n = 0; n < dimension; n++) {
                if(n != i) {
                    subVect[m-1][z] = vect[m][n];
                    z++;
                }
            }
        }

        //recursive call
        result = result + sign * vect[0][i] * getDeterminant(subVect);
        sign = -sign;
    }

    return result;
}

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) {
    if(getDeterminant(vect) == 0) {
        throw std::runtime_error("Determinant is 0");
    } 
    double d = 1.0/getDeterminant(vect);
    std::vector<std::vector<double>> solution(vect.size(), std::vector<double> (vect.size()));

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

    return getTranspose(getCofactor(solution));
}

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

int main() {
    std::vector<std::vector<double>> matrix(3, std::vector<double> (3));
    matrix = {
        {1,2,3},
        {4,5,6},
        {7,8,8}
    };

    printMatrix(getInverse(matrix));
    return 0;
}

The functions for calculating the determinant, the transpose- and the cofactor-matrix work correctly (as far as I can see), but the function for calculating the inverse-matrix doesn't. I searched the internet and found this, which uses the same function for calculating the inverse.

Is this formula incorrect, or do you have any other idea, why it doesnt work?


The matrix I am using is

and the inverse of it should be

Upvotes: 3

Views: 21904

Answers (2)

Jaden Park
Jaden Park

Reputation: 5

I know this is an old question but your code does not work when the input matrix's dimension is 1. Here is a workaround that I used:

vector<vector<double>> inverse(const vector<vector<double>> A) {
    double d = 1.0/det(A);
    vector<vector<double>> solution(A.size(), vector<double> (A.size()));

    if(A.size() == 1){
        vector<double> ans = {0};
        ans[0] = 1.0/det(A);
        solution[0] = ans;
        return solution;
    }

    for(size_t i = 0; i < A.size(); i++) {
        for(size_t j = 0; j < A.size(); j++) {
            solution[i][j] = A[i][j] * d; 
        }
    }

    return transpose(cofactor(solution));
}

Upvotes: 0

user12493956
user12493956

Reputation:

First of all, thanks for your comments.

The problem was the order of execution. The correct solution is:

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

    double d = 1.0/getDeterminant(vect);
    std::vector<std::vector<double>> solution(vect.size(), std::vector<double> (vect.size()));

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

    solution = getTranspose(getCofactor(solution));

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

    return solution;
}

Upvotes: 3

Related Questions