001:       SUBROUTINE ZTPTRI( 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:       COMPLEX*16         AP( * )
014: *     ..
015: *
016: *  Purpose
017: *  =======
018: *
019: *  ZTPTRI computes the inverse of a complex 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) COMPLEX*16 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:       COMPLEX*16         ONE, ZERO
072:       PARAMETER          ( ONE = ( 1.0D+0, 0.0D+0 ),
073:      $                   ZERO = ( 0.0D+0, 0.0D+0 ) )
074: *     ..
075: *     .. Local Scalars ..
076:       LOGICAL            NOUNIT, UPPER
077:       INTEGER            J, JC, JCLAST, JJ
078:       COMPLEX*16         AJJ
079: *     ..
080: *     .. External Functions ..
081:       LOGICAL            LSAME
082:       EXTERNAL           LSAME
083: *     ..
084: *     .. External Subroutines ..
085:       EXTERNAL           XERBLA, ZSCAL, ZTPMV
086: *     ..
087: *     .. Executable Statements ..
088: *
089: *     Test the input parameters.
090: *
091:       INFO = 0
092:       UPPER = LSAME( UPLO, 'U' )
093:       NOUNIT = LSAME( DIAG, 'N' )
094:       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
095:          INFO = -1
096:       ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN
097:          INFO = -2
098:       ELSE IF( N.LT.0 ) THEN
099:          INFO = -3
100:       END IF
101:       IF( INFO.NE.0 ) THEN
102:          CALL XERBLA( 'ZTPTRI', -INFO )
103:          RETURN
104:       END IF
105: *
106: *     Check for singularity if non-unit.
107: *
108:       IF( NOUNIT ) THEN
109:          IF( UPPER ) THEN
110:             JJ = 0
111:             DO 10 INFO = 1, N
112:                JJ = JJ + INFO
113:                IF( AP( JJ ).EQ.ZERO )
114:      $            RETURN
115:    10       CONTINUE
116:          ELSE
117:             JJ = 1
118:             DO 20 INFO = 1, N
119:                IF( AP( JJ ).EQ.ZERO )
120:      $            RETURN
121:                JJ = JJ + N - INFO + 1
122:    20       CONTINUE
123:          END IF
124:          INFO = 0
125:       END IF
126: *
127:       IF( UPPER ) THEN
128: *
129: *        Compute inverse of upper triangular matrix.
130: *
131:          JC = 1
132:          DO 30 J = 1, N
133:             IF( NOUNIT ) THEN
134:                AP( JC+J-1 ) = ONE / AP( JC+J-1 )
135:                AJJ = -AP( JC+J-1 )
136:             ELSE
137:                AJJ = -ONE
138:             END IF
139: *
140: *           Compute elements 1:j-1 of j-th column.
141: *
142:             CALL ZTPMV( 'Upper', 'No transpose', DIAG, J-1, AP,
143:      $                  AP( JC ), 1 )
144:             CALL ZSCAL( J-1, AJJ, AP( JC ), 1 )
145:             JC = JC + J
146:    30    CONTINUE
147: *
148:       ELSE
149: *
150: *        Compute inverse of lower triangular matrix.
151: *
152:          JC = N*( N+1 ) / 2
153:          DO 40 J = N, 1, -1
154:             IF( NOUNIT ) THEN
155:                AP( JC ) = ONE / AP( JC )
156:                AJJ = -AP( JC )
157:             ELSE
158:                AJJ = -ONE
159:             END IF
160:             IF( J.LT.N ) THEN
161: *
162: *              Compute elements j+1:n of j-th column.
163: *
164:                CALL ZTPMV( 'Lower', 'No transpose', DIAG, N-J,
165:      $                     AP( JCLAST ), AP( JC+1 ), 1 )
166:                CALL ZSCAL( N-J, AJJ, AP( JC+1 ), 1 )
167:             END IF
168:             JCLAST = JC
169:             JC = JC - N + J - 2
170:    40    CONTINUE
171:       END IF
172: *
173:       RETURN
174: *
175: *     End of ZTPTRI
176: *
177:       END
178: