Matrix inversion with MKL LAPACK

[swillmon]

Test driving my new Composer XE 2011 with a simple FORTRAN code that computes the inverse of a n x n matrix using the getrf and getri calls to LAPACK. Is there any rational reason why the results of the inverted matrix differ from what MATLAB is giving me?

[mecej4]

Are you sure that the original matrix being inverted is the same in your Matlab and Fortran runs? How big is the matrix? Is it possible to give an example where the inverse is noticeably different?

[swillmon]

Here are the arrays I've been working with / results:

Array
1.0658 -0.6313 0.0000 0.0000 0.0000 0.0000 0.0000
-0.6313 4.2632 -1.8940 0.0000 0.0000 0.0000 0.0000
0.0000 -1.8940 8.5265 -3.1567 0.0000 0.0000 0.0000
0.0000 0.0000 -3.1567 23.2603 -3.6892 0.0000 0.0000
0.0000 0.0000 0.0000 -3.6892 37.5680 -4.8192 0.0000
0.0000 0.0000 0.0000 0.0000 -4.8192 47.5467 -5.9492
0.0000 0.0000 0.0000 0.0000 0.0000 -5.9492 57.5255

GETRI Array
1.0401 0.1719 0.0402 0.0055 0.0006 0.0001 0.0000
0.1719 0.2902 0.0679 0.0094 0.0009 0.0001 0.0000
0.0402 0.0679 0.1395 0.0192 0.0019 0.0002 0.0000
0.0055 0.0094 0.0192 0.0463 0.0046 0.0005 0.0000
0.0006 0.0009 0.0019 0.0046 0.0274 0.0028 0.0003
0.0001 0.0001 0.0002 0.0005 0.0028 0.0216 0.0022
0.0000 0.0000 0.0000 0.0000 0.0003 0.0022 0.0176

MATLAB INV(ARRAY)
1.0042 0.1114 0.0182 0.0025 0.0003 0.0000 0.0000
0.1114 0.1881 0.0308 0.0042 0.0004 0.0000 0.0000
0.0182 0.0308 0.0915 0.0126 0.0013 0.0001 0.0000
0.0025 0.0042 0.0126 0.0454 0.0045 0.0005 0.0000
0.0003 0.0004 0.0013 0.0045 0.0274 0.0028 0.0003
0.0000 0.0000 0.0001 0.0005 0.0028 0.0216 0.0022
0.0000 0.0000 0.0000 0.0000 0.0003 0.0022 0.0176

Thanks for taking a look!

[styc]

[plain]$ matlab
Warning: No display specified.  You will not be able to display graphics on the screen.

                            < M A T L A B  >
                  Copyright 1984-2008 The MathWorks, Inc.
                         Version 7.6.0.324 (R2008a)
                             February 10, 2008


  To get started, type one of these: helpwin, helpdesk, or demo.
  For product information, visit www.mathworks.com.

>> A=[    1.0658   -0.6313    0.0000    0.0000    0.0000    0.0000    0.0000
   -0.6313    4.2632   -1.8940    0.0000    0.0000    0.0000    0.0000
    0.0000   -1.8940    8.5265   -3.1567    0.0000    0.0000    0.0000
    0.0000    0.0000   -3.1567   23.2603   -3.6892    0.0000    0.0000
    0.0000    0.0000    0.0000   -3.6892   37.5680   -4.8192    0.0000
    0.0000    0.0000    0.0000    0.0000   -4.8192   47.5467   -5.9492
    0.0000    0.0000    0.0000    0.0000    0.0000   -5.9492   57.5255
]
A =
    1.0658   -0.6313         0         0         0         0         0
   -0.6313    4.2632   -1.8940         0         0         0         0
         0   -1.8940    8.5265   -3.1567         0         0         0
         0         0   -3.1567   23.2603   -3.6892         0         0
         0         0         0   -3.6892   37.5680   -4.8192         0
         0         0         0         0   -4.8192   47.5467   -5.9492
         0         0         0         0         0   -5.9492   57.5255
>> inv(A)
ans =
    1.0401    0.1719    0.0402    0.0055    0.0006    0.0001    0.0000
    0.1719    0.2902    0.0679    0.0094    0.0009    0.0001    0.0000
    0.0402    0.0679    0.1395    0.0192    0.0019    0.0002    0.0000
    0.0055    0.0094    0.0192    0.0463    0.0046    0.0005    0.0000
    0.0006    0.0009    0.0019    0.0046    0.0274    0.0028    0.0003
    0.0001    0.0001    0.0002    0.0005    0.0028    0.0216    0.0022
    0.0000    0.0000    0.0000    0.0000    0.0003    0.0022    0.0176
>> quit
$
[/plain]

[mecej4 · Accepted solution]

With Matlab 7.5 (R2007b) I got results that agreed with your (and my) results from MKL Lapack95:

[bash]>> A=[   1.0658   -0.6313    0.0000    0.0000    0.0000    0.0000    0.0000
   -0.6313    4.2632   -1.8940    0.0000    0.0000    0.0000    0.0000
    0.0000   -1.8940    8.5265   -3.1567    0.0000    0.0000    0.0000
    0.0000    0.0000   -3.1567   23.2603   -3.6892    0.0000    0.0000
    0.0000    0.0000    0.0000   -3.6892   37.5680   -4.8192    0.0000
    0.0000    0.0000    0.0000    0.0000   -4.8192   47.5467   -5.9492
    0.0000    0.0000    0.0000    0.0000    0.0000   -5.9492   57.5255];
>> inv(A)

ans =

    1.0401    0.1719    0.0402    0.0055    0.0006    0.0001    0.0000
    0.1719    0.2902    0.0679    0.0094    0.0009    0.0001    0.0000
    0.0402    0.0679    0.1395    0.0192    0.0019    0.0002    0.0000
    0.0055    0.0094    0.0192    0.0463    0.0046    0.0005    0.0000
    0.0006    0.0009    0.0019    0.0046    0.0274    0.0028    0.0003
    0.0001    0.0001    0.0002    0.0005    0.0028    0.0216    0.0022
    0.0000    0.0000    0.0000    0.0000    0.0003    0.0022    0.0176

>> cond(A)

ans =

   65.3807

>> A-A'

ans =

     0     0     0     0     0     0     0
     0     0     0     0     0     0     0
     0     0     0     0     0     0     0
     0     0     0     0     0     0     0
     0     0     0     0     0     0     0
     0     0     0     0     0     0     0
     0     0     0     0     0     0     0[/bash]

[swillmon]

mecej4 - after your first reply, I went back and scrubbed my input data between what I was running w/ FORTRAN and what I put into MATLAB -- apples and oranges! Of course I get the same answer now that I'm using the same input data!

Thanks!