1 SUBROUTINE DSPTRS( UPLO, N, NRHS, AP, IPIV, B, LDB, 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, LDB, N, NRHS
11 * ..
12 * .. Array Arguments ..
13 INTEGER IPIV( * )
14 DOUBLE PRECISION AP( * ), B( LDB, * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * DSPTRS solves a system of linear equations A*X = B with a real
21 * symmetric matrix A stored in packed format using the factorization
22 * A = U*D*U**T or 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 * NRHS (input) INTEGER
37 * The number of right hand sides, i.e., the number of columns
38 * of the matrix B. NRHS >= 0.
39 *
40 * AP (input) DOUBLE PRECISION array, dimension (N*(N+1)/2)
41 * The block diagonal matrix D and the multipliers used to
42 * obtain the factor U or L as computed by DSPTRF, stored as a
43 * packed triangular matrix.
44 *
45 * IPIV (input) INTEGER array, dimension (N)
46 * Details of the interchanges and the block structure of D
47 * as determined by DSPTRF.
48 *
49 * B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
50 * On entry, the right hand side matrix B.
51 * On exit, the solution matrix X.
52 *
53 * LDB (input) INTEGER
54 * The leading dimension of the array B. LDB >= max(1,N).
55 *
56 * INFO (output) INTEGER
57 * = 0: successful exit
58 * < 0: if INFO = -i, the i-th argument had an illegal value
59 *
60 * =====================================================================
61 *
62 * .. Parameters ..
63 DOUBLE PRECISION ONE
64 PARAMETER ( ONE = 1.0D+0 )
65 * ..
66 * .. Local Scalars ..
67 LOGICAL UPPER
68 INTEGER J, K, KC, KP
69 DOUBLE PRECISION AK, AKM1, AKM1K, BK, BKM1, DENOM
70 * ..
71 * .. External Functions ..
72 LOGICAL LSAME
73 EXTERNAL LSAME
74 * ..
75 * .. External Subroutines ..
76 EXTERNAL DGEMV, DGER, DSCAL, DSWAP, XERBLA
77 * ..
78 * .. Intrinsic Functions ..
79 INTRINSIC MAX
80 * ..
81 * .. Executable Statements ..
82 *
83 INFO = 0
84 UPPER = LSAME( UPLO, 'U' )
85 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
86 INFO = -1
87 ELSE IF( N.LT.0 ) THEN
88 INFO = -2
89 ELSE IF( NRHS.LT.0 ) THEN
90 INFO = -3
91 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
92 INFO = -7
93 END IF
94 IF( INFO.NE.0 ) THEN
95 CALL XERBLA( 'DSPTRS', -INFO )
96 RETURN
97 END IF
98 *
99 * Quick return if possible
100 *
101 IF( N.EQ.0 .OR. NRHS.EQ.0 )
102 $ RETURN
103 *
104 IF( UPPER ) THEN
105 *
106 * Solve A*X = B, where A = U*D*U**T.
107 *
108 * First solve U*D*X = B, overwriting B with X.
109 *
110 * K is the main loop index, decreasing from N to 1 in steps of
111 * 1 or 2, depending on the size of the diagonal blocks.
112 *
113 K = N
114 KC = N*( N+1 ) / 2 + 1
115 10 CONTINUE
116 *
117 * If K < 1, exit from loop.
118 *
119 IF( K.LT.1 )
120 $ GO TO 30
121 *
122 KC = KC - K
123 IF( IPIV( K ).GT.0 ) THEN
124 *
125 * 1 x 1 diagonal block
126 *
127 * Interchange rows K and IPIV(K).
128 *
129 KP = IPIV( K )
130 IF( KP.NE.K )
131 $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
132 *
133 * Multiply by inv(U(K)), where U(K) is the transformation
134 * stored in column K of A.
135 *
136 CALL DGER( K-1, NRHS, -ONE, AP( KC ), 1, B( K, 1 ), LDB,
137 $ B( 1, 1 ), LDB )
138 *
139 * Multiply by the inverse of the diagonal block.
140 *
141 CALL DSCAL( NRHS, ONE / AP( KC+K-1 ), B( K, 1 ), LDB )
142 K = K - 1
143 ELSE
144 *
145 * 2 x 2 diagonal block
146 *
147 * Interchange rows K-1 and -IPIV(K).
148 *
149 KP = -IPIV( K )
150 IF( KP.NE.K-1 )
151 $ CALL DSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ), LDB )
152 *
153 * Multiply by inv(U(K)), where U(K) is the transformation
154 * stored in columns K-1 and K of A.
155 *
156 CALL DGER( K-2, NRHS, -ONE, AP( KC ), 1, B( K, 1 ), LDB,
157 $ B( 1, 1 ), LDB )
158 CALL DGER( K-2, NRHS, -ONE, AP( KC-( K-1 ) ), 1,
159 $ B( K-1, 1 ), LDB, B( 1, 1 ), LDB )
160 *
161 * Multiply by the inverse of the diagonal block.
162 *
163 AKM1K = AP( KC+K-2 )
164 AKM1 = AP( KC-1 ) / AKM1K
165 AK = AP( KC+K-1 ) / AKM1K
166 DENOM = AKM1*AK - ONE
167 DO 20 J = 1, NRHS
168 BKM1 = B( K-1, J ) / AKM1K
169 BK = B( K, J ) / AKM1K
170 B( K-1, J ) = ( AK*BKM1-BK ) / DENOM
171 B( K, J ) = ( AKM1*BK-BKM1 ) / DENOM
172 20 CONTINUE
173 KC = KC - K + 1
174 K = K - 2
175 END IF
176 *
177 GO TO 10
178 30 CONTINUE
179 *
180 * Next solve U**T*X = B, overwriting B with X.
181 *
182 * K is the main loop index, increasing from 1 to N in steps of
183 * 1 or 2, depending on the size of the diagonal blocks.
184 *
185 K = 1
186 KC = 1
187 40 CONTINUE
188 *
189 * If K > N, exit from loop.
190 *
191 IF( K.GT.N )
192 $ GO TO 50
193 *
194 IF( IPIV( K ).GT.0 ) THEN
195 *
196 * 1 x 1 diagonal block
197 *
198 * Multiply by inv(U**T(K)), where U(K) is the transformation
199 * stored in column K of A.
200 *
201 CALL DGEMV( 'Transpose', K-1, NRHS, -ONE, B, LDB, AP( KC ),
202 $ 1, ONE, B( K, 1 ), LDB )
203 *
204 * Interchange rows K and IPIV(K).
205 *
206 KP = IPIV( K )
207 IF( KP.NE.K )
208 $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
209 KC = KC + K
210 K = K + 1
211 ELSE
212 *
213 * 2 x 2 diagonal block
214 *
215 * Multiply by inv(U**T(K+1)), where U(K+1) is the transformation
216 * stored in columns K and K+1 of A.
217 *
218 CALL DGEMV( 'Transpose', K-1, NRHS, -ONE, B, LDB, AP( KC ),
219 $ 1, ONE, B( K, 1 ), LDB )
220 CALL DGEMV( 'Transpose', K-1, NRHS, -ONE, B, LDB,
221 $ AP( KC+K ), 1, ONE, B( K+1, 1 ), LDB )
222 *
223 * Interchange rows K and -IPIV(K).
224 *
225 KP = -IPIV( K )
226 IF( KP.NE.K )
227 $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
228 KC = KC + 2*K + 1
229 K = K + 2
230 END IF
231 *
232 GO TO 40
233 50 CONTINUE
234 *
235 ELSE
236 *
237 * Solve A*X = B, where A = L*D*L**T.
238 *
239 * First solve L*D*X = B, overwriting B with X.
240 *
241 * K is the main loop index, increasing from 1 to N in steps of
242 * 1 or 2, depending on the size of the diagonal blocks.
243 *
244 K = 1
245 KC = 1
246 60 CONTINUE
247 *
248 * If K > N, exit from loop.
249 *
250 IF( K.GT.N )
251 $ GO TO 80
252 *
253 IF( IPIV( K ).GT.0 ) THEN
254 *
255 * 1 x 1 diagonal block
256 *
257 * Interchange rows K and IPIV(K).
258 *
259 KP = IPIV( K )
260 IF( KP.NE.K )
261 $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
262 *
263 * Multiply by inv(L(K)), where L(K) is the transformation
264 * stored in column K of A.
265 *
266 IF( K.LT.N )
267 $ CALL DGER( N-K, NRHS, -ONE, AP( KC+1 ), 1, B( K, 1 ),
268 $ LDB, B( K+1, 1 ), LDB )
269 *
270 * Multiply by the inverse of the diagonal block.
271 *
272 CALL DSCAL( NRHS, ONE / AP( KC ), B( K, 1 ), LDB )
273 KC = KC + N - K + 1
274 K = K + 1
275 ELSE
276 *
277 * 2 x 2 diagonal block
278 *
279 * Interchange rows K+1 and -IPIV(K).
280 *
281 KP = -IPIV( K )
282 IF( KP.NE.K+1 )
283 $ CALL DSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ), LDB )
284 *
285 * Multiply by inv(L(K)), where L(K) is the transformation
286 * stored in columns K and K+1 of A.
287 *
288 IF( K.LT.N-1 ) THEN
289 CALL DGER( N-K-1, NRHS, -ONE, AP( KC+2 ), 1, B( K, 1 ),
290 $ LDB, B( K+2, 1 ), LDB )
291 CALL DGER( N-K-1, NRHS, -ONE, AP( KC+N-K+2 ), 1,
292 $ B( K+1, 1 ), LDB, B( K+2, 1 ), LDB )
293 END IF
294 *
295 * Multiply by the inverse of the diagonal block.
296 *
297 AKM1K = AP( KC+1 )
298 AKM1 = AP( KC ) / AKM1K
299 AK = AP( KC+N-K+1 ) / AKM1K
300 DENOM = AKM1*AK - ONE
301 DO 70 J = 1, NRHS
302 BKM1 = B( K, J ) / AKM1K
303 BK = B( K+1, J ) / AKM1K
304 B( K, J ) = ( AK*BKM1-BK ) / DENOM
305 B( K+1, J ) = ( AKM1*BK-BKM1 ) / DENOM
306 70 CONTINUE
307 KC = KC + 2*( N-K ) + 1
308 K = K + 2
309 END IF
310 *
311 GO TO 60
312 80 CONTINUE
313 *
314 * Next solve L**T*X = B, overwriting B with X.
315 *
316 * K is the main loop index, decreasing from N to 1 in steps of
317 * 1 or 2, depending on the size of the diagonal blocks.
318 *
319 K = N
320 KC = N*( N+1 ) / 2 + 1
321 90 CONTINUE
322 *
323 * If K < 1, exit from loop.
324 *
325 IF( K.LT.1 )
326 $ GO TO 100
327 *
328 KC = KC - ( N-K+1 )
329 IF( IPIV( K ).GT.0 ) THEN
330 *
331 * 1 x 1 diagonal block
332 *
333 * Multiply by inv(L**T(K)), where L(K) is the transformation
334 * stored in column K of A.
335 *
336 IF( K.LT.N )
337 $ CALL DGEMV( 'Transpose', N-K, NRHS, -ONE, B( K+1, 1 ),
338 $ LDB, AP( KC+1 ), 1, ONE, B( K, 1 ), LDB )
339 *
340 * Interchange rows K and IPIV(K).
341 *
342 KP = IPIV( K )
343 IF( KP.NE.K )
344 $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
345 K = K - 1
346 ELSE
347 *
348 * 2 x 2 diagonal block
349 *
350 * Multiply by inv(L**T(K-1)), where L(K-1) is the transformation
351 * stored in columns K-1 and K of A.
352 *
353 IF( K.LT.N ) THEN
354 CALL DGEMV( 'Transpose', N-K, NRHS, -ONE, B( K+1, 1 ),
355 $ LDB, AP( KC+1 ), 1, ONE, B( K, 1 ), LDB )
356 CALL DGEMV( 'Transpose', N-K, NRHS, -ONE, B( K+1, 1 ),
357 $ LDB, AP( KC-( N-K ) ), 1, ONE, B( K-1, 1 ),
358 $ LDB )
359 END IF
360 *
361 * Interchange rows K and -IPIV(K).
362 *
363 KP = -IPIV( K )
364 IF( KP.NE.K )
365 $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
366 KC = KC - ( N-K+2 )
367 K = K - 2
368 END IF
369 *
370 GO TO 90
371 100 CONTINUE
372 END IF
373 *
374 RETURN
375 *
376 * End of DSPTRS
377 *
378 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, LDB, N, NRHS
11 * ..
12 * .. Array Arguments ..
13 INTEGER IPIV( * )
14 DOUBLE PRECISION AP( * ), B( LDB, * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * DSPTRS solves a system of linear equations A*X = B with a real
21 * symmetric matrix A stored in packed format using the factorization
22 * A = U*D*U**T or 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 * NRHS (input) INTEGER
37 * The number of right hand sides, i.e., the number of columns
38 * of the matrix B. NRHS >= 0.
39 *
40 * AP (input) DOUBLE PRECISION array, dimension (N*(N+1)/2)
41 * The block diagonal matrix D and the multipliers used to
42 * obtain the factor U or L as computed by DSPTRF, stored as a
43 * packed triangular matrix.
44 *
45 * IPIV (input) INTEGER array, dimension (N)
46 * Details of the interchanges and the block structure of D
47 * as determined by DSPTRF.
48 *
49 * B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
50 * On entry, the right hand side matrix B.
51 * On exit, the solution matrix X.
52 *
53 * LDB (input) INTEGER
54 * The leading dimension of the array B. LDB >= max(1,N).
55 *
56 * INFO (output) INTEGER
57 * = 0: successful exit
58 * < 0: if INFO = -i, the i-th argument had an illegal value
59 *
60 * =====================================================================
61 *
62 * .. Parameters ..
63 DOUBLE PRECISION ONE
64 PARAMETER ( ONE = 1.0D+0 )
65 * ..
66 * .. Local Scalars ..
67 LOGICAL UPPER
68 INTEGER J, K, KC, KP
69 DOUBLE PRECISION AK, AKM1, AKM1K, BK, BKM1, DENOM
70 * ..
71 * .. External Functions ..
72 LOGICAL LSAME
73 EXTERNAL LSAME
74 * ..
75 * .. External Subroutines ..
76 EXTERNAL DGEMV, DGER, DSCAL, DSWAP, XERBLA
77 * ..
78 * .. Intrinsic Functions ..
79 INTRINSIC MAX
80 * ..
81 * .. Executable Statements ..
82 *
83 INFO = 0
84 UPPER = LSAME( UPLO, 'U' )
85 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
86 INFO = -1
87 ELSE IF( N.LT.0 ) THEN
88 INFO = -2
89 ELSE IF( NRHS.LT.0 ) THEN
90 INFO = -3
91 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
92 INFO = -7
93 END IF
94 IF( INFO.NE.0 ) THEN
95 CALL XERBLA( 'DSPTRS', -INFO )
96 RETURN
97 END IF
98 *
99 * Quick return if possible
100 *
101 IF( N.EQ.0 .OR. NRHS.EQ.0 )
102 $ RETURN
103 *
104 IF( UPPER ) THEN
105 *
106 * Solve A*X = B, where A = U*D*U**T.
107 *
108 * First solve U*D*X = B, overwriting B with X.
109 *
110 * K is the main loop index, decreasing from N to 1 in steps of
111 * 1 or 2, depending on the size of the diagonal blocks.
112 *
113 K = N
114 KC = N*( N+1 ) / 2 + 1
115 10 CONTINUE
116 *
117 * If K < 1, exit from loop.
118 *
119 IF( K.LT.1 )
120 $ GO TO 30
121 *
122 KC = KC - K
123 IF( IPIV( K ).GT.0 ) THEN
124 *
125 * 1 x 1 diagonal block
126 *
127 * Interchange rows K and IPIV(K).
128 *
129 KP = IPIV( K )
130 IF( KP.NE.K )
131 $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
132 *
133 * Multiply by inv(U(K)), where U(K) is the transformation
134 * stored in column K of A.
135 *
136 CALL DGER( K-1, NRHS, -ONE, AP( KC ), 1, B( K, 1 ), LDB,
137 $ B( 1, 1 ), LDB )
138 *
139 * Multiply by the inverse of the diagonal block.
140 *
141 CALL DSCAL( NRHS, ONE / AP( KC+K-1 ), B( K, 1 ), LDB )
142 K = K - 1
143 ELSE
144 *
145 * 2 x 2 diagonal block
146 *
147 * Interchange rows K-1 and -IPIV(K).
148 *
149 KP = -IPIV( K )
150 IF( KP.NE.K-1 )
151 $ CALL DSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ), LDB )
152 *
153 * Multiply by inv(U(K)), where U(K) is the transformation
154 * stored in columns K-1 and K of A.
155 *
156 CALL DGER( K-2, NRHS, -ONE, AP( KC ), 1, B( K, 1 ), LDB,
157 $ B( 1, 1 ), LDB )
158 CALL DGER( K-2, NRHS, -ONE, AP( KC-( K-1 ) ), 1,
159 $ B( K-1, 1 ), LDB, B( 1, 1 ), LDB )
160 *
161 * Multiply by the inverse of the diagonal block.
162 *
163 AKM1K = AP( KC+K-2 )
164 AKM1 = AP( KC-1 ) / AKM1K
165 AK = AP( KC+K-1 ) / AKM1K
166 DENOM = AKM1*AK - ONE
167 DO 20 J = 1, NRHS
168 BKM1 = B( K-1, J ) / AKM1K
169 BK = B( K, J ) / AKM1K
170 B( K-1, J ) = ( AK*BKM1-BK ) / DENOM
171 B( K, J ) = ( AKM1*BK-BKM1 ) / DENOM
172 20 CONTINUE
173 KC = KC - K + 1
174 K = K - 2
175 END IF
176 *
177 GO TO 10
178 30 CONTINUE
179 *
180 * Next solve U**T*X = B, overwriting B with X.
181 *
182 * K is the main loop index, increasing from 1 to N in steps of
183 * 1 or 2, depending on the size of the diagonal blocks.
184 *
185 K = 1
186 KC = 1
187 40 CONTINUE
188 *
189 * If K > N, exit from loop.
190 *
191 IF( K.GT.N )
192 $ GO TO 50
193 *
194 IF( IPIV( K ).GT.0 ) THEN
195 *
196 * 1 x 1 diagonal block
197 *
198 * Multiply by inv(U**T(K)), where U(K) is the transformation
199 * stored in column K of A.
200 *
201 CALL DGEMV( 'Transpose', K-1, NRHS, -ONE, B, LDB, AP( KC ),
202 $ 1, ONE, B( K, 1 ), LDB )
203 *
204 * Interchange rows K and IPIV(K).
205 *
206 KP = IPIV( K )
207 IF( KP.NE.K )
208 $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
209 KC = KC + K
210 K = K + 1
211 ELSE
212 *
213 * 2 x 2 diagonal block
214 *
215 * Multiply by inv(U**T(K+1)), where U(K+1) is the transformation
216 * stored in columns K and K+1 of A.
217 *
218 CALL DGEMV( 'Transpose', K-1, NRHS, -ONE, B, LDB, AP( KC ),
219 $ 1, ONE, B( K, 1 ), LDB )
220 CALL DGEMV( 'Transpose', K-1, NRHS, -ONE, B, LDB,
221 $ AP( KC+K ), 1, ONE, B( K+1, 1 ), LDB )
222 *
223 * Interchange rows K and -IPIV(K).
224 *
225 KP = -IPIV( K )
226 IF( KP.NE.K )
227 $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
228 KC = KC + 2*K + 1
229 K = K + 2
230 END IF
231 *
232 GO TO 40
233 50 CONTINUE
234 *
235 ELSE
236 *
237 * Solve A*X = B, where A = L*D*L**T.
238 *
239 * First solve L*D*X = B, overwriting B with X.
240 *
241 * K is the main loop index, increasing from 1 to N in steps of
242 * 1 or 2, depending on the size of the diagonal blocks.
243 *
244 K = 1
245 KC = 1
246 60 CONTINUE
247 *
248 * If K > N, exit from loop.
249 *
250 IF( K.GT.N )
251 $ GO TO 80
252 *
253 IF( IPIV( K ).GT.0 ) THEN
254 *
255 * 1 x 1 diagonal block
256 *
257 * Interchange rows K and IPIV(K).
258 *
259 KP = IPIV( K )
260 IF( KP.NE.K )
261 $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
262 *
263 * Multiply by inv(L(K)), where L(K) is the transformation
264 * stored in column K of A.
265 *
266 IF( K.LT.N )
267 $ CALL DGER( N-K, NRHS, -ONE, AP( KC+1 ), 1, B( K, 1 ),
268 $ LDB, B( K+1, 1 ), LDB )
269 *
270 * Multiply by the inverse of the diagonal block.
271 *
272 CALL DSCAL( NRHS, ONE / AP( KC ), B( K, 1 ), LDB )
273 KC = KC + N - K + 1
274 K = K + 1
275 ELSE
276 *
277 * 2 x 2 diagonal block
278 *
279 * Interchange rows K+1 and -IPIV(K).
280 *
281 KP = -IPIV( K )
282 IF( KP.NE.K+1 )
283 $ CALL DSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ), LDB )
284 *
285 * Multiply by inv(L(K)), where L(K) is the transformation
286 * stored in columns K and K+1 of A.
287 *
288 IF( K.LT.N-1 ) THEN
289 CALL DGER( N-K-1, NRHS, -ONE, AP( KC+2 ), 1, B( K, 1 ),
290 $ LDB, B( K+2, 1 ), LDB )
291 CALL DGER( N-K-1, NRHS, -ONE, AP( KC+N-K+2 ), 1,
292 $ B( K+1, 1 ), LDB, B( K+2, 1 ), LDB )
293 END IF
294 *
295 * Multiply by the inverse of the diagonal block.
296 *
297 AKM1K = AP( KC+1 )
298 AKM1 = AP( KC ) / AKM1K
299 AK = AP( KC+N-K+1 ) / AKM1K
300 DENOM = AKM1*AK - ONE
301 DO 70 J = 1, NRHS
302 BKM1 = B( K, J ) / AKM1K
303 BK = B( K+1, J ) / AKM1K
304 B( K, J ) = ( AK*BKM1-BK ) / DENOM
305 B( K+1, J ) = ( AKM1*BK-BKM1 ) / DENOM
306 70 CONTINUE
307 KC = KC + 2*( N-K ) + 1
308 K = K + 2
309 END IF
310 *
311 GO TO 60
312 80 CONTINUE
313 *
314 * Next solve L**T*X = B, overwriting B with X.
315 *
316 * K is the main loop index, decreasing from N to 1 in steps of
317 * 1 or 2, depending on the size of the diagonal blocks.
318 *
319 K = N
320 KC = N*( N+1 ) / 2 + 1
321 90 CONTINUE
322 *
323 * If K < 1, exit from loop.
324 *
325 IF( K.LT.1 )
326 $ GO TO 100
327 *
328 KC = KC - ( N-K+1 )
329 IF( IPIV( K ).GT.0 ) THEN
330 *
331 * 1 x 1 diagonal block
332 *
333 * Multiply by inv(L**T(K)), where L(K) is the transformation
334 * stored in column K of A.
335 *
336 IF( K.LT.N )
337 $ CALL DGEMV( 'Transpose', N-K, NRHS, -ONE, B( K+1, 1 ),
338 $ LDB, AP( KC+1 ), 1, ONE, B( K, 1 ), LDB )
339 *
340 * Interchange rows K and IPIV(K).
341 *
342 KP = IPIV( K )
343 IF( KP.NE.K )
344 $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
345 K = K - 1
346 ELSE
347 *
348 * 2 x 2 diagonal block
349 *
350 * Multiply by inv(L**T(K-1)), where L(K-1) is the transformation
351 * stored in columns K-1 and K of A.
352 *
353 IF( K.LT.N ) THEN
354 CALL DGEMV( 'Transpose', N-K, NRHS, -ONE, B( K+1, 1 ),
355 $ LDB, AP( KC+1 ), 1, ONE, B( K, 1 ), LDB )
356 CALL DGEMV( 'Transpose', N-K, NRHS, -ONE, B( K+1, 1 ),
357 $ LDB, AP( KC-( N-K ) ), 1, ONE, B( K-1, 1 ),
358 $ LDB )
359 END IF
360 *
361 * Interchange rows K and -IPIV(K).
362 *
363 KP = -IPIV( K )
364 IF( KP.NE.K )
365 $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
366 KC = KC - ( N-K+2 )
367 K = K - 2
368 END IF
369 *
370 GO TO 90
371 100 CONTINUE
372 END IF
373 *
374 RETURN
375 *
376 * End of DSPTRS
377 *
378 END