LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
zhetf2.f
Go to the documentation of this file.
00001 *> \brief \b ZHETF2
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download ZHETF2 + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhetf2.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhetf2.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhetf2.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE ZHETF2( UPLO, N, A, LDA, IPIV, INFO )
00022 * 
00023 *       .. Scalar Arguments ..
00024 *       CHARACTER          UPLO
00025 *       INTEGER            INFO, LDA, N
00026 *       ..
00027 *       .. Array Arguments ..
00028 *       INTEGER            IPIV( * )
00029 *       COMPLEX*16         A( LDA, * )
00030 *       ..
00031 *  
00032 *
00033 *> \par Purpose:
00034 *  =============
00035 *>
00036 *> \verbatim
00037 *>
00038 *> ZHETF2 computes the factorization of a complex Hermitian matrix A
00039 *> using the Bunch-Kaufman diagonal pivoting method:
00040 *>
00041 *>    A = U*D*U**H  or  A = L*D*L**H
00042 *>
00043 *> where U (or L) is a product of permutation and unit upper (lower)
00044 *> triangular matrices, U**H is the conjugate transpose of U, and D is
00045 *> Hermitian and block diagonal with 1-by-1 and 2-by-2 diagonal blocks.
00046 *>
00047 *> This is the unblocked version of the algorithm, calling Level 2 BLAS.
00048 *> \endverbatim
00049 *
00050 *  Arguments:
00051 *  ==========
00052 *
00053 *> \param[in] UPLO
00054 *> \verbatim
00055 *>          UPLO is CHARACTER*1
00056 *>          Specifies whether the upper or lower triangular part of the
00057 *>          Hermitian matrix A is stored:
00058 *>          = 'U':  Upper triangular
00059 *>          = 'L':  Lower triangular
00060 *> \endverbatim
00061 *>
00062 *> \param[in] N
00063 *> \verbatim
00064 *>          N is INTEGER
00065 *>          The order of the matrix A.  N >= 0.
00066 *> \endverbatim
00067 *>
00068 *> \param[in,out] A
00069 *> \verbatim
00070 *>          A is COMPLEX*16 array, dimension (LDA,N)
00071 *>          On entry, the Hermitian matrix A.  If UPLO = 'U', the leading
00072 *>          n-by-n upper triangular part of A contains the upper
00073 *>          triangular part of the matrix A, and the strictly lower
00074 *>          triangular part of A is not referenced.  If UPLO = 'L', the
00075 *>          leading n-by-n lower triangular part of A contains the lower
00076 *>          triangular part of the matrix A, and the strictly upper
00077 *>          triangular part of A is not referenced.
00078 *>
00079 *>          On exit, the block diagonal matrix D and the multipliers used
00080 *>          to obtain the factor U or L (see below for further details).
00081 *> \endverbatim
00082 *>
00083 *> \param[in] LDA
00084 *> \verbatim
00085 *>          LDA is INTEGER
00086 *>          The leading dimension of the array A.  LDA >= max(1,N).
00087 *> \endverbatim
00088 *>
00089 *> \param[out] IPIV
00090 *> \verbatim
00091 *>          IPIV is INTEGER array, dimension (N)
00092 *>          Details of the interchanges and the block structure of D.
00093 *>          If IPIV(k) > 0, then rows and columns k and IPIV(k) were
00094 *>          interchanged and D(k,k) is a 1-by-1 diagonal block.
00095 *>          If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and
00096 *>          columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
00097 *>          is a 2-by-2 diagonal block.  If UPLO = 'L' and IPIV(k) =
00098 *>          IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were
00099 *>          interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
00100 *> \endverbatim
00101 *>
00102 *> \param[out] INFO
00103 *> \verbatim
00104 *>          INFO is INTEGER
00105 *>          = 0: successful exit
00106 *>          < 0: if INFO = -k, the k-th argument had an illegal value
00107 *>          > 0: if INFO = k, D(k,k) is exactly zero.  The factorization
00108 *>               has been completed, but the block diagonal matrix D is
00109 *>               exactly singular, and division by zero will occur if it
00110 *>               is used to solve a system of equations.
00111 *> \endverbatim
00112 *
00113 *  Authors:
00114 *  ========
00115 *
00116 *> \author Univ. of Tennessee 
00117 *> \author Univ. of California Berkeley 
00118 *> \author Univ. of Colorado Denver 
00119 *> \author NAG Ltd. 
00120 *
00121 *> \date November 2011
00122 *
00123 *> \ingroup complex16HEcomputational
00124 *
00125 *> \par Further Details:
00126 *  =====================
00127 *>
00128 *> \verbatim
00129 *>
00130 *>  If UPLO = 'U', then A = U*D*U**H, where
00131 *>     U = P(n)*U(n)* ... *P(k)U(k)* ...,
00132 *>  i.e., U is a product of terms P(k)*U(k), where k decreases from n to
00133 *>  1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
00134 *>  and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
00135 *>  defined by IPIV(k), and U(k) is a unit upper triangular matrix, such
00136 *>  that if the diagonal block D(k) is of order s (s = 1 or 2), then
00137 *>
00138 *>             (   I    v    0   )   k-s
00139 *>     U(k) =  (   0    I    0   )   s
00140 *>             (   0    0    I   )   n-k
00141 *>                k-s   s   n-k
00142 *>
00143 *>  If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k).
00144 *>  If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k),
00145 *>  and A(k,k), and v overwrites A(1:k-2,k-1:k).
00146 *>
00147 *>  If UPLO = 'L', then A = L*D*L**H, where
00148 *>     L = P(1)*L(1)* ... *P(k)*L(k)* ...,
00149 *>  i.e., L is a product of terms P(k)*L(k), where k increases from 1 to
00150 *>  n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
00151 *>  and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
00152 *>  defined by IPIV(k), and L(k) is a unit lower triangular matrix, such
00153 *>  that if the diagonal block D(k) is of order s (s = 1 or 2), then
00154 *>
00155 *>             (   I    0     0   )  k-1
00156 *>     L(k) =  (   0    I     0   )  s
00157 *>             (   0    v     I   )  n-k-s+1
00158 *>                k-1   s  n-k-s+1
00159 *>
00160 *>  If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k).
00161 *>  If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k),
00162 *>  and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1).
00163 *> \endverbatim
00164 *
00165 *> \par Contributors:
00166 *  ==================
00167 *>
00168 *> \verbatim
00169 *>  09-29-06 - patch from
00170 *>    Bobby Cheng, MathWorks
00171 *>
00172 *>    Replace l.210 and l.393
00173 *>         IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
00174 *>    by
00175 *>         IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. DISNAN(ABSAKK) ) THEN
00176 *>
00177 *>  01-01-96 - Based on modifications by
00178 *>    J. Lewis, Boeing Computer Services Company
00179 *>    A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA
00180 *> \endverbatim
00181 *
00182 *  =====================================================================
00183       SUBROUTINE ZHETF2( UPLO, N, A, LDA, IPIV, INFO )
00184 *
00185 *  -- LAPACK computational routine (version 3.4.0) --
00186 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00187 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00188 *     November 2011
00189 *
00190 *     .. Scalar Arguments ..
00191       CHARACTER          UPLO
00192       INTEGER            INFO, LDA, N
00193 *     ..
00194 *     .. Array Arguments ..
00195       INTEGER            IPIV( * )
00196       COMPLEX*16         A( LDA, * )
00197 *     ..
00198 *
00199 *  =====================================================================
00200 *
00201 *     .. Parameters ..
00202       DOUBLE PRECISION   ZERO, ONE
00203       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
00204       DOUBLE PRECISION   EIGHT, SEVTEN
00205       PARAMETER          ( EIGHT = 8.0D+0, SEVTEN = 17.0D+0 )
00206 *     ..
00207 *     .. Local Scalars ..
00208       LOGICAL            UPPER
00209       INTEGER            I, IMAX, J, JMAX, K, KK, KP, KSTEP
00210       DOUBLE PRECISION   ABSAKK, ALPHA, COLMAX, D, D11, D22, R1, ROWMAX,
00211      $                   TT
00212       COMPLEX*16         D12, D21, T, WK, WKM1, WKP1, ZDUM
00213 *     ..
00214 *     .. External Functions ..
00215       LOGICAL            LSAME, DISNAN
00216       INTEGER            IZAMAX
00217       DOUBLE PRECISION   DLAPY2
00218       EXTERNAL           LSAME, IZAMAX, DLAPY2, DISNAN
00219 *     ..
00220 *     .. External Subroutines ..
00221       EXTERNAL           XERBLA, ZDSCAL, ZHER, ZSWAP
00222 *     ..
00223 *     .. Intrinsic Functions ..
00224       INTRINSIC          ABS, DBLE, DCMPLX, DCONJG, DIMAG, MAX, SQRT
00225 *     ..
00226 *     .. Statement Functions ..
00227       DOUBLE PRECISION   CABS1
00228 *     ..
00229 *     .. Statement Function definitions ..
00230       CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
00231 *     ..
00232 *     .. Executable Statements ..
00233 *
00234 *     Test the input parameters.
00235 *
00236       INFO = 0
00237       UPPER = LSAME( UPLO, 'U' )
00238       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
00239          INFO = -1
00240       ELSE IF( N.LT.0 ) THEN
00241          INFO = -2
00242       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00243          INFO = -4
00244       END IF
00245       IF( INFO.NE.0 ) THEN
00246          CALL XERBLA( 'ZHETF2', -INFO )
00247          RETURN
00248       END IF
00249 *
00250 *     Initialize ALPHA for use in choosing pivot block size.
00251 *
00252       ALPHA = ( ONE+SQRT( SEVTEN ) ) / EIGHT
00253 *
00254       IF( UPPER ) THEN
00255 *
00256 *        Factorize A as U*D*U**H using the upper triangle of A
00257 *
00258 *        K is the main loop index, decreasing from N to 1 in steps of
00259 *        1 or 2
00260 *
00261          K = N
00262    10    CONTINUE
00263 *
00264 *        If K < 1, exit from loop
00265 *
00266          IF( K.LT.1 )
00267      $      GO TO 90
00268          KSTEP = 1
00269 *
00270 *        Determine rows and columns to be interchanged and whether
00271 *        a 1-by-1 or 2-by-2 pivot block will be used
00272 *
00273          ABSAKK = ABS( DBLE( A( K, K ) ) )
00274 *
00275 *        IMAX is the row-index of the largest off-diagonal element in
00276 *        column K, and COLMAX is its absolute value
00277 *
00278          IF( K.GT.1 ) THEN
00279             IMAX = IZAMAX( K-1, A( 1, K ), 1 )
00280             COLMAX = CABS1( A( IMAX, K ) )
00281          ELSE
00282             COLMAX = ZERO
00283          END IF
00284 *
00285          IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. DISNAN(ABSAKK) ) THEN
00286 *
00287 *           Column K is zero or contains a NaN: set INFO and continue
00288 *
00289             IF( INFO.EQ.0 )
00290      $         INFO = K
00291             KP = K
00292             A( K, K ) = DBLE( A( K, K ) )
00293          ELSE
00294             IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
00295 *
00296 *              no interchange, use 1-by-1 pivot block
00297 *
00298                KP = K
00299             ELSE
00300 *
00301 *              JMAX is the column-index of the largest off-diagonal
00302 *              element in row IMAX, and ROWMAX is its absolute value
00303 *
00304                JMAX = IMAX + IZAMAX( K-IMAX, A( IMAX, IMAX+1 ), LDA )
00305                ROWMAX = CABS1( A( IMAX, JMAX ) )
00306                IF( IMAX.GT.1 ) THEN
00307                   JMAX = IZAMAX( IMAX-1, A( 1, IMAX ), 1 )
00308                   ROWMAX = MAX( ROWMAX, CABS1( A( JMAX, IMAX ) ) )
00309                END IF
00310 *
00311                IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
00312 *
00313 *                 no interchange, use 1-by-1 pivot block
00314 *
00315                   KP = K
00316                ELSE IF( ABS( DBLE( A( IMAX, IMAX ) ) ).GE.ALPHA*ROWMAX )
00317      $                   THEN
00318 *
00319 *                 interchange rows and columns K and IMAX, use 1-by-1
00320 *                 pivot block
00321 *
00322                   KP = IMAX
00323                ELSE
00324 *
00325 *                 interchange rows and columns K-1 and IMAX, use 2-by-2
00326 *                 pivot block
00327 *
00328                   KP = IMAX
00329                   KSTEP = 2
00330                END IF
00331             END IF
00332 *
00333             KK = K - KSTEP + 1
00334             IF( KP.NE.KK ) THEN
00335 *
00336 *              Interchange rows and columns KK and KP in the leading
00337 *              submatrix A(1:k,1:k)
00338 *
00339                CALL ZSWAP( KP-1, A( 1, KK ), 1, A( 1, KP ), 1 )
00340                DO 20 J = KP + 1, KK - 1
00341                   T = DCONJG( A( J, KK ) )
00342                   A( J, KK ) = DCONJG( A( KP, J ) )
00343                   A( KP, J ) = T
00344    20          CONTINUE
00345                A( KP, KK ) = DCONJG( A( KP, KK ) )
00346                R1 = DBLE( A( KK, KK ) )
00347                A( KK, KK ) = DBLE( A( KP, KP ) )
00348                A( KP, KP ) = R1
00349                IF( KSTEP.EQ.2 ) THEN
00350                   A( K, K ) = DBLE( A( K, K ) )
00351                   T = A( K-1, K )
00352                   A( K-1, K ) = A( KP, K )
00353                   A( KP, K ) = T
00354                END IF
00355             ELSE
00356                A( K, K ) = DBLE( A( K, K ) )
00357                IF( KSTEP.EQ.2 )
00358      $            A( K-1, K-1 ) = DBLE( A( K-1, K-1 ) )
00359             END IF
00360 *
00361 *           Update the leading submatrix
00362 *
00363             IF( KSTEP.EQ.1 ) THEN
00364 *
00365 *              1-by-1 pivot block D(k): column k now holds
00366 *
00367 *              W(k) = U(k)*D(k)
00368 *
00369 *              where U(k) is the k-th column of U
00370 *
00371 *              Perform a rank-1 update of A(1:k-1,1:k-1) as
00372 *
00373 *              A := A - U(k)*D(k)*U(k)**H = A - W(k)*1/D(k)*W(k)**H
00374 *
00375                R1 = ONE / DBLE( A( K, K ) )
00376                CALL ZHER( UPLO, K-1, -R1, A( 1, K ), 1, A, LDA )
00377 *
00378 *              Store U(k) in column k
00379 *
00380                CALL ZDSCAL( K-1, R1, A( 1, K ), 1 )
00381             ELSE
00382 *
00383 *              2-by-2 pivot block D(k): columns k and k-1 now hold
00384 *
00385 *              ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
00386 *
00387 *              where U(k) and U(k-1) are the k-th and (k-1)-th columns
00388 *              of U
00389 *
00390 *              Perform a rank-2 update of A(1:k-2,1:k-2) as
00391 *
00392 *              A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )**H
00393 *                 = A - ( W(k-1) W(k) )*inv(D(k))*( W(k-1) W(k) )**H
00394 *
00395                IF( K.GT.2 ) THEN
00396 *
00397                   D = DLAPY2( DBLE( A( K-1, K ) ),
00398      $                DIMAG( A( K-1, K ) ) )
00399                   D22 = DBLE( A( K-1, K-1 ) ) / D
00400                   D11 = DBLE( A( K, K ) ) / D
00401                   TT = ONE / ( D11*D22-ONE )
00402                   D12 = A( K-1, K ) / D
00403                   D = TT / D
00404 *
00405                   DO 40 J = K - 2, 1, -1
00406                      WKM1 = D*( D11*A( J, K-1 )-DCONJG( D12 )*
00407      $                      A( J, K ) )
00408                      WK = D*( D22*A( J, K )-D12*A( J, K-1 ) )
00409                      DO 30 I = J, 1, -1
00410                         A( I, J ) = A( I, J ) - A( I, K )*DCONJG( WK ) -
00411      $                              A( I, K-1 )*DCONJG( WKM1 )
00412    30                CONTINUE
00413                      A( J, K ) = WK
00414                      A( J, K-1 ) = WKM1
00415                      A( J, J ) = DCMPLX( DBLE( A( J, J ) ), 0.0D+0 )
00416    40             CONTINUE
00417 *
00418                END IF
00419 *
00420             END IF
00421          END IF
00422 *
00423 *        Store details of the interchanges in IPIV
00424 *
00425          IF( KSTEP.EQ.1 ) THEN
00426             IPIV( K ) = KP
00427          ELSE
00428             IPIV( K ) = -KP
00429             IPIV( K-1 ) = -KP
00430          END IF
00431 *
00432 *        Decrease K and return to the start of the main loop
00433 *
00434          K = K - KSTEP
00435          GO TO 10
00436 *
00437       ELSE
00438 *
00439 *        Factorize A as L*D*L**H using the lower triangle of A
00440 *
00441 *        K is the main loop index, increasing from 1 to N in steps of
00442 *        1 or 2
00443 *
00444          K = 1
00445    50    CONTINUE
00446 *
00447 *        If K > N, exit from loop
00448 *
00449          IF( K.GT.N )
00450      $      GO TO 90
00451          KSTEP = 1
00452 *
00453 *        Determine rows and columns to be interchanged and whether
00454 *        a 1-by-1 or 2-by-2 pivot block will be used
00455 *
00456          ABSAKK = ABS( DBLE( A( K, K ) ) )
00457 *
00458 *        IMAX is the row-index of the largest off-diagonal element in
00459 *        column K, and COLMAX is its absolute value
00460 *
00461          IF( K.LT.N ) THEN
00462             IMAX = K + IZAMAX( N-K, A( K+1, K ), 1 )
00463             COLMAX = CABS1( A( IMAX, K ) )
00464          ELSE
00465             COLMAX = ZERO
00466          END IF
00467 *
00468          IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. DISNAN(ABSAKK) ) THEN
00469 *
00470 *           Column K is zero or contains a NaN: set INFO and continue
00471 *
00472             IF( INFO.EQ.0 )
00473      $         INFO = K
00474             KP = K
00475             A( K, K ) = DBLE( A( K, K ) )
00476          ELSE
00477             IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
00478 *
00479 *              no interchange, use 1-by-1 pivot block
00480 *
00481                KP = K
00482             ELSE
00483 *
00484 *              JMAX is the column-index of the largest off-diagonal
00485 *              element in row IMAX, and ROWMAX is its absolute value
00486 *
00487                JMAX = K - 1 + IZAMAX( IMAX-K, A( IMAX, K ), LDA )
00488                ROWMAX = CABS1( A( IMAX, JMAX ) )
00489                IF( IMAX.LT.N ) THEN
00490                   JMAX = IMAX + IZAMAX( N-IMAX, A( IMAX+1, IMAX ), 1 )
00491                   ROWMAX = MAX( ROWMAX, CABS1( A( JMAX, IMAX ) ) )
00492                END IF
00493 *
00494                IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
00495 *
00496 *                 no interchange, use 1-by-1 pivot block
00497 *
00498                   KP = K
00499                ELSE IF( ABS( DBLE( A( IMAX, IMAX ) ) ).GE.ALPHA*ROWMAX )
00500      $                   THEN
00501 *
00502 *                 interchange rows and columns K and IMAX, use 1-by-1
00503 *                 pivot block
00504 *
00505                   KP = IMAX
00506                ELSE
00507 *
00508 *                 interchange rows and columns K+1 and IMAX, use 2-by-2
00509 *                 pivot block
00510 *
00511                   KP = IMAX
00512                   KSTEP = 2
00513                END IF
00514             END IF
00515 *
00516             KK = K + KSTEP - 1
00517             IF( KP.NE.KK ) THEN
00518 *
00519 *              Interchange rows and columns KK and KP in the trailing
00520 *              submatrix A(k:n,k:n)
00521 *
00522                IF( KP.LT.N )
00523      $            CALL ZSWAP( N-KP, A( KP+1, KK ), 1, A( KP+1, KP ), 1 )
00524                DO 60 J = KK + 1, KP - 1
00525                   T = DCONJG( A( J, KK ) )
00526                   A( J, KK ) = DCONJG( A( KP, J ) )
00527                   A( KP, J ) = T
00528    60          CONTINUE
00529                A( KP, KK ) = DCONJG( A( KP, KK ) )
00530                R1 = DBLE( A( KK, KK ) )
00531                A( KK, KK ) = DBLE( A( KP, KP ) )
00532                A( KP, KP ) = R1
00533                IF( KSTEP.EQ.2 ) THEN
00534                   A( K, K ) = DBLE( A( K, K ) )
00535                   T = A( K+1, K )
00536                   A( K+1, K ) = A( KP, K )
00537                   A( KP, K ) = T
00538                END IF
00539             ELSE
00540                A( K, K ) = DBLE( A( K, K ) )
00541                IF( KSTEP.EQ.2 )
00542      $            A( K+1, K+1 ) = DBLE( A( K+1, K+1 ) )
00543             END IF
00544 *
00545 *           Update the trailing submatrix
00546 *
00547             IF( KSTEP.EQ.1 ) THEN
00548 *
00549 *              1-by-1 pivot block D(k): column k now holds
00550 *
00551 *              W(k) = L(k)*D(k)
00552 *
00553 *              where L(k) is the k-th column of L
00554 *
00555                IF( K.LT.N ) THEN
00556 *
00557 *                 Perform a rank-1 update of A(k+1:n,k+1:n) as
00558 *
00559 *                 A := A - L(k)*D(k)*L(k)**H = A - W(k)*(1/D(k))*W(k)**H
00560 *
00561                   R1 = ONE / DBLE( A( K, K ) )
00562                   CALL ZHER( UPLO, N-K, -R1, A( K+1, K ), 1,
00563      $                       A( K+1, K+1 ), LDA )
00564 *
00565 *                 Store L(k) in column K
00566 *
00567                   CALL ZDSCAL( N-K, R1, A( K+1, K ), 1 )
00568                END IF
00569             ELSE
00570 *
00571 *              2-by-2 pivot block D(k)
00572 *
00573                IF( K.LT.N-1 ) THEN
00574 *
00575 *                 Perform a rank-2 update of A(k+2:n,k+2:n) as
00576 *
00577 *                 A := A - ( L(k) L(k+1) )*D(k)*( L(k) L(k+1) )**H
00578 *                    = A - ( W(k) W(k+1) )*inv(D(k))*( W(k) W(k+1) )**H
00579 *
00580 *                 where L(k) and L(k+1) are the k-th and (k+1)-th
00581 *                 columns of L
00582 *
00583                   D = DLAPY2( DBLE( A( K+1, K ) ),
00584      $                DIMAG( A( K+1, K ) ) )
00585                   D11 = DBLE( A( K+1, K+1 ) ) / D
00586                   D22 = DBLE( A( K, K ) ) / D
00587                   TT = ONE / ( D11*D22-ONE )
00588                   D21 = A( K+1, K ) / D
00589                   D = TT / D
00590 *
00591                   DO 80 J = K + 2, N
00592                      WK = D*( D11*A( J, K )-D21*A( J, K+1 ) )
00593                      WKP1 = D*( D22*A( J, K+1 )-DCONJG( D21 )*
00594      $                      A( J, K ) )
00595                      DO 70 I = J, N
00596                         A( I, J ) = A( I, J ) - A( I, K )*DCONJG( WK ) -
00597      $                              A( I, K+1 )*DCONJG( WKP1 )
00598    70                CONTINUE
00599                      A( J, K ) = WK
00600                      A( J, K+1 ) = WKP1
00601                      A( J, J ) = DCMPLX( DBLE( A( J, J ) ), 0.0D+0 )
00602    80             CONTINUE
00603                END IF
00604             END IF
00605          END IF
00606 *
00607 *        Store details of the interchanges in IPIV
00608 *
00609          IF( KSTEP.EQ.1 ) THEN
00610             IPIV( K ) = KP
00611          ELSE
00612             IPIV( K ) = -KP
00613             IPIV( K+1 ) = -KP
00614          END IF
00615 *
00616 *        Increase K and return to the start of the main loop
00617 *
00618          K = K + KSTEP
00619          GO TO 50
00620 *
00621       END IF
00622 *
00623    90 CONTINUE
00624       RETURN
00625 *
00626 *     End of ZHETF2
00627 *
00628       END
 All Files Functions