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