LAPACK: SRC/clasyf.f Source File

Go to the documentation of this file.
00001       SUBROUTINE CLASYF( 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            A( LDA, * ), W( LDW, * )
00015 *     ..
00016 *
00017 *  Purpose
00018 *  =======
00019 *
00020 *  CLASYF computes a partial factorization of a complex symmetric matrix
00021 *  A using the Bunch-Kaufman diagonal pivoting method. The partial
00022 *  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 transpose of U.
00033 *
00034 *  CLASYF is an auxiliary routine called by CSYTRF. 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 *          symmetric 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 array, dimension (LDA,N)
00060 *          On entry, the symmetric 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 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       REAL               ZERO, ONE
00100       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
00101       REAL               EIGHT, SEVTEN
00102       PARAMETER          ( EIGHT = 8.0E+0, SEVTEN = 17.0E+0 )
00103       COMPLEX            CONE
00104       PARAMETER          ( CONE = ( 1.0E+0, 0.0E+0 ) )
00105 *     ..
00106 *     .. Local Scalars ..
00107       INTEGER            IMAX, J, JB, JJ, JMAX, JP, K, KK, KKW, KP,
00108      $                   KSTEP, KW
00109       REAL               ABSAKK, ALPHA, COLMAX, ROWMAX
00110       COMPLEX            D11, D21, D22, R1, T, Z
00111 *     ..
00112 *     .. External Functions ..
00113       LOGICAL            LSAME
00114       INTEGER            ICAMAX
00115       EXTERNAL           LSAME, ICAMAX
00116 *     ..
00117 *     .. External Subroutines ..
00118       EXTERNAL           CCOPY, CGEMM, CGEMV, CSCAL, CSWAP
00119 *     ..
00120 *     .. Intrinsic Functions ..
00121       INTRINSIC          ABS, AIMAG, MAX, MIN, REAL, SQRT
00122 *     ..
00123 *     .. Statement Functions ..
00124       REAL               CABS1
00125 *     ..
00126 *     .. Statement Function definitions ..
00127       CABS1( Z ) = ABS( REAL( Z ) ) + ABS( AIMAG( 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
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 CCOPY( K, A( 1, K ), 1, W( 1, KW ), 1 )
00159          IF( K.LT.N )
00160      $      CALL CGEMV( 'No transpose', K, N-K, -CONE, A( 1, K+1 ), LDA,
00161      $                  W( K, KW+1 ), LDW, CONE, W( 1, KW ), 1 )
00162 *
00163          KSTEP = 1
00164 *
00165 *        Determine rows and columns to be interchanged and whether
00166 *        a 1-by-1 or 2-by-2 pivot block will be used
00167 *
00168          ABSAKK = CABS1( W( K, KW ) )
00169 *
00170 *        IMAX is the row-index of the largest off-diagonal element in
00171 *        column K, and COLMAX is its absolute value
00172 *
00173          IF( K.GT.1 ) THEN
00174             IMAX = ICAMAX( K-1, W( 1, KW ), 1 )
00175             COLMAX = CABS1( W( IMAX, KW ) )
00176          ELSE
00177             COLMAX = ZERO
00178          END IF
00179 *
00180          IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
00181 *
00182 *           Column K is zero: set INFO and continue
00183 *
00184             IF( INFO.EQ.0 )
00185      $         INFO = K
00186             KP = K
00187          ELSE
00188             IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
00189 *
00190 *              no interchange, use 1-by-1 pivot block
00191 *
00192                KP = K
00193             ELSE
00194 *
00195 *              Copy column IMAX to column KW-1 of W and update it
00196 *
00197                CALL CCOPY( IMAX, A( 1, IMAX ), 1, W( 1, KW-1 ), 1 )
00198                CALL CCOPY( K-IMAX, A( IMAX, IMAX+1 ), LDA,
00199      $                     W( IMAX+1, KW-1 ), 1 )
00200                IF( K.LT.N )
00201      $            CALL CGEMV( 'No transpose', K, N-K, -CONE,
00202      $                        A( 1, K+1 ), LDA, W( IMAX, KW+1 ), LDW,
00203      $                        CONE, W( 1, KW-1 ), 1 )
00204 *
00205 *              JMAX is the column-index of the largest off-diagonal
00206 *              element in row IMAX, and ROWMAX is its absolute value
00207 *
00208                JMAX = IMAX + ICAMAX( K-IMAX, W( IMAX+1, KW-1 ), 1 )
00209                ROWMAX = CABS1( W( JMAX, KW-1 ) )
00210                IF( IMAX.GT.1 ) THEN
00211                   JMAX = ICAMAX( IMAX-1, W( 1, KW-1 ), 1 )
00212                   ROWMAX = MAX( ROWMAX, CABS1( W( JMAX, KW-1 ) ) )
00213                END IF
00214 *
00215                IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
00216 *
00217 *                 no interchange, use 1-by-1 pivot block
00218 *
00219                   KP = K
00220                ELSE IF( CABS1( W( IMAX, KW-1 ) ).GE.ALPHA*ROWMAX ) THEN
00221 *
00222 *                 interchange rows and columns K and IMAX, use 1-by-1
00223 *                 pivot block
00224 *
00225                   KP = IMAX
00226 *
00227 *                 copy column KW-1 of W to column KW
00228 *
00229                   CALL CCOPY( K, W( 1, KW-1 ), 1, W( 1, KW ), 1 )
00230                ELSE
00231 *
00232 *                 interchange rows and columns K-1 and IMAX, use 2-by-2
00233 *                 pivot block
00234 *
00235                   KP = IMAX
00236                   KSTEP = 2
00237                END IF
00238             END IF
00239 *
00240             KK = K - KSTEP + 1
00241             KKW = NB + KK - N
00242 *
00243 *           Updated column KP is already stored in column KKW of W
00244 *
00245             IF( KP.NE.KK ) THEN
00246 *
00247 *              Copy non-updated column KK to column KP
00248 *
00249                A( KP, K ) = A( KK, K )
00250                CALL CCOPY( K-1-KP, A( KP+1, KK ), 1, A( KP, KP+1 ),
00251      $                     LDA )
00252                CALL CCOPY( KP, A( 1, KK ), 1, A( 1, KP ), 1 )
00253 *
00254 *              Interchange rows KK and KP in last KK columns of A and W
00255 *
00256                CALL CSWAP( N-KK+1, A( KK, KK ), LDA, A( KP, KK ), LDA )
00257                CALL CSWAP( N-KK+1, W( KK, KKW ), LDW, W( KP, KKW ),
00258      $                     LDW )
00259             END IF
00260 *
00261             IF( KSTEP.EQ.1 ) THEN
00262 *
00263 *              1-by-1 pivot block D(k): column KW of W now holds
00264 *
00265 *              W(k) = U(k)*D(k)
00266 *
00267 *              where U(k) is the k-th column of U
00268 *
00269 *              Store U(k) in column k of A
00270 *
00271                CALL CCOPY( K, W( 1, KW ), 1, A( 1, K ), 1 )
00272                R1 = CONE / A( K, K )
00273                CALL CSCAL( K-1, R1, A( 1, K ), 1 )
00274             ELSE
00275 *
00276 *              2-by-2 pivot block D(k): columns KW and KW-1 of W now
00277 *              hold
00278 *
00279 *              ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
00280 *
00281 *              where U(k) and U(k-1) are the k-th and (k-1)-th columns
00282 *              of U
00283 *
00284                IF( K.GT.2 ) THEN
00285 *
00286 *                 Store U(k) and U(k-1) in columns k and k-1 of A
00287 *
00288                   D21 = W( K-1, KW )
00289                   D11 = W( K, KW ) / D21
00290                   D22 = W( K-1, KW-1 ) / D21
00291                   T = CONE / ( D11*D22-CONE )
00292                   D21 = T / D21
00293                   DO 20 J = 1, K - 2
00294                      A( J, K-1 ) = D21*( D11*W( J, KW-1 )-W( J, KW ) )
00295                      A( J, K ) = D21*( D22*W( J, KW )-W( J, KW-1 ) )
00296    20             CONTINUE
00297                END IF
00298 *
00299 *              Copy D(k) to A
00300 *
00301                A( K-1, K-1 ) = W( K-1, KW-1 )
00302                A( K-1, K ) = W( K-1, KW )
00303                A( K, K ) = W( K, KW )
00304             END IF
00305          END IF
00306 *
00307 *        Store details of the interchanges in IPIV
00308 *
00309          IF( KSTEP.EQ.1 ) THEN
00310             IPIV( K ) = KP
00311          ELSE
00312             IPIV( K ) = -KP
00313             IPIV( K-1 ) = -KP
00314          END IF
00315 *
00316 *        Decrease K and return to the start of the main loop
00317 *
00318          K = K - KSTEP
00319          GO TO 10
00320 *
00321    30    CONTINUE
00322 *
00323 *        Update the upper triangle of A11 (= A(1:k,1:k)) as
00324 *
00325 *        A11 := A11 - U12*D*U12' = A11 - U12*W'
00326 *
00327 *        computing blocks of NB columns at a time
00328 *
00329          DO 50 J = ( ( K-1 ) / NB )*NB + 1, 1, -NB
00330             JB = MIN( NB, K-J+1 )
00331 *
00332 *           Update the upper triangle of the diagonal block
00333 *
00334             DO 40 JJ = J, J + JB - 1
00335                CALL CGEMV( 'No transpose', JJ-J+1, N-K, -CONE,
00336      $                     A( J, K+1 ), LDA, W( JJ, KW+1 ), LDW, CONE,
00337      $                     A( J, JJ ), 1 )
00338    40       CONTINUE
00339 *
00340 *           Update the rectangular superdiagonal block
00341 *
00342             CALL CGEMM( 'No transpose', 'Transpose', J-1, JB, N-K,
00343      $                  -CONE, A( 1, K+1 ), LDA, W( J, KW+1 ), LDW,
00344      $                  CONE, A( 1, J ), LDA )
00345    50    CONTINUE
00346 *
00347 *        Put U12 in standard form by partially undoing the interchanges
00348 *        in columns k+1:n
00349 *
00350          J = K + 1
00351    60    CONTINUE
00352          JJ = J
00353          JP = IPIV( J )
00354          IF( JP.LT.0 ) THEN
00355             JP = -JP
00356             J = J + 1
00357          END IF
00358          J = J + 1
00359          IF( JP.NE.JJ .AND. J.LE.N )
00360      $      CALL CSWAP( N-J+1, A( JP, J ), LDA, A( JJ, J ), LDA )
00361          IF( J.LE.N )
00362      $      GO TO 60
00363 *
00364 *        Set KB to the number of columns factorized
00365 *
00366          KB = N - K
00367 *
00368       ELSE
00369 *
00370 *        Factorize the leading columns of A using the lower triangle
00371 *        of A and working forwards, and compute the matrix W = L21*D
00372 *        for use in updating A22
00373 *
00374 *        K is the main loop index, increasing from 1 in steps of 1 or 2
00375 *
00376          K = 1
00377    70    CONTINUE
00378 *
00379 *        Exit from loop
00380 *
00381          IF( ( K.GE.NB .AND. NB.LT.N ) .OR. K.GT.N )
00382      $      GO TO 90
00383 *
00384 *        Copy column K of A to column K of W and update it
00385 *
00386          CALL CCOPY( N-K+1, A( K, K ), 1, W( K, K ), 1 )
00387          CALL CGEMV( 'No transpose', N-K+1, K-1, -CONE, A( K, 1 ), LDA,
00388      $               W( K, 1 ), LDW, CONE, W( K, K ), 1 )
00389 *
00390          KSTEP = 1
00391 *
00392 *        Determine rows and columns to be interchanged and whether
00393 *        a 1-by-1 or 2-by-2 pivot block will be used
00394 *
00395          ABSAKK = CABS1( W( K, K ) )
00396 *
00397 *        IMAX is the row-index of the largest off-diagonal element in
00398 *        column K, and COLMAX is its absolute value
00399 *
00400          IF( K.LT.N ) THEN
00401             IMAX = K + ICAMAX( N-K, W( K+1, K ), 1 )
00402             COLMAX = CABS1( W( IMAX, K ) )
00403          ELSE
00404             COLMAX = ZERO
00405          END IF
00406 *
00407          IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
00408 *
00409 *           Column K is zero: set INFO and continue
00410 *
00411             IF( INFO.EQ.0 )
00412      $         INFO = K
00413             KP = K
00414          ELSE
00415             IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
00416 *
00417 *              no interchange, use 1-by-1 pivot block
00418 *
00419                KP = K
00420             ELSE
00421 *
00422 *              Copy column IMAX to column K+1 of W and update it
00423 *
00424                CALL CCOPY( IMAX-K, A( IMAX, K ), LDA, W( K, K+1 ), 1 )
00425                CALL CCOPY( N-IMAX+1, A( IMAX, IMAX ), 1, W( IMAX, K+1 ),
00426      $                     1 )
00427                CALL CGEMV( 'No transpose', N-K+1, K-1, -CONE, A( K, 1 ),
00428      $                     LDA, W( IMAX, 1 ), LDW, CONE, W( K, K+1 ),
00429      $                     1 )
00430 *
00431 *              JMAX is the column-index of the largest off-diagonal
00432 *              element in row IMAX, and ROWMAX is its absolute value
00433 *
00434                JMAX = K - 1 + ICAMAX( IMAX-K, W( K, K+1 ), 1 )
00435                ROWMAX = CABS1( W( JMAX, K+1 ) )
00436                IF( IMAX.LT.N ) THEN
00437                   JMAX = IMAX + ICAMAX( N-IMAX, W( IMAX+1, K+1 ), 1 )
00438                   ROWMAX = MAX( ROWMAX, CABS1( W( JMAX, K+1 ) ) )
00439                END IF
00440 *
00441                IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
00442 *
00443 *                 no interchange, use 1-by-1 pivot block
00444 *
00445                   KP = K
00446                ELSE IF( CABS1( W( IMAX, K+1 ) ).GE.ALPHA*ROWMAX ) THEN
00447 *
00448 *                 interchange rows and columns K and IMAX, use 1-by-1
00449 *                 pivot block
00450 *
00451                   KP = IMAX
00452 *
00453 *                 copy column K+1 of W to column K
00454 *
00455                   CALL CCOPY( N-K+1, W( K, K+1 ), 1, W( K, K ), 1 )
00456                ELSE
00457 *
00458 *                 interchange rows and columns K+1 and IMAX, use 2-by-2
00459 *                 pivot block
00460 *
00461                   KP = IMAX
00462                   KSTEP = 2
00463                END IF
00464             END IF
00465 *
00466             KK = K + KSTEP - 1
00467 *
00468 *           Updated column KP is already stored in column KK of W
00469 *
00470             IF( KP.NE.KK ) THEN
00471 *
00472 *              Copy non-updated column KK to column KP
00473 *
00474                A( KP, K ) = A( KK, K )
00475                CALL CCOPY( KP-K-1, A( K+1, KK ), 1, A( KP, K+1 ), LDA )
00476                CALL CCOPY( N-KP+1, A( KP, KK ), 1, A( KP, KP ), 1 )
00477 *
00478 *              Interchange rows KK and KP in first KK columns of A and W
00479 *
00480                CALL CSWAP( KK, A( KK, 1 ), LDA, A( KP, 1 ), LDA )
00481                CALL CSWAP( KK, W( KK, 1 ), LDW, W( KP, 1 ), LDW )
00482             END IF
00483 *
00484             IF( KSTEP.EQ.1 ) THEN
00485 *
00486 *              1-by-1 pivot block D(k): column k of W now holds
00487 *
00488 *              W(k) = L(k)*D(k)
00489 *
00490 *              where L(k) is the k-th column of L
00491 *
00492 *              Store L(k) in column k of A
00493 *
00494                CALL CCOPY( N-K+1, W( K, K ), 1, A( K, K ), 1 )
00495                IF( K.LT.N ) THEN
00496                   R1 = CONE / A( K, K )
00497                   CALL CSCAL( N-K, R1, A( K+1, K ), 1 )
00498                END IF
00499             ELSE
00500 *
00501 *              2-by-2 pivot block D(k): columns k and k+1 of W now hold
00502 *
00503 *              ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k)
00504 *
00505 *              where L(k) and L(k+1) are the k-th and (k+1)-th columns
00506 *              of L
00507 *
00508                IF( K.LT.N-1 ) THEN
00509 *
00510 *                 Store L(k) and L(k+1) in columns k and k+1 of A
00511 *
00512                   D21 = W( K+1, K )
00513                   D11 = W( K+1, K+1 ) / D21
00514                   D22 = W( K, K ) / D21
00515                   T = CONE / ( D11*D22-CONE )
00516                   D21 = T / D21
00517                   DO 80 J = K + 2, N
00518                      A( J, K ) = D21*( D11*W( J, K )-W( J, K+1 ) )
00519                      A( J, K+1 ) = D21*( D22*W( J, K+1 )-W( J, K ) )
00520    80             CONTINUE
00521                END IF
00522 *
00523 *              Copy D(k) to A
00524 *
00525                A( K, K ) = W( K, K )
00526                A( K+1, K ) = W( K+1, K )
00527                A( K+1, K+1 ) = W( K+1, K+1 )
00528             END IF
00529          END IF
00530 *
00531 *        Store details of the interchanges in IPIV
00532 *
00533          IF( KSTEP.EQ.1 ) THEN
00534             IPIV( K ) = KP
00535          ELSE
00536             IPIV( K ) = -KP
00537             IPIV( K+1 ) = -KP
00538          END IF
00539 *
00540 *        Increase K and return to the start of the main loop
00541 *
00542          K = K + KSTEP
00543          GO TO 70
00544 *
00545    90    CONTINUE
00546 *
00547 *        Update the lower triangle of A22 (= A(k:n,k:n)) as
00548 *
00549 *        A22 := A22 - L21*D*L21' = A22 - L21*W'
00550 *
00551 *        computing blocks of NB columns at a time
00552 *
00553          DO 110 J = K, N, NB
00554             JB = MIN( NB, N-J+1 )
00555 *
00556 *           Update the lower triangle of the diagonal block
00557 *
00558             DO 100 JJ = J, J + JB - 1
00559                CALL CGEMV( 'No transpose', J+JB-JJ, K-1, -CONE,
00560      $                     A( JJ, 1 ), LDA, W( JJ, 1 ), LDW, CONE,
00561      $                     A( JJ, JJ ), 1 )
00562   100       CONTINUE
00563 *
00564 *           Update the rectangular subdiagonal block
00565 *
00566             IF( J+JB.LE.N )
00567      $         CALL CGEMM( 'No transpose', 'Transpose', N-J-JB+1, JB,
00568      $                     K-1, -CONE, A( J+JB, 1 ), LDA, W( J, 1 ),
00569      $                     LDW, CONE, A( J+JB, J ), LDA )
00570   110    CONTINUE
00571 *
00572 *        Put L21 in standard form by partially undoing the interchanges
00573 *        in columns 1:k-1
00574 *
00575          J = K - 1
00576   120    CONTINUE
00577          JJ = J
00578          JP = IPIV( J )
00579          IF( JP.LT.0 ) THEN
00580             JP = -JP
00581             J = J - 1
00582          END IF
00583          J = J - 1
00584          IF( JP.NE.JJ .AND. J.GE.1 )
00585      $      CALL CSWAP( J, A( JP, 1 ), LDA, A( JJ, 1 ), LDA )
00586          IF( J.GE.1 )
00587      $      GO TO 120
00588 *
00589 *        Set KB to the number of columns factorized
00590 *
00591          KB = K - 1
00592 *
00593       END IF
00594       RETURN
00595 *
00596 *     End of CLASYF
00597 *
00598       END