1 *> \brief \b CLASYF_ROOK computes a partial factorization of a complex symmetric matrix using the bounded Bunch-Kaufman ("rook") diagonal pivoting method.
3 * =========== DOCUMENTATION ===========
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
9 *> Download CLASYF_ROOK + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clasyf_rook.f">
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/clasyf_rook.f">
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/clasyf_rook.f">
21 * SUBROUTINE CLASYF_ROOK( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, INFO )
23 * .. Scalar Arguments ..
25 * INTEGER INFO, KB, LDA, LDW, N, NB
27 * .. Array Arguments ..
29 * COMPLEX A( LDA, * ), W( LDW, * )
38 *> CLASYF_ROOK computes a partial factorization of a complex symmetric
39 *> matrix A using the bounded Bunch-Kaufman ("rook") diagonal
40 *> pivoting method. The partial factorization has the form:
42 *> A = ( I U12 ) ( A11 0 ) ( I 0 ) if UPLO = 'U', or:
43 *> ( 0 U22 ) ( 0 D ) ( U12**T U22**T )
45 *> A = ( L11 0 ) ( D 0 ) ( L11**T L21**T ) if UPLO = 'L'
46 *> ( L21 I ) ( 0 A22 ) ( 0 I )
48 *> where the order of D is at most NB. The actual order is returned in
49 *> the argument KB, and is either NB or NB-1, or N if N <= NB.
51 *> CLASYF_ROOK is an auxiliary routine called by CSYTRF_ROOK. It uses
52 *> blocked code (calling Level 3 BLAS) to update the submatrix
53 *> A11 (if UPLO = 'U') or A22 (if UPLO = 'L').
61 *> UPLO is CHARACTER*1
62 *> Specifies whether the upper or lower triangular part of the
63 *> symmetric matrix A is stored:
64 *> = 'U': Upper triangular
65 *> = 'L': Lower triangular
71 *> The order of the matrix A. N >= 0.
77 *> The maximum number of columns of the matrix A that should be
78 *> factored. NB should be at least 2 to allow for 2-by-2 pivot
85 *> The number of columns of A that were actually factored.
86 *> KB is either NB-1 or NB, or N if N <= NB.
91 *> A is COMPLEX array, dimension (LDA,N)
92 *> On entry, the symmetric matrix A. If UPLO = 'U', the leading
93 *> n-by-n upper triangular part of A contains the upper
94 *> triangular part of the matrix A, and the strictly lower
95 *> triangular part of A is not referenced. If UPLO = 'L', the
96 *> leading n-by-n lower triangular part of A contains the lower
97 *> triangular part of the matrix A, and the strictly upper
98 *> triangular part of A is not referenced.
99 *> On exit, A contains details of the partial factorization.
105 *> The leading dimension of the array A. LDA >= max(1,N).
110 *> IPIV is INTEGER array, dimension (N)
111 *> Details of the interchanges and the block structure of D.
114 *> Only the last KB elements of IPIV are set.
116 *> If IPIV(k) > 0, then rows and columns k and IPIV(k) were
117 *> interchanged and D(k,k) is a 1-by-1 diagonal block.
119 *> If IPIV(k) < 0 and IPIV(k-1) < 0, then rows and
120 *> columns k and -IPIV(k) were interchanged and rows and
121 *> columns k-1 and -IPIV(k-1) were inerchaged,
122 *> D(k-1:k,k-1:k) is a 2-by-2 diagonal block.
125 *> Only the first KB elements of IPIV are set.
127 *> If IPIV(k) > 0, then rows and columns k and IPIV(k)
128 *> were interchanged and D(k,k) is a 1-by-1 diagonal block.
130 *> If IPIV(k) < 0 and IPIV(k+1) < 0, then rows and
131 *> columns k and -IPIV(k) were interchanged and rows and
132 *> columns k+1 and -IPIV(k+1) were inerchaged,
133 *> D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
138 *> W is COMPLEX array, dimension (LDW,NB)
144 *> The leading dimension of the array W. LDW >= max(1,N).
150 *> = 0: successful exit
151 *> > 0: if INFO = k, D(k,k) is exactly zero. The factorization
152 *> has been completed, but the block diagonal matrix D is
159 *> \author Univ. of Tennessee
160 *> \author Univ. of California Berkeley
161 *> \author Univ. of Colorado Denver
164 *> \date November 2013
166 *> \ingroup complexSYcomputational
168 *> \par Contributors:
173 *> November 2013, Igor Kozachenko,
174 *> Computer Science Division,
175 *> University of California, Berkeley
177 *> September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas,
178 *> School of Mathematics,
179 *> University of Manchester
183 * =====================================================================
184 SUBROUTINE CLASYF_ROOK( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW,
187 * -- LAPACK computational routine (version 3.5.0) --
188 * -- LAPACK is a software package provided by Univ. of Tennessee, --
189 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
192 * .. Scalar Arguments ..
194 INTEGER INFO, KB, LDA, LDW, N, NB
196 * .. Array Arguments ..
198 COMPLEX A( LDA, * ), W( LDW, * )
201 * =====================================================================
205 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 )
207 PARAMETER ( EIGHT = 8.0E+0, SEVTEN = 17.0E+0 )
209 PARAMETER ( CONE = ( 1.0E+0, 0.0E+0 ),
210 $ CZERO = ( 0.0E+0, 0.0E+0 ) )
212 * .. Local Scalars ..
214 INTEGER IMAX, ITEMP, J, JB, JJ, JMAX, JP1, JP2, K, KK,
215 $ KW, KKW, KP, KSTEP, P, II
216 REAL ABSAKK, ALPHA, COLMAX, ROWMAX, STEMP, SFMIN
217 COMPLEX D11, D12, D21, D22, R1, T, Z
219 * .. External Functions ..
223 EXTERNAL LSAME, ICAMAX, SLAMCH
225 * .. External Subroutines ..
226 EXTERNAL CCOPY, CGEMM, CGEMV, CSCAL, CSWAP
228 * .. Intrinsic Functions ..
229 INTRINSIC ABS, MAX, MIN, SQRT, AIMAG, REAL
231 * .. Statement Functions ..
234 * .. Statement Function definitions ..
235 CABS1( Z ) = ABS( REAL( Z ) ) + ABS( AIMAG( Z ) )
237 * .. Executable Statements ..
241 * Initialize ALPHA for use in choosing pivot block size.
243 ALPHA = ( ONE+SQRT( SEVTEN ) ) / EIGHT
245 * Compute machine safe minimum
247 SFMIN = SLAMCH( 'S' )
249 IF( LSAME( UPLO, 'U' ) ) THEN
251 * Factorize the trailing columns of A using the upper triangle
252 * of A and working backwards, and compute the matrix W = U12*D
253 * for use in updating A11
255 * K is the main loop index, decreasing from N in steps of 1 or 2
260 * KW is the column of W which corresponds to column K of A
266 IF( ( K.LE.N-NB+1 .AND. NB.LT.N ) .OR. K.LT.1 )
272 * Copy column K of A to column KW of W and update it
274 CALL CCOPY( K, A( 1, K ), 1, W( 1, KW ), 1 )
276 $ CALL CGEMV( 'No transpose', K, N-K, -CONE, A( 1, K+1 ),
277 $ LDA, W( K, KW+1 ), LDW, CONE, W( 1, KW ), 1 )
279 * Determine rows and columns to be interchanged and whether
280 * a 1-by-1 or 2-by-2 pivot block will be used
282 ABSAKK = CABS1( W( K, KW ) )
284 * IMAX is the row-index of the largest off-diagonal element in
285 * column K, and COLMAX is its absolute value.
286 * Determine both COLMAX and IMAX.
289 IMAX = ICAMAX( K-1, W( 1, KW ), 1 )
290 COLMAX = CABS1( W( IMAX, KW ) )
295 IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
297 * Column K is zero or underflow: set INFO and continue
302 CALL CCOPY( K, W( 1, KW ), 1, A( 1, K ), 1 )
305 * ============================================================
307 * Test for interchange
309 * Equivalent to testing for ABSAKK.GE.ALPHA*COLMAX
310 * (used to handle NaN and Inf)
312 IF( .NOT.( ABSAKK.LT.ALPHA*COLMAX ) ) THEN
314 * no interchange, use 1-by-1 pivot block
322 * Loop until pivot found
326 * Begin pivot search loop body
329 * Copy column IMAX to column KW-1 of W and update it
331 CALL CCOPY( IMAX, A( 1, IMAX ), 1, W( 1, KW-1 ), 1 )
332 CALL CCOPY( K-IMAX, A( IMAX, IMAX+1 ), LDA,
333 $ W( IMAX+1, KW-1 ), 1 )
336 $ CALL CGEMV( 'No transpose', K, N-K, -CONE,
337 $ A( 1, K+1 ), LDA, W( IMAX, KW+1 ), LDW,
338 $ CONE, W( 1, KW-1 ), 1 )
340 * JMAX is the column-index of the largest off-diagonal
341 * element in row IMAX, and ROWMAX is its absolute value.
342 * Determine both ROWMAX and JMAX.
345 JMAX = IMAX + ICAMAX( K-IMAX, W( IMAX+1, KW-1 ),
347 ROWMAX = CABS1( W( JMAX, KW-1 ) )
353 ITEMP = ICAMAX( IMAX-1, W( 1, KW-1 ), 1 )
354 STEMP = CABS1( W( ITEMP, KW-1 ) )
355 IF( STEMP.GT.ROWMAX ) THEN
361 * Equivalent to testing for
362 * CABS1( W( IMAX, KW-1 ) ).GE.ALPHA*ROWMAX
363 * (used to handle NaN and Inf)
365 IF( .NOT.(CABS1( W( IMAX, KW-1 ) ).LT.ALPHA*ROWMAX ) )
368 * interchange rows and columns K and IMAX,
369 * use 1-by-1 pivot block
373 * copy column KW-1 of W to column KW of W
375 CALL CCOPY( K, W( 1, KW-1 ), 1, W( 1, KW ), 1 )
379 * Equivalent to testing for ROWMAX.EQ.COLMAX,
380 * (used to handle NaN and Inf)
382 ELSE IF( ( P.EQ.JMAX ) .OR. ( ROWMAX.LE.COLMAX ) )
385 * interchange rows and columns K-1 and IMAX,
386 * use 2-by-2 pivot block
393 * Pivot not found: set params and repeat
399 * Copy updated JMAXth (next IMAXth) column to Kth of W
401 CALL CCOPY( K, W( 1, KW-1 ), 1, W( 1, KW ), 1 )
405 * End pivot search loop body
407 IF( .NOT. DONE ) GOTO 12
411 * ============================================================
415 * KKW is the column of W which corresponds to column KK of A
419 IF( ( KSTEP.EQ.2 ) .AND. ( P.NE.K ) ) THEN
421 * Copy non-updated column K to column P
423 CALL CCOPY( K-P, A( P+1, K ), 1, A( P, P+1 ), LDA )
424 CALL CCOPY( P, A( 1, K ), 1, A( 1, P ), 1 )
426 * Interchange rows K and P in last N-K+1 columns of A
427 * and last N-K+2 columns of W
429 CALL CSWAP( N-K+1, A( K, K ), LDA, A( P, K ), LDA )
430 CALL CSWAP( N-KK+1, W( K, KKW ), LDW, W( P, KKW ), LDW )
433 * Updated column KP is already stored in column KKW of W
437 * Copy non-updated column KK to column KP
439 A( KP, K ) = A( KK, K )
440 CALL CCOPY( K-1-KP, A( KP+1, KK ), 1, A( KP, KP+1 ),
442 CALL CCOPY( KP, A( 1, KK ), 1, A( 1, KP ), 1 )
444 * Interchange rows KK and KP in last N-KK+1 columns
447 CALL CSWAP( N-KK+1, A( KK, KK ), LDA, A( KP, KK ), LDA )
448 CALL CSWAP( N-KK+1, W( KK, KKW ), LDW, W( KP, KKW ),
452 IF( KSTEP.EQ.1 ) THEN
454 * 1-by-1 pivot block D(k): column KW of W now holds
458 * where U(k) is the k-th column of U
460 * Store U(k) in column k of A
462 CALL CCOPY( K, W( 1, KW ), 1, A( 1, K ), 1 )
464 IF( CABS1( A( K, K ) ).GE.SFMIN ) THEN
465 R1 = CONE / A( K, K )
466 CALL CSCAL( K-1, R1, A( 1, K ), 1 )
467 ELSE IF( A( K, K ).NE.CZERO ) THEN
469 A( II, K ) = A( II, K ) / A( K, K )
476 * 2-by-2 pivot block D(k): columns KW and KW-1 of W now
479 * ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
481 * where U(k) and U(k-1) are the k-th and (k-1)-th columns
486 * Store U(k) and U(k-1) in columns k and k-1 of A
489 D11 = W( K, KW ) / D12
490 D22 = W( K-1, KW-1 ) / D12
491 T = CONE / ( D11*D22-CONE )
493 A( J, K-1 ) = T*( (D11*W( J, KW-1 )-W( J, KW ) ) /
495 A( J, K ) = T*( ( D22*W( J, KW )-W( J, KW-1 ) ) /
502 A( K-1, K-1 ) = W( K-1, KW-1 )
503 A( K-1, K ) = W( K-1, KW )
504 A( K, K ) = W( K, KW )
508 * Store details of the interchanges in IPIV
510 IF( KSTEP.EQ.1 ) THEN
517 * Decrease K and return to the start of the main loop
524 * Update the upper triangle of A11 (= A(1:k,1:k)) as
526 * A11 := A11 - U12*D*U12**T = A11 - U12*W**T
528 * computing blocks of NB columns at a time
530 DO 50 J = ( ( K-1 ) / NB )*NB + 1, 1, -NB
531 JB = MIN( NB, K-J+1 )
533 * Update the upper triangle of the diagonal block
535 DO 40 JJ = J, J + JB - 1
536 CALL CGEMV( 'No transpose', JJ-J+1, N-K, -CONE,
537 $ A( J, K+1 ), LDA, W( JJ, KW+1 ), LDW, CONE,
541 * Update the rectangular superdiagonal block
544 $ CALL CGEMM( 'No transpose', 'Transpose', J-1, JB,
545 $ N-K, -CONE, A( 1, K+1 ), LDA, W( J, KW+1 ), LDW,
546 $ CONE, A( 1, J ), LDA )
549 * Put U12 in standard form by partially undoing the interchanges
567 IF( JP2.NE.JJ .AND. J.LE.N )
568 $ CALL CSWAP( N-J+1, A( JP2, J ), LDA, A( JJ, J ), LDA )
570 IF( JP1.NE.JJ .AND. KSTEP.EQ.2 )
571 $ CALL CSWAP( N-J+1, A( JP1, J ), LDA, A( JJ, J ), LDA )
575 * Set KB to the number of columns factorized
581 * Factorize the leading columns of A using the lower triangle
582 * of A and working forwards, and compute the matrix W = L21*D
583 * for use in updating A22
585 * K is the main loop index, increasing from 1 in steps of 1 or 2
592 IF( ( K.GE.NB .AND. NB.LT.N ) .OR. K.GT.N )
598 * Copy column K of A to column K of W and update it
600 CALL CCOPY( N-K+1, A( K, K ), 1, W( K, K ), 1 )
602 $ CALL CGEMV( 'No transpose', N-K+1, K-1, -CONE, A( K, 1 ),
603 $ LDA, W( K, 1 ), LDW, CONE, W( K, K ), 1 )
605 * Determine rows and columns to be interchanged and whether
606 * a 1-by-1 or 2-by-2 pivot block will be used
608 ABSAKK = CABS1( W( K, K ) )
610 * IMAX is the row-index of the largest off-diagonal element in
611 * column K, and COLMAX is its absolute value.
612 * Determine both COLMAX and IMAX.
615 IMAX = K + ICAMAX( N-K, W( K+1, K ), 1 )
616 COLMAX = CABS1( W( IMAX, K ) )
621 IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
623 * Column K is zero or underflow: set INFO and continue
628 CALL CCOPY( N-K+1, W( K, K ), 1, A( K, K ), 1 )
631 * ============================================================
633 * Test for interchange
635 * Equivalent to testing for ABSAKK.GE.ALPHA*COLMAX
636 * (used to handle NaN and Inf)
638 IF( .NOT.( ABSAKK.LT.ALPHA*COLMAX ) ) THEN
640 * no interchange, use 1-by-1 pivot block
648 * Loop until pivot found
652 * Begin pivot search loop body
655 * Copy column IMAX to column K+1 of W and update it
657 CALL CCOPY( IMAX-K, A( IMAX, K ), LDA, W( K, K+1 ), 1)
658 CALL CCOPY( N-IMAX+1, A( IMAX, IMAX ), 1,
659 $ W( IMAX, K+1 ), 1 )
661 $ CALL CGEMV( 'No transpose', N-K+1, K-1, -CONE,
662 $ A( K, 1 ), LDA, W( IMAX, 1 ), LDW,
663 $ CONE, W( K, K+1 ), 1 )
665 * JMAX is the column-index of the largest off-diagonal
666 * element in row IMAX, and ROWMAX is its absolute value.
667 * Determine both ROWMAX and JMAX.
670 JMAX = K - 1 + ICAMAX( IMAX-K, W( K, K+1 ), 1 )
671 ROWMAX = CABS1( W( JMAX, K+1 ) )
677 ITEMP = IMAX + ICAMAX( N-IMAX, W( IMAX+1, K+1 ), 1)
678 STEMP = CABS1( W( ITEMP, K+1 ) )
679 IF( STEMP.GT.ROWMAX ) THEN
685 * Equivalent to testing for
686 * CABS1( W( IMAX, K+1 ) ).GE.ALPHA*ROWMAX
687 * (used to handle NaN and Inf)
689 IF( .NOT.( CABS1( W( IMAX, K+1 ) ).LT.ALPHA*ROWMAX ) )
692 * interchange rows and columns K and IMAX,
693 * use 1-by-1 pivot block
697 * copy column K+1 of W to column K of W
699 CALL CCOPY( N-K+1, W( K, K+1 ), 1, W( K, K ), 1 )
703 * Equivalent to testing for ROWMAX.EQ.COLMAX,
704 * (used to handle NaN and Inf)
706 ELSE IF( ( P.EQ.JMAX ) .OR. ( ROWMAX.LE.COLMAX ) )
709 * interchange rows and columns K+1 and IMAX,
710 * use 2-by-2 pivot block
717 * Pivot not found: set params and repeat
723 * Copy updated JMAXth (next IMAXth) column to Kth of W
725 CALL CCOPY( N-K+1, W( K, K+1 ), 1, W( K, K ), 1 )
729 * End pivot search loop body
731 IF( .NOT. DONE ) GOTO 72
735 * ============================================================
739 IF( ( KSTEP.EQ.2 ) .AND. ( P.NE.K ) ) THEN
741 * Copy non-updated column K to column P
743 CALL CCOPY( P-K, A( K, K ), 1, A( P, K ), LDA )
744 CALL CCOPY( N-P+1, A( P, K ), 1, A( P, P ), 1 )
746 * Interchange rows K and P in first K columns of A
747 * and first K+1 columns of W
749 CALL CSWAP( K, A( K, 1 ), LDA, A( P, 1 ), LDA )
750 CALL CSWAP( KK, W( K, 1 ), LDW, W( P, 1 ), LDW )
753 * Updated column KP is already stored in column KK of W
757 * Copy non-updated column KK to column KP
759 A( KP, K ) = A( KK, K )
760 CALL CCOPY( KP-K-1, A( K+1, KK ), 1, A( KP, K+1 ), LDA )
761 CALL CCOPY( N-KP+1, A( KP, KK ), 1, A( KP, KP ), 1 )
763 * Interchange rows KK and KP in first KK columns of A and W
765 CALL CSWAP( KK, A( KK, 1 ), LDA, A( KP, 1 ), LDA )
766 CALL CSWAP( KK, W( KK, 1 ), LDW, W( KP, 1 ), LDW )
769 IF( KSTEP.EQ.1 ) THEN
771 * 1-by-1 pivot block D(k): column k of W now holds
775 * where L(k) is the k-th column of L
777 * Store L(k) in column k of A
779 CALL CCOPY( N-K+1, W( K, K ), 1, A( K, K ), 1 )
781 IF( CABS1( A( K, K ) ).GE.SFMIN ) THEN
782 R1 = CONE / A( K, K )
783 CALL CSCAL( N-K, R1, A( K+1, K ), 1 )
784 ELSE IF( A( K, K ).NE.CZERO ) THEN
786 A( II, K ) = A( II, K ) / A( K, K )
793 * 2-by-2 pivot block D(k): columns k and k+1 of W now hold
795 * ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k)
797 * where L(k) and L(k+1) are the k-th and (k+1)-th columns
802 * Store L(k) and L(k+1) in columns k and k+1 of A
805 D11 = W( K+1, K+1 ) / D21
806 D22 = W( K, K ) / D21
807 T = CONE / ( D11*D22-CONE )
809 A( J, K ) = T*( ( D11*W( J, K )-W( J, K+1 ) ) /
811 A( J, K+1 ) = T*( ( D22*W( J, K+1 )-W( J, K ) ) /
818 A( K, K ) = W( K, K )
819 A( K+1, K ) = W( K+1, K )
820 A( K+1, K+1 ) = W( K+1, K+1 )
824 * Store details of the interchanges in IPIV
826 IF( KSTEP.EQ.1 ) THEN
833 * Increase K and return to the start of the main loop
840 * Update the lower triangle of A22 (= A(k:n,k:n)) as
842 * A22 := A22 - L21*D*L21**T = A22 - L21*W**T
844 * computing blocks of NB columns at a time
847 JB = MIN( NB, N-J+1 )
849 * Update the lower triangle of the diagonal block
851 DO 100 JJ = J, J + JB - 1
852 CALL CGEMV( 'No transpose', J+JB-JJ, K-1, -CONE,
853 $ A( JJ, 1 ), LDA, W( JJ, 1 ), LDW, CONE,
857 * Update the rectangular subdiagonal block
860 $ CALL CGEMM( 'No transpose', 'Transpose', N-J-JB+1, JB,
861 $ K-1, -CONE, A( J+JB, 1 ), LDA, W( J, 1 ), LDW,
862 $ CONE, A( J+JB, J ), LDA )
865 * Put L21 in standard form by partially undoing the interchanges
883 IF( JP2.NE.JJ .AND. J.GE.1 )
884 $ CALL CSWAP( J, A( JP2, 1 ), LDA, A( JJ, 1 ), LDA )
886 IF( JP1.NE.JJ .AND. KSTEP.EQ.2 )
887 $ CALL CSWAP( J, A( JP1, 1 ), LDA, A( JJ, 1 ), LDA )
891 * Set KB to the number of columns factorized