1 SUBROUTINE ZLASYF( 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 * ZLASYF computes a partial factorization of a complex symmetric matrix
21 * A using the Bunch-Kaufman diagonal pivoting method. The partial
22 * factorization has the form:
23 *
24 * A = ( I U12 ) ( A11 0 ) ( I 0 ) if UPLO = 'U', or:
25 * ( 0 U22 ) ( 0 D ) ( U12**T U22**T )
26 *
27 * A = ( L11 0 ) ( D 0 ) ( L11**T L21**T ) 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**T denotes the transpose of U.
33 *
34 * ZLASYF is an auxiliary routine called by ZSYTRF. 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 * symmetric 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 symmetric 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 DOUBLE PRECISION EIGHT, SEVTEN
102 PARAMETER ( EIGHT = 8.0D+0, SEVTEN = 17.0D+0 )
103 COMPLEX*16 CONE
104 PARAMETER ( CONE = ( 1.0D+0, 0.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, ROWMAX
110 COMPLEX*16 D11, D21, D22, R1, T, Z
111 * ..
112 * .. External Functions ..
113 LOGICAL LSAME
114 INTEGER IZAMAX
115 EXTERNAL LSAME, IZAMAX
116 * ..
117 * .. External Subroutines ..
118 EXTERNAL ZCOPY, ZGEMM, ZGEMV, ZSCAL, ZSWAP
119 * ..
120 * .. Intrinsic Functions ..
121 INTRINSIC ABS, DBLE, 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
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, A( 1, K ), 1, W( 1, KW ), 1 )
159 IF( K.LT.N )
160 $ CALL ZGEMV( 'No transpose', K, N-K, -CONE, A( 1, K+1 ), LDA,
161 $ W( K, KW+1 ), LDW, CONE, W( 1, KW ), 1 )
162 *
163 KSTEP = 1
164 *
165 * Determine rows and columns to be interchanged and whether
166 * a 1-by-1 or 2-by-2 pivot block will be used
167 *
168 ABSAKK = CABS1( W( K, KW ) )
169 *
170 * IMAX is the row-index of the largest off-diagonal element in
171 * column K, and COLMAX is its absolute value
172 *
173 IF( K.GT.1 ) THEN
174 IMAX = IZAMAX( K-1, W( 1, KW ), 1 )
175 COLMAX = CABS1( W( IMAX, KW ) )
176 ELSE
177 COLMAX = ZERO
178 END IF
179 *
180 IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
181 *
182 * Column K is zero: set INFO and continue
183 *
184 IF( INFO.EQ.0 )
185 $ INFO = K
186 KP = K
187 ELSE
188 IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
189 *
190 * no interchange, use 1-by-1 pivot block
191 *
192 KP = K
193 ELSE
194 *
195 * Copy column IMAX to column KW-1 of W and update it
196 *
197 CALL ZCOPY( IMAX, A( 1, IMAX ), 1, W( 1, KW-1 ), 1 )
198 CALL ZCOPY( K-IMAX, A( IMAX, IMAX+1 ), LDA,
199 $ W( IMAX+1, KW-1 ), 1 )
200 IF( K.LT.N )
201 $ CALL ZGEMV( 'No transpose', K, N-K, -CONE,
202 $ A( 1, K+1 ), LDA, W( IMAX, KW+1 ), LDW,
203 $ CONE, W( 1, KW-1 ), 1 )
204 *
205 * JMAX is the column-index of the largest off-diagonal
206 * element in row IMAX, and ROWMAX is its absolute value
207 *
208 JMAX = IMAX + IZAMAX( K-IMAX, W( IMAX+1, KW-1 ), 1 )
209 ROWMAX = CABS1( W( JMAX, KW-1 ) )
210 IF( IMAX.GT.1 ) THEN
211 JMAX = IZAMAX( IMAX-1, W( 1, KW-1 ), 1 )
212 ROWMAX = MAX( ROWMAX, CABS1( W( JMAX, KW-1 ) ) )
213 END IF
214 *
215 IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
216 *
217 * no interchange, use 1-by-1 pivot block
218 *
219 KP = K
220 ELSE IF( CABS1( W( IMAX, KW-1 ) ).GE.ALPHA*ROWMAX ) THEN
221 *
222 * interchange rows and columns K and IMAX, use 1-by-1
223 * pivot block
224 *
225 KP = IMAX
226 *
227 * copy column KW-1 of W to column KW
228 *
229 CALL ZCOPY( K, W( 1, KW-1 ), 1, W( 1, KW ), 1 )
230 ELSE
231 *
232 * interchange rows and columns K-1 and IMAX, use 2-by-2
233 * pivot block
234 *
235 KP = IMAX
236 KSTEP = 2
237 END IF
238 END IF
239 *
240 KK = K - KSTEP + 1
241 KKW = NB + KK - N
242 *
243 * Updated column KP is already stored in column KKW of W
244 *
245 IF( KP.NE.KK ) THEN
246 *
247 * Copy non-updated column KK to column KP
248 *
249 A( KP, K ) = A( KK, K )
250 CALL ZCOPY( K-1-KP, A( KP+1, KK ), 1, A( KP, KP+1 ),
251 $ LDA )
252 CALL ZCOPY( KP, A( 1, KK ), 1, A( 1, KP ), 1 )
253 *
254 * Interchange rows KK and KP in last KK columns of A and W
255 *
256 CALL ZSWAP( N-KK+1, A( KK, KK ), LDA, A( KP, KK ), LDA )
257 CALL ZSWAP( N-KK+1, W( KK, KKW ), LDW, W( KP, KKW ),
258 $ LDW )
259 END IF
260 *
261 IF( KSTEP.EQ.1 ) THEN
262 *
263 * 1-by-1 pivot block D(k): column KW of W now holds
264 *
265 * W(k) = U(k)*D(k)
266 *
267 * where U(k) is the k-th column of U
268 *
269 * Store U(k) in column k of A
270 *
271 CALL ZCOPY( K, W( 1, KW ), 1, A( 1, K ), 1 )
272 R1 = CONE / A( K, K )
273 CALL ZSCAL( K-1, R1, A( 1, K ), 1 )
274 ELSE
275 *
276 * 2-by-2 pivot block D(k): columns KW and KW-1 of W now
277 * hold
278 *
279 * ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
280 *
281 * where U(k) and U(k-1) are the k-th and (k-1)-th columns
282 * of U
283 *
284 IF( K.GT.2 ) THEN
285 *
286 * Store U(k) and U(k-1) in columns k and k-1 of A
287 *
288 D21 = W( K-1, KW )
289 D11 = W( K, KW ) / D21
290 D22 = W( K-1, KW-1 ) / D21
291 T = CONE / ( D11*D22-CONE )
292 D21 = T / D21
293 DO 20 J = 1, K - 2
294 A( J, K-1 ) = D21*( D11*W( J, KW-1 )-W( J, KW ) )
295 A( J, K ) = D21*( D22*W( J, KW )-W( J, KW-1 ) )
296 20 CONTINUE
297 END IF
298 *
299 * Copy D(k) to A
300 *
301 A( K-1, K-1 ) = W( K-1, KW-1 )
302 A( K-1, K ) = W( K-1, KW )
303 A( K, K ) = W( K, KW )
304 END IF
305 END IF
306 *
307 * Store details of the interchanges in IPIV
308 *
309 IF( KSTEP.EQ.1 ) THEN
310 IPIV( K ) = KP
311 ELSE
312 IPIV( K ) = -KP
313 IPIV( K-1 ) = -KP
314 END IF
315 *
316 * Decrease K and return to the start of the main loop
317 *
318 K = K - KSTEP
319 GO TO 10
320 *
321 30 CONTINUE
322 *
323 * Update the upper triangle of A11 (= A(1:k,1:k)) as
324 *
325 * A11 := A11 - U12*D*U12**T = A11 - U12*W**T
326 *
327 * computing blocks of NB columns at a time
328 *
329 DO 50 J = ( ( K-1 ) / NB )*NB + 1, 1, -NB
330 JB = MIN( NB, K-J+1 )
331 *
332 * Update the upper triangle of the diagonal block
333 *
334 DO 40 JJ = J, J + JB - 1
335 CALL ZGEMV( 'No transpose', JJ-J+1, N-K, -CONE,
336 $ A( J, K+1 ), LDA, W( JJ, KW+1 ), LDW, CONE,
337 $ A( J, JJ ), 1 )
338 40 CONTINUE
339 *
340 * Update the rectangular superdiagonal block
341 *
342 CALL ZGEMM( 'No transpose', 'Transpose', J-1, JB, N-K,
343 $ -CONE, A( 1, K+1 ), LDA, W( J, KW+1 ), LDW,
344 $ CONE, A( 1, J ), LDA )
345 50 CONTINUE
346 *
347 * Put U12 in standard form by partially undoing the interchanges
348 * in columns k+1:n
349 *
350 J = K + 1
351 60 CONTINUE
352 JJ = J
353 JP = IPIV( J )
354 IF( JP.LT.0 ) THEN
355 JP = -JP
356 J = J + 1
357 END IF
358 J = J + 1
359 IF( JP.NE.JJ .AND. J.LE.N )
360 $ CALL ZSWAP( N-J+1, A( JP, J ), LDA, A( JJ, J ), LDA )
361 IF( J.LE.N )
362 $ GO TO 60
363 *
364 * Set KB to the number of columns factorized
365 *
366 KB = N - K
367 *
368 ELSE
369 *
370 * Factorize the leading columns of A using the lower triangle
371 * of A and working forwards, and compute the matrix W = L21*D
372 * for use in updating A22
373 *
374 * K is the main loop index, increasing from 1 in steps of 1 or 2
375 *
376 K = 1
377 70 CONTINUE
378 *
379 * Exit from loop
380 *
381 IF( ( K.GE.NB .AND. NB.LT.N ) .OR. K.GT.N )
382 $ GO TO 90
383 *
384 * Copy column K of A to column K of W and update it
385 *
386 CALL ZCOPY( N-K+1, A( K, K ), 1, W( K, K ), 1 )
387 CALL ZGEMV( 'No transpose', N-K+1, K-1, -CONE, A( K, 1 ), LDA,
388 $ W( K, 1 ), LDW, CONE, W( K, K ), 1 )
389 *
390 KSTEP = 1
391 *
392 * Determine rows and columns to be interchanged and whether
393 * a 1-by-1 or 2-by-2 pivot block will be used
394 *
395 ABSAKK = CABS1( W( K, K ) )
396 *
397 * IMAX is the row-index of the largest off-diagonal element in
398 * column K, and COLMAX is its absolute value
399 *
400 IF( K.LT.N ) THEN
401 IMAX = K + IZAMAX( N-K, W( K+1, K ), 1 )
402 COLMAX = CABS1( W( IMAX, K ) )
403 ELSE
404 COLMAX = ZERO
405 END IF
406 *
407 IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
408 *
409 * Column K is zero: set INFO and continue
410 *
411 IF( INFO.EQ.0 )
412 $ INFO = K
413 KP = K
414 ELSE
415 IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
416 *
417 * no interchange, use 1-by-1 pivot block
418 *
419 KP = K
420 ELSE
421 *
422 * Copy column IMAX to column K+1 of W and update it
423 *
424 CALL ZCOPY( IMAX-K, A( IMAX, K ), LDA, W( K, K+1 ), 1 )
425 CALL ZCOPY( N-IMAX+1, A( IMAX, IMAX ), 1, W( IMAX, K+1 ),
426 $ 1 )
427 CALL ZGEMV( 'No transpose', N-K+1, K-1, -CONE, A( K, 1 ),
428 $ LDA, W( IMAX, 1 ), LDW, CONE, W( K, K+1 ),
429 $ 1 )
430 *
431 * JMAX is the column-index of the largest off-diagonal
432 * element in row IMAX, and ROWMAX is its absolute value
433 *
434 JMAX = K - 1 + IZAMAX( IMAX-K, W( K, K+1 ), 1 )
435 ROWMAX = CABS1( W( JMAX, K+1 ) )
436 IF( IMAX.LT.N ) THEN
437 JMAX = IMAX + IZAMAX( N-IMAX, W( IMAX+1, K+1 ), 1 )
438 ROWMAX = MAX( ROWMAX, CABS1( W( JMAX, K+1 ) ) )
439 END IF
440 *
441 IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
442 *
443 * no interchange, use 1-by-1 pivot block
444 *
445 KP = K
446 ELSE IF( CABS1( W( IMAX, K+1 ) ).GE.ALPHA*ROWMAX ) THEN
447 *
448 * interchange rows and columns K and IMAX, use 1-by-1
449 * pivot block
450 *
451 KP = IMAX
452 *
453 * copy column K+1 of W to column K
454 *
455 CALL ZCOPY( N-K+1, W( K, K+1 ), 1, W( K, K ), 1 )
456 ELSE
457 *
458 * interchange rows and columns K+1 and IMAX, use 2-by-2
459 * pivot block
460 *
461 KP = IMAX
462 KSTEP = 2
463 END IF
464 END IF
465 *
466 KK = K + KSTEP - 1
467 *
468 * Updated column KP is already stored in column KK of W
469 *
470 IF( KP.NE.KK ) THEN
471 *
472 * Copy non-updated column KK to column KP
473 *
474 A( KP, K ) = A( KK, K )
475 CALL ZCOPY( KP-K-1, A( K+1, KK ), 1, A( KP, K+1 ), LDA )
476 CALL ZCOPY( N-KP+1, A( KP, KK ), 1, A( KP, KP ), 1 )
477 *
478 * Interchange rows KK and KP in first KK columns of A and W
479 *
480 CALL ZSWAP( KK, A( KK, 1 ), LDA, A( KP, 1 ), LDA )
481 CALL ZSWAP( KK, W( KK, 1 ), LDW, W( KP, 1 ), LDW )
482 END IF
483 *
484 IF( KSTEP.EQ.1 ) THEN
485 *
486 * 1-by-1 pivot block D(k): column k of W now holds
487 *
488 * W(k) = L(k)*D(k)
489 *
490 * where L(k) is the k-th column of L
491 *
492 * Store L(k) in column k of A
493 *
494 CALL ZCOPY( N-K+1, W( K, K ), 1, A( K, K ), 1 )
495 IF( K.LT.N ) THEN
496 R1 = CONE / A( K, K )
497 CALL ZSCAL( N-K, R1, A( K+1, K ), 1 )
498 END IF
499 ELSE
500 *
501 * 2-by-2 pivot block D(k): columns k and k+1 of W now hold
502 *
503 * ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k)
504 *
505 * where L(k) and L(k+1) are the k-th and (k+1)-th columns
506 * of L
507 *
508 IF( K.LT.N-1 ) THEN
509 *
510 * Store L(k) and L(k+1) in columns k and k+1 of A
511 *
512 D21 = W( K+1, K )
513 D11 = W( K+1, K+1 ) / D21
514 D22 = W( K, K ) / D21
515 T = CONE / ( D11*D22-CONE )
516 D21 = T / D21
517 DO 80 J = K + 2, N
518 A( J, K ) = D21*( D11*W( J, K )-W( J, K+1 ) )
519 A( J, K+1 ) = D21*( D22*W( J, K+1 )-W( J, K ) )
520 80 CONTINUE
521 END IF
522 *
523 * Copy D(k) to A
524 *
525 A( K, K ) = W( K, K )
526 A( K+1, K ) = W( K+1, K )
527 A( K+1, K+1 ) = W( K+1, K+1 )
528 END IF
529 END IF
530 *
531 * Store details of the interchanges in IPIV
532 *
533 IF( KSTEP.EQ.1 ) THEN
534 IPIV( K ) = KP
535 ELSE
536 IPIV( K ) = -KP
537 IPIV( K+1 ) = -KP
538 END IF
539 *
540 * Increase K and return to the start of the main loop
541 *
542 K = K + KSTEP
543 GO TO 70
544 *
545 90 CONTINUE
546 *
547 * Update the lower triangle of A22 (= A(k:n,k:n)) as
548 *
549 * A22 := A22 - L21*D*L21**T = A22 - L21*W**T
550 *
551 * computing blocks of NB columns at a time
552 *
553 DO 110 J = K, N, NB
554 JB = MIN( NB, N-J+1 )
555 *
556 * Update the lower triangle of the diagonal block
557 *
558 DO 100 JJ = J, J + JB - 1
559 CALL ZGEMV( 'No transpose', J+JB-JJ, K-1, -CONE,
560 $ A( JJ, 1 ), LDA, W( JJ, 1 ), LDW, CONE,
561 $ A( JJ, JJ ), 1 )
562 100 CONTINUE
563 *
564 * Update the rectangular subdiagonal block
565 *
566 IF( J+JB.LE.N )
567 $ CALL ZGEMM( 'No transpose', 'Transpose', N-J-JB+1, JB,
568 $ K-1, -CONE, A( J+JB, 1 ), LDA, W( J, 1 ),
569 $ LDW, CONE, A( J+JB, J ), LDA )
570 110 CONTINUE
571 *
572 * Put L21 in standard form by partially undoing the interchanges
573 * in columns 1:k-1
574 *
575 J = K - 1
576 120 CONTINUE
577 JJ = J
578 JP = IPIV( J )
579 IF( JP.LT.0 ) THEN
580 JP = -JP
581 J = J - 1
582 END IF
583 J = J - 1
584 IF( JP.NE.JJ .AND. J.GE.1 )
585 $ CALL ZSWAP( J, A( JP, 1 ), LDA, A( JJ, 1 ), LDA )
586 IF( J.GE.1 )
587 $ GO TO 120
588 *
589 * Set KB to the number of columns factorized
590 *
591 KB = K - 1
592 *
593 END IF
594 RETURN
595 *
596 * End of ZLASYF
597 *
598 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 * ZLASYF computes a partial factorization of a complex symmetric matrix
21 * A using the Bunch-Kaufman diagonal pivoting method. The partial
22 * factorization has the form:
23 *
24 * A = ( I U12 ) ( A11 0 ) ( I 0 ) if UPLO = 'U', or:
25 * ( 0 U22 ) ( 0 D ) ( U12**T U22**T )
26 *
27 * A = ( L11 0 ) ( D 0 ) ( L11**T L21**T ) 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**T denotes the transpose of U.
33 *
34 * ZLASYF is an auxiliary routine called by ZSYTRF. 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 * symmetric 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 symmetric 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 DOUBLE PRECISION EIGHT, SEVTEN
102 PARAMETER ( EIGHT = 8.0D+0, SEVTEN = 17.0D+0 )
103 COMPLEX*16 CONE
104 PARAMETER ( CONE = ( 1.0D+0, 0.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, ROWMAX
110 COMPLEX*16 D11, D21, D22, R1, T, Z
111 * ..
112 * .. External Functions ..
113 LOGICAL LSAME
114 INTEGER IZAMAX
115 EXTERNAL LSAME, IZAMAX
116 * ..
117 * .. External Subroutines ..
118 EXTERNAL ZCOPY, ZGEMM, ZGEMV, ZSCAL, ZSWAP
119 * ..
120 * .. Intrinsic Functions ..
121 INTRINSIC ABS, DBLE, 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
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, A( 1, K ), 1, W( 1, KW ), 1 )
159 IF( K.LT.N )
160 $ CALL ZGEMV( 'No transpose', K, N-K, -CONE, A( 1, K+1 ), LDA,
161 $ W( K, KW+1 ), LDW, CONE, W( 1, KW ), 1 )
162 *
163 KSTEP = 1
164 *
165 * Determine rows and columns to be interchanged and whether
166 * a 1-by-1 or 2-by-2 pivot block will be used
167 *
168 ABSAKK = CABS1( W( K, KW ) )
169 *
170 * IMAX is the row-index of the largest off-diagonal element in
171 * column K, and COLMAX is its absolute value
172 *
173 IF( K.GT.1 ) THEN
174 IMAX = IZAMAX( K-1, W( 1, KW ), 1 )
175 COLMAX = CABS1( W( IMAX, KW ) )
176 ELSE
177 COLMAX = ZERO
178 END IF
179 *
180 IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
181 *
182 * Column K is zero: set INFO and continue
183 *
184 IF( INFO.EQ.0 )
185 $ INFO = K
186 KP = K
187 ELSE
188 IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
189 *
190 * no interchange, use 1-by-1 pivot block
191 *
192 KP = K
193 ELSE
194 *
195 * Copy column IMAX to column KW-1 of W and update it
196 *
197 CALL ZCOPY( IMAX, A( 1, IMAX ), 1, W( 1, KW-1 ), 1 )
198 CALL ZCOPY( K-IMAX, A( IMAX, IMAX+1 ), LDA,
199 $ W( IMAX+1, KW-1 ), 1 )
200 IF( K.LT.N )
201 $ CALL ZGEMV( 'No transpose', K, N-K, -CONE,
202 $ A( 1, K+1 ), LDA, W( IMAX, KW+1 ), LDW,
203 $ CONE, W( 1, KW-1 ), 1 )
204 *
205 * JMAX is the column-index of the largest off-diagonal
206 * element in row IMAX, and ROWMAX is its absolute value
207 *
208 JMAX = IMAX + IZAMAX( K-IMAX, W( IMAX+1, KW-1 ), 1 )
209 ROWMAX = CABS1( W( JMAX, KW-1 ) )
210 IF( IMAX.GT.1 ) THEN
211 JMAX = IZAMAX( IMAX-1, W( 1, KW-1 ), 1 )
212 ROWMAX = MAX( ROWMAX, CABS1( W( JMAX, KW-1 ) ) )
213 END IF
214 *
215 IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
216 *
217 * no interchange, use 1-by-1 pivot block
218 *
219 KP = K
220 ELSE IF( CABS1( W( IMAX, KW-1 ) ).GE.ALPHA*ROWMAX ) THEN
221 *
222 * interchange rows and columns K and IMAX, use 1-by-1
223 * pivot block
224 *
225 KP = IMAX
226 *
227 * copy column KW-1 of W to column KW
228 *
229 CALL ZCOPY( K, W( 1, KW-1 ), 1, W( 1, KW ), 1 )
230 ELSE
231 *
232 * interchange rows and columns K-1 and IMAX, use 2-by-2
233 * pivot block
234 *
235 KP = IMAX
236 KSTEP = 2
237 END IF
238 END IF
239 *
240 KK = K - KSTEP + 1
241 KKW = NB + KK - N
242 *
243 * Updated column KP is already stored in column KKW of W
244 *
245 IF( KP.NE.KK ) THEN
246 *
247 * Copy non-updated column KK to column KP
248 *
249 A( KP, K ) = A( KK, K )
250 CALL ZCOPY( K-1-KP, A( KP+1, KK ), 1, A( KP, KP+1 ),
251 $ LDA )
252 CALL ZCOPY( KP, A( 1, KK ), 1, A( 1, KP ), 1 )
253 *
254 * Interchange rows KK and KP in last KK columns of A and W
255 *
256 CALL ZSWAP( N-KK+1, A( KK, KK ), LDA, A( KP, KK ), LDA )
257 CALL ZSWAP( N-KK+1, W( KK, KKW ), LDW, W( KP, KKW ),
258 $ LDW )
259 END IF
260 *
261 IF( KSTEP.EQ.1 ) THEN
262 *
263 * 1-by-1 pivot block D(k): column KW of W now holds
264 *
265 * W(k) = U(k)*D(k)
266 *
267 * where U(k) is the k-th column of U
268 *
269 * Store U(k) in column k of A
270 *
271 CALL ZCOPY( K, W( 1, KW ), 1, A( 1, K ), 1 )
272 R1 = CONE / A( K, K )
273 CALL ZSCAL( K-1, R1, A( 1, K ), 1 )
274 ELSE
275 *
276 * 2-by-2 pivot block D(k): columns KW and KW-1 of W now
277 * hold
278 *
279 * ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
280 *
281 * where U(k) and U(k-1) are the k-th and (k-1)-th columns
282 * of U
283 *
284 IF( K.GT.2 ) THEN
285 *
286 * Store U(k) and U(k-1) in columns k and k-1 of A
287 *
288 D21 = W( K-1, KW )
289 D11 = W( K, KW ) / D21
290 D22 = W( K-1, KW-1 ) / D21
291 T = CONE / ( D11*D22-CONE )
292 D21 = T / D21
293 DO 20 J = 1, K - 2
294 A( J, K-1 ) = D21*( D11*W( J, KW-1 )-W( J, KW ) )
295 A( J, K ) = D21*( D22*W( J, KW )-W( J, KW-1 ) )
296 20 CONTINUE
297 END IF
298 *
299 * Copy D(k) to A
300 *
301 A( K-1, K-1 ) = W( K-1, KW-1 )
302 A( K-1, K ) = W( K-1, KW )
303 A( K, K ) = W( K, KW )
304 END IF
305 END IF
306 *
307 * Store details of the interchanges in IPIV
308 *
309 IF( KSTEP.EQ.1 ) THEN
310 IPIV( K ) = KP
311 ELSE
312 IPIV( K ) = -KP
313 IPIV( K-1 ) = -KP
314 END IF
315 *
316 * Decrease K and return to the start of the main loop
317 *
318 K = K - KSTEP
319 GO TO 10
320 *
321 30 CONTINUE
322 *
323 * Update the upper triangle of A11 (= A(1:k,1:k)) as
324 *
325 * A11 := A11 - U12*D*U12**T = A11 - U12*W**T
326 *
327 * computing blocks of NB columns at a time
328 *
329 DO 50 J = ( ( K-1 ) / NB )*NB + 1, 1, -NB
330 JB = MIN( NB, K-J+1 )
331 *
332 * Update the upper triangle of the diagonal block
333 *
334 DO 40 JJ = J, J + JB - 1
335 CALL ZGEMV( 'No transpose', JJ-J+1, N-K, -CONE,
336 $ A( J, K+1 ), LDA, W( JJ, KW+1 ), LDW, CONE,
337 $ A( J, JJ ), 1 )
338 40 CONTINUE
339 *
340 * Update the rectangular superdiagonal block
341 *
342 CALL ZGEMM( 'No transpose', 'Transpose', J-1, JB, N-K,
343 $ -CONE, A( 1, K+1 ), LDA, W( J, KW+1 ), LDW,
344 $ CONE, A( 1, J ), LDA )
345 50 CONTINUE
346 *
347 * Put U12 in standard form by partially undoing the interchanges
348 * in columns k+1:n
349 *
350 J = K + 1
351 60 CONTINUE
352 JJ = J
353 JP = IPIV( J )
354 IF( JP.LT.0 ) THEN
355 JP = -JP
356 J = J + 1
357 END IF
358 J = J + 1
359 IF( JP.NE.JJ .AND. J.LE.N )
360 $ CALL ZSWAP( N-J+1, A( JP, J ), LDA, A( JJ, J ), LDA )
361 IF( J.LE.N )
362 $ GO TO 60
363 *
364 * Set KB to the number of columns factorized
365 *
366 KB = N - K
367 *
368 ELSE
369 *
370 * Factorize the leading columns of A using the lower triangle
371 * of A and working forwards, and compute the matrix W = L21*D
372 * for use in updating A22
373 *
374 * K is the main loop index, increasing from 1 in steps of 1 or 2
375 *
376 K = 1
377 70 CONTINUE
378 *
379 * Exit from loop
380 *
381 IF( ( K.GE.NB .AND. NB.LT.N ) .OR. K.GT.N )
382 $ GO TO 90
383 *
384 * Copy column K of A to column K of W and update it
385 *
386 CALL ZCOPY( N-K+1, A( K, K ), 1, W( K, K ), 1 )
387 CALL ZGEMV( 'No transpose', N-K+1, K-1, -CONE, A( K, 1 ), LDA,
388 $ W( K, 1 ), LDW, CONE, W( K, K ), 1 )
389 *
390 KSTEP = 1
391 *
392 * Determine rows and columns to be interchanged and whether
393 * a 1-by-1 or 2-by-2 pivot block will be used
394 *
395 ABSAKK = CABS1( W( K, K ) )
396 *
397 * IMAX is the row-index of the largest off-diagonal element in
398 * column K, and COLMAX is its absolute value
399 *
400 IF( K.LT.N ) THEN
401 IMAX = K + IZAMAX( N-K, W( K+1, K ), 1 )
402 COLMAX = CABS1( W( IMAX, K ) )
403 ELSE
404 COLMAX = ZERO
405 END IF
406 *
407 IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
408 *
409 * Column K is zero: set INFO and continue
410 *
411 IF( INFO.EQ.0 )
412 $ INFO = K
413 KP = K
414 ELSE
415 IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
416 *
417 * no interchange, use 1-by-1 pivot block
418 *
419 KP = K
420 ELSE
421 *
422 * Copy column IMAX to column K+1 of W and update it
423 *
424 CALL ZCOPY( IMAX-K, A( IMAX, K ), LDA, W( K, K+1 ), 1 )
425 CALL ZCOPY( N-IMAX+1, A( IMAX, IMAX ), 1, W( IMAX, K+1 ),
426 $ 1 )
427 CALL ZGEMV( 'No transpose', N-K+1, K-1, -CONE, A( K, 1 ),
428 $ LDA, W( IMAX, 1 ), LDW, CONE, W( K, K+1 ),
429 $ 1 )
430 *
431 * JMAX is the column-index of the largest off-diagonal
432 * element in row IMAX, and ROWMAX is its absolute value
433 *
434 JMAX = K - 1 + IZAMAX( IMAX-K, W( K, K+1 ), 1 )
435 ROWMAX = CABS1( W( JMAX, K+1 ) )
436 IF( IMAX.LT.N ) THEN
437 JMAX = IMAX + IZAMAX( N-IMAX, W( IMAX+1, K+1 ), 1 )
438 ROWMAX = MAX( ROWMAX, CABS1( W( JMAX, K+1 ) ) )
439 END IF
440 *
441 IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
442 *
443 * no interchange, use 1-by-1 pivot block
444 *
445 KP = K
446 ELSE IF( CABS1( W( IMAX, K+1 ) ).GE.ALPHA*ROWMAX ) THEN
447 *
448 * interchange rows and columns K and IMAX, use 1-by-1
449 * pivot block
450 *
451 KP = IMAX
452 *
453 * copy column K+1 of W to column K
454 *
455 CALL ZCOPY( N-K+1, W( K, K+1 ), 1, W( K, K ), 1 )
456 ELSE
457 *
458 * interchange rows and columns K+1 and IMAX, use 2-by-2
459 * pivot block
460 *
461 KP = IMAX
462 KSTEP = 2
463 END IF
464 END IF
465 *
466 KK = K + KSTEP - 1
467 *
468 * Updated column KP is already stored in column KK of W
469 *
470 IF( KP.NE.KK ) THEN
471 *
472 * Copy non-updated column KK to column KP
473 *
474 A( KP, K ) = A( KK, K )
475 CALL ZCOPY( KP-K-1, A( K+1, KK ), 1, A( KP, K+1 ), LDA )
476 CALL ZCOPY( N-KP+1, A( KP, KK ), 1, A( KP, KP ), 1 )
477 *
478 * Interchange rows KK and KP in first KK columns of A and W
479 *
480 CALL ZSWAP( KK, A( KK, 1 ), LDA, A( KP, 1 ), LDA )
481 CALL ZSWAP( KK, W( KK, 1 ), LDW, W( KP, 1 ), LDW )
482 END IF
483 *
484 IF( KSTEP.EQ.1 ) THEN
485 *
486 * 1-by-1 pivot block D(k): column k of W now holds
487 *
488 * W(k) = L(k)*D(k)
489 *
490 * where L(k) is the k-th column of L
491 *
492 * Store L(k) in column k of A
493 *
494 CALL ZCOPY( N-K+1, W( K, K ), 1, A( K, K ), 1 )
495 IF( K.LT.N ) THEN
496 R1 = CONE / A( K, K )
497 CALL ZSCAL( N-K, R1, A( K+1, K ), 1 )
498 END IF
499 ELSE
500 *
501 * 2-by-2 pivot block D(k): columns k and k+1 of W now hold
502 *
503 * ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k)
504 *
505 * where L(k) and L(k+1) are the k-th and (k+1)-th columns
506 * of L
507 *
508 IF( K.LT.N-1 ) THEN
509 *
510 * Store L(k) and L(k+1) in columns k and k+1 of A
511 *
512 D21 = W( K+1, K )
513 D11 = W( K+1, K+1 ) / D21
514 D22 = W( K, K ) / D21
515 T = CONE / ( D11*D22-CONE )
516 D21 = T / D21
517 DO 80 J = K + 2, N
518 A( J, K ) = D21*( D11*W( J, K )-W( J, K+1 ) )
519 A( J, K+1 ) = D21*( D22*W( J, K+1 )-W( J, K ) )
520 80 CONTINUE
521 END IF
522 *
523 * Copy D(k) to A
524 *
525 A( K, K ) = W( K, K )
526 A( K+1, K ) = W( K+1, K )
527 A( K+1, K+1 ) = W( K+1, K+1 )
528 END IF
529 END IF
530 *
531 * Store details of the interchanges in IPIV
532 *
533 IF( KSTEP.EQ.1 ) THEN
534 IPIV( K ) = KP
535 ELSE
536 IPIV( K ) = -KP
537 IPIV( K+1 ) = -KP
538 END IF
539 *
540 * Increase K and return to the start of the main loop
541 *
542 K = K + KSTEP
543 GO TO 70
544 *
545 90 CONTINUE
546 *
547 * Update the lower triangle of A22 (= A(k:n,k:n)) as
548 *
549 * A22 := A22 - L21*D*L21**T = A22 - L21*W**T
550 *
551 * computing blocks of NB columns at a time
552 *
553 DO 110 J = K, N, NB
554 JB = MIN( NB, N-J+1 )
555 *
556 * Update the lower triangle of the diagonal block
557 *
558 DO 100 JJ = J, J + JB - 1
559 CALL ZGEMV( 'No transpose', J+JB-JJ, K-1, -CONE,
560 $ A( JJ, 1 ), LDA, W( JJ, 1 ), LDW, CONE,
561 $ A( JJ, JJ ), 1 )
562 100 CONTINUE
563 *
564 * Update the rectangular subdiagonal block
565 *
566 IF( J+JB.LE.N )
567 $ CALL ZGEMM( 'No transpose', 'Transpose', N-J-JB+1, JB,
568 $ K-1, -CONE, A( J+JB, 1 ), LDA, W( J, 1 ),
569 $ LDW, CONE, A( J+JB, J ), LDA )
570 110 CONTINUE
571 *
572 * Put L21 in standard form by partially undoing the interchanges
573 * in columns 1:k-1
574 *
575 J = K - 1
576 120 CONTINUE
577 JJ = J
578 JP = IPIV( J )
579 IF( JP.LT.0 ) THEN
580 JP = -JP
581 J = J - 1
582 END IF
583 J = J - 1
584 IF( JP.NE.JJ .AND. J.GE.1 )
585 $ CALL ZSWAP( J, A( JP, 1 ), LDA, A( JJ, 1 ), LDA )
586 IF( J.GE.1 )
587 $ GO TO 120
588 *
589 * Set KB to the number of columns factorized
590 *
591 KB = K - 1
592 *
593 END IF
594 RETURN
595 *
596 * End of ZLASYF
597 *
598 END