c++matrix-inverse

inverse of a 2D array (matrix)


I've made a function to calculate the inverse of a matrix, but I'm having an issue getting the correct result when I run my program.

For example, when I insert the following matrix:

A1  A2  A3
(2  1   0)
(1  2   1)
(0  1   2)

I get this matrix:

0.1875    -0.125    0.0625
-0.125      0.25    -0.125
0.0625    -0.125    0.1875

While the correct matrix solution is this:

3/4 -1/2    1/4
-1/2    1   -1/2
1/4 -1/2    3/4

Here is my code:

void inverse(double A[][MAX], int n) {

    double B[MAX][MAX];
    double det = 1.0;
    int i, j, k;

    // Initializing B to identity matrix
    for(i = 0; i < n; i++)
        for(j = 0; j < n; j++)
            B[i][j] = (i == j) ? 1 : 0;

    // Performing row operations to convert A to identity matrix
    for(k = 0; k < n; k++) {
        // Find the pivot element
        int max_row = k;
        for (i = k + 1; i < n; i++) {
            if (abs(A[i][k]) > abs(A[max_row][k])) {
                max_row = i;
            }
        }
        // Swap the current row with the row with the largest pivot element
        if (max_row != k) {
            for (int col = 0; col < n; col++) {
                swap(A[k][col], A[max_row][col]);
                swap(B[k][col], B[max_row][col]);
            }
            det *= -1;
        }
        // Divide the current row by the pivot element to make it equal to 1
        double pivot = A[k][k];
        if (pivot == 0) {
            cout << "Error: Singular matrix." << endl;
            return;
        }
        det *= pivot;
        for(j = 0; j < n; j++) {
            A[k][j] /= pivot;
            B[k][j] /= pivot;
        }
        // Use the current row to eliminate the pivot elements in the other rows
        for(i = 0; i < n; i++) {
            if(i == k) continue;
            double factor = A[i][k];
            for(j = 0; j < n; j++) {
                A[i][j] -= A[k][j] * factor;
                B[i][j] -= B[k][j] * factor;
            }
        }
    }

    // Multiplying by reciprocal of the determinant
    if(det == 0) {
        cout << "Error: Singular matrix." << endl;
        return;
    }
    det = 1/det;

    for(i = 0; i < n; i++)
        for(j = 0; j < n; j++)
            B[i][j] *= det;

    // Copying the inverse matrix to A
    for(i = 0; i < n; i++)
        for(j = 0; j < n; j++)
            A[i][j] = B[i][j];
}

Solution

  • Note that the determinant is 4 and your result is 1/4th of the intended result

    // Multiplying by reciprocal of the determinant

    You should not multiply by the reciprocal of the determinant. If you had computed the adjugate of A, in that case yes, multiplying it by the reciprocal of the determinant would give you the inverse. But you already have the inverse: when A has been turned into the identity matrix, B must be the inverse.