1 SUBROUTINE ZLAHEF( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, INFO )
2 *
3 * -- LAPACK routine (version 3.3.1) --
4 * -- LAPACK is a software package provided by Univ. of Tennessee, --
5 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6 * -- April 2011 --
7 *
8 * .. Scalar Arguments ..
9 CHARACTER UPLO
10 INTEGER INFO, KB, LDA, LDW, N, NB
11 * ..
12 * .. Array Arguments ..
13 INTEGER IPIV( * )
14 COMPLEX*16 A( LDA, * ), W( LDW, * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * ZLAHEF computes a partial factorization of a complex Hermitian
21 * matrix A using the Bunch-Kaufman diagonal pivoting method. The
22 * partial factorization has the form:
23 *
24 * A = ( I U12 ) ( A11 0 ) ( I 0 ) if UPLO = 'U', or:
25 * ( 0 U22 ) ( 0 D ) ( U12**H U22**H )
26 *
27 * A = ( L11 0 ) ( D 0 ) ( L11**H L21**H ) if UPLO = 'L'
28 * ( L21 I ) ( 0 A22 ) ( 0 I )
29 *
30 * where the order of D is at most NB. The actual order is returned in
31 * the argument KB, and is either NB or NB-1, or N if N <= NB.
32 * Note that U**H denotes the conjugate transpose of U.
33 *
34 * ZLAHEF is an auxiliary routine called by ZHETRF. It uses blocked code
35 * (calling Level 3 BLAS) to update the submatrix A11 (if UPLO = 'U') or
36 * A22 (if UPLO = 'L').
37 *
38 * Arguments
39 * =========
40 *
41 * UPLO (input) CHARACTER*1
42 * Specifies whether the upper or lower triangular part of the
43 * Hermitian matrix A is stored:
44 * = 'U': Upper triangular
45 * = 'L': Lower triangular
46 *
47 * N (input) INTEGER
48 * The order of the matrix A. N >= 0.
49 *
50 * NB (input) INTEGER
51 * The maximum number of columns of the matrix A that should be
52 * factored. NB should be at least 2 to allow for 2-by-2 pivot
53 * blocks.
54 *
55 * KB (output) INTEGER
56 * The number of columns of A that were actually factored.
57 * KB is either NB-1 or NB, or N if N <= NB.
58 *
59 * A (input/output) COMPLEX*16 array, dimension (LDA,N)
60 * On entry, the Hermitian matrix A. If UPLO = 'U', the leading
61 * n-by-n upper triangular part of A contains the upper
62 * triangular part of the matrix A, and the strictly lower
63 * triangular part of A is not referenced. If UPLO = 'L', the
64 * leading n-by-n lower triangular part of A contains the lower
65 * triangular part of the matrix A, and the strictly upper
66 * triangular part of A is not referenced.
67 * On exit, A contains details of the partial factorization.
68 *
69 * LDA (input) INTEGER
70 * The leading dimension of the array A. LDA >= max(1,N).
71 *
72 * IPIV (output) INTEGER array, dimension (N)
73 * Details of the interchanges and the block structure of D.
74 * If UPLO = 'U', only the last KB elements of IPIV are set;
75 * if UPLO = 'L', only the first KB elements are set.
76 *
77 * If IPIV(k) > 0, then rows and columns k and IPIV(k) were
78 * interchanged and D(k,k) is a 1-by-1 diagonal block.
79 * If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and
80 * columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
81 * is a 2-by-2 diagonal block. If UPLO = 'L' and IPIV(k) =
82 * IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were
83 * interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
84 *
85 * W (workspace) COMPLEX*16 array, dimension (LDW,NB)
86 *
87 * LDW (input) INTEGER
88 * The leading dimension of the array W. LDW >= max(1,N).
89 *
90 * INFO (output) INTEGER
91 * = 0: successful exit
92 * > 0: if INFO = k, D(k,k) is exactly zero. The factorization
93 * has been completed, but the block diagonal matrix D is
94 * exactly singular.
95 *
96 * =====================================================================
97 *
98 * .. Parameters ..
99 DOUBLE PRECISION ZERO, ONE
100 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
101 COMPLEX*16 CONE
102 PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ) )
103 DOUBLE PRECISION EIGHT, SEVTEN
104 PARAMETER ( EIGHT = 8.0D+0, SEVTEN = 17.0D+0 )
105 * ..
106 * .. Local Scalars ..
107 INTEGER IMAX, J, JB, JJ, JMAX, JP, K, KK, KKW, KP,
108 $ KSTEP, KW
109 DOUBLE PRECISION ABSAKK, ALPHA, COLMAX, R1, ROWMAX, T
110 COMPLEX*16 D11, D21, D22, Z
111 * ..
112 * .. External Functions ..
113 LOGICAL LSAME
114 INTEGER IZAMAX
115 EXTERNAL LSAME, IZAMAX
116 * ..
117 * .. External Subroutines ..
118 EXTERNAL ZCOPY, ZDSCAL, ZGEMM, ZGEMV, ZLACGV, ZSWAP
119 * ..
120 * .. Intrinsic Functions ..
121 INTRINSIC ABS, DBLE, DCONJG, DIMAG, MAX, MIN, SQRT
122 * ..
123 * .. Statement Functions ..
124 DOUBLE PRECISION CABS1
125 * ..
126 * .. Statement Function definitions ..
127 CABS1( Z ) = ABS( DBLE( Z ) ) + ABS( DIMAG( Z ) )
128 * ..
129 * .. Executable Statements ..
130 *
131 INFO = 0
132 *
133 * Initialize ALPHA for use in choosing pivot block size.
134 *
135 ALPHA = ( ONE+SQRT( SEVTEN ) ) / EIGHT
136 *
137 IF( LSAME( UPLO, 'U' ) ) THEN
138 *
139 * Factorize the trailing columns of A using the upper triangle
140 * of A and working backwards, and compute the matrix W = U12*D
141 * for use in updating A11 (note that conjg(W) is actually stored)
142 *
143 * K is the main loop index, decreasing from N in steps of 1 or 2
144 *
145 * KW is the column of W which corresponds to column K of A
146 *
147 K = N
148 10 CONTINUE
149 KW = NB + K - N
150 *
151 * Exit from loop
152 *
153 IF( ( K.LE.N-NB+1 .AND. NB.LT.N ) .OR. K.LT.1 )
154 $ GO TO 30
155 *
156 * Copy column K of A to column KW of W and update it
157 *
158 CALL ZCOPY( K-1, A( 1, K ), 1, W( 1, KW ), 1 )
159 W( K, KW ) = DBLE( A( K, K ) )
160 IF( K.LT.N ) THEN
161 CALL ZGEMV( 'No transpose', K, N-K, -CONE, A( 1, K+1 ), LDA,
162 $ W( K, KW+1 ), LDW, CONE, W( 1, KW ), 1 )
163 W( K, KW ) = DBLE( W( K, KW ) )
164 END IF
165 *
166 KSTEP = 1
167 *
168 * Determine rows and columns to be interchanged and whether
169 * a 1-by-1 or 2-by-2 pivot block will be used
170 *
171 ABSAKK = ABS( DBLE( W( K, KW ) ) )
172 *
173 * IMAX is the row-index of the largest off-diagonal element in
174 * column K, and COLMAX is its absolute value
175 *
176 IF( K.GT.1 ) THEN
177 IMAX = IZAMAX( K-1, W( 1, KW ), 1 )
178 COLMAX = CABS1( W( IMAX, KW ) )
179 ELSE
180 COLMAX = ZERO
181 END IF
182 *
183 IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
184 *
185 * Column K is zero: set INFO and continue
186 *
187 IF( INFO.EQ.0 )
188 $ INFO = K
189 KP = K
190 A( K, K ) = DBLE( A( K, K ) )
191 ELSE
192 IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
193 *
194 * no interchange, use 1-by-1 pivot block
195 *
196 KP = K
197 ELSE
198 *
199 * Copy column IMAX to column KW-1 of W and update it
200 *
201 CALL ZCOPY( IMAX-1, A( 1, IMAX ), 1, W( 1, KW-1 ), 1 )
202 W( IMAX, KW-1 ) = DBLE( A( IMAX, IMAX ) )
203 CALL ZCOPY( K-IMAX, A( IMAX, IMAX+1 ), LDA,
204 $ W( IMAX+1, KW-1 ), 1 )
205 CALL ZLACGV( K-IMAX, W( IMAX+1, KW-1 ), 1 )
206 IF( K.LT.N ) THEN
207 CALL ZGEMV( 'No transpose', K, N-K, -CONE,
208 $ A( 1, K+1 ), LDA, W( IMAX, KW+1 ), LDW,
209 $ CONE, W( 1, KW-1 ), 1 )
210 W( IMAX, KW-1 ) = DBLE( W( IMAX, KW-1 ) )
211 END IF
212 *
213 * JMAX is the column-index of the largest off-diagonal
214 * element in row IMAX, and ROWMAX is its absolute value
215 *
216 JMAX = IMAX + IZAMAX( K-IMAX, W( IMAX+1, KW-1 ), 1 )
217 ROWMAX = CABS1( W( JMAX, KW-1 ) )
218 IF( IMAX.GT.1 ) THEN
219 JMAX = IZAMAX( IMAX-1, W( 1, KW-1 ), 1 )
220 ROWMAX = MAX( ROWMAX, CABS1( W( JMAX, KW-1 ) ) )
221 END IF
222 *
223 IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
224 *
225 * no interchange, use 1-by-1 pivot block
226 *
227 KP = K
228 ELSE IF( ABS( DBLE( W( IMAX, KW-1 ) ) ).GE.ALPHA*ROWMAX )
229 $ THEN
230 *
231 * interchange rows and columns K and IMAX, use 1-by-1
232 * pivot block
233 *
234 KP = IMAX
235 *
236 * copy column KW-1 of W to column KW
237 *
238 CALL ZCOPY( K, W( 1, KW-1 ), 1, W( 1, KW ), 1 )
239 ELSE
240 *
241 * interchange rows and columns K-1 and IMAX, use 2-by-2
242 * pivot block
243 *
244 KP = IMAX
245 KSTEP = 2
246 END IF
247 END IF
248 *
249 KK = K - KSTEP + 1
250 KKW = NB + KK - N
251 *
252 * Updated column KP is already stored in column KKW of W
253 *
254 IF( KP.NE.KK ) THEN
255 *
256 * Copy non-updated column KK to column KP
257 *
258 A( KP, KP ) = DBLE( A( KK, KK ) )
259 CALL ZCOPY( KK-1-KP, A( KP+1, KK ), 1, A( KP, KP+1 ),
260 $ LDA )
261 CALL ZLACGV( KK-1-KP, A( KP, KP+1 ), LDA )
262 CALL ZCOPY( KP-1, A( 1, KK ), 1, A( 1, KP ), 1 )
263 *
264 * Interchange rows KK and KP in last KK columns of A and W
265 *
266 IF( KK.LT.N )
267 $ CALL ZSWAP( N-KK, A( KK, KK+1 ), LDA, A( KP, KK+1 ),
268 $ LDA )
269 CALL ZSWAP( N-KK+1, W( KK, KKW ), LDW, W( KP, KKW ),
270 $ LDW )
271 END IF
272 *
273 IF( KSTEP.EQ.1 ) THEN
274 *
275 * 1-by-1 pivot block D(k): column KW of W now holds
276 *
277 * W(k) = U(k)*D(k)
278 *
279 * where U(k) is the k-th column of U
280 *
281 * Store U(k) in column k of A
282 *
283 CALL ZCOPY( K, W( 1, KW ), 1, A( 1, K ), 1 )
284 R1 = ONE / DBLE( A( K, K ) )
285 CALL ZDSCAL( K-1, R1, A( 1, K ), 1 )
286 *
287 * Conjugate W(k)
288 *
289 CALL ZLACGV( K-1, W( 1, KW ), 1 )
290 ELSE
291 *
292 * 2-by-2 pivot block D(k): columns KW and KW-1 of W now
293 * hold
294 *
295 * ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
296 *
297 * where U(k) and U(k-1) are the k-th and (k-1)-th columns
298 * of U
299 *
300 IF( K.GT.2 ) THEN
301 *
302 * Store U(k) and U(k-1) in columns k and k-1 of A
303 *
304 D21 = W( K-1, KW )
305 D11 = W( K, KW ) / DCONJG( D21 )
306 D22 = W( K-1, KW-1 ) / D21
307 T = ONE / ( DBLE( D11*D22 )-ONE )
308 D21 = T / D21
309 DO 20 J = 1, K - 2
310 A( J, K-1 ) = D21*( D11*W( J, KW-1 )-W( J, KW ) )
311 A( J, K ) = DCONJG( D21 )*
312 $ ( D22*W( J, KW )-W( J, KW-1 ) )
313 20 CONTINUE
314 END IF
315 *
316 * Copy D(k) to A
317 *
318 A( K-1, K-1 ) = W( K-1, KW-1 )
319 A( K-1, K ) = W( K-1, KW )
320 A( K, K ) = W( K, KW )
321 *
322 * Conjugate W(k) and W(k-1)
323 *
324 CALL ZLACGV( K-1, W( 1, KW ), 1 )
325 CALL ZLACGV( K-2, W( 1, KW-1 ), 1 )
326 END IF
327 END IF
328 *
329 * Store details of the interchanges in IPIV
330 *
331 IF( KSTEP.EQ.1 ) THEN
332 IPIV( K ) = KP
333 ELSE
334 IPIV( K ) = -KP
335 IPIV( K-1 ) = -KP
336 END IF
337 *
338 * Decrease K and return to the start of the main loop
339 *
340 K = K - KSTEP
341 GO TO 10
342 *
343 30 CONTINUE
344 *
345 * Update the upper triangle of A11 (= A(1:k,1:k)) as
346 *
347 * A11 := A11 - U12*D*U12**H = A11 - U12*W**H
348 *
349 * computing blocks of NB columns at a time (note that conjg(W) is
350 * actually stored)
351 *
352 DO 50 J = ( ( K-1 ) / NB )*NB + 1, 1, -NB
353 JB = MIN( NB, K-J+1 )
354 *
355 * Update the upper triangle of the diagonal block
356 *
357 DO 40 JJ = J, J + JB - 1
358 A( JJ, JJ ) = DBLE( A( JJ, JJ ) )
359 CALL ZGEMV( 'No transpose', JJ-J+1, N-K, -CONE,
360 $ A( J, K+1 ), LDA, W( JJ, KW+1 ), LDW, CONE,
361 $ A( J, JJ ), 1 )
362 A( JJ, JJ ) = DBLE( A( JJ, JJ ) )
363 40 CONTINUE
364 *
365 * Update the rectangular superdiagonal block
366 *
367 CALL ZGEMM( 'No transpose', 'Transpose', J-1, JB, N-K,
368 $ -CONE, A( 1, K+1 ), LDA, W( J, KW+1 ), LDW,
369 $ CONE, A( 1, J ), LDA )
370 50 CONTINUE
371 *
372 * Put U12 in standard form by partially undoing the interchanges
373 * in columns k+1:n
374 *
375 J = K + 1
376 60 CONTINUE
377 JJ = J
378 JP = IPIV( J )
379 IF( JP.LT.0 ) THEN
380 JP = -JP
381 J = J + 1
382 END IF
383 J = J + 1
384 IF( JP.NE.JJ .AND. J.LE.N )
385 $ CALL ZSWAP( N-J+1, A( JP, J ), LDA, A( JJ, J ), LDA )
386 IF( J.LE.N )
387 $ GO TO 60
388 *
389 * Set KB to the number of columns factorized
390 *
391 KB = N - K
392 *
393 ELSE
394 *
395 * Factorize the leading columns of A using the lower triangle
396 * of A and working forwards, and compute the matrix W = L21*D
397 * for use in updating A22 (note that conjg(W) is actually stored)
398 *
399 * K is the main loop index, increasing from 1 in steps of 1 or 2
400 *
401 K = 1
402 70 CONTINUE
403 *
404 * Exit from loop
405 *
406 IF( ( K.GE.NB .AND. NB.LT.N ) .OR. K.GT.N )
407 $ GO TO 90
408 *
409 * Copy column K of A to column K of W and update it
410 *
411 W( K, K ) = DBLE( A( K, K ) )
412 IF( K.LT.N )
413 $ CALL ZCOPY( N-K, A( K+1, K ), 1, W( K+1, K ), 1 )
414 CALL ZGEMV( 'No transpose', N-K+1, K-1, -CONE, A( K, 1 ), LDA,
415 $ W( K, 1 ), LDW, CONE, W( K, K ), 1 )
416 W( K, K ) = DBLE( W( K, K ) )
417 *
418 KSTEP = 1
419 *
420 * Determine rows and columns to be interchanged and whether
421 * a 1-by-1 or 2-by-2 pivot block will be used
422 *
423 ABSAKK = ABS( DBLE( W( K, K ) ) )
424 *
425 * IMAX is the row-index of the largest off-diagonal element in
426 * column K, and COLMAX is its absolute value
427 *
428 IF( K.LT.N ) THEN
429 IMAX = K + IZAMAX( N-K, W( K+1, K ), 1 )
430 COLMAX = CABS1( W( IMAX, K ) )
431 ELSE
432 COLMAX = ZERO
433 END IF
434 *
435 IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
436 *
437 * Column K is zero: set INFO and continue
438 *
439 IF( INFO.EQ.0 )
440 $ INFO = K
441 KP = K
442 A( K, K ) = DBLE( A( K, K ) )
443 ELSE
444 IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
445 *
446 * no interchange, use 1-by-1 pivot block
447 *
448 KP = K
449 ELSE
450 *
451 * Copy column IMAX to column K+1 of W and update it
452 *
453 CALL ZCOPY( IMAX-K, A( IMAX, K ), LDA, W( K, K+1 ), 1 )
454 CALL ZLACGV( IMAX-K, W( K, K+1 ), 1 )
455 W( IMAX, K+1 ) = DBLE( A( IMAX, IMAX ) )
456 IF( IMAX.LT.N )
457 $ CALL ZCOPY( N-IMAX, A( IMAX+1, IMAX ), 1,
458 $ W( IMAX+1, K+1 ), 1 )
459 CALL ZGEMV( 'No transpose', N-K+1, K-1, -CONE, A( K, 1 ),
460 $ LDA, W( IMAX, 1 ), LDW, CONE, W( K, K+1 ),
461 $ 1 )
462 W( IMAX, K+1 ) = DBLE( W( IMAX, K+1 ) )
463 *
464 * JMAX is the column-index of the largest off-diagonal
465 * element in row IMAX, and ROWMAX is its absolute value
466 *
467 JMAX = K - 1 + IZAMAX( IMAX-K, W( K, K+1 ), 1 )
468 ROWMAX = CABS1( W( JMAX, K+1 ) )
469 IF( IMAX.LT.N ) THEN
470 JMAX = IMAX + IZAMAX( N-IMAX, W( IMAX+1, K+1 ), 1 )
471 ROWMAX = MAX( ROWMAX, CABS1( W( JMAX, K+1 ) ) )
472 END IF
473 *
474 IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
475 *
476 * no interchange, use 1-by-1 pivot block
477 *
478 KP = K
479 ELSE IF( ABS( DBLE( W( IMAX, K+1 ) ) ).GE.ALPHA*ROWMAX )
480 $ THEN
481 *
482 * interchange rows and columns K and IMAX, use 1-by-1
483 * pivot block
484 *
485 KP = IMAX
486 *
487 * copy column K+1 of W to column K
488 *
489 CALL ZCOPY( N-K+1, W( K, K+1 ), 1, W( K, K ), 1 )
490 ELSE
491 *
492 * interchange rows and columns K+1 and IMAX, use 2-by-2
493 * pivot block
494 *
495 KP = IMAX
496 KSTEP = 2
497 END IF
498 END IF
499 *
500 KK = K + KSTEP - 1
501 *
502 * Updated column KP is already stored in column KK of W
503 *
504 IF( KP.NE.KK ) THEN
505 *
506 * Copy non-updated column KK to column KP
507 *
508 A( KP, KP ) = DBLE( A( KK, KK ) )
509 CALL ZCOPY( KP-KK-1, A( KK+1, KK ), 1, A( KP, KK+1 ),
510 $ LDA )
511 CALL ZLACGV( KP-KK-1, A( KP, KK+1 ), LDA )
512 IF( KP.LT.N )
513 $ CALL ZCOPY( N-KP, A( KP+1, KK ), 1, A( KP+1, KP ), 1 )
514 *
515 * Interchange rows KK and KP in first KK columns of A and W
516 *
517 CALL ZSWAP( KK-1, A( KK, 1 ), LDA, A( KP, 1 ), LDA )
518 CALL ZSWAP( KK, W( KK, 1 ), LDW, W( KP, 1 ), LDW )
519 END IF
520 *
521 IF( KSTEP.EQ.1 ) THEN
522 *
523 * 1-by-1 pivot block D(k): column k of W now holds
524 *
525 * W(k) = L(k)*D(k)
526 *
527 * where L(k) is the k-th column of L
528 *
529 * Store L(k) in column k of A
530 *
531 CALL ZCOPY( N-K+1, W( K, K ), 1, A( K, K ), 1 )
532 IF( K.LT.N ) THEN
533 R1 = ONE / DBLE( A( K, K ) )
534 CALL ZDSCAL( N-K, R1, A( K+1, K ), 1 )
535 *
536 * Conjugate W(k)
537 *
538 CALL ZLACGV( N-K, W( K+1, K ), 1 )
539 END IF
540 ELSE
541 *
542 * 2-by-2 pivot block D(k): columns k and k+1 of W now hold
543 *
544 * ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k)
545 *
546 * where L(k) and L(k+1) are the k-th and (k+1)-th columns
547 * of L
548 *
549 IF( K.LT.N-1 ) THEN
550 *
551 * Store L(k) and L(k+1) in columns k and k+1 of A
552 *
553 D21 = W( K+1, K )
554 D11 = W( K+1, K+1 ) / D21
555 D22 = W( K, K ) / DCONJG( D21 )
556 T = ONE / ( DBLE( D11*D22 )-ONE )
557 D21 = T / D21
558 DO 80 J = K + 2, N
559 A( J, K ) = DCONJG( D21 )*
560 $ ( D11*W( J, K )-W( J, K+1 ) )
561 A( J, K+1 ) = D21*( D22*W( J, K+1 )-W( J, K ) )
562 80 CONTINUE
563 END IF
564 *
565 * Copy D(k) to A
566 *
567 A( K, K ) = W( K, K )
568 A( K+1, K ) = W( K+1, K )
569 A( K+1, K+1 ) = W( K+1, K+1 )
570 *
571 * Conjugate W(k) and W(k+1)
572 *
573 CALL ZLACGV( N-K, W( K+1, K ), 1 )
574 CALL ZLACGV( N-K-1, W( K+2, K+1 ), 1 )
575 END IF
576 END IF
577 *
578 * Store details of the interchanges in IPIV
579 *
580 IF( KSTEP.EQ.1 ) THEN
581 IPIV( K ) = KP
582 ELSE
583 IPIV( K ) = -KP
584 IPIV( K+1 ) = -KP
585 END IF
586 *
587 * Increase K and return to the start of the main loop
588 *
589 K = K + KSTEP
590 GO TO 70
591 *
592 90 CONTINUE
593 *
594 * Update the lower triangle of A22 (= A(k:n,k:n)) as
595 *
596 * A22 := A22 - L21*D*L21**H = A22 - L21*W**H
597 *
598 * computing blocks of NB columns at a time (note that conjg(W) is
599 * actually stored)
600 *
601 DO 110 J = K, N, NB
602 JB = MIN( NB, N-J+1 )
603 *
604 * Update the lower triangle of the diagonal block
605 *
606 DO 100 JJ = J, J + JB - 1
607 A( JJ, JJ ) = DBLE( A( JJ, JJ ) )
608 CALL ZGEMV( 'No transpose', J+JB-JJ, K-1, -CONE,
609 $ A( JJ, 1 ), LDA, W( JJ, 1 ), LDW, CONE,
610 $ A( JJ, JJ ), 1 )
611 A( JJ, JJ ) = DBLE( A( JJ, JJ ) )
612 100 CONTINUE
613 *
614 * Update the rectangular subdiagonal block
615 *
616 IF( J+JB.LE.N )
617 $ CALL ZGEMM( 'No transpose', 'Transpose', N-J-JB+1, JB,
618 $ K-1, -CONE, A( J+JB, 1 ), LDA, W( J, 1 ),
619 $ LDW, CONE, A( J+JB, J ), LDA )
620 110 CONTINUE
621 *
622 * Put L21 in standard form by partially undoing the interchanges
623 * in columns 1:k-1
624 *
625 J = K - 1
626 120 CONTINUE
627 JJ = J
628 JP = IPIV( J )
629 IF( JP.LT.0 ) THEN
630 JP = -JP
631 J = J - 1
632 END IF
633 J = J - 1
634 IF( JP.NE.JJ .AND. J.GE.1 )
635 $ CALL ZSWAP( J, A( JP, 1 ), LDA, A( JJ, 1 ), LDA )
636 IF( J.GE.1 )
637 $ GO TO 120
638 *
639 * Set KB to the number of columns factorized
640 *
641 KB = K - 1
642 *
643 END IF
644 RETURN
645 *
646 * End of ZLAHEF
647 *
648 END
2 *
3 * -- LAPACK routine (version 3.3.1) --
4 * -- LAPACK is a software package provided by Univ. of Tennessee, --
5 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6 * -- April 2011 --
7 *
8 * .. Scalar Arguments ..
9 CHARACTER UPLO
10 INTEGER INFO, KB, LDA, LDW, N, NB
11 * ..
12 * .. Array Arguments ..
13 INTEGER IPIV( * )
14 COMPLEX*16 A( LDA, * ), W( LDW, * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * ZLAHEF computes a partial factorization of a complex Hermitian
21 * matrix A using the Bunch-Kaufman diagonal pivoting method. The
22 * partial factorization has the form:
23 *
24 * A = ( I U12 ) ( A11 0 ) ( I 0 ) if UPLO = 'U', or:
25 * ( 0 U22 ) ( 0 D ) ( U12**H U22**H )
26 *
27 * A = ( L11 0 ) ( D 0 ) ( L11**H L21**H ) if UPLO = 'L'
28 * ( L21 I ) ( 0 A22 ) ( 0 I )
29 *
30 * where the order of D is at most NB. The actual order is returned in
31 * the argument KB, and is either NB or NB-1, or N if N <= NB.
32 * Note that U**H denotes the conjugate transpose of U.
33 *
34 * ZLAHEF is an auxiliary routine called by ZHETRF. It uses blocked code
35 * (calling Level 3 BLAS) to update the submatrix A11 (if UPLO = 'U') or
36 * A22 (if UPLO = 'L').
37 *
38 * Arguments
39 * =========
40 *
41 * UPLO (input) CHARACTER*1
42 * Specifies whether the upper or lower triangular part of the
43 * Hermitian matrix A is stored:
44 * = 'U': Upper triangular
45 * = 'L': Lower triangular
46 *
47 * N (input) INTEGER
48 * The order of the matrix A. N >= 0.
49 *
50 * NB (input) INTEGER
51 * The maximum number of columns of the matrix A that should be
52 * factored. NB should be at least 2 to allow for 2-by-2 pivot
53 * blocks.
54 *
55 * KB (output) INTEGER
56 * The number of columns of A that were actually factored.
57 * KB is either NB-1 or NB, or N if N <= NB.
58 *
59 * A (input/output) COMPLEX*16 array, dimension (LDA,N)
60 * On entry, the Hermitian matrix A. If UPLO = 'U', the leading
61 * n-by-n upper triangular part of A contains the upper
62 * triangular part of the matrix A, and the strictly lower
63 * triangular part of A is not referenced. If UPLO = 'L', the
64 * leading n-by-n lower triangular part of A contains the lower
65 * triangular part of the matrix A, and the strictly upper
66 * triangular part of A is not referenced.
67 * On exit, A contains details of the partial factorization.
68 *
69 * LDA (input) INTEGER
70 * The leading dimension of the array A. LDA >= max(1,N).
71 *
72 * IPIV (output) INTEGER array, dimension (N)
73 * Details of the interchanges and the block structure of D.
74 * If UPLO = 'U', only the last KB elements of IPIV are set;
75 * if UPLO = 'L', only the first KB elements are set.
76 *
77 * If IPIV(k) > 0, then rows and columns k and IPIV(k) were
78 * interchanged and D(k,k) is a 1-by-1 diagonal block.
79 * If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and
80 * columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
81 * is a 2-by-2 diagonal block. If UPLO = 'L' and IPIV(k) =
82 * IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were
83 * interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
84 *
85 * W (workspace) COMPLEX*16 array, dimension (LDW,NB)
86 *
87 * LDW (input) INTEGER
88 * The leading dimension of the array W. LDW >= max(1,N).
89 *
90 * INFO (output) INTEGER
91 * = 0: successful exit
92 * > 0: if INFO = k, D(k,k) is exactly zero. The factorization
93 * has been completed, but the block diagonal matrix D is
94 * exactly singular.
95 *
96 * =====================================================================
97 *
98 * .. Parameters ..
99 DOUBLE PRECISION ZERO, ONE
100 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
101 COMPLEX*16 CONE
102 PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ) )
103 DOUBLE PRECISION EIGHT, SEVTEN
104 PARAMETER ( EIGHT = 8.0D+0, SEVTEN = 17.0D+0 )
105 * ..
106 * .. Local Scalars ..
107 INTEGER IMAX, J, JB, JJ, JMAX, JP, K, KK, KKW, KP,
108 $ KSTEP, KW
109 DOUBLE PRECISION ABSAKK, ALPHA, COLMAX, R1, ROWMAX, T
110 COMPLEX*16 D11, D21, D22, Z
111 * ..
112 * .. External Functions ..
113 LOGICAL LSAME
114 INTEGER IZAMAX
115 EXTERNAL LSAME, IZAMAX
116 * ..
117 * .. External Subroutines ..
118 EXTERNAL ZCOPY, ZDSCAL, ZGEMM, ZGEMV, ZLACGV, ZSWAP
119 * ..
120 * .. Intrinsic Functions ..
121 INTRINSIC ABS, DBLE, DCONJG, DIMAG, MAX, MIN, SQRT
122 * ..
123 * .. Statement Functions ..
124 DOUBLE PRECISION CABS1
125 * ..
126 * .. Statement Function definitions ..
127 CABS1( Z ) = ABS( DBLE( Z ) ) + ABS( DIMAG( Z ) )
128 * ..
129 * .. Executable Statements ..
130 *
131 INFO = 0
132 *
133 * Initialize ALPHA for use in choosing pivot block size.
134 *
135 ALPHA = ( ONE+SQRT( SEVTEN ) ) / EIGHT
136 *
137 IF( LSAME( UPLO, 'U' ) ) THEN
138 *
139 * Factorize the trailing columns of A using the upper triangle
140 * of A and working backwards, and compute the matrix W = U12*D
141 * for use in updating A11 (note that conjg(W) is actually stored)
142 *
143 * K is the main loop index, decreasing from N in steps of 1 or 2
144 *
145 * KW is the column of W which corresponds to column K of A
146 *
147 K = N
148 10 CONTINUE
149 KW = NB + K - N
150 *
151 * Exit from loop
152 *
153 IF( ( K.LE.N-NB+1 .AND. NB.LT.N ) .OR. K.LT.1 )
154 $ GO TO 30
155 *
156 * Copy column K of A to column KW of W and update it
157 *
158 CALL ZCOPY( K-1, A( 1, K ), 1, W( 1, KW ), 1 )
159 W( K, KW ) = DBLE( A( K, K ) )
160 IF( K.LT.N ) THEN
161 CALL ZGEMV( 'No transpose', K, N-K, -CONE, A( 1, K+1 ), LDA,
162 $ W( K, KW+1 ), LDW, CONE, W( 1, KW ), 1 )
163 W( K, KW ) = DBLE( W( K, KW ) )
164 END IF
165 *
166 KSTEP = 1
167 *
168 * Determine rows and columns to be interchanged and whether
169 * a 1-by-1 or 2-by-2 pivot block will be used
170 *
171 ABSAKK = ABS( DBLE( W( K, KW ) ) )
172 *
173 * IMAX is the row-index of the largest off-diagonal element in
174 * column K, and COLMAX is its absolute value
175 *
176 IF( K.GT.1 ) THEN
177 IMAX = IZAMAX( K-1, W( 1, KW ), 1 )
178 COLMAX = CABS1( W( IMAX, KW ) )
179 ELSE
180 COLMAX = ZERO
181 END IF
182 *
183 IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
184 *
185 * Column K is zero: set INFO and continue
186 *
187 IF( INFO.EQ.0 )
188 $ INFO = K
189 KP = K
190 A( K, K ) = DBLE( A( K, K ) )
191 ELSE
192 IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
193 *
194 * no interchange, use 1-by-1 pivot block
195 *
196 KP = K
197 ELSE
198 *
199 * Copy column IMAX to column KW-1 of W and update it
200 *
201 CALL ZCOPY( IMAX-1, A( 1, IMAX ), 1, W( 1, KW-1 ), 1 )
202 W( IMAX, KW-1 ) = DBLE( A( IMAX, IMAX ) )
203 CALL ZCOPY( K-IMAX, A( IMAX, IMAX+1 ), LDA,
204 $ W( IMAX+1, KW-1 ), 1 )
205 CALL ZLACGV( K-IMAX, W( IMAX+1, KW-1 ), 1 )
206 IF( K.LT.N ) THEN
207 CALL ZGEMV( 'No transpose', K, N-K, -CONE,
208 $ A( 1, K+1 ), LDA, W( IMAX, KW+1 ), LDW,
209 $ CONE, W( 1, KW-1 ), 1 )
210 W( IMAX, KW-1 ) = DBLE( W( IMAX, KW-1 ) )
211 END IF
212 *
213 * JMAX is the column-index of the largest off-diagonal
214 * element in row IMAX, and ROWMAX is its absolute value
215 *
216 JMAX = IMAX + IZAMAX( K-IMAX, W( IMAX+1, KW-1 ), 1 )
217 ROWMAX = CABS1( W( JMAX, KW-1 ) )
218 IF( IMAX.GT.1 ) THEN
219 JMAX = IZAMAX( IMAX-1, W( 1, KW-1 ), 1 )
220 ROWMAX = MAX( ROWMAX, CABS1( W( JMAX, KW-1 ) ) )
221 END IF
222 *
223 IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
224 *
225 * no interchange, use 1-by-1 pivot block
226 *
227 KP = K
228 ELSE IF( ABS( DBLE( W( IMAX, KW-1 ) ) ).GE.ALPHA*ROWMAX )
229 $ THEN
230 *
231 * interchange rows and columns K and IMAX, use 1-by-1
232 * pivot block
233 *
234 KP = IMAX
235 *
236 * copy column KW-1 of W to column KW
237 *
238 CALL ZCOPY( K, W( 1, KW-1 ), 1, W( 1, KW ), 1 )
239 ELSE
240 *
241 * interchange rows and columns K-1 and IMAX, use 2-by-2
242 * pivot block
243 *
244 KP = IMAX
245 KSTEP = 2
246 END IF
247 END IF
248 *
249 KK = K - KSTEP + 1
250 KKW = NB + KK - N
251 *
252 * Updated column KP is already stored in column KKW of W
253 *
254 IF( KP.NE.KK ) THEN
255 *
256 * Copy non-updated column KK to column KP
257 *
258 A( KP, KP ) = DBLE( A( KK, KK ) )
259 CALL ZCOPY( KK-1-KP, A( KP+1, KK ), 1, A( KP, KP+1 ),
260 $ LDA )
261 CALL ZLACGV( KK-1-KP, A( KP, KP+1 ), LDA )
262 CALL ZCOPY( KP-1, A( 1, KK ), 1, A( 1, KP ), 1 )
263 *
264 * Interchange rows KK and KP in last KK columns of A and W
265 *
266 IF( KK.LT.N )
267 $ CALL ZSWAP( N-KK, A( KK, KK+1 ), LDA, A( KP, KK+1 ),
268 $ LDA )
269 CALL ZSWAP( N-KK+1, W( KK, KKW ), LDW, W( KP, KKW ),
270 $ LDW )
271 END IF
272 *
273 IF( KSTEP.EQ.1 ) THEN
274 *
275 * 1-by-1 pivot block D(k): column KW of W now holds
276 *
277 * W(k) = U(k)*D(k)
278 *
279 * where U(k) is the k-th column of U
280 *
281 * Store U(k) in column k of A
282 *
283 CALL ZCOPY( K, W( 1, KW ), 1, A( 1, K ), 1 )
284 R1 = ONE / DBLE( A( K, K ) )
285 CALL ZDSCAL( K-1, R1, A( 1, K ), 1 )
286 *
287 * Conjugate W(k)
288 *
289 CALL ZLACGV( K-1, W( 1, KW ), 1 )
290 ELSE
291 *
292 * 2-by-2 pivot block D(k): columns KW and KW-1 of W now
293 * hold
294 *
295 * ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
296 *
297 * where U(k) and U(k-1) are the k-th and (k-1)-th columns
298 * of U
299 *
300 IF( K.GT.2 ) THEN
301 *
302 * Store U(k) and U(k-1) in columns k and k-1 of A
303 *
304 D21 = W( K-1, KW )
305 D11 = W( K, KW ) / DCONJG( D21 )
306 D22 = W( K-1, KW-1 ) / D21
307 T = ONE / ( DBLE( D11*D22 )-ONE )
308 D21 = T / D21
309 DO 20 J = 1, K - 2
310 A( J, K-1 ) = D21*( D11*W( J, KW-1 )-W( J, KW ) )
311 A( J, K ) = DCONJG( D21 )*
312 $ ( D22*W( J, KW )-W( J, KW-1 ) )
313 20 CONTINUE
314 END IF
315 *
316 * Copy D(k) to A
317 *
318 A( K-1, K-1 ) = W( K-1, KW-1 )
319 A( K-1, K ) = W( K-1, KW )
320 A( K, K ) = W( K, KW )
321 *
322 * Conjugate W(k) and W(k-1)
323 *
324 CALL ZLACGV( K-1, W( 1, KW ), 1 )
325 CALL ZLACGV( K-2, W( 1, KW-1 ), 1 )
326 END IF
327 END IF
328 *
329 * Store details of the interchanges in IPIV
330 *
331 IF( KSTEP.EQ.1 ) THEN
332 IPIV( K ) = KP
333 ELSE
334 IPIV( K ) = -KP
335 IPIV( K-1 ) = -KP
336 END IF
337 *
338 * Decrease K and return to the start of the main loop
339 *
340 K = K - KSTEP
341 GO TO 10
342 *
343 30 CONTINUE
344 *
345 * Update the upper triangle of A11 (= A(1:k,1:k)) as
346 *
347 * A11 := A11 - U12*D*U12**H = A11 - U12*W**H
348 *
349 * computing blocks of NB columns at a time (note that conjg(W) is
350 * actually stored)
351 *
352 DO 50 J = ( ( K-1 ) / NB )*NB + 1, 1, -NB
353 JB = MIN( NB, K-J+1 )
354 *
355 * Update the upper triangle of the diagonal block
356 *
357 DO 40 JJ = J, J + JB - 1
358 A( JJ, JJ ) = DBLE( A( JJ, JJ ) )
359 CALL ZGEMV( 'No transpose', JJ-J+1, N-K, -CONE,
360 $ A( J, K+1 ), LDA, W( JJ, KW+1 ), LDW, CONE,
361 $ A( J, JJ ), 1 )
362 A( JJ, JJ ) = DBLE( A( JJ, JJ ) )
363 40 CONTINUE
364 *
365 * Update the rectangular superdiagonal block
366 *
367 CALL ZGEMM( 'No transpose', 'Transpose', J-1, JB, N-K,
368 $ -CONE, A( 1, K+1 ), LDA, W( J, KW+1 ), LDW,
369 $ CONE, A( 1, J ), LDA )
370 50 CONTINUE
371 *
372 * Put U12 in standard form by partially undoing the interchanges
373 * in columns k+1:n
374 *
375 J = K + 1
376 60 CONTINUE
377 JJ = J
378 JP = IPIV( J )
379 IF( JP.LT.0 ) THEN
380 JP = -JP
381 J = J + 1
382 END IF
383 J = J + 1
384 IF( JP.NE.JJ .AND. J.LE.N )
385 $ CALL ZSWAP( N-J+1, A( JP, J ), LDA, A( JJ, J ), LDA )
386 IF( J.LE.N )
387 $ GO TO 60
388 *
389 * Set KB to the number of columns factorized
390 *
391 KB = N - K
392 *
393 ELSE
394 *
395 * Factorize the leading columns of A using the lower triangle
396 * of A and working forwards, and compute the matrix W = L21*D
397 * for use in updating A22 (note that conjg(W) is actually stored)
398 *
399 * K is the main loop index, increasing from 1 in steps of 1 or 2
400 *
401 K = 1
402 70 CONTINUE
403 *
404 * Exit from loop
405 *
406 IF( ( K.GE.NB .AND. NB.LT.N ) .OR. K.GT.N )
407 $ GO TO 90
408 *
409 * Copy column K of A to column K of W and update it
410 *
411 W( K, K ) = DBLE( A( K, K ) )
412 IF( K.LT.N )
413 $ CALL ZCOPY( N-K, A( K+1, K ), 1, W( K+1, K ), 1 )
414 CALL ZGEMV( 'No transpose', N-K+1, K-1, -CONE, A( K, 1 ), LDA,
415 $ W( K, 1 ), LDW, CONE, W( K, K ), 1 )
416 W( K, K ) = DBLE( W( K, K ) )
417 *
418 KSTEP = 1
419 *
420 * Determine rows and columns to be interchanged and whether
421 * a 1-by-1 or 2-by-2 pivot block will be used
422 *
423 ABSAKK = ABS( DBLE( W( K, K ) ) )
424 *
425 * IMAX is the row-index of the largest off-diagonal element in
426 * column K, and COLMAX is its absolute value
427 *
428 IF( K.LT.N ) THEN
429 IMAX = K + IZAMAX( N-K, W( K+1, K ), 1 )
430 COLMAX = CABS1( W( IMAX, K ) )
431 ELSE
432 COLMAX = ZERO
433 END IF
434 *
435 IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
436 *
437 * Column K is zero: set INFO and continue
438 *
439 IF( INFO.EQ.0 )
440 $ INFO = K
441 KP = K
442 A( K, K ) = DBLE( A( K, K ) )
443 ELSE
444 IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
445 *
446 * no interchange, use 1-by-1 pivot block
447 *
448 KP = K
449 ELSE
450 *
451 * Copy column IMAX to column K+1 of W and update it
452 *
453 CALL ZCOPY( IMAX-K, A( IMAX, K ), LDA, W( K, K+1 ), 1 )
454 CALL ZLACGV( IMAX-K, W( K, K+1 ), 1 )
455 W( IMAX, K+1 ) = DBLE( A( IMAX, IMAX ) )
456 IF( IMAX.LT.N )
457 $ CALL ZCOPY( N-IMAX, A( IMAX+1, IMAX ), 1,
458 $ W( IMAX+1, K+1 ), 1 )
459 CALL ZGEMV( 'No transpose', N-K+1, K-1, -CONE, A( K, 1 ),
460 $ LDA, W( IMAX, 1 ), LDW, CONE, W( K, K+1 ),
461 $ 1 )
462 W( IMAX, K+1 ) = DBLE( W( IMAX, K+1 ) )
463 *
464 * JMAX is the column-index of the largest off-diagonal
465 * element in row IMAX, and ROWMAX is its absolute value
466 *
467 JMAX = K - 1 + IZAMAX( IMAX-K, W( K, K+1 ), 1 )
468 ROWMAX = CABS1( W( JMAX, K+1 ) )
469 IF( IMAX.LT.N ) THEN
470 JMAX = IMAX + IZAMAX( N-IMAX, W( IMAX+1, K+1 ), 1 )
471 ROWMAX = MAX( ROWMAX, CABS1( W( JMAX, K+1 ) ) )
472 END IF
473 *
474 IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
475 *
476 * no interchange, use 1-by-1 pivot block
477 *
478 KP = K
479 ELSE IF( ABS( DBLE( W( IMAX, K+1 ) ) ).GE.ALPHA*ROWMAX )
480 $ THEN
481 *
482 * interchange rows and columns K and IMAX, use 1-by-1
483 * pivot block
484 *
485 KP = IMAX
486 *
487 * copy column K+1 of W to column K
488 *
489 CALL ZCOPY( N-K+1, W( K, K+1 ), 1, W( K, K ), 1 )
490 ELSE
491 *
492 * interchange rows and columns K+1 and IMAX, use 2-by-2
493 * pivot block
494 *
495 KP = IMAX
496 KSTEP = 2
497 END IF
498 END IF
499 *
500 KK = K + KSTEP - 1
501 *
502 * Updated column KP is already stored in column KK of W
503 *
504 IF( KP.NE.KK ) THEN
505 *
506 * Copy non-updated column KK to column KP
507 *
508 A( KP, KP ) = DBLE( A( KK, KK ) )
509 CALL ZCOPY( KP-KK-1, A( KK+1, KK ), 1, A( KP, KK+1 ),
510 $ LDA )
511 CALL ZLACGV( KP-KK-1, A( KP, KK+1 ), LDA )
512 IF( KP.LT.N )
513 $ CALL ZCOPY( N-KP, A( KP+1, KK ), 1, A( KP+1, KP ), 1 )
514 *
515 * Interchange rows KK and KP in first KK columns of A and W
516 *
517 CALL ZSWAP( KK-1, A( KK, 1 ), LDA, A( KP, 1 ), LDA )
518 CALL ZSWAP( KK, W( KK, 1 ), LDW, W( KP, 1 ), LDW )
519 END IF
520 *
521 IF( KSTEP.EQ.1 ) THEN
522 *
523 * 1-by-1 pivot block D(k): column k of W now holds
524 *
525 * W(k) = L(k)*D(k)
526 *
527 * where L(k) is the k-th column of L
528 *
529 * Store L(k) in column k of A
530 *
531 CALL ZCOPY( N-K+1, W( K, K ), 1, A( K, K ), 1 )
532 IF( K.LT.N ) THEN
533 R1 = ONE / DBLE( A( K, K ) )
534 CALL ZDSCAL( N-K, R1, A( K+1, K ), 1 )
535 *
536 * Conjugate W(k)
537 *
538 CALL ZLACGV( N-K, W( K+1, K ), 1 )
539 END IF
540 ELSE
541 *
542 * 2-by-2 pivot block D(k): columns k and k+1 of W now hold
543 *
544 * ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k)
545 *
546 * where L(k) and L(k+1) are the k-th and (k+1)-th columns
547 * of L
548 *
549 IF( K.LT.N-1 ) THEN
550 *
551 * Store L(k) and L(k+1) in columns k and k+1 of A
552 *
553 D21 = W( K+1, K )
554 D11 = W( K+1, K+1 ) / D21
555 D22 = W( K, K ) / DCONJG( D21 )
556 T = ONE / ( DBLE( D11*D22 )-ONE )
557 D21 = T / D21
558 DO 80 J = K + 2, N
559 A( J, K ) = DCONJG( D21 )*
560 $ ( D11*W( J, K )-W( J, K+1 ) )
561 A( J, K+1 ) = D21*( D22*W( J, K+1 )-W( J, K ) )
562 80 CONTINUE
563 END IF
564 *
565 * Copy D(k) to A
566 *
567 A( K, K ) = W( K, K )
568 A( K+1, K ) = W( K+1, K )
569 A( K+1, K+1 ) = W( K+1, K+1 )
570 *
571 * Conjugate W(k) and W(k+1)
572 *
573 CALL ZLACGV( N-K, W( K+1, K ), 1 )
574 CALL ZLACGV( N-K-1, W( K+2, K+1 ), 1 )
575 END IF
576 END IF
577 *
578 * Store details of the interchanges in IPIV
579 *
580 IF( KSTEP.EQ.1 ) THEN
581 IPIV( K ) = KP
582 ELSE
583 IPIV( K ) = -KP
584 IPIV( K+1 ) = -KP
585 END IF
586 *
587 * Increase K and return to the start of the main loop
588 *
589 K = K + KSTEP
590 GO TO 70
591 *
592 90 CONTINUE
593 *
594 * Update the lower triangle of A22 (= A(k:n,k:n)) as
595 *
596 * A22 := A22 - L21*D*L21**H = A22 - L21*W**H
597 *
598 * computing blocks of NB columns at a time (note that conjg(W) is
599 * actually stored)
600 *
601 DO 110 J = K, N, NB
602 JB = MIN( NB, N-J+1 )
603 *
604 * Update the lower triangle of the diagonal block
605 *
606 DO 100 JJ = J, J + JB - 1
607 A( JJ, JJ ) = DBLE( A( JJ, JJ ) )
608 CALL ZGEMV( 'No transpose', J+JB-JJ, K-1, -CONE,
609 $ A( JJ, 1 ), LDA, W( JJ, 1 ), LDW, CONE,
610 $ A( JJ, JJ ), 1 )
611 A( JJ, JJ ) = DBLE( A( JJ, JJ ) )
612 100 CONTINUE
613 *
614 * Update the rectangular subdiagonal block
615 *
616 IF( J+JB.LE.N )
617 $ CALL ZGEMM( 'No transpose', 'Transpose', N-J-JB+1, JB,
618 $ K-1, -CONE, A( J+JB, 1 ), LDA, W( J, 1 ),
619 $ LDW, CONE, A( J+JB, J ), LDA )
620 110 CONTINUE
621 *
622 * Put L21 in standard form by partially undoing the interchanges
623 * in columns 1:k-1
624 *
625 J = K - 1
626 120 CONTINUE
627 JJ = J
628 JP = IPIV( J )
629 IF( JP.LT.0 ) THEN
630 JP = -JP
631 J = J - 1
632 END IF
633 J = J - 1
634 IF( JP.NE.JJ .AND. J.GE.1 )
635 $ CALL ZSWAP( J, A( JP, 1 ), LDA, A( JJ, 1 ), LDA )
636 IF( J.GE.1 )
637 $ GO TO 120
638 *
639 * Set KB to the number of columns factorized
640 *
641 KB = K - 1
642 *
643 END IF
644 RETURN
645 *
646 * End of ZLAHEF
647 *
648 END