LAPACK Archives

[Lapack] Lapack's DGETRF+DGETRI for matrix inversion fails quite badly

Oupseeeee ...
Merci Antoine (Petitet) for catching my mistake.

In my previous email, please replace all occurences of DLANSY with DLANGE!

On Sat, 7 Nov 2009, Julien Langou wrote:

Welcome in the world of finite precision arithmetic. Operations in a computer 
are made with round-off errors and results are sometime quite not the same as 
in exact arithmetic.

So yes, if a matrix has two identical columns or two identical rows, it is 
singular so not invertible. From the LU factorization point of view this 
means that we should get an exact zero pivot during the factorization. And so 
the LU factorization algorithm breaks.

DGETRF does check for an exact zero as pivot. If the pivot is zero then you 
have an error message. However it is unlikely to obtain this zero pivot in 
finite precision computation. You'll get a very small number as a pivot but 
not zero. We are not bold enough to divide by zero, but we are bold enough to 
divide by a very small number. So the computation in general proceeds and, to 
your dismay, DGETRF+DGETRI completes succesfully.

During the LU factorization algorithm, the operations performed on the rows 
of A are often exactly the same. Therefore it is likely that if you start 
with two same rows, perform the same operations on their elements, you will 
obtain a zero pivot. Regarding columns, the operations during the LU 
factorization are completely different. Therefore, in this case, it is 
unlikely to obtain a zero pivot. This explain your observation on the INFO 
error messages.

If you really want to test for numerical singularity (non invertibility) of 
your matrix, you will need to use DGECON. Its use is "similar" to DGETRI in 
the sense that the input of DGECON are the L and U factors of the inital 
matrix A. So ou first need to run DGETRF to use it. A pseudo code would be:

NRMA = DLANSY ( )  % compute the norm of A (1, infty)
DGETRF             % compute the LU factors of A
DGECON             % compute the (one over the) condition number of A
DGETRI             % compute the inverse of A from the LU factors

* Note 1: DGECON needs NRMA and the LU factors of A, the LU factors of A
 are input only.

* Note 2: Since DGETRF overwrites A, you want to compute the norm of A

* Note 3: The cost of DLANSY and DGECON is O(n^2) so is neglictible over
 DGETRF and DGETRI O(n^3). Well, if n=4 ... not that neglictible.

* Note 4: So you run DGECON and if you have a large condition number, you
 quit and do not call DGETRI. You define large. Say 1e14. So
 (RCOND < 1e-14) for example.

Best wishes,

On Fri, 6 Nov 2009, Alex A. Schmidt wrote:

Dear Lapack developers team,

I have been using lapack's DGETRF+DGETRI to do matrix inversion in a linux 
running Fedora 7. Lapack and Blas libraries are from rpm packages 
and lapack-3.1.1-1.fc7.i386.rpm.

DGETRF+DGETRI seem to work just fine but since my (fortran-77) program use 
them for a
lot of matrixes I made a small program to test these two subroutines, by 
generating random
4x4 matrixes to calculate their inverses. Everything went fine until I 
decided to test the reaction
of? DGETRF+DGETRI when it finds two equal rows or two equal columns in the 
matrix to be
inverted. To my surprise the routines sometimes still try to do the 
inversion process. Eventually
I get a error flag from DGETRF+DGETRI when the matrix have two equal lines 
but not always and
I never get an error message when the matrix has two equal columns. I 
guess, this has something
to do about the LU decomposition behind the strategy of DGETRF+DGETRI to do 
the matrix inversion.

Some examples of the matrixes where this has happened are shown below.

Of course, one should never call DGETRF+DGETRI to do matrix inversion when 
the matrix has two
or more equal lines (or columns). Should I test the matrix myself before 
or should these routines find out this kind of situation by themselves?

Thanks in advance for you attention,

Alex A. Schmidt
Dep. of Mathematics/UFSM, Brazil

???????????????????????? 4x4 matrix ? 
??????????????????????????????????????????????????? 4x4 inverse ? ? ? ? ? ? 
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?
?? 4x4 product ?

-1.6130E+00 -1.1020E+00 -1.6130E+00? 2.2540E+00?? -1.7889E+13? 6.9052E+15? 
6.3183E+15? 1.3882E+16??? 7.1933E-01 -5.5129E-02 -1.3698E+00? 1.1371E+00
?2.4020E+00 -1.8960E+00? 2.4020E+00? 1.5350E+00??? 5.9788E-02 -3.6745E-02 
-4.5782E-02? 2.7993E-01??? 4.1761E-01 -5.2895E-01? 1.7921E-02? 2.9202E-01
-3.0010E+00 -4.2300E+00 -3.0010E+00 -7.7700E-01??? 1.7889E+13 -6.9052E+15 
-6.3183E+15 -1.3882E+16?? -5.2300E-01? 4.3789E-01 -3.2938E-01? 1.9254E-01
?1.6900E-01? 2.8670E+00? 1.6900E-01 -4.0700E-01??? 3.4841E-01 -1.1111E-01 
-2.9752E-01 -3.7852E-01??? 2.9773E-02? 1.2005E-01 -5.9453E-02? 1.0029E+00

?3.3920E+00 -7.2800E-01 -1.4400E+00? 3.3920E+00??? 1.5132E+16 -6.8311E+15? 
4.0857E+14 -1.4379E+16?? -1.8125E+00 -2.2559E+00 -1.4188E+00? 2.7344E-02
?1.4470E+00 -3.6190E+00 -3.7510E+00? 1.4470E+00?? -5.9252E-01? 1.5156E-01? 
1.6178E-01? 9.3682E-01??? 2.9590E+00 -2.1631E+00 -5.6195E-01 -5.6016E+00
-4.8280E+00? 1.1900E+00 -2.8370E+00 -4.8280E+00??? 2.7498E-01 -2.5928E-01 
-1.8376E-01 -5.2632E-01?? -3.6172E+00? 4.4102E+00? 3.2223E+00 -3.0625E+00
?2.7450E+00? 9.8700E-01? 1.8600E-01? 2.7450E+00?? -1.5132E+16? 6.8311E+15 
-4.0857E+14? 1.4379E+16?? -2.4023E+00 -4.2715E+00 -1.2907E+00 -3.3359E+00

?4.2280E+00? 3.4440E+00 -3.9640E+00? 3.4440E+00??? 7.8510E-01 -9.7194E-01 
-1.1424E+00 -9.5392E-01??? 6.1641E+00 -2.2148E+00 -5.5859E+00 -5.7617E+00
?1.7530E+00? 3.4530E+00? 4.2120E+00? 3.4530E+00?? -1.5059E+16? 2.0766E+16? 
2.5755E+16? 1.7998E+16??? 4.6445E+00 -3.3047E+00 -6.9297E+00? 8.7891E-01
?3.7260E+00 -1.7950E+00 -2.6670E+00 -1.7950E+00?? -1.2906E-01? 1.6487E-01? 
7.2103E-02? 3.1464E-03??? 4.9102E+00 -3.3594E-01 -6.8086E+00 -5.1992E+00
-3.8170E+00? 1.4660E+00 -4.3600E+00? 1.4660E+00??? 1.5059E+16 -2.0766E+16 
-2.5755E+16 -1.7998E+16?? -3.6875E+00? 4.1543E+00? 3.6836E+00? 3.3594E+00

<Prev in Thread] Current Thread [Next in Thread>

For additional information you may use the LAPACK/ScaLAPACK Forum.
Or one of the mailing lists, or