LAPACK 3.12.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches

## ◆ 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*
subroutine ccopy(n, cx, incx, cy, incy)
CCOPY
Definition ccopy.f:81
subroutine cgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
CGEMM
Definition cgemm.f:188
subroutine cgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
CGEMV
Definition cgemv.f:160
integer function icamax(n, cx, incx)
ICAMAX
Definition icamax.f:71
subroutine clacgv(n, x, incx)
CLACGV conjugates a complex vector.
Definition clacgv.f:74
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
subroutine csscal(n, sa, cx, incx)
CSSCAL
Definition csscal.f:78
subroutine cswap(n, cx, incx, cy, incy)
CSWAP
Definition cswap.f:81
Here is the call graph for this function:
Here is the caller graph for this function: