LAPACK  3.7.0
LAPACK: Linear Algebra PACKage
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.
Date
December 2016
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 243 of file zhetf2_rk.f.

243 *
244 * -- LAPACK computational routine (version 3.7.0) --
245 * -- LAPACK is a software package provided by Univ. of Tennessee, --
246 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
247 * December 2016
248 *
249 * .. Scalar Arguments ..
250  CHARACTER uplo
251  INTEGER info, lda, n
252 * ..
253 * .. Array Arguments ..
254  INTEGER ipiv( * )
255  COMPLEX*16 a( lda, * ), e( * )
256 * ..
257 *
258 * ======================================================================
259 *
260 * .. Parameters ..
261  DOUBLE PRECISION zero, one
262  parameter ( zero = 0.0d+0, one = 1.0d+0 )
263  DOUBLE PRECISION eight, sevten
264  parameter ( eight = 8.0d+0, sevten = 17.0d+0 )
265  COMPLEX*16 czero
266  parameter ( czero = ( 0.0d+0, 0.0d+0 ) )
267 * ..
268 * .. Local Scalars ..
269  LOGICAL done, upper
270  INTEGER i, ii, imax, itemp, j, jmax, k, kk, kp, kstep,
271  $ p
272  DOUBLE PRECISION absakk, alpha, colmax, d, d11, d22, r1, dtemp,
273  $ rowmax, tt, sfmin
274  COMPLEX*16 d12, d21, t, wk, wkm1, wkp1, z
275 * ..
276 * .. External Functions ..
277 *
278  LOGICAL lsame
279  INTEGER izamax
280  DOUBLE PRECISION dlamch, dlapy2
281  EXTERNAL lsame, izamax, dlamch, dlapy2
282 * ..
283 * .. External Subroutines ..
284  EXTERNAL xerbla, zdscal, zher, zswap
285 * ..
286 * .. Intrinsic Functions ..
287  INTRINSIC abs, dble, dcmplx, dconjg, dimag, max, sqrt
288 * ..
289 * .. Statement Functions ..
290  DOUBLE PRECISION cabs1
291 * ..
292 * .. Statement Function definitions ..
293  cabs1( z ) = abs( dble( z ) ) + abs( dimag( z ) )
294 * ..
295 * .. Executable Statements ..
296 *
297 * Test the input parameters.
298 *
299  info = 0
300  upper = lsame( uplo, 'U' )
301  IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
302  info = -1
303  ELSE IF( n.LT.0 ) THEN
304  info = -2
305  ELSE IF( lda.LT.max( 1, n ) ) THEN
306  info = -4
307  END IF
308  IF( info.NE.0 ) THEN
309  CALL xerbla( 'ZHETF2_RK', -info )
310  RETURN
311  END IF
312 *
313 * Initialize ALPHA for use in choosing pivot block size.
314 *
315  alpha = ( one+sqrt( sevten ) ) / eight
316 *
317 * Compute machine safe minimum
318 *
319  sfmin = dlamch( 'S' )
320 *
321  IF( upper ) THEN
322 *
323 * Factorize A as U*D*U**H using the upper triangle of A
324 *
325 * Initilize the first entry of array E, where superdiagonal
326 * elements of D are stored
327 *
328  e( 1 ) = czero
329 *
330 * K is the main loop index, decreasing from N to 1 in steps of
331 * 1 or 2
332 *
333  k = n
334  10 CONTINUE
335 *
336 * If K < 1, exit from loop
337 *
338  IF( k.LT.1 )
339  $ GO TO 34
340  kstep = 1
341  p = k
342 *
343 * Determine rows and columns to be interchanged and whether
344 * a 1-by-1 or 2-by-2 pivot block will be used
345 *
346  absakk = abs( dble( a( k, k ) ) )
347 *
348 * IMAX is the row-index of the largest off-diagonal element in
349 * column K, and COLMAX is its absolute value.
350 * Determine both COLMAX and IMAX.
351 *
352  IF( k.GT.1 ) THEN
353  imax = izamax( k-1, a( 1, k ), 1 )
354  colmax = cabs1( a( imax, k ) )
355  ELSE
356  colmax = zero
357  END IF
358 *
359  IF( ( max( absakk, colmax ).EQ.zero ) ) THEN
360 *
361 * Column K is zero or underflow: set INFO and continue
362 *
363  IF( info.EQ.0 )
364  $ info = k
365  kp = k
366  a( k, k ) = dble( a( k, k ) )
367 *
368 * Set E( K ) to zero
369 *
370  IF( k.GT.1 )
371  $ e( k ) = czero
372 *
373  ELSE
374 *
375 * ============================================================
376 *
377 * BEGIN pivot search
378 *
379 * Case(1)
380 * Equivalent to testing for ABSAKK.GE.ALPHA*COLMAX
381 * (used to handle NaN and Inf)
382 *
383  IF( .NOT.( absakk.LT.alpha*colmax ) ) THEN
384 *
385 * no interchange, use 1-by-1 pivot block
386 *
387  kp = k
388 *
389  ELSE
390 *
391  done = .false.
392 *
393 * Loop until pivot found
394 *
395  12 CONTINUE
396 *
397 * BEGIN pivot search loop body
398 *
399 *
400 * JMAX is the column-index of the largest off-diagonal
401 * element in row IMAX, and ROWMAX is its absolute value.
402 * Determine both ROWMAX and JMAX.
403 *
404  IF( imax.NE.k ) THEN
405  jmax = imax + izamax( k-imax, a( imax, imax+1 ),
406  $ lda )
407  rowmax = cabs1( a( imax, jmax ) )
408  ELSE
409  rowmax = zero
410  END IF
411 *
412  IF( imax.GT.1 ) THEN
413  itemp = izamax( imax-1, a( 1, imax ), 1 )
414  dtemp = cabs1( a( itemp, imax ) )
415  IF( dtemp.GT.rowmax ) THEN
416  rowmax = dtemp
417  jmax = itemp
418  END IF
419  END IF
420 *
421 * Case(2)
422 * Equivalent to testing for
423 * ABS( REAL( W( IMAX,KW-1 ) ) ).GE.ALPHA*ROWMAX
424 * (used to handle NaN and Inf)
425 *
426  IF( .NOT.( abs( dble( a( imax, imax ) ) )
427  $ .LT.alpha*rowmax ) ) THEN
428 *
429 * interchange rows and columns K and IMAX,
430 * use 1-by-1 pivot block
431 *
432  kp = imax
433  done = .true.
434 *
435 * Case(3)
436 * Equivalent to testing for ROWMAX.EQ.COLMAX,
437 * (used to handle NaN and Inf)
438 *
439  ELSE IF( ( p.EQ.jmax ) .OR. ( rowmax.LE.colmax ) )
440  $ THEN
441 *
442 * interchange rows and columns K-1 and IMAX,
443 * use 2-by-2 pivot block
444 *
445  kp = imax
446  kstep = 2
447  done = .true.
448 *
449 * Case(4)
450  ELSE
451 *
452 * Pivot not found: set params and repeat
453 *
454  p = imax
455  colmax = rowmax
456  imax = jmax
457  END IF
458 *
459 * END pivot search loop body
460 *
461  IF( .NOT.done ) GOTO 12
462 *
463  END IF
464 *
465 * END pivot search
466 *
467 * ============================================================
468 *
469 * KK is the column of A where pivoting step stopped
470 *
471  kk = k - kstep + 1
472 *
473 * For only a 2x2 pivot, interchange rows and columns K and P
474 * in the leading submatrix A(1:k,1:k)
475 *
476  IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) ) THEN
477 * (1) Swap columnar parts
478  IF( p.GT.1 )
479  $ CALL zswap( p-1, a( 1, k ), 1, a( 1, p ), 1 )
480 * (2) Swap and conjugate middle parts
481  DO 14 j = p + 1, k - 1
482  t = dconjg( a( j, k ) )
483  a( j, k ) = dconjg( a( p, j ) )
484  a( p, j ) = t
485  14 CONTINUE
486 * (3) Swap and conjugate corner elements at row-col interserction
487  a( p, k ) = dconjg( a( p, k ) )
488 * (4) Swap diagonal elements at row-col intersection
489  r1 = dble( a( k, k ) )
490  a( k, k ) = dble( a( p, p ) )
491  a( p, p ) = r1
492 *
493 * Convert upper triangle of A into U form by applying
494 * the interchanges in columns k+1:N.
495 *
496  IF( k.LT.n )
497  $ CALL zswap( n-k, a( k, k+1 ), lda, a( p, k+1 ), lda )
498 *
499  END IF
500 *
501 * For both 1x1 and 2x2 pivots, interchange rows and
502 * columns KK and KP in the leading submatrix A(1:k,1:k)
503 *
504  IF( kp.NE.kk ) THEN
505 * (1) Swap columnar parts
506  IF( kp.GT.1 )
507  $ CALL zswap( kp-1, a( 1, kk ), 1, a( 1, kp ), 1 )
508 * (2) Swap and conjugate middle parts
509  DO 15 j = kp + 1, kk - 1
510  t = dconjg( a( j, kk ) )
511  a( j, kk ) = dconjg( a( kp, j ) )
512  a( kp, j ) = t
513  15 CONTINUE
514 * (3) Swap and conjugate corner elements at row-col interserction
515  a( kp, kk ) = dconjg( a( kp, kk ) )
516 * (4) Swap diagonal elements at row-col intersection
517  r1 = dble( a( kk, kk ) )
518  a( kk, kk ) = dble( a( kp, kp ) )
519  a( kp, kp ) = r1
520 *
521  IF( kstep.EQ.2 ) THEN
522 * (*) Make sure that diagonal element of pivot is real
523  a( k, k ) = dble( a( k, k ) )
524 * (5) Swap row elements
525  t = a( k-1, k )
526  a( k-1, k ) = a( kp, k )
527  a( kp, k ) = t
528  END IF
529 *
530 * Convert upper triangle of A into U form by applying
531 * the interchanges in columns k+1:N.
532 *
533  IF( k.LT.n )
534  $ CALL zswap( n-k, a( kk, k+1 ), lda, a( kp, k+1 ),
535  $ lda )
536 *
537  ELSE
538 * (*) Make sure that diagonal element of pivot is real
539  a( k, k ) = dble( a( k, k ) )
540  IF( kstep.EQ.2 )
541  $ a( k-1, k-1 ) = dble( a( k-1, k-1 ) )
542  END IF
543 *
544 * Update the leading submatrix
545 *
546  IF( kstep.EQ.1 ) THEN
547 *
548 * 1-by-1 pivot block D(k): column k now holds
549 *
550 * W(k) = U(k)*D(k)
551 *
552 * where U(k) is the k-th column of U
553 *
554  IF( k.GT.1 ) THEN
555 *
556 * Perform a rank-1 update of A(1:k-1,1:k-1) and
557 * store U(k) in column k
558 *
559  IF( abs( dble( a( k, k ) ) ).GE.sfmin ) THEN
560 *
561 * Perform a rank-1 update of A(1:k-1,1:k-1) as
562 * A := A - U(k)*D(k)*U(k)**T
563 * = A - W(k)*1/D(k)*W(k)**T
564 *
565  d11 = one / dble( a( k, k ) )
566  CALL zher( uplo, k-1, -d11, a( 1, k ), 1, a, lda )
567 *
568 * Store U(k) in column k
569 *
570  CALL zdscal( k-1, d11, a( 1, k ), 1 )
571  ELSE
572 *
573 * Store L(k) in column K
574 *
575  d11 = dble( a( k, k ) )
576  DO 16 ii = 1, k - 1
577  a( ii, k ) = a( ii, k ) / d11
578  16 CONTINUE
579 *
580 * Perform a rank-1 update of A(k+1:n,k+1:n) as
581 * A := A - U(k)*D(k)*U(k)**T
582 * = A - W(k)*(1/D(k))*W(k)**T
583 * = A - (W(k)/D(k))*(D(k))*(W(k)/D(K))**T
584 *
585  CALL zher( uplo, k-1, -d11, a( 1, k ), 1, a, lda )
586  END IF
587 *
588 * Store the superdiagonal element of D in array E
589 *
590  e( k ) = czero
591 *
592  END IF
593 *
594  ELSE
595 *
596 * 2-by-2 pivot block D(k): columns k and k-1 now hold
597 *
598 * ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
599 *
600 * where U(k) and U(k-1) are the k-th and (k-1)-th columns
601 * of U
602 *
603 * Perform a rank-2 update of A(1:k-2,1:k-2) as
604 *
605 * A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )**T
606 * = A - ( ( A(k-1)A(k) )*inv(D(k)) ) * ( A(k-1)A(k) )**T
607 *
608 * and store L(k) and L(k+1) in columns k and k+1
609 *
610  IF( k.GT.2 ) THEN
611 * D = |A12|
612  d = dlapy2( dble( a( k-1, k ) ),
613  $ dimag( a( k-1, k ) ) )
614  d11 = a( k, k ) / d
615  d22 = a( k-1, k-1 ) / d
616  d12 = a( k-1, k ) / d
617  tt = one / ( d11*d22-one )
618 *
619  DO 30 j = k - 2, 1, -1
620 *
621 * Compute D21 * ( W(k)W(k+1) ) * inv(D(k)) for row J
622 *
623  wkm1 = tt*( d11*a( j, k-1 )-dconjg( d12 )*
624  $ a( j, k ) )
625  wk = tt*( d22*a( j, k )-d12*a( j, k-1 ) )
626 *
627 * Perform a rank-2 update of A(1:k-2,1:k-2)
628 *
629  DO 20 i = j, 1, -1
630  a( i, j ) = a( i, j ) -
631  $ ( a( i, k ) / d )*dconjg( wk ) -
632  $ ( a( i, k-1 ) / d )*dconjg( wkm1 )
633  20 CONTINUE
634 *
635 * Store U(k) and U(k-1) in cols k and k-1 for row J
636 *
637  a( j, k ) = wk / d
638  a( j, k-1 ) = wkm1 / d
639 * (*) Make sure that diagonal element of pivot is real
640  a( j, j ) = dcmplx( dble( a( j, j ) ), zero )
641 *
642  30 CONTINUE
643 *
644  END IF
645 *
646 * Copy superdiagonal elements of D(K) to E(K) and
647 * ZERO out superdiagonal entry of A
648 *
649  e( k ) = a( k-1, k )
650  e( k-1 ) = czero
651  a( k-1, k ) = czero
652 *
653  END IF
654 *
655 * End column K is nonsingular
656 *
657  END IF
658 *
659 * Store details of the interchanges in IPIV
660 *
661  IF( kstep.EQ.1 ) THEN
662  ipiv( k ) = kp
663  ELSE
664  ipiv( k ) = -p
665  ipiv( k-1 ) = -kp
666  END IF
667 *
668 * Decrease K and return to the start of the main loop
669 *
670  k = k - kstep
671  GO TO 10
672 *
673  34 CONTINUE
674 *
675  ELSE
676 *
677 * Factorize A as L*D*L**H using the lower triangle of A
678 *
679 * Initilize the unused last entry of the subdiagonal array E.
680 *
681  e( n ) = czero
682 *
683 * K is the main loop index, increasing from 1 to N in steps of
684 * 1 or 2
685 *
686  k = 1
687  40 CONTINUE
688 *
689 * If K > N, exit from loop
690 *
691  IF( k.GT.n )
692  $ GO TO 64
693  kstep = 1
694  p = k
695 *
696 * Determine rows and columns to be interchanged and whether
697 * a 1-by-1 or 2-by-2 pivot block will be used
698 *
699  absakk = abs( dble( a( k, k ) ) )
700 *
701 * IMAX is the row-index of the largest off-diagonal element in
702 * column K, and COLMAX is its absolute value.
703 * Determine both COLMAX and IMAX.
704 *
705  IF( k.LT.n ) THEN
706  imax = k + izamax( n-k, a( k+1, k ), 1 )
707  colmax = cabs1( a( imax, k ) )
708  ELSE
709  colmax = zero
710  END IF
711 *
712  IF( max( absakk, colmax ).EQ.zero ) THEN
713 *
714 * Column K is zero or underflow: set INFO and continue
715 *
716  IF( info.EQ.0 )
717  $ info = k
718  kp = k
719  a( k, k ) = dble( a( k, k ) )
720 *
721 * Set E( K ) to zero
722 *
723  IF( k.LT.n )
724  $ e( k ) = czero
725 *
726  ELSE
727 *
728 * ============================================================
729 *
730 * BEGIN pivot search
731 *
732 * Case(1)
733 * Equivalent to testing for ABSAKK.GE.ALPHA*COLMAX
734 * (used to handle NaN and Inf)
735 *
736  IF( .NOT.( absakk.LT.alpha*colmax ) ) THEN
737 *
738 * no interchange, use 1-by-1 pivot block
739 *
740  kp = k
741 *
742  ELSE
743 *
744  done = .false.
745 *
746 * Loop until pivot found
747 *
748  42 CONTINUE
749 *
750 * BEGIN pivot search loop body
751 *
752 *
753 * JMAX is the column-index of the largest off-diagonal
754 * element in row IMAX, and ROWMAX is its absolute value.
755 * Determine both ROWMAX and JMAX.
756 *
757  IF( imax.NE.k ) THEN
758  jmax = k - 1 + izamax( imax-k, a( imax, k ), lda )
759  rowmax = cabs1( a( imax, jmax ) )
760  ELSE
761  rowmax = zero
762  END IF
763 *
764  IF( imax.LT.n ) THEN
765  itemp = imax + izamax( n-imax, a( imax+1, imax ),
766  $ 1 )
767  dtemp = cabs1( a( itemp, imax ) )
768  IF( dtemp.GT.rowmax ) THEN
769  rowmax = dtemp
770  jmax = itemp
771  END IF
772  END IF
773 *
774 * Case(2)
775 * Equivalent to testing for
776 * ABS( REAL( W( IMAX,KW-1 ) ) ).GE.ALPHA*ROWMAX
777 * (used to handle NaN and Inf)
778 *
779  IF( .NOT.( abs( dble( a( imax, imax ) ) )
780  $ .LT.alpha*rowmax ) ) THEN
781 *
782 * interchange rows and columns K and IMAX,
783 * use 1-by-1 pivot block
784 *
785  kp = imax
786  done = .true.
787 *
788 * Case(3)
789 * Equivalent to testing for ROWMAX.EQ.COLMAX,
790 * (used to handle NaN and Inf)
791 *
792  ELSE IF( ( p.EQ.jmax ) .OR. ( rowmax.LE.colmax ) )
793  $ THEN
794 *
795 * interchange rows and columns K+1 and IMAX,
796 * use 2-by-2 pivot block
797 *
798  kp = imax
799  kstep = 2
800  done = .true.
801 *
802 * Case(4)
803  ELSE
804 *
805 * Pivot not found: set params and repeat
806 *
807  p = imax
808  colmax = rowmax
809  imax = jmax
810  END IF
811 *
812 *
813 * END pivot search loop body
814 *
815  IF( .NOT.done ) GOTO 42
816 *
817  END IF
818 *
819 * END pivot search
820 *
821 * ============================================================
822 *
823 * KK is the column of A where pivoting step stopped
824 *
825  kk = k + kstep - 1
826 *
827 * For only a 2x2 pivot, interchange rows and columns K and P
828 * in the trailing submatrix A(k:n,k:n)
829 *
830  IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) ) THEN
831 * (1) Swap columnar parts
832  IF( p.LT.n )
833  $ CALL zswap( n-p, a( p+1, k ), 1, a( p+1, p ), 1 )
834 * (2) Swap and conjugate middle parts
835  DO 44 j = k + 1, p - 1
836  t = dconjg( a( j, k ) )
837  a( j, k ) = dconjg( a( p, j ) )
838  a( p, j ) = t
839  44 CONTINUE
840 * (3) Swap and conjugate corner elements at row-col interserction
841  a( p, k ) = dconjg( a( p, k ) )
842 * (4) Swap diagonal elements at row-col intersection
843  r1 = dble( a( k, k ) )
844  a( k, k ) = dble( a( p, p ) )
845  a( p, p ) = r1
846 *
847 * Convert lower triangle of A into L form by applying
848 * the interchanges in columns 1:k-1.
849 *
850  IF ( k.GT.1 )
851  $ CALL zswap( k-1, a( k, 1 ), lda, a( p, 1 ), lda )
852 *
853  END IF
854 *
855 * For both 1x1 and 2x2 pivots, interchange rows and
856 * columns KK and KP in the trailing submatrix A(k:n,k:n)
857 *
858  IF( kp.NE.kk ) THEN
859 * (1) Swap columnar parts
860  IF( kp.LT.n )
861  $ CALL zswap( n-kp, a( kp+1, kk ), 1, a( kp+1, kp ), 1 )
862 * (2) Swap and conjugate middle parts
863  DO 45 j = kk + 1, kp - 1
864  t = dconjg( a( j, kk ) )
865  a( j, kk ) = dconjg( a( kp, j ) )
866  a( kp, j ) = t
867  45 CONTINUE
868 * (3) Swap and conjugate corner elements at row-col interserction
869  a( kp, kk ) = dconjg( a( kp, kk ) )
870 * (4) Swap diagonal elements at row-col intersection
871  r1 = dble( a( kk, kk ) )
872  a( kk, kk ) = dble( a( kp, kp ) )
873  a( kp, kp ) = r1
874 *
875  IF( kstep.EQ.2 ) THEN
876 * (*) Make sure that diagonal element of pivot is real
877  a( k, k ) = dble( a( k, k ) )
878 * (5) Swap row elements
879  t = a( k+1, k )
880  a( k+1, k ) = a( kp, k )
881  a( kp, k ) = t
882  END IF
883 *
884 * Convert lower triangle of A into L form by applying
885 * the interchanges in columns 1:k-1.
886 *
887  IF ( k.GT.1 )
888  $ CALL zswap( k-1, a( kk, 1 ), lda, a( kp, 1 ), lda )
889 *
890  ELSE
891 * (*) Make sure that diagonal element of pivot is real
892  a( k, k ) = dble( a( k, k ) )
893  IF( kstep.EQ.2 )
894  $ a( k+1, k+1 ) = dble( a( k+1, k+1 ) )
895  END IF
896 *
897 * Update the trailing submatrix
898 *
899  IF( kstep.EQ.1 ) THEN
900 *
901 * 1-by-1 pivot block D(k): column k of A now holds
902 *
903 * W(k) = L(k)*D(k),
904 *
905 * where L(k) is the k-th column of L
906 *
907  IF( k.LT.n ) THEN
908 *
909 * Perform a rank-1 update of A(k+1:n,k+1:n) and
910 * store L(k) in column k
911 *
912 * Handle division by a small number
913 *
914  IF( abs( dble( a( k, k ) ) ).GE.sfmin ) THEN
915 *
916 * Perform a rank-1 update of A(k+1:n,k+1:n) as
917 * A := A - L(k)*D(k)*L(k)**T
918 * = A - W(k)*(1/D(k))*W(k)**T
919 *
920  d11 = one / dble( a( k, k ) )
921  CALL zher( uplo, n-k, -d11, a( k+1, k ), 1,
922  $ a( k+1, k+1 ), lda )
923 *
924 * Store L(k) in column k
925 *
926  CALL zdscal( n-k, d11, a( k+1, k ), 1 )
927  ELSE
928 *
929 * Store L(k) in column k
930 *
931  d11 = dble( a( k, k ) )
932  DO 46 ii = k + 1, n
933  a( ii, k ) = a( ii, k ) / d11
934  46 CONTINUE
935 *
936 * Perform a rank-1 update of A(k+1:n,k+1:n) as
937 * A := A - L(k)*D(k)*L(k)**T
938 * = A - W(k)*(1/D(k))*W(k)**T
939 * = A - (W(k)/D(k))*(D(k))*(W(k)/D(K))**T
940 *
941  CALL zher( uplo, n-k, -d11, a( k+1, k ), 1,
942  $ a( k+1, k+1 ), lda )
943  END IF
944 *
945 * Store the subdiagonal element of D in array E
946 *
947  e( k ) = czero
948 *
949  END IF
950 *
951  ELSE
952 *
953 * 2-by-2 pivot block D(k): columns k and k+1 now hold
954 *
955 * ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k)
956 *
957 * where L(k) and L(k+1) are the k-th and (k+1)-th columns
958 * of L
959 *
960 *
961 * Perform a rank-2 update of A(k+2:n,k+2:n) as
962 *
963 * A := A - ( L(k) L(k+1) ) * D(k) * ( L(k) L(k+1) )**T
964 * = A - ( ( A(k)A(k+1) )*inv(D(k) ) * ( A(k)A(k+1) )**T
965 *
966 * and store L(k) and L(k+1) in columns k and k+1
967 *
968  IF( k.LT.n-1 ) THEN
969 * D = |A21|
970  d = dlapy2( dble( a( k+1, k ) ),
971  $ dimag( a( k+1, k ) ) )
972  d11 = dble( a( k+1, k+1 ) ) / d
973  d22 = dble( a( k, k ) ) / d
974  d21 = a( k+1, k ) / d
975  tt = one / ( d11*d22-one )
976 *
977  DO 60 j = k + 2, n
978 *
979 * Compute D21 * ( W(k)W(k+1) ) * inv(D(k)) for row J
980 *
981  wk = tt*( d11*a( j, k )-d21*a( j, k+1 ) )
982  wkp1 = tt*( d22*a( j, k+1 )-dconjg( d21 )*
983  $ a( j, k ) )
984 *
985 * Perform a rank-2 update of A(k+2:n,k+2:n)
986 *
987  DO 50 i = j, n
988  a( i, j ) = a( i, j ) -
989  $ ( a( i, k ) / d )*dconjg( wk ) -
990  $ ( a( i, k+1 ) / d )*dconjg( wkp1 )
991  50 CONTINUE
992 *
993 * Store L(k) and L(k+1) in cols k and k+1 for row J
994 *
995  a( j, k ) = wk / d
996  a( j, k+1 ) = wkp1 / d
997 * (*) Make sure that diagonal element of pivot is real
998  a( j, j ) = dcmplx( dble( a( j, j ) ), zero )
999 *
1000  60 CONTINUE
1001 *
1002  END IF
1003 *
1004 * Copy subdiagonal elements of D(K) to E(K) and
1005 * ZERO out subdiagonal entry of A
1006 *
1007  e( k ) = a( k+1, k )
1008  e( k+1 ) = czero
1009  a( k+1, k ) = czero
1010 *
1011  END IF
1012 *
1013 * End column K is nonsingular
1014 *
1015  END IF
1016 *
1017 * Store details of the interchanges in IPIV
1018 *
1019  IF( kstep.EQ.1 ) THEN
1020  ipiv( k ) = kp
1021  ELSE
1022  ipiv( k ) = -p
1023  ipiv( k+1 ) = -kp
1024  END IF
1025 *
1026 * Increase K and return to the start of the main loop
1027 *
1028  k = k + kstep
1029  GO TO 40
1030 *
1031  64 CONTINUE
1032 *
1033  END IF
1034 *
1035  RETURN
1036 *
1037 * End of ZHETF2_RK
1038 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
Definition: zswap.f:52
integer function izamax(N, ZX, INCX)
IZAMAX
Definition: izamax.f:53
subroutine zher(UPLO, N, ALPHA, X, INCX, A, LDA)
ZHER
Definition: zher.f:137
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine zdscal(N, DA, ZX, INCX)
ZDSCAL
Definition: zdscal.f:54
double precision function dlapy2(X, Y)
DLAPY2 returns sqrt(x2+y2).
Definition: dlapy2.f:65

Here is the call graph for this function:

Here is the caller graph for this function: