001:       SUBROUTINE DSPTRI( 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:       DOUBLE PRECISION   AP( * ), WORK( * )
014: *     ..
015: *
016: *  Purpose
017: *  =======
018: *
019: *  DSPTRI computes the inverse of a real 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 DSPTRF.
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) DOUBLE PRECISION 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 DSPTRF,
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 DSPTRF.
050: *
051: *  WORK    (workspace) DOUBLE PRECISION 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:       DOUBLE PRECISION   ONE, ZERO
063:       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
064: *     ..
065: *     .. Local Scalars ..
066:       LOGICAL            UPPER
067:       INTEGER            J, K, KC, KCNEXT, KP, KPC, KSTEP, KX, NPP
068:       DOUBLE PRECISION   AK, AKKP1, AKP1, D, T, TEMP
069: *     ..
070: *     .. External Functions ..
071:       LOGICAL            LSAME
072:       DOUBLE PRECISION   DDOT
073:       EXTERNAL           LSAME, DDOT
074: *     ..
075: *     .. External Subroutines ..
076:       EXTERNAL           DCOPY, DSPMV, DSWAP, XERBLA
077: *     ..
078: *     .. Intrinsic Functions ..
079:       INTRINSIC          ABS
080: *     ..
081: *     .. Executable Statements ..
082: *
083: *     Test the input parameters.
084: *
085:       INFO = 0
086:       UPPER = LSAME( UPLO, 'U' )
087:       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
088:          INFO = -1
089:       ELSE IF( N.LT.0 ) THEN
090:          INFO = -2
091:       END IF
092:       IF( INFO.NE.0 ) THEN
093:          CALL XERBLA( 'DSPTRI', -INFO )
094:          RETURN
095:       END IF
096: *
097: *     Quick return if possible
098: *
099:       IF( N.EQ.0 )
100:      $   RETURN
101: *
102: *     Check that the diagonal matrix D is nonsingular.
103: *
104:       IF( UPPER ) THEN
105: *
106: *        Upper triangular storage: examine D from bottom to top
107: *
108:          KP = N*( N+1 ) / 2
109:          DO 10 INFO = N, 1, -1
110:             IF( IPIV( INFO ).GT.0 .AND. AP( KP ).EQ.ZERO )
111:      $         RETURN
112:             KP = KP - INFO
113:    10    CONTINUE
114:       ELSE
115: *
116: *        Lower triangular storage: examine D from top to bottom.
117: *
118:          KP = 1
119:          DO 20 INFO = 1, N
120:             IF( IPIV( INFO ).GT.0 .AND. AP( KP ).EQ.ZERO )
121:      $         RETURN
122:             KP = KP + N - INFO + 1
123:    20    CONTINUE
124:       END IF
125:       INFO = 0
126: *
127:       IF( UPPER ) THEN
128: *
129: *        Compute inv(A) from the factorization A = U*D*U'.
130: *
131: *        K is the main loop index, increasing from 1 to N in steps of
132: *        1 or 2, depending on the size of the diagonal blocks.
133: *
134:          K = 1
135:          KC = 1
136:    30    CONTINUE
137: *
138: *        If K > N, exit from loop.
139: *
140:          IF( K.GT.N )
141:      $      GO TO 50
142: *
143:          KCNEXT = KC + K
144:          IF( IPIV( K ).GT.0 ) THEN
145: *
146: *           1 x 1 diagonal block
147: *
148: *           Invert the diagonal block.
149: *
150:             AP( KC+K-1 ) = ONE / AP( KC+K-1 )
151: *
152: *           Compute column K of the inverse.
153: *
154:             IF( K.GT.1 ) THEN
155:                CALL DCOPY( K-1, AP( KC ), 1, WORK, 1 )
156:                CALL DSPMV( UPLO, K-1, -ONE, AP, WORK, 1, ZERO, AP( KC ),
157:      $                     1 )
158:                AP( KC+K-1 ) = AP( KC+K-1 ) -
159:      $                        DDOT( K-1, WORK, 1, AP( KC ), 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 = ABS( AP( KCNEXT+K-1 ) )
169:             AK = AP( KC+K-1 ) / T
170:             AKP1 = AP( KCNEXT+K ) / T
171:             AKKP1 = AP( KCNEXT+K-1 ) / T
172:             D = T*( AK*AKP1-ONE )
173:             AP( KC+K-1 ) = AKP1 / D
174:             AP( KCNEXT+K ) = AK / D
175:             AP( KCNEXT+K-1 ) = -AKKP1 / D
176: *
177: *           Compute columns K and K+1 of the inverse.
178: *
179:             IF( K.GT.1 ) THEN
180:                CALL DCOPY( K-1, AP( KC ), 1, WORK, 1 )
181:                CALL DSPMV( UPLO, K-1, -ONE, AP, WORK, 1, ZERO, AP( KC ),
182:      $                     1 )
183:                AP( KC+K-1 ) = AP( KC+K-1 ) -
184:      $                        DDOT( K-1, WORK, 1, AP( KC ), 1 )
185:                AP( KCNEXT+K-1 ) = AP( KCNEXT+K-1 ) -
186:      $                            DDOT( K-1, AP( KC ), 1, AP( KCNEXT ),
187:      $                            1 )
188:                CALL DCOPY( K-1, AP( KCNEXT ), 1, WORK, 1 )
189:                CALL DSPMV( UPLO, K-1, -ONE, AP, WORK, 1, ZERO,
190:      $                     AP( KCNEXT ), 1 )
191:                AP( KCNEXT+K ) = AP( KCNEXT+K ) -
192:      $                          DDOT( K-1, WORK, 1, AP( KCNEXT ), 1 )
193:             END IF
194:             KSTEP = 2
195:             KCNEXT = KCNEXT + K + 1
196:          END IF
197: *
198:          KP = ABS( IPIV( K ) )
199:          IF( KP.NE.K ) THEN
200: *
201: *           Interchange rows and columns K and KP in the leading
202: *           submatrix A(1:k+1,1:k+1)
203: *
204:             KPC = ( KP-1 )*KP / 2 + 1
205:             CALL DSWAP( KP-1, AP( KC ), 1, AP( KPC ), 1 )
206:             KX = KPC + KP - 1
207:             DO 40 J = KP + 1, K - 1
208:                KX = KX + J - 1
209:                TEMP = AP( KC+J-1 )
210:                AP( KC+J-1 ) = AP( KX )
211:                AP( KX ) = TEMP
212:    40       CONTINUE
213:             TEMP = AP( KC+K-1 )
214:             AP( KC+K-1 ) = AP( KPC+KP-1 )
215:             AP( KPC+KP-1 ) = TEMP
216:             IF( KSTEP.EQ.2 ) THEN
217:                TEMP = AP( KC+K+K-1 )
218:                AP( KC+K+K-1 ) = AP( KC+K+KP-1 )
219:                AP( KC+K+KP-1 ) = TEMP
220:             END IF
221:          END IF
222: *
223:          K = K + KSTEP
224:          KC = KCNEXT
225:          GO TO 30
226:    50    CONTINUE
227: *
228:       ELSE
229: *
230: *        Compute inv(A) from the factorization A = L*D*L'.
231: *
232: *        K is the main loop index, increasing from 1 to N in steps of
233: *        1 or 2, depending on the size of the diagonal blocks.
234: *
235:          NPP = N*( N+1 ) / 2
236:          K = N
237:          KC = NPP
238:    60    CONTINUE
239: *
240: *        If K < 1, exit from loop.
241: *
242:          IF( K.LT.1 )
243:      $      GO TO 80
244: *
245:          KCNEXT = KC - ( N-K+2 )
246:          IF( IPIV( K ).GT.0 ) THEN
247: *
248: *           1 x 1 diagonal block
249: *
250: *           Invert the diagonal block.
251: *
252:             AP( KC ) = ONE / AP( KC )
253: *
254: *           Compute column K of the inverse.
255: *
256:             IF( K.LT.N ) THEN
257:                CALL DCOPY( N-K, AP( KC+1 ), 1, WORK, 1 )
258:                CALL DSPMV( UPLO, N-K, -ONE, AP( KC+N-K+1 ), WORK, 1,
259:      $                     ZERO, AP( KC+1 ), 1 )
260:                AP( KC ) = AP( KC ) - DDOT( N-K, WORK, 1, AP( KC+1 ), 1 )
261:             END IF
262:             KSTEP = 1
263:          ELSE
264: *
265: *           2 x 2 diagonal block
266: *
267: *           Invert the diagonal block.
268: *
269:             T = ABS( AP( KCNEXT+1 ) )
270:             AK = AP( KCNEXT ) / T
271:             AKP1 = AP( KC ) / T
272:             AKKP1 = AP( KCNEXT+1 ) / T
273:             D = T*( AK*AKP1-ONE )
274:             AP( KCNEXT ) = AKP1 / D
275:             AP( KC ) = AK / D
276:             AP( KCNEXT+1 ) = -AKKP1 / D
277: *
278: *           Compute columns K-1 and K of the inverse.
279: *
280:             IF( K.LT.N ) THEN
281:                CALL DCOPY( N-K, AP( KC+1 ), 1, WORK, 1 )
282:                CALL DSPMV( UPLO, N-K, -ONE, AP( KC+( N-K+1 ) ), WORK, 1,
283:      $                     ZERO, AP( KC+1 ), 1 )
284:                AP( KC ) = AP( KC ) - DDOT( N-K, WORK, 1, AP( KC+1 ), 1 )
285:                AP( KCNEXT+1 ) = AP( KCNEXT+1 ) -
286:      $                          DDOT( N-K, AP( KC+1 ), 1,
287:      $                          AP( KCNEXT+2 ), 1 )
288:                CALL DCOPY( N-K, AP( KCNEXT+2 ), 1, WORK, 1 )
289:                CALL DSPMV( UPLO, N-K, -ONE, AP( KC+( N-K+1 ) ), WORK, 1,
290:      $                     ZERO, AP( KCNEXT+2 ), 1 )
291:                AP( KCNEXT ) = AP( KCNEXT ) -
292:      $                        DDOT( N-K, WORK, 1, AP( KCNEXT+2 ), 1 )
293:             END IF
294:             KSTEP = 2
295:             KCNEXT = KCNEXT - ( N-K+3 )
296:          END IF
297: *
298:          KP = ABS( IPIV( K ) )
299:          IF( KP.NE.K ) THEN
300: *
301: *           Interchange rows and columns K and KP in the trailing
302: *           submatrix A(k-1:n,k-1:n)
303: *
304:             KPC = NPP - ( N-KP+1 )*( N-KP+2 ) / 2 + 1
305:             IF( KP.LT.N )
306:      $         CALL DSWAP( N-KP, AP( KC+KP-K+1 ), 1, AP( KPC+1 ), 1 )
307:             KX = KC + KP - K
308:             DO 70 J = K + 1, KP - 1
309:                KX = KX + N - J + 1
310:                TEMP = AP( KC+J-K )
311:                AP( KC+J-K ) = AP( KX )
312:                AP( KX ) = TEMP
313:    70       CONTINUE
314:             TEMP = AP( KC )
315:             AP( KC ) = AP( KPC )
316:             AP( KPC ) = TEMP
317:             IF( KSTEP.EQ.2 ) THEN
318:                TEMP = AP( KC-N+K-1 )
319:                AP( KC-N+K-1 ) = AP( KC-N+KP-1 )
320:                AP( KC-N+KP-1 ) = TEMP
321:             END IF
322:          END IF
323: *
324:          K = K - KSTEP
325:          KC = KCNEXT
326:          GO TO 60
327:    80    CONTINUE
328:       END IF
329: *
330:       RETURN
331: *
332: *     End of DSPTRI
333: *
334:       END
335: