File:  [local] / rpl / lapack / lapack / zhetf2.f
Revision 1.10: download - view: text, annotated - select for diffs - revision graph
Mon Nov 21 22:19:48 2011 UTC (12 years, 5 months ago) by bertrand
Branches: MAIN
CVS tags: rpl-4_1_8, rpl-4_1_7, rpl-4_1_6, rpl-4_1_5, rpl-4_1_4, HEAD
Cohérence

    1: *> \brief \b ZHETF2
    2: *
    3: *  =========== DOCUMENTATION ===========
    4: *
    5: * Online html documentation available at 
    6: *            http://www.netlib.org/lapack/explore-html/ 
    7: *
    8: *> \htmlonly
    9: *> Download ZHETF2 + dependencies 
   10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhetf2.f"> 
   11: *> [TGZ]</a> 
   12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhetf2.f"> 
   13: *> [ZIP]</a> 
   14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhetf2.f"> 
   15: *> [TXT]</a>
   16: *> \endhtmlonly 
   17: *
   18: *  Definition:
   19: *  ===========
   20: *
   21: *       SUBROUTINE ZHETF2( UPLO, N, A, LDA, IPIV, INFO )
   22:    23: *       .. Scalar Arguments ..
   24: *       CHARACTER          UPLO
   25: *       INTEGER            INFO, LDA, N
   26: *       ..
   27: *       .. Array Arguments ..
   28: *       INTEGER            IPIV( * )
   29: *       COMPLEX*16         A( LDA, * )
   30: *       ..
   31: *  
   32: *
   33: *> \par Purpose:
   34: *  =============
   35: *>
   36: *> \verbatim
   37: *>
   38: *> ZHETF2 computes the factorization of a complex Hermitian matrix A
   39: *> using the Bunch-Kaufman diagonal pivoting method:
   40: *>
   41: *>    A = U*D*U**H  or  A = L*D*L**H
   42: *>
   43: *> where U (or L) is a product of permutation and unit upper (lower)
   44: *> triangular matrices, U**H is the conjugate transpose of U, and D is
   45: *> Hermitian and block diagonal with 1-by-1 and 2-by-2 diagonal blocks.
   46: *>
   47: *> This is the unblocked version of the algorithm, calling Level 2 BLAS.
   48: *> \endverbatim
   49: *
   50: *  Arguments:
   51: *  ==========
   52: *
   53: *> \param[in] UPLO
   54: *> \verbatim
   55: *>          UPLO is CHARACTER*1
   56: *>          Specifies whether the upper or lower triangular part of the
   57: *>          Hermitian matrix A is stored:
   58: *>          = 'U':  Upper triangular
   59: *>          = 'L':  Lower triangular
   60: *> \endverbatim
   61: *>
   62: *> \param[in] N
   63: *> \verbatim
   64: *>          N is INTEGER
   65: *>          The order of the matrix A.  N >= 0.
   66: *> \endverbatim
   67: *>
   68: *> \param[in,out] A
   69: *> \verbatim
   70: *>          A is COMPLEX*16 array, dimension (LDA,N)
   71: *>          On entry, the Hermitian matrix A.  If UPLO = 'U', the leading
   72: *>          n-by-n upper triangular part of A contains the upper
   73: *>          triangular part of the matrix A, and the strictly lower
   74: *>          triangular part of A is not referenced.  If UPLO = 'L', the
   75: *>          leading n-by-n lower triangular part of A contains the lower
   76: *>          triangular part of the matrix A, and the strictly upper
   77: *>          triangular part of A is not referenced.
   78: *>
   79: *>          On exit, the block diagonal matrix D and the multipliers used
   80: *>          to obtain the factor U or L (see below for further details).
   81: *> \endverbatim
   82: *>
   83: *> \param[in] LDA
   84: *> \verbatim
   85: *>          LDA is INTEGER
   86: *>          The leading dimension of the array A.  LDA >= max(1,N).
   87: *> \endverbatim
   88: *>
   89: *> \param[out] IPIV
   90: *> \verbatim
   91: *>          IPIV is INTEGER array, dimension (N)
   92: *>          Details of the interchanges and the block structure of D.
   93: *>          If IPIV(k) > 0, then rows and columns k and IPIV(k) were
   94: *>          interchanged and D(k,k) is a 1-by-1 diagonal block.
   95: *>          If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and
   96: *>          columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
   97: *>          is a 2-by-2 diagonal block.  If UPLO = 'L' and IPIV(k) =
   98: *>          IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were
   99: *>          interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
  100: *> \endverbatim
  101: *>
  102: *> \param[out] INFO
  103: *> \verbatim
  104: *>          INFO is INTEGER
  105: *>          = 0: successful exit
  106: *>          < 0: if INFO = -k, the k-th argument had an illegal value
  107: *>          > 0: if INFO = k, D(k,k) is exactly zero.  The factorization
  108: *>               has been completed, but the block diagonal matrix D is
  109: *>               exactly singular, and division by zero will occur if it
  110: *>               is used to solve a system of equations.
  111: *> \endverbatim
  112: *
  113: *  Authors:
  114: *  ========
  115: *
  116: *> \author Univ. of Tennessee 
  117: *> \author Univ. of California Berkeley 
  118: *> \author Univ. of Colorado Denver 
  119: *> \author NAG Ltd. 
  120: *
  121: *> \date November 2011
  122: *
  123: *> \ingroup complex16HEcomputational
  124: *
  125: *> \par Further Details:
  126: *  =====================
  127: *>
  128: *> \verbatim
  129: *>
  130: *>  If UPLO = 'U', then A = U*D*U**H, where
  131: *>     U = P(n)*U(n)* ... *P(k)U(k)* ...,
  132: *>  i.e., U is a product of terms P(k)*U(k), where k decreases from n to
  133: *>  1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
  134: *>  and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
  135: *>  defined by IPIV(k), and U(k) is a unit upper triangular matrix, such
  136: *>  that if the diagonal block D(k) is of order s (s = 1 or 2), then
  137: *>
  138: *>             (   I    v    0   )   k-s
  139: *>     U(k) =  (   0    I    0   )   s
  140: *>             (   0    0    I   )   n-k
  141: *>                k-s   s   n-k
  142: *>
  143: *>  If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k).
  144: *>  If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k),
  145: *>  and A(k,k), and v overwrites A(1:k-2,k-1:k).
  146: *>
  147: *>  If UPLO = 'L', then A = L*D*L**H, where
  148: *>     L = P(1)*L(1)* ... *P(k)*L(k)* ...,
  149: *>  i.e., L is a product of terms P(k)*L(k), where k increases from 1 to
  150: *>  n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
  151: *>  and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
  152: *>  defined by IPIV(k), and L(k) is a unit lower triangular matrix, such
  153: *>  that if the diagonal block D(k) is of order s (s = 1 or 2), then
  154: *>
  155: *>             (   I    0     0   )  k-1
  156: *>     L(k) =  (   0    I     0   )  s
  157: *>             (   0    v     I   )  n-k-s+1
  158: *>                k-1   s  n-k-s+1
  159: *>
  160: *>  If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k).
  161: *>  If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k),
  162: *>  and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1).
  163: *> \endverbatim
  164: *
  165: *> \par Contributors:
  166: *  ==================
  167: *>
  168: *> \verbatim
  169: *>  09-29-06 - patch from
  170: *>    Bobby Cheng, MathWorks
  171: *>
  172: *>    Replace l.210 and l.393
  173: *>         IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
  174: *>    by
  175: *>         IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. DISNAN(ABSAKK) ) THEN
  176: *>
  177: *>  01-01-96 - Based on modifications by
  178: *>    J. Lewis, Boeing Computer Services Company
  179: *>    A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA
  180: *> \endverbatim
  181: *
  182: *  =====================================================================
  183:       SUBROUTINE ZHETF2( UPLO, N, A, LDA, IPIV, INFO )
  184: *
  185: *  -- LAPACK computational routine (version 3.4.0) --
  186: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  187: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  188: *     November 2011
  189: *
  190: *     .. Scalar Arguments ..
  191:       CHARACTER          UPLO
  192:       INTEGER            INFO, LDA, N
  193: *     ..
  194: *     .. Array Arguments ..
  195:       INTEGER            IPIV( * )
  196:       COMPLEX*16         A( LDA, * )
  197: *     ..
  198: *
  199: *  =====================================================================
  200: *
  201: *     .. Parameters ..
  202:       DOUBLE PRECISION   ZERO, ONE
  203:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
  204:       DOUBLE PRECISION   EIGHT, SEVTEN
  205:       PARAMETER          ( EIGHT = 8.0D+0, SEVTEN = 17.0D+0 )
  206: *     ..
  207: *     .. Local Scalars ..
  208:       LOGICAL            UPPER
  209:       INTEGER            I, IMAX, J, JMAX, K, KK, KP, KSTEP
  210:       DOUBLE PRECISION   ABSAKK, ALPHA, COLMAX, D, D11, D22, R1, ROWMAX,
  211:      $                   TT
  212:       COMPLEX*16         D12, D21, T, WK, WKM1, WKP1, ZDUM
  213: *     ..
  214: *     .. External Functions ..
  215:       LOGICAL            LSAME, DISNAN
  216:       INTEGER            IZAMAX
  217:       DOUBLE PRECISION   DLAPY2
  218:       EXTERNAL           LSAME, IZAMAX, DLAPY2, DISNAN
  219: *     ..
  220: *     .. External Subroutines ..
  221:       EXTERNAL           XERBLA, ZDSCAL, ZHER, ZSWAP
  222: *     ..
  223: *     .. Intrinsic Functions ..
  224:       INTRINSIC          ABS, DBLE, DCMPLX, DCONJG, DIMAG, MAX, SQRT
  225: *     ..
  226: *     .. Statement Functions ..
  227:       DOUBLE PRECISION   CABS1
  228: *     ..
  229: *     .. Statement Function definitions ..
  230:       CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
  231: *     ..
  232: *     .. Executable Statements ..
  233: *
  234: *     Test the input parameters.
  235: *
  236:       INFO = 0
  237:       UPPER = LSAME( UPLO, 'U' )
  238:       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
  239:          INFO = -1
  240:       ELSE IF( N.LT.0 ) THEN
  241:          INFO = -2
  242:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
  243:          INFO = -4
  244:       END IF
  245:       IF( INFO.NE.0 ) THEN
  246:          CALL XERBLA( 'ZHETF2', -INFO )
  247:          RETURN
  248:       END IF
  249: *
  250: *     Initialize ALPHA for use in choosing pivot block size.
  251: *
  252:       ALPHA = ( ONE+SQRT( SEVTEN ) ) / EIGHT
  253: *
  254:       IF( UPPER ) THEN
  255: *
  256: *        Factorize A as U*D*U**H using the upper triangle of A
  257: *
  258: *        K is the main loop index, decreasing from N to 1 in steps of
  259: *        1 or 2
  260: *
  261:          K = N
  262:    10    CONTINUE
  263: *
  264: *        If K < 1, exit from loop
  265: *
  266:          IF( K.LT.1 )
  267:      $      GO TO 90
  268:          KSTEP = 1
  269: *
  270: *        Determine rows and columns to be interchanged and whether
  271: *        a 1-by-1 or 2-by-2 pivot block will be used
  272: *
  273:          ABSAKK = ABS( DBLE( A( K, K ) ) )
  274: *
  275: *        IMAX is the row-index of the largest off-diagonal element in
  276: *        column K, and COLMAX is its absolute value
  277: *
  278:          IF( K.GT.1 ) THEN
  279:             IMAX = IZAMAX( K-1, A( 1, K ), 1 )
  280:             COLMAX = CABS1( A( IMAX, K ) )
  281:          ELSE
  282:             COLMAX = ZERO
  283:          END IF
  284: *
  285:          IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. DISNAN(ABSAKK) ) THEN
  286: *
  287: *           Column K is zero or contains a NaN: set INFO and continue
  288: *
  289:             IF( INFO.EQ.0 )
  290:      $         INFO = K
  291:             KP = K
  292:             A( K, K ) = DBLE( A( K, K ) )
  293:          ELSE
  294:             IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
  295: *
  296: *              no interchange, use 1-by-1 pivot block
  297: *
  298:                KP = K
  299:             ELSE
  300: *
  301: *              JMAX is the column-index of the largest off-diagonal
  302: *              element in row IMAX, and ROWMAX is its absolute value
  303: *
  304:                JMAX = IMAX + IZAMAX( K-IMAX, A( IMAX, IMAX+1 ), LDA )
  305:                ROWMAX = CABS1( A( IMAX, JMAX ) )
  306:                IF( IMAX.GT.1 ) THEN
  307:                   JMAX = IZAMAX( IMAX-1, A( 1, IMAX ), 1 )
  308:                   ROWMAX = MAX( ROWMAX, CABS1( A( JMAX, IMAX ) ) )
  309:                END IF
  310: *
  311:                IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
  312: *
  313: *                 no interchange, use 1-by-1 pivot block
  314: *
  315:                   KP = K
  316:                ELSE IF( ABS( DBLE( A( IMAX, IMAX ) ) ).GE.ALPHA*ROWMAX )
  317:      $                   THEN
  318: *
  319: *                 interchange rows and columns K and IMAX, use 1-by-1
  320: *                 pivot block
  321: *
  322:                   KP = IMAX
  323:                ELSE
  324: *
  325: *                 interchange rows and columns K-1 and IMAX, use 2-by-2
  326: *                 pivot block
  327: *
  328:                   KP = IMAX
  329:                   KSTEP = 2
  330:                END IF
  331:             END IF
  332: *
  333:             KK = K - KSTEP + 1
  334:             IF( KP.NE.KK ) THEN
  335: *
  336: *              Interchange rows and columns KK and KP in the leading
  337: *              submatrix A(1:k,1:k)
  338: *
  339:                CALL ZSWAP( KP-1, A( 1, KK ), 1, A( 1, KP ), 1 )
  340:                DO 20 J = KP + 1, KK - 1
  341:                   T = DCONJG( A( J, KK ) )
  342:                   A( J, KK ) = DCONJG( A( KP, J ) )
  343:                   A( KP, J ) = T
  344:    20          CONTINUE
  345:                A( KP, KK ) = DCONJG( A( KP, KK ) )
  346:                R1 = DBLE( A( KK, KK ) )
  347:                A( KK, KK ) = DBLE( A( KP, KP ) )
  348:                A( KP, KP ) = R1
  349:                IF( KSTEP.EQ.2 ) THEN
  350:                   A( K, K ) = DBLE( A( K, K ) )
  351:                   T = A( K-1, K )
  352:                   A( K-1, K ) = A( KP, K )
  353:                   A( KP, K ) = T
  354:                END IF
  355:             ELSE
  356:                A( K, K ) = DBLE( A( K, K ) )
  357:                IF( KSTEP.EQ.2 )
  358:      $            A( K-1, K-1 ) = DBLE( A( K-1, K-1 ) )
  359:             END IF
  360: *
  361: *           Update the leading submatrix
  362: *
  363:             IF( KSTEP.EQ.1 ) THEN
  364: *
  365: *              1-by-1 pivot block D(k): column k now holds
  366: *
  367: *              W(k) = U(k)*D(k)
  368: *
  369: *              where U(k) is the k-th column of U
  370: *
  371: *              Perform a rank-1 update of A(1:k-1,1:k-1) as
  372: *
  373: *              A := A - U(k)*D(k)*U(k)**H = A - W(k)*1/D(k)*W(k)**H
  374: *
  375:                R1 = ONE / DBLE( A( K, K ) )
  376:                CALL ZHER( UPLO, K-1, -R1, A( 1, K ), 1, A, LDA )
  377: *
  378: *              Store U(k) in column k
  379: *
  380:                CALL ZDSCAL( K-1, R1, A( 1, K ), 1 )
  381:             ELSE
  382: *
  383: *              2-by-2 pivot block D(k): columns k and k-1 now hold
  384: *
  385: *              ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
  386: *
  387: *              where U(k) and U(k-1) are the k-th and (k-1)-th columns
  388: *              of U
  389: *
  390: *              Perform a rank-2 update of A(1:k-2,1:k-2) as
  391: *
  392: *              A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )**H
  393: *                 = A - ( W(k-1) W(k) )*inv(D(k))*( W(k-1) W(k) )**H
  394: *
  395:                IF( K.GT.2 ) THEN
  396: *
  397:                   D = DLAPY2( DBLE( A( K-1, K ) ),
  398:      $                DIMAG( A( K-1, K ) ) )
  399:                   D22 = DBLE( A( K-1, K-1 ) ) / D
  400:                   D11 = DBLE( A( K, K ) ) / D
  401:                   TT = ONE / ( D11*D22-ONE )
  402:                   D12 = A( K-1, K ) / D
  403:                   D = TT / D
  404: *
  405:                   DO 40 J = K - 2, 1, -1
  406:                      WKM1 = D*( D11*A( J, K-1 )-DCONJG( D12 )*
  407:      $                      A( J, K ) )
  408:                      WK = D*( D22*A( J, K )-D12*A( J, K-1 ) )
  409:                      DO 30 I = J, 1, -1
  410:                         A( I, J ) = A( I, J ) - A( I, K )*DCONJG( WK ) -
  411:      $                              A( I, K-1 )*DCONJG( WKM1 )
  412:    30                CONTINUE
  413:                      A( J, K ) = WK
  414:                      A( J, K-1 ) = WKM1
  415:                      A( J, J ) = DCMPLX( DBLE( A( J, J ) ), 0.0D+0 )
  416:    40             CONTINUE
  417: *
  418:                END IF
  419: *
  420:             END IF
  421:          END IF
  422: *
  423: *        Store details of the interchanges in IPIV
  424: *
  425:          IF( KSTEP.EQ.1 ) THEN
  426:             IPIV( K ) = KP
  427:          ELSE
  428:             IPIV( K ) = -KP
  429:             IPIV( K-1 ) = -KP
  430:          END IF
  431: *
  432: *        Decrease K and return to the start of the main loop
  433: *
  434:          K = K - KSTEP
  435:          GO TO 10
  436: *
  437:       ELSE
  438: *
  439: *        Factorize A as L*D*L**H using the lower triangle of A
  440: *
  441: *        K is the main loop index, increasing from 1 to N in steps of
  442: *        1 or 2
  443: *
  444:          K = 1
  445:    50    CONTINUE
  446: *
  447: *        If K > N, exit from loop
  448: *
  449:          IF( K.GT.N )
  450:      $      GO TO 90
  451:          KSTEP = 1
  452: *
  453: *        Determine rows and columns to be interchanged and whether
  454: *        a 1-by-1 or 2-by-2 pivot block will be used
  455: *
  456:          ABSAKK = ABS( DBLE( A( K, K ) ) )
  457: *
  458: *        IMAX is the row-index of the largest off-diagonal element in
  459: *        column K, and COLMAX is its absolute value
  460: *
  461:          IF( K.LT.N ) THEN
  462:             IMAX = K + IZAMAX( N-K, A( K+1, K ), 1 )
  463:             COLMAX = CABS1( A( IMAX, K ) )
  464:          ELSE
  465:             COLMAX = ZERO
  466:          END IF
  467: *
  468:          IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. DISNAN(ABSAKK) ) THEN
  469: *
  470: *           Column K is zero or contains a NaN: set INFO and continue
  471: *
  472:             IF( INFO.EQ.0 )
  473:      $         INFO = K
  474:             KP = K
  475:             A( K, K ) = DBLE( A( K, K ) )
  476:          ELSE
  477:             IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
  478: *
  479: *              no interchange, use 1-by-1 pivot block
  480: *
  481:                KP = K
  482:             ELSE
  483: *
  484: *              JMAX is the column-index of the largest off-diagonal
  485: *              element in row IMAX, and ROWMAX is its absolute value
  486: *
  487:                JMAX = K - 1 + IZAMAX( IMAX-K, A( IMAX, K ), LDA )
  488:                ROWMAX = CABS1( A( IMAX, JMAX ) )
  489:                IF( IMAX.LT.N ) THEN
  490:                   JMAX = IMAX + IZAMAX( N-IMAX, A( IMAX+1, IMAX ), 1 )
  491:                   ROWMAX = MAX( ROWMAX, CABS1( A( JMAX, IMAX ) ) )
  492:                END IF
  493: *
  494:                IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
  495: *
  496: *                 no interchange, use 1-by-1 pivot block
  497: *
  498:                   KP = K
  499:                ELSE IF( ABS( DBLE( A( IMAX, IMAX ) ) ).GE.ALPHA*ROWMAX )
  500:      $                   THEN
  501: *
  502: *                 interchange rows and columns K and IMAX, use 1-by-1
  503: *                 pivot block
  504: *
  505:                   KP = IMAX
  506:                ELSE
  507: *
  508: *                 interchange rows and columns K+1 and IMAX, use 2-by-2
  509: *                 pivot block
  510: *
  511:                   KP = IMAX
  512:                   KSTEP = 2
  513:                END IF
  514:             END IF
  515: *
  516:             KK = K + KSTEP - 1
  517:             IF( KP.NE.KK ) THEN
  518: *
  519: *              Interchange rows and columns KK and KP in the trailing
  520: *              submatrix A(k:n,k:n)
  521: *
  522:                IF( KP.LT.N )
  523:      $            CALL ZSWAP( N-KP, A( KP+1, KK ), 1, A( KP+1, KP ), 1 )
  524:                DO 60 J = KK + 1, KP - 1
  525:                   T = DCONJG( A( J, KK ) )
  526:                   A( J, KK ) = DCONJG( A( KP, J ) )
  527:                   A( KP, J ) = T
  528:    60          CONTINUE
  529:                A( KP, KK ) = DCONJG( A( KP, KK ) )
  530:                R1 = DBLE( A( KK, KK ) )
  531:                A( KK, KK ) = DBLE( A( KP, KP ) )
  532:                A( KP, KP ) = R1
  533:                IF( KSTEP.EQ.2 ) THEN
  534:                   A( K, K ) = DBLE( A( K, K ) )
  535:                   T = A( K+1, K )
  536:                   A( K+1, K ) = A( KP, K )
  537:                   A( KP, K ) = T
  538:                END IF
  539:             ELSE
  540:                A( K, K ) = DBLE( A( K, K ) )
  541:                IF( KSTEP.EQ.2 )
  542:      $            A( K+1, K+1 ) = DBLE( A( K+1, K+1 ) )
  543:             END IF
  544: *
  545: *           Update the trailing submatrix
  546: *
  547:             IF( KSTEP.EQ.1 ) THEN
  548: *
  549: *              1-by-1 pivot block D(k): column k now holds
  550: *
  551: *              W(k) = L(k)*D(k)
  552: *
  553: *              where L(k) is the k-th column of L
  554: *
  555:                IF( K.LT.N ) THEN
  556: *
  557: *                 Perform a rank-1 update of A(k+1:n,k+1:n) as
  558: *
  559: *                 A := A - L(k)*D(k)*L(k)**H = A - W(k)*(1/D(k))*W(k)**H
  560: *
  561:                   R1 = ONE / DBLE( A( K, K ) )
  562:                   CALL ZHER( UPLO, N-K, -R1, A( K+1, K ), 1,
  563:      $                       A( K+1, K+1 ), LDA )
  564: *
  565: *                 Store L(k) in column K
  566: *
  567:                   CALL ZDSCAL( N-K, R1, A( K+1, K ), 1 )
  568:                END IF
  569:             ELSE
  570: *
  571: *              2-by-2 pivot block D(k)
  572: *
  573:                IF( K.LT.N-1 ) THEN
  574: *
  575: *                 Perform a rank-2 update of A(k+2:n,k+2:n) as
  576: *
  577: *                 A := A - ( L(k) L(k+1) )*D(k)*( L(k) L(k+1) )**H
  578: *                    = A - ( W(k) W(k+1) )*inv(D(k))*( W(k) W(k+1) )**H
  579: *
  580: *                 where L(k) and L(k+1) are the k-th and (k+1)-th
  581: *                 columns of L
  582: *
  583:                   D = DLAPY2( DBLE( A( K+1, K ) ),
  584:      $                DIMAG( A( K+1, K ) ) )
  585:                   D11 = DBLE( A( K+1, K+1 ) ) / D
  586:                   D22 = DBLE( A( K, K ) ) / D
  587:                   TT = ONE / ( D11*D22-ONE )
  588:                   D21 = A( K+1, K ) / D
  589:                   D = TT / D
  590: *
  591:                   DO 80 J = K + 2, N
  592:                      WK = D*( D11*A( J, K )-D21*A( J, K+1 ) )
  593:                      WKP1 = D*( D22*A( J, K+1 )-DCONJG( D21 )*
  594:      $                      A( J, K ) )
  595:                      DO 70 I = J, N
  596:                         A( I, J ) = A( I, J ) - A( I, K )*DCONJG( WK ) -
  597:      $                              A( I, K+1 )*DCONJG( WKP1 )
  598:    70                CONTINUE
  599:                      A( J, K ) = WK
  600:                      A( J, K+1 ) = WKP1
  601:                      A( J, J ) = DCMPLX( DBLE( A( J, J ) ), 0.0D+0 )
  602:    80             CONTINUE
  603:                END IF
  604:             END IF
  605:          END IF
  606: *
  607: *        Store details of the interchanges in IPIV
  608: *
  609:          IF( KSTEP.EQ.1 ) THEN
  610:             IPIV( K ) = KP
  611:          ELSE
  612:             IPIV( K ) = -KP
  613:             IPIV( K+1 ) = -KP
  614:          END IF
  615: *
  616: *        Increase K and return to the start of the main loop
  617: *
  618:          K = K + KSTEP
  619:          GO TO 50
  620: *
  621:       END IF
  622: *
  623:    90 CONTINUE
  624:       RETURN
  625: *
  626: *     End of ZHETF2
  627: *
  628:       END

CVSweb interface <joel.bertrand@systella.fr>