1 SUBROUTINE DLAVSP( UPLO, TRANS, DIAG, N, NRHS, A, IPIV, B, LDB,
2 $ INFO )
3 *
4 * -- LAPACK auxiliary routine (version 3.1) --
5 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
6 * November 2006
7 *
8 * .. Scalar Arguments ..
9 CHARACTER DIAG, TRANS, UPLO
10 INTEGER INFO, LDB, N, NRHS
11 * ..
12 * .. Array Arguments ..
13 INTEGER IPIV( * )
14 DOUBLE PRECISION A( * ), B( LDB, * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * DLAVSP performs one of the matrix-vector operations
21 * x := A*x or x := A'*x,
22 * where x is an N element vector and A is one of the factors
23 * from the block U*D*U' or L*D*L' factorization computed by DSPTRF.
24 *
25 * If TRANS = 'N', multiplies by U or U * D (or L or L * D)
26 * If TRANS = 'T', multiplies by U' or D * U' (or L' or D * L' )
27 * If TRANS = 'C', multiplies by U' or D * U' (or L' or D * L' )
28 *
29 * Arguments
30 * ==========
31 *
32 * UPLO (input) CHARACTER*1
33 * Specifies whether the factor stored in A is upper or lower
34 * triangular.
35 * = 'U': Upper triangular
36 * = 'L': Lower triangular
37 *
38 * TRANS (input) CHARACTER*1
39 * Specifies the operation to be performed:
40 * = 'N': x := A*x
41 * = 'T': x := A'*x
42 * = 'C': x := A'*x
43 *
44 * DIAG (input) CHARACTER*1
45 * Specifies whether or not the diagonal blocks are unit
46 * matrices. If the diagonal blocks are assumed to be unit,
47 * then A = U or A = L, otherwise A = U*D or A = L*D.
48 * = 'U': Diagonal blocks are assumed to be unit matrices.
49 * = 'N': Diagonal blocks are assumed to be non-unit matrices.
50 *
51 * N (input) INTEGER
52 * The number of rows and columns of the matrix A. N >= 0.
53 *
54 * NRHS (input) INTEGER
55 * The number of right hand sides, i.e., the number of vectors
56 * x to be multiplied by A. NRHS >= 0.
57 *
58 * A (input) DOUBLE PRECISION array, dimension (N*(N+1)/2)
59 * The block diagonal matrix D and the multipliers used to
60 * obtain the factor U or L, stored as a packed triangular
61 * matrix as computed by DSPTRF.
62 *
63 * IPIV (input) INTEGER array, dimension (N)
64 * The pivot indices from DSPTRF.
65 *
66 * B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
67 * On entry, B contains NRHS vectors of length N.
68 * On exit, B is overwritten with the product A * B.
69 *
70 * LDB (input) INTEGER
71 * The leading dimension of the array B. LDB >= max(1,N).
72 *
73 * INFO (output) INTEGER
74 * = 0: successful exit
75 * < 0: if INFO = -k, the k-th argument had an illegal value
76 *
77 * =====================================================================
78 *
79 * .. Parameters ..
80 DOUBLE PRECISION ONE
81 PARAMETER ( ONE = 1.0D+0 )
82 * ..
83 * .. Local Scalars ..
84 LOGICAL NOUNIT
85 INTEGER J, K, KC, KCNEXT, KP
86 DOUBLE PRECISION D11, D12, D21, D22, T1, T2
87 * ..
88 * .. External Functions ..
89 LOGICAL LSAME
90 EXTERNAL LSAME
91 * ..
92 * .. External Subroutines ..
93 EXTERNAL DGEMV, DGER, DSCAL, DSWAP, XERBLA
94 * ..
95 * .. Intrinsic Functions ..
96 INTRINSIC ABS, MAX
97 * ..
98 * .. Executable Statements ..
99 *
100 * Test the input parameters.
101 *
102 INFO = 0
103 IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
104 INFO = -1
105 ELSE IF( .NOT.LSAME( TRANS, 'N' ) .AND. .NOT.
106 $ LSAME( TRANS, 'T' ) .AND. .NOT.LSAME( TRANS, 'C' ) ) THEN
107 INFO = -2
108 ELSE IF( .NOT.LSAME( DIAG, 'U' ) .AND. .NOT.LSAME( DIAG, 'N' ) )
109 $ THEN
110 INFO = -3
111 ELSE IF( N.LT.0 ) THEN
112 INFO = -4
113 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
114 INFO = -8
115 END IF
116 IF( INFO.NE.0 ) THEN
117 CALL XERBLA( 'DLAVSP ', -INFO )
118 RETURN
119 END IF
120 *
121 * Quick return if possible.
122 *
123 IF( N.EQ.0 )
124 $ RETURN
125 *
126 NOUNIT = LSAME( DIAG, 'N' )
127 *------------------------------------------
128 *
129 * Compute B := A * B (No transpose)
130 *
131 *------------------------------------------
132 IF( LSAME( TRANS, 'N' ) ) THEN
133 *
134 * Compute B := U*B
135 * where U = P(m)*inv(U(m))* ... *P(1)*inv(U(1))
136 *
137 IF( LSAME( UPLO, 'U' ) ) THEN
138 *
139 * Loop forward applying the transformations.
140 *
141 K = 1
142 KC = 1
143 10 CONTINUE
144 IF( K.GT.N )
145 $ GO TO 30
146 *
147 * 1 x 1 pivot block
148 *
149 IF( IPIV( K ).GT.0 ) THEN
150 *
151 * Multiply by the diagonal element if forming U * D.
152 *
153 IF( NOUNIT )
154 $ CALL DSCAL( NRHS, A( KC+K-1 ), B( K, 1 ), LDB )
155 *
156 * Multiply by P(K) * inv(U(K)) if K > 1.
157 *
158 IF( K.GT.1 ) THEN
159 *
160 * Apply the transformation.
161 *
162 CALL DGER( K-1, NRHS, ONE, A( KC ), 1, B( K, 1 ), LDB,
163 $ B( 1, 1 ), LDB )
164 *
165 * Interchange if P(K) != I.
166 *
167 KP = IPIV( K )
168 IF( KP.NE.K )
169 $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
170 END IF
171 KC = KC + K
172 K = K + 1
173 ELSE
174 *
175 * 2 x 2 pivot block
176 *
177 KCNEXT = KC + K
178 *
179 * Multiply by the diagonal block if forming U * D.
180 *
181 IF( NOUNIT ) THEN
182 D11 = A( KCNEXT-1 )
183 D22 = A( KCNEXT+K )
184 D12 = A( KCNEXT+K-1 )
185 D21 = D12
186 DO 20 J = 1, NRHS
187 T1 = B( K, J )
188 T2 = B( K+1, J )
189 B( K, J ) = D11*T1 + D12*T2
190 B( K+1, J ) = D21*T1 + D22*T2
191 20 CONTINUE
192 END IF
193 *
194 * Multiply by P(K) * inv(U(K)) if K > 1.
195 *
196 IF( K.GT.1 ) THEN
197 *
198 * Apply the transformations.
199 *
200 CALL DGER( K-1, NRHS, ONE, A( KC ), 1, B( K, 1 ), LDB,
201 $ B( 1, 1 ), LDB )
202 CALL DGER( K-1, NRHS, ONE, A( KCNEXT ), 1,
203 $ B( K+1, 1 ), LDB, B( 1, 1 ), LDB )
204 *
205 * Interchange if P(K) != I.
206 *
207 KP = ABS( IPIV( K ) )
208 IF( KP.NE.K )
209 $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
210 END IF
211 KC = KCNEXT + K + 1
212 K = K + 2
213 END IF
214 GO TO 10
215 30 CONTINUE
216 *
217 * Compute B := L*B
218 * where L = P(1)*inv(L(1))* ... *P(m)*inv(L(m)) .
219 *
220 ELSE
221 *
222 * Loop backward applying the transformations to B.
223 *
224 K = N
225 KC = N*( N+1 ) / 2 + 1
226 40 CONTINUE
227 IF( K.LT.1 )
228 $ GO TO 60
229 KC = KC - ( N-K+1 )
230 *
231 * Test the pivot index. If greater than zero, a 1 x 1
232 * pivot was used, otherwise a 2 x 2 pivot was used.
233 *
234 IF( IPIV( K ).GT.0 ) THEN
235 *
236 * 1 x 1 pivot block:
237 *
238 * Multiply by the diagonal element if forming L * D.
239 *
240 IF( NOUNIT )
241 $ CALL DSCAL( NRHS, A( KC ), B( K, 1 ), LDB )
242 *
243 * Multiply by P(K) * inv(L(K)) if K < N.
244 *
245 IF( K.NE.N ) THEN
246 KP = IPIV( K )
247 *
248 * Apply the transformation.
249 *
250 CALL DGER( N-K, NRHS, ONE, A( KC+1 ), 1, B( K, 1 ),
251 $ LDB, B( K+1, 1 ), LDB )
252 *
253 * Interchange if a permutation was applied at the
254 * K-th step of the factorization.
255 *
256 IF( KP.NE.K )
257 $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
258 END IF
259 K = K - 1
260 *
261 ELSE
262 *
263 * 2 x 2 pivot block:
264 *
265 KCNEXT = KC - ( N-K+2 )
266 *
267 * Multiply by the diagonal block if forming L * D.
268 *
269 IF( NOUNIT ) THEN
270 D11 = A( KCNEXT )
271 D22 = A( KC )
272 D21 = A( KCNEXT+1 )
273 D12 = D21
274 DO 50 J = 1, NRHS
275 T1 = B( K-1, J )
276 T2 = B( K, J )
277 B( K-1, J ) = D11*T1 + D12*T2
278 B( K, J ) = D21*T1 + D22*T2
279 50 CONTINUE
280 END IF
281 *
282 * Multiply by P(K) * inv(L(K)) if K < N.
283 *
284 IF( K.NE.N ) THEN
285 *
286 * Apply the transformation.
287 *
288 CALL DGER( N-K, NRHS, ONE, A( KC+1 ), 1, B( K, 1 ),
289 $ LDB, B( K+1, 1 ), LDB )
290 CALL DGER( N-K, NRHS, ONE, A( KCNEXT+2 ), 1,
291 $ B( K-1, 1 ), LDB, B( K+1, 1 ), LDB )
292 *
293 * Interchange if a permutation was applied at the
294 * K-th step of the factorization.
295 *
296 KP = ABS( IPIV( K ) )
297 IF( KP.NE.K )
298 $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
299 END IF
300 KC = KCNEXT
301 K = K - 2
302 END IF
303 GO TO 40
304 60 CONTINUE
305 END IF
306 *----------------------------------------
307 *
308 * Compute B := A' * B (transpose)
309 *
310 *----------------------------------------
311 ELSE
312 *
313 * Form B := U'*B
314 * where U = P(m)*inv(U(m))* ... *P(1)*inv(U(1))
315 * and U' = inv(U'(1))*P(1)* ... *inv(U'(m))*P(m)
316 *
317 IF( LSAME( UPLO, 'U' ) ) THEN
318 *
319 * Loop backward applying the transformations.
320 *
321 K = N
322 KC = N*( N+1 ) / 2 + 1
323 70 CONTINUE
324 IF( K.LT.1 )
325 $ GO TO 90
326 KC = KC - K
327 *
328 * 1 x 1 pivot block.
329 *
330 IF( IPIV( K ).GT.0 ) THEN
331 IF( K.GT.1 ) THEN
332 *
333 * Interchange if P(K) != I.
334 *
335 KP = IPIV( K )
336 IF( KP.NE.K )
337 $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
338 *
339 * Apply the transformation
340 *
341 CALL DGEMV( 'Transpose', K-1, NRHS, ONE, B, LDB,
342 $ A( KC ), 1, ONE, B( K, 1 ), LDB )
343 END IF
344 IF( NOUNIT )
345 $ CALL DSCAL( NRHS, A( KC+K-1 ), B( K, 1 ), LDB )
346 K = K - 1
347 *
348 * 2 x 2 pivot block.
349 *
350 ELSE
351 KCNEXT = KC - ( K-1 )
352 IF( K.GT.2 ) THEN
353 *
354 * Interchange if P(K) != I.
355 *
356 KP = ABS( IPIV( K ) )
357 IF( KP.NE.K-1 )
358 $ CALL DSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ),
359 $ LDB )
360 *
361 * Apply the transformations
362 *
363 CALL DGEMV( 'Transpose', K-2, NRHS, ONE, B, LDB,
364 $ A( KC ), 1, ONE, B( K, 1 ), LDB )
365 CALL DGEMV( 'Transpose', K-2, NRHS, ONE, B, LDB,
366 $ A( KCNEXT ), 1, ONE, B( K-1, 1 ), LDB )
367 END IF
368 *
369 * Multiply by the diagonal block if non-unit.
370 *
371 IF( NOUNIT ) THEN
372 D11 = A( KC-1 )
373 D22 = A( KC+K-1 )
374 D12 = A( KC+K-2 )
375 D21 = D12
376 DO 80 J = 1, NRHS
377 T1 = B( K-1, J )
378 T2 = B( K, J )
379 B( K-1, J ) = D11*T1 + D12*T2
380 B( K, J ) = D21*T1 + D22*T2
381 80 CONTINUE
382 END IF
383 KC = KCNEXT
384 K = K - 2
385 END IF
386 GO TO 70
387 90 CONTINUE
388 *
389 * Form B := L'*B
390 * where L = P(1)*inv(L(1))* ... *P(m)*inv(L(m))
391 * and L' = inv(L(m))*P(m)* ... *inv(L(1))*P(1)
392 *
393 ELSE
394 *
395 * Loop forward applying the L-transformations.
396 *
397 K = 1
398 KC = 1
399 100 CONTINUE
400 IF( K.GT.N )
401 $ GO TO 120
402 *
403 * 1 x 1 pivot block
404 *
405 IF( IPIV( K ).GT.0 ) THEN
406 IF( K.LT.N ) THEN
407 *
408 * Interchange if P(K) != I.
409 *
410 KP = IPIV( K )
411 IF( KP.NE.K )
412 $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
413 *
414 * Apply the transformation
415 *
416 CALL DGEMV( 'Transpose', N-K, NRHS, ONE, B( K+1, 1 ),
417 $ LDB, A( KC+1 ), 1, ONE, B( K, 1 ), LDB )
418 END IF
419 IF( NOUNIT )
420 $ CALL DSCAL( NRHS, A( KC ), B( K, 1 ), LDB )
421 KC = KC + N - K + 1
422 K = K + 1
423 *
424 * 2 x 2 pivot block.
425 *
426 ELSE
427 KCNEXT = KC + N - K + 1
428 IF( K.LT.N-1 ) THEN
429 *
430 * Interchange if P(K) != I.
431 *
432 KP = ABS( IPIV( K ) )
433 IF( KP.NE.K+1 )
434 $ CALL DSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ),
435 $ LDB )
436 *
437 * Apply the transformation
438 *
439 CALL DGEMV( 'Transpose', N-K-1, NRHS, ONE,
440 $ B( K+2, 1 ), LDB, A( KCNEXT+1 ), 1, ONE,
441 $ B( K+1, 1 ), LDB )
442 CALL DGEMV( 'Transpose', N-K-1, NRHS, ONE,
443 $ B( K+2, 1 ), LDB, A( KC+2 ), 1, ONE,
444 $ B( K, 1 ), LDB )
445 END IF
446 *
447 * Multiply by the diagonal block if non-unit.
448 *
449 IF( NOUNIT ) THEN
450 D11 = A( KC )
451 D22 = A( KCNEXT )
452 D21 = A( KC+1 )
453 D12 = D21
454 DO 110 J = 1, NRHS
455 T1 = B( K, J )
456 T2 = B( K+1, J )
457 B( K, J ) = D11*T1 + D12*T2
458 B( K+1, J ) = D21*T1 + D22*T2
459 110 CONTINUE
460 END IF
461 KC = KCNEXT + ( N-K )
462 K = K + 2
463 END IF
464 GO TO 100
465 120 CONTINUE
466 END IF
467 *
468 END IF
469 RETURN
470 *
471 * End of DLAVSP
472 *
473 END
2 $ INFO )
3 *
4 * -- LAPACK auxiliary routine (version 3.1) --
5 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
6 * November 2006
7 *
8 * .. Scalar Arguments ..
9 CHARACTER DIAG, TRANS, UPLO
10 INTEGER INFO, LDB, N, NRHS
11 * ..
12 * .. Array Arguments ..
13 INTEGER IPIV( * )
14 DOUBLE PRECISION A( * ), B( LDB, * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * DLAVSP performs one of the matrix-vector operations
21 * x := A*x or x := A'*x,
22 * where x is an N element vector and A is one of the factors
23 * from the block U*D*U' or L*D*L' factorization computed by DSPTRF.
24 *
25 * If TRANS = 'N', multiplies by U or U * D (or L or L * D)
26 * If TRANS = 'T', multiplies by U' or D * U' (or L' or D * L' )
27 * If TRANS = 'C', multiplies by U' or D * U' (or L' or D * L' )
28 *
29 * Arguments
30 * ==========
31 *
32 * UPLO (input) CHARACTER*1
33 * Specifies whether the factor stored in A is upper or lower
34 * triangular.
35 * = 'U': Upper triangular
36 * = 'L': Lower triangular
37 *
38 * TRANS (input) CHARACTER*1
39 * Specifies the operation to be performed:
40 * = 'N': x := A*x
41 * = 'T': x := A'*x
42 * = 'C': x := A'*x
43 *
44 * DIAG (input) CHARACTER*1
45 * Specifies whether or not the diagonal blocks are unit
46 * matrices. If the diagonal blocks are assumed to be unit,
47 * then A = U or A = L, otherwise A = U*D or A = L*D.
48 * = 'U': Diagonal blocks are assumed to be unit matrices.
49 * = 'N': Diagonal blocks are assumed to be non-unit matrices.
50 *
51 * N (input) INTEGER
52 * The number of rows and columns of the matrix A. N >= 0.
53 *
54 * NRHS (input) INTEGER
55 * The number of right hand sides, i.e., the number of vectors
56 * x to be multiplied by A. NRHS >= 0.
57 *
58 * A (input) DOUBLE PRECISION array, dimension (N*(N+1)/2)
59 * The block diagonal matrix D and the multipliers used to
60 * obtain the factor U or L, stored as a packed triangular
61 * matrix as computed by DSPTRF.
62 *
63 * IPIV (input) INTEGER array, dimension (N)
64 * The pivot indices from DSPTRF.
65 *
66 * B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
67 * On entry, B contains NRHS vectors of length N.
68 * On exit, B is overwritten with the product A * B.
69 *
70 * LDB (input) INTEGER
71 * The leading dimension of the array B. LDB >= max(1,N).
72 *
73 * INFO (output) INTEGER
74 * = 0: successful exit
75 * < 0: if INFO = -k, the k-th argument had an illegal value
76 *
77 * =====================================================================
78 *
79 * .. Parameters ..
80 DOUBLE PRECISION ONE
81 PARAMETER ( ONE = 1.0D+0 )
82 * ..
83 * .. Local Scalars ..
84 LOGICAL NOUNIT
85 INTEGER J, K, KC, KCNEXT, KP
86 DOUBLE PRECISION D11, D12, D21, D22, T1, T2
87 * ..
88 * .. External Functions ..
89 LOGICAL LSAME
90 EXTERNAL LSAME
91 * ..
92 * .. External Subroutines ..
93 EXTERNAL DGEMV, DGER, DSCAL, DSWAP, XERBLA
94 * ..
95 * .. Intrinsic Functions ..
96 INTRINSIC ABS, MAX
97 * ..
98 * .. Executable Statements ..
99 *
100 * Test the input parameters.
101 *
102 INFO = 0
103 IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
104 INFO = -1
105 ELSE IF( .NOT.LSAME( TRANS, 'N' ) .AND. .NOT.
106 $ LSAME( TRANS, 'T' ) .AND. .NOT.LSAME( TRANS, 'C' ) ) THEN
107 INFO = -2
108 ELSE IF( .NOT.LSAME( DIAG, 'U' ) .AND. .NOT.LSAME( DIAG, 'N' ) )
109 $ THEN
110 INFO = -3
111 ELSE IF( N.LT.0 ) THEN
112 INFO = -4
113 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
114 INFO = -8
115 END IF
116 IF( INFO.NE.0 ) THEN
117 CALL XERBLA( 'DLAVSP ', -INFO )
118 RETURN
119 END IF
120 *
121 * Quick return if possible.
122 *
123 IF( N.EQ.0 )
124 $ RETURN
125 *
126 NOUNIT = LSAME( DIAG, 'N' )
127 *------------------------------------------
128 *
129 * Compute B := A * B (No transpose)
130 *
131 *------------------------------------------
132 IF( LSAME( TRANS, 'N' ) ) THEN
133 *
134 * Compute B := U*B
135 * where U = P(m)*inv(U(m))* ... *P(1)*inv(U(1))
136 *
137 IF( LSAME( UPLO, 'U' ) ) THEN
138 *
139 * Loop forward applying the transformations.
140 *
141 K = 1
142 KC = 1
143 10 CONTINUE
144 IF( K.GT.N )
145 $ GO TO 30
146 *
147 * 1 x 1 pivot block
148 *
149 IF( IPIV( K ).GT.0 ) THEN
150 *
151 * Multiply by the diagonal element if forming U * D.
152 *
153 IF( NOUNIT )
154 $ CALL DSCAL( NRHS, A( KC+K-1 ), B( K, 1 ), LDB )
155 *
156 * Multiply by P(K) * inv(U(K)) if K > 1.
157 *
158 IF( K.GT.1 ) THEN
159 *
160 * Apply the transformation.
161 *
162 CALL DGER( K-1, NRHS, ONE, A( KC ), 1, B( K, 1 ), LDB,
163 $ B( 1, 1 ), LDB )
164 *
165 * Interchange if P(K) != I.
166 *
167 KP = IPIV( K )
168 IF( KP.NE.K )
169 $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
170 END IF
171 KC = KC + K
172 K = K + 1
173 ELSE
174 *
175 * 2 x 2 pivot block
176 *
177 KCNEXT = KC + K
178 *
179 * Multiply by the diagonal block if forming U * D.
180 *
181 IF( NOUNIT ) THEN
182 D11 = A( KCNEXT-1 )
183 D22 = A( KCNEXT+K )
184 D12 = A( KCNEXT+K-1 )
185 D21 = D12
186 DO 20 J = 1, NRHS
187 T1 = B( K, J )
188 T2 = B( K+1, J )
189 B( K, J ) = D11*T1 + D12*T2
190 B( K+1, J ) = D21*T1 + D22*T2
191 20 CONTINUE
192 END IF
193 *
194 * Multiply by P(K) * inv(U(K)) if K > 1.
195 *
196 IF( K.GT.1 ) THEN
197 *
198 * Apply the transformations.
199 *
200 CALL DGER( K-1, NRHS, ONE, A( KC ), 1, B( K, 1 ), LDB,
201 $ B( 1, 1 ), LDB )
202 CALL DGER( K-1, NRHS, ONE, A( KCNEXT ), 1,
203 $ B( K+1, 1 ), LDB, B( 1, 1 ), LDB )
204 *
205 * Interchange if P(K) != I.
206 *
207 KP = ABS( IPIV( K ) )
208 IF( KP.NE.K )
209 $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
210 END IF
211 KC = KCNEXT + K + 1
212 K = K + 2
213 END IF
214 GO TO 10
215 30 CONTINUE
216 *
217 * Compute B := L*B
218 * where L = P(1)*inv(L(1))* ... *P(m)*inv(L(m)) .
219 *
220 ELSE
221 *
222 * Loop backward applying the transformations to B.
223 *
224 K = N
225 KC = N*( N+1 ) / 2 + 1
226 40 CONTINUE
227 IF( K.LT.1 )
228 $ GO TO 60
229 KC = KC - ( N-K+1 )
230 *
231 * Test the pivot index. If greater than zero, a 1 x 1
232 * pivot was used, otherwise a 2 x 2 pivot was used.
233 *
234 IF( IPIV( K ).GT.0 ) THEN
235 *
236 * 1 x 1 pivot block:
237 *
238 * Multiply by the diagonal element if forming L * D.
239 *
240 IF( NOUNIT )
241 $ CALL DSCAL( NRHS, A( KC ), B( K, 1 ), LDB )
242 *
243 * Multiply by P(K) * inv(L(K)) if K < N.
244 *
245 IF( K.NE.N ) THEN
246 KP = IPIV( K )
247 *
248 * Apply the transformation.
249 *
250 CALL DGER( N-K, NRHS, ONE, A( KC+1 ), 1, B( K, 1 ),
251 $ LDB, B( K+1, 1 ), LDB )
252 *
253 * Interchange if a permutation was applied at the
254 * K-th step of the factorization.
255 *
256 IF( KP.NE.K )
257 $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
258 END IF
259 K = K - 1
260 *
261 ELSE
262 *
263 * 2 x 2 pivot block:
264 *
265 KCNEXT = KC - ( N-K+2 )
266 *
267 * Multiply by the diagonal block if forming L * D.
268 *
269 IF( NOUNIT ) THEN
270 D11 = A( KCNEXT )
271 D22 = A( KC )
272 D21 = A( KCNEXT+1 )
273 D12 = D21
274 DO 50 J = 1, NRHS
275 T1 = B( K-1, J )
276 T2 = B( K, J )
277 B( K-1, J ) = D11*T1 + D12*T2
278 B( K, J ) = D21*T1 + D22*T2
279 50 CONTINUE
280 END IF
281 *
282 * Multiply by P(K) * inv(L(K)) if K < N.
283 *
284 IF( K.NE.N ) THEN
285 *
286 * Apply the transformation.
287 *
288 CALL DGER( N-K, NRHS, ONE, A( KC+1 ), 1, B( K, 1 ),
289 $ LDB, B( K+1, 1 ), LDB )
290 CALL DGER( N-K, NRHS, ONE, A( KCNEXT+2 ), 1,
291 $ B( K-1, 1 ), LDB, B( K+1, 1 ), LDB )
292 *
293 * Interchange if a permutation was applied at the
294 * K-th step of the factorization.
295 *
296 KP = ABS( IPIV( K ) )
297 IF( KP.NE.K )
298 $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
299 END IF
300 KC = KCNEXT
301 K = K - 2
302 END IF
303 GO TO 40
304 60 CONTINUE
305 END IF
306 *----------------------------------------
307 *
308 * Compute B := A' * B (transpose)
309 *
310 *----------------------------------------
311 ELSE
312 *
313 * Form B := U'*B
314 * where U = P(m)*inv(U(m))* ... *P(1)*inv(U(1))
315 * and U' = inv(U'(1))*P(1)* ... *inv(U'(m))*P(m)
316 *
317 IF( LSAME( UPLO, 'U' ) ) THEN
318 *
319 * Loop backward applying the transformations.
320 *
321 K = N
322 KC = N*( N+1 ) / 2 + 1
323 70 CONTINUE
324 IF( K.LT.1 )
325 $ GO TO 90
326 KC = KC - K
327 *
328 * 1 x 1 pivot block.
329 *
330 IF( IPIV( K ).GT.0 ) THEN
331 IF( K.GT.1 ) THEN
332 *
333 * Interchange if P(K) != I.
334 *
335 KP = IPIV( K )
336 IF( KP.NE.K )
337 $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
338 *
339 * Apply the transformation
340 *
341 CALL DGEMV( 'Transpose', K-1, NRHS, ONE, B, LDB,
342 $ A( KC ), 1, ONE, B( K, 1 ), LDB )
343 END IF
344 IF( NOUNIT )
345 $ CALL DSCAL( NRHS, A( KC+K-1 ), B( K, 1 ), LDB )
346 K = K - 1
347 *
348 * 2 x 2 pivot block.
349 *
350 ELSE
351 KCNEXT = KC - ( K-1 )
352 IF( K.GT.2 ) THEN
353 *
354 * Interchange if P(K) != I.
355 *
356 KP = ABS( IPIV( K ) )
357 IF( KP.NE.K-1 )
358 $ CALL DSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ),
359 $ LDB )
360 *
361 * Apply the transformations
362 *
363 CALL DGEMV( 'Transpose', K-2, NRHS, ONE, B, LDB,
364 $ A( KC ), 1, ONE, B( K, 1 ), LDB )
365 CALL DGEMV( 'Transpose', K-2, NRHS, ONE, B, LDB,
366 $ A( KCNEXT ), 1, ONE, B( K-1, 1 ), LDB )
367 END IF
368 *
369 * Multiply by the diagonal block if non-unit.
370 *
371 IF( NOUNIT ) THEN
372 D11 = A( KC-1 )
373 D22 = A( KC+K-1 )
374 D12 = A( KC+K-2 )
375 D21 = D12
376 DO 80 J = 1, NRHS
377 T1 = B( K-1, J )
378 T2 = B( K, J )
379 B( K-1, J ) = D11*T1 + D12*T2
380 B( K, J ) = D21*T1 + D22*T2
381 80 CONTINUE
382 END IF
383 KC = KCNEXT
384 K = K - 2
385 END IF
386 GO TO 70
387 90 CONTINUE
388 *
389 * Form B := L'*B
390 * where L = P(1)*inv(L(1))* ... *P(m)*inv(L(m))
391 * and L' = inv(L(m))*P(m)* ... *inv(L(1))*P(1)
392 *
393 ELSE
394 *
395 * Loop forward applying the L-transformations.
396 *
397 K = 1
398 KC = 1
399 100 CONTINUE
400 IF( K.GT.N )
401 $ GO TO 120
402 *
403 * 1 x 1 pivot block
404 *
405 IF( IPIV( K ).GT.0 ) THEN
406 IF( K.LT.N ) THEN
407 *
408 * Interchange if P(K) != I.
409 *
410 KP = IPIV( K )
411 IF( KP.NE.K )
412 $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
413 *
414 * Apply the transformation
415 *
416 CALL DGEMV( 'Transpose', N-K, NRHS, ONE, B( K+1, 1 ),
417 $ LDB, A( KC+1 ), 1, ONE, B( K, 1 ), LDB )
418 END IF
419 IF( NOUNIT )
420 $ CALL DSCAL( NRHS, A( KC ), B( K, 1 ), LDB )
421 KC = KC + N - K + 1
422 K = K + 1
423 *
424 * 2 x 2 pivot block.
425 *
426 ELSE
427 KCNEXT = KC + N - K + 1
428 IF( K.LT.N-1 ) THEN
429 *
430 * Interchange if P(K) != I.
431 *
432 KP = ABS( IPIV( K ) )
433 IF( KP.NE.K+1 )
434 $ CALL DSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ),
435 $ LDB )
436 *
437 * Apply the transformation
438 *
439 CALL DGEMV( 'Transpose', N-K-1, NRHS, ONE,
440 $ B( K+2, 1 ), LDB, A( KCNEXT+1 ), 1, ONE,
441 $ B( K+1, 1 ), LDB )
442 CALL DGEMV( 'Transpose', N-K-1, NRHS, ONE,
443 $ B( K+2, 1 ), LDB, A( KC+2 ), 1, ONE,
444 $ B( K, 1 ), LDB )
445 END IF
446 *
447 * Multiply by the diagonal block if non-unit.
448 *
449 IF( NOUNIT ) THEN
450 D11 = A( KC )
451 D22 = A( KCNEXT )
452 D21 = A( KC+1 )
453 D12 = D21
454 DO 110 J = 1, NRHS
455 T1 = B( K, J )
456 T2 = B( K+1, J )
457 B( K, J ) = D11*T1 + D12*T2
458 B( K+1, J ) = D21*T1 + D22*T2
459 110 CONTINUE
460 END IF
461 KC = KCNEXT + ( N-K )
462 K = K + 2
463 END IF
464 GO TO 100
465 120 CONTINUE
466 END IF
467 *
468 END IF
469 RETURN
470 *
471 * End of DLAVSP
472 *
473 END