```001:       SUBROUTINE DSYTRI( 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:       DOUBLE PRECISION   A( LDA, * ), WORK( * )
015: *     ..
016: *
017: *  Purpose
018: *  =======
019: *
020: *  DSYTRI computes the inverse of a real symmetric indefinite matrix
021: *  A using the factorization A = U*D*U**T or A = L*D*L**T computed by
022: *  DSYTRF.
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) DOUBLE PRECISION 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 DSYTRF.
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 DSYTRF.
053: *
054: *  WORK    (workspace) DOUBLE PRECISION array, dimension (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:       DOUBLE PRECISION   ONE, ZERO
066:       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
067: *     ..
068: *     .. Local Scalars ..
069:       LOGICAL            UPPER
070:       INTEGER            K, KP, KSTEP
071:       DOUBLE PRECISION   AK, AKKP1, AKP1, D, T, TEMP
072: *     ..
073: *     .. External Functions ..
074:       LOGICAL            LSAME
075:       DOUBLE PRECISION   DDOT
076:       EXTERNAL           LSAME, DDOT
077: *     ..
078: *     .. External Subroutines ..
079:       EXTERNAL           DCOPY, DSWAP, DSYMV, XERBLA
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( 'DSYTRI', -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 DCOPY( K-1, A( 1, K ), 1, WORK, 1 )
155:                CALL DSYMV( UPLO, K-1, -ONE, A, LDA, WORK, 1, ZERO,
156:      \$                     A( 1, K ), 1 )
157:                A( K, K ) = A( K, K ) - DDOT( 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 = ABS( 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 DCOPY( K-1, A( 1, K ), 1, WORK, 1 )
180:                CALL DSYMV( UPLO, K-1, -ONE, A, LDA, WORK, 1, ZERO,
181:      \$                     A( 1, K ), 1 )
182:                A( K, K ) = A( K, K ) - DDOT( K-1, WORK, 1, A( 1, K ),
183:      \$                     1 )
184:                A( K, K+1 ) = A( K, K+1 ) -
185:      \$                       DDOT( K-1, A( 1, K ), 1, A( 1, K+1 ), 1 )
186:                CALL DCOPY( K-1, A( 1, K+1 ), 1, WORK, 1 )
187:                CALL DSYMV( 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:      \$                         DDOT( 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 DSWAP( KP-1, A( 1, K ), 1, A( 1, KP ), 1 )
202:             CALL DSWAP( 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 DCOPY( N-K, A( K+1, K ), 1, WORK, 1 )
244:                CALL DSYMV( 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 ) - DDOT( 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 = ABS( 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 DCOPY( N-K, A( K+1, K ), 1, WORK, 1 )
269:                CALL DSYMV( 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 ) - DDOT( N-K, WORK, 1, A( K+1, K ),
272:      \$                     1 )
273:                A( K, K-1 ) = A( K, K-1 ) -
274:      \$                       DDOT( N-K, A( K+1, K ), 1, A( K+1, K-1 ),
275:      \$                       1 )
276:                CALL DCOPY( N-K, A( K+1, K-1 ), 1, WORK, 1 )
277:                CALL DSYMV( 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:      \$                         DDOT( 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 DSWAP( N-KP, A( KP+1, K ), 1, A( KP+1, KP ), 1 )
293:             CALL DSWAP( 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 DSYTRI
312: *
313:       END
314: ```