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