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

CVSweb interface <joel.bertrand@systella.fr>