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