LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
clasyf.f
Go to the documentation of this file.
00001 *> \brief \b CLASYF
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download CLASYF + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clasyf.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/clasyf.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/clasyf.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE CLASYF( 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 *> CLASYF computes a partial factorization of a complex symmetric matrix
00039 *> A 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 *> Note that U**T denotes the transpose of U.
00051 *>
00052 *> CLASYF is an auxiliary routine called by CSYTRF. 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 *>          symmetric 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 symmetric 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 complexSYcomputational
00156 *
00157 *  =====================================================================
00158       SUBROUTINE CLASYF( 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       REAL               EIGHT, SEVTEN
00180       PARAMETER          ( EIGHT = 8.0E+0, SEVTEN = 17.0E+0 )
00181       COMPLEX            CONE
00182       PARAMETER          ( CONE = ( 1.0E+0, 0.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, ROWMAX
00188       COMPLEX            D11, D21, D22, R1, T, Z
00189 *     ..
00190 *     .. External Functions ..
00191       LOGICAL            LSAME
00192       INTEGER            ICAMAX
00193       EXTERNAL           LSAME, ICAMAX
00194 *     ..
00195 *     .. External Subroutines ..
00196       EXTERNAL           CCOPY, CGEMM, CGEMV, CSCAL, CSWAP
00197 *     ..
00198 *     .. Intrinsic Functions ..
00199       INTRINSIC          ABS, AIMAG, 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
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, A( 1, K ), 1, W( 1, KW ), 1 )
00237          IF( K.LT.N )
00238      $      CALL CGEMV( 'No transpose', K, N-K, -CONE, A( 1, K+1 ), LDA,
00239      $                  W( K, KW+1 ), LDW, CONE, W( 1, KW ), 1 )
00240 *
00241          KSTEP = 1
00242 *
00243 *        Determine rows and columns to be interchanged and whether
00244 *        a 1-by-1 or 2-by-2 pivot block will be used
00245 *
00246          ABSAKK = CABS1( W( K, KW ) )
00247 *
00248 *        IMAX is the row-index of the largest off-diagonal element in
00249 *        column K, and COLMAX is its absolute value
00250 *
00251          IF( K.GT.1 ) THEN
00252             IMAX = ICAMAX( K-1, W( 1, KW ), 1 )
00253             COLMAX = CABS1( W( IMAX, KW ) )
00254          ELSE
00255             COLMAX = ZERO
00256          END IF
00257 *
00258          IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
00259 *
00260 *           Column K is zero: set INFO and continue
00261 *
00262             IF( INFO.EQ.0 )
00263      $         INFO = K
00264             KP = K
00265          ELSE
00266             IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
00267 *
00268 *              no interchange, use 1-by-1 pivot block
00269 *
00270                KP = K
00271             ELSE
00272 *
00273 *              Copy column IMAX to column KW-1 of W and update it
00274 *
00275                CALL CCOPY( IMAX, A( 1, IMAX ), 1, W( 1, KW-1 ), 1 )
00276                CALL CCOPY( K-IMAX, A( IMAX, IMAX+1 ), LDA,
00277      $                     W( IMAX+1, KW-1 ), 1 )
00278                IF( K.LT.N )
00279      $            CALL CGEMV( 'No transpose', K, N-K, -CONE,
00280      $                        A( 1, K+1 ), LDA, W( IMAX, KW+1 ), LDW,
00281      $                        CONE, W( 1, KW-1 ), 1 )
00282 *
00283 *              JMAX is the column-index of the largest off-diagonal
00284 *              element in row IMAX, and ROWMAX is its absolute value
00285 *
00286                JMAX = IMAX + ICAMAX( K-IMAX, W( IMAX+1, KW-1 ), 1 )
00287                ROWMAX = CABS1( W( JMAX, KW-1 ) )
00288                IF( IMAX.GT.1 ) THEN
00289                   JMAX = ICAMAX( IMAX-1, W( 1, KW-1 ), 1 )
00290                   ROWMAX = MAX( ROWMAX, CABS1( W( JMAX, KW-1 ) ) )
00291                END IF
00292 *
00293                IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
00294 *
00295 *                 no interchange, use 1-by-1 pivot block
00296 *
00297                   KP = K
00298                ELSE IF( CABS1( W( IMAX, KW-1 ) ).GE.ALPHA*ROWMAX ) THEN
00299 *
00300 *                 interchange rows and columns K and IMAX, use 1-by-1
00301 *                 pivot block
00302 *
00303                   KP = IMAX
00304 *
00305 *                 copy column KW-1 of W to column KW
00306 *
00307                   CALL CCOPY( K, W( 1, KW-1 ), 1, W( 1, KW ), 1 )
00308                ELSE
00309 *
00310 *                 interchange rows and columns K-1 and IMAX, use 2-by-2
00311 *                 pivot block
00312 *
00313                   KP = IMAX
00314                   KSTEP = 2
00315                END IF
00316             END IF
00317 *
00318             KK = K - KSTEP + 1
00319             KKW = NB + KK - N
00320 *
00321 *           Updated column KP is already stored in column KKW of W
00322 *
00323             IF( KP.NE.KK ) THEN
00324 *
00325 *              Copy non-updated column KK to column KP
00326 *
00327                A( KP, K ) = A( KK, K )
00328                CALL CCOPY( K-1-KP, A( KP+1, KK ), 1, A( KP, KP+1 ),
00329      $                     LDA )
00330                CALL CCOPY( KP, A( 1, KK ), 1, A( 1, KP ), 1 )
00331 *
00332 *              Interchange rows KK and KP in last KK columns of A and W
00333 *
00334                CALL CSWAP( N-KK+1, A( KK, KK ), LDA, A( KP, KK ), LDA )
00335                CALL CSWAP( N-KK+1, W( KK, KKW ), LDW, W( KP, KKW ),
00336      $                     LDW )
00337             END IF
00338 *
00339             IF( KSTEP.EQ.1 ) THEN
00340 *
00341 *              1-by-1 pivot block D(k): column KW of W now holds
00342 *
00343 *              W(k) = U(k)*D(k)
00344 *
00345 *              where U(k) is the k-th column of U
00346 *
00347 *              Store U(k) in column k of A
00348 *
00349                CALL CCOPY( K, W( 1, KW ), 1, A( 1, K ), 1 )
00350                R1 = CONE / A( K, K )
00351                CALL CSCAL( K-1, R1, A( 1, K ), 1 )
00352             ELSE
00353 *
00354 *              2-by-2 pivot block D(k): columns KW and KW-1 of W now
00355 *              hold
00356 *
00357 *              ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
00358 *
00359 *              where U(k) and U(k-1) are the k-th and (k-1)-th columns
00360 *              of U
00361 *
00362                IF( K.GT.2 ) THEN
00363 *
00364 *                 Store U(k) and U(k-1) in columns k and k-1 of A
00365 *
00366                   D21 = W( K-1, KW )
00367                   D11 = W( K, KW ) / D21
00368                   D22 = W( K-1, KW-1 ) / D21
00369                   T = CONE / ( D11*D22-CONE )
00370                   D21 = T / D21
00371                   DO 20 J = 1, K - 2
00372                      A( J, K-1 ) = D21*( D11*W( J, KW-1 )-W( J, KW ) )
00373                      A( J, K ) = D21*( D22*W( J, KW )-W( J, KW-1 ) )
00374    20             CONTINUE
00375                END IF
00376 *
00377 *              Copy D(k) to A
00378 *
00379                A( K-1, K-1 ) = W( K-1, KW-1 )
00380                A( K-1, K ) = W( K-1, KW )
00381                A( K, K ) = W( K, KW )
00382             END IF
00383          END IF
00384 *
00385 *        Store details of the interchanges in IPIV
00386 *
00387          IF( KSTEP.EQ.1 ) THEN
00388             IPIV( K ) = KP
00389          ELSE
00390             IPIV( K ) = -KP
00391             IPIV( K-1 ) = -KP
00392          END IF
00393 *
00394 *        Decrease K and return to the start of the main loop
00395 *
00396          K = K - KSTEP
00397          GO TO 10
00398 *
00399    30    CONTINUE
00400 *
00401 *        Update the upper triangle of A11 (= A(1:k,1:k)) as
00402 *
00403 *        A11 := A11 - U12*D*U12**T = A11 - U12*W**T
00404 *
00405 *        computing blocks of NB columns at a time
00406 *
00407          DO 50 J = ( ( K-1 ) / NB )*NB + 1, 1, -NB
00408             JB = MIN( NB, K-J+1 )
00409 *
00410 *           Update the upper triangle of the diagonal block
00411 *
00412             DO 40 JJ = J, J + JB - 1
00413                CALL CGEMV( 'No transpose', JJ-J+1, N-K, -CONE,
00414      $                     A( J, K+1 ), LDA, W( JJ, KW+1 ), LDW, CONE,
00415      $                     A( J, JJ ), 1 )
00416    40       CONTINUE
00417 *
00418 *           Update the rectangular superdiagonal block
00419 *
00420             CALL CGEMM( 'No transpose', 'Transpose', J-1, JB, N-K,
00421      $                  -CONE, A( 1, K+1 ), LDA, W( J, KW+1 ), LDW,
00422      $                  CONE, A( 1, J ), LDA )
00423    50    CONTINUE
00424 *
00425 *        Put U12 in standard form by partially undoing the interchanges
00426 *        in columns k+1:n
00427 *
00428          J = K + 1
00429    60    CONTINUE
00430          JJ = J
00431          JP = IPIV( J )
00432          IF( JP.LT.0 ) THEN
00433             JP = -JP
00434             J = J + 1
00435          END IF
00436          J = J + 1
00437          IF( JP.NE.JJ .AND. J.LE.N )
00438      $      CALL CSWAP( N-J+1, A( JP, J ), LDA, A( JJ, J ), LDA )
00439          IF( J.LE.N )
00440      $      GO TO 60
00441 *
00442 *        Set KB to the number of columns factorized
00443 *
00444          KB = N - K
00445 *
00446       ELSE
00447 *
00448 *        Factorize the leading columns of A using the lower triangle
00449 *        of A and working forwards, and compute the matrix W = L21*D
00450 *        for use in updating A22
00451 *
00452 *        K is the main loop index, increasing from 1 in steps of 1 or 2
00453 *
00454          K = 1
00455    70    CONTINUE
00456 *
00457 *        Exit from loop
00458 *
00459          IF( ( K.GE.NB .AND. NB.LT.N ) .OR. K.GT.N )
00460      $      GO TO 90
00461 *
00462 *        Copy column K of A to column K of W and update it
00463 *
00464          CALL CCOPY( N-K+1, A( K, K ), 1, W( K, K ), 1 )
00465          CALL CGEMV( 'No transpose', N-K+1, K-1, -CONE, A( K, 1 ), LDA,
00466      $               W( K, 1 ), LDW, CONE, W( K, K ), 1 )
00467 *
00468          KSTEP = 1
00469 *
00470 *        Determine rows and columns to be interchanged and whether
00471 *        a 1-by-1 or 2-by-2 pivot block will be used
00472 *
00473          ABSAKK = CABS1( W( K, K ) )
00474 *
00475 *        IMAX is the row-index of the largest off-diagonal element in
00476 *        column K, and COLMAX is its absolute value
00477 *
00478          IF( K.LT.N ) THEN
00479             IMAX = K + ICAMAX( N-K, W( K+1, K ), 1 )
00480             COLMAX = CABS1( W( IMAX, K ) )
00481          ELSE
00482             COLMAX = ZERO
00483          END IF
00484 *
00485          IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
00486 *
00487 *           Column K is zero: set INFO and continue
00488 *
00489             IF( INFO.EQ.0 )
00490      $         INFO = K
00491             KP = K
00492          ELSE
00493             IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
00494 *
00495 *              no interchange, use 1-by-1 pivot block
00496 *
00497                KP = K
00498             ELSE
00499 *
00500 *              Copy column IMAX to column K+1 of W and update it
00501 *
00502                CALL CCOPY( IMAX-K, A( IMAX, K ), LDA, W( K, K+1 ), 1 )
00503                CALL CCOPY( N-IMAX+1, A( IMAX, IMAX ), 1, W( IMAX, K+1 ),
00504      $                     1 )
00505                CALL CGEMV( 'No transpose', N-K+1, K-1, -CONE, A( K, 1 ),
00506      $                     LDA, W( IMAX, 1 ), LDW, CONE, W( K, K+1 ),
00507      $                     1 )
00508 *
00509 *              JMAX is the column-index of the largest off-diagonal
00510 *              element in row IMAX, and ROWMAX is its absolute value
00511 *
00512                JMAX = K - 1 + ICAMAX( IMAX-K, W( K, K+1 ), 1 )
00513                ROWMAX = CABS1( W( JMAX, K+1 ) )
00514                IF( IMAX.LT.N ) THEN
00515                   JMAX = IMAX + ICAMAX( N-IMAX, W( IMAX+1, K+1 ), 1 )
00516                   ROWMAX = MAX( ROWMAX, CABS1( W( JMAX, K+1 ) ) )
00517                END IF
00518 *
00519                IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
00520 *
00521 *                 no interchange, use 1-by-1 pivot block
00522 *
00523                   KP = K
00524                ELSE IF( CABS1( W( IMAX, K+1 ) ).GE.ALPHA*ROWMAX ) THEN
00525 *
00526 *                 interchange rows and columns K and IMAX, use 1-by-1
00527 *                 pivot block
00528 *
00529                   KP = IMAX
00530 *
00531 *                 copy column K+1 of W to column K
00532 *
00533                   CALL CCOPY( N-K+1, W( K, K+1 ), 1, W( K, K ), 1 )
00534                ELSE
00535 *
00536 *                 interchange rows and columns K+1 and IMAX, use 2-by-2
00537 *                 pivot block
00538 *
00539                   KP = IMAX
00540                   KSTEP = 2
00541                END IF
00542             END IF
00543 *
00544             KK = K + KSTEP - 1
00545 *
00546 *           Updated column KP is already stored in column KK of W
00547 *
00548             IF( KP.NE.KK ) THEN
00549 *
00550 *              Copy non-updated column KK to column KP
00551 *
00552                A( KP, K ) = A( KK, K )
00553                CALL CCOPY( KP-K-1, A( K+1, KK ), 1, A( KP, K+1 ), LDA )
00554                CALL CCOPY( N-KP+1, A( KP, KK ), 1, A( KP, KP ), 1 )
00555 *
00556 *              Interchange rows KK and KP in first KK columns of A and W
00557 *
00558                CALL CSWAP( KK, A( KK, 1 ), LDA, A( KP, 1 ), LDA )
00559                CALL CSWAP( KK, W( KK, 1 ), LDW, W( KP, 1 ), LDW )
00560             END IF
00561 *
00562             IF( KSTEP.EQ.1 ) THEN
00563 *
00564 *              1-by-1 pivot block D(k): column k of W now holds
00565 *
00566 *              W(k) = L(k)*D(k)
00567 *
00568 *              where L(k) is the k-th column of L
00569 *
00570 *              Store L(k) in column k of A
00571 *
00572                CALL CCOPY( N-K+1, W( K, K ), 1, A( K, K ), 1 )
00573                IF( K.LT.N ) THEN
00574                   R1 = CONE / A( K, K )
00575                   CALL CSCAL( N-K, R1, A( K+1, K ), 1 )
00576                END IF
00577             ELSE
00578 *
00579 *              2-by-2 pivot block D(k): columns k and k+1 of W now hold
00580 *
00581 *              ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k)
00582 *
00583 *              where L(k) and L(k+1) are the k-th and (k+1)-th columns
00584 *              of L
00585 *
00586                IF( K.LT.N-1 ) THEN
00587 *
00588 *                 Store L(k) and L(k+1) in columns k and k+1 of A
00589 *
00590                   D21 = W( K+1, K )
00591                   D11 = W( K+1, K+1 ) / D21
00592                   D22 = W( K, K ) / D21
00593                   T = CONE / ( D11*D22-CONE )
00594                   D21 = T / D21
00595                   DO 80 J = K + 2, N
00596                      A( J, K ) = D21*( D11*W( J, K )-W( J, K+1 ) )
00597                      A( J, K+1 ) = D21*( D22*W( J, K+1 )-W( J, K ) )
00598    80             CONTINUE
00599                END IF
00600 *
00601 *              Copy D(k) to A
00602 *
00603                A( K, K ) = W( K, K )
00604                A( K+1, K ) = W( K+1, K )
00605                A( K+1, K+1 ) = W( K+1, K+1 )
00606             END IF
00607          END IF
00608 *
00609 *        Store details of the interchanges in IPIV
00610 *
00611          IF( KSTEP.EQ.1 ) THEN
00612             IPIV( K ) = KP
00613          ELSE
00614             IPIV( K ) = -KP
00615             IPIV( K+1 ) = -KP
00616          END IF
00617 *
00618 *        Increase K and return to the start of the main loop
00619 *
00620          K = K + KSTEP
00621          GO TO 70
00622 *
00623    90    CONTINUE
00624 *
00625 *        Update the lower triangle of A22 (= A(k:n,k:n)) as
00626 *
00627 *        A22 := A22 - L21*D*L21**T = A22 - L21*W**T
00628 *
00629 *        computing blocks of NB columns at a time
00630 *
00631          DO 110 J = K, N, NB
00632             JB = MIN( NB, N-J+1 )
00633 *
00634 *           Update the lower triangle of the diagonal block
00635 *
00636             DO 100 JJ = J, J + JB - 1
00637                CALL CGEMV( 'No transpose', J+JB-JJ, K-1, -CONE,
00638      $                     A( JJ, 1 ), LDA, W( JJ, 1 ), LDW, CONE,
00639      $                     A( JJ, JJ ), 1 )
00640   100       CONTINUE
00641 *
00642 *           Update the rectangular subdiagonal block
00643 *
00644             IF( J+JB.LE.N )
00645      $         CALL CGEMM( 'No transpose', 'Transpose', N-J-JB+1, JB,
00646      $                     K-1, -CONE, A( J+JB, 1 ), LDA, W( J, 1 ),
00647      $                     LDW, CONE, A( J+JB, J ), LDA )
00648   110    CONTINUE
00649 *
00650 *        Put L21 in standard form by partially undoing the interchanges
00651 *        in columns 1:k-1
00652 *
00653          J = K - 1
00654   120    CONTINUE
00655          JJ = J
00656          JP = IPIV( J )
00657          IF( JP.LT.0 ) THEN
00658             JP = -JP
00659             J = J - 1
00660          END IF
00661          J = J - 1
00662          IF( JP.NE.JJ .AND. J.GE.1 )
00663      $      CALL CSWAP( J, A( JP, 1 ), LDA, A( JJ, 1 ), LDA )
00664          IF( J.GE.1 )
00665      $      GO TO 120
00666 *
00667 *        Set KB to the number of columns factorized
00668 *
00669          KB = K - 1
00670 *
00671       END IF
00672       RETURN
00673 *
00674 *     End of CLASYF
00675 *
00676       END
 All Files Functions