1 SUBROUTINE DSPTRI( UPLO, N, AP, IPIV, WORK, INFO )
2 *
3 * -- LAPACK routine (version 3.3.1) --
4 * -- LAPACK is a software package provided by Univ. of Tennessee, --
5 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6 * -- April 2011 --
7 *
8 * .. Scalar Arguments ..
9 CHARACTER UPLO
10 INTEGER INFO, N
11 * ..
12 * .. Array Arguments ..
13 INTEGER IPIV( * )
14 DOUBLE PRECISION AP( * ), WORK( * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * DSPTRI computes the inverse of a real symmetric indefinite matrix
21 * A in packed storage using the factorization A = U*D*U**T or
22 * A = L*D*L**T computed by DSPTRF.
23 *
24 * Arguments
25 * =========
26 *
27 * UPLO (input) CHARACTER*1
28 * Specifies whether the details of the factorization are stored
29 * as an upper or lower triangular matrix.
30 * = 'U': Upper triangular, form is A = U*D*U**T;
31 * = 'L': Lower triangular, form is A = L*D*L**T.
32 *
33 * N (input) INTEGER
34 * The order of the matrix A. N >= 0.
35 *
36 * AP (input/output) DOUBLE PRECISION array, dimension (N*(N+1)/2)
37 * On entry, the block diagonal matrix D and the multipliers
38 * used to obtain the factor U or L as computed by DSPTRF,
39 * stored as a packed triangular matrix.
40 *
41 * On exit, if INFO = 0, the (symmetric) inverse of the original
42 * matrix, stored as a packed triangular matrix. The j-th column
43 * of inv(A) is stored in the array AP as follows:
44 * if UPLO = 'U', AP(i + (j-1)*j/2) = inv(A)(i,j) for 1<=i<=j;
45 * if UPLO = 'L',
46 * AP(i + (j-1)*(2n-j)/2) = inv(A)(i,j) for j<=i<=n.
47 *
48 * IPIV (input) INTEGER array, dimension (N)
49 * Details of the interchanges and the block structure of D
50 * as determined by DSPTRF.
51 *
52 * WORK (workspace) DOUBLE PRECISION array, dimension (N)
53 *
54 * INFO (output) INTEGER
55 * = 0: successful exit
56 * < 0: if INFO = -i, the i-th argument had an illegal value
57 * > 0: if INFO = i, D(i,i) = 0; the matrix is singular and its
58 * inverse could not be computed.
59 *
60 * =====================================================================
61 *
62 * .. Parameters ..
63 DOUBLE PRECISION ONE, ZERO
64 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
65 * ..
66 * .. Local Scalars ..
67 LOGICAL UPPER
68 INTEGER J, K, KC, KCNEXT, KP, KPC, KSTEP, KX, NPP
69 DOUBLE PRECISION AK, AKKP1, AKP1, D, T, TEMP
70 * ..
71 * .. External Functions ..
72 LOGICAL LSAME
73 DOUBLE PRECISION DDOT
74 EXTERNAL LSAME, DDOT
75 * ..
76 * .. External Subroutines ..
77 EXTERNAL DCOPY, DSPMV, DSWAP, XERBLA
78 * ..
79 * .. Intrinsic Functions ..
80 INTRINSIC ABS
81 * ..
82 * .. Executable Statements ..
83 *
84 * Test the input parameters.
85 *
86 INFO = 0
87 UPPER = LSAME( UPLO, 'U' )
88 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
89 INFO = -1
90 ELSE IF( N.LT.0 ) THEN
91 INFO = -2
92 END IF
93 IF( INFO.NE.0 ) THEN
94 CALL XERBLA( 'DSPTRI', -INFO )
95 RETURN
96 END IF
97 *
98 * Quick return if possible
99 *
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**T.
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 DCOPY( K-1, AP( KC ), 1, WORK, 1 )
157 CALL DSPMV( UPLO, K-1, -ONE, AP, WORK, 1, ZERO, AP( KC ),
158 $ 1 )
159 AP( KC+K-1 ) = AP( KC+K-1 ) -
160 $ DDOT( 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 = ABS( 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 DCOPY( K-1, AP( KC ), 1, WORK, 1 )
182 CALL DSPMV( UPLO, K-1, -ONE, AP, WORK, 1, ZERO, AP( KC ),
183 $ 1 )
184 AP( KC+K-1 ) = AP( KC+K-1 ) -
185 $ DDOT( K-1, WORK, 1, AP( KC ), 1 )
186 AP( KCNEXT+K-1 ) = AP( KCNEXT+K-1 ) -
187 $ DDOT( K-1, AP( KC ), 1, AP( KCNEXT ),
188 $ 1 )
189 CALL DCOPY( K-1, AP( KCNEXT ), 1, WORK, 1 )
190 CALL DSPMV( UPLO, K-1, -ONE, AP, WORK, 1, ZERO,
191 $ AP( KCNEXT ), 1 )
192 AP( KCNEXT+K ) = AP( KCNEXT+K ) -
193 $ DDOT( 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 DSWAP( 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**T.
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 DCOPY( N-K, AP( KC+1 ), 1, WORK, 1 )
259 CALL DSPMV( UPLO, N-K, -ONE, AP( KC+N-K+1 ), WORK, 1,
260 $ ZERO, AP( KC+1 ), 1 )
261 AP( KC ) = AP( KC ) - DDOT( N-K, WORK, 1, AP( KC+1 ), 1 )
262 END IF
263 KSTEP = 1
264 ELSE
265 *
266 * 2 x 2 diagonal block
267 *
268 * Invert the diagonal block.
269 *
270 T = ABS( AP( KCNEXT+1 ) )
271 AK = AP( KCNEXT ) / T
272 AKP1 = AP( KC ) / T
273 AKKP1 = AP( KCNEXT+1 ) / T
274 D = T*( AK*AKP1-ONE )
275 AP( KCNEXT ) = AKP1 / D
276 AP( KC ) = AK / D
277 AP( KCNEXT+1 ) = -AKKP1 / D
278 *
279 * Compute columns K-1 and K of the inverse.
280 *
281 IF( K.LT.N ) THEN
282 CALL DCOPY( N-K, AP( KC+1 ), 1, WORK, 1 )
283 CALL DSPMV( UPLO, N-K, -ONE, AP( KC+( N-K+1 ) ), WORK, 1,
284 $ ZERO, AP( KC+1 ), 1 )
285 AP( KC ) = AP( KC ) - DDOT( N-K, WORK, 1, AP( KC+1 ), 1 )
286 AP( KCNEXT+1 ) = AP( KCNEXT+1 ) -
287 $ DDOT( N-K, AP( KC+1 ), 1,
288 $ AP( KCNEXT+2 ), 1 )
289 CALL DCOPY( N-K, AP( KCNEXT+2 ), 1, WORK, 1 )
290 CALL DSPMV( UPLO, N-K, -ONE, AP( KC+( N-K+1 ) ), WORK, 1,
291 $ ZERO, AP( KCNEXT+2 ), 1 )
292 AP( KCNEXT ) = AP( KCNEXT ) -
293 $ DDOT( N-K, WORK, 1, AP( KCNEXT+2 ), 1 )
294 END IF
295 KSTEP = 2
296 KCNEXT = KCNEXT - ( N-K+3 )
297 END IF
298 *
299 KP = ABS( IPIV( K ) )
300 IF( KP.NE.K ) THEN
301 *
302 * Interchange rows and columns K and KP in the trailing
303 * submatrix A(k-1:n,k-1:n)
304 *
305 KPC = NPP - ( N-KP+1 )*( N-KP+2 ) / 2 + 1
306 IF( KP.LT.N )
307 $ CALL DSWAP( N-KP, AP( KC+KP-K+1 ), 1, AP( KPC+1 ), 1 )
308 KX = KC + KP - K
309 DO 70 J = K + 1, KP - 1
310 KX = KX + N - J + 1
311 TEMP = AP( KC+J-K )
312 AP( KC+J-K ) = AP( KX )
313 AP( KX ) = TEMP
314 70 CONTINUE
315 TEMP = AP( KC )
316 AP( KC ) = AP( KPC )
317 AP( KPC ) = TEMP
318 IF( KSTEP.EQ.2 ) THEN
319 TEMP = AP( KC-N+K-1 )
320 AP( KC-N+K-1 ) = AP( KC-N+KP-1 )
321 AP( KC-N+KP-1 ) = TEMP
322 END IF
323 END IF
324 *
325 K = K - KSTEP
326 KC = KCNEXT
327 GO TO 60
328 80 CONTINUE
329 END IF
330 *
331 RETURN
332 *
333 * End of DSPTRI
334 *
335 END
2 *
3 * -- LAPACK routine (version 3.3.1) --
4 * -- LAPACK is a software package provided by Univ. of Tennessee, --
5 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6 * -- April 2011 --
7 *
8 * .. Scalar Arguments ..
9 CHARACTER UPLO
10 INTEGER INFO, N
11 * ..
12 * .. Array Arguments ..
13 INTEGER IPIV( * )
14 DOUBLE PRECISION AP( * ), WORK( * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * DSPTRI computes the inverse of a real symmetric indefinite matrix
21 * A in packed storage using the factorization A = U*D*U**T or
22 * A = L*D*L**T computed by DSPTRF.
23 *
24 * Arguments
25 * =========
26 *
27 * UPLO (input) CHARACTER*1
28 * Specifies whether the details of the factorization are stored
29 * as an upper or lower triangular matrix.
30 * = 'U': Upper triangular, form is A = U*D*U**T;
31 * = 'L': Lower triangular, form is A = L*D*L**T.
32 *
33 * N (input) INTEGER
34 * The order of the matrix A. N >= 0.
35 *
36 * AP (input/output) DOUBLE PRECISION array, dimension (N*(N+1)/2)
37 * On entry, the block diagonal matrix D and the multipliers
38 * used to obtain the factor U or L as computed by DSPTRF,
39 * stored as a packed triangular matrix.
40 *
41 * On exit, if INFO = 0, the (symmetric) inverse of the original
42 * matrix, stored as a packed triangular matrix. The j-th column
43 * of inv(A) is stored in the array AP as follows:
44 * if UPLO = 'U', AP(i + (j-1)*j/2) = inv(A)(i,j) for 1<=i<=j;
45 * if UPLO = 'L',
46 * AP(i + (j-1)*(2n-j)/2) = inv(A)(i,j) for j<=i<=n.
47 *
48 * IPIV (input) INTEGER array, dimension (N)
49 * Details of the interchanges and the block structure of D
50 * as determined by DSPTRF.
51 *
52 * WORK (workspace) DOUBLE PRECISION array, dimension (N)
53 *
54 * INFO (output) INTEGER
55 * = 0: successful exit
56 * < 0: if INFO = -i, the i-th argument had an illegal value
57 * > 0: if INFO = i, D(i,i) = 0; the matrix is singular and its
58 * inverse could not be computed.
59 *
60 * =====================================================================
61 *
62 * .. Parameters ..
63 DOUBLE PRECISION ONE, ZERO
64 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
65 * ..
66 * .. Local Scalars ..
67 LOGICAL UPPER
68 INTEGER J, K, KC, KCNEXT, KP, KPC, KSTEP, KX, NPP
69 DOUBLE PRECISION AK, AKKP1, AKP1, D, T, TEMP
70 * ..
71 * .. External Functions ..
72 LOGICAL LSAME
73 DOUBLE PRECISION DDOT
74 EXTERNAL LSAME, DDOT
75 * ..
76 * .. External Subroutines ..
77 EXTERNAL DCOPY, DSPMV, DSWAP, XERBLA
78 * ..
79 * .. Intrinsic Functions ..
80 INTRINSIC ABS
81 * ..
82 * .. Executable Statements ..
83 *
84 * Test the input parameters.
85 *
86 INFO = 0
87 UPPER = LSAME( UPLO, 'U' )
88 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
89 INFO = -1
90 ELSE IF( N.LT.0 ) THEN
91 INFO = -2
92 END IF
93 IF( INFO.NE.0 ) THEN
94 CALL XERBLA( 'DSPTRI', -INFO )
95 RETURN
96 END IF
97 *
98 * Quick return if possible
99 *
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**T.
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 DCOPY( K-1, AP( KC ), 1, WORK, 1 )
157 CALL DSPMV( UPLO, K-1, -ONE, AP, WORK, 1, ZERO, AP( KC ),
158 $ 1 )
159 AP( KC+K-1 ) = AP( KC+K-1 ) -
160 $ DDOT( 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 = ABS( 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 DCOPY( K-1, AP( KC ), 1, WORK, 1 )
182 CALL DSPMV( UPLO, K-1, -ONE, AP, WORK, 1, ZERO, AP( KC ),
183 $ 1 )
184 AP( KC+K-1 ) = AP( KC+K-1 ) -
185 $ DDOT( K-1, WORK, 1, AP( KC ), 1 )
186 AP( KCNEXT+K-1 ) = AP( KCNEXT+K-1 ) -
187 $ DDOT( K-1, AP( KC ), 1, AP( KCNEXT ),
188 $ 1 )
189 CALL DCOPY( K-1, AP( KCNEXT ), 1, WORK, 1 )
190 CALL DSPMV( UPLO, K-1, -ONE, AP, WORK, 1, ZERO,
191 $ AP( KCNEXT ), 1 )
192 AP( KCNEXT+K ) = AP( KCNEXT+K ) -
193 $ DDOT( 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 DSWAP( 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**T.
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 DCOPY( N-K, AP( KC+1 ), 1, WORK, 1 )
259 CALL DSPMV( UPLO, N-K, -ONE, AP( KC+N-K+1 ), WORK, 1,
260 $ ZERO, AP( KC+1 ), 1 )
261 AP( KC ) = AP( KC ) - DDOT( N-K, WORK, 1, AP( KC+1 ), 1 )
262 END IF
263 KSTEP = 1
264 ELSE
265 *
266 * 2 x 2 diagonal block
267 *
268 * Invert the diagonal block.
269 *
270 T = ABS( AP( KCNEXT+1 ) )
271 AK = AP( KCNEXT ) / T
272 AKP1 = AP( KC ) / T
273 AKKP1 = AP( KCNEXT+1 ) / T
274 D = T*( AK*AKP1-ONE )
275 AP( KCNEXT ) = AKP1 / D
276 AP( KC ) = AK / D
277 AP( KCNEXT+1 ) = -AKKP1 / D
278 *
279 * Compute columns K-1 and K of the inverse.
280 *
281 IF( K.LT.N ) THEN
282 CALL DCOPY( N-K, AP( KC+1 ), 1, WORK, 1 )
283 CALL DSPMV( UPLO, N-K, -ONE, AP( KC+( N-K+1 ) ), WORK, 1,
284 $ ZERO, AP( KC+1 ), 1 )
285 AP( KC ) = AP( KC ) - DDOT( N-K, WORK, 1, AP( KC+1 ), 1 )
286 AP( KCNEXT+1 ) = AP( KCNEXT+1 ) -
287 $ DDOT( N-K, AP( KC+1 ), 1,
288 $ AP( KCNEXT+2 ), 1 )
289 CALL DCOPY( N-K, AP( KCNEXT+2 ), 1, WORK, 1 )
290 CALL DSPMV( UPLO, N-K, -ONE, AP( KC+( N-K+1 ) ), WORK, 1,
291 $ ZERO, AP( KCNEXT+2 ), 1 )
292 AP( KCNEXT ) = AP( KCNEXT ) -
293 $ DDOT( N-K, WORK, 1, AP( KCNEXT+2 ), 1 )
294 END IF
295 KSTEP = 2
296 KCNEXT = KCNEXT - ( N-K+3 )
297 END IF
298 *
299 KP = ABS( IPIV( K ) )
300 IF( KP.NE.K ) THEN
301 *
302 * Interchange rows and columns K and KP in the trailing
303 * submatrix A(k-1:n,k-1:n)
304 *
305 KPC = NPP - ( N-KP+1 )*( N-KP+2 ) / 2 + 1
306 IF( KP.LT.N )
307 $ CALL DSWAP( N-KP, AP( KC+KP-K+1 ), 1, AP( KPC+1 ), 1 )
308 KX = KC + KP - K
309 DO 70 J = K + 1, KP - 1
310 KX = KX + N - J + 1
311 TEMP = AP( KC+J-K )
312 AP( KC+J-K ) = AP( KX )
313 AP( KX ) = TEMP
314 70 CONTINUE
315 TEMP = AP( KC )
316 AP( KC ) = AP( KPC )
317 AP( KPC ) = TEMP
318 IF( KSTEP.EQ.2 ) THEN
319 TEMP = AP( KC-N+K-1 )
320 AP( KC-N+K-1 ) = AP( KC-N+KP-1 )
321 AP( KC-N+KP-1 ) = TEMP
322 END IF
323 END IF
324 *
325 K = K - KSTEP
326 KC = KCNEXT
327 GO TO 60
328 80 CONTINUE
329 END IF
330 *
331 RETURN
332 *
333 * End of DSPTRI
334 *
335 END