src/LAPACK/dtrtri.f

   1       SUBROUTINE DTRTRI( UPLO, DIAG, N, A, LDA, INFO )
   2 *
   3 *  -- LAPACK routine (version 3.1) --
   4 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
   5 *     November 2006
   6 *
   7 *     .. Scalar Arguments ..
   8       CHARACTER          DIAG, UPLO
   9       INTEGER            INFO, LDA, N
  10 *     ..
  11 *     .. Array Arguments ..
  12       DOUBLE PRECISION   A( LDA, * )
  13 *     ..
  14 *
  15 *  Purpose
  16 *  =======
  17 *
  18 *  DTRTRI computes the inverse of a real upper or lower triangular
  19 *  matrix A.
  20 *
  21 *  This is the Level 3 BLAS version of the algorithm.
  22 *
  23 *  Arguments
  24 *  =========
  25 *
  26 *  UPLO    (input) CHARACTER*1
  27 *          = 'U':  A is upper triangular;
  28 *          = 'L':  A is lower triangular.
  29 *
  30 *  DIAG    (input) CHARACTER*1
  31 *          = 'N':  A is non-unit triangular;
  32 *          = 'U':  A is unit triangular.
  33 *
  34 *  N       (input) INTEGER
  35 *          The order of the matrix A.  N >= 0.
  36 *
  37 *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
  38 *          On entry, the triangular matrix A.  If UPLO = 'U', the
  39 *          leading N-by-N upper triangular part of the array A contains
  40 *          the upper triangular matrix, and the strictly lower
  41 *          triangular part of A is not referenced.  If UPLO = 'L', the
  42 *          leading N-by-N lower triangular part of the array A contains
  43 *          the lower triangular matrix, and the strictly upper
  44 *          triangular part of A is not referenced.  If DIAG = 'U', the
  45 *          diagonal elements of A are also not referenced and are
  46 *          assumed to be 1.
  47 *          On exit, the (triangular) inverse of the original matrix, in
  48 *          the same storage format.
  49 *
  50 *  LDA     (input) INTEGER
  51 *          The leading dimension of the array A.  LDA >= max(1,N).
  52 *
  53 *  INFO    (output) INTEGER
  54 *          = 0: successful exit
  55 *          < 0: if INFO = -i, the i-th argument had an illegal value
  56 *          > 0: if INFO = i, A(i,i) is exactly zero.  The triangular
  57 *               matrix is singular and its inverse can not be computed.
  58 *
  59 *  =====================================================================
  60 *
  61 *     .. Parameters ..
  62       DOUBLE PRECISION   ONE, ZERO
  63       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
  64 *     ..
  65 *     .. Local Scalars ..
  66       LOGICAL            NOUNIT, UPPER
  67       INTEGER            J, JB, NB, NN
  68 *     ..
  69 *     .. External Functions ..
  70       LOGICAL            LSAME
  71       INTEGER            ILAENV
  72       EXTERNAL           LSAME, ILAENV
  73 *     ..
  74 *     .. External Subroutines ..
  75       EXTERNAL           DTRMM, DTRSM, DTRTI2, XERBLA
  76 *     ..
  77 *     .. Intrinsic Functions ..
  78       INTRINSIC          MAX, MIN
  79 *     ..
  80 *     .. Executable Statements ..
  81 *
  82 *     Test the input parameters.
  83 *
  84       INFO = 0
  85       UPPER = LSAME( UPLO, 'U' )
  86       NOUNIT = LSAME( DIAG, 'N' )
  87       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
  88          INFO = -1
  89       ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN
  90          INFO = -2
  91       ELSE IF( N.LT.0 ) THEN
  92          INFO = -3
  93       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
  94          INFO = -5
  95       END IF
  96       IF( INFO.NE.0 ) THEN
  97          CALL XERBLA( 'DTRTRI', -INFO )
  98          RETURN
  99       END IF
 100 *
 101 *     Quick return if possible
 102 *
 103       IF( N.EQ.0 )
 104      $   RETURN
 105 *
 106 *     Check for singularity if non-unit.
 107 *
 108       IF( NOUNIT ) THEN
 109          DO 10 INFO = 1, N
 110             IF( A( INFO, INFO ).EQ.ZERO )
 111      $         RETURN
 112    10    CONTINUE
 113          INFO = 0
 114       END IF
 115 *
 116 *     Determine the block size for this environment.
 117 *
 118       NB = ILAENV( 1, 'DTRTRI', UPLO // DIAG, N, -1, -1, -1 )
 119       IF( NB.LE.1 .OR. NB.GE.N ) THEN
 120 *
 121 *        Use unblocked code
 122 *
 123          CALL DTRTI2( UPLO, DIAG, N, A, LDA, INFO )
 124       ELSE
 125 *
 126 *        Use blocked code
 127 *
 128          IF( UPPER ) THEN
 129 *
 130 *           Compute inverse of upper triangular matrix
 131 *
 132             DO 20 J = 1, N, NB
 133                JB = MIN( NB, N-J+1 )
 134 *
 135 *              Compute rows 1:j-1 of current block column
 136 *
 137                CALL DTRMM( 'Left', 'Upper', 'No transpose', DIAG, J-1,
 138      $                     JB, ONE, A, LDA, A( 1, J ), LDA )
 139                CALL DTRSM( 'Right', 'Upper', 'No transpose', DIAG, J-1,
 140      $                     JB, -ONE, A( J, J ), LDA, A( 1, J ), LDA )
 141 *
 142 *              Compute inverse of current diagonal block
 143 *
 144                CALL DTRTI2( 'Upper', DIAG, JB, A( J, J ), LDA, INFO )
 145    20       CONTINUE
 146          ELSE
 147 *
 148 *           Compute inverse of lower triangular matrix
 149 *
 150             NN = ( ( N-1 ) / NB )*NB + 1
 151             DO 30 J = NN, 1, -NB
 152                JB = MIN( NB, N-J+1 )
 153                IF( J+JB.LE.N ) THEN
 154 *
 155 *                 Compute rows j+jb:n of current block column
 156 *
 157                   CALL DTRMM( 'Left', 'Lower', 'No transpose', DIAG,
 158      $                        N-J-JB+1, JB, ONE, A( J+JB, J+JB ), LDA,
 159      $                        A( J+JB, J ), LDA )
 160                   CALL DTRSM( 'Right', 'Lower', 'No transpose', DIAG,
 161      $                        N-J-JB+1, JB, -ONE, A( J, J ), LDA,
 162      $                        A( J+JB, J ), LDA )
 163                END IF
 164 *
 165 *              Compute inverse of current diagonal block
 166 *
 167                CALL DTRTI2( 'Lower', DIAG, JB, A( J, J ), LDA, INFO )
 168    30       CONTINUE
 169          END IF
 170       END IF
 171 *
 172       RETURN
 173 *
 174 *     End of DTRTRI
 175 *
 176       END