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