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