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