LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ zhetf2_rk()

subroutine zhetf2_rk ( character  uplo,
integer  n,
complex*16, dimension( lda, * )  a,
integer  lda,
complex*16, dimension( * )  e,
integer, dimension( * )  ipiv,
integer  info 
)

ZHETF2_RK computes the factorization of a complex Hermitian indefinite matrix using the bounded Bunch-Kaufman (rook) diagonal pivoting method (BLAS2 unblocked algorithm).

Download ZHETF2_RK + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 ZHETF2_RK computes the factorization of a complex Hermitian matrix A
 using the bounded Bunch-Kaufman (rook) diagonal pivoting method:

    A = P*U*D*(U**H)*(P**T) or A = P*L*D*(L**H)*(P**T),

 where U (or L) is unit upper (or lower) triangular matrix,
 U**H (or L**H) is the conjugate of U (or L), P is a permutation
 matrix, P**T is the transpose of P, and D is Hermitian and block
 diagonal with 1-by-1 and 2-by-2 diagonal blocks.

 This is the unblocked version of the algorithm, calling Level 2 BLAS.
 For more information see Further Details section.
Parameters
[in]UPLO
          UPLO is CHARACTER*1
          Specifies whether the upper or lower triangular part of the
          Hermitian matrix A is stored:
          = 'U':  Upper triangular
          = 'L':  Lower triangular
[in]N
          N is INTEGER
          The order of the matrix A.  N >= 0.
[in,out]A
          A is COMPLEX*16 array, dimension (LDA,N)
          On entry, the Hermitian matrix A.
            If UPLO = 'U': the leading N-by-N upper triangular part
            of A contains the upper triangular part of the matrix A,
            and the strictly lower triangular part of A is not
            referenced.

            If UPLO = 'L': the leading N-by-N lower triangular part
            of A contains the lower triangular part of the matrix A,
            and the strictly upper triangular part of A is not
            referenced.

          On exit, contains:
            a) ONLY diagonal elements of the Hermitian block diagonal
               matrix D on the diagonal of A, i.e. D(k,k) = A(k,k);
               (superdiagonal (or subdiagonal) elements of D
                are stored on exit in array E), and
            b) If UPLO = 'U': factor U in the superdiagonal part of A.
               If UPLO = 'L': factor L in the subdiagonal part of A.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,N).
[out]E
          E is COMPLEX*16 array, dimension (N)
          On exit, contains the superdiagonal (or subdiagonal)
          elements of the Hermitian block diagonal matrix D
          with 1-by-1 or 2-by-2 diagonal blocks, where
          If UPLO = 'U': E(i) = D(i-1,i), i=2:N, E(1) is set to 0;
          If UPLO = 'L': E(i) = D(i+1,i), i=1:N-1, E(N) is set to 0.

          NOTE: For 1-by-1 diagonal block D(k), where
          1 <= k <= N, the element E(k) is set to 0 in both
          UPLO = 'U' or UPLO = 'L' cases.
[out]IPIV
          IPIV is INTEGER array, dimension (N)
          IPIV describes the permutation matrix P in the factorization
          of matrix A as follows. The absolute value of IPIV(k)
          represents the index of row and column that were
          interchanged with the k-th row and column. The value of UPLO
          describes the order in which the interchanges were applied.
          Also, the sign of IPIV represents the block structure of
          the Hermitian block diagonal matrix D with 1-by-1 or 2-by-2
          diagonal blocks which correspond to 1 or 2 interchanges
          at each factorization step. For more info see Further
          Details section.

          If UPLO = 'U',
          ( in factorization order, k decreases from N to 1 ):
            a) A single positive entry IPIV(k) > 0 means:
               D(k,k) is a 1-by-1 diagonal block.
               If IPIV(k) != k, rows and columns k and IPIV(k) were
               interchanged in the matrix A(1:N,1:N);
               If IPIV(k) = k, no interchange occurred.

            b) A pair of consecutive negative entries
               IPIV(k) < 0 and IPIV(k-1) < 0 means:
               D(k-1:k,k-1:k) is a 2-by-2 diagonal block.
               (NOTE: negative entries in IPIV appear ONLY in pairs).
               1) If -IPIV(k) != k, rows and columns
                  k and -IPIV(k) were interchanged
                  in the matrix A(1:N,1:N).
                  If -IPIV(k) = k, no interchange occurred.
               2) If -IPIV(k-1) != k-1, rows and columns
                  k-1 and -IPIV(k-1) were interchanged
                  in the matrix A(1:N,1:N).
                  If -IPIV(k-1) = k-1, no interchange occurred.

            c) In both cases a) and b), always ABS( IPIV(k) ) <= k.

            d) NOTE: Any entry IPIV(k) is always NONZERO on output.

          If UPLO = 'L',
          ( in factorization order, k increases from 1 to N ):
            a) A single positive entry IPIV(k) > 0 means:
               D(k,k) is a 1-by-1 diagonal block.
               If IPIV(k) != k, rows and columns k and IPIV(k) were
               interchanged in the matrix A(1:N,1:N).
               If IPIV(k) = k, no interchange occurred.

            b) A pair of consecutive negative entries
               IPIV(k) < 0 and IPIV(k+1) < 0 means:
               D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
               (NOTE: negative entries in IPIV appear ONLY in pairs).
               1) If -IPIV(k) != k, rows and columns
                  k and -IPIV(k) were interchanged
                  in the matrix A(1:N,1:N).
                  If -IPIV(k) = k, no interchange occurred.
               2) If -IPIV(k+1) != k+1, rows and columns
                  k-1 and -IPIV(k-1) were interchanged
                  in the matrix A(1:N,1:N).
                  If -IPIV(k+1) = k+1, no interchange occurred.

            c) In both cases a) and b), always ABS( IPIV(k) ) >= k.

            d) NOTE: Any entry IPIV(k) is always NONZERO on output.
[out]INFO
          INFO is INTEGER
          = 0: successful exit

          < 0: If INFO = -k, the k-th argument had an illegal value

          > 0: If INFO = k, the matrix A is singular, because:
                 If UPLO = 'U': column k in the upper
                 triangular part of A contains all zeros.
                 If UPLO = 'L': column k in the lower
                 triangular part of A contains all zeros.

               Therefore D(k,k) is exactly zero, and superdiagonal
               elements of column k of U (or subdiagonal elements of
               column k of L ) are all zeros. The factorization has
               been completed, but the block diagonal matrix D is
               exactly singular, and division by zero will occur if
               it is used to solve a system of equations.

               NOTE: INFO only stores the first occurrence of
               a singularity, any subsequent occurrence of singularity
               is not stored in INFO even though the factorization
               always completes.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
 TODO: put further details
Contributors:
  December 2016,  Igor Kozachenko,
                  Computer Science Division,
                  University of California, Berkeley

  September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas,
                  School of Mathematics,
                  University of Manchester

  01-01-96 - Based on modifications by
    J. Lewis, Boeing Computer Services Company
    A. Petitet, Computer Science Dept.,
                Univ. of Tenn., Knoxville abd , USA

Definition at line 240 of file zhetf2_rk.f.

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*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine zher(uplo, n, alpha, x, incx, a, lda)
ZHER
Definition zher.f:135
integer function izamax(n, zx, incx)
IZAMAX
Definition izamax.f:71
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
double precision function dlapy2(x, y)
DLAPY2 returns sqrt(x2+y2).
Definition dlapy2.f:63
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
subroutine zdscal(n, da, zx, incx)
ZDSCAL
Definition zdscal.f:78
subroutine zswap(n, zx, incx, zy, incy)
ZSWAP
Definition zswap.f:81
Here is the call graph for this function:
Here is the caller graph for this function: