1 *> \brief \b CHETF2 computes the factorization of a complex Hermitian matrix, using the diagonal pivoting method (unblocked algorithm calling Level 2 BLAS).
3 * =========== DOCUMENTATION ===========
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
9 *> Download CHETF2 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/chetf2.f">
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/chetf2.f">
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/chetf2.f">
21 * SUBROUTINE CHETF2( UPLO, N, A, LDA, IPIV, INFO )
23 * .. Scalar Arguments ..
25 * INTEGER INFO, LDA, N
27 * .. Array Arguments ..
38 *> CHETF2 computes the factorization of a complex Hermitian matrix A
39 *> using the Bunch-Kaufman diagonal pivoting method:
41 *> A = U*D*U**H or A = L*D*L**H
43 *> where U (or L) is a product of permutation and unit upper (lower)
44 *> triangular matrices, U**H is the conjugate transpose of U, and D is
45 *> Hermitian and block diagonal with 1-by-1 and 2-by-2 diagonal blocks.
47 *> This is the unblocked version of the algorithm, calling Level 2 BLAS.
55 *> UPLO is CHARACTER*1
56 *> Specifies whether the upper or lower triangular part of the
57 *> Hermitian matrix A is stored:
58 *> = 'U': Upper triangular
59 *> = 'L': Lower triangular
65 *> The order of the matrix A. N >= 0.
70 *> A is COMPLEX array, dimension (LDA,N)
71 *> On entry, the Hermitian matrix A. If UPLO = 'U', the leading
72 *> n-by-n upper triangular part of A contains the upper
73 *> triangular part of the matrix A, and the strictly lower
74 *> triangular part of A is not referenced. If UPLO = 'L', the
75 *> leading n-by-n lower triangular part of A contains the lower
76 *> triangular part of the matrix A, and the strictly upper
77 *> triangular part of A is not referenced.
79 *> On exit, the block diagonal matrix D and the multipliers used
80 *> to obtain the factor U or L (see below for further details).
86 *> The leading dimension of the array A. LDA >= max(1,N).
91 *> IPIV is INTEGER array, dimension (N)
92 *> Details of the interchanges and the block structure of D.
95 *> If IPIV(k) > 0, then rows and columns k and IPIV(k) were
96 *> interchanged and D(k,k) is a 1-by-1 diagonal block.
98 *> If IPIV(k) = IPIV(k-1) < 0, then rows and columns
99 *> k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
100 *> is a 2-by-2 diagonal block.
103 *> If IPIV(k) > 0, then rows and columns k and IPIV(k) were
104 *> interchanged and D(k,k) is a 1-by-1 diagonal block.
106 *> If IPIV(k) = IPIV(k+1) < 0, then rows and columns
107 *> k+1 and -IPIV(k) were interchanged and D(k:k+1,k:k+1)
108 *> is a 2-by-2 diagonal block.
114 *> = 0: successful exit
115 *> < 0: if INFO = -k, the k-th argument had an illegal value
116 *> > 0: if INFO = k, D(k,k) is exactly zero. The factorization
117 *> has been completed, but the block diagonal matrix D is
118 *> exactly singular, and division by zero will occur if it
119 *> is used to solve a system of equations.
125 *> \author Univ. of Tennessee
126 *> \author Univ. of California Berkeley
127 *> \author Univ. of Colorado Denver
130 *> \date November 2013
132 *> \ingroup complexHEcomputational
134 *> \par Further Details:
135 * =====================
139 *> 09-29-06 - patch from
140 *> Bobby Cheng, MathWorks
142 *> Replace l.210 and l.392
143 *> IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
145 *> IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. SISNAN(ABSAKK) ) THEN
147 *> 01-01-96 - Based on modifications by
148 *> J. Lewis, Boeing Computer Services Company
149 *> A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA
151 *> If UPLO = 'U', then A = U*D*U**H, where
152 *> U = P(n)*U(n)* ... *P(k)U(k)* ...,
153 *> i.e., U is a product of terms P(k)*U(k), where k decreases from n to
154 *> 1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
155 *> and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as
156 *> defined by IPIV(k), and U(k) is a unit upper triangular matrix, such
157 *> that if the diagonal block D(k) is of order s (s = 1 or 2), then
160 *> U(k) = ( 0 I 0 ) s
164 *> If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k).
165 *> If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k),
166 *> and A(k,k), and v overwrites A(1:k-2,k-1:k).
168 *> If UPLO = 'L', then A = L*D*L**H, where
169 *> L = P(1)*L(1)* ... *P(k)*L(k)* ...,
170 *> i.e., L is a product of terms P(k)*L(k), where k increases from 1 to
171 *> n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
172 *> and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as
173 *> defined by IPIV(k), and L(k) is a unit lower triangular matrix, such
174 *> that if the diagonal block D(k) is of order s (s = 1 or 2), then
177 *> L(k) = ( 0 I 0 ) s
181 *> If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k).
182 *> If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k),
183 *> and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1).
186 * =====================================================================
187 SUBROUTINE CHETF2( UPLO, N, A, LDA, IPIV, INFO )
189 * -- LAPACK computational routine (version 3.5.0) --
190 * -- LAPACK is a software package provided by Univ. of Tennessee, --
191 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
194 * .. Scalar Arguments ..
198 * .. Array Arguments ..
203 * =====================================================================
207 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 )
209 PARAMETER ( EIGHT = 8.0E+0, SEVTEN = 17.0E+0 )
211 * .. Local Scalars ..
213 INTEGER I, IMAX, J, JMAX, K, KK, KP, KSTEP
214 REAL ABSAKK, ALPHA, COLMAX, D, D11, D22, R1, ROWMAX,
216 COMPLEX D12, D21, T, WK, WKM1, WKP1, ZDUM
218 * .. External Functions ..
219 LOGICAL LSAME, SISNAN
222 EXTERNAL LSAME, ICAMAX, SLAPY2, SISNAN
224 * .. External Subroutines ..
225 EXTERNAL CHER, CSSCAL, CSWAP, XERBLA
227 * .. Intrinsic Functions ..
228 INTRINSIC ABS, AIMAG, CMPLX, CONJG, MAX, REAL, SQRT
230 * .. Statement Functions ..
233 * .. Statement Function definitions ..
234 CABS1( ZDUM ) = ABS( REAL( ZDUM ) ) + ABS( AIMAG( ZDUM ) )
236 * .. Executable Statements ..
238 * Test the input parameters.
241 UPPER = LSAME( UPLO, 'U' )
242 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
244 ELSE IF( N.LT.0 ) THEN
246 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
250 CALL XERBLA( 'CHETF2', -INFO )
254 * Initialize ALPHA for use in choosing pivot block size.
256 ALPHA = ( ONE+SQRT( SEVTEN ) ) / EIGHT
260 * Factorize A as U*D*U**H using the upper triangle of A
262 * K is the main loop index, decreasing from N to 1 in steps of
268 * If K < 1, exit from loop
274 * Determine rows and columns to be interchanged and whether
275 * a 1-by-1 or 2-by-2 pivot block will be used
277 ABSAKK = ABS( REAL( A( K, K ) ) )
279 * IMAX is the row-index of the largest off-diagonal element in
280 * column K, and COLMAX is its absolute value.
281 * Determine both COLMAX and IMAX.
284 IMAX = ICAMAX( K-1, A( 1, K ), 1 )
285 COLMAX = CABS1( A( IMAX, K ) )
290 IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. SISNAN(ABSAKK) ) THEN
292 * Column K is or underflow, or contains a NaN:
293 * set INFO and continue
298 A( K, K ) = REAL( A( K, K ) )
300 IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
302 * no interchange, use 1-by-1 pivot block
307 * JMAX is the column-index of the largest off-diagonal
308 * element in row IMAX, and ROWMAX is its absolute value
310 JMAX = IMAX + ICAMAX( K-IMAX, A( IMAX, IMAX+1 ), LDA )
311 ROWMAX = CABS1( A( IMAX, JMAX ) )
313 JMAX = ICAMAX( IMAX-1, A( 1, IMAX ), 1 )
314 ROWMAX = MAX( ROWMAX, CABS1( A( JMAX, IMAX ) ) )
317 IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
319 * no interchange, use 1-by-1 pivot block
322 ELSE IF( ABS( REAL( A( IMAX, IMAX ) ) ).GE.ALPHA*ROWMAX )
325 * interchange rows and columns K and IMAX, use 1-by-1
331 * interchange rows and columns K-1 and IMAX, use 2-by-2
342 * Interchange rows and columns KK and KP in the leading
343 * submatrix A(1:k,1:k)
345 CALL CSWAP( KP-1, A( 1, KK ), 1, A( 1, KP ), 1 )
346 DO 20 J = KP + 1, KK - 1
347 T = CONJG( A( J, KK ) )
348 A( J, KK ) = CONJG( A( KP, J ) )
351 A( KP, KK ) = CONJG( A( KP, KK ) )
352 R1 = REAL( A( KK, KK ) )
353 A( KK, KK ) = REAL( A( KP, KP ) )
355 IF( KSTEP.EQ.2 ) THEN
356 A( K, K ) = REAL( A( K, K ) )
358 A( K-1, K ) = A( KP, K )
362 A( K, K ) = REAL( A( K, K ) )
364 $ A( K-1, K-1 ) = REAL( A( K-1, K-1 ) )
367 * Update the leading submatrix
369 IF( KSTEP.EQ.1 ) THEN
371 * 1-by-1 pivot block D(k): column k now holds
375 * where U(k) is the k-th column of U
377 * Perform a rank-1 update of A(1:k-1,1:k-1) as
379 * A := A - U(k)*D(k)*U(k)**H = A - W(k)*1/D(k)*W(k)**H
381 R1 = ONE / REAL( A( K, K ) )
382 CALL CHER( UPLO, K-1, -R1, A( 1, K ), 1, A, LDA )
384 * Store U(k) in column k
386 CALL CSSCAL( K-1, R1, A( 1, K ), 1 )
389 * 2-by-2 pivot block D(k): columns k and k-1 now hold
391 * ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
393 * where U(k) and U(k-1) are the k-th and (k-1)-th columns
396 * Perform a rank-2 update of A(1:k-2,1:k-2) as
398 * A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )**H
399 * = A - ( W(k-1) W(k) )*inv(D(k))*( W(k-1) W(k) )**H
403 D = SLAPY2( REAL( A( K-1, K ) ),
404 $ AIMAG( A( K-1, K ) ) )
405 D22 = REAL( A( K-1, K-1 ) ) / D
406 D11 = REAL( A( K, K ) ) / D
407 TT = ONE / ( D11*D22-ONE )
408 D12 = A( K-1, K ) / D
411 DO 40 J = K - 2, 1, -1
412 WKM1 = D*( D11*A( J, K-1 )-CONJG( D12 )*A( J, K ) )
413 WK = D*( D22*A( J, K )-D12*A( J, K-1 ) )
415 A( I, J ) = A( I, J ) - A( I, K )*CONJG( WK ) -
416 $ A( I, K-1 )*CONJG( WKM1 )
420 A( J, J ) = CMPLX( REAL( A( J, J ) ), 0.0E+0 )
428 * Store details of the interchanges in IPIV
430 IF( KSTEP.EQ.1 ) THEN
437 * Decrease K and return to the start of the main loop
444 * Factorize A as L*D*L**H using the lower triangle of A
446 * K is the main loop index, increasing from 1 to N in steps of
452 * If K > N, exit from loop
458 * Determine rows and columns to be interchanged and whether
459 * a 1-by-1 or 2-by-2 pivot block will be used
461 ABSAKK = ABS( REAL( A( K, K ) ) )
463 * IMAX is the row-index of the largest off-diagonal element in
464 * column K, and COLMAX is its absolute value.
465 * Determine both COLMAX and IMAX.
468 IMAX = K + ICAMAX( N-K, A( K+1, K ), 1 )
469 COLMAX = CABS1( A( IMAX, K ) )
474 IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. SISNAN(ABSAKK) ) THEN
476 * Column K is zero or underflow, contains a NaN:
477 * set INFO and continue
482 A( K, K ) = REAL( A( K, K ) )
484 IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
486 * no interchange, use 1-by-1 pivot block
491 * JMAX is the column-index of the largest off-diagonal
492 * element in row IMAX, and ROWMAX is its absolute value
494 JMAX = K - 1 + ICAMAX( IMAX-K, A( IMAX, K ), LDA )
495 ROWMAX = CABS1( A( IMAX, JMAX ) )
497 JMAX = IMAX + ICAMAX( N-IMAX, A( IMAX+1, IMAX ), 1 )
498 ROWMAX = MAX( ROWMAX, CABS1( A( JMAX, IMAX ) ) )
501 IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
503 * no interchange, use 1-by-1 pivot block
506 ELSE IF( ABS( REAL( A( IMAX, IMAX ) ) ).GE.ALPHA*ROWMAX )
509 * interchange rows and columns K and IMAX, use 1-by-1
515 * interchange rows and columns K+1 and IMAX, use 2-by-2
526 * Interchange rows and columns KK and KP in the trailing
527 * submatrix A(k:n,k:n)
530 $ CALL CSWAP( N-KP, A( KP+1, KK ), 1, A( KP+1, KP ), 1 )
531 DO 60 J = KK + 1, KP - 1
532 T = CONJG( A( J, KK ) )
533 A( J, KK ) = CONJG( A( KP, J ) )
536 A( KP, KK ) = CONJG( A( KP, KK ) )
537 R1 = REAL( A( KK, KK ) )
538 A( KK, KK ) = REAL( A( KP, KP ) )
540 IF( KSTEP.EQ.2 ) THEN
541 A( K, K ) = REAL( A( K, K ) )
543 A( K+1, K ) = A( KP, K )
547 A( K, K ) = REAL( A( K, K ) )
549 $ A( K+1, K+1 ) = REAL( A( K+1, K+1 ) )
552 * Update the trailing submatrix
554 IF( KSTEP.EQ.1 ) THEN
556 * 1-by-1 pivot block D(k): column k now holds
560 * where L(k) is the k-th column of L
564 * Perform a rank-1 update of A(k+1:n,k+1:n) as
566 * A := A - L(k)*D(k)*L(k)**H = A - W(k)*(1/D(k))*W(k)**H
568 R1 = ONE / REAL( A( K, K ) )
569 CALL CHER( UPLO, N-K, -R1, A( K+1, K ), 1,
570 $ A( K+1, K+1 ), LDA )
572 * Store L(k) in column K
574 CALL CSSCAL( N-K, R1, A( K+1, K ), 1 )
578 * 2-by-2 pivot block D(k)
582 * Perform a rank-2 update of A(k+2:n,k+2:n) as
584 * A := A - ( L(k) L(k+1) )*D(k)*( L(k) L(k+1) )**H
585 * = A - ( W(k) W(k+1) )*inv(D(k))*( W(k) W(k+1) )**H
587 * where L(k) and L(k+1) are the k-th and (k+1)-th
590 D = SLAPY2( REAL( A( K+1, K ) ),
591 $ AIMAG( A( K+1, K ) ) )
592 D11 = REAL( A( K+1, K+1 ) ) / D
593 D22 = REAL( A( K, K ) ) / D
594 TT = ONE / ( D11*D22-ONE )
595 D21 = A( K+1, K ) / D
599 WK = D*( D11*A( J, K )-D21*A( J, K+1 ) )
600 WKP1 = D*( D22*A( J, K+1 )-CONJG( D21 )*A( J, K ) )
602 A( I, J ) = A( I, J ) - A( I, K )*CONJG( WK ) -
603 $ A( I, K+1 )*CONJG( WKP1 )
607 A( J, J ) = CMPLX( REAL( A( J, J ) ), 0.0E+0 )
613 * Store details of the interchanges in IPIV
615 IF( KSTEP.EQ.1 ) THEN
622 * Increase K and return to the start of the main loop