001:       SUBROUTINE DTPTRI( UPLO, DIAG, N, AP, INFO )
002: *
003: *  -- LAPACK routine (version 3.2) --
004: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
005: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
006: *     November 2006
007: *
008: *     .. Scalar Arguments ..
009:       CHARACTER          DIAG, UPLO
010:       INTEGER            INFO, N
011: *     ..
012: *     .. Array Arguments ..
013:       DOUBLE PRECISION   AP( * )
014: *     ..
015: *
016: *  Purpose
017: *  =======
018: *
019: *  DTPTRI computes the inverse of a real upper or lower triangular
020: *  matrix A stored in packed format.
021: *
022: *  Arguments
023: *  =========
024: *
025: *  UPLO    (input) CHARACTER*1
026: *          = 'U':  A is upper triangular;
027: *          = 'L':  A is lower triangular.
028: *
029: *  DIAG    (input) CHARACTER*1
030: *          = 'N':  A is non-unit triangular;
031: *          = 'U':  A is unit triangular.
032: *
033: *  N       (input) INTEGER
034: *          The order of the matrix A.  N >= 0.
035: *
036: *  AP      (input/output) DOUBLE PRECISION array, dimension (N*(N+1)/2)
037: *          On entry, the upper or lower triangular matrix A, stored
038: *          columnwise in a linear array.  The j-th column of A is stored
039: *          in the array AP as follows:
040: *          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
041: *          if UPLO = 'L', AP(i + (j-1)*((2*n-j)/2) = A(i,j) for j<=i<=n.
042: *          See below for further details.
043: *          On exit, the (triangular) inverse of the original matrix, in
044: *          the same packed storage format.
045: *
046: *  INFO    (output) INTEGER
047: *          = 0:  successful exit
048: *          < 0:  if INFO = -i, the i-th argument had an illegal value
049: *          > 0:  if INFO = i, A(i,i) is exactly zero.  The triangular
050: *                matrix is singular and its inverse can not be computed.
051: *
052: *  Further Details
053: *  ===============
054: *
055: *  A triangular matrix A can be transferred to packed storage using one
056: *  of the following program segments:
057: *
058: *  UPLO = 'U':                      UPLO = 'L':
059: *
060: *        JC = 1                           JC = 1
061: *        DO 2 J = 1, N                    DO 2 J = 1, N
062: *           DO 1 I = 1, J                    DO 1 I = J, N
063: *              AP(JC+I-1) = A(I,J)              AP(JC+I-J) = A(I,J)
064: *      1    CONTINUE                    1    CONTINUE
065: *           JC = JC + J                      JC = JC + N - J + 1
066: *      2 CONTINUE                       2 CONTINUE
067: *
068: *  =====================================================================
069: *
070: *     .. Parameters ..
071:       DOUBLE PRECISION   ONE, ZERO
072:       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
073: *     ..
074: *     .. Local Scalars ..
075:       LOGICAL            NOUNIT, UPPER
076:       INTEGER            J, JC, JCLAST, JJ
077:       DOUBLE PRECISION   AJJ
078: *     ..
079: *     .. External Functions ..
080:       LOGICAL            LSAME
081:       EXTERNAL           LSAME
082: *     ..
083: *     .. External Subroutines ..
084:       EXTERNAL           DSCAL, DTPMV, XERBLA
085: *     ..
086: *     .. Executable Statements ..
087: *
088: *     Test the input parameters.
089: *
090:       INFO = 0
091:       UPPER = LSAME( UPLO, 'U' )
092:       NOUNIT = LSAME( DIAG, 'N' )
093:       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
094:          INFO = -1
095:       ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN
096:          INFO = -2
097:       ELSE IF( N.LT.0 ) THEN
098:          INFO = -3
099:       END IF
100:       IF( INFO.NE.0 ) THEN
101:          CALL XERBLA( 'DTPTRI', -INFO )
102:          RETURN
103:       END IF
104: *
105: *     Check for singularity if non-unit.
106: *
107:       IF( NOUNIT ) THEN
108:          IF( UPPER ) THEN
109:             JJ = 0
110:             DO 10 INFO = 1, N
111:                JJ = JJ + INFO
112:                IF( AP( JJ ).EQ.ZERO )
113:      $            RETURN
114:    10       CONTINUE
115:          ELSE
116:             JJ = 1
117:             DO 20 INFO = 1, N
118:                IF( AP( JJ ).EQ.ZERO )
119:      $            RETURN
120:                JJ = JJ + N - INFO + 1
121:    20       CONTINUE
122:          END IF
123:          INFO = 0
124:       END IF
125: *
126:       IF( UPPER ) THEN
127: *
128: *        Compute inverse of upper triangular matrix.
129: *
130:          JC = 1
131:          DO 30 J = 1, N
132:             IF( NOUNIT ) THEN
133:                AP( JC+J-1 ) = ONE / AP( JC+J-1 )
134:                AJJ = -AP( JC+J-1 )
135:             ELSE
136:                AJJ = -ONE
137:             END IF
138: *
139: *           Compute elements 1:j-1 of j-th column.
140: *
141:             CALL DTPMV( 'Upper', 'No transpose', DIAG, J-1, AP,
142:      $                  AP( JC ), 1 )
143:             CALL DSCAL( J-1, AJJ, AP( JC ), 1 )
144:             JC = JC + J
145:    30    CONTINUE
146: *
147:       ELSE
148: *
149: *        Compute inverse of lower triangular matrix.
150: *
151:          JC = N*( N+1 ) / 2
152:          DO 40 J = N, 1, -1
153:             IF( NOUNIT ) THEN
154:                AP( JC ) = ONE / AP( JC )
155:                AJJ = -AP( JC )
156:             ELSE
157:                AJJ = -ONE
158:             END IF
159:             IF( J.LT.N ) THEN
160: *
161: *              Compute elements j+1:n of j-th column.
162: *
163:                CALL DTPMV( 'Lower', 'No transpose', DIAG, N-J,
164:      $                     AP( JCLAST ), AP( JC+1 ), 1 )
165:                CALL DSCAL( N-J, AJJ, AP( JC+1 ), 1 )
166:             END IF
167:             JCLAST = JC
168:             JC = JC - N + J - 2
169:    40    CONTINUE
170:       END IF
171: *
172:       RETURN
173: *
174: *     End of DTPTRI
175: *
176:       END
177: