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
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
(DLANSY) before DGETRF
* 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.
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
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
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