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