topic Matrix inversion with MKL LAPACK in Intel® oneAPI Math Kernel Library

Matrix inversion with MKL LAPACK

swillmon — Fri, 02 Sep 2011 21:43:45 GMT

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?

Matrix inversion with MKL LAPACK

mecej4 — Sat, 03 Sep 2011 02:37:02 GMT

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?

Matrix inversion with MKL LAPACK

swillmon — Sat, 03 Sep 2011 12:56:48 GMT

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!

Matrix inversion with MKL LAPACK

styc — Sat, 03 Sep 2011 16:59:08 GMT

[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]

Matrix inversion with MKL LAPACK

mecej4 — Sat, 03 Sep 2011 17:13:54 GMT

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]

Matrix inversion with MKL LAPACK

swillmon — Sat, 03 Sep 2011 17:52:02 GMT

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!