1        2        3        4        5        6        7        8        9       10       11       12       13       14       15       16       17       18       19       20       21       22       23       24       25       26       27       28       29       30       31       32       33       34       35       36       37       38       39       40       41       42       43       44       45       46       47       48       49       50       51       52       53       54       55       56       57       58       59       60       61       62       63       64       65       66       67       68       69       70       71       72       73       74       75       76       77       78       79       80       81       82       83       84       85       86       87       88       89       90       91       92       93       94       95       96       97       98       99      100      101      102      103      104      105      106      107      108      109      110      111      112      113      114      115      116      117      118      119      120      121      122      123      124      125      126      127      128      129      130      131      132      133      134      135      136      137      138      139      140      141      142      143      144      145      146      147      148      149      150      151      152      153      154      155      156      157      158      159      160      161      162      163      164      165      166      167      168      169      170      171      172      173      174      175      176      177      178      179      180      181      182      183      184      185      186      187      188      189      190      191      192      193 SUBROUTINE DGETRI( N, A, LDA, IPIV, WORK, LWORK, INFO ) * *  -- LAPACK routine (version 3.2) -- *  -- LAPACK is a software package provided by Univ. of Tennessee,    -- *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- *     November 2006 * *     .. Scalar Arguments ..       INTEGER            INFO, LDA, LWORK, N *     .. *     .. Array Arguments ..       INTEGER            IPIV( * )       DOUBLE PRECISION   A( LDA, * ), WORK( * ) *     .. * *  Purpose *  ======= * *  DGETRI computes the inverse of a matrix using the LU factorization *  computed by DGETRF. * *  This method inverts U and then computes inv(A) by solving the system *  inv(A)*L = inv(U) for inv(A). * *  Arguments *  ========= * *  N       (input) INTEGER *          The order of the matrix A.  N >= 0. * *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N) *          On entry, the factors L and U from the factorization *          A = P*L*U as computed by DGETRF. *          On exit, if INFO = 0, the inverse of the original matrix A. * *  LDA     (input) INTEGER *          The leading dimension of the array A.  LDA >= max(1,N). * *  IPIV    (input) INTEGER array, dimension (N) *          The pivot indices from DGETRF; for 1<=i<=N, row i of the *          matrix was interchanged with row IPIV(i). * *  WORK    (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) *          On exit, if INFO=0, then WORK(1) returns the optimal LWORK. * *  LWORK   (input) INTEGER *          The dimension of the array WORK.  LWORK >= max(1,N). *          For optimal performance LWORK >= N*NB, where NB is *          the optimal blocksize returned by ILAENV. * *          If LWORK = -1, then a workspace query is assumed; the routine *          only calculates the optimal size of the WORK array, returns *          this value as the first entry of the WORK array, and no error *          message related to LWORK is issued by XERBLA. * *  INFO    (output) INTEGER *          = 0:  successful exit *          < 0:  if INFO = -i, the i-th argument had an illegal value *          > 0:  if INFO = i, U(i,i) is exactly zero; the matrix is *                singular and its inverse could not be computed. * *  ===================================================================== * *     .. Parameters ..       DOUBLE PRECISION   ZERO, ONE       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 ) *     .. *     .. Local Scalars ..       LOGICAL            LQUERY       INTEGER            I, IWS, J, JB, JJ, JP, LDWORK, LWKOPT, NB,      $NBMIN, NN * .. * .. External Functions .. INTEGER ILAENV EXTERNAL ILAENV * .. * .. External Subroutines .. EXTERNAL DGEMM, DGEMV, DSWAP, DTRSM, DTRTRI, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX, MIN * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 NB = ILAENV( 1, 'DGETRI', ' ', N, -1, -1, -1 ) LWKOPT = N*NB WORK( 1 ) = LWKOPT LQUERY = ( LWORK.EQ.-1 ) IF( N.LT.0 ) THEN INFO = -1 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN INFO = -3 ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN INFO = -6 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'DGETRI', -INFO ) RETURN ELSE IF( LQUERY ) THEN RETURN END IF * * Quick return if possible * IF( N.EQ.0 )$   RETURN * *     Form inv(U).  If INFO > 0 from DTRTRI, then U is singular, *     and the inverse is not computed. *       CALL DTRTRI( 'Upper', 'Non-unit', N, A, LDA, INFO )       IF( INFO.GT.0 )      $RETURN * NBMIN = 2 LDWORK = N IF( NB.GT.1 .AND. NB.LT.N ) THEN IWS = MAX( LDWORK*NB, 1 ) IF( LWORK.LT.IWS ) THEN NB = LWORK / LDWORK NBMIN = MAX( 2, ILAENV( 2, 'DGETRI', ' ', N, -1, -1, -1 ) ) END IF ELSE IWS = N END IF * * Solve the equation inv(A)*L = inv(U) for inv(A). * IF( NB.LT.NBMIN .OR. NB.GE.N ) THEN * * Use unblocked code. * DO 20 J = N, 1, -1 * * Copy current column of L to WORK and replace with zeros. * DO 10 I = J + 1, N WORK( I ) = A( I, J ) A( I, J ) = ZERO 10 CONTINUE * * Compute current column of inv(A). * IF( J.LT.N )$         CALL DGEMV( 'No transpose', N, N-J, -ONE, A( 1, J+1 ),      $LDA, WORK( J+1 ), 1, ONE, A( 1, J ), 1 ) 20 CONTINUE ELSE * * Use blocked code. * NN = ( ( N-1 ) / NB )*NB + 1 DO 50 J = NN, 1, -NB JB = MIN( NB, N-J+1 ) * * Copy current block column of L to WORK and replace with * zeros. * DO 40 JJ = J, J + JB - 1 DO 30 I = JJ + 1, N WORK( I+( JJ-J )*LDWORK ) = A( I, JJ ) A( I, JJ ) = ZERO 30 CONTINUE 40 CONTINUE * * Compute current block column of inv(A). * IF( J+JB.LE.N )$         CALL DGEMM( 'No transpose', 'No transpose', N, JB,      $N-J-JB+1, -ONE, A( 1, J+JB ), LDA,$                     WORK( J+JB ), LDWORK, ONE, A( 1, J ), LDA )             CALL DTRSM( 'Right', 'Lower', 'No transpose', 'Unit', N, JB,      $ONE, WORK( J ), LDWORK, A( 1, J ), LDA ) 50 CONTINUE END IF * * Apply column interchanges. * DO 60 J = N - 1, 1, -1 JP = IPIV( J ) IF( JP.NE.J )$      CALL DSWAP( N, A( 1, J ), 1, A( 1, JP ), 1 )    60 CONTINUE *       WORK( 1 ) = IWS       RETURN * *     End of DGETRI *       END