001:       SUBROUTINE ZSYTRI( UPLO, N, A, LDA, IPIV, WORK, 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          UPLO
010:       INTEGER            INFO, LDA, N
011: *     ..
012: *     .. Array Arguments ..
013:       INTEGER            IPIV( * )
014:       COMPLEX*16         A( LDA, * ), WORK( * )
015: *     ..
016: *
017: *  Purpose
018: *  =======
019: *
020: *  ZSYTRI computes the inverse of a complex symmetric indefinite matrix
021: *  A using the factorization A = U*D*U**T or A = L*D*L**T computed by
022: *  ZSYTRF.
023: *
024: *  Arguments
025: *  =========
026: *
027: *  UPLO    (input) CHARACTER*1
028: *          Specifies whether the details of the factorization are stored
029: *          as an upper or lower triangular matrix.
030: *          = 'U':  Upper triangular, form is A = U*D*U**T;
031: *          = 'L':  Lower triangular, form is A = L*D*L**T.
032: *
033: *  N       (input) INTEGER
034: *          The order of the matrix A.  N >= 0.
035: *
036: *  A       (input/output) COMPLEX*16 array, dimension (LDA,N)
037: *          On entry, the block diagonal matrix D and the multipliers
038: *          used to obtain the factor U or L as computed by ZSYTRF.
039: *
040: *          On exit, if INFO = 0, the (symmetric) inverse of the original
041: *          matrix.  If UPLO = 'U', the upper triangular part of the
042: *          inverse is formed and the part of A below the diagonal is not
043: *          referenced; if UPLO = 'L' the lower triangular part of the
044: *          inverse is formed and the part of A above the diagonal is
045: *          not referenced.
046: *
047: *  LDA     (input) INTEGER
048: *          The leading dimension of the array A.  LDA >= max(1,N).
049: *
050: *  IPIV    (input) INTEGER array, dimension (N)
051: *          Details of the interchanges and the block structure of D
052: *          as determined by ZSYTRF.
053: *
054: *  WORK    (workspace) COMPLEX*16 array, dimension (2*N)
055: *
056: *  INFO    (output) INTEGER
057: *          = 0: successful exit
058: *          < 0: if INFO = -i, the i-th argument had an illegal value
059: *          > 0: if INFO = i, D(i,i) = 0; the matrix is singular and its
060: *               inverse could not be computed.
061: *
062: *  =====================================================================
063: *
064: *     .. Parameters ..
065:       COMPLEX*16         ONE, ZERO
066:       PARAMETER          ( ONE = ( 1.0D+0, 0.0D+0 ),
067:      $                   ZERO = ( 0.0D+0, 0.0D+0 ) )
068: *     ..
069: *     .. Local Scalars ..
070:       LOGICAL            UPPER
071:       INTEGER            K, KP, KSTEP
072:       COMPLEX*16         AK, AKKP1, AKP1, D, T, TEMP
073: *     ..
074: *     .. External Functions ..
075:       LOGICAL            LSAME
076:       COMPLEX*16         ZDOTU
077:       EXTERNAL           LSAME, ZDOTU
078: *     ..
079: *     .. External Subroutines ..
080:       EXTERNAL           XERBLA, ZCOPY, ZSWAP, ZSYMV
081: *     ..
082: *     .. Intrinsic Functions ..
083:       INTRINSIC          ABS, MAX
084: *     ..
085: *     .. Executable Statements ..
086: *
087: *     Test the input parameters.
088: *
089:       INFO = 0
090:       UPPER = LSAME( UPLO, 'U' )
091:       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
092:          INFO = -1
093:       ELSE IF( N.LT.0 ) THEN
094:          INFO = -2
095:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
096:          INFO = -4
097:       END IF
098:       IF( INFO.NE.0 ) THEN
099:          CALL XERBLA( 'ZSYTRI', -INFO )
100:          RETURN
101:       END IF
102: *
103: *     Quick return if possible
104: *
105:       IF( N.EQ.0 )
106:      $   RETURN
107: *
108: *     Check that the diagonal matrix D is nonsingular.
109: *
110:       IF( UPPER ) THEN
111: *
112: *        Upper triangular storage: examine D from bottom to top
113: *
114:          DO 10 INFO = N, 1, -1
115:             IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO )
116:      $         RETURN
117:    10    CONTINUE
118:       ELSE
119: *
120: *        Lower triangular storage: examine D from top to bottom.
121: *
122:          DO 20 INFO = 1, N
123:             IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO )
124:      $         RETURN
125:    20    CONTINUE
126:       END IF
127:       INFO = 0
128: *
129:       IF( UPPER ) THEN
130: *
131: *        Compute inv(A) from the factorization A = U*D*U'.
132: *
133: *        K is the main loop index, increasing from 1 to N in steps of
134: *        1 or 2, depending on the size of the diagonal blocks.
135: *
136:          K = 1
137:    30    CONTINUE
138: *
139: *        If K > N, exit from loop.
140: *
141:          IF( K.GT.N )
142:      $      GO TO 40
143: *
144:          IF( IPIV( K ).GT.0 ) THEN
145: *
146: *           1 x 1 diagonal block
147: *
148: *           Invert the diagonal block.
149: *
150:             A( K, K ) = ONE / A( K, K )
151: *
152: *           Compute column K of the inverse.
153: *
154:             IF( K.GT.1 ) THEN
155:                CALL ZCOPY( K-1, A( 1, K ), 1, WORK, 1 )
156:                CALL ZSYMV( UPLO, K-1, -ONE, A, LDA, WORK, 1, ZERO,
157:      $                     A( 1, K ), 1 )
158:                A( K, K ) = A( K, K ) - ZDOTU( K-1, WORK, 1, A( 1, K ),
159:      $                     1 )
160:             END IF
161:             KSTEP = 1
162:          ELSE
163: *
164: *           2 x 2 diagonal block
165: *
166: *           Invert the diagonal block.
167: *
168:             T = A( K, K+1 )
169:             AK = A( K, K ) / T
170:             AKP1 = A( K+1, K+1 ) / T
171:             AKKP1 = A( K, K+1 ) / T
172:             D = T*( AK*AKP1-ONE )
173:             A( K, K ) = AKP1 / D
174:             A( K+1, K+1 ) = AK / D
175:             A( K, K+1 ) = -AKKP1 / D
176: *
177: *           Compute columns K and K+1 of the inverse.
178: *
179:             IF( K.GT.1 ) THEN
180:                CALL ZCOPY( K-1, A( 1, K ), 1, WORK, 1 )
181:                CALL ZSYMV( UPLO, K-1, -ONE, A, LDA, WORK, 1, ZERO,
182:      $                     A( 1, K ), 1 )
183:                A( K, K ) = A( K, K ) - ZDOTU( K-1, WORK, 1, A( 1, K ),
184:      $                     1 )
185:                A( K, K+1 ) = A( K, K+1 ) -
186:      $                       ZDOTU( K-1, A( 1, K ), 1, A( 1, K+1 ), 1 )
187:                CALL ZCOPY( K-1, A( 1, K+1 ), 1, WORK, 1 )
188:                CALL ZSYMV( UPLO, K-1, -ONE, A, LDA, WORK, 1, ZERO,
189:      $                     A( 1, K+1 ), 1 )
190:                A( K+1, K+1 ) = A( K+1, K+1 ) -
191:      $                         ZDOTU( K-1, WORK, 1, A( 1, K+1 ), 1 )
192:             END IF
193:             KSTEP = 2
194:          END IF
195: *
196:          KP = ABS( IPIV( K ) )
197:          IF( KP.NE.K ) THEN
198: *
199: *           Interchange rows and columns K and KP in the leading
200: *           submatrix A(1:k+1,1:k+1)
201: *
202:             CALL ZSWAP( KP-1, A( 1, K ), 1, A( 1, KP ), 1 )
203:             CALL ZSWAP( K-KP-1, A( KP+1, K ), 1, A( KP, KP+1 ), LDA )
204:             TEMP = A( K, K )
205:             A( K, K ) = A( KP, KP )
206:             A( KP, KP ) = TEMP
207:             IF( KSTEP.EQ.2 ) THEN
208:                TEMP = A( K, K+1 )
209:                A( K, K+1 ) = A( KP, K+1 )
210:                A( KP, K+1 ) = TEMP
211:             END IF
212:          END IF
213: *
214:          K = K + KSTEP
215:          GO TO 30
216:    40    CONTINUE
217: *
218:       ELSE
219: *
220: *        Compute inv(A) from the factorization A = L*D*L'.
221: *
222: *        K is the main loop index, increasing from 1 to N in steps of
223: *        1 or 2, depending on the size of the diagonal blocks.
224: *
225:          K = N
226:    50    CONTINUE
227: *
228: *        If K < 1, exit from loop.
229: *
230:          IF( K.LT.1 )
231:      $      GO TO 60
232: *
233:          IF( IPIV( K ).GT.0 ) THEN
234: *
235: *           1 x 1 diagonal block
236: *
237: *           Invert the diagonal block.
238: *
239:             A( K, K ) = ONE / A( K, K )
240: *
241: *           Compute column K of the inverse.
242: *
243:             IF( K.LT.N ) THEN
244:                CALL ZCOPY( N-K, A( K+1, K ), 1, WORK, 1 )
245:                CALL ZSYMV( UPLO, N-K, -ONE, A( K+1, K+1 ), LDA, WORK, 1,
246:      $                     ZERO, A( K+1, K ), 1 )
247:                A( K, K ) = A( K, K ) - ZDOTU( N-K, WORK, 1, A( K+1, K ),
248:      $                     1 )
249:             END IF
250:             KSTEP = 1
251:          ELSE
252: *
253: *           2 x 2 diagonal block
254: *
255: *           Invert the diagonal block.
256: *
257:             T = A( K, K-1 )
258:             AK = A( K-1, K-1 ) / T
259:             AKP1 = A( K, K ) / T
260:             AKKP1 = A( K, K-1 ) / T
261:             D = T*( AK*AKP1-ONE )
262:             A( K-1, K-1 ) = AKP1 / D
263:             A( K, K ) = AK / D
264:             A( K, K-1 ) = -AKKP1 / D
265: *
266: *           Compute columns K-1 and K of the inverse.
267: *
268:             IF( K.LT.N ) THEN
269:                CALL ZCOPY( N-K, A( K+1, K ), 1, WORK, 1 )
270:                CALL ZSYMV( UPLO, N-K, -ONE, A( K+1, K+1 ), LDA, WORK, 1,
271:      $                     ZERO, A( K+1, K ), 1 )
272:                A( K, K ) = A( K, K ) - ZDOTU( N-K, WORK, 1, A( K+1, K ),
273:      $                     1 )
274:                A( K, K-1 ) = A( K, K-1 ) -
275:      $                       ZDOTU( N-K, A( K+1, K ), 1, A( K+1, K-1 ),
276:      $                       1 )
277:                CALL ZCOPY( N-K, A( K+1, K-1 ), 1, WORK, 1 )
278:                CALL ZSYMV( UPLO, N-K, -ONE, A( K+1, K+1 ), LDA, WORK, 1,
279:      $                     ZERO, A( K+1, K-1 ), 1 )
280:                A( K-1, K-1 ) = A( K-1, K-1 ) -
281:      $                         ZDOTU( N-K, WORK, 1, A( K+1, K-1 ), 1 )
282:             END IF
283:             KSTEP = 2
284:          END IF
285: *
286:          KP = ABS( IPIV( K ) )
287:          IF( KP.NE.K ) THEN
288: *
289: *           Interchange rows and columns K and KP in the trailing
290: *           submatrix A(k-1:n,k-1:n)
291: *
292:             IF( KP.LT.N )
293:      $         CALL ZSWAP( N-KP, A( KP+1, K ), 1, A( KP+1, KP ), 1 )
294:             CALL ZSWAP( KP-K-1, A( K+1, K ), 1, A( KP, K+1 ), LDA )
295:             TEMP = A( K, K )
296:             A( K, K ) = A( KP, KP )
297:             A( KP, KP ) = TEMP
298:             IF( KSTEP.EQ.2 ) THEN
299:                TEMP = A( K, K-1 )
300:                A( K, K-1 ) = A( KP, K-1 )
301:                A( KP, K-1 ) = TEMP
302:             END IF
303:          END IF
304: *
305:          K = K - KSTEP
306:          GO TO 50
307:    60    CONTINUE
308:       END IF
309: *
310:       RETURN
311: *
312: *     End of ZSYTRI
313: *
314:       END
315: