Welcome in the world of finite precision arithmetic. Operations in a
computer are made with roundoff 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
(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 < 1e14) for example.
Best wishes,
Julien
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 box
running Fedora 7. Lapack and Blas libraries are from rpm packages
blas3.1.11.fc7.i386.rpm
and lapack3.1.11.fc7.i386.rpm.
DGETRF+DGETRI seem to work just fine but since my (fortran77) 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
calling DGETRF+DGETRI
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.1933E01 5.5129E02 1.3698E+00? 1.1371E+00
?2.4020E+00 1.8960E+00? 2.4020E+00? 1.5350E+00??? 5.9788E02 3.6745E02
4.5782E02? 2.7993E01??? 4.1761E01 5.2895E01? 1.7921E02? 2.9202E01
3.0010E+00 4.2300E+00 3.0010E+00 7.7700E01??? 1.7889E+13 6.9052E+15
6.3183E+15 1.3882E+16?? 5.2300E01? 4.3789E01 3.2938E01? 1.9254E01
?1.6900E01? 2.8670E+00? 1.6900E01 4.0700E01??? 3.4841E01 1.1111E01
2.9752E01 3.7852E01??? 2.9773E02? 1.2005E01 5.9453E02? 1.0029E+00
?3.3920E+00 7.2800E01 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.7344E02
?1.4470E+00 3.6190E+00 3.7510E+00? 1.4470E+00?? 5.9252E01? 1.5156E01?
1.6178E01? 9.3682E01??? 2.9590E+00 2.1631E+00 5.6195E01 5.6016E+00
4.8280E+00? 1.1900E+00 2.8370E+00 4.8280E+00??? 2.7498E01 2.5928E01
1.8376E01 5.2632E01?? 3.6172E+00? 4.4102E+00? 3.2223E+00 3.0625E+00
?2.7450E+00? 9.8700E01? 1.8600E01? 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.8510E01 9.7194E01
1.1424E+00 9.5392E01??? 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.7891E01
?3.7260E+00 1.7950E+00 2.6670E+00 1.7950E+00?? 1.2906E01? 1.6487E01?
7.2103E02? 3.1464E03??? 4.9102E+00 3.3594E01 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
?
