3 * =========== DOCUMENTATION ===========
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
9 *> Download ZHPTRF + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhptrf.f">
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhptrf.f">
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhptrf.f">
21 * SUBROUTINE ZHPTRF( UPLO, N, AP, IPIV, INFO )
23 * .. Scalar Arguments ..
27 * .. Array Arguments ..
38 *> ZHPTRF computes the factorization of a complex Hermitian packed
39 *> matrix A 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, and D is Hermitian and block diagonal with
45 *> 1-by-1 and 2-by-2 diagonal blocks.
53 *> UPLO is CHARACTER*1
54 *> = 'U': Upper triangle of A is stored;
55 *> = 'L': Lower triangle of A is stored.
61 *> The order of the matrix A. N >= 0.
66 *> AP is COMPLEX*16 array, dimension (N*(N+1)/2)
67 *> On entry, the upper or lower triangle of the Hermitian matrix
68 *> A, packed columnwise in a linear array. The j-th column of A
69 *> is stored in the array AP as follows:
70 *> if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
71 *> if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
73 *> On exit, the block diagonal matrix D and the multipliers used
74 *> to obtain the factor U or L, stored as a packed triangular
75 *> matrix overwriting A (see below for further details).
80 *> IPIV is INTEGER array, dimension (N)
81 *> Details of the interchanges and the block structure of D.
82 *> If IPIV(k) > 0, then rows and columns k and IPIV(k) were
83 *> interchanged and D(k,k) is a 1-by-1 diagonal block.
84 *> If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and
85 *> columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
86 *> is a 2-by-2 diagonal block. If UPLO = 'L' and IPIV(k) =
87 *> IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were
88 *> interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
94 *> = 0: successful exit
95 *> < 0: if INFO = -i, the i-th argument had an illegal value
96 *> > 0: if INFO = i, D(i,i) is exactly zero. The factorization
97 *> has been completed, but the block diagonal matrix D is
98 *> exactly singular, and division by zero will occur if it
99 *> is used to solve a system of equations.
105 *> \author Univ. of Tennessee
106 *> \author Univ. of California Berkeley
107 *> \author Univ. of Colorado Denver
110 *> \date November 2011
112 *> \ingroup complex16OTHERcomputational
114 *> \par Further Details:
115 * =====================
119 *> If UPLO = 'U', then A = U*D*U**H, where
120 *> U = P(n)*U(n)* ... *P(k)U(k)* ...,
121 *> i.e., U is a product of terms P(k)*U(k), where k decreases from n to
122 *> 1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
123 *> and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as
124 *> defined by IPIV(k), and U(k) is a unit upper triangular matrix, such
125 *> that if the diagonal block D(k) is of order s (s = 1 or 2), then
128 *> U(k) = ( 0 I 0 ) s
132 *> If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k).
133 *> If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k),
134 *> and A(k,k), and v overwrites A(1:k-2,k-1:k).
136 *> If UPLO = 'L', then A = L*D*L**H, where
137 *> L = P(1)*L(1)* ... *P(k)*L(k)* ...,
138 *> i.e., L is a product of terms P(k)*L(k), where k increases from 1 to
139 *> n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
140 *> and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as
141 *> defined by IPIV(k), and L(k) is a unit lower triangular matrix, such
142 *> that if the diagonal block D(k) is of order s (s = 1 or 2), then
145 *> L(k) = ( 0 I 0 ) s
149 *> If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k).
150 *> If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k),
151 *> and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1).
154 *> \par Contributors:
157 *> J. Lewis, Boeing Computer Services Company
159 * =====================================================================
160 SUBROUTINE ZHPTRF( UPLO, N, AP, IPIV, INFO )
162 * -- LAPACK computational routine (version 3.4.0) --
163 * -- LAPACK is a software package provided by Univ. of Tennessee, --
164 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
167 * .. Scalar Arguments ..
171 * .. Array Arguments ..
176 * =====================================================================
179 DOUBLE PRECISION ZERO, ONE
180 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
181 DOUBLE PRECISION EIGHT, SEVTEN
182 PARAMETER ( EIGHT = 8.0D+0, SEVTEN = 17.0D+0 )
184 * .. Local Scalars ..
186 INTEGER I, IMAX, J, JMAX, K, KC, KK, KNC, KP, KPC,
188 DOUBLE PRECISION ABSAKK, ALPHA, COLMAX, D, D11, D22, R1, ROWMAX,
190 COMPLEX*16 D12, D21, T, WK, WKM1, WKP1, ZDUM
192 * .. External Functions ..
195 DOUBLE PRECISION DLAPY2
196 EXTERNAL LSAME, IZAMAX, DLAPY2
198 * .. External Subroutines ..
199 EXTERNAL XERBLA, ZDSCAL, ZHPR, ZSWAP
201 * .. Intrinsic Functions ..
202 INTRINSIC ABS, DBLE, DCMPLX, DCONJG, DIMAG, MAX, SQRT
204 * .. Statement Functions ..
205 DOUBLE PRECISION CABS1
207 * .. Statement Function definitions ..
208 CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
210 * .. Executable Statements ..
212 * Test the input parameters.
215 UPPER = LSAME( UPLO, 'U' )
216 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
218 ELSE IF( N.LT.0 ) THEN
222 CALL XERBLA( 'ZHPTRF', -INFO )
226 * Initialize ALPHA for use in choosing pivot block size.
228 ALPHA = ( ONE+SQRT( SEVTEN ) ) / EIGHT
232 * Factorize A as U*D*U**H using the upper triangle of A
234 * K is the main loop index, decreasing from N to 1 in steps of
238 KC = ( N-1 )*N / 2 + 1
242 * If K < 1, exit from loop
248 * Determine rows and columns to be interchanged and whether
249 * a 1-by-1 or 2-by-2 pivot block will be used
251 ABSAKK = ABS( DBLE( AP( KC+K-1 ) ) )
253 * IMAX is the row-index of the largest off-diagonal element in
254 * column K, and COLMAX is its absolute value
257 IMAX = IZAMAX( K-1, AP( KC ), 1 )
258 COLMAX = CABS1( AP( KC+IMAX-1 ) )
263 IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
265 * Column K is zero: set INFO and continue
270 AP( KC+K-1 ) = DBLE( AP( KC+K-1 ) )
272 IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
274 * no interchange, use 1-by-1 pivot block
279 * JMAX is the column-index of the largest off-diagonal
280 * element in row IMAX, and ROWMAX is its absolute value
284 KX = IMAX*( IMAX+1 ) / 2 + IMAX
285 DO 20 J = IMAX + 1, K
286 IF( CABS1( AP( KX ) ).GT.ROWMAX ) THEN
287 ROWMAX = CABS1( AP( KX ) )
292 KPC = ( IMAX-1 )*IMAX / 2 + 1
294 JMAX = IZAMAX( IMAX-1, AP( KPC ), 1 )
295 ROWMAX = MAX( ROWMAX, CABS1( AP( KPC+JMAX-1 ) ) )
298 IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
300 * no interchange, use 1-by-1 pivot block
303 ELSE IF( ABS( DBLE( AP( KPC+IMAX-1 ) ) ).GE.ALPHA*
306 * interchange rows and columns K and IMAX, use 1-by-1
312 * interchange rows and columns K-1 and IMAX, use 2-by-2
325 * Interchange rows and columns KK and KP in the leading
326 * submatrix A(1:k,1:k)
328 CALL ZSWAP( KP-1, AP( KNC ), 1, AP( KPC ), 1 )
330 DO 30 J = KP + 1, KK - 1
332 T = DCONJG( AP( KNC+J-1 ) )
333 AP( KNC+J-1 ) = DCONJG( AP( KX ) )
336 AP( KX+KK-1 ) = DCONJG( AP( KX+KK-1 ) )
337 R1 = DBLE( AP( KNC+KK-1 ) )
338 AP( KNC+KK-1 ) = DBLE( AP( KPC+KP-1 ) )
340 IF( KSTEP.EQ.2 ) THEN
341 AP( KC+K-1 ) = DBLE( AP( KC+K-1 ) )
343 AP( KC+K-2 ) = AP( KC+KP-1 )
347 AP( KC+K-1 ) = DBLE( AP( KC+K-1 ) )
349 $ AP( KC-1 ) = DBLE( AP( KC-1 ) )
352 * Update the leading submatrix
354 IF( KSTEP.EQ.1 ) THEN
356 * 1-by-1 pivot block D(k): column k now holds
360 * where U(k) is the k-th column of U
362 * Perform a rank-1 update of A(1:k-1,1:k-1) as
364 * A := A - U(k)*D(k)*U(k)**H = A - W(k)*1/D(k)*W(k)**H
366 R1 = ONE / DBLE( AP( KC+K-1 ) )
367 CALL ZHPR( UPLO, K-1, -R1, AP( KC ), 1, AP )
369 * Store U(k) in column k
371 CALL ZDSCAL( K-1, R1, AP( KC ), 1 )
374 * 2-by-2 pivot block D(k): columns k and k-1 now hold
376 * ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
378 * where U(k) and U(k-1) are the k-th and (k-1)-th columns
381 * Perform a rank-2 update of A(1:k-2,1:k-2) as
383 * A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )**H
384 * = A - ( W(k-1) W(k) )*inv(D(k))*( W(k-1) W(k) )**H
388 D = DLAPY2( DBLE( AP( K-1+( K-1 )*K / 2 ) ),
389 $ DIMAG( AP( K-1+( K-1 )*K / 2 ) ) )
390 D22 = DBLE( AP( K-1+( K-2 )*( K-1 ) / 2 ) ) / D
391 D11 = DBLE( AP( K+( K-1 )*K / 2 ) ) / D
392 TT = ONE / ( D11*D22-ONE )
393 D12 = AP( K-1+( K-1 )*K / 2 ) / D
396 DO 50 J = K - 2, 1, -1
397 WKM1 = D*( D11*AP( J+( K-2 )*( K-1 ) / 2 )-
398 $ DCONJG( D12 )*AP( J+( K-1 )*K / 2 ) )
399 WK = D*( D22*AP( J+( K-1 )*K / 2 )-D12*
400 $ AP( J+( K-2 )*( K-1 ) / 2 ) )
402 AP( I+( J-1 )*J / 2 ) = AP( I+( J-1 )*J / 2 ) -
403 $ AP( I+( K-1 )*K / 2 )*DCONJG( WK ) -
404 $ AP( I+( K-2 )*( K-1 ) / 2 )*DCONJG( WKM1 )
406 AP( J+( K-1 )*K / 2 ) = WK
407 AP( J+( K-2 )*( K-1 ) / 2 ) = WKM1
408 AP( J+( J-1 )*J / 2 ) = DCMPLX( DBLE( AP( J+( J-
409 $ 1 )*J / 2 ) ), 0.0D+0 )
417 * Store details of the interchanges in IPIV
419 IF( KSTEP.EQ.1 ) THEN
426 * Decrease K and return to the start of the main loop
434 * Factorize A as L*D*L**H using the lower triangle of A
436 * K is the main loop index, increasing from 1 to N in steps of
445 * If K > N, exit from loop
451 * Determine rows and columns to be interchanged and whether
452 * a 1-by-1 or 2-by-2 pivot block will be used
454 ABSAKK = ABS( DBLE( AP( KC ) ) )
456 * IMAX is the row-index of the largest off-diagonal element in
457 * column K, and COLMAX is its absolute value
460 IMAX = K + IZAMAX( N-K, AP( KC+1 ), 1 )
461 COLMAX = CABS1( AP( KC+IMAX-K ) )
466 IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
468 * Column K is zero: set INFO and continue
473 AP( KC ) = DBLE( AP( KC ) )
475 IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
477 * no interchange, use 1-by-1 pivot block
482 * JMAX is the column-index of the largest off-diagonal
483 * element in row IMAX, and ROWMAX is its absolute value
487 DO 70 J = K, IMAX - 1
488 IF( CABS1( AP( KX ) ).GT.ROWMAX ) THEN
489 ROWMAX = CABS1( AP( KX ) )
494 KPC = NPP - ( N-IMAX+1 )*( N-IMAX+2 ) / 2 + 1
496 JMAX = IMAX + IZAMAX( N-IMAX, AP( KPC+1 ), 1 )
497 ROWMAX = MAX( ROWMAX, CABS1( AP( KPC+JMAX-IMAX ) ) )
500 IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
502 * no interchange, use 1-by-1 pivot block
505 ELSE IF( ABS( DBLE( AP( KPC ) ) ).GE.ALPHA*ROWMAX ) THEN
507 * interchange rows and columns K and IMAX, use 1-by-1
513 * interchange rows and columns K+1 and IMAX, use 2-by-2
523 $ KNC = KNC + N - K + 1
526 * Interchange rows and columns KK and KP in the trailing
527 * submatrix A(k:n,k:n)
530 $ CALL ZSWAP( N-KP, AP( KNC+KP-KK+1 ), 1, AP( KPC+1 ),
533 DO 80 J = KK + 1, KP - 1
535 T = DCONJG( AP( KNC+J-KK ) )
536 AP( KNC+J-KK ) = DCONJG( AP( KX ) )
539 AP( KNC+KP-KK ) = DCONJG( AP( KNC+KP-KK ) )
540 R1 = DBLE( AP( KNC ) )
541 AP( KNC ) = DBLE( AP( KPC ) )
543 IF( KSTEP.EQ.2 ) THEN
544 AP( KC ) = DBLE( AP( KC ) )
546 AP( KC+1 ) = AP( KC+KP-K )
550 AP( KC ) = DBLE( AP( KC ) )
552 $ AP( KNC ) = DBLE( AP( KNC ) )
555 * Update the trailing submatrix
557 IF( KSTEP.EQ.1 ) THEN
559 * 1-by-1 pivot block D(k): column k now holds
563 * where L(k) is the k-th column of L
567 * Perform a rank-1 update of A(k+1:n,k+1:n) as
569 * A := A - L(k)*D(k)*L(k)**H = A - W(k)*(1/D(k))*W(k)**H
571 R1 = ONE / DBLE( AP( KC ) )
572 CALL ZHPR( UPLO, N-K, -R1, AP( KC+1 ), 1,
575 * Store L(k) in column K
577 CALL ZDSCAL( N-K, R1, AP( KC+1 ), 1 )
581 * 2-by-2 pivot block D(k): columns K and K+1 now hold
583 * ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k)
585 * where L(k) and L(k+1) are the k-th and (k+1)-th columns
590 * Perform a rank-2 update of A(k+2:n,k+2:n) as
592 * A := A - ( L(k) L(k+1) )*D(k)*( L(k) L(k+1) )**H
593 * = A - ( W(k) W(k+1) )*inv(D(k))*( W(k) W(k+1) )**H
595 * where L(k) and L(k+1) are the k-th and (k+1)-th
598 D = DLAPY2( DBLE( AP( K+1+( K-1 )*( 2*N-K ) / 2 ) ),
599 $ DIMAG( AP( K+1+( K-1 )*( 2*N-K ) / 2 ) ) )
600 D11 = DBLE( AP( K+1+K*( 2*N-K-1 ) / 2 ) ) / D
601 D22 = DBLE( AP( K+( K-1 )*( 2*N-K ) / 2 ) ) / D
602 TT = ONE / ( D11*D22-ONE )
603 D21 = AP( K+1+( K-1 )*( 2*N-K ) / 2 ) / D
607 WK = D*( D11*AP( J+( K-1 )*( 2*N-K ) / 2 )-D21*
608 $ AP( J+K*( 2*N-K-1 ) / 2 ) )
609 WKP1 = D*( D22*AP( J+K*( 2*N-K-1 ) / 2 )-
610 $ DCONJG( D21 )*AP( J+( K-1 )*( 2*N-K ) /
613 AP( I+( J-1 )*( 2*N-J ) / 2 ) = AP( I+( J-1 )*
614 $ ( 2*N-J ) / 2 ) - AP( I+( K-1 )*( 2*N-K ) /
615 $ 2 )*DCONJG( WK ) - AP( I+K*( 2*N-K-1 ) / 2 )*
618 AP( J+( K-1 )*( 2*N-K ) / 2 ) = WK
619 AP( J+K*( 2*N-K-1 ) / 2 ) = WKP1
620 AP( J+( J-1 )*( 2*N-J ) / 2 )
621 $ = DCMPLX( DBLE( AP( J+( J-1 )*( 2*N-J ) / 2 ) ),
628 * Store details of the interchanges in IPIV
630 IF( KSTEP.EQ.1 ) THEN
637 * Increase K and return to the start of the main loop