LAPACK 3.3.1 Linear Algebra PACKage

# slasyf.f

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