LAPACK: SRC/zlahef.f Source File

Go to the documentation of this file.
00001       SUBROUTINE ZLAHEF( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, INFO )
00002 *
00003 *  -- LAPACK routine (version 3.2) --
00004 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00005 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00006 *     November 2006
00007 *
00008 *     .. Scalar Arguments ..
00009       CHARACTER          UPLO
00010       INTEGER            INFO, KB, LDA, LDW, N, NB
00011 *     ..
00012 *     .. Array Arguments ..
00013       INTEGER            IPIV( * )
00014       COMPLEX*16         A( LDA, * ), W( LDW, * )
00015 *     ..
00016 *
00017 *  Purpose
00018 *  =======
00019 *
00020 *  ZLAHEF computes a partial factorization of a complex Hermitian
00021 *  matrix A using the Bunch-Kaufman diagonal pivoting method. The
00022 *  partial factorization has the form:
00023 *
00024 *  A  =  ( I  U12 ) ( A11  0  ) (  I    0   )  if UPLO = 'U', or:
00025 *        ( 0  U22 ) (  0   D  ) ( U12' U22' )
00026 *
00027 *  A  =  ( L11  0 ) (  D   0  ) ( L11' L21' )  if UPLO = 'L'
00028 *        ( L21  I ) (  0  A22 ) (  0    I   )
00029 *
00030 *  where the order of D is at most NB. The actual order is returned in
00031 *  the argument KB, and is either NB or NB-1, or N if N <= NB.
00032 *  Note that U' denotes the conjugate transpose of U.
00033 *
00034 *  ZLAHEF is an auxiliary routine called by ZHETRF. It uses blocked code
00035 *  (calling Level 3 BLAS) to update the submatrix A11 (if UPLO = 'U') or
00036 *  A22 (if UPLO = 'L').
00037 *
00038 *  Arguments
00039 *  =========
00040 *
00041 *  UPLO    (input) CHARACTER*1
00042 *          Specifies whether the upper or lower triangular part of the
00043 *          Hermitian matrix A is stored:
00044 *          = 'U':  Upper triangular
00045 *          = 'L':  Lower triangular
00046 *
00047 *  N       (input) INTEGER
00048 *          The order of the matrix A.  N >= 0.
00049 *
00050 *  NB      (input) INTEGER
00051 *          The maximum number of columns of the matrix A that should be
00052 *          factored.  NB should be at least 2 to allow for 2-by-2 pivot
00053 *          blocks.
00054 *
00055 *  KB      (output) INTEGER
00056 *          The number of columns of A that were actually factored.
00057 *          KB is either NB-1 or NB, or N if N <= NB.
00058 *
00059 *  A       (input/output) COMPLEX*16 array, dimension (LDA,N)
00060 *          On entry, the Hermitian matrix A.  If UPLO = 'U', the leading
00061 *          n-by-n upper triangular part of A contains the upper
00062 *          triangular part of the matrix A, and the strictly lower
00063 *          triangular part of A is not referenced.  If UPLO = 'L', the
00064 *          leading n-by-n lower triangular part of A contains the lower
00065 *          triangular part of the matrix A, and the strictly upper
00066 *          triangular part of A is not referenced.
00067 *          On exit, A contains details of the partial factorization.
00068 *
00069 *  LDA     (input) INTEGER
00070 *          The leading dimension of the array A.  LDA >= max(1,N).
00071 *
00072 *  IPIV    (output) INTEGER array, dimension (N)
00073 *          Details of the interchanges and the block structure of D.
00074 *          If UPLO = 'U', only the last KB elements of IPIV are set;
00075 *          if UPLO = 'L', only the first KB elements are set.
00076 *
00077 *          If IPIV(k) > 0, then rows and columns k and IPIV(k) were
00078 *          interchanged and D(k,k) is a 1-by-1 diagonal block.
00079 *          If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and
00080 *          columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
00081 *          is a 2-by-2 diagonal block.  If UPLO = 'L' and IPIV(k) =
00082 *          IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were
00083 *          interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
00084 *
00085 *  W       (workspace) COMPLEX*16 array, dimension (LDW,NB)
00086 *
00087 *  LDW     (input) INTEGER
00088 *          The leading dimension of the array W.  LDW >= max(1,N).
00089 *
00090 *  INFO    (output) INTEGER
00091 *          = 0: successful exit
00092 *          > 0: if INFO = k, D(k,k) is exactly zero.  The factorization
00093 *               has been completed, but the block diagonal matrix D is
00094 *               exactly singular.
00095 *
00096 *  =====================================================================
00097 *
00098 *     .. Parameters ..
00099       DOUBLE PRECISION   ZERO, ONE
00100       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
00101       COMPLEX*16         CONE
00102       PARAMETER          ( CONE = ( 1.0D+0, 0.0D+0 ) )
00103       DOUBLE PRECISION   EIGHT, SEVTEN
00104       PARAMETER          ( EIGHT = 8.0D+0, SEVTEN = 17.0D+0 )
00105 *     ..
00106 *     .. Local Scalars ..
00107       INTEGER            IMAX, J, JB, JJ, JMAX, JP, K, KK, KKW, KP,
00108      $                   KSTEP, KW
00109       DOUBLE PRECISION   ABSAKK, ALPHA, COLMAX, R1, ROWMAX, T
00110       COMPLEX*16         D11, D21, D22, Z
00111 *     ..
00112 *     .. External Functions ..
00113       LOGICAL            LSAME
00114       INTEGER            IZAMAX
00115       EXTERNAL           LSAME, IZAMAX
00116 *     ..
00117 *     .. External Subroutines ..
00118       EXTERNAL           ZCOPY, ZDSCAL, ZGEMM, ZGEMV, ZLACGV, ZSWAP
00119 *     ..
00120 *     .. Intrinsic Functions ..
00121       INTRINSIC          ABS, DBLE, DCONJG, DIMAG, MAX, MIN, SQRT
00122 *     ..
00123 *     .. Statement Functions ..
00124       DOUBLE PRECISION   CABS1
00125 *     ..
00126 *     .. Statement Function definitions ..
00127       CABS1( Z ) = ABS( DBLE( Z ) ) + ABS( DIMAG( Z ) )
00128 *     ..
00129 *     .. Executable Statements ..
00130 *
00131       INFO = 0
00132 *
00133 *     Initialize ALPHA for use in choosing pivot block size.
00134 *
00135       ALPHA = ( ONE+SQRT( SEVTEN ) ) / EIGHT
00136 *
00137       IF( LSAME( UPLO, 'U' ) ) THEN
00138 *
00139 *        Factorize the trailing columns of A using the upper triangle
00140 *        of A and working backwards, and compute the matrix W = U12*D
00141 *        for use in updating A11 (note that conjg(W) is actually stored)
00142 *
00143 *        K is the main loop index, decreasing from N in steps of 1 or 2
00144 *
00145 *        KW is the column of W which corresponds to column K of A
00146 *
00147          K = N
00148    10    CONTINUE
00149          KW = NB + K - N
00150 *
00151 *        Exit from loop
00152 *
00153          IF( ( K.LE.N-NB+1 .AND. NB.LT.N ) .OR. K.LT.1 )
00154      $      GO TO 30
00155 *
00156 *        Copy column K of A to column KW of W and update it
00157 *
00158          CALL ZCOPY( K-1, A( 1, K ), 1, W( 1, KW ), 1 )
00159          W( K, KW ) = DBLE( A( K, K ) )
00160          IF( K.LT.N ) THEN
00161             CALL ZGEMV( 'No transpose', K, N-K, -CONE, A( 1, K+1 ), LDA,
00162      $                  W( K, KW+1 ), LDW, CONE, W( 1, KW ), 1 )
00163             W( K, KW ) = DBLE( W( K, KW ) )
00164          END IF
00165 *
00166          KSTEP = 1
00167 *
00168 *        Determine rows and columns to be interchanged and whether
00169 *        a 1-by-1 or 2-by-2 pivot block will be used
00170 *
00171          ABSAKK = ABS( DBLE( W( K, KW ) ) )
00172 *
00173 *        IMAX is the row-index of the largest off-diagonal element in
00174 *        column K, and COLMAX is its absolute value
00175 *
00176          IF( K.GT.1 ) THEN
00177             IMAX = IZAMAX( K-1, W( 1, KW ), 1 )
00178             COLMAX = CABS1( W( IMAX, KW ) )
00179          ELSE
00180             COLMAX = ZERO
00181          END IF
00182 *
00183          IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
00184 *
00185 *           Column K is zero: set INFO and continue
00186 *
00187             IF( INFO.EQ.0 )
00188      $         INFO = K
00189             KP = K
00190             A( K, K ) = DBLE( A( K, K ) )
00191          ELSE
00192             IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
00193 *
00194 *              no interchange, use 1-by-1 pivot block
00195 *
00196                KP = K
00197             ELSE
00198 *
00199 *              Copy column IMAX to column KW-1 of W and update it
00200 *
00201                CALL ZCOPY( IMAX-1, A( 1, IMAX ), 1, W( 1, KW-1 ), 1 )
00202                W( IMAX, KW-1 ) = DBLE( A( IMAX, IMAX ) )
00203                CALL ZCOPY( K-IMAX, A( IMAX, IMAX+1 ), LDA,
00204      $                     W( IMAX+1, KW-1 ), 1 )
00205                CALL ZLACGV( K-IMAX, W( IMAX+1, KW-1 ), 1 )
00206                IF( K.LT.N ) THEN
00207                   CALL ZGEMV( 'No transpose', K, N-K, -CONE,
00208      $                        A( 1, K+1 ), LDA, W( IMAX, KW+1 ), LDW,
00209      $                        CONE, W( 1, KW-1 ), 1 )
00210                   W( IMAX, KW-1 ) = DBLE( W( IMAX, KW-1 ) )
00211                END IF
00212 *
00213 *              JMAX is the column-index of the largest off-diagonal
00214 *              element in row IMAX, and ROWMAX is its absolute value
00215 *
00216                JMAX = IMAX + IZAMAX( K-IMAX, W( IMAX+1, KW-1 ), 1 )
00217                ROWMAX = CABS1( W( JMAX, KW-1 ) )
00218                IF( IMAX.GT.1 ) THEN
00219                   JMAX = IZAMAX( IMAX-1, W( 1, KW-1 ), 1 )
00220                   ROWMAX = MAX( ROWMAX, CABS1( W( JMAX, KW-1 ) ) )
00221                END IF
00222 *
00223                IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
00224 *
00225 *                 no interchange, use 1-by-1 pivot block
00226 *
00227                   KP = K
00228                ELSE IF( ABS( DBLE( W( IMAX, KW-1 ) ) ).GE.ALPHA*ROWMAX )
00229      $                   THEN
00230 *
00231 *                 interchange rows and columns K and IMAX, use 1-by-1
00232 *                 pivot block
00233 *
00234                   KP = IMAX
00235 *
00236 *                 copy column KW-1 of W to column KW
00237 *
00238                   CALL ZCOPY( K, W( 1, KW-1 ), 1, W( 1, KW ), 1 )
00239                ELSE
00240 *
00241 *                 interchange rows and columns K-1 and IMAX, use 2-by-2
00242 *                 pivot block
00243 *
00244                   KP = IMAX
00245                   KSTEP = 2
00246                END IF
00247             END IF
00248 *
00249             KK = K - KSTEP + 1
00250             KKW = NB + KK - N
00251 *
00252 *           Updated column KP is already stored in column KKW of W
00253 *
00254             IF( KP.NE.KK ) THEN
00255 *
00256 *              Copy non-updated column KK to column KP
00257 *
00258                A( KP, KP ) = DBLE( A( KK, KK ) )
00259                CALL ZCOPY( KK-1-KP, A( KP+1, KK ), 1, A( KP, KP+1 ),
00260      $                     LDA )
00261                CALL ZLACGV( KK-1-KP, A( KP, KP+1 ), LDA )
00262                CALL ZCOPY( KP-1, A( 1, KK ), 1, A( 1, KP ), 1 )
00263 *
00264 *              Interchange rows KK and KP in last KK columns of A and W
00265 *
00266                IF( KK.LT.N )
00267      $            CALL ZSWAP( N-KK, A( KK, KK+1 ), LDA, A( KP, KK+1 ),
00268      $                        LDA )
00269                CALL ZSWAP( N-KK+1, W( KK, KKW ), LDW, W( KP, KKW ),
00270      $                     LDW )
00271             END IF
00272 *
00273             IF( KSTEP.EQ.1 ) THEN
00274 *
00275 *              1-by-1 pivot block D(k): column KW of W now holds
00276 *
00277 *              W(k) = U(k)*D(k)
00278 *
00279 *              where U(k) is the k-th column of U
00280 *
00281 *              Store U(k) in column k of A
00282 *
00283                CALL ZCOPY( K, W( 1, KW ), 1, A( 1, K ), 1 )
00284                R1 = ONE / DBLE( A( K, K ) )
00285                CALL ZDSCAL( K-1, R1, A( 1, K ), 1 )
00286 *
00287 *              Conjugate W(k)
00288 *
00289                CALL ZLACGV( K-1, W( 1, KW ), 1 )
00290             ELSE
00291 *
00292 *              2-by-2 pivot block D(k): columns KW and KW-1 of W now
00293 *              hold
00294 *
00295 *              ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
00296 *
00297 *              where U(k) and U(k-1) are the k-th and (k-1)-th columns
00298 *              of U
00299 *
00300                IF( K.GT.2 ) THEN
00301 *
00302 *                 Store U(k) and U(k-1) in columns k and k-1 of A
00303 *
00304                   D21 = W( K-1, KW )
00305                   D11 = W( K, KW ) / DCONJG( D21 )
00306                   D22 = W( K-1, KW-1 ) / D21
00307                   T = ONE / ( DBLE( D11*D22 )-ONE )
00308                   D21 = T / D21
00309                   DO 20 J = 1, K - 2
00310                      A( J, K-1 ) = D21*( D11*W( J, KW-1 )-W( J, KW ) )
00311                      A( J, K ) = DCONJG( D21 )*
00312      $                           ( D22*W( J, KW )-W( J, KW-1 ) )
00313    20             CONTINUE
00314                END IF
00315 *
00316 *              Copy D(k) to A
00317 *
00318                A( K-1, K-1 ) = W( K-1, KW-1 )
00319                A( K-1, K ) = W( K-1, KW )
00320                A( K, K ) = W( K, KW )
00321 *
00322 *              Conjugate W(k) and W(k-1)
00323 *
00324                CALL ZLACGV( K-1, W( 1, KW ), 1 )
00325                CALL ZLACGV( K-2, W( 1, KW-1 ), 1 )
00326             END IF
00327          END IF
00328 *
00329 *        Store details of the interchanges in IPIV
00330 *
00331          IF( KSTEP.EQ.1 ) THEN
00332             IPIV( K ) = KP
00333          ELSE
00334             IPIV( K ) = -KP
00335             IPIV( K-1 ) = -KP
00336          END IF
00337 *
00338 *        Decrease K and return to the start of the main loop
00339 *
00340          K = K - KSTEP
00341          GO TO 10
00342 *
00343    30    CONTINUE
00344 *
00345 *        Update the upper triangle of A11 (= A(1:k,1:k)) as
00346 *
00347 *        A11 := A11 - U12*D*U12' = A11 - U12*W'
00348 *
00349 *        computing blocks of NB columns at a time (note that conjg(W) is
00350 *        actually stored)
00351 *
00352          DO 50 J = ( ( K-1 ) / NB )*NB + 1, 1, -NB
00353             JB = MIN( NB, K-J+1 )
00354 *
00355 *           Update the upper triangle of the diagonal block
00356 *
00357             DO 40 JJ = J, J + JB - 1
00358                A( JJ, JJ ) = DBLE( A( JJ, JJ ) )
00359                CALL ZGEMV( 'No transpose', JJ-J+1, N-K, -CONE,
00360      $                     A( J, K+1 ), LDA, W( JJ, KW+1 ), LDW, CONE,
00361      $                     A( J, JJ ), 1 )
00362                A( JJ, JJ ) = DBLE( A( JJ, JJ ) )
00363    40       CONTINUE
00364 *
00365 *           Update the rectangular superdiagonal block
00366 *
00367             CALL ZGEMM( 'No transpose', 'Transpose', J-1, JB, N-K,
00368      $                  -CONE, A( 1, K+1 ), LDA, W( J, KW+1 ), LDW,
00369      $                  CONE, A( 1, J ), LDA )
00370    50    CONTINUE
00371 *
00372 *        Put U12 in standard form by partially undoing the interchanges
00373 *        in columns k+1:n
00374 *
00375          J = K + 1
00376    60    CONTINUE
00377          JJ = J
00378          JP = IPIV( J )
00379          IF( JP.LT.0 ) THEN
00380             JP = -JP
00381             J = J + 1
00382          END IF
00383          J = J + 1
00384          IF( JP.NE.JJ .AND. J.LE.N )
00385      $      CALL ZSWAP( N-J+1, A( JP, J ), LDA, A( JJ, J ), LDA )
00386          IF( J.LE.N )
00387      $      GO TO 60
00388 *
00389 *        Set KB to the number of columns factorized
00390 *
00391          KB = N - K
00392 *
00393       ELSE
00394 *
00395 *        Factorize the leading columns of A using the lower triangle
00396 *        of A and working forwards, and compute the matrix W = L21*D
00397 *        for use in updating A22 (note that conjg(W) is actually stored)
00398 *
00399 *        K is the main loop index, increasing from 1 in steps of 1 or 2
00400 *
00401          K = 1
00402    70    CONTINUE
00403 *
00404 *        Exit from loop
00405 *
00406          IF( ( K.GE.NB .AND. NB.LT.N ) .OR. K.GT.N )
00407      $      GO TO 90
00408 *
00409 *        Copy column K of A to column K of W and update it
00410 *
00411          W( K, K ) = DBLE( A( K, K ) )
00412          IF( K.LT.N )
00413      $      CALL ZCOPY( N-K, A( K+1, K ), 1, W( K+1, K ), 1 )
00414          CALL ZGEMV( 'No transpose', N-K+1, K-1, -CONE, A( K, 1 ), LDA,
00415      $               W( K, 1 ), LDW, CONE, W( K, K ), 1 )
00416          W( K, K ) = DBLE( W( K, K ) )
00417 *
00418          KSTEP = 1
00419 *
00420 *        Determine rows and columns to be interchanged and whether
00421 *        a 1-by-1 or 2-by-2 pivot block will be used
00422 *
00423          ABSAKK = ABS( DBLE( W( K, K ) ) )
00424 *
00425 *        IMAX is the row-index of the largest off-diagonal element in
00426 *        column K, and COLMAX is its absolute value
00427 *
00428          IF( K.LT.N ) THEN
00429             IMAX = K + IZAMAX( N-K, W( K+1, K ), 1 )
00430             COLMAX = CABS1( W( IMAX, K ) )
00431          ELSE
00432             COLMAX = ZERO
00433          END IF
00434 *
00435          IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
00436 *
00437 *           Column K is zero: set INFO and continue
00438 *
00439             IF( INFO.EQ.0 )
00440      $         INFO = K
00441             KP = K
00442             A( K, K ) = DBLE( A( K, K ) )
00443          ELSE
00444             IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
00445 *
00446 *              no interchange, use 1-by-1 pivot block
00447 *
00448                KP = K
00449             ELSE
00450 *
00451 *              Copy column IMAX to column K+1 of W and update it
00452 *
00453                CALL ZCOPY( IMAX-K, A( IMAX, K ), LDA, W( K, K+1 ), 1 )
00454                CALL ZLACGV( IMAX-K, W( K, K+1 ), 1 )
00455                W( IMAX, K+1 ) = DBLE( A( IMAX, IMAX ) )
00456                IF( IMAX.LT.N )
00457      $            CALL ZCOPY( N-IMAX, A( IMAX+1, IMAX ), 1,
00458      $                        W( IMAX+1, K+1 ), 1 )
00459                CALL ZGEMV( 'No transpose', N-K+1, K-1, -CONE, A( K, 1 ),
00460      $                     LDA, W( IMAX, 1 ), LDW, CONE, W( K, K+1 ),
00461      $                     1 )
00462                W( IMAX, K+1 ) = DBLE( W( IMAX, K+1 ) )
00463 *
00464 *              JMAX is the column-index of the largest off-diagonal
00465 *              element in row IMAX, and ROWMAX is its absolute value
00466 *
00467                JMAX = K - 1 + IZAMAX( IMAX-K, W( K, K+1 ), 1 )
00468                ROWMAX = CABS1( W( JMAX, K+1 ) )
00469                IF( IMAX.LT.N ) THEN
00470                   JMAX = IMAX + IZAMAX( N-IMAX, W( IMAX+1, K+1 ), 1 )
00471                   ROWMAX = MAX( ROWMAX, CABS1( W( JMAX, K+1 ) ) )
00472                END IF
00473 *
00474                IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
00475 *
00476 *                 no interchange, use 1-by-1 pivot block
00477 *
00478                   KP = K
00479                ELSE IF( ABS( DBLE( W( IMAX, K+1 ) ) ).GE.ALPHA*ROWMAX )
00480      $                   THEN
00481 *
00482 *                 interchange rows and columns K and IMAX, use 1-by-1
00483 *                 pivot block
00484 *
00485                   KP = IMAX
00486 *
00487 *                 copy column K+1 of W to column K
00488 *
00489                   CALL ZCOPY( N-K+1, W( K, K+1 ), 1, W( K, K ), 1 )
00490                ELSE
00491 *
00492 *                 interchange rows and columns K+1 and IMAX, use 2-by-2
00493 *                 pivot block
00494 *
00495                   KP = IMAX
00496                   KSTEP = 2
00497                END IF
00498             END IF
00499 *
00500             KK = K + KSTEP - 1
00501 *
00502 *           Updated column KP is already stored in column KK of W
00503 *
00504             IF( KP.NE.KK ) THEN
00505 *
00506 *              Copy non-updated column KK to column KP
00507 *
00508                A( KP, KP ) = DBLE( A( KK, KK ) )
00509                CALL ZCOPY( KP-KK-1, A( KK+1, KK ), 1, A( KP, KK+1 ),
00510      $                     LDA )
00511                CALL ZLACGV( KP-KK-1, A( KP, KK+1 ), LDA )
00512                IF( KP.LT.N )
00513      $            CALL ZCOPY( N-KP, A( KP+1, KK ), 1, A( KP+1, KP ), 1 )
00514 *
00515 *              Interchange rows KK and KP in first KK columns of A and W
00516 *
00517                CALL ZSWAP( KK-1, A( KK, 1 ), LDA, A( KP, 1 ), LDA )
00518                CALL ZSWAP( KK, W( KK, 1 ), LDW, W( KP, 1 ), LDW )
00519             END IF
00520 *
00521             IF( KSTEP.EQ.1 ) THEN
00522 *
00523 *              1-by-1 pivot block D(k): column k of W now holds
00524 *
00525 *              W(k) = L(k)*D(k)
00526 *
00527 *              where L(k) is the k-th column of L
00528 *
00529 *              Store L(k) in column k of A
00530 *
00531                CALL ZCOPY( N-K+1, W( K, K ), 1, A( K, K ), 1 )
00532                IF( K.LT.N ) THEN
00533                   R1 = ONE / DBLE( A( K, K ) )
00534                   CALL ZDSCAL( N-K, R1, A( K+1, K ), 1 )
00535 *
00536 *                 Conjugate W(k)
00537 *
00538                   CALL ZLACGV( N-K, W( K+1, K ), 1 )
00539                END IF
00540             ELSE
00541 *
00542 *              2-by-2 pivot block D(k): columns k and k+1 of W now hold
00543 *
00544 *              ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k)
00545 *
00546 *              where L(k) and L(k+1) are the k-th and (k+1)-th columns
00547 *              of L
00548 *
00549                IF( K.LT.N-1 ) THEN
00550 *
00551 *                 Store L(k) and L(k+1) in columns k and k+1 of A
00552 *
00553                   D21 = W( K+1, K )
00554                   D11 = W( K+1, K+1 ) / D21
00555                   D22 = W( K, K ) / DCONJG( D21 )
00556                   T = ONE / ( DBLE( D11*D22 )-ONE )
00557                   D21 = T / D21
00558                   DO 80 J = K + 2, N
00559                      A( J, K ) = DCONJG( D21 )*
00560      $                           ( D11*W( J, K )-W( J, K+1 ) )
00561                      A( J, K+1 ) = D21*( D22*W( J, K+1 )-W( J, K ) )
00562    80             CONTINUE
00563                END IF
00564 *
00565 *              Copy D(k) to A
00566 *
00567                A( K, K ) = W( K, K )
00568                A( K+1, K ) = W( K+1, K )
00569                A( K+1, K+1 ) = W( K+1, K+1 )
00570 *
00571 *              Conjugate W(k) and W(k+1)
00572 *
00573                CALL ZLACGV( N-K, W( K+1, K ), 1 )
00574                CALL ZLACGV( N-K-1, W( K+2, K+1 ), 1 )
00575             END IF
00576          END IF
00577 *
00578 *        Store details of the interchanges in IPIV
00579 *
00580          IF( KSTEP.EQ.1 ) THEN
00581             IPIV( K ) = KP
00582          ELSE
00583             IPIV( K ) = -KP
00584             IPIV( K+1 ) = -KP
00585          END IF
00586 *
00587 *        Increase K and return to the start of the main loop
00588 *
00589          K = K + KSTEP
00590          GO TO 70
00591 *
00592    90    CONTINUE
00593 *
00594 *        Update the lower triangle of A22 (= A(k:n,k:n)) as
00595 *
00596 *        A22 := A22 - L21*D*L21' = A22 - L21*W'
00597 *
00598 *        computing blocks of NB columns at a time (note that conjg(W) is
00599 *        actually stored)
00600 *
00601          DO 110 J = K, N, NB
00602             JB = MIN( NB, N-J+1 )
00603 *
00604 *           Update the lower triangle of the diagonal block
00605 *
00606             DO 100 JJ = J, J + JB - 1
00607                A( JJ, JJ ) = DBLE( A( JJ, JJ ) )
00608                CALL ZGEMV( 'No transpose', J+JB-JJ, K-1, -CONE,
00609      $                     A( JJ, 1 ), LDA, W( JJ, 1 ), LDW, CONE,
00610      $                     A( JJ, JJ ), 1 )
00611                A( JJ, JJ ) = DBLE( A( JJ, JJ ) )
00612   100       CONTINUE
00613 *
00614 *           Update the rectangular subdiagonal block
00615 *
00616             IF( J+JB.LE.N )
00617      $         CALL ZGEMM( 'No transpose', 'Transpose', N-J-JB+1, JB,
00618      $                     K-1, -CONE, A( J+JB, 1 ), LDA, W( J, 1 ),
00619      $                     LDW, CONE, A( J+JB, J ), LDA )
00620   110    CONTINUE
00621 *
00622 *        Put L21 in standard form by partially undoing the interchanges
00623 *        in columns 1:k-1
00624 *
00625          J = K - 1
00626   120    CONTINUE
00627          JJ = J
00628          JP = IPIV( J )
00629          IF( JP.LT.0 ) THEN
00630             JP = -JP
00631             J = J - 1
00632          END IF
00633          J = J - 1
00634          IF( JP.NE.JJ .AND. J.GE.1 )
00635      $      CALL ZSWAP( J, A( JP, 1 ), LDA, A( JJ, 1 ), LDA )
00636          IF( J.GE.1 )
00637      $      GO TO 120
00638 *
00639 *        Set KB to the number of columns factorized
00640 *
00641          KB = K - 1
00642 *
00643       END IF
00644       RETURN
00645 *
00646 *     End of ZLAHEF
00647 *
00648       END