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