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

◆ 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*
subroutine zcopy(n, zx, incx, zy, incy)
ZCOPY
Definition zcopy.f:81
subroutine zgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
ZGEMM
Definition zgemm.f:188
subroutine zgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
ZGEMV
Definition zgemv.f:160
integer function izamax(n, zx, incx)
IZAMAX
Definition izamax.f:71
subroutine zlacgv(n, x, incx)
ZLACGV conjugates a complex vector.
Definition zlacgv.f:74
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
subroutine zdscal(n, da, zx, incx)
ZDSCAL
Definition zdscal.f:78
subroutine zswap(n, zx, incx, zy, incy)
ZSWAP
Definition zswap.f:81
Here is the call graph for this function:
Here is the caller graph for this function: