LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ zlahef_rk()

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.
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 260 of file zlahef_rk.f.

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