LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
zhetf2_rk.f
Go to the documentation of this file.
1*> \brief \b ZHETF2_RK computes the factorization of a complex Hermitian indefinite matrix using the bounded Bunch-Kaufman (rook) diagonal pivoting method (BLAS2 unblocked algorithm).
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download ZHETF2_RK + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhetf2_rk.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhetf2_rk.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhetf2_rk.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE ZHETF2_RK( UPLO, N, A, LDA, E, IPIV, INFO )
22*
23* .. Scalar Arguments ..
24* CHARACTER UPLO
25* INTEGER INFO, LDA, N
26* ..
27* .. Array Arguments ..
28* INTEGER IPIV( * )
29* COMPLEX*16 A( LDA, * ), E ( * )
30* ..
31*
32*
33*> \par Purpose:
34* =============
35*>
36*> \verbatim
37*> ZHETF2_RK computes the factorization of a complex Hermitian matrix A
38*> using the bounded Bunch-Kaufman (rook) diagonal pivoting method:
39*>
40*> A = P*U*D*(U**H)*(P**T) or A = P*L*D*(L**H)*(P**T),
41*>
42*> where U (or L) is unit upper (or lower) triangular matrix,
43*> U**H (or L**H) is the conjugate of U (or L), P is a permutation
44*> matrix, P**T is the transpose of P, and D is Hermitian and block
45*> diagonal with 1-by-1 and 2-by-2 diagonal blocks.
46*>
47*> This is the unblocked version of the algorithm, calling Level 2 BLAS.
48*> For more information see Further Details section.
49*> \endverbatim
50*
51* Arguments:
52* ==========
53*
54*> \param[in] UPLO
55*> \verbatim
56*> UPLO is CHARACTER*1
57*> Specifies whether the upper or lower triangular part of the
58*> Hermitian matrix A is stored:
59*> = 'U': Upper triangular
60*> = 'L': Lower triangular
61*> \endverbatim
62*>
63*> \param[in] N
64*> \verbatim
65*> N is INTEGER
66*> The order of the matrix A. N >= 0.
67*> \endverbatim
68*>
69*> \param[in,out] A
70*> \verbatim
71*> A is COMPLEX*16 array, dimension (LDA,N)
72*> On entry, the Hermitian matrix A.
73*> If UPLO = 'U': the leading N-by-N upper triangular part
74*> of A contains the upper triangular part of the matrix A,
75*> and the strictly lower triangular part of A is not
76*> referenced.
77*>
78*> If UPLO = 'L': the leading N-by-N lower triangular part
79*> of A contains the lower triangular part of the matrix A,
80*> and the strictly upper triangular part of A is not
81*> referenced.
82*>
83*> On exit, contains:
84*> a) ONLY diagonal elements of the Hermitian block diagonal
85*> matrix D on the diagonal of A, i.e. D(k,k) = A(k,k);
86*> (superdiagonal (or subdiagonal) elements of D
87*> are stored on exit in array E), and
88*> b) If UPLO = 'U': factor U in the superdiagonal part of A.
89*> If UPLO = 'L': factor L in the subdiagonal part of A.
90*> \endverbatim
91*>
92*> \param[in] LDA
93*> \verbatim
94*> LDA is INTEGER
95*> The leading dimension of the array A. LDA >= max(1,N).
96*> \endverbatim
97*>
98*> \param[out] E
99*> \verbatim
100*> E is COMPLEX*16 array, dimension (N)
101*> On exit, contains the superdiagonal (or subdiagonal)
102*> elements of the Hermitian block diagonal matrix D
103*> with 1-by-1 or 2-by-2 diagonal blocks, where
104*> If UPLO = 'U': E(i) = D(i-1,i), i=2:N, E(1) is set to 0;
105*> If UPLO = 'L': E(i) = D(i+1,i), i=1:N-1, E(N) is set to 0.
106*>
107*> NOTE: For 1-by-1 diagonal block D(k), where
108*> 1 <= k <= N, the element E(k) is set to 0 in both
109*> UPLO = 'U' or UPLO = 'L' cases.
110*> \endverbatim
111*>
112*> \param[out] IPIV
113*> \verbatim
114*> IPIV is INTEGER array, dimension (N)
115*> IPIV describes the permutation matrix P in the factorization
116*> of matrix A as follows. The absolute value of IPIV(k)
117*> represents the index of row and column that were
118*> interchanged with the k-th row and column. The value of UPLO
119*> describes the order in which the interchanges were applied.
120*> Also, the sign of IPIV represents the block structure of
121*> the Hermitian block diagonal matrix D with 1-by-1 or 2-by-2
122*> diagonal blocks which correspond to 1 or 2 interchanges
123*> at each factorization step. For more info see Further
124*> Details section.
125*>
126*> If UPLO = 'U',
127*> ( in factorization order, k decreases from N to 1 ):
128*> a) A single positive entry IPIV(k) > 0 means:
129*> D(k,k) is a 1-by-1 diagonal block.
130*> If IPIV(k) != k, rows and columns k and IPIV(k) were
131*> interchanged in the matrix A(1:N,1:N);
132*> If IPIV(k) = k, no interchange occurred.
133*>
134*> b) A pair of consecutive negative entries
135*> IPIV(k) < 0 and IPIV(k-1) < 0 means:
136*> D(k-1:k,k-1:k) is a 2-by-2 diagonal block.
137*> (NOTE: negative entries in IPIV appear ONLY in pairs).
138*> 1) If -IPIV(k) != k, rows and columns
139*> k and -IPIV(k) were interchanged
140*> in the matrix A(1:N,1:N).
141*> If -IPIV(k) = k, no interchange occurred.
142*> 2) If -IPIV(k-1) != k-1, rows and columns
143*> k-1 and -IPIV(k-1) were interchanged
144*> in the matrix A(1:N,1:N).
145*> If -IPIV(k-1) = k-1, no interchange occurred.
146*>
147*> c) In both cases a) and b), always ABS( IPIV(k) ) <= k.
148*>
149*> d) NOTE: Any entry IPIV(k) is always NONZERO on output.
150*>
151*> If UPLO = 'L',
152*> ( in factorization order, k increases from 1 to N ):
153*> a) A single positive entry IPIV(k) > 0 means:
154*> D(k,k) is a 1-by-1 diagonal block.
155*> If IPIV(k) != k, rows and columns k and IPIV(k) were
156*> interchanged in the matrix A(1:N,1:N).
157*> If IPIV(k) = k, no interchange occurred.
158*>
159*> b) A pair of consecutive negative entries
160*> IPIV(k) < 0 and IPIV(k+1) < 0 means:
161*> D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
162*> (NOTE: negative entries in IPIV appear ONLY in pairs).
163*> 1) If -IPIV(k) != k, rows and columns
164*> k and -IPIV(k) were interchanged
165*> in the matrix A(1:N,1:N).
166*> If -IPIV(k) = k, no interchange occurred.
167*> 2) If -IPIV(k+1) != k+1, rows and columns
168*> k-1 and -IPIV(k-1) were interchanged
169*> in the matrix A(1:N,1:N).
170*> If -IPIV(k+1) = k+1, no interchange occurred.
171*>
172*> c) In both cases a) and b), always ABS( IPIV(k) ) >= k.
173*>
174*> d) NOTE: Any entry IPIV(k) is always NONZERO on output.
175*> \endverbatim
176*>
177*> \param[out] INFO
178*> \verbatim
179*> INFO is INTEGER
180*> = 0: successful exit
181*>
182*> < 0: If INFO = -k, the k-th argument had an illegal value
183*>
184*> > 0: If INFO = k, the matrix A is singular, because:
185*> If UPLO = 'U': column k in the upper
186*> triangular part of A contains all zeros.
187*> If UPLO = 'L': column k in the lower
188*> triangular part of A contains all zeros.
189*>
190*> Therefore D(k,k) is exactly zero, and superdiagonal
191*> elements of column k of U (or subdiagonal elements of
192*> column k of L ) are all zeros. The factorization has
193*> been completed, but the block diagonal matrix D is
194*> exactly singular, and division by zero will occur if
195*> it is used to solve a system of equations.
196*>
197*> NOTE: INFO only stores the first occurrence of
198*> a singularity, any subsequent occurrence of singularity
199*> is not stored in INFO even though the factorization
200*> always completes.
201*> \endverbatim
202*
203* Authors:
204* ========
205*
206*> \author Univ. of Tennessee
207*> \author Univ. of California Berkeley
208*> \author Univ. of Colorado Denver
209*> \author NAG Ltd.
210*
211*> \ingroup hetf2_rk
212*
213*> \par Further Details:
214* =====================
215*>
216*> \verbatim
217*> TODO: put further details
218*> \endverbatim
219*
220*> \par Contributors:
221* ==================
222*>
223*> \verbatim
224*>
225*> December 2016, Igor Kozachenko,
226*> Computer Science Division,
227*> University of California, Berkeley
228*>
229*> September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas,
230*> School of Mathematics,
231*> University of Manchester
232*>
233*> 01-01-96 - Based on modifications by
234*> J. Lewis, Boeing Computer Services Company
235*> A. Petitet, Computer Science Dept.,
236*> Univ. of Tenn., Knoxville abd , USA
237*> \endverbatim
238*
239* =====================================================================
240 SUBROUTINE zhetf2_rk( UPLO, N, A, LDA, E, IPIV, INFO )
241*
242* -- LAPACK computational routine --
243* -- LAPACK is a software package provided by Univ. of Tennessee, --
244* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
245*
246* .. Scalar Arguments ..
247 CHARACTER UPLO
248 INTEGER INFO, LDA, N
249* ..
250* .. Array Arguments ..
251 INTEGER IPIV( * )
252 COMPLEX*16 A( LDA, * ), E( * )
253* ..
254*
255* ======================================================================
256*
257* .. Parameters ..
258 DOUBLE PRECISION ZERO, ONE
259 parameter( zero = 0.0d+0, one = 1.0d+0 )
260 DOUBLE PRECISION EIGHT, SEVTEN
261 parameter( eight = 8.0d+0, sevten = 17.0d+0 )
262 COMPLEX*16 CZERO
263 parameter( czero = ( 0.0d+0, 0.0d+0 ) )
264* ..
265* .. Local Scalars ..
266 LOGICAL DONE, UPPER
267 INTEGER I, II, IMAX, ITEMP, J, JMAX, K, KK, KP, KSTEP,
268 $ P
269 DOUBLE PRECISION ABSAKK, ALPHA, COLMAX, D, D11, D22, R1, DTEMP,
270 $ ROWMAX, TT, SFMIN
271 COMPLEX*16 D12, D21, T, WK, WKM1, WKP1, Z
272* ..
273* .. External Functions ..
274*
275 LOGICAL LSAME
276 INTEGER IZAMAX
277 DOUBLE PRECISION DLAMCH, DLAPY2
278 EXTERNAL lsame, izamax, dlamch, dlapy2
279* ..
280* .. External Subroutines ..
281 EXTERNAL xerbla, zdscal, zher, zswap
282* ..
283* .. Intrinsic Functions ..
284 INTRINSIC abs, dble, dcmplx, dconjg, dimag, max, sqrt
285* ..
286* .. Statement Functions ..
287 DOUBLE PRECISION CABS1
288* ..
289* .. Statement Function definitions ..
290 cabs1( z ) = abs( dble( z ) ) + abs( dimag( z ) )
291* ..
292* .. Executable Statements ..
293*
294* Test the input parameters.
295*
296 info = 0
297 upper = lsame( uplo, 'U' )
298 IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
299 info = -1
300 ELSE IF( n.LT.0 ) THEN
301 info = -2
302 ELSE IF( lda.LT.max( 1, n ) ) THEN
303 info = -4
304 END IF
305 IF( info.NE.0 ) THEN
306 CALL xerbla( 'ZHETF2_RK', -info )
307 RETURN
308 END IF
309*
310* Initialize ALPHA for use in choosing pivot block size.
311*
312 alpha = ( one+sqrt( sevten ) ) / eight
313*
314* Compute machine safe minimum
315*
316 sfmin = dlamch( 'S' )
317*
318 IF( upper ) THEN
319*
320* Factorize A as U*D*U**H using the upper triangle of A
321*
322* Initialize the first entry of array E, where superdiagonal
323* elements of D are stored
324*
325 e( 1 ) = czero
326*
327* K is the main loop index, decreasing from N to 1 in steps of
328* 1 or 2
329*
330 k = n
331 10 CONTINUE
332*
333* If K < 1, exit from loop
334*
335 IF( k.LT.1 )
336 $ GO TO 34
337 kstep = 1
338 p = k
339*
340* Determine rows and columns to be interchanged and whether
341* a 1-by-1 or 2-by-2 pivot block will be used
342*
343 absakk = abs( dble( a( k, k ) ) )
344*
345* IMAX is the row-index of the largest off-diagonal element in
346* column K, and COLMAX is its absolute value.
347* Determine both COLMAX and IMAX.
348*
349 IF( k.GT.1 ) THEN
350 imax = izamax( k-1, a( 1, k ), 1 )
351 colmax = cabs1( a( imax, k ) )
352 ELSE
353 colmax = zero
354 END IF
355*
356 IF( ( max( absakk, colmax ).EQ.zero ) ) THEN
357*
358* Column K is zero or underflow: set INFO and continue
359*
360 IF( info.EQ.0 )
361 $ info = k
362 kp = k
363 a( k, k ) = dble( a( k, k ) )
364*
365* Set E( K ) to zero
366*
367 IF( k.GT.1 )
368 $ e( k ) = czero
369*
370 ELSE
371*
372* ============================================================
373*
374* BEGIN pivot search
375*
376* Case(1)
377* Equivalent to testing for ABSAKK.GE.ALPHA*COLMAX
378* (used to handle NaN and Inf)
379*
380 IF( .NOT.( absakk.LT.alpha*colmax ) ) THEN
381*
382* no interchange, use 1-by-1 pivot block
383*
384 kp = k
385*
386 ELSE
387*
388 done = .false.
389*
390* Loop until pivot found
391*
392 12 CONTINUE
393*
394* BEGIN pivot search loop body
395*
396*
397* JMAX is the column-index of the largest off-diagonal
398* element in row IMAX, and ROWMAX is its absolute value.
399* Determine both ROWMAX and JMAX.
400*
401 IF( imax.NE.k ) THEN
402 jmax = imax + izamax( k-imax, a( imax, imax+1 ),
403 $ lda )
404 rowmax = cabs1( a( imax, jmax ) )
405 ELSE
406 rowmax = zero
407 END IF
408*
409 IF( imax.GT.1 ) THEN
410 itemp = izamax( imax-1, a( 1, imax ), 1 )
411 dtemp = cabs1( a( itemp, imax ) )
412 IF( dtemp.GT.rowmax ) THEN
413 rowmax = dtemp
414 jmax = itemp
415 END IF
416 END IF
417*
418* Case(2)
419* Equivalent to testing for
420* ABS( DBLE( W( IMAX,KW-1 ) ) ).GE.ALPHA*ROWMAX
421* (used to handle NaN and Inf)
422*
423 IF( .NOT.( abs( dble( a( imax, imax ) ) )
424 $ .LT.alpha*rowmax ) ) THEN
425*
426* interchange rows and columns K and IMAX,
427* use 1-by-1 pivot block
428*
429 kp = imax
430 done = .true.
431*
432* Case(3)
433* Equivalent to testing for ROWMAX.EQ.COLMAX,
434* (used to handle NaN and Inf)
435*
436 ELSE IF( ( p.EQ.jmax ) .OR. ( rowmax.LE.colmax ) )
437 $ THEN
438*
439* interchange rows and columns K-1 and IMAX,
440* use 2-by-2 pivot block
441*
442 kp = imax
443 kstep = 2
444 done = .true.
445*
446* Case(4)
447 ELSE
448*
449* Pivot not found: set params and repeat
450*
451 p = imax
452 colmax = rowmax
453 imax = jmax
454 END IF
455*
456* END pivot search loop body
457*
458 IF( .NOT.done ) GOTO 12
459*
460 END IF
461*
462* END pivot search
463*
464* ============================================================
465*
466* KK is the column of A where pivoting step stopped
467*
468 kk = k - kstep + 1
469*
470* For only a 2x2 pivot, interchange rows and columns K and P
471* in the leading submatrix A(1:k,1:k)
472*
473 IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) ) THEN
474* (1) Swap columnar parts
475 IF( p.GT.1 )
476 $ CALL zswap( p-1, a( 1, k ), 1, a( 1, p ), 1 )
477* (2) Swap and conjugate middle parts
478 DO 14 j = p + 1, k - 1
479 t = dconjg( a( j, k ) )
480 a( j, k ) = dconjg( a( p, j ) )
481 a( p, j ) = t
482 14 CONTINUE
483* (3) Swap and conjugate corner elements at row-col intersection
484 a( p, k ) = dconjg( a( p, k ) )
485* (4) Swap diagonal elements at row-col intersection
486 r1 = dble( a( k, k ) )
487 a( k, k ) = dble( a( p, p ) )
488 a( p, p ) = r1
489*
490* Convert upper triangle of A into U form by applying
491* the interchanges in columns k+1:N.
492*
493 IF( k.LT.n )
494 $ CALL zswap( n-k, a( k, k+1 ), lda, a( p, k+1 ), lda )
495*
496 END IF
497*
498* For both 1x1 and 2x2 pivots, interchange rows and
499* columns KK and KP in the leading submatrix A(1:k,1:k)
500*
501 IF( kp.NE.kk ) THEN
502* (1) Swap columnar parts
503 IF( kp.GT.1 )
504 $ CALL zswap( kp-1, a( 1, kk ), 1, a( 1, kp ), 1 )
505* (2) Swap and conjugate middle parts
506 DO 15 j = kp + 1, kk - 1
507 t = dconjg( a( j, kk ) )
508 a( j, kk ) = dconjg( a( kp, j ) )
509 a( kp, j ) = t
510 15 CONTINUE
511* (3) Swap and conjugate corner elements at row-col intersection
512 a( kp, kk ) = dconjg( a( kp, kk ) )
513* (4) Swap diagonal elements at row-col intersection
514 r1 = dble( a( kk, kk ) )
515 a( kk, kk ) = dble( a( kp, kp ) )
516 a( kp, kp ) = r1
517*
518 IF( kstep.EQ.2 ) THEN
519* (*) Make sure that diagonal element of pivot is real
520 a( k, k ) = dble( a( k, k ) )
521* (5) Swap row elements
522 t = a( k-1, k )
523 a( k-1, k ) = a( kp, k )
524 a( kp, k ) = t
525 END IF
526*
527* Convert upper triangle of A into U form by applying
528* the interchanges in columns k+1:N.
529*
530 IF( k.LT.n )
531 $ CALL zswap( n-k, a( kk, k+1 ), lda, a( kp, k+1 ),
532 $ lda )
533*
534 ELSE
535* (*) Make sure that diagonal element of pivot is real
536 a( k, k ) = dble( a( k, k ) )
537 IF( kstep.EQ.2 )
538 $ a( k-1, k-1 ) = dble( a( k-1, k-1 ) )
539 END IF
540*
541* Update the leading submatrix
542*
543 IF( kstep.EQ.1 ) THEN
544*
545* 1-by-1 pivot block D(k): column k now holds
546*
547* W(k) = U(k)*D(k)
548*
549* where U(k) is the k-th column of U
550*
551 IF( k.GT.1 ) THEN
552*
553* Perform a rank-1 update of A(1:k-1,1:k-1) and
554* store U(k) in column k
555*
556 IF( abs( dble( a( k, k ) ) ).GE.sfmin ) THEN
557*
558* Perform a rank-1 update of A(1:k-1,1:k-1) as
559* A := A - U(k)*D(k)*U(k)**T
560* = A - W(k)*1/D(k)*W(k)**T
561*
562 d11 = one / dble( a( k, k ) )
563 CALL zher( uplo, k-1, -d11, a( 1, k ), 1, a, lda )
564*
565* Store U(k) in column k
566*
567 CALL zdscal( k-1, d11, a( 1, k ), 1 )
568 ELSE
569*
570* Store L(k) in column K
571*
572 d11 = dble( a( k, k ) )
573 DO 16 ii = 1, k - 1
574 a( ii, k ) = a( ii, k ) / d11
575 16 CONTINUE
576*
577* Perform a rank-1 update of A(k+1:n,k+1:n) as
578* A := A - U(k)*D(k)*U(k)**T
579* = A - W(k)*(1/D(k))*W(k)**T
580* = A - (W(k)/D(k))*(D(k))*(W(k)/D(K))**T
581*
582 CALL zher( uplo, k-1, -d11, a( 1, k ), 1, a, lda )
583 END IF
584*
585* Store the superdiagonal element of D in array E
586*
587 e( k ) = czero
588*
589 END IF
590*
591 ELSE
592*
593* 2-by-2 pivot block D(k): columns k and k-1 now hold
594*
595* ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
596*
597* where U(k) and U(k-1) are the k-th and (k-1)-th columns
598* of U
599*
600* Perform a rank-2 update of A(1:k-2,1:k-2) as
601*
602* A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )**T
603* = A - ( ( A(k-1)A(k) )*inv(D(k)) ) * ( A(k-1)A(k) )**T
604*
605* and store L(k) and L(k+1) in columns k and k+1
606*
607 IF( k.GT.2 ) THEN
608* D = |A12|
609 d = dlapy2( dble( a( k-1, k ) ),
610 $ dimag( a( k-1, k ) ) )
611 d11 = dble( a( k, k ) / d )
612 d22 = dble( a( k-1, k-1 ) / d )
613 d12 = a( k-1, k ) / d
614 tt = one / ( d11*d22-one )
615*
616 DO 30 j = k - 2, 1, -1
617*
618* Compute D21 * ( W(k)W(k+1) ) * inv(D(k)) for row J
619*
620 wkm1 = tt*( d11*a( j, k-1 )-dconjg( d12 )*
621 $ a( j, k ) )
622 wk = tt*( d22*a( j, k )-d12*a( j, k-1 ) )
623*
624* Perform a rank-2 update of A(1:k-2,1:k-2)
625*
626 DO 20 i = j, 1, -1
627 a( i, j ) = a( i, j ) -
628 $ ( a( i, k ) / d )*dconjg( wk ) -
629 $ ( a( i, k-1 ) / d )*dconjg( wkm1 )
630 20 CONTINUE
631*
632* Store U(k) and U(k-1) in cols k and k-1 for row J
633*
634 a( j, k ) = wk / d
635 a( j, k-1 ) = wkm1 / d
636* (*) Make sure that diagonal element of pivot is real
637 a( j, j ) = dcmplx( dble( a( j, j ) ), zero )
638*
639 30 CONTINUE
640*
641 END IF
642*
643* Copy superdiagonal elements of D(K) to E(K) and
644* ZERO out superdiagonal entry of A
645*
646 e( k ) = a( k-1, k )
647 e( k-1 ) = czero
648 a( k-1, k ) = czero
649*
650 END IF
651*
652* End column K is nonsingular
653*
654 END IF
655*
656* Store details of the interchanges in IPIV
657*
658 IF( kstep.EQ.1 ) THEN
659 ipiv( k ) = kp
660 ELSE
661 ipiv( k ) = -p
662 ipiv( k-1 ) = -kp
663 END IF
664*
665* Decrease K and return to the start of the main loop
666*
667 k = k - kstep
668 GO TO 10
669*
670 34 CONTINUE
671*
672 ELSE
673*
674* Factorize A as L*D*L**H using the lower triangle of A
675*
676* Initialize the unused last entry of the subdiagonal array E.
677*
678 e( n ) = czero
679*
680* K is the main loop index, increasing from 1 to N in steps of
681* 1 or 2
682*
683 k = 1
684 40 CONTINUE
685*
686* If K > N, exit from loop
687*
688 IF( k.GT.n )
689 $ GO TO 64
690 kstep = 1
691 p = k
692*
693* Determine rows and columns to be interchanged and whether
694* a 1-by-1 or 2-by-2 pivot block will be used
695*
696 absakk = abs( dble( a( k, k ) ) )
697*
698* IMAX is the row-index of the largest off-diagonal element in
699* column K, and COLMAX is its absolute value.
700* Determine both COLMAX and IMAX.
701*
702 IF( k.LT.n ) THEN
703 imax = k + izamax( n-k, a( k+1, k ), 1 )
704 colmax = cabs1( a( imax, k ) )
705 ELSE
706 colmax = zero
707 END IF
708*
709 IF( max( absakk, colmax ).EQ.zero ) THEN
710*
711* Column K is zero or underflow: set INFO and continue
712*
713 IF( info.EQ.0 )
714 $ info = k
715 kp = k
716 a( k, k ) = dble( a( k, k ) )
717*
718* Set E( K ) to zero
719*
720 IF( k.LT.n )
721 $ e( k ) = czero
722*
723 ELSE
724*
725* ============================================================
726*
727* BEGIN pivot search
728*
729* Case(1)
730* Equivalent to testing for ABSAKK.GE.ALPHA*COLMAX
731* (used to handle NaN and Inf)
732*
733 IF( .NOT.( absakk.LT.alpha*colmax ) ) THEN
734*
735* no interchange, use 1-by-1 pivot block
736*
737 kp = k
738*
739 ELSE
740*
741 done = .false.
742*
743* Loop until pivot found
744*
745 42 CONTINUE
746*
747* BEGIN pivot search loop body
748*
749*
750* JMAX is the column-index of the largest off-diagonal
751* element in row IMAX, and ROWMAX is its absolute value.
752* Determine both ROWMAX and JMAX.
753*
754 IF( imax.NE.k ) THEN
755 jmax = k - 1 + izamax( imax-k, a( imax, k ), lda )
756 rowmax = cabs1( a( imax, jmax ) )
757 ELSE
758 rowmax = zero
759 END IF
760*
761 IF( imax.LT.n ) THEN
762 itemp = imax + izamax( n-imax, a( imax+1, imax ),
763 $ 1 )
764 dtemp = cabs1( a( itemp, imax ) )
765 IF( dtemp.GT.rowmax ) THEN
766 rowmax = dtemp
767 jmax = itemp
768 END IF
769 END IF
770*
771* Case(2)
772* Equivalent to testing for
773* ABS( DBLE( W( IMAX,KW-1 ) ) ).GE.ALPHA*ROWMAX
774* (used to handle NaN and Inf)
775*
776 IF( .NOT.( abs( dble( a( imax, imax ) ) )
777 $ .LT.alpha*rowmax ) ) THEN
778*
779* interchange rows and columns K and IMAX,
780* use 1-by-1 pivot block
781*
782 kp = imax
783 done = .true.
784*
785* Case(3)
786* Equivalent to testing for ROWMAX.EQ.COLMAX,
787* (used to handle NaN and Inf)
788*
789 ELSE IF( ( p.EQ.jmax ) .OR. ( rowmax.LE.colmax ) )
790 $ THEN
791*
792* interchange rows and columns K+1 and IMAX,
793* use 2-by-2 pivot block
794*
795 kp = imax
796 kstep = 2
797 done = .true.
798*
799* Case(4)
800 ELSE
801*
802* Pivot not found: set params and repeat
803*
804 p = imax
805 colmax = rowmax
806 imax = jmax
807 END IF
808*
809*
810* END pivot search loop body
811*
812 IF( .NOT.done ) GOTO 42
813*
814 END IF
815*
816* END pivot search
817*
818* ============================================================
819*
820* KK is the column of A where pivoting step stopped
821*
822 kk = k + kstep - 1
823*
824* For only a 2x2 pivot, interchange rows and columns K and P
825* in the trailing submatrix A(k:n,k:n)
826*
827 IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) ) THEN
828* (1) Swap columnar parts
829 IF( p.LT.n )
830 $ CALL zswap( n-p, a( p+1, k ), 1, a( p+1, p ), 1 )
831* (2) Swap and conjugate middle parts
832 DO 44 j = k + 1, p - 1
833 t = dconjg( a( j, k ) )
834 a( j, k ) = dconjg( a( p, j ) )
835 a( p, j ) = t
836 44 CONTINUE
837* (3) Swap and conjugate corner elements at row-col intersection
838 a( p, k ) = dconjg( a( p, k ) )
839* (4) Swap diagonal elements at row-col intersection
840 r1 = dble( a( k, k ) )
841 a( k, k ) = dble( a( p, p ) )
842 a( p, p ) = r1
843*
844* Convert lower triangle of A into L form by applying
845* the interchanges in columns 1:k-1.
846*
847 IF ( k.GT.1 )
848 $ CALL zswap( k-1, a( k, 1 ), lda, a( p, 1 ), lda )
849*
850 END IF
851*
852* For both 1x1 and 2x2 pivots, interchange rows and
853* columns KK and KP in the trailing submatrix A(k:n,k:n)
854*
855 IF( kp.NE.kk ) THEN
856* (1) Swap columnar parts
857 IF( kp.LT.n )
858 $ CALL zswap( n-kp, a( kp+1, kk ), 1, a( kp+1, kp ), 1 )
859* (2) Swap and conjugate middle parts
860 DO 45 j = kk + 1, kp - 1
861 t = dconjg( a( j, kk ) )
862 a( j, kk ) = dconjg( a( kp, j ) )
863 a( kp, j ) = t
864 45 CONTINUE
865* (3) Swap and conjugate corner elements at row-col intersection
866 a( kp, kk ) = dconjg( a( kp, kk ) )
867* (4) Swap diagonal elements at row-col intersection
868 r1 = dble( a( kk, kk ) )
869 a( kk, kk ) = dble( a( kp, kp ) )
870 a( kp, kp ) = r1
871*
872 IF( kstep.EQ.2 ) THEN
873* (*) Make sure that diagonal element of pivot is real
874 a( k, k ) = dble( a( k, k ) )
875* (5) Swap row elements
876 t = a( k+1, k )
877 a( k+1, k ) = a( kp, k )
878 a( kp, k ) = t
879 END IF
880*
881* Convert lower triangle of A into L form by applying
882* the interchanges in columns 1:k-1.
883*
884 IF ( k.GT.1 )
885 $ CALL zswap( k-1, a( kk, 1 ), lda, a( kp, 1 ), lda )
886*
887 ELSE
888* (*) Make sure that diagonal element of pivot is real
889 a( k, k ) = dble( a( k, k ) )
890 IF( kstep.EQ.2 )
891 $ a( k+1, k+1 ) = dble( a( k+1, k+1 ) )
892 END IF
893*
894* Update the trailing submatrix
895*
896 IF( kstep.EQ.1 ) THEN
897*
898* 1-by-1 pivot block D(k): column k of A now holds
899*
900* W(k) = L(k)*D(k),
901*
902* where L(k) is the k-th column of L
903*
904 IF( k.LT.n ) THEN
905*
906* Perform a rank-1 update of A(k+1:n,k+1:n) and
907* store L(k) in column k
908*
909* Handle division by a small number
910*
911 IF( abs( dble( a( k, k ) ) ).GE.sfmin ) THEN
912*
913* Perform a rank-1 update of A(k+1:n,k+1:n) as
914* A := A - L(k)*D(k)*L(k)**T
915* = A - W(k)*(1/D(k))*W(k)**T
916*
917 d11 = one / dble( a( k, k ) )
918 CALL zher( uplo, n-k, -d11, a( k+1, k ), 1,
919 $ a( k+1, k+1 ), lda )
920*
921* Store L(k) in column k
922*
923 CALL zdscal( n-k, d11, a( k+1, k ), 1 )
924 ELSE
925*
926* Store L(k) in column k
927*
928 d11 = dble( a( k, k ) )
929 DO 46 ii = k + 1, n
930 a( ii, k ) = a( ii, k ) / d11
931 46 CONTINUE
932*
933* Perform a rank-1 update of A(k+1:n,k+1:n) as
934* A := A - L(k)*D(k)*L(k)**T
935* = A - W(k)*(1/D(k))*W(k)**T
936* = A - (W(k)/D(k))*(D(k))*(W(k)/D(K))**T
937*
938 CALL zher( uplo, n-k, -d11, a( k+1, k ), 1,
939 $ a( k+1, k+1 ), lda )
940 END IF
941*
942* Store the subdiagonal element of D in array E
943*
944 e( k ) = czero
945*
946 END IF
947*
948 ELSE
949*
950* 2-by-2 pivot block D(k): columns k and k+1 now hold
951*
952* ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k)
953*
954* where L(k) and L(k+1) are the k-th and (k+1)-th columns
955* of L
956*
957*
958* Perform a rank-2 update of A(k+2:n,k+2:n) as
959*
960* A := A - ( L(k) L(k+1) ) * D(k) * ( L(k) L(k+1) )**T
961* = A - ( ( A(k)A(k+1) )*inv(D(k) ) * ( A(k)A(k+1) )**T
962*
963* and store L(k) and L(k+1) in columns k and k+1
964*
965 IF( k.LT.n-1 ) THEN
966* D = |A21|
967 d = dlapy2( dble( a( k+1, k ) ),
968 $ dimag( a( k+1, k ) ) )
969 d11 = dble( a( k+1, k+1 ) ) / d
970 d22 = dble( a( k, k ) ) / d
971 d21 = a( k+1, k ) / d
972 tt = one / ( d11*d22-one )
973*
974 DO 60 j = k + 2, n
975*
976* Compute D21 * ( W(k)W(k+1) ) * inv(D(k)) for row J
977*
978 wk = tt*( d11*a( j, k )-d21*a( j, k+1 ) )
979 wkp1 = tt*( d22*a( j, k+1 )-dconjg( d21 )*
980 $ a( j, k ) )
981*
982* Perform a rank-2 update of A(k+2:n,k+2:n)
983*
984 DO 50 i = j, n
985 a( i, j ) = a( i, j ) -
986 $ ( a( i, k ) / d )*dconjg( wk ) -
987 $ ( a( i, k+1 ) / d )*dconjg( wkp1 )
988 50 CONTINUE
989*
990* Store L(k) and L(k+1) in cols k and k+1 for row J
991*
992 a( j, k ) = wk / d
993 a( j, k+1 ) = wkp1 / d
994* (*) Make sure that diagonal element of pivot is real
995 a( j, j ) = dcmplx( dble( a( j, j ) ), zero )
996*
997 60 CONTINUE
998*
999 END IF
1000*
1001* Copy subdiagonal elements of D(K) to E(K) and
1002* ZERO out subdiagonal entry of A
1003*
1004 e( k ) = a( k+1, k )
1005 e( k+1 ) = czero
1006 a( k+1, k ) = czero
1007*
1008 END IF
1009*
1010* End column K is nonsingular
1011*
1012 END IF
1013*
1014* Store details of the interchanges in IPIV
1015*
1016 IF( kstep.EQ.1 ) THEN
1017 ipiv( k ) = kp
1018 ELSE
1019 ipiv( k ) = -p
1020 ipiv( k+1 ) = -kp
1021 END IF
1022*
1023* Increase K and return to the start of the main loop
1024*
1025 k = k + kstep
1026 GO TO 40
1027*
1028 64 CONTINUE
1029*
1030 END IF
1031*
1032 RETURN
1033*
1034* End of ZHETF2_RK
1035*
1036 END
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
zher
subroutine zher(uplo, n, alpha, x, incx, a, lda)
ZHER
Definition zher.f:135
zhetf2_rk
subroutine zhetf2_rk(uplo, n, a, lda, e, ipiv, info)
ZHETF2_RK computes the factorization of a complex Hermitian indefinite matrix using the bounded Bunch...
Definition zhetf2_rk.f:241
zdscal
subroutine zdscal(n, da, zx, incx)
ZDSCAL
Definition zdscal.f:78
zswap
subroutine zswap(n, zx, incx, zy, incy)
ZSWAP
Definition zswap.f:81