inverse of a 2D array (matrix)

Question

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

Are you trying to calculate the inverse of a function, or are you trying to calculate the inverse of a matrix? — Jeremy Friesner, Mar 12 '23 at 01:28
[What is a debugger and how can it help me diagnose problems?](https://stackoverflow.com/questions/25385173/what-is-a-debugger-and-how-can-it-help-me-diagnose-problems) — Jesper Juhl, Mar 12 '23 at 05:05

score 2 · Accepted Answer · answered Mar 12 '23 at 01:30

2

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.

answered Mar 12 '23 at 01:30

harold

61,398
6
86
164

Here it is the original implementation with @harold's correction: https://godbolt.org/z/rhf9869zP – alfC Mar 12 '23 at 02:23

inverse of a 2D array (matrix)

1 Answers1