0

Hello all,

Inverting an upper (or lower) triangular matrix is a trivial algorithm, due to the nature of the matrix. I am having an issue getting a part of my upper-triangular matrix inversion function to work, and I would like to get it working soon for a personal project. From reading a verbal description about this special inverse case using back-substitution, I have written some C++ code. I give it an input example 6x6 matrix, it inverts the diagonal correctly, but produces incorrect results above the diagonal. I am checking the true answer using the 'Online Matrix Calculator'. The code is as follows. If you are familiar with this technique, do you have any advice.


Thank-you,
Sean Michnowski

// Upper-Triangular Matrix Inversion Test
double R[36] = {4., 8., 2., 9., 0., 3., 0., 8., 1., 8., 4., 8., 0., 0., 5., 1., 2., 7.,
                0., 0., 0., 1., 4., 5., 0., 0., 0., 0., 6., 4., 0., 0., 0., 0., 0., 2.};
double Rinv[36];
int i, j, k;
//
for(i=0; i<36; i++)
    Rinv[i] = 0.0;
//
for(j=0; j<6; j++)
{
    Rinv[6*j+j] = 1./R[6*j+j];
    //
    for(i=0; i<(j-1); i++)
    {
        for(k=0; k<(j-1); k++)
        {
            Rinv[6*i+j] += Rinv[6*i+k]*R[6*k+j];
        }
    }
    //
    for(k=0; k<(j-1); k++)
    {
        Rinv[6*k+j] /= -R[6*j+j];
    }
}
//
displayMatrix("Matrix:", Rinv, 6, 6);
//
return 0;

The expected output is:
0.250 -0.250 -0.050 -0.200 0.317 0.667
0.000 0.125 -0.025 -0.975 0.575 0.875
0.000 0.000 0.200 -0.200 0.067 -0.333
0.000 0.000 0.000 1.000 -0.667 -1.167
0.000 0.000 0.000 0.000 0.167 -0.333
0.000 0.000 0.000 0.000 0.000 0.500

Instead, my code produces:
0.25 0 -0.1 -2.25 0.0333333 5.6
0 0.125 0 -1 -0.0833333 2
0 0 0.2 0 -0.0666667 -0.7
0 0 0 1 0 -2.5
0 0 0 0 0.166667 0
0 0 0 0 0 0.5


The diagonals work correctly, but the top is incorrect. Anyone familiar?

3
Contributors
3
Replies
8
Views
6 Years
Discussion Span
Last Post by Ivan_12
0

Looping up to j-2 sure looks dubious.

That solves it! Thank you so much!

0

Could you be more explicit about j-2? I can not associate it with a particular part of the code.

This question has already been answered. Start a new discussion instead.
Have something to contribute to this discussion? Please be thoughtful, detailed and courteous, and be sure to adhere to our posting rules.