 LAPACK  3.10.1 LAPACK: Linear Algebra PACKage

## ◆ clahef_rk()

 subroutine clahef_rk ( character UPLO, integer N, integer NB, integer KB, complex, dimension( lda, * ) A, integer LDA, complex, dimension( * ) E, integer, dimension( * ) IPIV, complex, dimension( ldw, * ) W, integer LDW, integer INFO )

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

Purpose:
``` CLAHEF_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.

CLAHEF_RK is an auxiliary routine called by CHETRF_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 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 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 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.```
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 clahef_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 A( LDA, * ), W( LDW, * ), E( * )
274 * ..
275 *
276 * =====================================================================
277 *
278 * .. Parameters ..
279  REAL ZERO, ONE
280  parameter( zero = 0.0e+0, one = 1.0e+0 )
281  REAL EIGHT, SEVTEN
282  parameter( eight = 8.0e+0, sevten = 17.0e+0 )
283  COMPLEX CONE, CZERO
284  parameter( cone = ( 1.0e+0, 0.0e+0 ),
285  \$ czero = ( 0.0e+0, 0.0e+0 ) )
286 * ..
287 * .. Local Scalars ..
288  LOGICAL DONE
289  INTEGER IMAX, ITEMP, II, J, JB, JJ, JMAX, K, KK, KKW,
290  \$ KP, KSTEP, KW, P
291  REAL ABSAKK, ALPHA, COLMAX, STEMP, R1, ROWMAX, T,
292  \$ SFMIN
293  COMPLEX D11, D21, D22, Z
294 * ..
295 * .. External Functions ..
296  LOGICAL LSAME
297  INTEGER ICAMAX
298  REAL SLAMCH
299  EXTERNAL lsame, icamax, slamch
300 * ..
301 * .. External Subroutines ..
302  EXTERNAL ccopy, csscal, cgemm, cgemv, clacgv, cswap
303 * ..
304 * .. Intrinsic Functions ..
305  INTRINSIC abs, conjg, aimag, max, min, real, sqrt
306 * ..
307 * .. Statement Functions ..
308  REAL CABS1
309 * ..
310 * .. Statement Function definitions ..
311  cabs1( z ) = abs( real( z ) ) + abs( aimag( z ) )
312 * ..
313 * .. Executable Statements ..
314 *
315  info = 0
316 *
317 * Initialize ALPHA for use in choosing pivot block size.
318 *
319  alpha = ( one+sqrt( sevten ) ) / eight
320 *
321 * Compute machine safe minimum
322 *
323  sfmin = slamch( 'S' )
324 *
325  IF( lsame( uplo, 'U' ) ) THEN
326 *
327 * Factorize the trailing columns of A using the upper triangle
328 * of A and working backwards, and compute the matrix W = U12*D
329 * for use in updating A11 (note that conjg(W) is actually stored)
330 *
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 ccopy( k-1, a( 1, k ), 1, w( 1, kw ), 1 )
357  w( k, kw ) = real( a( k, k ) )
358  IF( k.LT.n ) THEN
359  CALL cgemv( '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 ) = real( 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( real( 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 = icamax( 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 ) = real( w( k, kw ) )
388  IF( k.GT.1 )
389  \$ CALL ccopy( 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 ccopy( imax-1, a( 1, imax ), 1, w( 1, kw-1 ),
426  \$ 1 )
427  w( imax, kw-1 ) = real( a( imax, imax ) )
428 *
429  CALL ccopy( k-imax, a( imax, imax+1 ), lda,
430  \$ w( imax+1, kw-1 ), 1 )
431  CALL clacgv( k-imax, w( imax+1, kw-1 ), 1 )
432 *
433  IF( k.LT.n ) THEN
434  CALL cgemv( '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 ) = real( 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 + icamax( 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 = icamax( imax-1, w( 1, kw-1 ), 1 )
454  stemp = cabs1( w( itemp, kw-1 ) )
455  IF( stemp.GT.rowmax ) THEN
456  rowmax = stemp
457  jmax = itemp
458  END IF
459  END IF
460 *
461 * Case(2)
462 * Equivalent to testing for
463 * ABS( REAL( W( IMAX,KW-1 ) ) ).GE.ALPHA*ROWMAX
464 * (used to handle NaN and Inf)
465 *
466  IF( .NOT.( abs( real( 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 ccopy( 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 *
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 ccopy( 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 ) = real( a( k, k ) )
539  CALL ccopy( k-1-p, a( p+1, k ), 1, a( p, p+1 ),
540  \$ lda )
541  CALL clacgv( k-1-p, a( p, p+1 ), lda )
542  IF( p.GT.1 )
543  \$ CALL ccopy( 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 cswap( n-k, a( k, k+1 ), lda, a( p, k+1 ),
552  \$ lda )
553  CALL cswap( 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 ) = real( a( kk, kk ) )
568  CALL ccopy( kk-1-kp, a( kp+1, kk ), 1, a( kp, kp+1 ),
569  \$ lda )
570  CALL clacgv( kk-1-kp, a( kp, kp+1 ), lda )
571  IF( kp.GT.1 )
572  \$ CALL ccopy( 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 cswap( n-k, a( kk, k+1 ), lda, a( kp, k+1 ),
581  \$ lda )
582  CALL cswap( 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 ) = REAL( W( K, K) ) to separately copy diagonal
603 * element D(k,k) from W (potentially saves only one load))
604  CALL ccopy( 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 = real( a( k, k ) )
614  IF( abs( t ).GE.sfmin ) THEN
615  r1 = one / t
616  CALL csscal( 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 clacgv( 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 ) / conjg( d21 )
699  d22 = w( k-1, kw-1 ) / d21
700  t = one / ( real( 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  \$ conjg( 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 clacgv( k-1, w( 1, kw ), 1 )
727  CALL clacgv( 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 ) = real( a( jj, jj ) )
765  CALL cgemv( '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 ) = real( a( jj, jj ) )
769  40 CONTINUE
770 *
771 * Update the rectangular superdiagonal block
772 *
773  IF( j.GE.2 )
774  \$ CALL cgemm( '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 ) = real( a( k, k ) )
809  IF( k.LT.n )
810  \$ CALL ccopy( n-k, a( k+1, k ), 1, w( k+1, k ), 1 )
811  IF( k.GT.1 ) THEN
812  CALL cgemv( '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 ) = real( 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( real( 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 + icamax( 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 ) = real( w( k, k ) )
841  IF( k.LT.n )
842  \$ CALL ccopy( 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 ccopy( imax-k, a( imax, k ), lda, w( k, k+1 ), 1)
879  CALL clacgv( imax-k, w( k, k+1 ), 1 )
880  w( imax, k+1 ) = real( a( imax, imax ) )
881 *
882  IF( imax.LT.n )
883  \$ CALL ccopy( n-imax, a( imax+1, imax ), 1,
884  \$ w( imax+1, k+1 ), 1 )
885 *
886  IF( k.GT.1 ) THEN
887  CALL cgemv( '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 ) = real( 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 + icamax( 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 + icamax( n-imax, w( imax+1, k+1 ), 1)
906  stemp = cabs1( w( itemp, k+1 ) )
907  IF( stemp.GT.rowmax ) THEN
908  rowmax = stemp
909  jmax = itemp
910  END IF
911  END IF
912 *
913 * Case(2)
914 * Equivalent to testing for
915 * ABS( REAL( W( IMAX,K+1 ) ) ).GE.ALPHA*ROWMAX
916 * (used to handle NaN and Inf)
917 *
918  IF( .NOT.( abs( real( 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 ccopy( 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 *
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 ccopy( 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 ) = real( a( k, k ) )
987  CALL ccopy( p-k-1, a( k+1, k ), 1, a( p, k+1 ), lda )
988  CALL clacgv( p-k-1, a( p, k+1 ), lda )
989  IF( p.LT.n )
990  \$ CALL ccopy( 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 cswap( k-1, a( k, 1 ), lda, a( p, 1 ), lda )
999  CALL cswap( 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 ) = real( a( kk, kk ) )
1013  CALL ccopy( kp-kk-1, a( kk+1, kk ), 1, a( kp, kk+1 ),
1014  \$ lda )
1015  CALL clacgv( kp-kk-1, a( kp, kk+1 ), lda )
1016  IF( kp.LT.n )
1017  \$ CALL ccopy( 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 cswap( k-1, a( kk, 1 ), lda, a( kp, 1 ), lda )
1026  CALL cswap( 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 ) = REAL( W( K, K) ) to separately copy diagonal
1046 * element D(k,k) from W (potentially saves only one load))
1047  CALL ccopy( 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 = real( a( k, k ) )
1057  IF( abs( t ).GE.sfmin ) THEN
1058  r1 = one / t
1059  CALL csscal( 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 clacgv( 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 ) / conjg( d21 )
1143  t = one / ( real( 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  \$ conjg( 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 clacgv( n-k, w( k+1, k ), 1 )
1170  CALL clacgv( 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 ) = real( a( jj, jj ) )
1208  CALL cgemv( '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 ) = real( a( jj, jj ) )
1212  100 CONTINUE
1213 *
1214 * Update the rectangular subdiagonal block
1215 *
1216  IF( j+jb.LE.n )
1217  \$ CALL cgemm( '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 CLAHEF_RK
1230 *
integer function icamax(N, CX, INCX)
ICAMAX
Definition: icamax.f:71
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:81
subroutine csscal(N, SA, CX, INCX)
CSSCAL
Definition: csscal.f:78
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
Definition: cswap.f:81
subroutine cgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
CGEMV
Definition: cgemv.f:158
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
Definition: cgemm.f:187
subroutine clacgv(N, X, INCX)
CLACGV conjugates a complex vector.
Definition: clacgv.f:74
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: