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 (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. 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 blas-3.1.1-1.fc7.i386.rpm 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 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.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 ? `````` `````` ```
 Current Thread [Lapack] Lapack's DGETRF+DGETRI for matrix inversion fails quite badly, Alex A. Schmidt [Lapack] Lapack's DGETRF+DGETRI for matrix inversion fails quite badly, Julien Langou [Lapack] Lapack's DGETRF+DGETRI for matrix inversion fails quite badly, Julien Langou <=

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