Example: 3x3 Matrix Inverse in C++

// 3x3 Matrix Inverse
#include <iostream>
#include <vector>
#include <iomanip> // For std::fixed and std::setprecision

// Function to calculate the determinant of a 2x2 matrix
double determinant2x2(double a, double b, double c, double d) {
    return a * d - b * c;
}

// Function to calculate the determinant of a 3x3 matrix
double determinant3x3(const std::vector<std::vector<double>>& matrix) {
    double det = 0;
    det += matrix[0][0] * determinant2x2(matrix[1][1], matrix[1][2], matrix[2][1], matrix[2][2]);
    det -= matrix[0][1] * determinant2x2(matrix[1][0], matrix[1][2], matrix[2][0], matrix[2][2]);
    det += matrix[0][2] * determinant2x2(matrix[1][0], matrix[1][1], matrix[2][0], matrix[2][1]);
    return det;
}

// Function to calculate the adjoint of a 3x3 matrix
std::vector<std::vector<double>> adjoint3x3(const std::vector<std::vector<double>>& matrix) {
    std::vector<std::vector<double>> adj(3, std::vector<double>(3));
    
    // Calculate cofactor matrix and then transpose it to get adjoint
    adj[0][0] = determinant2x2(matrix[1][1], matrix[1][2], matrix[2][1], matrix[2][2]);
    adj[0][1] = -determinant2x2(matrix[0][1], matrix[0][2], matrix[2][1], matrix[2][2]); // Transposed from (1,0) cofactor
    adj[0][2] = determinant2x2(matrix[0][1], matrix[0][2], matrix[1][1], matrix[1][2]);  // Transposed from (2,0) cofactor

    adj[1][0] = -determinant2x2(matrix[1][0], matrix[1][2], matrix[2][0], matrix[2][2]); // Transposed from (0,1) cofactor
    adj[1][1] = determinant2x2(matrix[0][0], matrix[0][2], matrix[2][0], matrix[2][2]);
    adj[1][2] = -determinant2x2(matrix[0][0], matrix[0][2], matrix[1][0], matrix[1][2]); // Transposed from (2,1) cofactor

    adj[2][0] = determinant2x2(matrix[1][0], matrix[1][1], matrix[2][0], matrix[2][1]);  // Transposed from (0,2) cofactor
    adj[2][1] = -determinant2x2(matrix[0][0], matrix[0][1], matrix[2][0], matrix[2][1]); // Transposed from (1,2) cofactor
    adj[2][2] = determinant2x2(matrix[0][0], matrix[0][1], matrix[1][0], matrix[1][1]);
    
    return adj;
}

// Function to find the inverse of a 3x3 matrix
std::vector<std::vector<double>> inverse3x3(const std::vector<std::vector<double>>& matrix) {
    double det = determinant3x3(matrix);
    if (det == 0) {
        std::cerr << "Determinant is zero, inverse does not exist." << std::endl;
        // Return an empty or zero matrix to signify failure
        return std::vector<std::vector<double>>(); 
    }

    std::vector<std::vector<double>> adj = adjoint3x3(matrix);
    std::vector<std::vector<double>> inv(3, std::vector<double>(3));
    
    for (int i = 0; i < 3; ++i) {
        for (int j = 0; j < 3; ++j) {
            inv[i][j] = adj[i][j] / det;
        }
    }
    return inv;
}

// Function to print a 3x3 matrix
void printMatrix(const std::vector<std::vector<double>>& matrix) {
    if (matrix.empty()) {
        std::cout << "Matrix is empty or invalid." << std::endl;
        return;
    }
    for (int i = 0; i < 3; ++i) {
        for (int j = 0; j < 3; ++j) {
            std::cout << std::fixed << std::setprecision(4) << matrix[i][j] << "\t";
        }
        std::cout << std::endl;
    }
}

int main() {
    // Step 1: Define the 3x3 matrix
    std::vector<std::vector<double>> A = {
        {1, 2, 3},
        {0, 1, 4},
        {5, 6, 0}
    };

    std::cout << "Original Matrix A:" << std::endl;
    printMatrix(A);
    std::cout << std::endl;

    // Step 2: Calculate the determinant
    double det_A = determinant3x3(A);
    std::cout << "Determinant of A: " << det_A << std::endl;
    std::cout << std::endl;

    // Step 3: Find the inverse matrix
    if (det_A != 0) {
        std::vector<std::vector<double>> A_inv = inverse3x3(A);
        std::cout << "Inverse of Matrix A (A^-1):" << std::endl;
        printMatrix(A_inv);
        std::cout << std::endl;

        // Optional: Verify A * A^-1 = I
        std::cout << "Verification (A * A^-1):" << std::endl;
        std::vector<std::vector<double>> product(3, std::vector<double>(3, 0.0));
        for (int i = 0; i < 3; ++i) {
            for (int j = 0; j < 3; ++j) {
                for (int k = 0; k < 3; ++k) {
                    product[i][j] += A[i][k] * A_inv[k][j];
                }
            }
        }
        printMatrix(product);
    } else {
        std::cout << "Cannot compute inverse as determinant is zero." << std::endl;
    }

    return 0;
}