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