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