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