File:  [local] / rpl / lapack / lapack / zlasyf.f
Revision 1.10: download - view: text, annotated - select for diffs - revision graph
Mon Nov 21 22:19:54 2011 UTC (12 years, 5 months ago) by bertrand
Branches: MAIN
CVS tags: rpl-4_1_8, rpl-4_1_7, rpl-4_1_6, rpl-4_1_5, rpl-4_1_4, HEAD
Cohérence

    1: *> \brief \b ZLASYF
    2: *
    3: *  =========== DOCUMENTATION ===========
    4: *
    5: * Online html documentation available at 
    6: *            http://www.netlib.org/lapack/explore-html/ 
    7: *
    8: *> \htmlonly
    9: *> Download ZLASYF + dependencies 
   10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlasyf.f"> 
   11: *> [TGZ]</a> 
   12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlasyf.f"> 
   13: *> [ZIP]</a> 
   14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlasyf.f"> 
   15: *> [TXT]</a>
   16: *> \endhtmlonly 
   17: *
   18: *  Definition:
   19: *  ===========
   20: *
   21: *       SUBROUTINE ZLASYF( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, INFO )
   22:    23: *       .. Scalar Arguments ..
   24: *       CHARACTER          UPLO
   25: *       INTEGER            INFO, KB, LDA, LDW, N, NB
   26: *       ..
   27: *       .. Array Arguments ..
   28: *       INTEGER            IPIV( * )
   29: *       COMPLEX*16         A( LDA, * ), W( LDW, * )
   30: *       ..
   31: *  
   32: *
   33: *> \par Purpose:
   34: *  =============
   35: *>
   36: *> \verbatim
   37: *>
   38: *> ZLASYF computes a partial factorization of a complex symmetric matrix
   39: *> A using the Bunch-Kaufman diagonal pivoting method. The partial
   40: *> factorization has the form:
   41: *>
   42: *> A  =  ( I  U12 ) ( A11  0  ) (  I       0    )  if UPLO = 'U', or:
   43: *>       ( 0  U22 ) (  0   D  ) ( U12**T U22**T )
   44: *>
   45: *> A  =  ( L11  0 ) ( D    0  ) ( L11**T L21**T )  if UPLO = 'L'
   46: *>       ( L21  I ) ( 0   A22 ) (  0       I    )
   47: *>
   48: *> where the order of D is at most NB. The actual order is returned in
   49: *> the argument KB, and is either NB or NB-1, or N if N <= NB.
   50: *> Note that U**T denotes the transpose of U.
   51: *>
   52: *> ZLASYF is an auxiliary routine called by ZSYTRF. It uses blocked code
   53: *> (calling Level 3 BLAS) to update the submatrix A11 (if UPLO = 'U') or
   54: *> A22 (if UPLO = 'L').
   55: *> \endverbatim
   56: *
   57: *  Arguments:
   58: *  ==========
   59: *
   60: *> \param[in] UPLO
   61: *> \verbatim
   62: *>          UPLO is CHARACTER*1
   63: *>          Specifies whether the upper or lower triangular part of the
   64: *>          symmetric matrix A is stored:
   65: *>          = 'U':  Upper triangular
   66: *>          = 'L':  Lower triangular
   67: *> \endverbatim
   68: *>
   69: *> \param[in] N
   70: *> \verbatim
   71: *>          N is INTEGER
   72: *>          The order of the matrix A.  N >= 0.
   73: *> \endverbatim
   74: *>
   75: *> \param[in] NB
   76: *> \verbatim
   77: *>          NB is INTEGER
   78: *>          The maximum number of columns of the matrix A that should be
   79: *>          factored.  NB should be at least 2 to allow for 2-by-2 pivot
   80: *>          blocks.
   81: *> \endverbatim
   82: *>
   83: *> \param[out] KB
   84: *> \verbatim
   85: *>          KB is INTEGER
   86: *>          The number of columns of A that were actually factored.
   87: *>          KB is either NB-1 or NB, or N if N <= NB.
   88: *> \endverbatim
   89: *>
   90: *> \param[in,out] A
   91: *> \verbatim
   92: *>          A is COMPLEX*16 array, dimension (LDA,N)
   93: *>          On entry, the symmetric matrix A.  If UPLO = 'U', the leading
   94: *>          n-by-n upper triangular part of A contains the upper
   95: *>          triangular part of the matrix A, and the strictly lower
   96: *>          triangular part of A is not referenced.  If UPLO = 'L', the
   97: *>          leading n-by-n lower triangular part of A contains the lower
   98: *>          triangular part of the matrix A, and the strictly upper
   99: *>          triangular part of A is not referenced.
  100: *>          On exit, A contains details of the partial factorization.
  101: *> \endverbatim
  102: *>
  103: *> \param[in] LDA
  104: *> \verbatim
  105: *>          LDA is INTEGER
  106: *>          The leading dimension of the array A.  LDA >= max(1,N).
  107: *> \endverbatim
  108: *>
  109: *> \param[out] IPIV
  110: *> \verbatim
  111: *>          IPIV is INTEGER array, dimension (N)
  112: *>          Details of the interchanges and the block structure of D.
  113: *>          If UPLO = 'U', only the last KB elements of IPIV are set;
  114: *>          if UPLO = 'L', only the first KB elements are set.
  115: *>
  116: *>          If IPIV(k) > 0, then rows and columns k and IPIV(k) were
  117: *>          interchanged and D(k,k) is a 1-by-1 diagonal block.
  118: *>          If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and
  119: *>          columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
  120: *>          is a 2-by-2 diagonal block.  If UPLO = 'L' and IPIV(k) =
  121: *>          IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were
  122: *>          interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
  123: *> \endverbatim
  124: *>
  125: *> \param[out] W
  126: *> \verbatim
  127: *>          W is COMPLEX*16 array, dimension (LDW,NB)
  128: *> \endverbatim
  129: *>
  130: *> \param[in] LDW
  131: *> \verbatim
  132: *>          LDW is INTEGER
  133: *>          The leading dimension of the array W.  LDW >= max(1,N).
  134: *> \endverbatim
  135: *>
  136: *> \param[out] INFO
  137: *> \verbatim
  138: *>          INFO is INTEGER
  139: *>          = 0: successful exit
  140: *>          > 0: if INFO = k, D(k,k) is exactly zero.  The factorization
  141: *>               has been completed, but the block diagonal matrix D is
  142: *>               exactly singular.
  143: *> \endverbatim
  144: *
  145: *  Authors:
  146: *  ========
  147: *
  148: *> \author Univ. of Tennessee 
  149: *> \author Univ. of California Berkeley 
  150: *> \author Univ. of Colorado Denver 
  151: *> \author NAG Ltd. 
  152: *
  153: *> \date November 2011
  154: *
  155: *> \ingroup complex16SYcomputational
  156: *
  157: *  =====================================================================
  158:       SUBROUTINE ZLASYF( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, INFO )
  159: *
  160: *  -- LAPACK computational routine (version 3.4.0) --
  161: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  162: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  163: *     November 2011
  164: *
  165: *     .. Scalar Arguments ..
  166:       CHARACTER          UPLO
  167:       INTEGER            INFO, KB, LDA, LDW, N, NB
  168: *     ..
  169: *     .. Array Arguments ..
  170:       INTEGER            IPIV( * )
  171:       COMPLEX*16         A( LDA, * ), W( LDW, * )
  172: *     ..
  173: *
  174: *  =====================================================================
  175: *
  176: *     .. Parameters ..
  177:       DOUBLE PRECISION   ZERO, ONE
  178:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
  179:       DOUBLE PRECISION   EIGHT, SEVTEN
  180:       PARAMETER          ( EIGHT = 8.0D+0, SEVTEN = 17.0D+0 )
  181:       COMPLEX*16         CONE
  182:       PARAMETER          ( CONE = ( 1.0D+0, 0.0D+0 ) )
  183: *     ..
  184: *     .. Local Scalars ..
  185:       INTEGER            IMAX, J, JB, JJ, JMAX, JP, K, KK, KKW, KP,
  186:      $                   KSTEP, KW
  187:       DOUBLE PRECISION   ABSAKK, ALPHA, COLMAX, ROWMAX
  188:       COMPLEX*16         D11, D21, D22, R1, T, Z
  189: *     ..
  190: *     .. External Functions ..
  191:       LOGICAL            LSAME
  192:       INTEGER            IZAMAX
  193:       EXTERNAL           LSAME, IZAMAX
  194: *     ..
  195: *     .. External Subroutines ..
  196:       EXTERNAL           ZCOPY, ZGEMM, ZGEMV, ZSCAL, ZSWAP
  197: *     ..
  198: *     .. Intrinsic Functions ..
  199:       INTRINSIC          ABS, DBLE, DIMAG, MAX, MIN, SQRT
  200: *     ..
  201: *     .. Statement Functions ..
  202:       DOUBLE PRECISION   CABS1
  203: *     ..
  204: *     .. Statement Function definitions ..
  205:       CABS1( Z ) = ABS( DBLE( Z ) ) + ABS( DIMAG( Z ) )
  206: *     ..
  207: *     .. Executable Statements ..
  208: *
  209:       INFO = 0
  210: *
  211: *     Initialize ALPHA for use in choosing pivot block size.
  212: *
  213:       ALPHA = ( ONE+SQRT( SEVTEN ) ) / EIGHT
  214: *
  215:       IF( LSAME( UPLO, 'U' ) ) THEN
  216: *
  217: *        Factorize the trailing columns of A using the upper triangle
  218: *        of A and working backwards, and compute the matrix W = U12*D
  219: *        for use in updating A11
  220: *
  221: *        K is the main loop index, decreasing from N in steps of 1 or 2
  222: *
  223: *        KW is the column of W which corresponds to column K of A
  224: *
  225:          K = N
  226:    10    CONTINUE
  227:          KW = NB + K - N
  228: *
  229: *        Exit from loop
  230: *
  231:          IF( ( K.LE.N-NB+1 .AND. NB.LT.N ) .OR. K.LT.1 )
  232:      $      GO TO 30
  233: *
  234: *        Copy column K of A to column KW of W and update it
  235: *
  236:          CALL ZCOPY( K, A( 1, K ), 1, W( 1, KW ), 1 )
  237:          IF( K.LT.N )
  238:      $      CALL ZGEMV( 'No transpose', K, N-K, -CONE, A( 1, K+1 ), LDA,
  239:      $                  W( K, KW+1 ), LDW, CONE, W( 1, KW ), 1 )
  240: *
  241:          KSTEP = 1
  242: *
  243: *        Determine rows and columns to be interchanged and whether
  244: *        a 1-by-1 or 2-by-2 pivot block will be used
  245: *
  246:          ABSAKK = CABS1( W( K, KW ) )
  247: *
  248: *        IMAX is the row-index of the largest off-diagonal element in
  249: *        column K, and COLMAX is its absolute value
  250: *
  251:          IF( K.GT.1 ) THEN
  252:             IMAX = IZAMAX( K-1, W( 1, KW ), 1 )
  253:             COLMAX = CABS1( W( IMAX, KW ) )
  254:          ELSE
  255:             COLMAX = ZERO
  256:          END IF
  257: *
  258:          IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
  259: *
  260: *           Column K is zero: set INFO and continue
  261: *
  262:             IF( INFO.EQ.0 )
  263:      $         INFO = K
  264:             KP = K
  265:          ELSE
  266:             IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
  267: *
  268: *              no interchange, use 1-by-1 pivot block
  269: *
  270:                KP = K
  271:             ELSE
  272: *
  273: *              Copy column IMAX to column KW-1 of W and update it
  274: *
  275:                CALL ZCOPY( IMAX, A( 1, IMAX ), 1, W( 1, KW-1 ), 1 )
  276:                CALL ZCOPY( K-IMAX, A( IMAX, IMAX+1 ), LDA,
  277:      $                     W( IMAX+1, KW-1 ), 1 )
  278:                IF( K.LT.N )
  279:      $            CALL ZGEMV( 'No transpose', K, N-K, -CONE,
  280:      $                        A( 1, K+1 ), LDA, W( IMAX, KW+1 ), LDW,
  281:      $                        CONE, W( 1, KW-1 ), 1 )
  282: *
  283: *              JMAX is the column-index of the largest off-diagonal
  284: *              element in row IMAX, and ROWMAX is its absolute value
  285: *
  286:                JMAX = IMAX + IZAMAX( K-IMAX, W( IMAX+1, KW-1 ), 1 )
  287:                ROWMAX = CABS1( W( JMAX, KW-1 ) )
  288:                IF( IMAX.GT.1 ) THEN
  289:                   JMAX = IZAMAX( IMAX-1, W( 1, KW-1 ), 1 )
  290:                   ROWMAX = MAX( ROWMAX, CABS1( W( JMAX, KW-1 ) ) )
  291:                END IF
  292: *
  293:                IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
  294: *
  295: *                 no interchange, use 1-by-1 pivot block
  296: *
  297:                   KP = K
  298:                ELSE IF( CABS1( W( IMAX, KW-1 ) ).GE.ALPHA*ROWMAX ) THEN
  299: *
  300: *                 interchange rows and columns K and IMAX, use 1-by-1
  301: *                 pivot block
  302: *
  303:                   KP = IMAX
  304: *
  305: *                 copy column KW-1 of W to column KW
  306: *
  307:                   CALL ZCOPY( K, W( 1, KW-1 ), 1, W( 1, KW ), 1 )
  308:                ELSE
  309: *
  310: *                 interchange rows and columns K-1 and IMAX, use 2-by-2
  311: *                 pivot block
  312: *
  313:                   KP = IMAX
  314:                   KSTEP = 2
  315:                END IF
  316:             END IF
  317: *
  318:             KK = K - KSTEP + 1
  319:             KKW = NB + KK - N
  320: *
  321: *           Updated column KP is already stored in column KKW of W
  322: *
  323:             IF( KP.NE.KK ) THEN
  324: *
  325: *              Copy non-updated column KK to column KP
  326: *
  327:                A( KP, K ) = A( KK, K )
  328:                CALL ZCOPY( K-1-KP, A( KP+1, KK ), 1, A( KP, KP+1 ),
  329:      $                     LDA )
  330:                CALL ZCOPY( KP, A( 1, KK ), 1, A( 1, KP ), 1 )
  331: *
  332: *              Interchange rows KK and KP in last KK columns of A and W
  333: *
  334:                CALL ZSWAP( N-KK+1, A( KK, KK ), LDA, A( KP, KK ), LDA )
  335:                CALL ZSWAP( N-KK+1, W( KK, KKW ), LDW, W( KP, KKW ),
  336:      $                     LDW )
  337:             END IF
  338: *
  339:             IF( KSTEP.EQ.1 ) THEN
  340: *
  341: *              1-by-1 pivot block D(k): column KW of W now holds
  342: *
  343: *              W(k) = U(k)*D(k)
  344: *
  345: *              where U(k) is the k-th column of U
  346: *
  347: *              Store U(k) in column k of A
  348: *
  349:                CALL ZCOPY( K, W( 1, KW ), 1, A( 1, K ), 1 )
  350:                R1 = CONE / A( K, K )
  351:                CALL ZSCAL( K-1, R1, A( 1, K ), 1 )
  352:             ELSE
  353: *
  354: *              2-by-2 pivot block D(k): columns KW and KW-1 of W now
  355: *              hold
  356: *
  357: *              ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
  358: *
  359: *              where U(k) and U(k-1) are the k-th and (k-1)-th columns
  360: *              of U
  361: *
  362:                IF( K.GT.2 ) THEN
  363: *
  364: *                 Store U(k) and U(k-1) in columns k and k-1 of A
  365: *
  366:                   D21 = W( K-1, KW )
  367:                   D11 = W( K, KW ) / D21
  368:                   D22 = W( K-1, KW-1 ) / D21
  369:                   T = CONE / ( D11*D22-CONE )
  370:                   D21 = T / D21
  371:                   DO 20 J = 1, K - 2
  372:                      A( J, K-1 ) = D21*( D11*W( J, KW-1 )-W( J, KW ) )
  373:                      A( J, K ) = D21*( D22*W( J, KW )-W( J, KW-1 ) )
  374:    20             CONTINUE
  375:                END IF
  376: *
  377: *              Copy D(k) to A
  378: *
  379:                A( K-1, K-1 ) = W( K-1, KW-1 )
  380:                A( K-1, K ) = W( K-1, KW )
  381:                A( K, K ) = W( K, KW )
  382:             END IF
  383:          END IF
  384: *
  385: *        Store details of the interchanges in IPIV
  386: *
  387:          IF( KSTEP.EQ.1 ) THEN
  388:             IPIV( K ) = KP
  389:          ELSE
  390:             IPIV( K ) = -KP
  391:             IPIV( K-1 ) = -KP
  392:          END IF
  393: *
  394: *        Decrease K and return to the start of the main loop
  395: *
  396:          K = K - KSTEP
  397:          GO TO 10
  398: *
  399:    30    CONTINUE
  400: *
  401: *        Update the upper triangle of A11 (= A(1:k,1:k)) as
  402: *
  403: *        A11 := A11 - U12*D*U12**T = A11 - U12*W**T
  404: *
  405: *        computing blocks of NB columns at a time
  406: *
  407:          DO 50 J = ( ( K-1 ) / NB )*NB + 1, 1, -NB
  408:             JB = MIN( NB, K-J+1 )
  409: *
  410: *           Update the upper triangle of the diagonal block
  411: *
  412:             DO 40 JJ = J, J + JB - 1
  413:                CALL ZGEMV( 'No transpose', JJ-J+1, N-K, -CONE,
  414:      $                     A( J, K+1 ), LDA, W( JJ, KW+1 ), LDW, CONE,
  415:      $                     A( J, JJ ), 1 )
  416:    40       CONTINUE
  417: *
  418: *           Update the rectangular superdiagonal block
  419: *
  420:             CALL ZGEMM( 'No transpose', 'Transpose', J-1, JB, N-K,
  421:      $                  -CONE, A( 1, K+1 ), LDA, W( J, KW+1 ), LDW,
  422:      $                  CONE, A( 1, J ), LDA )
  423:    50    CONTINUE
  424: *
  425: *        Put U12 in standard form by partially undoing the interchanges
  426: *        in columns k+1:n
  427: *
  428:          J = K + 1
  429:    60    CONTINUE
  430:          JJ = J
  431:          JP = IPIV( J )
  432:          IF( JP.LT.0 ) THEN
  433:             JP = -JP
  434:             J = J + 1
  435:          END IF
  436:          J = J + 1
  437:          IF( JP.NE.JJ .AND. J.LE.N )
  438:      $      CALL ZSWAP( N-J+1, A( JP, J ), LDA, A( JJ, J ), LDA )
  439:          IF( J.LE.N )
  440:      $      GO TO 60
  441: *
  442: *        Set KB to the number of columns factorized
  443: *
  444:          KB = N - K
  445: *
  446:       ELSE
  447: *
  448: *        Factorize the leading columns of A using the lower triangle
  449: *        of A and working forwards, and compute the matrix W = L21*D
  450: *        for use in updating A22
  451: *
  452: *        K is the main loop index, increasing from 1 in steps of 1 or 2
  453: *
  454:          K = 1
  455:    70    CONTINUE
  456: *
  457: *        Exit from loop
  458: *
  459:          IF( ( K.GE.NB .AND. NB.LT.N ) .OR. K.GT.N )
  460:      $      GO TO 90
  461: *
  462: *        Copy column K of A to column K of W and update it
  463: *
  464:          CALL ZCOPY( N-K+1, A( K, K ), 1, W( K, K ), 1 )
  465:          CALL ZGEMV( 'No transpose', N-K+1, K-1, -CONE, A( K, 1 ), LDA,
  466:      $               W( K, 1 ), LDW, CONE, W( K, K ), 1 )
  467: *
  468:          KSTEP = 1
  469: *
  470: *        Determine rows and columns to be interchanged and whether
  471: *        a 1-by-1 or 2-by-2 pivot block will be used
  472: *
  473:          ABSAKK = CABS1( W( K, K ) )
  474: *
  475: *        IMAX is the row-index of the largest off-diagonal element in
  476: *        column K, and COLMAX is its absolute value
  477: *
  478:          IF( K.LT.N ) THEN
  479:             IMAX = K + IZAMAX( N-K, W( K+1, K ), 1 )
  480:             COLMAX = CABS1( W( IMAX, K ) )
  481:          ELSE
  482:             COLMAX = ZERO
  483:          END IF
  484: *
  485:          IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
  486: *
  487: *           Column K is zero: set INFO and continue
  488: *
  489:             IF( INFO.EQ.0 )
  490:      $         INFO = K
  491:             KP = K
  492:          ELSE
  493:             IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
  494: *
  495: *              no interchange, use 1-by-1 pivot block
  496: *
  497:                KP = K
  498:             ELSE
  499: *
  500: *              Copy column IMAX to column K+1 of W and update it
  501: *
  502:                CALL ZCOPY( IMAX-K, A( IMAX, K ), LDA, W( K, K+1 ), 1 )
  503:                CALL ZCOPY( N-IMAX+1, A( IMAX, IMAX ), 1, W( IMAX, K+1 ),
  504:      $                     1 )
  505:                CALL ZGEMV( 'No transpose', N-K+1, K-1, -CONE, A( K, 1 ),
  506:      $                     LDA, W( IMAX, 1 ), LDW, CONE, W( K, K+1 ),
  507:      $                     1 )
  508: *
  509: *              JMAX is the column-index of the largest off-diagonal
  510: *              element in row IMAX, and ROWMAX is its absolute value
  511: *
  512:                JMAX = K - 1 + IZAMAX( IMAX-K, W( K, K+1 ), 1 )
  513:                ROWMAX = CABS1( W( JMAX, K+1 ) )
  514:                IF( IMAX.LT.N ) THEN
  515:                   JMAX = IMAX + IZAMAX( N-IMAX, W( IMAX+1, K+1 ), 1 )
  516:                   ROWMAX = MAX( ROWMAX, CABS1( W( JMAX, K+1 ) ) )
  517:                END IF
  518: *
  519:                IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
  520: *
  521: *                 no interchange, use 1-by-1 pivot block
  522: *
  523:                   KP = K
  524:                ELSE IF( CABS1( W( IMAX, K+1 ) ).GE.ALPHA*ROWMAX ) THEN
  525: *
  526: *                 interchange rows and columns K and IMAX, use 1-by-1
  527: *                 pivot block
  528: *
  529:                   KP = IMAX
  530: *
  531: *                 copy column K+1 of W to column K
  532: *
  533:                   CALL ZCOPY( N-K+1, W( K, K+1 ), 1, W( K, K ), 1 )
  534:                ELSE
  535: *
  536: *                 interchange rows and columns K+1 and IMAX, use 2-by-2
  537: *                 pivot block
  538: *
  539:                   KP = IMAX
  540:                   KSTEP = 2
  541:                END IF
  542:             END IF
  543: *
  544:             KK = K + KSTEP - 1
  545: *
  546: *           Updated column KP is already stored in column KK of W
  547: *
  548:             IF( KP.NE.KK ) THEN
  549: *
  550: *              Copy non-updated column KK to column KP
  551: *
  552:                A( KP, K ) = A( KK, K )
  553:                CALL ZCOPY( KP-K-1, A( K+1, KK ), 1, A( KP, K+1 ), LDA )
  554:                CALL ZCOPY( N-KP+1, A( KP, KK ), 1, A( KP, KP ), 1 )
  555: *
  556: *              Interchange rows KK and KP in first KK columns of A and W
  557: *
  558:                CALL ZSWAP( KK, A( KK, 1 ), LDA, A( KP, 1 ), LDA )
  559:                CALL ZSWAP( KK, W( KK, 1 ), LDW, W( KP, 1 ), LDW )
  560:             END IF
  561: *
  562:             IF( KSTEP.EQ.1 ) THEN
  563: *
  564: *              1-by-1 pivot block D(k): column k of W now holds
  565: *
  566: *              W(k) = L(k)*D(k)
  567: *
  568: *              where L(k) is the k-th column of L
  569: *
  570: *              Store L(k) in column k of A
  571: *
  572:                CALL ZCOPY( N-K+1, W( K, K ), 1, A( K, K ), 1 )
  573:                IF( K.LT.N ) THEN
  574:                   R1 = CONE / A( K, K )
  575:                   CALL ZSCAL( N-K, R1, A( K+1, K ), 1 )
  576:                END IF
  577:             ELSE
  578: *
  579: *              2-by-2 pivot block D(k): columns k and k+1 of W now hold
  580: *
  581: *              ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k)
  582: *
  583: *              where L(k) and L(k+1) are the k-th and (k+1)-th columns
  584: *              of L
  585: *
  586:                IF( K.LT.N-1 ) THEN
  587: *
  588: *                 Store L(k) and L(k+1) in columns k and k+1 of A
  589: *
  590:                   D21 = W( K+1, K )
  591:                   D11 = W( K+1, K+1 ) / D21
  592:                   D22 = W( K, K ) / D21
  593:                   T = CONE / ( D11*D22-CONE )
  594:                   D21 = T / D21
  595:                   DO 80 J = K + 2, N
  596:                      A( J, K ) = D21*( D11*W( J, K )-W( J, K+1 ) )
  597:                      A( J, K+1 ) = D21*( D22*W( J, K+1 )-W( J, K ) )
  598:    80             CONTINUE
  599:                END IF
  600: *
  601: *              Copy D(k) to A
  602: *
  603:                A( K, K ) = W( K, K )
  604:                A( K+1, K ) = W( K+1, K )
  605:                A( K+1, K+1 ) = W( K+1, K+1 )
  606:             END IF
  607:          END IF
  608: *
  609: *        Store details of the interchanges in IPIV
  610: *
  611:          IF( KSTEP.EQ.1 ) THEN
  612:             IPIV( K ) = KP
  613:          ELSE
  614:             IPIV( K ) = -KP
  615:             IPIV( K+1 ) = -KP
  616:          END IF
  617: *
  618: *        Increase K and return to the start of the main loop
  619: *
  620:          K = K + KSTEP
  621:          GO TO 70
  622: *
  623:    90    CONTINUE
  624: *
  625: *        Update the lower triangle of A22 (= A(k:n,k:n)) as
  626: *
  627: *        A22 := A22 - L21*D*L21**T = A22 - L21*W**T
  628: *
  629: *        computing blocks of NB columns at a time
  630: *
  631:          DO 110 J = K, N, NB
  632:             JB = MIN( NB, N-J+1 )
  633: *
  634: *           Update the lower triangle of the diagonal block
  635: *
  636:             DO 100 JJ = J, J + JB - 1
  637:                CALL ZGEMV( 'No transpose', J+JB-JJ, K-1, -CONE,
  638:      $                     A( JJ, 1 ), LDA, W( JJ, 1 ), LDW, CONE,
  639:      $                     A( JJ, JJ ), 1 )
  640:   100       CONTINUE
  641: *
  642: *           Update the rectangular subdiagonal block
  643: *
  644:             IF( J+JB.LE.N )
  645:      $         CALL ZGEMM( 'No transpose', 'Transpose', N-J-JB+1, JB,
  646:      $                     K-1, -CONE, A( J+JB, 1 ), LDA, W( J, 1 ),
  647:      $                     LDW, CONE, A( J+JB, J ), LDA )
  648:   110    CONTINUE
  649: *
  650: *        Put L21 in standard form by partially undoing the interchanges
  651: *        in columns 1:k-1
  652: *
  653:          J = K - 1
  654:   120    CONTINUE
  655:          JJ = J
  656:          JP = IPIV( J )
  657:          IF( JP.LT.0 ) THEN
  658:             JP = -JP
  659:             J = J - 1
  660:          END IF
  661:          J = J - 1
  662:          IF( JP.NE.JJ .AND. J.GE.1 )
  663:      $      CALL ZSWAP( J, A( JP, 1 ), LDA, A( JJ, 1 ), LDA )
  664:          IF( J.GE.1 )
  665:      $      GO TO 120
  666: *
  667: *        Set KB to the number of columns factorized
  668: *
  669:          KB = K - 1
  670: *
  671:       END IF
  672:       RETURN
  673: *
  674: *     End of ZLASYF
  675: *
  676:       END

CVSweb interface <joel.bertrand@systella.fr>