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