Matrix inversion via Gaussian Elimination

Question

I'm focusing on 3x3 matrices for now as my code is volatile. I read the matrix from a text file and print to the console, based on its dimensions I generate the identity matrix.

const int m = 3;
const int n = 3;

int ID[m][n] = {};

for (i = 1; i <= n; ++i){

        ID[i][i] = 1;
}

For some reason ID(2)[3] gets printed as 4227276 so I have to force it to zero manually after the fact.

Aside from other elementary row operations like swapping rows based on leading entry position, the main chunk of my code consists of the following:

float matrix[m][n];

int i,j,k,p,s;

for(s = 1;s <= m;++s){
        j = s;
        k = j + 1;
        p = j;

        for(i = n;i >= j;--i){                    // makes leading entries 1

                ID[j][i] = ID[j][i]/matrix[j][j];
                matrix[j][i] = matrix[j][i]/matrix[j][j];

            }

            for(j = k;j <= m;++j){             //converts to upper triangular

                for(i = n;i >= 1;--i){

                    ID[j][i] = ID[j][i] - matrix[j][i]*matrix[p][i];
                    matrix[j][i] = matrix[j][i] - matrix[j][i]*matrix[p][i];

                }

            }
    }
            for(j = (m-1);j >= 1;--j){    //makes entries above diagonal zero

                for(i = n;i > j;--i){

                    ID[j][i] = ID[j][i] - matrix[j][i]*matrix[i][i];
                    matrix[j][i] = matrix[j][i] - matrix[j][i]*matrix[i][i];

                }

            }

I'm basically doing to the identity matrix whatever I do to matrix[m][n] to reduce it to row echelon form as you would with the augmented matrix. The row operations are pretty haphazard as I was just doing whatever worked to make matrix[m][n] an identity matrix. Afterwards, I just slotted ID[m][n] in there... not really sure what's happening but the result is half right.

my result right answer

I realize that I the term I subtract from ID might need to be a multiple of ID but that makes it even worse. What mistakes have I made?

`for (i = 1; i <= n; ++i)` arrays in C++ are 0 based. So all your array accesses (and loops) ignore the 1st element and go past the end of the array resulting in Undefined behaviour. — Richard Critten, Apr 03 '20 at 10:00
use `std::array,3>` instead. it has a `at` method to access elements that would have told you what is wrong. Once you are completely certain to access valid indices you can go back to the cheaper `operator[]` access — 463035818_is_not_an_ai, Apr 03 '20 at 10:09
not really a duplicate, but the answers might help to understand what is going on here: https://stackoverflow.com/questions/1239938/accessing-an-array-out-of-bounds-gives-no-error-why — 463035818_is_not_an_ai, Apr 03 '20 at 10:11

score 0 · Answer 1 · answered Apr 03 '20 at 10:47

In C++ the indexes of a n-dimensional array start from 0 to n-1: so the first element of the array a is a[0], the second element is a[1], ..., the n-th element is a[n-1].

When you use the for

for (i = 1; i <= n; ++i){
        ID[i][i] = 1;
}

you are discarding the first elements of each row and each column, accessing moreover to memory positions that do not belong to ID (e.g. ID[n][n]) which contain some unknown values.

You have to iterate over your arrays using for cycles such as

for (i = 0; i < n; ++i){
        ID[i][i] = 1;
}

or if you desire

for (i = 1; i <= n; ++i){
        ID[i-1][i-1] = 1;
}

but I found last solution quite confusing.

Matrix inversion via Gaussian Elimination

1 Answers1