001:       SUBROUTINE DTRTRI( UPLO, DIAG, N, A, LDA, INFO )
002: *
003: *  -- LAPACK routine (version 3.2) --
004: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
005: *     November 2006
006: *
007: *     .. Scalar Arguments ..
008:       CHARACTER          DIAG, UPLO
009:       INTEGER            INFO, LDA, N
010: *     ..
011: *     .. Array Arguments ..
012:       DOUBLE PRECISION   A( LDA, * )
013: *     ..
014: *
015: *  Purpose
016: *  =======
017: *
018: *  DTRTRI computes the inverse of a real upper or lower triangular
019: *  matrix A.
020: *
021: *  This is the Level 3 BLAS version of the algorithm.
022: *
023: *  Arguments
024: *  =========
025: *
026: *  UPLO    (input) CHARACTER*1
027: *          = 'U':  A is upper triangular;
028: *          = 'L':  A is lower triangular.
029: *
030: *  DIAG    (input) CHARACTER*1
031: *          = 'N':  A is non-unit triangular;
032: *          = 'U':  A is unit triangular.
033: *
034: *  N       (input) INTEGER
035: *          The order of the matrix A.  N >= 0.
036: *
037: *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
038: *          On entry, the triangular matrix A.  If UPLO = 'U', the
039: *          leading N-by-N upper triangular part of the array A contains
040: *          the upper triangular matrix, and the strictly lower
041: *          triangular part of A is not referenced.  If UPLO = 'L', the
042: *          leading N-by-N lower triangular part of the array A contains
043: *          the lower triangular matrix, and the strictly upper
044: *          triangular part of A is not referenced.  If DIAG = 'U', the
045: *          diagonal elements of A are also not referenced and are
046: *          assumed to be 1.
047: *          On exit, the (triangular) inverse of the original matrix, in
048: *          the same storage format.
049: *
050: *  LDA     (input) INTEGER
051: *          The leading dimension of the array A.  LDA >= max(1,N).
052: *
053: *  INFO    (output) INTEGER
054: *          = 0: successful exit
055: *          < 0: if INFO = -i, the i-th argument had an illegal value
056: *          > 0: if INFO = i, A(i,i) is exactly zero.  The triangular
057: *               matrix is singular and its inverse can not be computed.
058: *
059: *  =====================================================================
060: *
061: *     .. Parameters ..
062:       DOUBLE PRECISION   ONE, ZERO
063:       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
064: *     ..
065: *     .. Local Scalars ..
066:       LOGICAL            NOUNIT, UPPER
067:       INTEGER            J, JB, NB, NN
068: *     ..
069: *     .. External Functions ..
070:       LOGICAL            LSAME
071:       INTEGER            ILAENV
072:       EXTERNAL           LSAME, ILAENV
073: *     ..
074: *     .. External Subroutines ..
075:       EXTERNAL           DTRMM, DTRSM, DTRTI2, XERBLA
076: *     ..
077: *     .. Intrinsic Functions ..
078:       INTRINSIC          MAX, MIN
079: *     ..
080: *     .. Executable Statements ..
081: *
082: *     Test the input parameters.
083: *
084:       INFO = 0
085:       UPPER = LSAME( UPLO, 'U' )
086:       NOUNIT = LSAME( DIAG, 'N' )
087:       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
088:          INFO = -1
089:       ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN
090:          INFO = -2
091:       ELSE IF( N.LT.0 ) THEN
092:          INFO = -3
093:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
094:          INFO = -5
095:       END IF
096:       IF( INFO.NE.0 ) THEN
097:          CALL XERBLA( 'DTRTRI', -INFO )
098:          RETURN
099:       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
177: