LAPACK  3.7.0
LAPACK: Linear Algebra PACKage
subroutine zlahef_rk ( character  UPLO,
integer  N,
integer  NB,
integer  KB,
complex*16, dimension( lda, * )  A,
integer  LDA,
complex*16, dimension( * )  E,
integer, dimension( * )  IPIV,
complex*16, dimension( ldw, * )  W,
integer  LDW,
integer  INFO 
)

ZLAHEF_RK computes a partial factorization of a complex Hermitian indefinite matrix using bounded Bunch-Kaufman (rook) diagonal pivoting method.

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

Purpose:
 ZLAHEF_RK computes a partial factorization of a complex Hermitian
 matrix A using the bounded Bunch-Kaufman (rook) diagonal
 pivoting method. The partial factorization has the form:

 A  =  ( I  U12 ) ( A11  0  ) (  I       0    )  if UPLO = 'U', or:
       ( 0  U22 ) (  0   D  ) ( U12**H U22**H )

 A  =  ( L11  0 ) (  D   0  ) ( L11**H L21**H )  if UPLO = 'L',
       ( L21  I ) (  0  A22 ) (  0       I    )

 where the order of D is at most NB. The actual order is returned in
 the argument KB, and is either NB or NB-1, or N if N <= NB.

 ZLAHEF_RK is an auxiliary routine called by ZHETRF_RK. It uses
 blocked code (calling Level 3 BLAS) to update the submatrix
 A11 (if UPLO = 'U') or A22 (if UPLO = 'L').
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]NB
          NB is INTEGER
          The maximum number of columns of the matrix A that should be
          factored.  NB should be at least 2 to allow for 2-by-2 pivot
          blocks.
[out]KB
          KB is INTEGER
          The number of columns of A that were actually factored.
          KB is either NB-1 or NB, or N if N <= NB.
[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.

          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 submatrix A(1:N,N-KB+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,N-KB+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 submatrix A(1:N,N-KB+1:N).
                  If -IPIV(k-1) = k-1, no interchange occurred.

            c) In both cases a) and b) is 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 submatrix A(1:N,1:KB).
               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 submatrix A(1:N,1:KB).
                  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 submatrix A(1:N,1:KB).
                  If -IPIV(k+1) = k+1, no interchange occurred.

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

            d) NOTE: Any entry IPIV(k) is always NONZERO on output.
[out]W
          W is COMPLEX*16 array, dimension (LDW,NB)
[in]LDW
          LDW is INTEGER
          The leading dimension of the array W.  LDW >= max(1,N).
[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
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

Definition at line 264 of file zlahef_rk.f.

264 *
265 * -- LAPACK computational routine (version 3.7.0) --
266 * -- LAPACK is a software package provided by Univ. of Tennessee, --
267 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
268 * December 2016
269 *
270 * .. Scalar Arguments ..
271  CHARACTER uplo
272  INTEGER info, kb, lda, ldw, n, nb
273 * ..
274 * .. Array Arguments ..
275  INTEGER ipiv( * )
276  COMPLEX*16 a( lda, * ), w( ldw, * ), e( * )
277 * ..
278 *
279 * =====================================================================
280 *
281 * .. Parameters ..
282  DOUBLE PRECISION zero, one
283  parameter ( zero = 0.0d+0, one = 1.0d+0 )
284  COMPLEX*16 cone
285  parameter ( cone = ( 1.0d+0, 0.0d+0 ) )
286  DOUBLE PRECISION eight, sevten
287  parameter ( eight = 8.0d+0, sevten = 17.0d+0 )
288  COMPLEX*16 czero
289  parameter ( czero = ( 0.0d+0, 0.0d+0 ) )
290 * ..
291 * .. Local Scalars ..
292  LOGICAL done
293  INTEGER imax, itemp, ii, j, jb, jj, jmax, k, kk, kkw,
294  $ kp, kstep, kw, p
295  DOUBLE PRECISION absakk, alpha, colmax, dtemp, r1, rowmax, t,
296  $ sfmin
297  COMPLEX*16 d11, d21, d22, z
298 * ..
299 * .. External Functions ..
300  LOGICAL lsame
301  INTEGER izamax
302  DOUBLE PRECISION dlamch
303  EXTERNAL lsame, izamax, dlamch
304 * ..
305 * .. External Subroutines ..
306  EXTERNAL zcopy, zdscal, zgemm, zgemv, zlacgv, zswap
307 * ..
308 * .. Intrinsic Functions ..
309  INTRINSIC abs, dble, dconjg, dimag, max, min, sqrt
310 * ..
311 * .. Statement Functions ..
312  DOUBLE PRECISION cabs1
313 * ..
314 * .. Statement Function definitions ..
315  cabs1( z ) = abs( dble( z ) ) + abs( dimag( z ) )
316 * ..
317 * .. Executable Statements ..
318 *
319  info = 0
320 *
321 * Initialize ALPHA for use in choosing pivot block size.
322 *
323  alpha = ( one+sqrt( sevten ) ) / eight
324 *
325 * Compute machine safe minimum
326 *
327  sfmin = dlamch( 'S' )
328 *
329  IF( lsame( uplo, 'U' ) ) THEN
330 *
331 * Factorize the trailing columns of A using the upper triangle
332 * of A and working backwards, and compute the matrix W = U12*D
333 * for use in updating A11 (note that conjg(W) is actually stored)
334 * Initilize the first entry of array E, where superdiagonal
335 * elements of D are stored
336 *
337  e( 1 ) = czero
338 *
339 * K is the main loop index, decreasing from N in steps of 1 or 2
340 *
341  k = n
342  10 CONTINUE
343 *
344 * KW is the column of W which corresponds to column K of A
345 *
346  kw = nb + k - n
347 *
348 * Exit from loop
349 *
350  IF( ( k.LE.n-nb+1 .AND. nb.LT.n ) .OR. k.LT.1 )
351  $ GO TO 30
352 *
353  kstep = 1
354  p = k
355 *
356 * Copy column K of A to column KW of W and update it
357 *
358  IF( k.GT.1 )
359  $ CALL zcopy( k-1, a( 1, k ), 1, w( 1, kw ), 1 )
360  w( k, kw ) = dble( a( k, k ) )
361  IF( k.LT.n ) THEN
362  CALL zgemv( 'No transpose', k, n-k, -cone, a( 1, k+1 ), lda,
363  $ w( k, kw+1 ), ldw, cone, w( 1, kw ), 1 )
364  w( k, kw ) = dble( w( k, kw ) )
365  END IF
366 *
367 * Determine rows and columns to be interchanged and whether
368 * a 1-by-1 or 2-by-2 pivot block will be used
369 *
370  absakk = abs( dble( w( k, kw ) ) )
371 *
372 * IMAX is the row-index of the largest off-diagonal element in
373 * column K, and COLMAX is its absolute value.
374 * Determine both COLMAX and IMAX.
375 *
376  IF( k.GT.1 ) THEN
377  imax = izamax( k-1, w( 1, kw ), 1 )
378  colmax = cabs1( w( imax, kw ) )
379  ELSE
380  colmax = zero
381  END IF
382 *
383  IF( max( absakk, colmax ).EQ.zero ) THEN
384 *
385 * Column K is zero or underflow: set INFO and continue
386 *
387  IF( info.EQ.0 )
388  $ info = k
389  kp = k
390  a( k, k ) = dble( w( k, kw ) )
391  IF( k.GT.1 )
392  $ CALL zcopy( k-1, w( 1, kw ), 1, a( 1, k ), 1 )
393 *
394 * Set E( K ) to zero
395 *
396  IF( k.GT.1 )
397  $ e( k ) = czero
398 *
399  ELSE
400 *
401 * ============================================================
402 *
403 * BEGIN pivot search
404 *
405 * Case(1)
406 * Equivalent to testing for ABSAKK.GE.ALPHA*COLMAX
407 * (used to handle NaN and Inf)
408  IF( .NOT.( absakk.LT.alpha*colmax ) ) THEN
409 *
410 * no interchange, use 1-by-1 pivot block
411 *
412  kp = k
413 *
414  ELSE
415 *
416 * Lop until pivot found
417 *
418  done = .false.
419 *
420  12 CONTINUE
421 *
422 * BEGIN pivot search loop body
423 *
424 *
425 * Copy column IMAX to column KW-1 of W and update it
426 *
427  IF( imax.GT.1 )
428  $ CALL zcopy( imax-1, a( 1, imax ), 1, w( 1, kw-1 ),
429  $ 1 )
430  w( imax, kw-1 ) = dble( a( imax, imax ) )
431 *
432  CALL zcopy( k-imax, a( imax, imax+1 ), lda,
433  $ w( imax+1, kw-1 ), 1 )
434  CALL zlacgv( k-imax, w( imax+1, kw-1 ), 1 )
435 *
436  IF( k.LT.n ) THEN
437  CALL zgemv( 'No transpose', k, n-k, -cone,
438  $ a( 1, k+1 ), lda, w( imax, kw+1 ), ldw,
439  $ cone, w( 1, kw-1 ), 1 )
440  w( imax, kw-1 ) = dble( w( imax, kw-1 ) )
441  END IF
442 *
443 * JMAX is the column-index of the largest off-diagonal
444 * element in row IMAX, and ROWMAX is its absolute value.
445 * Determine both ROWMAX and JMAX.
446 *
447  IF( imax.NE.k ) THEN
448  jmax = imax + izamax( k-imax, w( imax+1, kw-1 ),
449  $ 1 )
450  rowmax = cabs1( w( jmax, kw-1 ) )
451  ELSE
452  rowmax = zero
453  END IF
454 *
455  IF( imax.GT.1 ) THEN
456  itemp = izamax( imax-1, w( 1, kw-1 ), 1 )
457  dtemp = cabs1( w( itemp, kw-1 ) )
458  IF( dtemp.GT.rowmax ) THEN
459  rowmax = dtemp
460  jmax = itemp
461  END IF
462  END IF
463 *
464 * Case(2)
465 * Equivalent to testing for
466 * ABS( REAL( W( IMAX,KW-1 ) ) ).GE.ALPHA*ROWMAX
467 * (used to handle NaN and Inf)
468 *
469  IF( .NOT.( abs( dble( w( imax,kw-1 ) ) )
470  $ .LT.alpha*rowmax ) ) THEN
471 *
472 * interchange rows and columns K and IMAX,
473 * use 1-by-1 pivot block
474 *
475  kp = imax
476 *
477 * copy column KW-1 of W to column KW of W
478 *
479  CALL zcopy( k, w( 1, kw-1 ), 1, w( 1, kw ), 1 )
480 *
481  done = .true.
482 *
483 * Case(3)
484 * Equivalent to testing for ROWMAX.EQ.COLMAX,
485 * (used to handle NaN and Inf)
486 *
487  ELSE IF( ( p.EQ.jmax ) .OR. ( rowmax.LE.colmax ) )
488  $ THEN
489 *
490 * interchange rows and columns K-1 and IMAX,
491 * use 2-by-2 pivot block
492 *
493  kp = imax
494  kstep = 2
495  done = .true.
496 *
497 * Case(4)
498  ELSE
499 *
500 * Pivot not found: set params and repeat
501 *
502  p = imax
503  colmax = rowmax
504  imax = jmax
505 *
506 * Copy updated JMAXth (next IMAXth) column to Kth of W
507 *
508  CALL zcopy( k, w( 1, kw-1 ), 1, w( 1, kw ), 1 )
509 *
510  END IF
511 *
512 *
513 * END pivot search loop body
514 *
515  IF( .NOT.done ) GOTO 12
516 *
517  END IF
518 *
519 * END pivot search
520 *
521 * ============================================================
522 *
523 * KK is the column of A where pivoting step stopped
524 *
525  kk = k - kstep + 1
526 *
527 * KKW is the column of W which corresponds to column KK of A
528 *
529  kkw = nb + kk - n
530 *
531 * Interchange rows and columns P and K.
532 * Updated column P is already stored in column KW of W.
533 *
534  IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) ) THEN
535 *
536 * Copy non-updated column K to column P of submatrix A
537 * at step K. No need to copy element into columns
538 * K and K-1 of A for 2-by-2 pivot, since these columns
539 * will be later overwritten.
540 *
541  a( p, p ) = dble( a( k, k ) )
542  CALL zcopy( k-1-p, a( p+1, k ), 1, a( p, p+1 ),
543  $ lda )
544  CALL zlacgv( k-1-p, a( p, p+1 ), lda )
545  IF( p.GT.1 )
546  $ CALL zcopy( p-1, a( 1, k ), 1, a( 1, p ), 1 )
547 *
548 * Interchange rows K and P in the last K+1 to N columns of A
549 * (columns K and K-1 of A for 2-by-2 pivot will be
550 * later overwritten). Interchange rows K and P
551 * in last KKW to NB columns of W.
552 *
553  IF( k.LT.n )
554  $ CALL zswap( n-k, a( k, k+1 ), lda, a( p, k+1 ),
555  $ lda )
556  CALL zswap( n-kk+1, w( k, kkw ), ldw, w( p, kkw ),
557  $ ldw )
558  END IF
559 *
560 * Interchange rows and columns KP and KK.
561 * Updated column KP is already stored in column KKW of W.
562 *
563  IF( kp.NE.kk ) THEN
564 *
565 * Copy non-updated column KK to column KP of submatrix A
566 * at step K. No need to copy element into column K
567 * (or K and K-1 for 2-by-2 pivot) of A, since these columns
568 * will be later overwritten.
569 *
570  a( kp, kp ) = dble( a( kk, kk ) )
571  CALL zcopy( kk-1-kp, a( kp+1, kk ), 1, a( kp, kp+1 ),
572  $ lda )
573  CALL zlacgv( kk-1-kp, a( kp, kp+1 ), lda )
574  IF( kp.GT.1 )
575  $ CALL zcopy( kp-1, a( 1, kk ), 1, a( 1, kp ), 1 )
576 *
577 * Interchange rows KK and KP in last K+1 to N columns of A
578 * (columns K (or K and K-1 for 2-by-2 pivot) of A will be
579 * later overwritten). Interchange rows KK and KP
580 * in last KKW to NB columns of W.
581 *
582  IF( k.LT.n )
583  $ CALL zswap( n-k, a( kk, k+1 ), lda, a( kp, k+1 ),
584  $ lda )
585  CALL zswap( n-kk+1, w( kk, kkw ), ldw, w( kp, kkw ),
586  $ ldw )
587  END IF
588 *
589  IF( kstep.EQ.1 ) THEN
590 *
591 * 1-by-1 pivot block D(k): column kw of W now holds
592 *
593 * W(kw) = U(k)*D(k),
594 *
595 * where U(k) is the k-th column of U
596 *
597 * (1) Store subdiag. elements of column U(k)
598 * and 1-by-1 block D(k) in column k of A.
599 * (NOTE: Diagonal element U(k,k) is a UNIT element
600 * and not stored)
601 * A(k,k) := D(k,k) = W(k,kw)
602 * A(1:k-1,k) := U(1:k-1,k) = W(1:k-1,kw)/D(k,k)
603 *
604 * (NOTE: No need to use for Hermitian matrix
605 * A( K, K ) = REAL( W( K, K) ) to separately copy diagonal
606 * element D(k,k) from W (potentially saves only one load))
607  CALL zcopy( k, w( 1, kw ), 1, a( 1, k ), 1 )
608  IF( k.GT.1 ) THEN
609 *
610 * (NOTE: No need to check if A(k,k) is NOT ZERO,
611 * since that was ensured earlier in pivot search:
612 * case A(k,k) = 0 falls into 2x2 pivot case(3))
613 *
614 * Handle division by a small number
615 *
616  t = dble( a( k, k ) )
617  IF( abs( t ).GE.sfmin ) THEN
618  r1 = one / t
619  CALL zdscal( k-1, r1, a( 1, k ), 1 )
620  ELSE
621  DO 14 ii = 1, k-1
622  a( ii, k ) = a( ii, k ) / t
623  14 CONTINUE
624  END IF
625 *
626 * (2) Conjugate column W(kw)
627 *
628  CALL zlacgv( k-1, w( 1, kw ), 1 )
629 *
630 * Store the superdiagonal element of D in array E
631 *
632  e( k ) = czero
633 *
634  END IF
635 *
636  ELSE
637 *
638 * 2-by-2 pivot block D(k): columns kw and kw-1 of W now hold
639 *
640 * ( W(kw-1) W(kw) ) = ( U(k-1) U(k) )*D(k)
641 *
642 * where U(k) and U(k-1) are the k-th and (k-1)-th columns
643 * of U
644 *
645 * (1) Store U(1:k-2,k-1) and U(1:k-2,k) and 2-by-2
646 * block D(k-1:k,k-1:k) in columns k-1 and k of A.
647 * (NOTE: 2-by-2 diagonal block U(k-1:k,k-1:k) is a UNIT
648 * block and not stored)
649 * A(k-1:k,k-1:k) := D(k-1:k,k-1:k) = W(k-1:k,kw-1:kw)
650 * A(1:k-2,k-1:k) := U(1:k-2,k:k-1:k) =
651 * = W(1:k-2,kw-1:kw) * ( D(k-1:k,k-1:k)**(-1) )
652 *
653  IF( k.GT.2 ) THEN
654 *
655 * Factor out the columns of the inverse of 2-by-2 pivot
656 * block D, so that each column contains 1, to reduce the
657 * number of FLOPS when we multiply panel
658 * ( W(kw-1) W(kw) ) by this inverse, i.e. by D**(-1).
659 *
660 * D**(-1) = ( d11 cj(d21) )**(-1) =
661 * ( d21 d22 )
662 *
663 * = 1/(d11*d22-|d21|**2) * ( ( d22) (-cj(d21) ) ) =
664 * ( (-d21) ( d11 ) )
665 *
666 * = 1/(|d21|**2) * 1/((d11/cj(d21))*(d22/d21)-1) *
667 *
668 * * ( d21*( d22/d21 ) conj(d21)*( - 1 ) ) =
669 * ( ( -1 ) ( d11/conj(d21) ) )
670 *
671 * = 1/(|d21|**2) * 1/(D22*D11-1) *
672 *
673 * * ( d21*( D11 ) conj(d21)*( -1 ) ) =
674 * ( ( -1 ) ( D22 ) )
675 *
676 * = (1/|d21|**2) * T * ( d21*( D11 ) conj(d21)*( -1 ) ) =
677 * ( ( -1 ) ( D22 ) )
678 *
679 * = ( (T/conj(d21))*( D11 ) (T/d21)*( -1 ) ) =
680 * ( ( -1 ) ( D22 ) )
681 *
682 * Handle division by a small number. (NOTE: order of
683 * operations is important)
684 *
685 * = ( T*(( D11 )/conj(D21)) T*(( -1 )/D21 ) )
686 * ( (( -1 ) ) (( D22 ) ) ),
687 *
688 * where D11 = d22/d21,
689 * D22 = d11/conj(d21),
690 * D21 = d21,
691 * T = 1/(D22*D11-1).
692 *
693 * (NOTE: No need to check for division by ZERO,
694 * since that was ensured earlier in pivot search:
695 * (a) d21 != 0 in 2x2 pivot case(4),
696 * since |d21| should be larger than |d11| and |d22|;
697 * (b) (D22*D11 - 1) != 0, since from (a),
698 * both |D11| < 1, |D22| < 1, hence |D22*D11| << 1.)
699 *
700  d21 = w( k-1, kw )
701  d11 = w( k, kw ) / dconjg( d21 )
702  d22 = w( k-1, kw-1 ) / d21
703  t = one / ( dble( d11*d22 )-one )
704 *
705 * Update elements in columns A(k-1) and A(k) as
706 * dot products of rows of ( W(kw-1) W(kw) ) and columns
707 * of D**(-1)
708 *
709  DO 20 j = 1, k - 2
710  a( j, k-1 ) = t*( ( d11*w( j, kw-1 )-w( j, kw ) ) /
711  $ d21 )
712  a( j, k ) = t*( ( d22*w( j, kw )-w( j, kw-1 ) ) /
713  $ dconjg( d21 ) )
714  20 CONTINUE
715  END IF
716 *
717 * Copy diagonal elements of D(K) to A,
718 * copy superdiagonal element of D(K) to E(K) and
719 * ZERO out superdiagonal entry of A
720 *
721  a( k-1, k-1 ) = w( k-1, kw-1 )
722  a( k-1, k ) = czero
723  a( k, k ) = w( k, kw )
724  e( k ) = w( k-1, kw )
725  e( k-1 ) = czero
726 *
727 * (2) Conjugate columns W(kw) and W(kw-1)
728 *
729  CALL zlacgv( k-1, w( 1, kw ), 1 )
730  CALL zlacgv( k-2, w( 1, kw-1 ), 1 )
731 *
732  END IF
733 *
734 * End column K is nonsingular
735 *
736  END IF
737 *
738 * Store details of the interchanges in IPIV
739 *
740  IF( kstep.EQ.1 ) THEN
741  ipiv( k ) = kp
742  ELSE
743  ipiv( k ) = -p
744  ipiv( k-1 ) = -kp
745  END IF
746 *
747 * Decrease K and return to the start of the main loop
748 *
749  k = k - kstep
750  GO TO 10
751 *
752  30 CONTINUE
753 *
754 * Update the upper triangle of A11 (= A(1:k,1:k)) as
755 *
756 * A11 := A11 - U12*D*U12**H = A11 - U12*W**H
757 *
758 * computing blocks of NB columns at a time (note that conjg(W) is
759 * actually stored)
760 *
761  DO 50 j = ( ( k-1 ) / nb )*nb + 1, 1, -nb
762  jb = min( nb, k-j+1 )
763 *
764 * Update the upper triangle of the diagonal block
765 *
766  DO 40 jj = j, j + jb - 1
767  a( jj, jj ) = dble( a( jj, jj ) )
768  CALL zgemv( 'No transpose', jj-j+1, n-k, -cone,
769  $ a( j, k+1 ), lda, w( jj, kw+1 ), ldw, cone,
770  $ a( j, jj ), 1 )
771  a( jj, jj ) = dble( a( jj, jj ) )
772  40 CONTINUE
773 *
774 * Update the rectangular superdiagonal block
775 *
776  IF( j.GE.2 )
777  $ CALL zgemm( 'No transpose', 'Transpose', j-1, jb, n-k,
778  $ -cone, a( 1, k+1 ), lda, w( j, kw+1 ), ldw,
779  $ cone, a( 1, j ), lda )
780  50 CONTINUE
781 *
782 * Set KB to the number of columns factorized
783 *
784  kb = n - k
785 *
786  ELSE
787 *
788 * Factorize the leading columns of A using the lower triangle
789 * of A and working forwards, and compute the matrix W = L21*D
790 * for use in updating A22 (note that conjg(W) is actually stored)
791 *
792 * Initilize the unused last entry of the subdiagonal array E.
793 *
794  e( n ) = czero
795 *
796 * K is the main loop index, increasing from 1 in steps of 1 or 2
797 *
798  k = 1
799  70 CONTINUE
800 *
801 * Exit from loop
802 *
803  IF( ( k.GE.nb .AND. nb.LT.n ) .OR. k.GT.n )
804  $ GO TO 90
805 *
806  kstep = 1
807  p = k
808 *
809 * Copy column K of A to column K of W and update column K of W
810 *
811  w( k, k ) = dble( a( k, k ) )
812  IF( k.LT.n )
813  $ CALL zcopy( n-k, a( k+1, k ), 1, w( k+1, k ), 1 )
814  IF( k.GT.1 ) THEN
815  CALL zgemv( 'No transpose', n-k+1, k-1, -cone, a( k, 1 ),
816  $ lda, w( k, 1 ), ldw, cone, w( k, k ), 1 )
817  w( k, k ) = dble( w( k, k ) )
818  END IF
819 *
820 * Determine rows and columns to be interchanged and whether
821 * a 1-by-1 or 2-by-2 pivot block will be used
822 *
823  absakk = abs( dble( w( k, k ) ) )
824 *
825 * IMAX is the row-index of the largest off-diagonal element in
826 * column K, and COLMAX is its absolute value.
827 * Determine both COLMAX and IMAX.
828 *
829  IF( k.LT.n ) THEN
830  imax = k + izamax( n-k, w( k+1, k ), 1 )
831  colmax = cabs1( w( imax, k ) )
832  ELSE
833  colmax = zero
834  END IF
835 *
836  IF( max( absakk, colmax ).EQ.zero ) THEN
837 *
838 * Column K is zero or underflow: set INFO and continue
839 *
840  IF( info.EQ.0 )
841  $ info = k
842  kp = k
843  a( k, k ) = dble( w( k, k ) )
844  IF( k.LT.n )
845  $ CALL zcopy( n-k, w( k+1, k ), 1, a( k+1, k ), 1 )
846 *
847 * Set E( K ) to zero
848 *
849  IF( k.LT.n )
850  $ e( k ) = czero
851 *
852  ELSE
853 *
854 * ============================================================
855 *
856 * BEGIN pivot search
857 *
858 * Case(1)
859 * Equivalent to testing for ABSAKK.GE.ALPHA*COLMAX
860 * (used to handle NaN and Inf)
861 *
862  IF( .NOT.( absakk.LT.alpha*colmax ) ) THEN
863 *
864 * no interchange, use 1-by-1 pivot block
865 *
866  kp = k
867 *
868  ELSE
869 *
870  done = .false.
871 *
872 * Loop until pivot found
873 *
874  72 CONTINUE
875 *
876 * BEGIN pivot search loop body
877 *
878 *
879 * Copy column IMAX to column k+1 of W and update it
880 *
881  CALL zcopy( imax-k, a( imax, k ), lda, w( k, k+1 ), 1)
882  CALL zlacgv( imax-k, w( k, k+1 ), 1 )
883  w( imax, k+1 ) = dble( a( imax, imax ) )
884 *
885  IF( imax.LT.n )
886  $ CALL zcopy( n-imax, a( imax+1, imax ), 1,
887  $ w( imax+1, k+1 ), 1 )
888 *
889  IF( k.GT.1 ) THEN
890  CALL zgemv( 'No transpose', n-k+1, k-1, -cone,
891  $ a( k, 1 ), lda, w( imax, 1 ), ldw,
892  $ cone, w( k, k+1 ), 1 )
893  w( imax, k+1 ) = dble( w( imax, k+1 ) )
894  END IF
895 *
896 * JMAX is the column-index of the largest off-diagonal
897 * element in row IMAX, and ROWMAX is its absolute value.
898 * Determine both ROWMAX and JMAX.
899 *
900  IF( imax.NE.k ) THEN
901  jmax = k - 1 + izamax( imax-k, w( k, k+1 ), 1 )
902  rowmax = cabs1( w( jmax, k+1 ) )
903  ELSE
904  rowmax = zero
905  END IF
906 *
907  IF( imax.LT.n ) THEN
908  itemp = imax + izamax( n-imax, w( imax+1, k+1 ), 1)
909  dtemp = cabs1( w( itemp, k+1 ) )
910  IF( dtemp.GT.rowmax ) THEN
911  rowmax = dtemp
912  jmax = itemp
913  END IF
914  END IF
915 *
916 * Case(2)
917 * Equivalent to testing for
918 * ABS( REAL( W( IMAX,K+1 ) ) ).GE.ALPHA*ROWMAX
919 * (used to handle NaN and Inf)
920 *
921  IF( .NOT.( abs( dble( w( imax,k+1 ) ) )
922  $ .LT.alpha*rowmax ) ) THEN
923 *
924 * interchange rows and columns K and IMAX,
925 * use 1-by-1 pivot block
926 *
927  kp = imax
928 *
929 * copy column K+1 of W to column K of W
930 *
931  CALL zcopy( n-k+1, w( k, k+1 ), 1, w( k, k ), 1 )
932 *
933  done = .true.
934 *
935 * Case(3)
936 * Equivalent to testing for ROWMAX.EQ.COLMAX,
937 * (used to handle NaN and Inf)
938 *
939  ELSE IF( ( p.EQ.jmax ) .OR. ( rowmax.LE.colmax ) )
940  $ THEN
941 *
942 * interchange rows and columns K+1 and IMAX,
943 * use 2-by-2 pivot block
944 *
945  kp = imax
946  kstep = 2
947  done = .true.
948 *
949 * Case(4)
950  ELSE
951 *
952 * Pivot not found: set params and repeat
953 *
954  p = imax
955  colmax = rowmax
956  imax = jmax
957 *
958 * Copy updated JMAXth (next IMAXth) column to Kth of W
959 *
960  CALL zcopy( n-k+1, w( k, k+1 ), 1, w( k, k ), 1 )
961 *
962  END IF
963 *
964 *
965 * End pivot search loop body
966 *
967  IF( .NOT.done ) GOTO 72
968 *
969  END IF
970 *
971 * END pivot search
972 *
973 * ============================================================
974 *
975 * KK is the column of A where pivoting step stopped
976 *
977  kk = k + kstep - 1
978 *
979 * Interchange rows and columns P and K (only for 2-by-2 pivot).
980 * Updated column P is already stored in column K of W.
981 *
982  IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) ) THEN
983 *
984 * Copy non-updated column KK-1 to column P of submatrix A
985 * at step K. No need to copy element into columns
986 * K and K+1 of A for 2-by-2 pivot, since these columns
987 * will be later overwritten.
988 *
989  a( p, p ) = dble( a( k, k ) )
990  CALL zcopy( p-k-1, a( k+1, k ), 1, a( p, k+1 ), lda )
991  CALL zlacgv( p-k-1, a( p, k+1 ), lda )
992  IF( p.LT.n )
993  $ CALL zcopy( n-p, a( p+1, k ), 1, a( p+1, p ), 1 )
994 *
995 * Interchange rows K and P in first K-1 columns of A
996 * (columns K and K+1 of A for 2-by-2 pivot will be
997 * later overwritten). Interchange rows K and P
998 * in first KK columns of W.
999 *
1000  IF( k.GT.1 )
1001  $ CALL zswap( k-1, a( k, 1 ), lda, a( p, 1 ), lda )
1002  CALL zswap( kk, w( k, 1 ), ldw, w( p, 1 ), ldw )
1003  END IF
1004 *
1005 * Interchange rows and columns KP and KK.
1006 * Updated column KP is already stored in column KK of W.
1007 *
1008  IF( kp.NE.kk ) THEN
1009 *
1010 * Copy non-updated column KK to column KP of submatrix A
1011 * at step K. No need to copy element into column K
1012 * (or K and K+1 for 2-by-2 pivot) of A, since these columns
1013 * will be later overwritten.
1014 *
1015  a( kp, kp ) = dble( a( kk, kk ) )
1016  CALL zcopy( kp-kk-1, a( kk+1, kk ), 1, a( kp, kk+1 ),
1017  $ lda )
1018  CALL zlacgv( kp-kk-1, a( kp, kk+1 ), lda )
1019  IF( kp.LT.n )
1020  $ CALL zcopy( n-kp, a( kp+1, kk ), 1, a( kp+1, kp ), 1 )
1021 *
1022 * Interchange rows KK and KP in first K-1 columns of A
1023 * (column K (or K and K+1 for 2-by-2 pivot) of A will be
1024 * later overwritten). Interchange rows KK and KP
1025 * in first KK columns of W.
1026 *
1027  IF( k.GT.1 )
1028  $ CALL zswap( k-1, a( kk, 1 ), lda, a( kp, 1 ), lda )
1029  CALL zswap( kk, w( kk, 1 ), ldw, w( kp, 1 ), ldw )
1030  END IF
1031 *
1032  IF( kstep.EQ.1 ) THEN
1033 *
1034 * 1-by-1 pivot block D(k): column k of W now holds
1035 *
1036 * W(k) = L(k)*D(k),
1037 *
1038 * where L(k) is the k-th column of L
1039 *
1040 * (1) Store subdiag. elements of column L(k)
1041 * and 1-by-1 block D(k) in column k of A.
1042 * (NOTE: Diagonal element L(k,k) is a UNIT element
1043 * and not stored)
1044 * A(k,k) := D(k,k) = W(k,k)
1045 * A(k+1:N,k) := L(k+1:N,k) = W(k+1:N,k)/D(k,k)
1046 *
1047 * (NOTE: No need to use for Hermitian matrix
1048 * A( K, K ) = REAL( W( K, K) ) to separately copy diagonal
1049 * element D(k,k) from W (potentially saves only one load))
1050  CALL zcopy( n-k+1, w( k, k ), 1, a( k, k ), 1 )
1051  IF( k.LT.n ) THEN
1052 *
1053 * (NOTE: No need to check if A(k,k) is NOT ZERO,
1054 * since that was ensured earlier in pivot search:
1055 * case A(k,k) = 0 falls into 2x2 pivot case(3))
1056 *
1057 * Handle division by a small number
1058 *
1059  t = dble( a( k, k ) )
1060  IF( abs( t ).GE.sfmin ) THEN
1061  r1 = one / t
1062  CALL zdscal( n-k, r1, a( k+1, k ), 1 )
1063  ELSE
1064  DO 74 ii = k + 1, n
1065  a( ii, k ) = a( ii, k ) / t
1066  74 CONTINUE
1067  END IF
1068 *
1069 * (2) Conjugate column W(k)
1070 *
1071  CALL zlacgv( n-k, w( k+1, k ), 1 )
1072 *
1073 * Store the subdiagonal element of D in array E
1074 *
1075  e( k ) = czero
1076 *
1077  END IF
1078 *
1079  ELSE
1080 *
1081 * 2-by-2 pivot block D(k): columns k and k+1 of W now hold
1082 *
1083 * ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k)
1084 *
1085 * where L(k) and L(k+1) are the k-th and (k+1)-th columns
1086 * of L
1087 *
1088 * (1) Store L(k+2:N,k) and L(k+2:N,k+1) and 2-by-2
1089 * block D(k:k+1,k:k+1) in columns k and k+1 of A.
1090 * NOTE: 2-by-2 diagonal block L(k:k+1,k:k+1) is a UNIT
1091 * block and not stored.
1092 * A(k:k+1,k:k+1) := D(k:k+1,k:k+1) = W(k:k+1,k:k+1)
1093 * A(k+2:N,k:k+1) := L(k+2:N,k:k+1) =
1094 * = W(k+2:N,k:k+1) * ( D(k:k+1,k:k+1)**(-1) )
1095 *
1096  IF( k.LT.n-1 ) THEN
1097 *
1098 * Factor out the columns of the inverse of 2-by-2 pivot
1099 * block D, so that each column contains 1, to reduce the
1100 * number of FLOPS when we multiply panel
1101 * ( W(kw-1) W(kw) ) by this inverse, i.e. by D**(-1).
1102 *
1103 * D**(-1) = ( d11 cj(d21) )**(-1) =
1104 * ( d21 d22 )
1105 *
1106 * = 1/(d11*d22-|d21|**2) * ( ( d22) (-cj(d21) ) ) =
1107 * ( (-d21) ( d11 ) )
1108 *
1109 * = 1/(|d21|**2) * 1/((d11/cj(d21))*(d22/d21)-1) *
1110 *
1111 * * ( d21*( d22/d21 ) conj(d21)*( - 1 ) ) =
1112 * ( ( -1 ) ( d11/conj(d21) ) )
1113 *
1114 * = 1/(|d21|**2) * 1/(D22*D11-1) *
1115 *
1116 * * ( d21*( D11 ) conj(d21)*( -1 ) ) =
1117 * ( ( -1 ) ( D22 ) )
1118 *
1119 * = (1/|d21|**2) * T * ( d21*( D11 ) conj(d21)*( -1 ) ) =
1120 * ( ( -1 ) ( D22 ) )
1121 *
1122 * = ( (T/conj(d21))*( D11 ) (T/d21)*( -1 ) ) =
1123 * ( ( -1 ) ( D22 ) )
1124 *
1125 * Handle division by a small number. (NOTE: order of
1126 * operations is important)
1127 *
1128 * = ( T*(( D11 )/conj(D21)) T*(( -1 )/D21 ) )
1129 * ( (( -1 ) ) (( D22 ) ) ),
1130 *
1131 * where D11 = d22/d21,
1132 * D22 = d11/conj(d21),
1133 * D21 = d21,
1134 * T = 1/(D22*D11-1).
1135 *
1136 * (NOTE: No need to check for division by ZERO,
1137 * since that was ensured earlier in pivot search:
1138 * (a) d21 != 0 in 2x2 pivot case(4),
1139 * since |d21| should be larger than |d11| and |d22|;
1140 * (b) (D22*D11 - 1) != 0, since from (a),
1141 * both |D11| < 1, |D22| < 1, hence |D22*D11| << 1.)
1142 *
1143  d21 = w( k+1, k )
1144  d11 = w( k+1, k+1 ) / d21
1145  d22 = w( k, k ) / dconjg( d21 )
1146  t = one / ( dble( d11*d22 )-one )
1147 *
1148 * Update elements in columns A(k) and A(k+1) as
1149 * dot products of rows of ( W(k) W(k+1) ) and columns
1150 * of D**(-1)
1151 *
1152  DO 80 j = k + 2, n
1153  a( j, k ) = t*( ( d11*w( j, k )-w( j, k+1 ) ) /
1154  $ dconjg( d21 ) )
1155  a( j, k+1 ) = t*( ( d22*w( j, k+1 )-w( j, k ) ) /
1156  $ d21 )
1157  80 CONTINUE
1158  END IF
1159 *
1160 * Copy diagonal elements of D(K) to A,
1161 * copy subdiagonal element of D(K) to E(K) and
1162 * ZERO out subdiagonal entry of A
1163 *
1164  a( k, k ) = w( k, k )
1165  a( k+1, k ) = czero
1166  a( k+1, k+1 ) = w( k+1, k+1 )
1167  e( k ) = w( k+1, k )
1168  e( k+1 ) = czero
1169 *
1170 * (2) Conjugate columns W(k) and W(k+1)
1171 *
1172  CALL zlacgv( n-k, w( k+1, k ), 1 )
1173  CALL zlacgv( n-k-1, w( k+2, k+1 ), 1 )
1174 *
1175  END IF
1176 *
1177 * End column K is nonsingular
1178 *
1179  END IF
1180 *
1181 * Store details of the interchanges in IPIV
1182 *
1183  IF( kstep.EQ.1 ) THEN
1184  ipiv( k ) = kp
1185  ELSE
1186  ipiv( k ) = -p
1187  ipiv( k+1 ) = -kp
1188  END IF
1189 *
1190 * Increase K and return to the start of the main loop
1191 *
1192  k = k + kstep
1193  GO TO 70
1194 *
1195  90 CONTINUE
1196 *
1197 * Update the lower triangle of A22 (= A(k:n,k:n)) as
1198 *
1199 * A22 := A22 - L21*D*L21**H = A22 - L21*W**H
1200 *
1201 * computing blocks of NB columns at a time (note that conjg(W) is
1202 * actually stored)
1203 *
1204  DO 110 j = k, n, nb
1205  jb = min( nb, n-j+1 )
1206 *
1207 * Update the lower triangle of the diagonal block
1208 *
1209  DO 100 jj = j, j + jb - 1
1210  a( jj, jj ) = dble( a( jj, jj ) )
1211  CALL zgemv( 'No transpose', j+jb-jj, k-1, -cone,
1212  $ a( jj, 1 ), lda, w( jj, 1 ), ldw, cone,
1213  $ a( jj, jj ), 1 )
1214  a( jj, jj ) = dble( a( jj, jj ) )
1215  100 CONTINUE
1216 *
1217 * Update the rectangular subdiagonal block
1218 *
1219  IF( j+jb.LE.n )
1220  $ CALL zgemm( 'No transpose', 'Transpose', n-j-jb+1, jb,
1221  $ k-1, -cone, a( j+jb, 1 ), lda, w( j, 1 ),
1222  $ ldw, cone, a( j+jb, j ), lda )
1223  110 CONTINUE
1224 *
1225 * Set KB to the number of columns factorized
1226 *
1227  kb = k - 1
1228 *
1229  END IF
1230  RETURN
1231 *
1232 * End of ZLAHEF_RK
1233 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:52
subroutine zgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
ZGEMV
Definition: zgemv.f:160
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
Definition: zswap.f:52
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
Definition: zgemm.f:189
integer function izamax(N, ZX, INCX)
IZAMAX
Definition: izamax.f:53
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine zlacgv(N, X, INCX)
ZLACGV conjugates a complex vector.
Definition: zlacgv.f:76
subroutine zdscal(N, DA, ZX, INCX)
ZDSCAL
Definition: zdscal.f:54

Here is the call graph for this function:

Here is the caller graph for this function: