1 SUBROUTINE ZLAVSY( 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*16 A( LDA, * ), B( LDB, * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * ZLAVSY 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 ZSYTRF.
24 * ZSYTRF 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', ZLAVSY multiplies either by U or U * D
35 * (or L or L * D).
36 * If TRANS = 'T' or 't', ZLAVSY 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*16 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 ZSYTRF or ZHETRF.
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*16 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*16 ONE
117 PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ) )
118 * ..
119 * .. Local Scalars ..
120 LOGICAL NOUNIT
121 INTEGER J, K, KP
122 COMPLEX*16 D11, D12, D21, D22, T1, T2
123 * ..
124 * .. External Functions ..
125 LOGICAL LSAME
126 EXTERNAL LSAME
127 * ..
128 * .. External Subroutines ..
129 EXTERNAL XERBLA, ZGEMV, ZGERU, ZSCAL, ZSWAP
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( 'ZLAVSY ', -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 ZSCAL( 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 ZGERU( 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 ZSWAP( 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 ZGERU( K-1, NRHS, ONE, A( 1, K ), 1, B( K, 1 ),
234 $ LDB, B( 1, 1 ), LDB )
235 CALL ZGERU( 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 ZSWAP( 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 ZSCAL( 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 ZGERU( N-K, NRHS, ONE, A( K+1, K ), 1, B( K, 1 ),
281 $ 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 ZSWAP( 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 ZGERU( N-K, NRHS, ONE, A( K+1, K ), 1, B( K, 1 ),
317 $ LDB, B( K+1, 1 ), LDB )
318 CALL ZGERU( 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 ZSWAP( 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 CONTINUE
350 IF( K.LT.1 )
351 $ GO TO 90
352 *
353 * 1 x 1 pivot block.
354 *
355 IF( IPIV( K ).GT.0 ) THEN
356 IF( K.GT.1 ) THEN
357 *
358 * Interchange if P(K) != I.
359 *
360 KP = IPIV( K )
361 IF( KP.NE.K )
362 $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
363 *
364 * Apply the transformation
365 *
366 CALL ZGEMV( 'Transpose', K-1, NRHS, ONE, B, LDB,
367 $ A( 1, K ), 1, ONE, B( K, 1 ), LDB )
368 END IF
369 IF( NOUNIT )
370 $ CALL ZSCAL( NRHS, A( K, K ), B( K, 1 ), LDB )
371 K = K - 1
372 *
373 * 2 x 2 pivot block.
374 *
375 ELSE
376 IF( K.GT.2 ) THEN
377 *
378 * Interchange if P(K) != I.
379 *
380 KP = ABS( IPIV( K ) )
381 IF( KP.NE.K-1 )
382 $ CALL ZSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ),
383 $ LDB )
384 *
385 * Apply the transformations
386 *
387 CALL ZGEMV( 'Transpose', K-2, NRHS, ONE, B, LDB,
388 $ A( 1, K ), 1, ONE, B( K, 1 ), LDB )
389 CALL ZGEMV( 'Transpose', K-2, NRHS, ONE, B, LDB,
390 $ A( 1, K-1 ), 1, ONE, B( K-1, 1 ), LDB )
391 END IF
392 *
393 * Multiply by the diagonal block if non-unit.
394 *
395 IF( NOUNIT ) THEN
396 D11 = A( K-1, K-1 )
397 D22 = A( K, K )
398 D12 = A( K-1, K )
399 D21 = D12
400 DO 80 J = 1, NRHS
401 T1 = B( K-1, J )
402 T2 = B( K, J )
403 B( K-1, J ) = D11*T1 + D12*T2
404 B( K, J ) = D21*T1 + D22*T2
405 80 CONTINUE
406 END IF
407 K = K - 2
408 END IF
409 GO TO 70
410 90 CONTINUE
411 *
412 * Form B := L'*B
413 * where L = P(1)*inv(L(1))* ... *P(m)*inv(L(m))
414 * and L' = inv(L'(m))*P(m)* ... *inv(L'(1))*P(1)
415 *
416 ELSE
417 *
418 * Loop forward applying the L-transformations.
419 *
420 K = 1
421 100 CONTINUE
422 IF( K.GT.N )
423 $ GO TO 120
424 *
425 * 1 x 1 pivot block
426 *
427 IF( IPIV( K ).GT.0 ) THEN
428 IF( K.LT.N ) THEN
429 *
430 * Interchange if P(K) != I.
431 *
432 KP = IPIV( K )
433 IF( KP.NE.K )
434 $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
435 *
436 * Apply the transformation
437 *
438 CALL ZGEMV( 'Transpose', N-K, NRHS, ONE, B( K+1, 1 ),
439 $ LDB, A( K+1, K ), 1, ONE, B( K, 1 ), LDB )
440 END IF
441 IF( NOUNIT )
442 $ CALL ZSCAL( NRHS, A( K, K ), B( K, 1 ), LDB )
443 K = K + 1
444 *
445 * 2 x 2 pivot block.
446 *
447 ELSE
448 IF( K.LT.N-1 ) THEN
449 *
450 * Interchange if P(K) != I.
451 *
452 KP = ABS( IPIV( K ) )
453 IF( KP.NE.K+1 )
454 $ CALL ZSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ),
455 $ LDB )
456 *
457 * Apply the transformation
458 *
459 CALL ZGEMV( 'Transpose', N-K-1, NRHS, ONE,
460 $ B( K+2, 1 ), LDB, A( K+2, K+1 ), 1, ONE,
461 $ B( K+1, 1 ), LDB )
462 CALL ZGEMV( 'Transpose', N-K-1, NRHS, ONE,
463 $ B( K+2, 1 ), LDB, A( K+2, K ), 1, ONE,
464 $ B( K, 1 ), LDB )
465 END IF
466 *
467 * Multiply by the diagonal block if non-unit.
468 *
469 IF( NOUNIT ) THEN
470 D11 = A( K, K )
471 D22 = A( K+1, K+1 )
472 D21 = A( K+1, K )
473 D12 = D21
474 DO 110 J = 1, NRHS
475 T1 = B( K, J )
476 T2 = B( K+1, J )
477 B( K, J ) = D11*T1 + D12*T2
478 B( K+1, J ) = D21*T1 + D22*T2
479 110 CONTINUE
480 END IF
481 K = K + 2
482 END IF
483 GO TO 100
484 120 CONTINUE
485 END IF
486 END IF
487 RETURN
488 *
489 * End of ZLAVSY
490 *
491 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*16 A( LDA, * ), B( LDB, * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * ZLAVSY 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 ZSYTRF.
24 * ZSYTRF 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', ZLAVSY multiplies either by U or U * D
35 * (or L or L * D).
36 * If TRANS = 'T' or 't', ZLAVSY 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*16 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 ZSYTRF or ZHETRF.
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*16 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*16 ONE
117 PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ) )
118 * ..
119 * .. Local Scalars ..
120 LOGICAL NOUNIT
121 INTEGER J, K, KP
122 COMPLEX*16 D11, D12, D21, D22, T1, T2
123 * ..
124 * .. External Functions ..
125 LOGICAL LSAME
126 EXTERNAL LSAME
127 * ..
128 * .. External Subroutines ..
129 EXTERNAL XERBLA, ZGEMV, ZGERU, ZSCAL, ZSWAP
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( 'ZLAVSY ', -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 ZSCAL( 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 ZGERU( 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 ZSWAP( 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 ZGERU( K-1, NRHS, ONE, A( 1, K ), 1, B( K, 1 ),
234 $ LDB, B( 1, 1 ), LDB )
235 CALL ZGERU( 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 ZSWAP( 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 ZSCAL( 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 ZGERU( N-K, NRHS, ONE, A( K+1, K ), 1, B( K, 1 ),
281 $ 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 ZSWAP( 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 ZGERU( N-K, NRHS, ONE, A( K+1, K ), 1, B( K, 1 ),
317 $ LDB, B( K+1, 1 ), LDB )
318 CALL ZGERU( 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 ZSWAP( 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 CONTINUE
350 IF( K.LT.1 )
351 $ GO TO 90
352 *
353 * 1 x 1 pivot block.
354 *
355 IF( IPIV( K ).GT.0 ) THEN
356 IF( K.GT.1 ) THEN
357 *
358 * Interchange if P(K) != I.
359 *
360 KP = IPIV( K )
361 IF( KP.NE.K )
362 $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
363 *
364 * Apply the transformation
365 *
366 CALL ZGEMV( 'Transpose', K-1, NRHS, ONE, B, LDB,
367 $ A( 1, K ), 1, ONE, B( K, 1 ), LDB )
368 END IF
369 IF( NOUNIT )
370 $ CALL ZSCAL( NRHS, A( K, K ), B( K, 1 ), LDB )
371 K = K - 1
372 *
373 * 2 x 2 pivot block.
374 *
375 ELSE
376 IF( K.GT.2 ) THEN
377 *
378 * Interchange if P(K) != I.
379 *
380 KP = ABS( IPIV( K ) )
381 IF( KP.NE.K-1 )
382 $ CALL ZSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ),
383 $ LDB )
384 *
385 * Apply the transformations
386 *
387 CALL ZGEMV( 'Transpose', K-2, NRHS, ONE, B, LDB,
388 $ A( 1, K ), 1, ONE, B( K, 1 ), LDB )
389 CALL ZGEMV( 'Transpose', K-2, NRHS, ONE, B, LDB,
390 $ A( 1, K-1 ), 1, ONE, B( K-1, 1 ), LDB )
391 END IF
392 *
393 * Multiply by the diagonal block if non-unit.
394 *
395 IF( NOUNIT ) THEN
396 D11 = A( K-1, K-1 )
397 D22 = A( K, K )
398 D12 = A( K-1, K )
399 D21 = D12
400 DO 80 J = 1, NRHS
401 T1 = B( K-1, J )
402 T2 = B( K, J )
403 B( K-1, J ) = D11*T1 + D12*T2
404 B( K, J ) = D21*T1 + D22*T2
405 80 CONTINUE
406 END IF
407 K = K - 2
408 END IF
409 GO TO 70
410 90 CONTINUE
411 *
412 * Form B := L'*B
413 * where L = P(1)*inv(L(1))* ... *P(m)*inv(L(m))
414 * and L' = inv(L'(m))*P(m)* ... *inv(L'(1))*P(1)
415 *
416 ELSE
417 *
418 * Loop forward applying the L-transformations.
419 *
420 K = 1
421 100 CONTINUE
422 IF( K.GT.N )
423 $ GO TO 120
424 *
425 * 1 x 1 pivot block
426 *
427 IF( IPIV( K ).GT.0 ) THEN
428 IF( K.LT.N ) THEN
429 *
430 * Interchange if P(K) != I.
431 *
432 KP = IPIV( K )
433 IF( KP.NE.K )
434 $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
435 *
436 * Apply the transformation
437 *
438 CALL ZGEMV( 'Transpose', N-K, NRHS, ONE, B( K+1, 1 ),
439 $ LDB, A( K+1, K ), 1, ONE, B( K, 1 ), LDB )
440 END IF
441 IF( NOUNIT )
442 $ CALL ZSCAL( NRHS, A( K, K ), B( K, 1 ), LDB )
443 K = K + 1
444 *
445 * 2 x 2 pivot block.
446 *
447 ELSE
448 IF( K.LT.N-1 ) THEN
449 *
450 * Interchange if P(K) != I.
451 *
452 KP = ABS( IPIV( K ) )
453 IF( KP.NE.K+1 )
454 $ CALL ZSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ),
455 $ LDB )
456 *
457 * Apply the transformation
458 *
459 CALL ZGEMV( 'Transpose', N-K-1, NRHS, ONE,
460 $ B( K+2, 1 ), LDB, A( K+2, K+1 ), 1, ONE,
461 $ B( K+1, 1 ), LDB )
462 CALL ZGEMV( 'Transpose', N-K-1, NRHS, ONE,
463 $ B( K+2, 1 ), LDB, A( K+2, K ), 1, ONE,
464 $ B( K, 1 ), LDB )
465 END IF
466 *
467 * Multiply by the diagonal block if non-unit.
468 *
469 IF( NOUNIT ) THEN
470 D11 = A( K, K )
471 D22 = A( K+1, K+1 )
472 D21 = A( K+1, K )
473 D12 = D21
474 DO 110 J = 1, NRHS
475 T1 = B( K, J )
476 T2 = B( K+1, J )
477 B( K, J ) = D11*T1 + D12*T2
478 B( K+1, J ) = D21*T1 + D22*T2
479 110 CONTINUE
480 END IF
481 K = K + 2
482 END IF
483 GO TO 100
484 120 CONTINUE
485 END IF
486 END IF
487 RETURN
488 *
489 * End of ZLAVSY
490 *
491 END