LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ 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 interserction
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 interserction
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 interserction
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 interserction
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 *
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
integer function izamax(N, ZX, INCX)
IZAMAX
Definition: izamax.f:71
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
Definition: zswap.f:81
subroutine zdscal(N, DA, ZX, INCX)
ZDSCAL
Definition: zdscal.f:78
subroutine zher(UPLO, N, ALPHA, X, INCX, A, LDA)
ZHER
Definition: zher.f:135
Here is the call graph for this function:
Here is the caller graph for this function: