added standard precision (COMPLEX) routines for rook pivoting \\nalgorithm for symmetric indefinite matrices: slasyf_rook.f ssytri_rook.f ssytrs_rook.f ssytrf_rook.f ssysv_rook.f ssycon_rook.f ssytf2_rook.f

7 files changed, 3599 insertions, 0 deletions

diff --git a/SRC/clasyf_rook.f b/SRC/clasyf_rook.f
new file mode 100644
index 00000000..b46fb987
--- /dev/null
+++ b/SRC/clasyf_rook.f
@@ -0,0 +1,894 @@
+*> \brief \b CLASYF_ROOK
+*
+*  =========== DOCUMENTATION ===========
+*
+* Online html documentation available at 
+*            http://www.netlib.org/lapack/explore-html/ 
+*
+*> \htmlonly
+*> Download CLASYF_ROOK + dependencies 
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clasyf_rook.f"> 
+*> [TGZ]</a> 
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/clasyf_rook.f"> 
+*> [ZIP]</a> 
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/clasyf_rook.f"> 
+*> [TXT]</a>
+*> \endhtmlonly 
+*
+*  Definition:
+*  ===========
+*
+*       SUBROUTINE CLASYF_ROOK( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, INFO )
+* 
+*       .. Scalar Arguments ..
+*       CHARACTER          UPLO
+*       INTEGER            INFO, KB, LDA, LDW, N, NB
+*       ..
+*       .. Array Arguments ..
+*       INTEGER            IPIV( * )
+*       COMPLEX            A( LDA, * ), W( LDW, * )
+*       ..
+*  
+*
+*> \par Purpose:
+*  =============
+*>
+*> \verbatim
+*>
+*> 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').
+*> \endverbatim
+*
+*  Arguments:
+*  ==========
+*
+*> \param[in] UPLO
+*> \verbatim
+*>          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
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*>          N is INTEGER
+*>          The order of the matrix A.  N >= 0.
+*> \endverbatim
+*>
+*> \param[in] NB
+*> \verbatim
+*>          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.
+*> \endverbatim
+*>
+*> \param[out] KB
+*> \verbatim
+*>          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.
+*> \endverbatim
+*>
+*> \param[in,out] A
+*> \verbatim
+*>          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.
+*> \endverbatim
+*>
+*> \param[in] LDA
+*> \verbatim
+*>          LDA is INTEGER
+*>          The leading dimension of the array A.  LDA >= max(1,N).
+*> \endverbatim
+*>
+*> \param[out] IPIV
+*> \verbatim
+*>          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 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.
+*> \endverbatim
+*>
+*> \param[out] W
+*> \verbatim
+*>          W is COMPLEX array, dimension (LDW,NB)
+*> \endverbatim
+*>
+*> \param[in] LDW
+*> \verbatim
+*>          LDW is INTEGER
+*>          The leading dimension of the array W.  LDW >= max(1,N).
+*> \endverbatim
+*>
+*> \param[out] INFO
+*> \verbatim
+*>          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.
+*> \endverbatim
+*
+*  Authors:
+*  ========
+*
+*> \author Univ. of Tennessee 
+*> \author Univ. of California Berkeley 
+*> \author Univ. of Colorado Denver 
+*> \author NAG Ltd. 
+*
+*> \date November 2011
+*
+*> \ingroup complexSYcomputational
+*
+*> \par Contributors:
+*  ==================
+*>
+*> \verbatim
+*>
+*>   November 2011, Igor Kozachenko,
+*>                  Computer Science Division,
+*>                  University of California, Berkeley
+*>
+*>  September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas,
+*>                  School of Mathematics,
+*>                  University of Manchester
+*>
+*> \endverbatim
+*
+*  =====================================================================
+      SUBROUTINE CLASYF_ROOK( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW,
+     $                        INFO )
+*
+*  -- LAPACK computational routine (version 3.3.1) --
+*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*     November 2011
+*
+*     .. 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 (used to handle NaN and Inf)
+*           ABSAKK.GE.ALPHA*COLMAX
+*           
+            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 (used to handle NaN and Inf)
+*                 CABS1( W( IMAX, KW-1 ) ).GE.ALPHA*ROWMAX
+*
+                  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
+*
+                     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 (used to handle NaN and Inf)
+*           ABSAKK.GE.ALPHA*COLMAX
+*           
+            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 (used to handle NaN and Inf)
+*                 CABS1( W( IMAX, K+1 ) ).GE.ALPHA*ROWMAX
+*
+                  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
+*
+                     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
+*
+      END
diff --git a/SRC/csycon_rook.f b/SRC/csycon_rook.f
new file mode 100644
index 00000000..50690dfd
--- /dev/null
+++ b/SRC/csycon_rook.f
@@ -0,0 +1,260 @@
+*> \brief \b CSYCON_ROOK
+*
+*  =========== DOCUMENTATION ===========
+*
+* Online html documentation available at 
+*            http://www.netlib.org/lapack/explore-html/ 
+*
+*> \htmlonly
+*> Download CSYCON_ROOK + dependencies 
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/csycon_rook.f"> 
+*> [TGZ]</a> 
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/csycon_rook.f"> 
+*> [ZIP]</a> 
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/csycon_rook.f"> 
+*> [TXT]</a>
+*> \endhtmlonly 
+*
+*  Definition:
+*  ===========
+*
+*       SUBROUTINE CSYCON_ROOK( UPLO, N, A, LDA, IPIV, ANORM, RCOND,
+*                               WORK, IWORK, INFO )
+* 
+*       .. Scalar Arguments ..
+*       CHARACTER          UPLO
+*       INTEGER            INFO, LDA, N
+*       REAL               ANORM, RCOND
+*       ..
+*       .. Array Arguments ..
+*       INTEGER            IPIV( * ), IWORK( * )
+*       COMPLEX            A( LDA, * ), WORK( * )
+*       ..
+*  
+*
+*> \par Purpose:
+*  =============
+*>
+*> \verbatim
+*>
+*> CSYCON_ROOK estimates the reciprocal of the condition number (in the
+*> 1-norm) of a complex symmetric matrix A using the factorization
+*> A = U*D*U**T or A = L*D*L**T computed by CSYTRF_ROOK.
+*>
+*> An estimate is obtained for norm(inv(A)), and the reciprocal of the
+*> condition number is computed as RCOND = 1 / (ANORM * norm(inv(A))).
+*> \endverbatim
+*
+*  Arguments:
+*  ==========
+*
+*> \param[in] UPLO
+*> \verbatim
+*>          UPLO is CHARACTER*1
+*>          Specifies whether the details of the factorization are stored
+*>          as an upper or lower triangular matrix.
+*>          = 'U':  Upper triangular, form is A = U*D*U**T;
+*>          = 'L':  Lower triangular, form is A = L*D*L**T.
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*>          N is INTEGER
+*>          The order of the matrix A.  N >= 0.
+*> \endverbatim
+*>
+*> \param[in] A
+*> \verbatim
+*>          A is COMPLEX array, dimension (LDA,N)
+*>          The block diagonal matrix D and the multipliers used to
+*>          obtain the factor U or L as computed by CSYTRF_ROOK.
+*> \endverbatim
+*>
+*> \param[in] LDA
+*> \verbatim
+*>          LDA is INTEGER
+*>          The leading dimension of the array A.  LDA >= max(1,N).
+*> \endverbatim
+*>
+*> \param[in] IPIV
+*> \verbatim
+*>          IPIV is INTEGER array, dimension (N)
+*>          Details of the interchanges and the block structure of D
+*>          as determined by CSYTRF_ROOK.
+*> \endverbatim
+*>
+*> \param[in] ANORM
+*> \verbatim
+*>          ANORM is REAL
+*>          The 1-norm of the original matrix A.
+*> \endverbatim
+*>
+*> \param[out] RCOND
+*> \verbatim
+*>          RCOND is REAL
+*>          The reciprocal of the condition number of the matrix A,
+*>          computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is an
+*>          estimate of the 1-norm of inv(A) computed in this routine.
+*> \endverbatim
+*>
+*> \param[out] WORK
+*> \verbatim
+*>          WORK is COMPLEX array, dimension (2*N)
+*> \endverbatim
+*>
+*> \param[out] IWORK
+*> \verbatim
+*>          IWORK is INTEGER array, dimension (N)
+*> \endverbatim
+*>
+*> \param[out] INFO
+*> \verbatim
+*>          INFO is INTEGER
+*>          = 0:  successful exit
+*>          < 0:  if INFO = -i, the i-th argument had an illegal value
+*> \endverbatim
+*
+*  Authors:
+*  ========
+*
+*> \author Univ. of Tennessee 
+*> \author Univ. of California Berkeley 
+*> \author Univ. of Colorado Denver 
+*> \author NAG Ltd. 
+*
+*> \date November 2011
+*
+*> \ingroup complexSYcomputational
+*
+*> \par Contributors:
+*  ==================
+*> \verbatim
+*>
+*>   November 2011, Igor Kozachenko,
+*>                  Computer Science Division,
+*>                  University of California, Berkeley
+*>
+*>  September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas,
+*>                  School of Mathematics,
+*>                  University of Manchester
+*>
+*> \endverbatim
+*
+*  =====================================================================
+      SUBROUTINE CSYCON_ROOK( UPLO, N, A, LDA, IPIV, ANORM, RCOND, WORK,
+     $                   IWORK, INFO )
+*
+*  -- LAPACK computational routine (version 3.4.0) --
+*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*     November 2011
+*
+*     .. Scalar Arguments ..
+      CHARACTER          UPLO
+      INTEGER            INFO, LDA, N
+      REAL               ANORM, RCOND
+*     ..
+*     .. Array Arguments ..
+      INTEGER            IPIV( * ), IWORK( * )
+      COMPLEX            A( LDA, * ), WORK( * )
+*     ..
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      REAL               ONE, ZERO
+      PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
+      COMPLEX            CZERO
+      PARAMETER          ( CZERO = ( 0.0E+0, 0.0E+0 ) )
+*     ..
+*     .. Local Scalars ..
+      LOGICAL            UPPER
+      INTEGER            I, KASE
+      REAL               AINVNM
+*     ..
+*     .. Local Arrays ..
+      INTEGER            ISAVE( 3 )
+*     ..
+*     .. External Functions ..
+      LOGICAL            LSAME
+      EXTERNAL           LSAME
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           CLACN2, CSYTRS_ROOK, XERBLA
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          MAX
+*     ..
+*     .. Executable Statements ..
+*
+*     Test the input parameters.
+*
+      INFO = 0
+      UPPER = LSAME( UPLO, 'U' )
+      IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
+         INFO = -1
+      ELSE IF( N.LT.0 ) THEN
+         INFO = -2
+      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
+         INFO = -4
+      ELSE IF( ANORM.LT.ZERO ) THEN
+         INFO = -6
+      END IF
+      IF( INFO.NE.0 ) THEN
+         CALL XERBLA( 'CSYCON_ROOK', -INFO )
+         RETURN
+      END IF
+*
+*     Quick return if possible
+*
+      RCOND = ZERO
+      IF( N.EQ.0 ) THEN
+         RCOND = ONE
+         RETURN
+      ELSE IF( ANORM.LE.ZERO ) THEN
+         RETURN
+      END IF
+*
+*     Check that the diagonal matrix D is nonsingular.
+*
+      IF( UPPER ) THEN
+*
+*        Upper triangular storage: examine D from bottom to top
+*
+         DO 10 I = N, 1, -1
+            IF( IPIV( I ).GT.0 .AND. A( I, I ).EQ.CZERO )
+     $         RETURN
+   10    CONTINUE
+      ELSE
+*
+*        Lower triangular storage: examine D from top to bottom.
+*
+         DO 20 I = 1, N
+            IF( IPIV( I ).GT.0 .AND. A( I, I ).EQ.CZERO )
+     $         RETURN
+   20    CONTINUE
+      END IF
+*
+*     Estimate the 1-norm of the inverse.
+*
+      KASE = 0
+   30 CONTINUE
+      CALL CLACN2( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE, ISAVE )
+      IF( KASE.NE.0 ) THEN
+*
+*        Multiply by inv(L*D*L**T) or inv(U*D*U**T).
+*
+         CALL CSYTRS_ROOK( UPLO, N, 1, A, LDA, IPIV, WORK, N, INFO )
+         GO TO 30
+      END IF
+*
+*     Compute the estimate of the reciprocal condition number.
+*
+      IF( AINVNM.NE.ZERO )
+     $   RCOND = ( ONE / AINVNM ) / ANORM
+*
+      RETURN
+*
+*     End of CSYCON_ROOK
+*
+      END
diff --git a/SRC/csysv_rook.f b/SRC/csysv_rook.f
new file mode 100644
index 00000000..efac6684
--- /dev/null
+++ b/SRC/csysv_rook.f
@@ -0,0 +1,293 @@
+*> \brief <b> CSYSV_ROOK computes the solution to system of linear equations A * X = B for SY matrices</b>
+*
+*  =========== DOCUMENTATION ===========
+*
+* Online html documentation available at 
+*            http://www.netlib.org/lapack/explore-html/ 
+*
+*> \htmlonly
+*> Download CSYSV_ROOK + dependencies 
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/csysv_rook.f"> 
+*> [TGZ]</a> 
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/csysv_rook.f"> 
+*> [ZIP]</a> 
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/csysv_rook.f"> 
+*> [TXT]</a>
+*> \endhtmlonly 
+*
+*  Definition:
+*  ===========
+*
+*       SUBROUTINE CSYSV_ROOK( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, WORK,
+*                         LWORK, INFO )
+* 
+*       .. Scalar Arguments ..
+*       CHARACTER          UPLO
+*       INTEGER            INFO, LDA, LDB, LWORK, N, NRHS
+*       ..
+*       .. Array Arguments ..
+*       INTEGER            IPIV( * )
+*       COMPLEX            A( LDA, * ), B( LDB, * ), WORK( * )
+*       ..
+*  
+*
+*> \par Purpose:
+*  =============
+*>
+*> \verbatim
+*>
+*> CSYSV_ROOK computes the solution to a complex system of linear
+*> equations
+*>    A * X = B,
+*> where A is an N-by-N symmetric matrix and X and B are N-by-NRHS
+*> matrices.
+*>
+*> The diagonal pivoting method is used to factor A as
+*>    A = U * D * U**T,  if UPLO = 'U', or
+*>    A = L * D * L**T,  if UPLO = 'L',
+*> where U (or L) is a product of permutation and unit upper (lower)
+*> triangular matrices, and D is symmetric and block diagonal with
+*> 1-by-1 and 2-by-2 diagonal blocks.  
+*>
+*> CSYTRF_ROOK is called to compute the factorization of a complex
+*> symmetric matrix A using the bounded Bunch-Kaufman ("rook") diagonal
+*> pivoting method.
+*>
+*> The factored form of A is then used to solve the system 
+*> of equations A * X = B by calling CSYTRS_ROOK.
+*> \endverbatim
+*
+*  Arguments:
+*  ==========
+*
+*> \param[in] UPLO
+*> \verbatim
+*>          UPLO is CHARACTER*1
+*>          = 'U':  Upper triangle of A is stored;
+*>          = 'L':  Lower triangle of A is stored.
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*>          N is INTEGER
+*>          The number of linear equations, i.e., the order of the
+*>          matrix A.  N >= 0.
+*> \endverbatim
+*>
+*> \param[in] NRHS
+*> \verbatim
+*>          NRHS is INTEGER
+*>          The number of right hand sides, i.e., the number of columns
+*>          of the matrix B.  NRHS >= 0.
+*> \endverbatim
+*>
+*> \param[in,out] A
+*> \verbatim
+*>          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, if INFO = 0, the block diagonal matrix D and the
+*>          multipliers used to obtain the factor U or L from the
+*>          factorization A = U*D*U**T or A = L*D*L**T as computed by
+*>          CSYTRF_ROOK.
+*> \endverbatim
+*>
+*> \param[in] LDA
+*> \verbatim
+*>          LDA is INTEGER
+*>          The leading dimension of the array A.  LDA >= max(1,N).
+*> \endverbatim
+*>
+*> \param[out] IPIV
+*> \verbatim
+*>          IPIV is INTEGER array, dimension (N)
+*>          Details of the interchanges and the block structure of D,
+*>          as determined by CSYTRF_ROOK.
+*>
+*>          If UPLO = 'U':
+*>               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':
+*>               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.
+*> \endverbatim
+*>
+*> \param[in,out] B
+*> \verbatim
+*>          B is COMPLEX array, dimension (LDB,NRHS)
+*>          On entry, the N-by-NRHS right hand side matrix B.
+*>          On exit, if INFO = 0, the N-by-NRHS solution matrix X.
+*> \endverbatim
+*>
+*> \param[in] LDB
+*> \verbatim
+*>          LDB is INTEGER
+*>          The leading dimension of the array B.  LDB >= max(1,N).
+*> \endverbatim
+*>
+*> \param[out] WORK
+*> \verbatim
+*>          WORK is COMPLEX array, dimension (MAX(1,LWORK))
+*>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
+*> \endverbatim
+*>
+*> \param[in] LWORK
+*> \verbatim
+*>          LWORK is INTEGER
+*>          The length of WORK.  LWORK >= 1, and for best performance
+*>          LWORK >= max(1,N*NB), where NB is the optimal blocksize for
+*>          CSYTRF_ROOK.
+*>          
+*>          TRS will be done with Level 2 BLAS
+*>
+*>          If LWORK = -1, then a workspace query is assumed; the routine
+*>          only calculates the optimal size of the WORK array, returns
+*>          this value as the first entry of the WORK array, and no error
+*>          message related to LWORK is issued by XERBLA.
+*> \endverbatim
+*>
+*> \param[out] INFO
+*> \verbatim
+*>          INFO is INTEGER
+*>          = 0: successful exit
+*>          < 0: if INFO = -i, the i-th argument had an illegal value
+*>          > 0: if INFO = i, D(i,i) is exactly zero.  The factorization
+*>               has been completed, but the block diagonal matrix D is
+*>               exactly singular, so the solution could not be computed.
+*> \endverbatim
+*
+*  Authors:
+*  ========
+*
+*> \author Univ. of Tennessee 
+*> \author Univ. of California Berkeley 
+*> \author Univ. of Colorado Denver 
+*> \author NAG Ltd. 
+*
+*> \date November 2011
+*
+*> \ingroup complexSYsolve
+*
+*> \par Contributors:
+*  ==================
+*>
+*> \verbatim
+*>
+*>   November 2011, Igor Kozachenko,
+*>                  Computer Science Division,
+*>                  University of California, Berkeley
+*>
+*>  September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas,
+*>                  School of Mathematics,
+*>                  University of Manchester
+*>
+*> \endverbatim
+*
+*  =====================================================================
+      SUBROUTINE CSYSV_ROOK( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, WORK,
+     $                  LWORK, INFO )
+*
+*  -- LAPACK driver routine (version 3.4.0) --
+*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*     November 2011
+*
+*     .. Scalar Arguments ..
+      CHARACTER          UPLO
+      INTEGER            INFO, LDA, LDB, LWORK, N, NRHS
+*     ..
+*     .. Array Arguments ..
+      INTEGER            IPIV( * )
+      COMPLEX            A( LDA, * ), B( LDB, * ), WORK( * )
+*     ..
+*
+*  =====================================================================
+*
+*     .. Local Scalars ..
+      LOGICAL            LQUERY
+      INTEGER            LWKOPT
+*     ..
+*     .. External Functions ..
+      LOGICAL            LSAME
+      EXTERNAL           LSAME
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           XERBLA, CSYTRF_ROOK, CSYTRS_ROOK
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          MAX
+*     ..
+*     .. Executable Statements ..
+*
+*     Test the input parameters.
+*
+      INFO = 0
+      LQUERY = ( LWORK.EQ.-1 )
+      IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
+         INFO = -1
+      ELSE IF( N.LT.0 ) THEN
+         INFO = -2
+      ELSE IF( NRHS.LT.0 ) THEN
+         INFO = -3
+      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
+         INFO = -5
+      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
+         INFO = -8
+      ELSE IF( LWORK.LT.1 .AND. .NOT.LQUERY ) THEN
+         INFO = -10
+      END IF
+*
+      IF( INFO.EQ.0 ) THEN
+         IF( N.EQ.0 ) THEN
+            LWKOPT = 1
+         ELSE
+            CALL DSYTRF( UPLO, N, A, LDA, IPIV, WORK, -1, INFO )
+            LWKOPT = WORK(1)
+         END IF
+         WORK( 1 ) = LWKOPT
+      END IF
+*
+      IF( INFO.NE.0 ) THEN
+         CALL XERBLA( 'CSYSV_ROOK ', -INFO )
+         RETURN
+      ELSE IF( LQUERY ) THEN
+         RETURN
+      END IF
+*
+*     Compute the factorization A = U*D*U**T or A = L*D*L**T.
+*
+      CALL CSYTRF_ROOK( UPLO, N, A, LDA, IPIV, WORK, LWORK, INFO )
+      IF( INFO.EQ.0 ) THEN
+*
+*        Solve the system A*X = B, overwriting B with X.
+*
+*        Solve with TRS_ROOK ( Use Level 2 BLAS)
+*
+         CALL CSYTRS_ROOK( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, INFO )
+*
+      END IF
+*
+      WORK( 1 ) = LWKOPT
+*
+      RETURN
+*
+*     End of CSYSV_ROOK
+*
+      END
diff --git a/SRC/csytf2_rook.f b/SRC/csytf2_rook.f
new file mode 100644
index 00000000..50b9c952
--- /dev/null
+++ b/SRC/csytf2_rook.f
@@ -0,0 +1,824 @@
+*> \brief \b CSYTF2_ROOK
+*
+*  =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+*            http://www.netlib.org/lapack/explore-html/
+*
+*> \htmlonly
+*> Download CSYTF2_ROOK + dependencies
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/csytf2_rook.f">
+*> [TGZ]</a>
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/csytf2_rook.f">
+*> [ZIP]</a>
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/csytf2_rook.f">
+*> [TXT]</a>
+*> \endhtmlonly
+*
+*  Definition:
+*  ===========
+*
+*       SUBROUTINE CSYTF2_ROOK( UPLO, N, A, LDA, IPIV, INFO )
+*
+*       .. Scalar Arguments ..
+*       CHARACTER          UPLO
+*       INTEGER            INFO, LDA, N
+*       ..
+*       .. Array Arguments ..
+*       INTEGER            IPIV( * )
+*       COMPLEX            A( LDA, * )
+*       ..
+*
+*
+*> \par Purpose:
+*  =============
+*>
+*> \verbatim
+*>
+*> CSYTF2_ROOK computes the factorization of a complex symmetric matrix A
+*> using the bounded Bunch-Kaufman ("rook") diagonal pivoting method:
+*>
+*>    A = U*D*U**T  or  A = L*D*L**T
+*>
+*> where U (or L) is a product of permutation and unit upper (lower)
+*> triangular matrices, U**T is the transpose of U, and D is symmetric and
+*> block diagonal with 1-by-1 and 2-by-2 diagonal blocks.
+*>
+*> This is the unblocked version of the algorithm, calling Level 2 BLAS.
+*> \endverbatim
+*
+*  Arguments:
+*  ==========
+*
+*> \param[in] UPLO
+*> \verbatim
+*>          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
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*>          N is INTEGER
+*>          The order of the matrix A.  N >= 0.
+*> \endverbatim
+*>
+*> \param[in,out] A
+*> \verbatim
+*>          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, the block diagonal matrix D and the multipliers used
+*>          to obtain the factor U or L (see below for further details).
+*> \endverbatim
+*>
+*> \param[in] LDA
+*> \verbatim
+*>          LDA is INTEGER
+*>          The leading dimension of the array A.  LDA >= max(1,N).
+*> \endverbatim
+*>
+*> \param[out] IPIV
+*> \verbatim
+*>          IPIV is INTEGER array, dimension (N)
+*>          Details of the interchanges and the block structure of D.
+*>
+*>          If UPLO = 'U':
+*>               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':
+*>               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.
+*> \endverbatim
+*>
+*> \param[out] INFO
+*> \verbatim
+*>          INFO is INTEGER
+*>          = 0: successful exit
+*>          < 0: if INFO = -k, the k-th argument had an illegal value
+*>          > 0: if INFO = k, D(k,k) is exactly zero.  The factorization
+*>               has been completed, but the block diagonal matrix D is
+*>               exactly singular, and division by zero will occur if it
+*>               is used to solve a system of equations.
+*> \endverbatim
+*
+*  Authors:
+*  ========
+*
+*> \author Univ. of Tennessee
+*> \author Univ. of California Berkeley
+*> \author Univ. of Colorado Denver
+*> \author NAG Ltd.
+*
+*> \date November 2011
+*
+*> \ingroup complexSYcomputational
+*
+*> \par Further Details:
+*  =====================
+*>
+*> \verbatim
+*>
+*>  If UPLO = 'U', then A = U*D*U**T, where
+*>     U = P(n)*U(n)* ... *P(k)U(k)* ...,
+*>  i.e., U is a product of terms P(k)*U(k), where k decreases from n to
+*>  1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
+*>  and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
+*>  defined by IPIV(k), and U(k) is a unit upper triangular matrix, such
+*>  that if the diagonal block D(k) is of order s (s = 1 or 2), then
+*>
+*>             (   I    v    0   )   k-s
+*>     U(k) =  (   0    I    0   )   s
+*>             (   0    0    I   )   n-k
+*>                k-s   s   n-k
+*>
+*>  If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k).
+*>  If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k),
+*>  and A(k,k), and v overwrites A(1:k-2,k-1:k).
+*>
+*>  If UPLO = 'L', then A = L*D*L**T, where
+*>     L = P(1)*L(1)* ... *P(k)*L(k)* ...,
+*>  i.e., L is a product of terms P(k)*L(k), where k increases from 1 to
+*>  n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
+*>  and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
+*>  defined by IPIV(k), and L(k) is a unit lower triangular matrix, such
+*>  that if the diagonal block D(k) is of order s (s = 1 or 2), then
+*>
+*>             (   I    0     0   )  k-1
+*>     L(k) =  (   0    I     0   )  s
+*>             (   0    v     I   )  n-k-s+1
+*>                k-1   s  n-k-s+1
+*>
+*>  If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k).
+*>  If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k),
+*>  and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1).
+*> \endverbatim
+*
+*> \par Contributors:
+*  ==================
+*>
+*> \verbatim
+*>
+*>   November 2011, Igor Kozachenko,
+*>                  Computer Science Division,
+*>                  University of California, Berkeley
+*>
+*>  September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas,
+*>                  School of Mathematics,
+*>                  University of Manchester
+*>
+*>  01-01-96 - Based on modifications by
+*>    J. Lewis, Boeing Computer Services Company
+*>    A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA
+*>  1-96 - Based on modifications by J. Lewis, Boeing Computer Services
+*>         Company
+*> \endverbatim
+*
+*  =====================================================================
+      SUBROUTINE CSYTF2_ROOK( UPLO, N, A, LDA, IPIV, INFO )
+*
+*  -- LAPACK computational routine (version 3.3.1) --
+*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*     November 2011
+*
+*     .. Scalar Arguments ..
+      CHARACTER          UPLO
+      INTEGER            INFO, LDA, N
+*     ..
+*     .. Array Arguments ..
+      INTEGER            IPIV( * )
+      COMPLEX            A( LDA, * )
+*     ..
+*
+*  =====================================================================
+*
+*     .. 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
+      PARAMETER          ( CONE = ( 1.0E+0, 0.0E+0 ) )
+*     ..
+*     .. Local Scalars ..
+      LOGICAL            UPPER, DONE
+      INTEGER            I, IMAX, J, JMAX, ITEMP, K, KK, KP, KSTEP,
+     $                   P, II
+      REAL               ABSAKK, ALPHA, COLMAX, ROWMAX, STEMP, SFMIN
+      COMPLEX            D11, D12, D21, D22,
+     $                   T, WK, WKM1, WKP1, Z
+*     ..
+*     .. External Functions ..
+      LOGICAL            LSAME
+      INTEGER            ICAMAX
+      REAL               SLAMCH
+      EXTERNAL           LSAME, ICAMAX, SLAMCH
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           CSCAL, CSWAP, CSYR, XERBLA
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          ABS, MAX, SQRT, AIMAG, REAL
+*     ..
+*     .. Statement Functions ..
+      REAL               CABS1
+*     ..
+*     .. Statement Function definitions ..
+      CABS1( Z ) = ABS( REAL( Z ) ) + ABS( AIMAG( Z ) )
+*     ..
+*     .. Executable Statements ..
+*
+*     Test the input parameters.
+*
+      INFO = 0
+      UPPER = LSAME( UPLO, 'U' )
+      IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
+         INFO = -1
+      ELSE IF( N.LT.0 ) THEN
+         INFO = -2
+      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
+         INFO = -4
+      END IF
+      IF( INFO.NE.0 ) THEN
+         CALL XERBLA( 'CSYTF2_ROOK', -INFO )
+         RETURN
+      END IF
+*
+*     Initialize ALPHA for use in choosing pivot block size.
+*
+      ALPHA = ( ONE+SQRT( SEVTEN ) ) / EIGHT
+*
+*     Compute machine safe minimum
+*
+      SFMIN = SLAMCH( 'S' )
+*
+      IF( UPPER ) THEN
+*
+*        Factorize A as U*D*U**T using the upper triangle of A
+*
+*        K is the main loop index, decreasing from N to 1 in steps of
+*        1 or 2
+*
+         K = N
+   10    CONTINUE
+*
+*        If K < 1, exit from loop
+*
+         IF( K.LT.1 )
+     $      GO TO 70
+         KSTEP = 1
+         P = K
+*
+*        Determine rows and columns to be interchanged and whether
+*        a 1-by-1 or 2-by-2 pivot block will be used
+*
+         ABSAKK = CABS1( A( 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.GT.1 ) THEN
+            IMAX = ICAMAX( K-1, A( 1, K ), 1 )
+            COLMAX = CABS1( A( IMAX, K ) )
+         ELSE
+            COLMAX = ZERO
+         END IF
+*
+         IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) ) THEN
+*
+*           Column K is zero: set INFO and continue
+*
+            IF( INFO.EQ.0 )
+     $         INFO = K
+            KP = K
+         ELSE
+*
+*           Test for interchange
+*
+*           Equivalent to testing for (used to handle NaN and Inf)
+*           ABSAKK.GE.ALPHA*COLMAX
+*
+            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
+*
+*                 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, A( IMAX, IMAX+1 ),
+     $                                    LDA )
+                     ROWMAX = CABS1( A( IMAX, JMAX ) )
+                  ELSE
+                     ROWMAX = ZERO
+                  END IF
+*
+                  IF( IMAX.GT.1 ) THEN
+                     ITEMP = ICAMAX( IMAX-1, A( 1, IMAX ), 1 )
+                     STEMP = CABS1( A( ITEMP, IMAX ) )
+                     IF( STEMP.GT.ROWMAX ) THEN
+                        ROWMAX = STEMP
+                        JMAX = ITEMP
+                     END IF
+                  END IF
+*
+*                 Equivalent to testing for (used to handle NaN and Inf)
+*                 CABS1( A( IMAX, IMAX ) ).GE.ALPHA*ROWMAX
+*
+                  IF( .NOT.( CABS1(A( IMAX, IMAX )).LT.ALPHA*ROWMAX ) )
+     $            THEN
+*
+*                    interchange rows and columns K and IMAX,
+*                    use 1-by-1 pivot block
+*
+                     KP = IMAX
+                     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 variables and repeat
+*
+                     P = IMAX
+                     COLMAX = ROWMAX
+                     IMAX = JMAX
+                  END IF
+*
+*                 End pivot search loop body
+*
+               IF( .NOT. DONE ) GOTO 12
+*
+            END IF
+*
+*           Swap TWO rows and TWO columns
+*
+*           First swap
+*
+            IF( ( KSTEP.EQ.2 ) .AND. ( P.NE.K ) ) THEN
+*
+*              Interchange rows and column K and P in the leading
+*              submatrix A(1:k,1:k) if we have a 2-by-2 pivot
+*
+               IF( P.GT.1 )
+     $            CALL CSWAP( P-1, A( 1, K ), 1, A( 1, P ), 1 )
+               IF( P.LT.(K-1) )
+     $            CALL CSWAP( K-P-1, A( P+1, K ), 1, A( P, P+1 ),
+     $                     LDA )
+               T = A( K, K )
+               A( K, K ) = A( P, P )
+               A( P, P ) = T
+            END IF
+*
+*           Second swap
+*
+            KK = K - KSTEP + 1
+            IF( KP.NE.KK ) THEN
+*
+*              Interchange rows and columns KK and KP in the leading
+*              submatrix A(1:k,1:k)
+*
+               IF( KP.GT.1 )
+     $            CALL CSWAP( KP-1, A( 1, KK ), 1, A( 1, KP ), 1 )
+               IF( ( KK.GT.1 ) .AND. ( KP.LT.(KK-1) ) )
+     $            CALL CSWAP( KK-KP-1, A( KP+1, KK ), 1, A( KP, KP+1 ),
+     $                     LDA )
+               T = A( KK, KK )
+               A( KK, KK ) = A( KP, KP )
+               A( KP, KP ) = T
+               IF( KSTEP.EQ.2 ) THEN
+                  T = A( K-1, K )
+                  A( K-1, K ) = A( KP, K )
+                  A( KP, K ) = T
+               END IF
+            END IF
+*
+*           Update the leading submatrix
+*
+            IF( KSTEP.EQ.1 ) THEN
+*
+*              1-by-1 pivot block D(k): column k now holds
+*
+*              W(k) = U(k)*D(k)
+*
+*              where U(k) is the k-th column of U
+*
+               IF( K.GT.1 ) THEN
+*
+*                 Perform a rank-1 update of A(1:k-1,1:k-1) and
+*                 store U(k) in column k
+*
+                  IF( CABS1( A( K, K ) ).GE.SFMIN ) THEN
+*
+*                    Perform a rank-1 update of A(1:k-1,1:k-1) as
+*                    A := A - U(k)*D(k)*U(k)**T 
+*                       = A - W(k)*1/D(k)*W(k)**T
+*
+                     D11 = CONE / A( K, K )
+                     CALL CSYR( UPLO, K-1, -D11, A( 1, K ), 1, A, LDA )
+*
+*                    Store U(k) in column k
+*
+                     CALL CSCAL( K-1, D11, A( 1, K ), 1 )
+                  ELSE
+*
+*                    Store L(k) in column K
+*
+                     D11 = A( K, K )
+                     DO 16 II = 1, K - 1
+                        A( II, K ) = A( II, K ) / D11
+   16                CONTINUE
+*
+*                    Perform a rank-1 update of A(k+1:n,k+1:n) as
+*                    A := A - U(k)*D(k)*U(k)**T
+*                       = A - W(k)*(1/D(k))*W(k)**T
+*                       = A - (W(k)/D(k))*(D(k))*(W(k)/D(K))**T
+*
+                     CALL CSYR( UPLO, K-1, -D11, A( 1, K ), 1, A, LDA )
+                  END IF
+               END IF
+*
+            ELSE
+*
+*              2-by-2 pivot block D(k): columns k and k-1 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
+*
+*              Perform a rank-2 update of A(1:k-2,1:k-2) as
+*
+*              A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )**T
+*                 = A - ( ( A(k-1)A(k) )*inv(D(k)) ) * ( A(k-1)A(k) )**T
+*
+*              and store L(k) and L(k+1) in columns k and k+1
+*
+               IF( K.GT.2 ) THEN
+*
+                  D12 = A( K-1, K )
+                  D22 = A( K-1, K-1 ) / D12
+                  D11 = A( K, K ) / D12
+                  T = CONE / ( D11*D22-CONE )
+*
+                  DO 30 J = K - 2, 1, -1
+*
+                     WKM1 = T*( D11*A( J, K-1 )-A( J, K ) )
+                     WK = T*( D22*A( J, K )-A( J, K-1 ) )
+*
+                     DO 20 I = J, 1, -1
+                        A( I, J ) = A( I, J ) - (A( I, K ) / D12 )*WK -
+     $                              ( A( I, K-1 ) / D12 )*WKM1
+   20                CONTINUE
+*
+*                    Store U(k) and U(k-1) in cols k and k-1 for row J
+*
+                     A( J, K ) = WK / D12
+                     A( J, K-1 ) = WKM1 / D12
+*
+   30             CONTINUE
+*
+               END IF
+*
+            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
+*
+      ELSE
+*
+*        Factorize A as L*D*L**T using the lower triangle of A
+*
+*        K is the main loop index, increasing from 1 to N in steps of
+*        1 or 2
+*
+         K = 1
+   40    CONTINUE
+*
+*        If K > N, exit from loop
+*
+         IF( K.GT.N )
+     $      GO TO 70
+         KSTEP = 1
+         P = K
+*
+*        Determine rows and columns to be interchanged and whether
+*        a 1-by-1 or 2-by-2 pivot block will be used
+*
+         ABSAKK = CABS1( A( 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, A( K+1, K ), 1 )
+            COLMAX = CABS1( A( IMAX, K ) )
+         ELSE
+            COLMAX = ZERO
+         END IF
+*
+         IF( ( MAX( ABSAKK, COLMAX ).EQ.ZERO ) ) THEN
+*
+*           Column K is zero: set INFO and continue
+*
+            IF( INFO.EQ.0 )
+     $         INFO = K
+            KP = K
+         ELSE
+*
+*           Test for interchange
+*
+*           Equivalent to testing for (used to handle NaN and Inf)
+*           ABSAKK.GE.ALPHA*COLMAX
+*
+            IF( .NOT.( ABSAKK.LT.ALPHA*COLMAX ) ) THEN
+*
+*              no interchange, use 1-by-1 pivot block
+*
+               KP = K
+            ELSE
+*
+               DONE = .FALSE.
+*
+*              Loop until pivot found
+*
+   42          CONTINUE
+*
+*                 Begin pivot search loop body
+*
+*                 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, A( IMAX, K ), LDA )
+                     ROWMAX = CABS1( A( IMAX, JMAX ) )
+                  ELSE
+                     ROWMAX = ZERO
+                  END IF
+*
+                  IF( IMAX.LT.N ) THEN
+                     ITEMP = IMAX + ICAMAX( N-IMAX, A( IMAX+1, IMAX ),
+     $                                     1 )
+                     STEMP = CABS1( A( ITEMP, IMAX ) )
+                     IF( STEMP.GT.ROWMAX ) THEN
+                        ROWMAX = STEMP
+                        JMAX = ITEMP
+                     END IF
+                  END IF
+*
+*                 Equivalent to testing for (used to handle NaN and Inf) 
+*                 CABS1( A( IMAX, IMAX ) ).GE.ALPHA*ROWMAX
+*
+                  IF( .NOT.( CABS1(A( IMAX, IMAX )).LT.ALPHA*ROWMAX ) )
+     $            THEN
+*
+*                    interchange rows and columns K and IMAX,
+*                    use 1-by-1 pivot block
+*
+                     KP = IMAX
+                     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 variables and repeat
+*
+                     P = IMAX
+                     COLMAX = ROWMAX
+                     IMAX = JMAX
+                  END IF
+*
+*                 End pivot search loop body
+*
+               IF( .NOT. DONE ) GOTO 42
+*
+            END IF
+*
+*           Swap TWO rows and TWO columns
+*
+*           First swap
+*
+            IF( ( KSTEP.EQ.2 ) .AND. ( P.NE.K ) ) THEN
+*
+*              Interchange rows and column K and P in the trailing
+*              submatrix A(k:n,k:n) if we have a 2-by-2 pivot
+*
+               IF( P.LT.N )
+     $            CALL CSWAP( N-P, A( P+1, K ), 1, A( P+1, P ), 1 )
+               IF( P.GT.(K+1) )
+     $            CALL CSWAP( P-K-1, A( K+1, K ), 1, A( P, K+1 ), LDA )
+               T = A( K, K )
+               A( K, K ) = A( P, P )
+               A( P, P ) = T
+            END IF
+*
+*           Second swap
+*
+            KK = K + KSTEP - 1
+            IF( KP.NE.KK ) THEN
+*
+*              Interchange rows and columns KK and KP in the trailing
+*              submatrix A(k:n,k:n)
+*
+               IF( KP.LT.N )
+     $            CALL CSWAP( N-KP, A( KP+1, KK ), 1, A( KP+1, KP ), 1 )
+               IF( ( KK.LT.N ) .AND. ( KP.GT.(KK+1) ) )
+     $            CALL CSWAP( KP-KK-1, A( KK+1, KK ), 1, A( KP, KK+1 ),
+     $                     LDA )
+               T = A( KK, KK )
+               A( KK, KK ) = A( KP, KP )
+               A( KP, KP ) = T
+               IF( KSTEP.EQ.2 ) THEN
+                  T = A( K+1, K )
+                  A( K+1, K ) = A( KP, K )
+                  A( KP, K ) = T
+               END IF
+            END IF
+*
+*           Update the trailing submatrix
+*
+            IF( KSTEP.EQ.1 ) THEN
+*
+*              1-by-1 pivot block D(k): column k now holds
+*
+*              W(k) = L(k)*D(k)
+*
+*              where L(k) is the k-th column of L
+*
+               IF( K.LT.N ) THEN
+*
+*              Perform a rank-1 update of A(k+1:n,k+1:n) and
+*              store L(k) in column k
+*
+                  IF( CABS1( A( K, K ) ).GE.SFMIN ) THEN
+*
+*                    Perform a rank-1 update of A(k+1:n,k+1:n) as
+*                    A := A - L(k)*D(k)*L(k)**T 
+*                       = A - W(k)*(1/D(k))*W(k)**T
+*
+                     D11 = CONE / A( K, K )
+                     CALL CSYR( UPLO, N-K, -D11, A( K+1, K ), 1,
+     $                          A( K+1, K+1 ), LDA )
+*
+*                    Store L(k) in column k
+*
+                     CALL CSCAL( N-K, D11, A( K+1, K ), 1 )
+                  ELSE
+*
+*                    Store L(k) in column k
+*
+                     D11 = A( K, K )
+                     DO 46 II = K + 1, N
+                        A( II, K ) = A( II, K ) / D11
+   46                CONTINUE
+*
+*                    Perform a rank-1 update of A(k+1:n,k+1:n) as
+*                    A := A - L(k)*D(k)*L(k)**T
+*                       = A - W(k)*(1/D(k))*W(k)**T
+*                       = A - (W(k)/D(k))*(D(k))*(W(k)/D(K))**T
+*
+                     CALL CSYR( UPLO, N-K, -D11, A( K+1, K ), 1,
+     $                          A( K+1, K+1 ), LDA )
+                  END IF
+               END IF
+*
+            ELSE
+*
+*              2-by-2 pivot block D(k): columns k and k+1 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
+*
+*
+*              Perform a rank-2 update of A(k+2:n,k+2:n) as
+*
+*              A := A - ( L(k) L(k+1) ) * D(k) * ( L(k) L(k+1) )**T
+*                 = A - ( ( A(k)A(k+1) )*inv(D(k) ) * ( A(k)A(k+1) )**T
+*
+*              and store L(k) and L(k+1) in columns k and k+1
+*
+               IF( K.LT.N-1 ) THEN
+*
+                  D21 = A( K+1, K )
+                  D11 = A( K+1, K+1 ) / D21
+                  D22 = A( K, K ) / D21
+                  T = CONE / ( D11*D22-CONE )
+*
+                  DO 60 J = K + 2, N
+*
+*                    Compute  D21 * ( W(k)W(k+1) ) * inv(D(k)) for row J
+*
+                     WK = T*( D11*A( J, K )-A( J, K+1 ) )
+                     WKP1 = T*( D22*A( J, K+1 )-A( J, K ) )
+*
+*                    Perform a rank-2 update of A(k+2:n,k+2:n)
+*
+                     DO 50 I = J, N
+                        A( I, J ) = A( I, J ) - ( A( I, K ) / D21 )*WK -
+     $                              ( A( I, K+1 ) / D21 )*WKP1
+   50                CONTINUE
+*
+*                    Store L(k) and L(k+1) in cols k and k+1 for row J
+*
+                     A( J, K ) = WK / D21
+                     A( J, K+1 ) = WKP1 / D21
+*
+   60             CONTINUE
+*
+               END IF
+*
+            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 40
+*
+      END IF
+*
+   70 CONTINUE
+*
+      RETURN
+*
+*     End of CSYTF2_ROOK
+*
+      END
diff --git a/SRC/csytrf_rook.f b/SRC/csytrf_rook.f
new file mode 100644
index 00000000..d4169d74
--- /dev/null
+++ b/SRC/csytrf_rook.f
@@ -0,0 +1,393 @@
+*> \brief \b CSYTRF_ROOK
+*
+*  =========== DOCUMENTATION ===========
+*
+* Online html documentation available at 
+*            http://www.netlib.org/lapack/explore-html/ 
+*
+*> \htmlonly
+*> Download CSYTRF_ROOK + dependencies 
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/csytrf_rook.f"> 
+*> [TGZ]</a> 
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/csytrf_rook.f"> 
+*> [ZIP]</a> 
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/csytrf_rook.f"> 
+*> [TXT]</a>
+*> \endhtmlonly 
+*
+*  Definition:
+*  ===========
+*
+*       SUBROUTINE CSYTRF_ROOK( UPLO, N, A, LDA, IPIV, WORK, LWORK, INFO )
+* 
+*       .. Scalar Arguments ..
+*       CHARACTER          UPLO
+*       INTEGER            INFO, LDA, LWORK, N
+*       ..
+*       .. Array Arguments ..
+*       INTEGER            IPIV( * )
+*       COMPLEX            A( LDA, * ), WORK( * )
+*       ..
+*  
+*
+*> \par Purpose:
+*  =============
+*>
+*> \verbatim
+*>
+*> CSYTRF_ROOK computes the factorization of a complex symmetric matrix A
+*> using the bounded Bunch-Kaufman ("rook") diagonal pivoting method.
+*> The form of the factorization is
+*>
+*>    A = U*D*U**T  or  A = L*D*L**T
+*>
+*> where U (or L) is a product of permutation and unit upper (lower)
+*> triangular matrices, and D is symmetric and block diagonal with
+*> 1-by-1 and 2-by-2 diagonal blocks.
+*>
+*> This is the blocked version of the algorithm, calling Level 3 BLAS.
+*> \endverbatim
+*
+*  Arguments:
+*  ==========
+*
+*> \param[in] UPLO
+*> \verbatim
+*>          UPLO is CHARACTER*1
+*>          = 'U':  Upper triangle of A is stored;
+*>          = 'L':  Lower triangle of A is stored.
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*>          N is INTEGER
+*>          The order of the matrix A.  N >= 0.
+*> \endverbatim
+*>
+*> \param[in,out] A
+*> \verbatim
+*>          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, the block diagonal matrix D and the multipliers used
+*>          to obtain the factor U or L (see below for further details).
+*> \endverbatim
+*>
+*> \param[in] LDA
+*> \verbatim
+*>          LDA is INTEGER
+*>          The leading dimension of the array A.  LDA >= max(1,N).
+*> \endverbatim
+*>
+*> \param[out] IPIV
+*> \verbatim
+*>          IPIV is INTEGER array, dimension (N)
+*>          Details of the interchanges and the block structure of D.
+*>
+*>          If UPLO = 'U':
+*>               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':
+*>               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.
+*> \endverbatim
+*>
+*> \param[out] WORK
+*> \verbatim
+*>          WORK is COMPLEX array, dimension (MAX(1,LWORK)).
+*>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
+*> \endverbatim
+*>
+*> \param[in] LWORK
+*> \verbatim
+*>          LWORK is INTEGER
+*>          The length of WORK.  LWORK >=1.  For best performance
+*>          LWORK >= N*NB, where NB is the block size returned by ILAENV.
+*>
+*>          If LWORK = -1, then a workspace query is assumed; the routine
+*>          only calculates the optimal size of the WORK array, returns
+*>          this value as the first entry of the WORK array, and no error
+*>          message related to LWORK is issued by XERBLA.
+*> \endverbatim
+*>
+*> \param[out] INFO
+*> \verbatim
+*>          INFO is INTEGER
+*>          = 0:  successful exit
+*>          < 0:  if INFO = -i, the i-th argument had an illegal value
+*>          > 0:  if INFO = i, D(i,i) is exactly zero.  The factorization
+*>                has been completed, but the block diagonal matrix D is
+*>                exactly singular, and division by zero will occur if it
+*>                is used to solve a system of equations.
+*> \endverbatim
+*
+*  Authors:
+*  ========
+*
+*> \author Univ. of Tennessee 
+*> \author Univ. of California Berkeley 
+*> \author Univ. of Colorado Denver 
+*> \author NAG Ltd. 
+*
+*> \date November 2011
+*
+*> \ingroup complexSYcomputational
+*
+*> \par Further Details:
+*  =====================
+*>
+*> \verbatim
+*>
+*>  If UPLO = 'U', then A = U*D*U**T, where
+*>     U = P(n)*U(n)* ... *P(k)U(k)* ...,
+*>  i.e., U is a product of terms P(k)*U(k), where k decreases from n to
+*>  1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
+*>  and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
+*>  defined by IPIV(k), and U(k) is a unit upper triangular matrix, such
+*>  that if the diagonal block D(k) is of order s (s = 1 or 2), then
+*>
+*>             (   I    v    0   )   k-s
+*>     U(k) =  (   0    I    0   )   s
+*>             (   0    0    I   )   n-k
+*>                k-s   s   n-k
+*>
+*>  If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k).
+*>  If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k),
+*>  and A(k,k), and v overwrites A(1:k-2,k-1:k).
+*>
+*>  If UPLO = 'L', then A = L*D*L**T, where
+*>     L = P(1)*L(1)* ... *P(k)*L(k)* ...,
+*>  i.e., L is a product of terms P(k)*L(k), where k increases from 1 to
+*>  n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
+*>  and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
+*>  defined by IPIV(k), and L(k) is a unit lower triangular matrix, such
+*>  that if the diagonal block D(k) is of order s (s = 1 or 2), then
+*>
+*>             (   I    0     0   )  k-1
+*>     L(k) =  (   0    I     0   )  s
+*>             (   0    v     I   )  n-k-s+1
+*>                k-1   s  n-k-s+1
+*>
+*>  If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k).
+*>  If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k),
+*>  and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1).
+*> \endverbatim
+*
+*> \par Contributors:
+*  ==================
+*>
+*> \verbatim
+*>
+*>   November 2011, Igor Kozachenko,
+*>                  Computer Science Division,
+*>                  University of California, Berkeley
+*>
+*>  September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas,
+*>                  School of Mathematics,
+*>                  University of Manchester
+*>
+*> \endverbatim
+*
+*  =====================================================================
+      SUBROUTINE CSYTRF_ROOK( UPLO, N, A, LDA, IPIV, WORK, LWORK, INFO )
+*
+*  -- LAPACK computational routine (version 3.4.0) --
+*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*     November 2011
+*
+*     .. Scalar Arguments ..
+      CHARACTER          UPLO
+      INTEGER            INFO, LDA, LWORK, N
+*     ..
+*     .. Array Arguments ..
+      INTEGER            IPIV( * )
+      COMPLEX            A( LDA, * ), WORK( * )
+*     ..
+*
+*  =====================================================================
+*
+*     .. Local Scalars ..
+      LOGICAL            LQUERY, UPPER
+      INTEGER            IINFO, IWS, J, K, KB, LDWORK, LWKOPT, NB, NBMIN
+*     ..
+*     .. External Functions ..
+      LOGICAL            LSAME
+      INTEGER            ILAENV
+      EXTERNAL           LSAME, ILAENV
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           CLASYF_ROOK, CSYTF2_ROOK, XERBLA
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          MAX
+*     ..
+*     .. Executable Statements ..
+*
+*     Test the input parameters.
+*
+      INFO = 0
+      UPPER = LSAME( UPLO, 'U' )
+      LQUERY = ( LWORK.EQ.-1 )
+      IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
+         INFO = -1
+      ELSE IF( N.LT.0 ) THEN
+         INFO = -2
+      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
+         INFO = -4
+      ELSE IF( LWORK.LT.1 .AND. .NOT.LQUERY ) THEN
+         INFO = -7
+      END IF
+*
+      IF( INFO.EQ.0 ) THEN
+*
+*        Determine the block size
+*
+         NB = ILAENV( 1, 'CSYTRF_ROOK', UPLO, N, -1, -1, -1 )
+         LWKOPT = N*NB
+         WORK( 1 ) = LWKOPT
+      END IF
+*
+      IF( INFO.NE.0 ) THEN
+         CALL XERBLA( 'CSYTRF_ROOK', -INFO )
+         RETURN
+      ELSE IF( LQUERY ) THEN
+         RETURN
+      END IF
+*
+      NBMIN = 2
+      LDWORK = N
+      IF( NB.GT.1 .AND. NB.LT.N ) THEN
+         IWS = LDWORK*NB
+         IF( LWORK.LT.IWS ) THEN
+            NB = MAX( LWORK / LDWORK, 1 )
+            NBMIN = MAX( 2, ILAENV( 2, 'CSYTRF_ROOK',
+     $                              UPLO, N, -1, -1, -1 ) )
+         END IF
+      ELSE
+         IWS = 1
+      END IF
+      IF( NB.LT.NBMIN )
+     $   NB = N
+*
+      IF( UPPER ) THEN
+*
+*        Factorize A as U*D*U**T using the upper triangle of A
+*
+*        K is the main loop index, decreasing from N to 1 in steps of
+*        KB, where KB is the number of columns factorized by CLASYF_ROOK;
+*        KB is either NB or NB-1, or K for the last block
+*
+         K = N
+   10    CONTINUE
+*
+*        If K < 1, exit from loop
+*
+         IF( K.LT.1 )
+     $      GO TO 40
+*
+         IF( K.GT.NB ) THEN
+*
+*           Factorize columns k-kb+1:k of A and use blocked code to
+*           update columns 1:k-kb
+*
+            CALL CLASYF_ROOK( UPLO, K, NB, KB, A, LDA,
+     $                        IPIV, WORK, LDWORK, IINFO )
+         ELSE
+*
+*           Use unblocked code to factorize columns 1:k of A
+*
+            CALL CSYTF2_ROOK( UPLO, K, A, LDA, IPIV, IINFO )
+            KB = K
+         END IF
+*
+*        Set INFO on the first occurrence of a zero pivot
+*
+         IF( INFO.EQ.0 .AND. IINFO.GT.0 )
+     $      INFO = IINFO     
+*
+*        No need to adjust IPIV
+*
+*        Decrease K and return to the start of the main loop
+*
+         K = K - KB
+         GO TO 10
+*
+      ELSE
+*
+*        Factorize A as L*D*L**T using the lower triangle of A
+*
+*        K is the main loop index, increasing from 1 to N in steps of
+*        KB, where KB is the number of columns factorized by CLASYF_ROOK;
+*        KB is either NB or NB-1, or N-K+1 for the last block
+*
+         K = 1
+   20    CONTINUE
+*
+*        If K > N, exit from loop
+*
+         IF( K.GT.N )
+     $      GO TO 40
+*
+         IF( K.LE.N-NB ) THEN
+*
+*           Factorize columns k:k+kb-1 of A and use blocked code to
+*           update columns k+kb:n
+*
+            CALL CLASYF_ROOK( UPLO, N-K+1, NB, KB, A( K, K ), LDA,
+     $                        IPIV( K ), WORK, LDWORK, IINFO )
+         ELSE
+*
+*           Use unblocked code to factorize columns k:n of A
+*
+            CALL CSYTF2_ROOK( UPLO, N-K+1, A( K, K ), LDA, IPIV( K ),
+     $                   IINFO )
+            KB = N - K + 1
+         END IF
+*
+*        Set INFO on the first occurrence of a zero pivot
+*
+         IF( INFO.EQ.0 .AND. IINFO.GT.0 )
+     $      INFO = IINFO + K - 1
+*
+*        Adjust IPIV
+*
+         DO 30 J = K, K + KB - 1
+            IF( IPIV( J ).GT.0 ) THEN
+               IPIV( J ) = IPIV( J ) + K - 1
+            ELSE
+               IPIV( J ) = IPIV( J ) - K + 1
+            END IF
+   30    CONTINUE
+*
+*        Increase K and return to the start of the main loop
+*
+         K = K + KB
+         GO TO 20
+*
+      END IF
+*
+   40 CONTINUE
+      WORK( 1 ) = LWKOPT
+      RETURN
+*
+*     End of CSYTRF_ROOK
+*
+      END
diff --git a/SRC/csytri_rook.f b/SRC/csytri_rook.f
new file mode 100644
index 00000000..0b57a713
--- /dev/null
+++ b/SRC/csytri_rook.f
@@ -0,0 +1,451 @@
+*> \brief \b CSYTRI_ROOK
+*
+*  =========== DOCUMENTATION ===========
+*
+* Online html documentation available at 
+*            http://www.netlib.org/lapack/explore-html/ 
+*
+*> \htmlonly
+*> Download CSYTRI_ROOK + dependencies 
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/csytri_rook.f"> 
+*> [TGZ]</a> 
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/csytri_rook.f"> 
+*> [ZIP]</a> 
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/csytri_rook.f"> 
+*> [TXT]</a>
+*> \endhtmlonly 
+*
+*  Definition:
+*  ===========
+*
+*       SUBROUTINE CSYTRI_ROOK( UPLO, N, A, LDA, IPIV, WORK, INFO )
+* 
+*       .. Scalar Arguments ..
+*       CHARACTER          UPLO
+*       INTEGER            INFO, LDA, N
+*       ..
+*       .. Array Arguments ..
+*       INTEGER            IPIV( * )
+*       COMPLEX            A( LDA, * ), WORK( * )
+*       ..
+*  
+*
+*> \par Purpose:
+*  =============
+*>
+*> \verbatim
+*>
+*> CSYTRI_ROOK computes the inverse of a complex symmetric
+*> matrix A using the factorization A = U*D*U**T or A = L*D*L**T
+*> computed by CSYTRF_ROOK.
+*> \endverbatim
+*
+*  Arguments:
+*  ==========
+*
+*> \param[in] UPLO
+*> \verbatim
+*>          UPLO is CHARACTER*1
+*>          Specifies whether the details of the factorization are stored
+*>          as an upper or lower triangular matrix.
+*>          = 'U':  Upper triangular, form is A = U*D*U**T;
+*>          = 'L':  Lower triangular, form is A = L*D*L**T.
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*>          N is INTEGER
+*>          The order of the matrix A.  N >= 0.
+*> \endverbatim
+*>
+*> \param[in,out] A
+*> \verbatim
+*>          A is COMPLEX array, dimension (LDA,N)
+*>          On entry, the block diagonal matrix D and the multipliers
+*>          used to obtain the factor U or L as computed by CSYTRF_ROOK.
+*>
+*>          On exit, if INFO = 0, the (symmetric) inverse of the original
+*>          matrix.  If UPLO = 'U', the upper triangular part of the
+*>          inverse is formed and the part of A below the diagonal is not
+*>          referenced; if UPLO = 'L' the lower triangular part of the
+*>          inverse is formed and the part of A above the diagonal is
+*>          not referenced.
+*> \endverbatim
+*>
+*> \param[in] LDA
+*> \verbatim
+*>          LDA is INTEGER
+*>          The leading dimension of the array A.  LDA >= max(1,N).
+*> \endverbatim
+*>
+*> \param[in] IPIV
+*> \verbatim
+*>          IPIV is INTEGER array, dimension (N)
+*>          Details of the interchanges and the block structure of D
+*>          as determined by CSYTRF_ROOK.
+*> \endverbatim
+*>
+*> \param[out] WORK
+*> \verbatim
+*>          WORK is COMPLEX array, dimension (N)
+*> \endverbatim
+*>
+*> \param[out] INFO
+*> \verbatim
+*>          INFO is INTEGER
+*>          = 0: successful exit
+*>          < 0: if INFO = -i, the i-th argument had an illegal value
+*>          > 0: if INFO = i, D(i,i) = 0; the matrix is singular and its
+*>               inverse could not be computed.
+*> \endverbatim
+*
+*  Authors:
+*  ========
+*
+*> \author Univ. of Tennessee 
+*> \author Univ. of California Berkeley 
+*> \author Univ. of Colorado Denver 
+*> \author NAG Ltd. 
+*
+*> \date November 2011
+*
+*> \ingroup complexSYcomputational
+*
+*> \par Contributors:
+*  ==================
+*>
+*> \verbatim
+*>
+*>   November 2011, Igor Kozachenko,
+*>                  Computer Science Division,
+*>                  University of California, Berkeley
+*>
+*>  September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas,
+*>                  School of Mathematics,
+*>                  University of Manchester
+*>
+*> \endverbatim
+*
+*  =====================================================================
+      SUBROUTINE CSYTRI_ROOK( UPLO, N, A, LDA, IPIV, WORK, INFO )
+*
+*  -- LAPACK computational routine (version 3.4.0) --
+*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*     November 2011
+*
+*     .. Scalar Arguments ..
+      CHARACTER          UPLO
+      INTEGER            INFO, LDA, N
+*     ..
+*     .. Array Arguments ..
+      INTEGER            IPIV( * )
+      COMPLEX            A( LDA, * ), WORK( * )
+*     ..
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      COMPLEX            CONE, CZERO
+      PARAMETER          ( CONE = ( 1.0E+0, 0.0E+0 ),
+     $                   CZERO = ( 0.0E+0, 0.0E+0 ) )
+*     ..
+*     .. Local Scalars ..
+      LOGICAL            UPPER
+      INTEGER            K, KP, KSTEP
+      COMPLEX            AK, AKKP1, AKP1, D, T, TEMP
+*     ..
+*     .. External Functions ..
+      LOGICAL            LSAME
+      COMPLEX            CDOTU
+      EXTERNAL           LSAME, CDOTU
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           CCOPY, CSWAP, CSYMV, XERBLA
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          MAX
+*     ..
+*     .. Executable Statements ..
+*
+*     Test the input parameters.
+*
+      INFO = 0
+      UPPER = LSAME( UPLO, 'U' )
+      IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
+         INFO = -1
+      ELSE IF( N.LT.0 ) THEN
+         INFO = -2
+      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
+         INFO = -4
+      END IF
+      IF( INFO.NE.0 ) THEN
+         CALL XERBLA( 'CSYTRI_ROOK', -INFO )
+         RETURN
+      END IF
+*
+*     Quick return if possible
+*
+      IF( N.EQ.0 )
+     $   RETURN
+*
+*     Check that the diagonal matrix D is nonsingular.
+*
+      IF( UPPER ) THEN
+*
+*        Upper triangular storage: examine D from bottom to top
+*
+         DO 10 INFO = N, 1, -1
+            IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.CZERO )
+     $         RETURN
+   10    CONTINUE
+      ELSE
+*
+*        Lower triangular storage: examine D from top to bottom.
+*
+         DO 20 INFO = 1, N
+            IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.CZERO )
+     $         RETURN
+   20    CONTINUE
+      END IF
+      INFO = 0
+*
+      IF( UPPER ) THEN
+*
+*        Compute inv(A) from the factorization A = U*D*U**T.
+*
+*        K is the main loop index, increasing from 1 to N in steps of
+*        1 or 2, depending on the size of the diagonal blocks.
+*
+         K = 1
+   30    CONTINUE
+*
+*        If K > N, exit from loop.
+*
+         IF( K.GT.N )
+     $      GO TO 40
+*
+         IF( IPIV( K ).GT.0 ) THEN
+*
+*           1 x 1 diagonal block
+*
+*           Invert the diagonal block.
+*
+            A( K, K ) = CONE / A( K, K )
+*
+*           Compute column K of the inverse.
+*
+            IF( K.GT.1 ) THEN
+               CALL CCOPY( K-1, A( 1, K ), 1, WORK, 1 )
+               CALL CSYMV( UPLO, K-1, -CONE, A, LDA, WORK, 1, CZERO,
+     $                     A( 1, K ), 1 )
+               A( K, K ) = A( K, K ) - CDOTU( K-1, WORK, 1, A( 1, K ),
+     $                     1 )
+            END IF
+            KSTEP = 1
+         ELSE
+*
+*           2 x 2 diagonal block
+*
+*           Invert the diagonal block.
+*
+            T = A( K, K+1 )
+            AK = A( K, K ) / T
+            AKP1 = A( K+1, K+1 ) / T
+            AKKP1 = A( K, K+1 ) / T
+            D = T*( AK*AKP1-CONE )
+            A( K, K ) = AKP1 / D
+            A( K+1, K+1 ) = AK / D
+            A( K, K+1 ) = -AKKP1 / D
+*
+*           Compute columns K and K+1 of the inverse.
+*
+            IF( K.GT.1 ) THEN
+               CALL CCOPY( K-1, A( 1, K ), 1, WORK, 1 )
+               CALL CSYMV( UPLO, K-1, -CONE, A, LDA, WORK, 1, CZERO,
+     $                     A( 1, K ), 1 )
+               A( K, K ) = A( K, K ) - CDOTU( K-1, WORK, 1, A( 1, K ),
+     $                     1 )
+               A( K, K+1 ) = A( K, K+1 ) -
+     $                       CDOTU( K-1, A( 1, K ), 1, A( 1, K+1 ), 1 )
+               CALL CCOPY( K-1, A( 1, K+1 ), 1, WORK, 1 )
+               CALL CSYMV( UPLO, K-1, -CONE, A, LDA, WORK, 1, CZERO,
+     $                     A( 1, K+1 ), 1 )
+               A( K+1, K+1 ) = A( K+1, K+1 ) -
+     $                         CDOTU( K-1, WORK, 1, A( 1, K+1 ), 1 )
+            END IF
+            KSTEP = 2
+         END IF
+*
+         IF( KSTEP.EQ.1 ) THEN
+*
+*           Interchange rows and columns K and IPIV(K) in the leading
+*           submatrix A(1:k+1,1:k+1)
+*
+            KP = IPIV( K )
+            IF( KP.NE.K ) THEN
+               IF( KP.GT.1 )
+     $             CALL CSWAP( KP-1, A( 1, K ), 1, A( 1, KP ), 1 )
+               CALL CSWAP( K-KP-1, A( KP+1, K ), 1, A( KP, KP+1 ), LDA )
+               TEMP = A( K, K )
+               A( K, K ) = A( KP, KP )
+               A( KP, KP ) = TEMP
+            END IF
+         ELSE
+*
+*           Interchange rows and columns K and K+1 with -IPIV(K) and
+*           -IPIV(K+1)in the leading submatrix A(1:k+1,1:k+1)
+*
+            KP = -IPIV( K )
+            IF( KP.NE.K ) THEN
+               IF( KP.GT.1 )
+     $            CALL CSWAP( KP-1, A( 1, K ), 1, A( 1, KP ), 1 )
+               CALL CSWAP( K-KP-1, A( KP+1, K ), 1, A( KP, KP+1 ), LDA )
+*               
+               TEMP = A( K, K )
+               A( K, K ) = A( KP, KP )
+               A( KP, KP ) = TEMP
+               TEMP = A( K, K+1 )
+               A( K, K+1 ) = A( KP, K+1 )
+               A( KP, K+1 ) = TEMP
+            END IF
+*
+            K = K + 1
+            KP = -IPIV( K )
+            IF( KP.NE.K ) THEN
+               IF( KP.GT.1 )
+     $            CALL CSWAP( KP-1, A( 1, K ), 1, A( 1, KP ), 1 )
+               CALL CSWAP( K-KP-1, A( KP+1, K ), 1, A( KP, KP+1 ), LDA )
+               TEMP = A( K, K )
+               A( K, K ) = A( KP, KP )
+               A( KP, KP ) = TEMP
+            END IF
+         END IF
+*
+         K = K + 1
+         GO TO 30
+   40    CONTINUE
+*
+      ELSE
+*
+*        Compute inv(A) from the factorization A = L*D*L**T.
+*
+*        K is the main loop index, increasing from 1 to N in steps of
+*        1 or 2, depending on the size of the diagonal blocks.
+*
+         K = N
+   50    CONTINUE
+*
+*        If K < 1, exit from loop.
+*
+         IF( K.LT.1 )
+     $      GO TO 60
+*
+         IF( IPIV( K ).GT.0 ) THEN
+*
+*           1 x 1 diagonal block
+*
+*           Invert the diagonal block.
+*
+            A( K, K ) = CONE / A( K, K )
+*
+*           Compute column K of the inverse.
+*
+            IF( K.LT.N ) THEN
+               CALL CCOPY( N-K, A( K+1, K ), 1, WORK, 1 )
+               CALL CSYMV( UPLO, N-K,-CONE, A( K+1, K+1 ), LDA, WORK, 1,
+     $                     CZERO, A( K+1, K ), 1 )
+               A( K, K ) = A( K, K ) - CDOTU( N-K, WORK, 1, A( K+1, K ),
+     $                     1 )
+            END IF
+            KSTEP = 1
+         ELSE
+*
+*           2 x 2 diagonal block
+*
+*           Invert the diagonal block.
+*
+            T = A( K, K-1 )
+            AK = A( K-1, K-1 ) / T
+            AKP1 = A( K, K ) / T
+            AKKP1 = A( K, K-1 ) / T
+            D = T*( AK*AKP1-CONE )
+            A( K-1, K-1 ) = AKP1 / D
+            A( K, K ) = AK / D
+            A( K, K-1 ) = -AKKP1 / D
+*
+*           Compute columns K-1 and K of the inverse.
+*
+            IF( K.LT.N ) THEN
+               CALL CCOPY( N-K, A( K+1, K ), 1, WORK, 1 )
+               CALL CSYMV( UPLO, N-K,-CONE, A( K+1, K+1 ), LDA, WORK, 1,
+     $                     CZERO, A( K+1, K ), 1 )
+               A( K, K ) = A( K, K ) - CDOTU( N-K, WORK, 1, A( K+1, K ),
+     $                     1 )
+               A( K, K-1 ) = A( K, K-1 ) -
+     $                       CDOTU( N-K, A( K+1, K ), 1, A( K+1, K-1 ),
+     $                       1 )
+               CALL CCOPY( N-K, A( K+1, K-1 ), 1, WORK, 1 )
+               CALL CSYMV( UPLO, N-K,-CONE, A( K+1, K+1 ), LDA, WORK, 1,
+     $                     CZERO, A( K+1, K-1 ), 1 )
+               A( K-1, K-1 ) = A( K-1, K-1 ) -
+     $                         CDOTU( N-K, WORK, 1, A( K+1, K-1 ), 1 )
+            END IF
+            KSTEP = 2
+         END IF  
+*
+         IF( KSTEP.EQ.1 ) THEN
+*
+*           Interchange rows and columns K and IPIV(K) in the trailing
+*           submatrix A(k-1:n,k-1:n)
+*
+            KP = IPIV( K )
+            IF( KP.NE.K ) THEN
+               IF( KP.LT.N )
+     $            CALL CSWAP( N-KP, A( KP+1, K ), 1, A( KP+1, KP ), 1 )
+               CALL CSWAP( KP-K-1, A( K+1, K ), 1, A( KP, K+1 ), LDA )
+               TEMP = A( K, K )
+               A( K, K ) = A( KP, KP )
+               A( KP, KP ) = TEMP
+            END IF
+         ELSE
+*
+*           Interchange rows and columns K and K-1 with -IPIV(K) and
+*           -IPIV(K-1) in the trailing submatrix A(k-1:n,k-1:n)
+*
+            KP = -IPIV( K )
+            IF( KP.NE.K ) THEN
+               IF( KP.LT.N )
+     $            CALL CSWAP( N-KP, A( KP+1, K ), 1, A( KP+1, KP ), 1 )
+               CALL CSWAP( KP-K-1, A( K+1, K ), 1, A( KP, K+1 ), LDA )
+*
+               TEMP = A( K, K )
+               A( K, K ) = A( KP, KP )
+               A( KP, KP ) = TEMP
+               TEMP = A( K, K-1 )
+               A( K, K-1 ) = A( KP, K-1 )
+               A( KP, K-1 ) = TEMP
+            END IF
+*
+            K = K - 1
+            KP = -IPIV( K )
+            IF( KP.NE.K ) THEN
+               IF( KP.LT.N )
+     $            CALL CSWAP( N-KP, A( KP+1, K ), 1, A( KP+1, KP ), 1 )
+               CALL CSWAP( KP-K-1, A( K+1, K ), 1, A( KP, K+1 ), LDA )
+               TEMP = A( K, K )
+               A( K, K ) = A( KP, KP )
+               A( KP, KP ) = TEMP
+            END IF
+         END IF
+*
+         K = K - 1
+         GO TO 50
+   60    CONTINUE
+      END IF
+*
+      RETURN
+*
+*     End of CSYTRI_ROOK
+*
+      END
diff --git a/SRC/csytrs_rook.f b/SRC/csytrs_rook.f
new file mode 100644
index 00000000..44727212
--- /dev/null
+++ b/SRC/csytrs_rook.f
@@ -0,0 +1,484 @@
+*> \brief \b CSYTRS_ROOK
+*
+*  =========== DOCUMENTATION ===========
+*
+* Online html documentation available at 
+*            http://www.netlib.org/lapack/explore-html/ 
+*
+*> \htmlonly
+*> Download CSYTRS_ROOK + dependencies 
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/csytrs_rook.f"> 
+*> [TGZ]</a> 
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/csytrs_rook.f"> 
+*> [ZIP]</a> 
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/csytrs_rook.f"> 
+*> [TXT]</a>
+*> \endhtmlonly 
+*
+*  Definition:
+*  ===========
+*
+*       SUBROUTINE CSYTRS_ROOK( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, INFO )
+* 
+*       .. Scalar Arguments ..
+*       CHARACTER          UPLO
+*       INTEGER            INFO, LDA, LDB, N, NRHS
+*       ..
+*       .. Array Arguments ..
+*       INTEGER            IPIV( * )
+*       COMPLEX            A( LDA, * ), B( LDB, * )
+*       ..
+*  
+*
+*> \par Purpose:
+*  =============
+*>
+*> \verbatim
+*>
+*> CSYTRS_ROOK solves a system of linear equations A*X = B with
+*> a complex symmetric matrix A using the factorization A = U*D*U**T or
+*> A = L*D*L**T computed by CSYTRF_ROOK.
+*> \endverbatim
+*
+*  Arguments:
+*  ==========
+*
+*> \param[in] UPLO
+*> \verbatim
+*>          UPLO is CHARACTER*1
+*>          Specifies whether the details of the factorization are stored
+*>          as an upper or lower triangular matrix.
+*>          = 'U':  Upper triangular, form is A = U*D*U**T;
+*>          = 'L':  Lower triangular, form is A = L*D*L**T.
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*>          N is INTEGER
+*>          The order of the matrix A.  N >= 0.
+*> \endverbatim
+*>
+*> \param[in] NRHS
+*> \verbatim
+*>          NRHS is INTEGER
+*>          The number of right hand sides, i.e., the number of columns
+*>          of the matrix B.  NRHS >= 0.
+*> \endverbatim
+*>
+*> \param[in] A
+*> \verbatim
+*>          A is COMPLEX array, dimension (LDA,N)
+*>          The block diagonal matrix D and the multipliers used to
+*>          obtain the factor U or L as computed by CSYTRF_ROOK.
+*> \endverbatim
+*>
+*> \param[in] LDA
+*> \verbatim
+*>          LDA is INTEGER
+*>          The leading dimension of the array A.  LDA >= max(1,N).
+*> \endverbatim
+*>
+*> \param[in] IPIV
+*> \verbatim
+*>          IPIV is INTEGER array, dimension (N)
+*>          Details of the interchanges and the block structure of D
+*>          as determined by CSYTRF_ROOK.
+*> \endverbatim
+*>
+*> \param[in,out] B
+*> \verbatim
+*>          B is COMPLEX array, dimension (LDB,NRHS)
+*>          On entry, the right hand side matrix B.
+*>          On exit, the solution matrix X.
+*> \endverbatim
+*>
+*> \param[in] LDB
+*> \verbatim
+*>          LDB is INTEGER
+*>          The leading dimension of the array B.  LDB >= max(1,N).
+*> \endverbatim
+*>
+*> \param[out] INFO
+*> \verbatim
+*>          INFO is INTEGER
+*>          = 0:  successful exit
+*>          < 0:  if INFO = -i, the i-th argument had an illegal value
+*> \endverbatim
+*
+*  Authors:
+*  ========
+*
+*> \author Univ. of Tennessee 
+*> \author Univ. of California Berkeley 
+*> \author Univ. of Colorado Denver 
+*> \author NAG Ltd. 
+*
+*> \date November 2011
+*
+*> \ingroup complexSYcomputational
+*
+*> \par Contributors:
+*  ==================
+*>
+*> \verbatim
+*>
+*>   November 2011, Igor Kozachenko,
+*>                  Computer Science Division,
+*>                  University of California, Berkeley
+*>
+*>  September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas,
+*>                  School of Mathematics,
+*>                  University of Manchester
+*>
+*> \endverbatim
+*
+*  =====================================================================
+      SUBROUTINE CSYTRS_ROOK( UPLO, N, NRHS, A, LDA, IPIV, B, LDB,
+     $                        INFO )
+*
+*  -- LAPACK computational routine (version 3.4.0) --
+*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*     November 2011
+*
+*     .. Scalar Arguments ..
+      CHARACTER          UPLO
+      INTEGER            INFO, LDA, LDB, N, NRHS
+*     ..
+*     .. Array Arguments ..
+      INTEGER            IPIV( * )
+      COMPLEX            A( LDA, * ), B( LDB, * )
+*     ..
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      COMPLEX            CONE
+      PARAMETER          ( CONE = ( 1.0E+0, 0.0E+0 ) )
+*     ..
+*     .. Local Scalars ..
+      LOGICAL            UPPER
+      INTEGER            J, K, KP
+      COMPLEX            AK, AKM1, AKM1K, BK, BKM1, DENOM
+*     ..
+*     .. External Functions ..
+      LOGICAL            LSAME
+      EXTERNAL           LSAME
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           CGEMV, CGERU, CSCAL, CSWAP, XERBLA
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          MAX
+*     ..
+*     .. Executable Statements ..
+*
+      INFO = 0
+      UPPER = LSAME( UPLO, 'U' )
+      IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
+         INFO = -1
+      ELSE IF( N.LT.0 ) THEN
+         INFO = -2
+      ELSE IF( NRHS.LT.0 ) THEN
+         INFO = -3
+      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
+         INFO = -5
+      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
+         INFO = -8
+      END IF
+      IF( INFO.NE.0 ) THEN
+         CALL XERBLA( 'CSYTRS_ROOK', -INFO )
+         RETURN
+      END IF
+*
+*     Quick return if possible
+*
+      IF( N.EQ.0 .OR. NRHS.EQ.0 )
+     $   RETURN
+*
+      IF( UPPER ) THEN
+*
+*        Solve A*X = B, where A = U*D*U**T.
+*
+*        First solve U*D*X = B, overwriting B with X.
+*
+*        K is the main loop index, decreasing from N to 1 in steps of
+*        1 or 2, depending on the size of the diagonal blocks.
+*
+         K = N
+   10    CONTINUE
+*
+*        If K < 1, exit from loop.
+*
+         IF( K.LT.1 )
+     $      GO TO 30
+*
+         IF( IPIV( K ).GT.0 ) THEN
+*
+*           1 x 1 diagonal block
+*
+*           Interchange rows K and IPIV(K).
+*
+            KP = IPIV( K )
+            IF( KP.NE.K )
+     $         CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+*
+*           Multiply by inv(U(K)), where U(K) is the transformation
+*           stored in column K of A.
+*
+            CALL CGERU( K-1, NRHS, -CONE, A( 1, K ), 1, B( K, 1 ), LDB,
+     $                 B( 1, 1 ), LDB )
+*
+*           Multiply by the inverse of the diagonal block.
+*
+            CALL CSCAL( NRHS, CONE / A( K, K ), B( K, 1 ), LDB )
+            K = K - 1
+         ELSE
+*
+*           2 x 2 diagonal block
+*
+*           Interchange rows K and -IPIV(K) THEN K-1 and -IPIV(K-1)
+*
+            KP = -IPIV( K )
+            IF( KP.NE.K )
+     $         CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+*
+            KP = -IPIV( K-1 )
+            IF( KP.NE.K-1 )
+     $         CALL CSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ), LDB )
+*
+*           Multiply by inv(U(K)), where U(K) is the transformation
+*           stored in columns K-1 and K of A.
+*
+            IF( K.GT.2 ) THEN
+               CALL CGERU( K-2, NRHS,-CONE, A( 1, K ), 1, B( K, 1 ),
+     $                    LDB, B( 1, 1 ), LDB )
+               CALL CGERU( K-2, NRHS,-CONE, A( 1, K-1 ), 1, B( K-1, 1 ),
+     $                    LDB, B( 1, 1 ), LDB )
+            END IF
+*
+*           Multiply by the inverse of the diagonal block.
+*
+            AKM1K = A( K-1, K )
+            AKM1 = A( K-1, K-1 ) / AKM1K
+            AK = A( K, K ) / AKM1K
+            DENOM = AKM1*AK - CONE
+            DO 20 J = 1, NRHS
+               BKM1 = B( K-1, J ) / AKM1K
+               BK = B( K, J ) / AKM1K
+               B( K-1, J ) = ( AK*BKM1-BK ) / DENOM
+               B( K, J ) = ( AKM1*BK-BKM1 ) / DENOM
+   20       CONTINUE
+            K = K - 2
+         END IF
+*
+         GO TO 10
+   30    CONTINUE
+*
+*        Next solve U**T *X = B, overwriting B with X.
+*
+*        K is the main loop index, increasing from 1 to N in steps of
+*        1 or 2, depending on the size of the diagonal blocks.
+*
+         K = 1
+   40    CONTINUE
+*
+*        If K > N, exit from loop.
+*
+         IF( K.GT.N )
+     $      GO TO 50
+*
+         IF( IPIV( K ).GT.0 ) THEN
+*
+*           1 x 1 diagonal block
+*
+*           Multiply by inv(U**T(K)), where U(K) is the transformation
+*           stored in column K of A.
+*
+            IF( K.GT.1 )
+     $         CALL CGEMV( 'Transpose', K-1, NRHS, -CONE, B,
+     $                     LDB, A( 1, K ), 1, CONE, B( K, 1 ), LDB )
+*
+*           Interchange rows K and IPIV(K).
+*
+            KP = IPIV( K )
+            IF( KP.NE.K )
+     $         CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+            K = K + 1
+         ELSE
+*
+*           2 x 2 diagonal block
+*
+*           Multiply by inv(U**T(K+1)), where U(K+1) is the transformation
+*           stored in columns K and K+1 of A.
+*
+            IF( K.GT.1 ) THEN
+               CALL CGEMV( 'Transpose', K-1, NRHS, -CONE, B,
+     $                     LDB, A( 1, K ), 1, CONE, B( K, 1 ), LDB )
+               CALL CGEMV( 'Transpose', K-1, NRHS, -CONE, B,
+     $                     LDB, A( 1, K+1 ), 1, CONE, B( K+1, 1 ), LDB )
+            END IF
+*
+*           Interchange rows K and -IPIV(K) THEN K+1 and -IPIV(K+1).
+*
+            KP = -IPIV( K )
+            IF( KP.NE.K )
+     $         CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+*
+            KP = -IPIV( K+1 )
+            IF( KP.NE.K+1 )
+     $         CALL CSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ), LDB )
+*
+            K = K + 2
+         END IF
+*
+         GO TO 40
+   50    CONTINUE
+*
+      ELSE
+*
+*        Solve A*X = B, where A = L*D*L**T.
+*
+*        First solve L*D*X = B, overwriting B with X.
+*
+*        K is the main loop index, increasing from 1 to N in steps of
+*        1 or 2, depending on the size of the diagonal blocks.
+*
+         K = 1
+   60    CONTINUE
+*
+*        If K > N, exit from loop.
+*
+         IF( K.GT.N )
+     $      GO TO 80
+*
+         IF( IPIV( K ).GT.0 ) THEN
+*
+*           1 x 1 diagonal block
+*
+*           Interchange rows K and IPIV(K).
+*
+            KP = IPIV( K )
+            IF( KP.NE.K )
+     $         CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+*
+*           Multiply by inv(L(K)), where L(K) is the transformation
+*           stored in column K of A.
+*
+            IF( K.LT.N )
+     $         CALL CGERU( N-K, NRHS, -CONE, A( K+1, K ), 1, B( K, 1 ),
+     $                    LDB, B( K+1, 1 ), LDB )
+*
+*           Multiply by the inverse of the diagonal block.
+*
+            CALL CSCAL( NRHS, CONE / A( K, K ), B( K, 1 ), LDB )
+            K = K + 1
+         ELSE
+*
+*           2 x 2 diagonal block
+*
+*           Interchange rows K and -IPIV(K) THEN K+1 and -IPIV(K+1)
+*
+            KP = -IPIV( K )
+            IF( KP.NE.K )
+     $         CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+*
+            KP = -IPIV( K+1 )
+            IF( KP.NE.K+1 )
+     $         CALL CSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ), LDB )
+*
+*           Multiply by inv(L(K)), where L(K) is the transformation
+*           stored in columns K and K+1 of A.
+*
+            IF( K.LT.N-1 ) THEN
+               CALL CGERU( N-K-1, NRHS,-CONE, A( K+2, K ), 1, B( K, 1 ),
+     $                    LDB, B( K+2, 1 ), LDB )
+               CALL CGERU( N-K-1, NRHS,-CONE, A( K+2, K+1 ), 1,
+     $                    B( K+1, 1 ), LDB, B( K+2, 1 ), LDB )
+            END IF
+*
+*           Multiply by the inverse of the diagonal block.
+*
+            AKM1K = A( K+1, K )
+            AKM1 = A( K, K ) / AKM1K
+            AK = A( K+1, K+1 ) / AKM1K
+            DENOM = AKM1*AK - CONE
+            DO 70 J = 1, NRHS
+               BKM1 = B( K, J ) / AKM1K
+               BK = B( K+1, J ) / AKM1K
+               B( K, J ) = ( AK*BKM1-BK ) / DENOM
+               B( K+1, J ) = ( AKM1*BK-BKM1 ) / DENOM
+   70       CONTINUE
+            K = K + 2
+         END IF
+*
+         GO TO 60
+   80    CONTINUE
+*
+*        Next solve L**T *X = B, overwriting B with X.
+*
+*        K is the main loop index, decreasing from N to 1 in steps of
+*        1 or 2, depending on the size of the diagonal blocks.
+*
+         K = N
+   90    CONTINUE
+*
+*        If K < 1, exit from loop.
+*
+         IF( K.LT.1 )
+     $      GO TO 100
+*
+         IF( IPIV( K ).GT.0 ) THEN
+*
+*           1 x 1 diagonal block
+*
+*           Multiply by inv(L**T(K)), where L(K) is the transformation
+*           stored in column K of A.
+*
+            IF( K.LT.N )
+     $         CALL CGEMV( 'Transpose', N-K, NRHS, -CONE, B( K+1, 1 ),
+     $                     LDB, A( K+1, K ), 1, CONE, B( K, 1 ), LDB )
+*
+*           Interchange rows K and IPIV(K).
+*
+            KP = IPIV( K )
+            IF( KP.NE.K )
+     $         CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+            K = K - 1
+         ELSE
+*
+*           2 x 2 diagonal block
+*
+*           Multiply by inv(L**T(K-1)), where L(K-1) is the transformation
+*           stored in columns K-1 and K of A.
+*
+            IF( K.LT.N ) THEN
+               CALL CGEMV( 'Transpose', N-K, NRHS, -CONE, B( K+1, 1 ),
+     $                     LDB, A( K+1, K ), 1, CONE, B( K, 1 ), LDB )
+               CALL CGEMV( 'Transpose', N-K, NRHS, -CONE, B( K+1, 1 ),
+     $                     LDB, A( K+1, K-1 ), 1, CONE, B( K-1, 1 ),
+     $                     LDB )
+            END IF
+*
+*           Interchange rows K and -IPIV(K) THEN K-1 and -IPIV(K-1)
+*
+            KP = -IPIV( K )
+            IF( KP.NE.K )
+     $         CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+*
+            KP = -IPIV( K-1 )
+            IF( KP.NE.K-1 )
+     $         CALL CSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ), LDB )
+*
+            K = K - 2
+         END IF
+*
+         GO TO 90
+  100    CONTINUE
+      END IF
+*
+      RETURN
+*
+*     End of CSYTRS_ROOK
+*
+      END