Hi,

I'm using gaussj to obtain the inverse of a matrix a(1:3,1:3).
According to the comments in the routine,
1) a(1:n,1:n) is an input matrix stored in an array of physical dimensions np by np.
2) b(1:n,1:m) is an input matrix containing the m right-hand side vectors, stored in an array of physical dimensions np by mp.
3) On output, a(1:n,1:n) is replaced by its matrix inverse, and b(1:n,1:m) is replaced by the corresponding set of solution vectors.

However, I find that on output, a(1:3,1:3) does not store its inverse matrix, even though b(1:n,1:m) contains the correct solution! Here's an example:

initial a(1:3,1:3)
2.000000 2.000000 0.0000000E+00
0.0000000E+00 4.000000 0.0000000E+00
0.0000000E+00 0.0000000E+00 6.000000

initial b(1:3,1:3) (identity matrix)
1.000000 0.0000000E+00 0.0000000E+00
0.0000000E+00 1.000000 0.0000000E+00
0.0000000E+00 0.0000000E+00 1.000000


a(1:3,1:3) on output (should be inverse matrix but it is not)
0.0000000E+00 0.5000000 -0.2500000
0.0000000E+00 0.0000000E+00 0.2500000
1.000000 0.0000000E+00 0.0000000E+00

b(1:3,1:3) on output (in this case, it is the inverse matrix because b(1:3,1:3)=identity matrix initially)
0.5000000 -0.2500000 0.0000000E+00
0.0000000E+00 0.2500000 0.0000000E+00
0.0000000E+00 0.0000000E+00 0.1666667

Can anyone explain what I am doing wrong?

Thanks,

Ana

Hi Ana,

Sorry, I am not able to reproduce your results. The original matrix is replaced by its inverse when I use your test matrix.

Saul Teukolsky

Thanks for the quick reply and for having taken the time to test gaussj. After your message I had a look (for the 100th time...) at my subroutine and I found the error. I had accidentally deleted one line...

Thank you for helping the newbie-with-dumb-question. Before your reply I was wondering whether I had misunderstood the comments that came with the subroutine (I'm not a native English speaker), because I had already looked for the error without success. Your reply gave me the certainty I needed to keep looking.

Ana Vila Verde