slasyf_rook()

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

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

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

Purpose:

 SLASYF_ROOK computes a partial factorization of a real 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.

 SLASYF_ROOK is an auxiliary routine called by SSYTRF_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 REAL 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 REAL 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 slasyf_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( * )
      REAL               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 )
*     ..
*     .. 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, D11, D12, D21, D22,
     $                   STEMP, R1, ROWMAX, T, SFMIN
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      INTEGER            ISAMAX
      REAL               SLAMCH
      EXTERNAL           lsame, isamax, slamch
*     ..
*     .. External Subroutines ..
      EXTERNAL           scopy, sgemm, sgemv, sscal, sswap
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, max, min, sqrt
*     ..
*     .. 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 scopy( k, a( 1, k ), 1, w( 1, kw ), 1 )
         IF( k.LT.n )
     $      CALL sgemv( 'No transpose', k, n-k, -one, a( 1, k+1 ),
     $                  lda, w( k, kw+1 ), ldw, one, 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 = abs( 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 = isamax( k-1, w( 1, kw ), 1 )
            colmax = abs( 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 scopy( 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 scopy( imax, a( 1, imax ), 1, w( 1, kw-1 ), 1 )
                  CALL scopy( k-imax, a( imax, imax+1 ), lda,
     $                        w( imax+1, kw-1 ), 1 )
*
                  IF( k.LT.n )
     $               CALL sgemv( 'No transpose', k, n-k, -one,
     $                           a( 1, k+1 ), lda, w( imax, kw+1 ), ldw,
     $                           one, 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 + isamax( k-imax, w( imax+1, kw-1 ),
     $                                     1 )
                     rowmax = abs( w( jmax, kw-1 ) )
                  ELSE
                     rowmax = zero
                  END IF
*
                  IF( imax.GT.1 ) THEN
                     itemp = isamax( imax-1, w( 1, kw-1 ), 1 )
                     stemp = abs( w( itemp, kw-1 ) )
                     IF( stemp.GT.rowmax ) THEN
                        rowmax = stemp
                        jmax = itemp
                     END IF
                  END IF
*
*                 Equivalent to testing for
*                 ABS( W( IMAX, KW-1 ) ).GE.ALPHA*ROWMAX
*                 (used to handle NaN and Inf)
*
                  IF( .NOT.(abs( 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 scopy( 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 scopy( 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 scopy( k-p, a( p+1, k ), 1, a( p, p+1 ), lda )
               CALL scopy( 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 sswap( n-k+1, a( k, k ), lda, a( p, k ), lda )
               CALL sswap( 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 scopy( k-1-kp, a( kp+1, kk ), 1, a( kp, kp+1 ),
     $                     lda )
               CALL scopy( 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 sswap( n-kk+1, a( kk, kk ), lda, a( kp, kk ), lda )
               CALL sswap( 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 scopy( k, w( 1, kw ), 1, a( 1, k ), 1 )
               IF( k.GT.1 ) THEN
                  IF( abs( a( k, k ) ).GE.sfmin ) THEN
                     r1 = one / a( k, k )
                     CALL sscal( k-1, r1, a( 1, k ), 1 )
                  ELSE IF( a( k, k ).NE.zero ) 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 = one / ( d11*d22-one )
                  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 sgemv( 'No transpose', jj-j+1, n-k, -one,
     $                     a( j, k+1 ), lda, w( jj, kw+1 ), ldw, one,
     $                     a( j, jj ), 1 )
   40       CONTINUE
*
*           Update the rectangular superdiagonal block
*
            IF( j.GE.2 )
     $         CALL sgemm( 'No transpose', 'Transpose', j-1, jb,
     $                  n-k, -one, a( 1, k+1 ), lda, w( j, kw+1 ), ldw,
     $                  one, 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 sswap( n-j+1, a( jp2, j ), lda, a( jj, j ), lda )
            jj = j - 1
            IF( jp1.NE.jj .AND. kstep.EQ.2 )
     $         CALL sswap( 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 scopy( n-k+1, a( k, k ), 1, w( k, k ), 1 )
         IF( k.GT.1 )
     $      CALL sgemv( 'No transpose', n-k+1, k-1, -one, a( k, 1 ),
     $                  lda, w( k, 1 ), ldw, one, 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 = abs( 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 + isamax( n-k, w( k+1, k ), 1 )
            colmax = abs( 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 scopy( 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 scopy( imax-k, a( imax, k ), lda, w( k, k+1 ), 1)
                  CALL scopy( n-imax+1, a( imax, imax ), 1,
     $                        w( imax, k+1 ), 1 )
                  IF( k.GT.1 )
     $               CALL sgemv( 'No transpose', n-k+1, k-1, -one,
     $                           a( k, 1 ), lda, w( imax, 1 ), ldw,
     $                           one, 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 + isamax( imax-k, w( k, k+1 ), 1 )
                     rowmax = abs( w( jmax, k+1 ) )
                  ELSE
                     rowmax = zero
                  END IF
*
                  IF( imax.LT.n ) THEN
                     itemp = imax + isamax( n-imax, w( imax+1, k+1 ), 1)
                     stemp = abs( w( itemp, k+1 ) )
                     IF( stemp.GT.rowmax ) THEN
                        rowmax = stemp
                        jmax = itemp
                     END IF
                  END IF
*
*                 Equivalent to testing for
*                 ABS( W( IMAX, K+1 ) ).GE.ALPHA*ROWMAX
*                 (used to handle NaN and Inf)
*
                  IF( .NOT.( abs( 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 scopy( 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 scopy( 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 scopy( p-k, a( k, k ), 1, a( p, k ), lda )
               CALL scopy( 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 sswap( k, a( k, 1 ), lda, a( p, 1 ), lda )
               CALL sswap( 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 scopy( kp-k-1, a( k+1, kk ), 1, a( kp, k+1 ), lda )
               CALL scopy( 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 sswap( kk, a( kk, 1 ), lda, a( kp, 1 ), lda )
               CALL sswap( 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 scopy( n-k+1, w( k, k ), 1, a( k, k ), 1 )
               IF( k.LT.n ) THEN
                  IF( abs( a( k, k ) ).GE.sfmin ) THEN
                     r1 = one / a( k, k )
                     CALL sscal( n-k, r1, a( k+1, k ), 1 )
                  ELSE IF( a( k, k ).NE.zero ) 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 = one / ( d11*d22-one )
                  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 sgemv( 'No transpose', j+jb-jj, k-1, -one,
     $                     a( jj, 1 ), lda, w( jj, 1 ), ldw, one,
     $                     a( jj, jj ), 1 )
  100       CONTINUE
*
*           Update the rectangular subdiagonal block
*
            IF( j+jb.LE.n )
     $         CALL sgemm( 'No transpose', 'Transpose', n-j-jb+1, jb,
     $                     k-1, -one, a( j+jb, 1 ), lda, w( j, 1 ), ldw,
     $                     one, 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 sswap( j, a( jp2, 1 ), lda, a( jj, 1 ), lda )
            jj = j + 1
            IF( jp1.NE.jj .AND. kstep.EQ.2 )
     $         CALL sswap( 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 SLASYF_ROOK
*

Here is the call graph for this function:

Here is the caller graph for this function: