clasyf_rook()

subroutine clasyf_rook	(	character	uplo,
		integer	n,
		integer	nb,
		integer	kb,
		complex, dimension( lda, * )	a,
		integer	lda,
		integer, dimension( * )	ipiv,
		complex, dimension( ldw, * )	w,
		integer	ldw,
		integer	info
	)

CLASYF_ROOK computes a partial factorization of a complex symmetric matrix using the bounded Bunch-Kaufman ("rook") diagonal pivoting method.

Download CLASYF_ROOK + dependencies [TGZ] [ZIP] [TXT]

Purpose:

 CLASYF_ROOK computes a partial factorization of a complex symmetric
 matrix A using the bounded Bunch-Kaufman ("rook") diagonal
 pivoting method. The partial factorization has the form:

 A  =  ( I  U12 ) ( A11  0  ) (  I       0    )  if UPLO = 'U', or:
       ( 0  U22 ) (  0   D  ) ( U12**T U22**T )

 A  =  ( L11  0 ) (  D   0  ) ( L11**T L21**T )  if UPLO = 'L'
       ( L21  I ) (  0  A22 ) (  0       I    )

 where the order of D is at most NB. The actual order is returned in
 the argument KB, and is either NB or NB-1, or N if N <= NB.

 CLASYF_ROOK is an auxiliary routine called by CSYTRF_ROOK. It uses
 blocked code (calling Level 3 BLAS) to update the submatrix
 A11 (if UPLO = 'U') or A22 (if UPLO = 'L').

Parameters

[in]	UPLO	UPLO is CHARACTER*1 Specifies whether the upper or lower triangular part of the symmetric matrix A is stored: = 'U': Upper triangular = 'L': Lower triangular
[in]	N	N is INTEGER The order of the matrix A. N >= 0.
[in]	NB	NB is INTEGER The maximum number of columns of the matrix A that should be factored. NB should be at least 2 to allow for 2-by-2 pivot blocks.
[out]	KB	KB is INTEGER The number of columns of A that were actually factored. KB is either NB-1 or NB, or N if N <= NB.
[in,out]	A	A is COMPLEX array, dimension (LDA,N) On entry, the symmetric matrix A. If UPLO = 'U', the leading n-by-n upper triangular part of A contains the upper triangular part of the matrix A, and the strictly lower triangular part of A is not referenced. If UPLO = 'L', the leading n-by-n lower triangular part of A contains the lower triangular part of the matrix A, and the strictly upper triangular part of A is not referenced. On exit, A contains details of the partial factorization.
[in]	LDA	LDA is INTEGER The leading dimension of the array A. LDA >= max(1,N).
[out]	IPIV	IPIV is INTEGER array, dimension (N) Details of the interchanges and the block structure of D. If UPLO = 'U': Only the last KB elements of IPIV are set. If IPIV(k) > 0, then rows and columns k and IPIV(k) were interchanged and D(k,k) is a 1-by-1 diagonal block. If IPIV(k) < 0 and IPIV(k-1) < 0, then rows and columns k and -IPIV(k) were interchanged and rows and columns k-1 and -IPIV(k-1) were inerchaged, D(k-1:k,k-1:k) is a 2-by-2 diagonal block. If UPLO = 'L': Only the first KB elements of IPIV are set. If IPIV(k) > 0, then rows and columns k and IPIV(k) were interchanged and D(k,k) is a 1-by-1 diagonal block. If IPIV(k) < 0 and IPIV(k+1) < 0, then rows and columns k and -IPIV(k) were interchanged and rows and columns k+1 and -IPIV(k+1) were inerchaged, D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
[out]	W	W is COMPLEX array, dimension (LDW,NB)
[in]	LDW	LDW is INTEGER The leading dimension of the array W. LDW >= max(1,N).
[out]	INFO	INFO is INTEGER = 0: successful exit > 0: if INFO = k, D(k,k) is exactly zero. The factorization has been completed, but the block diagonal matrix D is exactly singular.

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Contributors:

  November 2013,     Igor Kozachenko,
                  Computer Science Division,
                  University of California, Berkeley

  September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas,
                  School of Mathematics,
                  University of Manchester

Definition at line 182 of file clasyf_rook.f.

*
*  -- LAPACK computational routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      CHARACTER          UPLO
      INTEGER            INFO, KB, LDA, LDW, N, NB
*     ..
*     .. Array Arguments ..
      INTEGER            IPIV( * )
      COMPLEX            A( LDA, * ), W( LDW, * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      REAL               ZERO, ONE
      parameter( zero = 0.0e+0, one = 1.0e+0 )
      REAL               EIGHT, SEVTEN
      parameter( eight = 8.0e+0, sevten = 17.0e+0 )
      COMPLEX            CONE, CZERO
      parameter( cone = ( 1.0e+0, 0.0e+0 ),
     $                   czero = ( 0.0e+0, 0.0e+0 ) )
*     ..
*     .. Local Scalars ..
      LOGICAL            DONE
      INTEGER            IMAX, ITEMP, J, JB, JJ, JMAX, JP1, JP2, K, KK,
     $                   KW, KKW, KP, KSTEP, P, II
      REAL               ABSAKK, ALPHA, COLMAX, ROWMAX, STEMP, SFMIN
      COMPLEX            D11, D12, D21, D22, R1, T, Z
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      INTEGER            ICAMAX
      REAL               SLAMCH
      EXTERNAL           lsame, icamax, slamch
*     ..
*     .. External Subroutines ..
      EXTERNAL           ccopy, cgemm, cgemv, cscal, cswap
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, max, min, sqrt, aimag, real
*     ..
*     .. Statement Functions ..
      REAL               CABS1
*     ..
*     .. Statement Function definitions ..
      cabs1( z ) = abs( real( z ) ) + abs( aimag( z ) )
*     ..
*     .. Executable Statements ..
*
      info = 0
*
*     Initialize ALPHA for use in choosing pivot block size.
*
      alpha = ( one+sqrt( sevten ) ) / eight
*
*     Compute machine safe minimum
*
      sfmin = slamch( 'S' )
*
      IF( lsame( uplo, 'U' ) ) THEN
*
*        Factorize the trailing columns of A using the upper triangle
*        of A and working backwards, and compute the matrix W = U12*D
*        for use in updating A11
*
*        K is the main loop index, decreasing from N in steps of 1 or 2
*
         k = n
   10    CONTINUE
*
*        KW is the column of W which corresponds to column K of A
*
         kw = nb + k - n
*
*        Exit from loop
*
         IF( ( k.LE.n-nb+1 .AND. nb.LT.n ) .OR. k.LT.1 )
     $      GO TO 30
*
         kstep = 1
         p = k
*
*        Copy column K of A to column KW of W and update it
*
         CALL ccopy( k, a( 1, k ), 1, w( 1, kw ), 1 )
         IF( k.LT.n )
     $      CALL cgemv( 'No transpose', k, n-k, -cone, a( 1, k+1 ),
     $                  lda, w( k, kw+1 ), ldw, cone, w( 1, kw ), 1 )
*
*        Determine rows and columns to be interchanged and whether
*        a 1-by-1 or 2-by-2 pivot block will be used
*
         absakk = cabs1( w( k, kw ) )
*
*        IMAX is the row-index of the largest off-diagonal element in
*        column K, and COLMAX is its absolute value.
*        Determine both COLMAX and IMAX.
*
         IF( k.GT.1 ) THEN
            imax = icamax( k-1, w( 1, kw ), 1 )
            colmax = cabs1( w( imax, kw ) )
         ELSE
            colmax = zero
         END IF
*
         IF( max( absakk, colmax ).EQ.zero ) THEN
*
*           Column K is zero or underflow: set INFO and continue
*
            IF( info.EQ.0 )
     $         info = k
            kp = k
            CALL ccopy( k, w( 1, kw ), 1, a( 1, k ), 1 )
         ELSE
*
*           ============================================================
*
*           Test for interchange
*
*           Equivalent to testing for ABSAKK.GE.ALPHA*COLMAX
*           (used to handle NaN and Inf)
*
            IF( .NOT.( absakk.LT.alpha*colmax ) ) THEN
*
*              no interchange, use 1-by-1 pivot block
*
               kp = k
*
            ELSE
*
               done = .false.
*
*              Loop until pivot found
*
   12          CONTINUE
*
*                 Begin pivot search loop body
*
*
*                 Copy column IMAX to column KW-1 of W and update it
*
                  CALL ccopy( imax, a( 1, imax ), 1, w( 1, kw-1 ), 1 )
                  CALL ccopy( k-imax, a( imax, imax+1 ), lda,
     $                        w( imax+1, kw-1 ), 1 )
*
                  IF( k.LT.n )
     $               CALL cgemv( 'No transpose', k, n-k, -cone,
     $                           a( 1, k+1 ), lda, w( imax, kw+1 ), ldw,
     $                           cone, w( 1, kw-1 ), 1 )
*
*                 JMAX is the column-index of the largest off-diagonal
*                 element in row IMAX, and ROWMAX is its absolute value.
*                 Determine both ROWMAX and JMAX.
*
                  IF( imax.NE.k ) THEN
                     jmax = imax + icamax( k-imax, w( imax+1, kw-1 ),
     $                                     1 )
                     rowmax = cabs1( w( jmax, kw-1 ) )
                  ELSE
                     rowmax = zero
                  END IF
*
                  IF( imax.GT.1 ) THEN
                     itemp = icamax( imax-1, w( 1, kw-1 ), 1 )
                     stemp = cabs1( w( itemp, kw-1 ) )
                     IF( stemp.GT.rowmax ) THEN
                        rowmax = stemp
                        jmax = itemp
                     END IF
                  END IF
*
*                 Equivalent to testing for
*                 CABS1( W( IMAX, KW-1 ) ).GE.ALPHA*ROWMAX
*                 (used to handle NaN and Inf)
*
                  IF( .NOT.(cabs1( w( imax, kw-1 ) ).LT.alpha*rowmax ) )
     $            THEN
*
*                    interchange rows and columns K and IMAX,
*                    use 1-by-1 pivot block
*
                     kp = imax
*
*                    copy column KW-1 of W to column KW of W
*
                     CALL ccopy( k, w( 1, kw-1 ), 1, w( 1, kw ), 1 )
*
                     done = .true.
*
*                 Equivalent to testing for ROWMAX.EQ.COLMAX,
*                 (used to handle NaN and Inf)
*
                  ELSE IF( ( p.EQ.jmax ) .OR. ( rowmax.LE.colmax ) )
     $            THEN
*
*                    interchange rows and columns K-1 and IMAX,
*                    use 2-by-2 pivot block
*
                     kp = imax
                     kstep = 2
                     done = .true.
                  ELSE
*
*                    Pivot not found: set params and repeat
*
                     p = imax
                     colmax = rowmax
                     imax = jmax
*
*                    Copy updated JMAXth (next IMAXth) column to Kth of W
*
                     CALL ccopy( k, w( 1, kw-1 ), 1, w( 1, kw ), 1 )
*
                  END IF
*
*                 End pivot search loop body
*
               IF( .NOT. done ) GOTO 12
*
            END IF
*
*           ============================================================
*
            kk = k - kstep + 1
*
*           KKW is the column of W which corresponds to column KK of A
*
            kkw = nb + kk - n
*
            IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) ) THEN
*
*              Copy non-updated column K to column P
*
               CALL ccopy( k-p, a( p+1, k ), 1, a( p, p+1 ), lda )
               CALL ccopy( p, a( 1, k ), 1, a( 1, p ), 1 )
*
*              Interchange rows K and P in last N-K+1 columns of A
*              and last N-K+2 columns of W
*
               CALL cswap( n-k+1, a( k, k ), lda, a( p, k ), lda )
               CALL cswap( n-kk+1, w( k, kkw ), ldw, w( p, kkw ), ldw )
            END IF
*
*           Updated column KP is already stored in column KKW of W
*
            IF( kp.NE.kk ) THEN
*
*              Copy non-updated column KK to column KP
*
               a( kp, k ) = a( kk, k )
               CALL ccopy( k-1-kp, a( kp+1, kk ), 1, a( kp, kp+1 ),
     $                     lda )
               CALL ccopy( kp, a( 1, kk ), 1, a( 1, kp ), 1 )
*
*              Interchange rows KK and KP in last N-KK+1 columns
*              of A and W
*
               CALL cswap( n-kk+1, a( kk, kk ), lda, a( kp, kk ), lda )
               CALL cswap( n-kk+1, w( kk, kkw ), ldw, w( kp, kkw ),
     $                     ldw )
            END IF
*
            IF( kstep.EQ.1 ) THEN
*
*              1-by-1 pivot block D(k): column KW of W now holds
*
*              W(k) = U(k)*D(k)
*
*              where U(k) is the k-th column of U
*
*              Store U(k) in column k of A
*
               CALL ccopy( k, w( 1, kw ), 1, a( 1, k ), 1 )
               IF( k.GT.1 ) THEN
                  IF( cabs1( a( k, k ) ).GE.sfmin ) THEN
                     r1 = cone / a( k, k )
                     CALL cscal( k-1, r1, a( 1, k ), 1 )
                  ELSE IF( a( k, k ).NE.czero ) THEN
                     DO 14 ii = 1, k - 1
                        a( ii, k ) = a( ii, k ) / a( k, k )
   14                CONTINUE
                  END IF
               END IF
*
            ELSE
*
*              2-by-2 pivot block D(k): columns KW and KW-1 of W now
*              hold
*
*              ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
*
*              where U(k) and U(k-1) are the k-th and (k-1)-th columns
*              of U
*
               IF( k.GT.2 ) THEN
*
*                 Store U(k) and U(k-1) in columns k and k-1 of A
*
                  d12 = w( k-1, kw )
                  d11 = w( k, kw ) / d12
                  d22 = w( k-1, kw-1 ) / d12
                  t = cone / ( d11*d22-cone )
                  DO 20 j = 1, k - 2
                     a( j, k-1 ) = t*( (d11*w( j, kw-1 )-w( j, kw ) ) /
     $                             d12 )
                     a( j, k ) = t*( ( d22*w( j, kw )-w( j, kw-1 ) ) /
     $                           d12 )
   20             CONTINUE
               END IF
*
*              Copy D(k) to A
*
               a( k-1, k-1 ) = w( k-1, kw-1 )
               a( k-1, k ) = w( k-1, kw )
               a( k, k ) = w( k, kw )
            END IF
         END IF
*
*        Store details of the interchanges in IPIV
*
         IF( kstep.EQ.1 ) THEN
            ipiv( k ) = kp
         ELSE
            ipiv( k ) = -p
            ipiv( k-1 ) = -kp
         END IF
*
*        Decrease K and return to the start of the main loop
*
         k = k - kstep
         GO TO 10
*
   30    CONTINUE
*
*        Update the upper triangle of A11 (= A(1:k,1:k)) as
*
*        A11 := A11 - U12*D*U12**T = A11 - U12*W**T
*
*        computing blocks of NB columns at a time
*
         DO 50 j = ( ( k-1 ) / nb )*nb + 1, 1, -nb
            jb = min( nb, k-j+1 )
*
*           Update the upper triangle of the diagonal block
*
            DO 40 jj = j, j + jb - 1
               CALL cgemv( 'No transpose', jj-j+1, n-k, -cone,
     $                     a( j, k+1 ), lda, w( jj, kw+1 ), ldw, cone,
     $                     a( j, jj ), 1 )
   40       CONTINUE
*
*           Update the rectangular superdiagonal block
*
            IF( j.GE.2 )
     $         CALL cgemm( 'No transpose', 'Transpose', j-1, jb,
     $                  n-k, -cone, a( 1, k+1 ), lda, w( j, kw+1 ), ldw,
     $                  cone, a( 1, j ), lda )
   50    CONTINUE
*
*        Put U12 in standard form by partially undoing the interchanges
*        in columns k+1:n
*
         j = k + 1
   60    CONTINUE
*
            kstep = 1
            jp1 = 1
            jj = j
            jp2 = ipiv( j )
            IF( jp2.LT.0 ) THEN
               jp2 = -jp2
               j = j + 1
               jp1 = -ipiv( j )
               kstep = 2
            END IF
*
            j = j + 1
            IF( jp2.NE.jj .AND. j.LE.n )
     $         CALL cswap( n-j+1, a( jp2, j ), lda, a( jj, j ), lda )
            jj = j - 1
            IF( jp1.NE.jj .AND. kstep.EQ.2 )
     $         CALL cswap( n-j+1, a( jp1, j ), lda, a( jj, j ), lda )
         IF( j.LE.n )
     $      GO TO 60
*
*        Set KB to the number of columns factorized
*
         kb = n - k
*
      ELSE
*
*        Factorize the leading columns of A using the lower triangle
*        of A and working forwards, and compute the matrix W = L21*D
*        for use in updating A22
*
*        K is the main loop index, increasing from 1 in steps of 1 or 2
*
         k = 1
   70   CONTINUE
*
*        Exit from loop
*
         IF( ( k.GE.nb .AND. nb.LT.n ) .OR. k.GT.n )
     $      GO TO 90
*
         kstep = 1
         p = k
*
*        Copy column K of A to column K of W and update it
*
         CALL ccopy( n-k+1, a( k, k ), 1, w( k, k ), 1 )
         IF( k.GT.1 )
     $      CALL cgemv( 'No transpose', n-k+1, k-1, -cone, a( k, 1 ),
     $                  lda, w( k, 1 ), ldw, cone, w( k, k ), 1 )
*
*        Determine rows and columns to be interchanged and whether
*        a 1-by-1 or 2-by-2 pivot block will be used
*
         absakk = cabs1( w( k, k ) )
*
*        IMAX is the row-index of the largest off-diagonal element in
*        column K, and COLMAX is its absolute value.
*        Determine both COLMAX and IMAX.
*
         IF( k.LT.n ) THEN
            imax = k + icamax( n-k, w( k+1, k ), 1 )
            colmax = cabs1( w( imax, k ) )
         ELSE
            colmax = zero
         END IF
*
         IF( max( absakk, colmax ).EQ.zero ) THEN
*
*           Column K is zero or underflow: set INFO and continue
*
            IF( info.EQ.0 )
     $         info = k
            kp = k
            CALL ccopy( n-k+1, w( k, k ), 1, a( k, k ), 1 )
         ELSE
*
*           ============================================================
*
*           Test for interchange
*
*           Equivalent to testing for ABSAKK.GE.ALPHA*COLMAX
*           (used to handle NaN and Inf)
*
            IF( .NOT.( absakk.LT.alpha*colmax ) ) THEN
*
*              no interchange, use 1-by-1 pivot block
*
               kp = k
*
            ELSE
*
               done = .false.
*
*              Loop until pivot found
*
   72          CONTINUE
*
*                 Begin pivot search loop body
*
*
*                 Copy column IMAX to column K+1 of W and update it
*
                  CALL ccopy( imax-k, a( imax, k ), lda, w( k, k+1 ), 1)
                  CALL ccopy( n-imax+1, a( imax, imax ), 1,
     $                        w( imax, k+1 ), 1 )
                  IF( k.GT.1 )
     $               CALL cgemv( 'No transpose', n-k+1, k-1, -cone,
     $                           a( k, 1 ), lda, w( imax, 1 ), ldw,
     $                           cone, w( k, k+1 ), 1 )
*
*                 JMAX is the column-index of the largest off-diagonal
*                 element in row IMAX, and ROWMAX is its absolute value.
*                 Determine both ROWMAX and JMAX.
*
                  IF( imax.NE.k ) THEN
                     jmax = k - 1 + icamax( imax-k, w( k, k+1 ), 1 )
                     rowmax = cabs1( w( jmax, k+1 ) )
                  ELSE
                     rowmax = zero
                  END IF
*
                  IF( imax.LT.n ) THEN
                     itemp = imax + icamax( n-imax, w( imax+1, k+1 ), 1)
                     stemp = cabs1( w( itemp, k+1 ) )
                     IF( stemp.GT.rowmax ) THEN
                        rowmax = stemp
                        jmax = itemp
                     END IF
                  END IF
*
*                 Equivalent to testing for
*                 CABS1( W( IMAX, K+1 ) ).GE.ALPHA*ROWMAX
*                 (used to handle NaN and Inf)
*
                  IF( .NOT.( cabs1( w( imax, k+1 ) ).LT.alpha*rowmax ) )
     $            THEN
*
*                    interchange rows and columns K and IMAX,
*                    use 1-by-1 pivot block
*
                     kp = imax
*
*                    copy column K+1 of W to column K of W
*
                     CALL ccopy( n-k+1, w( k, k+1 ), 1, w( k, k ), 1 )
*
                     done = .true.
*
*                 Equivalent to testing for ROWMAX.EQ.COLMAX,
*                 (used to handle NaN and Inf)
*
                  ELSE IF( ( p.EQ.jmax ) .OR. ( rowmax.LE.colmax ) )
     $            THEN
*
*                    interchange rows and columns K+1 and IMAX,
*                    use 2-by-2 pivot block
*
                     kp = imax
                     kstep = 2
                     done = .true.
                  ELSE
*
*                    Pivot not found: set params and repeat
*
                     p = imax
                     colmax = rowmax
                     imax = jmax
*
*                    Copy updated JMAXth (next IMAXth) column to Kth of W
*
                     CALL ccopy( n-k+1, w( k, k+1 ), 1, w( k, k ), 1 )
*
                  END IF
*
*                 End pivot search loop body
*
               IF( .NOT. done ) GOTO 72
*
            END IF
*
*           ============================================================
*
            kk = k + kstep - 1
*
            IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) ) THEN
*
*              Copy non-updated column K to column P
*
               CALL ccopy( p-k, a( k, k ), 1, a( p, k ), lda )
               CALL ccopy( n-p+1, a( p, k ), 1, a( p, p ), 1 )
*
*              Interchange rows K and P in first K columns of A
*              and first K+1 columns of W
*
               CALL cswap( k, a( k, 1 ), lda, a( p, 1 ), lda )
               CALL cswap( kk, w( k, 1 ), ldw, w( p, 1 ), ldw )
            END IF
*
*           Updated column KP is already stored in column KK of W
*
            IF( kp.NE.kk ) THEN
*
*              Copy non-updated column KK to column KP
*
               a( kp, k ) = a( kk, k )
               CALL ccopy( kp-k-1, a( k+1, kk ), 1, a( kp, k+1 ), lda )
               CALL ccopy( n-kp+1, a( kp, kk ), 1, a( kp, kp ), 1 )
*
*              Interchange rows KK and KP in first KK columns of A and W
*
               CALL cswap( kk, a( kk, 1 ), lda, a( kp, 1 ), lda )
               CALL cswap( kk, w( kk, 1 ), ldw, w( kp, 1 ), ldw )
            END IF
*
            IF( kstep.EQ.1 ) THEN
*
*              1-by-1 pivot block D(k): column k of W now holds
*
*              W(k) = L(k)*D(k)
*
*              where L(k) is the k-th column of L
*
*              Store L(k) in column k of A
*
               CALL ccopy( n-k+1, w( k, k ), 1, a( k, k ), 1 )
               IF( k.LT.n ) THEN
                  IF( cabs1( a( k, k ) ).GE.sfmin ) THEN
                     r1 = cone / a( k, k )
                     CALL cscal( n-k, r1, a( k+1, k ), 1 )
                  ELSE IF( a( k, k ).NE.czero ) THEN
                     DO 74 ii = k + 1, n
                        a( ii, k ) = a( ii, k ) / a( k, k )
   74                CONTINUE
                  END IF
               END IF
*
            ELSE
*
*              2-by-2 pivot block D(k): columns k and k+1 of W now hold
*
*              ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k)
*
*              where L(k) and L(k+1) are the k-th and (k+1)-th columns
*              of L
*
               IF( k.LT.n-1 ) THEN
*
*                 Store L(k) and L(k+1) in columns k and k+1 of A
*
                  d21 = w( k+1, k )
                  d11 = w( k+1, k+1 ) / d21
                  d22 = w( k, k ) / d21
                  t = cone / ( d11*d22-cone )
                  DO 80 j = k + 2, n
                     a( j, k ) = t*( ( d11*w( j, k )-w( j, k+1 ) ) /
     $                           d21 )
                     a( j, k+1 ) = t*( ( d22*w( j, k+1 )-w( j, k ) ) /
     $                             d21 )
   80             CONTINUE
               END IF
*
*              Copy D(k) to A
*
               a( k, k ) = w( k, k )
               a( k+1, k ) = w( k+1, k )
               a( k+1, k+1 ) = w( k+1, k+1 )
            END IF
         END IF
*
*        Store details of the interchanges in IPIV
*
         IF( kstep.EQ.1 ) THEN
            ipiv( k ) = kp
         ELSE
            ipiv( k ) = -p
            ipiv( k+1 ) = -kp
         END IF
*
*        Increase K and return to the start of the main loop
*
         k = k + kstep
         GO TO 70
*
   90    CONTINUE
*
*        Update the lower triangle of A22 (= A(k:n,k:n)) as
*
*        A22 := A22 - L21*D*L21**T = A22 - L21*W**T
*
*        computing blocks of NB columns at a time
*
         DO 110 j = k, n, nb
            jb = min( nb, n-j+1 )
*
*           Update the lower triangle of the diagonal block
*
            DO 100 jj = j, j + jb - 1
               CALL cgemv( 'No transpose', j+jb-jj, k-1, -cone,
     $                     a( jj, 1 ), lda, w( jj, 1 ), ldw, cone,
     $                     a( jj, jj ), 1 )
  100       CONTINUE
*
*           Update the rectangular subdiagonal block
*
            IF( j+jb.LE.n )
     $         CALL cgemm( 'No transpose', 'Transpose', n-j-jb+1, jb,
     $                    k-1, -cone, a( j+jb, 1 ), lda, w( j, 1 ), ldw,
     $                    cone, a( j+jb, j ), lda )
  110    CONTINUE
*
*        Put L21 in standard form by partially undoing the interchanges
*        in columns 1:k-1
*
         j = k - 1
  120    CONTINUE
*
            kstep = 1
            jp1 = 1
            jj = j
            jp2 = ipiv( j )
            IF( jp2.LT.0 ) THEN
               jp2 = -jp2
               j = j - 1
               jp1 = -ipiv( j )
               kstep = 2
            END IF
*
            j = j - 1
            IF( jp2.NE.jj .AND. j.GE.1 )
     $         CALL cswap( j, a( jp2, 1 ), lda, a( jj, 1 ), lda )
            jj = j + 1
            IF( jp1.NE.jj .AND. kstep.EQ.2 )
     $         CALL cswap( j, a( jp1, 1 ), lda, a( jj, 1 ), lda )
         IF( j.GE.1 )
     $      GO TO 120
*
*        Set KB to the number of columns factorized
*
         kb = k - 1
*
      END IF
      RETURN
*
*     End of CLASYF_ROOK
*

Here is the call graph for this function:

Here is the caller graph for this function: