added test routines (c,z)chkhe_rook.f and (c,z)drvhe_rook.f for Hermitian factorization routines with rook pivoting algorithm

6 files changed, 3404 insertions, 0 deletions

diff --git a/TESTING/LIN/cchkhe_rook.f b/TESTING/LIN/cchkhe_rook.f
new file mode 100644
index 00000000..948660a2
--- /dev/null
+++ b/TESTING/LIN/cchkhe_rook.f
@@ -0,0 +1,849 @@
+*> \brief \b CCHKHE_ROOK
+*
+*  =========== DOCUMENTATION ===========
+*
+* Online html documentation available at 
+*            http://www.netlib.org/lapack/explore-html/ 
+*
+*  Definition:
+*  ===========
+*
+*       SUBROUTINE CCHKHE_ROOK( DOTYPE, NN, NVAL, NNB, NBVAL, NNS, NSVAL,
+*                               THRESH, TSTERR, NMAX, A, AFAC, AINV, B, X,
+*                               XACT, WORK, RWORK, IWORK, NOUT )
+* 
+*       .. Scalar Arguments ..
+*       LOGICAL            TSTERR
+*       INTEGER            NMAX, NN, NNB, NNS, NOUT
+*       REAL               THRESH
+*       ..
+*       .. Array Arguments ..
+*       LOGICAL            DOTYPE( * )
+*       INTEGER            IWORK( * ), NBVAL( * ), NSVAL( * ), NVAL( * )
+*       REAL               RWORK( * )
+*       COMPLEX            A( * ), AFAC( * ), AINV( * ), B( * ),
+*      $                   WORK( * ), X( * ), XACT( * )
+*       ..
+*  
+*
+*> \par Purpose:
+*  =============
+*>
+*> \verbatim
+*>
+*> CCHKHE_ROOK tests CHETRF_ROOK, -TRI_ROOK, -TRS_ROOK,
+*> and -CON_ROOK.
+*> \endverbatim
+*
+*  Arguments:
+*  ==========
+*
+*> \param[in] DOTYPE
+*> \verbatim
+*>          DOTYPE is LOGICAL array, dimension (NTYPES)
+*>          The matrix types to be used for testing.  Matrices of type j
+*>          (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) =
+*>          .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used.
+*> \endverbatim
+*>
+*> \param[in] NN
+*> \verbatim
+*>          NN is INTEGER
+*>          The number of values of N contained in the vector NVAL.
+*> \endverbatim
+*>
+*> \param[in] NVAL
+*> \verbatim
+*>          NVAL is INTEGER array, dimension (NN)
+*>          The values of the matrix dimension N.
+*> \endverbatim
+*>
+*> \param[in] NNB
+*> \verbatim
+*>          NNB is INTEGER
+*>          The number of values of NB contained in the vector NBVAL.
+*> \endverbatim
+*>
+*> \param[in] NBVAL
+*> \verbatim
+*>          NBVAL is INTEGER array, dimension (NBVAL)
+*>          The values of the blocksize NB.
+*> \endverbatim
+*>
+*> \param[in] NNS
+*> \verbatim
+*>          NNS is INTEGER
+*>          The number of values of NRHS contained in the vector NSVAL.
+*> \endverbatim
+*>
+*> \param[in] NSVAL
+*> \verbatim
+*>          NSVAL is INTEGER array, dimension (NNS)
+*>          The values of the number of right hand sides NRHS.
+*> \endverbatim
+*>
+*> \param[in] THRESH
+*> \verbatim
+*>          THRESH is REAL
+*>          The threshold value for the test ratios.  A result is
+*>          included in the output file if RESULT >= THRESH.  To have
+*>          every test ratio printed, use THRESH = 0.
+*> \endverbatim
+*>
+*> \param[in] TSTERR
+*> \verbatim
+*>          TSTERR is LOGICAL
+*>          Flag that indicates whether error exits are to be tested.
+*> \endverbatim
+*>
+*> \param[in] NMAX
+*> \verbatim
+*>          NMAX is INTEGER
+*>          The maximum value permitted for N, used in dimensioning the
+*>          work arrays.
+*> \endverbatim
+*>
+*> \param[out] A
+*> \verbatim
+*>          A is COMPLEX array, dimension (NMAX*NMAX)
+*> \endverbatim
+*>
+*> \param[out] AFAC
+*> \verbatim
+*>          AFAC is COMPLEX array, dimension (NMAX*NMAX)
+*> \endverbatim
+*>
+*> \param[out] AINV
+*> \verbatim
+*>          AINV is COMPLEX array, dimension (NMAX*NMAX)
+*> \endverbatim
+*>
+*> \param[out] B
+*> \verbatim
+*>          B is COMPLEX array, dimension (NMAX*NSMAX)
+*>          where NSMAX is the largest entry in NSVAL.
+*> \endverbatim
+*>
+*> \param[out] X
+*> \verbatim
+*>          X is COMPLEX array, dimension (NMAX*NSMAX)
+*> \endverbatim
+*>
+*> \param[out] XACT
+*> \verbatim
+*>          XACT is COMPLEX array, dimension (NMAX*NSMAX)
+*> \endverbatim
+*>
+*> \param[out] WORK
+*> \verbatim
+*>          WORK is COMPLEX array, dimension
+*>                      (NMAX*max(3,NSMAX))
+*> \endverbatim
+*>
+*> \param[out] RWORK
+*> \verbatim
+*>          RWORK is REAL array, dimension
+*>                      (max(NMAX,2*NSMAX))
+*> \endverbatim
+*>
+*> \param[out] IWORK
+*> \verbatim
+*>          IWORK is INTEGER array, dimension (2*NMAX)
+*> \endverbatim
+*>
+*> \param[in] NOUT
+*> \verbatim
+*>          NOUT is INTEGER
+*>          The unit number for output.
+*> \endverbatim
+*
+*  Authors:
+*  ========
+*
+*> \author Univ. of Tennessee 
+*> \author Univ. of California Berkeley 
+*> \author Univ. of Colorado Denver 
+*> \author NAG Ltd. 
+*
+*> \date April 2013
+*
+*> \ingroup complex_lin
+*
+*  =====================================================================
+      SUBROUTINE CCHKHE_ROOK( DOTYPE, NN, NVAL, NNB, NBVAL, NNS, NSVAL,
+     $                        THRESH, TSTERR, NMAX, A, AFAC, AINV, B, X,
+     $                        XACT, WORK, RWORK, IWORK, NOUT )
+*
+*  -- LAPACK test routine (version 3.4.1) --
+*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*     April 2012
+*
+*     .. Scalar Arguments ..
+      LOGICAL            TSTERR
+      INTEGER            NMAX, NN, NNB, NNS, NOUT
+      REAL               THRESH
+*     ..
+*     .. Array Arguments ..
+      LOGICAL            DOTYPE( * )
+      INTEGER            IWORK( * ), NBVAL( * ), NSVAL( * ), NVAL( * )
+      REAL               RWORK( * )
+      COMPLEX            A( * ), AFAC( * ), AINV( * ), B( * ),
+     $                   WORK( * ), X( * ), XACT( * )
+*     ..
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      REAL               ZERO, ONE
+      PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
+      REAL               ONEHALF
+      PARAMETER          ( ONEHALF = 0.5E+0 )
+      REAL               EIGHT, SEVTEN
+      PARAMETER          ( EIGHT = 8.0E+0, SEVTEN = 17.0E+0 )
+      COMPLEX            CZERO
+      PARAMETER          ( CZERO = ( 0.0E+0, 0.0E+0 )  )
+      INTEGER            NTYPES
+      PARAMETER          ( NTYPES = 10 )
+      INTEGER            NTESTS
+      PARAMETER          ( NTESTS = 7 )
+*     ..
+*     .. Local Scalars ..
+      LOGICAL            TRFCON, ZEROT
+      CHARACTER          DIST, TYPE, UPLO, XTYPE
+      CHARACTER*3        PATH, MATPATH
+      INTEGER            I, I1, I2, IMAT, IN, INB, INFO, IOFF, IRHS,
+     $                   ITEMP, ITEMP2, IUPLO, IZERO, J, K, KL, KU, LDA,
+     $                   LWORK, MODE, N, NB, NERRS, NFAIL, NIMAT, NRHS,
+     $                   NRUN, NT
+      REAL               ALPHA, ANORM, CNDNUM, CONST, LAM_MAX, LAM_MIN,
+     $                   RCOND, RCONDC, STEMP
+*     ..
+*     .. Local Arrays ..
+      CHARACTER          UPLOS( 2 )
+      INTEGER            ISEED( 4 ), ISEEDY( 4 )
+      REAL               RESULT( NTESTS )
+      COMPLEX            BLOCK( 2, 2 ), CDUMMY( 1 )
+*     ..
+*     .. External Functions ..
+      REAL               CLANGE, CLANHE, SGET06
+      EXTERNAL           CLANGE, CLANHE, SGET06
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           ALAERH, ALAHD, ALASUM, CERRHE, CHEEVX, CGET04,
+     $                   CLACPY, CLARHS, CLATB4, CLATMS, CPOT02,
+     $                   CPOT03, CHECON_ROOK, CHET01_ROOK, CHETRF_ROOK,
+     $                   CHETRI_ROOK, CHETRS_ROOK, XLAENV
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          ABS, MAX, MIN, SQRT
+*     ..
+*     .. Scalars in Common ..
+      LOGICAL            LERR, OK
+      CHARACTER*32       SRNAMT
+      INTEGER            INFOT, NUNIT
+*     ..
+*     .. Common blocks ..
+      COMMON             / INFOC / INFOT, NUNIT, OK, LERR
+      COMMON             / SRNAMC / SRNAMT
+*     ..
+*     .. Data statements ..
+      DATA               ISEEDY / 1988, 1989, 1990, 1991 /
+      DATA               UPLOS / 'U', 'L' /
+*     ..
+*     .. Executable Statements ..
+*
+*     Initialize constants and the random number seed.
+*
+      ALPHA = ( ONE+SQRT( SEVTEN ) ) / EIGHT
+*
+*     Test path
+*
+      PATH( 1: 1 ) = 'Complex precision'
+      PATH( 2: 3 ) = 'HR'
+*
+*     Path to generate matrices
+*
+      MATPATH( 1: 1 ) = 'Complex precision'
+      MATPATH( 2: 3 ) = 'HE'
+*
+      NRUN = 0
+      NFAIL = 0
+      NERRS = 0
+      DO 10 I = 1, 4
+         ISEED( I ) = ISEEDY( I )
+   10 CONTINUE
+*
+*     Test the error exits
+*
+      IF( TSTERR )
+     $   CALL CERRHE( PATH, NOUT )
+      INFOT = 0
+*
+*     Set the minimum block size for which the block routine should
+*     be used, which will be later returned by ILAENV
+*
+      CALL XLAENV( 2, 2 )
+*
+*     Do for each value of N in NVAL
+*
+      DO 270 IN = 1, NN
+         N = NVAL( IN )
+         LDA = MAX( N, 1 )
+         XTYPE = 'N'
+         NIMAT = NTYPES
+         IF( N.LE.0 )
+     $      NIMAT = 1
+*
+         IZERO = 0
+*
+*        Do for each value of matrix type IMAT
+*
+         DO 260 IMAT = 1, NIMAT
+*
+*           Do the tests only if DOTYPE( IMAT ) is true.
+*
+            IF( .NOT.DOTYPE( IMAT ) )
+     $         GO TO 260
+*
+*           Skip types 3, 4, 5, or 6 if the matrix size is too small.
+*
+            ZEROT = IMAT.GE.3 .AND. IMAT.LE.6
+            IF( ZEROT .AND. N.LT.IMAT-2 )
+     $         GO TO 260
+*
+*           Do first for UPLO = 'U', then for UPLO = 'L'
+*
+            DO 250 IUPLO = 1, 2
+               UPLO = UPLOS( IUPLO )
+*
+*                 Begin generate the test matrix A.
+*
+*                 Set up parameters with CLATB4 for the matrix generator
+*                 based on the type of matrix to be generated.
+*
+                  CALL CLATB4( MATPATH, IMAT, N, N, TYPE, KL, KU, ANORM,
+     $                         MODE, CNDNUM, DIST )
+*
+*                 Generate a matrix with CLATMS.
+*
+                  SRNAMT = 'CLATMS'
+                  CALL CLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE,
+     $                         CNDNUM, ANORM, KL, KU, UPLO, A, LDA,
+     $                         WORK, INFO )
+*
+*                 Check error code from CLATMS and handle error.
+*
+                  IF( INFO.NE.0 ) THEN
+                     CALL ALAERH( PATH, 'CLATMS', INFO, 0, UPLO, N, N,
+     $                            -1, -1, -1, IMAT, NFAIL, NERRS, NOUT )
+*
+*                    Skip all tests for this generated matrix
+*
+                     GO TO 250
+                  END IF
+*
+*                 For matrix types 3-6, zero one or more rows and
+*                 columns of the matrix to test that INFO is returned
+*                 correctly.
+*
+                  IF( ZEROT ) THEN
+                     IF( IMAT.EQ.3 ) THEN
+                        IZERO = 1
+                     ELSE IF( IMAT.EQ.4 ) THEN
+                        IZERO = N
+                     ELSE
+                        IZERO = N / 2 + 1
+                     END IF
+*
+                     IF( IMAT.LT.6 ) THEN
+*
+*                       Set row and column IZERO to zero.
+*
+                        IF( IUPLO.EQ.1 ) THEN
+                           IOFF = ( IZERO-1 )*LDA
+                           DO 20 I = 1, IZERO - 1
+                              A( IOFF+I ) = CZERO
+   20                      CONTINUE
+                           IOFF = IOFF + IZERO
+                           DO 30 I = IZERO, N
+                              A( IOFF ) = CZERO
+                              IOFF = IOFF + LDA
+   30                      CONTINUE
+                        ELSE
+                           IOFF = IZERO
+                           DO 40 I = 1, IZERO - 1
+                              A( IOFF ) = CZERO
+                              IOFF = IOFF + LDA
+   40                      CONTINUE
+                           IOFF = IOFF - IZERO
+                           DO 50 I = IZERO, N
+                              A( IOFF+I ) = CZERO
+   50                      CONTINUE
+                        END IF
+                     ELSE
+                        IOFF = 0
+                        IF( IUPLO.EQ.1 ) THEN
+*
+*                          Set the first IZERO rows and columns to zero.
+*
+                           DO 70 J = 1, N
+                              I2 = MIN( J, IZERO )
+                              DO 60 I = 1, I2
+                                 A( IOFF+I ) = CZERO
+   60                         CONTINUE
+                              IOFF = IOFF + LDA
+   70                      CONTINUE
+                        ELSE
+*
+*                          Set the last IZERO rows and columns to zero.
+*
+                           DO 90 J = 1, N
+                              I1 = MAX( J, IZERO )
+                              DO 80 I = I1, N
+                                 A( IOFF+I ) = CZERO
+   80                         CONTINUE
+                              IOFF = IOFF + LDA
+   90                      CONTINUE
+                        END IF
+                     END IF
+                  ELSE
+                     IZERO = 0
+                  END IF
+*
+*                 End generate the test matrix A.
+*
+*              Do for each value of NB in NBVAL
+*
+               DO 240 INB = 1, NNB
+*
+*                 Set the optimal blocksize, which will be later
+*                 returned by ILAENV.
+*
+                  NB = NBVAL( INB )
+                  CALL XLAENV( 1, NB )
+*
+*                 Copy the test matrix A into matrix AFAC which
+*                 will be factorized in place. This is needed to
+*                 preserve the test matrix A for subsequent tests.
+*
+                  CALL CLACPY( UPLO, N, N, A, LDA, AFAC, LDA )
+*
+*                 Compute the L*D*L**T or U*D*U**T factorization of the
+*                 matrix. IWORK stores details of the interchanges and
+*                 the block structure of D. AINV is a work array for
+*                 block factorization, LWORK is the length of AINV.
+*
+                  LWORK = MAX( 2, NB )*LDA
+                  SRNAMT = 'CHETRF_ROOK'
+                  CALL CHETRF_ROOK( UPLO, N, AFAC, LDA, IWORK, AINV,
+     $                              LWORK, INFO )
+*
+*                 Adjust the expected value of INFO to account for
+*                 pivoting.
+*
+                  K = IZERO
+                  IF( K.GT.0 ) THEN
+  100                CONTINUE
+                     IF( IWORK( K ).LT.0 ) THEN
+                        IF( IWORK( K ).NE.-K ) THEN
+                           K = -IWORK( K )
+                           GO TO 100
+                        END IF
+                     ELSE IF( IWORK( K ).NE.K ) THEN
+                        K = IWORK( K )
+                        GO TO 100
+                     END IF
+                  END IF
+*
+*                 Check error code from CHETRF_ROOK and handle error.
+*
+                  IF( INFO.NE.K)
+     $               CALL ALAERH( PATH, 'CHETRF_ROOK', INFO, K,
+     $                            UPLO, N, N, -1, -1, NB, IMAT,
+     $                            NFAIL, NERRS, NOUT )
+*
+*                 Set the condition estimate flag if the INFO is not 0.
+*
+                  IF( INFO.NE.0 ) THEN
+                     TRFCON = .TRUE.
+                  ELSE
+                     TRFCON = .FALSE.
+                  END IF
+*
+*+    TEST 1
+*                 Reconstruct matrix from factors and compute residual.
+*
+                  CALL CHET01_ROOK( UPLO, N, A, LDA, AFAC, LDA, IWORK,
+     $                              AINV, LDA, RWORK, RESULT( 1 ) )
+                  NT = 1
+*
+*+    TEST 2
+*                 Form the inverse and compute the residual,
+*                 if the factorization was competed without INFO > 0
+*                 (i.e. there is no zero rows and columns).
+*                 Do it only for the first block size.
+*
+                  IF( INB.EQ.1 .AND. .NOT.TRFCON ) THEN
+                     CALL CLACPY( UPLO, N, N, AFAC, LDA, AINV, LDA )
+                     SRNAMT = 'CHETRI_ROOK'
+                     CALL CHETRI_ROOK( UPLO, N, AINV, LDA, IWORK, WORK,
+     $                                 INFO )
+*
+*                    Check error code from CHETRI_ROOK and handle error.
+*
+                     IF( INFO.NE.0 )
+     $                  CALL ALAERH( PATH, 'CHETRI_ROOK', INFO, -1,
+     $                               UPLO, N, N, -1, -1, -1, IMAT,
+     $                               NFAIL, NERRS, NOUT )
+*
+*                    Compute the residual for a Hermitian matrix times
+*                    its inverse.
+*
+                     CALL CPOT03( UPLO, N, A, LDA, AINV, LDA, WORK, LDA,
+     $                            RWORK, RCONDC, RESULT( 2 ) )
+                     NT = 2
+                  END IF
+*
+*                 Print information about the tests that did not pass
+*                 the threshold.
+*
+                  DO 110 K = 1, NT
+                     IF( RESULT( K ).GE.THRESH ) THEN
+                        IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
+     $                     CALL ALAHD( NOUT, PATH )
+                        WRITE( NOUT, FMT = 9999 )UPLO, N, NB, IMAT, K,
+     $                     RESULT( K )
+                        NFAIL = NFAIL + 1
+                     END IF
+  110             CONTINUE
+                  NRUN = NRUN + NT
+*
+*+    TEST 3
+*                 Compute largest element in U or L
+*
+                  RESULT( 3 ) = ZERO
+                  STEMP = ZERO
+*
+                  CONST = ( ( ALPHA**2-ONE ) / ( ALPHA**2-ONEHALF ) ) /
+     $                    ( ONE-ALPHA )
+*
+                  IF( IUPLO.EQ.1 ) THEN
+*
+*                 Compute largest element in U
+*
+                     K = N
+  120                CONTINUE
+                     IF( K.LE.1 )
+     $                  GO TO 130
+*
+                     IF( IWORK( K ).GT.ZERO ) THEN
+*
+*                       Get max absolute value from elements
+*                       in column k in in U
+*
+                        STEMP = CLANGE( 'M', K-1, 1,
+     $                          AFAC( ( K-1 )*LDA+1 ), LDA, RWORK )
+                     ELSE
+*
+*                       Get max absolute value from elements
+*                       in columns k and k-1 in U
+*
+                        STEMP = CLANGE( 'M', K-2, 2,
+     $                          AFAC( ( K-2 )*LDA+1 ), LDA, RWORK )
+                        K = K - 1
+*
+                     END IF
+*
+*                    STEMP should be bounded by CONST
+*
+                     STEMP = STEMP - CONST + THRESH
+                     IF( STEMP.GT.RESULT( 3 ) )
+     $                  RESULT( 3 ) = STEMP
+*
+                     K = K - 1
+*
+                     GO TO 120
+  130                CONTINUE
+*
+                  ELSE
+*
+*                 Compute largest element in L
+*
+                     K = 1
+  140                CONTINUE
+                     IF( K.GE.N )
+     $                  GO TO 150
+*
+                     IF( IWORK( K ).GT.ZERO ) THEN
+*
+*                       Get max absolute value from elements
+*                       in column k in L
+*
+                        STEMP = CLANGE( 'M', N-K, 1,
+     $                          AFAC( ( K-1 )*LDA+K+1 ), LDA, RWORK )
+                     ELSE
+*
+*                       Get max absolute value from elements
+*                       in columns k and k+1 in L
+*
+                        STEMP = CLANGE( 'M', N-K-1, 2,
+     $                          AFAC( ( K-1 )*LDA+K+2 ), LDA, RWORK )
+                        K = K + 1
+*
+                     END IF
+*
+*                    STEMP should be bounded by CONST
+*
+                     STEMP = STEMP - CONST + THRESH
+                     IF( STEMP.GT.RESULT( 3 ) )
+     $                  RESULT( 3 ) = STEMP
+*
+                     K = K + 1
+*
+                     GO TO 140
+  150                CONTINUE
+                  END IF
+*
+*
+*+    TEST 4
+*                 Compute largest 2-Norm of 2-by-2 diag blocks
+*
+                  RESULT( 4 ) = ZERO
+                  STEMP = ZERO
+*
+                  CONST = ( ( ALPHA**2-ONE ) / ( ALPHA**2-ONEHALF ) )*
+     $                    ( ( ONE + ALPHA ) / ( ONE - ALPHA ) )
+*
+                  IF( IUPLO.EQ.1 ) THEN
+*
+*                    Loop backward for UPLO = 'U'
+*
+                     K = N
+  160                CONTINUE
+                     IF( K.LE.1 )
+     $                  GO TO 170
+*
+                     IF( IWORK( K ).LT.ZERO ) THEN
+*
+*                       Get the two eigenvalues of a 2-by-2 block,
+*                       store them in WORK array
+*
+                        BLOCK( 1, 1 ) = AFAC( ( K-2 )*LDA+K-1 )
+                        BLOCK( 2, 1 ) = AFAC( ( K-2 )*LDA+K )
+                        BLOCK( 1, 2 ) = BLOCK( 2, 1 )
+                        BLOCK( 2, 2 ) = AFAC( (K-1)*LDA+K )
+*
+                        CALL CHEEVX( 'N', 'N', 'N', 'N', 2, BLOCK,
+     $                               2, WORK, CDUMMY, 1, CDUMMY, 1,
+     $                               ITEMP, ITEMP2, RWORK, STEMP,
+     $                               RWORK( 3 ), RWORK( 5 ), WORK( 3 ),
+     $                               4, RWORK( 7 ), INFO )
+*
+                        LAM_MAX = MAX( ABS( WORK( 1 ) ),
+     $                            ABS( WORK( 2 ) ) )
+                        LAM_MIN = MIN( ABS( WORK( 1 ) ),
+     $                            ABS( WORK( 2 ) ) )
+*
+                        STEMP = LAM_MAX / LAM_MIN
+*
+*                       STEMP should be bounded by CONST
+*
+                        STEMP = ABS( STEMP ) - CONST + THRESH
+                        IF( STEMP.GT.RESULT( 4 ) )
+     $                     RESULT( 4 ) = STEMP
+                        K = K - 1
+*
+                     END IF
+*
+                     K = K - 1
+*
+                     GO TO 160
+  170                CONTINUE
+*
+                  ELSE
+*
+*                    Loop forward for UPLO = 'L'
+*
+                     K = 1
+  180                CONTINUE
+                     IF( K.GE.N )
+     $                  GO TO 190
+*
+                     IF( IWORK( K ).LT.ZERO ) THEN
+*
+*                       Get the two eigenvalues of a 2-by-2 block,
+*                       store them in WORK array
+*
+                        BLOCK( 1, 1 ) = AFAC( ( K-1 )*LDA+K )
+                        BLOCK( 2, 1 ) = AFAC( ( K-1 )*LDA+K+1 )
+                        BLOCK( 1, 2 ) = BLOCK( 2, 1 )
+                        BLOCK( 2, 2 ) = AFAC( K*LDA+K+1 )
+*
+                        CALL CHEEVX( 'N', 'N', 'N', 'N', 2, BLOCK,
+     $                               2, WORK, CDUMMY, 1, CDUMMY, 1,
+     $                               ITEMP, ITEMP2, RWORK, STEMP,
+     $                               RWORK( 3 ), RWORK( 5 ), WORK( 3 ),
+     $                               4, RWORK( 7 ), INFO )
+*
+                        LAM_MAX = MAX( ABS( WORK( 1 ) ),
+     $                            ABS( WORK( 2 ) ) )
+                        LAM_MIN = MIN( ABS( WORK( 1 ) ),
+     $                            ABS( WORK( 2 ) ) )
+*
+                        STEMP = LAM_MAX / LAM_MIN
+*
+*                       STEMP should be bounded by CONST
+*
+                        STEMP = ABS( STEMP ) - CONST + THRESH
+                        IF( STEMP.GT.RESULT( 4 ) )
+     $                     RESULT( 4 ) = STEMP
+                        K = K + 1
+*
+                     END IF
+*
+                     K = K + 1
+*
+                     GO TO 180
+  190                CONTINUE
+                  END IF
+*
+*                 Print information about the tests that did not pass
+*                 the threshold.
+*
+                  DO 200 K = 3, 4
+                     IF( RESULT( K ).GE.THRESH ) THEN
+                        IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
+     $                     CALL ALAHD( NOUT, PATH )
+                        WRITE( NOUT, FMT = 9999 )UPLO, N, NB, IMAT, K,
+     $                     RESULT( K )
+                        NFAIL = NFAIL + 1
+                     END IF
+  200             CONTINUE
+                  NRUN = NRUN + NT
+*
+*                 Skip the other tests if this is not the first block
+*                 size.
+*
+                  IF( INB.GT.1 )
+     $               GO TO 240
+*
+*                 Do only the condition estimate if INFO is not 0.
+*
+                  IF( TRFCON ) THEN
+                     RCONDC = ZERO
+                     GO TO 230
+                  END IF
+*
+                  DO 220 IRHS = 1, NNS
+                     NRHS = NSVAL( IRHS )
+*
+*                    Begin loop over NRHS values
+*
+*
+*+    TEST 5 ( Using TRS_ROOK)
+*                 Solve and compute residual for  A * X = B.
+*
+*                    Choose a set of NRHS random solution vectors
+*                    stored in XACT and set up the right hand side B
+*
+                     SRNAMT = 'CLARHS'
+                     CALL CLARHS( MATPATH, XTYPE, UPLO, ' ', N, N,
+     $                            KL, KU, NRHS, A, LDA, XACT, LDA,
+     $                            B, LDA, ISEED, INFO )
+                     CALL CLACPY( 'Full', N, NRHS, B, LDA, X, LDA )
+*
+                     SRNAMT = 'CHETRS_ROOK'
+                     CALL CHETRS_ROOK( UPLO, N, NRHS, AFAC, LDA, IWORK,
+     $                                 X, LDA, INFO )
+*
+*                    Check error code from CHETRS_ROOK and handle error.
+*
+                     IF( INFO.NE.0 )
+     $                  CALL ALAERH( PATH, 'CHETRS_ROOK', INFO, 0,
+     $                               UPLO, N, N, -1, -1, NRHS, IMAT,
+     $                               NFAIL, NERRS, NOUT )
+*
+                     CALL CLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA )
+*
+*                    Compute the residual for the solution
+*
+                     CALL CPOT02( UPLO, N, NRHS, A, LDA, X, LDA, WORK,
+     $                            LDA, RWORK, RESULT( 5 ) )
+*
+*+    TEST 6
+*                 Check solution from generated exact solution.
+*
+                     CALL CGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
+     $                            RESULT( 6 ) )
+*
+*                    Print information about the tests that did not pass
+*                    the threshold.
+*
+                     DO 210 K = 5, 6
+                        IF( RESULT( K ).GE.THRESH ) THEN
+                           IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
+     $                        CALL ALAHD( NOUT, PATH )
+                           WRITE( NOUT, FMT = 9998 )UPLO, N, NRHS,
+     $                        IMAT, K, RESULT( K )
+                           NFAIL = NFAIL + 1
+                        END IF
+  210                CONTINUE
+                     NRUN = NRUN + 2
+*
+*                    End loop over NRHS values
+*
+  220             CONTINUE
+*
+*+    TEST 7
+*                 Get an estimate of RCOND = 1/CNDNUM.
+*
+  230             CONTINUE
+                  ANORM = CLANHE( '1', UPLO, N, A, LDA, RWORK )
+                  SRNAMT = 'CHECON_ROOK'
+                  CALL CHECON_ROOK( UPLO, N, AFAC, LDA, IWORK, ANORM,
+     $                              RCOND, WORK, INFO )
+*
+*                 Check error code from CHECON_ROOK and handle error.
+*
+                  IF( INFO.NE.0 )
+     $               CALL ALAERH( PATH, 'CHECON_ROOK', INFO, 0,
+     $                             UPLO, N, N, -1, -1, -1, IMAT,
+     $                             NFAIL, NERRS, NOUT )
+*
+*                 Compute the test ratio to compare to values of RCOND
+*
+                  RESULT( 7 ) = SGET06( RCOND, RCONDC )
+*
+*                 Print information about the tests that did not pass
+*                 the threshold.
+*
+                  IF( RESULT( 7 ).GE.THRESH ) THEN
+                     IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
+     $                  CALL ALAHD( NOUT, PATH )
+                     WRITE( NOUT, FMT = 9997 )UPLO, N, IMAT, 7,
+     $                  RESULT( 7 )
+                     NFAIL = NFAIL + 1
+                  END IF
+                  NRUN = NRUN + 1
+  240          CONTINUE
+*
+  250       CONTINUE
+  260    CONTINUE
+  270 CONTINUE
+*
+*     Print a summary of the results.
+*
+      CALL ALASUM( PATH, NOUT, NFAIL, NRUN, NERRS )
+*
+ 9999 FORMAT( ' UPLO = ''', A1, ''', N =', I5, ', NB =', I4, ', type ',
+     $      I2, ', test ', I2, ', ratio =', G12.5 )
+ 9998 FORMAT( ' UPLO = ''', A1, ''', N =', I5, ', NRHS=', I3, ', type ',
+     $      I2, ', test(', I2, ') =', G12.5 )
+ 9997 FORMAT( ' UPLO = ''', A1, ''', N =', I5, ',', 10X, ' type ', I2,
+     $      ', test(', I2, ') =', G12.5 )
+      RETURN
+*
+*     End of CCHKHE_ROOK
+*
+      END
diff --git a/TESTING/LIN/chet01_rook.f b/TESTING/LIN/chet01_rook.f
new file mode 100644
index 00000000..b40f2fb1
--- /dev/null
+++ b/TESTING/LIN/chet01_rook.f
@@ -0,0 +1,239 @@
+*> \brief \b CHET01_ROOK
+*
+*  =========== DOCUMENTATION ===========
+*
+* Online html documentation available at 
+*            http://www.netlib.org/lapack/explore-html/ 
+*
+*  Definition:
+*  ===========
+*
+*       SUBROUTINE CHET01_ROOK( UPLO, N, A, LDA, AFAC, LDAFAC, IPIV, C, LDC,
+*                               RWORK, RESID )
+* 
+*       .. Scalar Arguments ..
+*       CHARACTER          UPLO
+*       INTEGER            LDA, LDAFAC, LDC, N
+*       REAL               RESID
+*       ..
+*       .. Array Arguments ..
+*       INTEGER            IPIV( * )
+*       REAL               RWORK( * )
+*       COMPLEX            A( LDA, * ), AFAC( LDAFAC, * ), C( LDC, * )
+*       ..
+*  
+*
+*> \par Purpose:
+*  =============
+*>
+*> \verbatim
+*>
+*> CHET01_ROOK reconstructs a complex Hermitian indefinite matrix A from its
+*> block L*D*L' or U*D*U' factorization and computes the residual
+*>    norm( C - A ) / ( N * norm(A) * EPS ),
+*> where C is the reconstructed matrix, EPS is the machine epsilon,
+*> L' is the transpose of L, and U' is the transpose of U.
+*> \endverbatim
+*
+*  Arguments:
+*  ==========
+*
+*> \param[in] UPLO
+*> \verbatim
+*>          UPLO is CHARACTER*1
+*>          Specifies whether the upper or lower triangular part of the
+*>          complex Hermitian matrix A is stored:
+*>          = 'U':  Upper triangular
+*>          = 'L':  Lower triangular
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*>          N is INTEGER
+*>          The number of rows and columns of the matrix A.  N >= 0.
+*> \endverbatim
+*>
+*> \param[in] A
+*> \verbatim
+*>          A is COMPLEX array, dimension (LDA,N)
+*>          The original complex Hermitian matrix A.
+*> \endverbatim
+*>
+*> \param[in] LDA
+*> \verbatim
+*>          LDA is INTEGER
+*>          The leading dimension of the array A.  LDA >= max(1,N)
+*> \endverbatim
+*>
+*> \param[in] AFAC
+*> \verbatim
+*>          AFAC is COMPLEX array, dimension (LDAFAC,N)
+*>          The factored form of the matrix A.  AFAC contains the block
+*>          diagonal matrix D and the multipliers used to obtain the
+*>          factor L or U from the block L*D*L' or U*D*U' factorization
+*>          as computed by CSYTRF_ROOK.
+*> \endverbatim
+*>
+*> \param[in] LDAFAC
+*> \verbatim
+*>          LDAFAC is INTEGER
+*>          The leading dimension of the array AFAC.  LDAFAC >= max(1,N).
+*> \endverbatim
+*>
+*> \param[in] IPIV
+*> \verbatim
+*>          IPIV is INTEGER array, dimension (N)
+*>          The pivot indices from CSYTRF_ROOK.
+*> \endverbatim
+*>
+*> \param[out] C
+*> \verbatim
+*>          C is COMPLEX array, dimension (LDC,N)
+*> \endverbatim
+*>
+*> \param[in] LDC
+*> \verbatim
+*>          LDC is INTEGER
+*>          The leading dimension of the array C.  LDC >= max(1,N).
+*> \endverbatim
+*>
+*> \param[out] RWORK
+*> \verbatim
+*>          RWORK is REAL array, dimension (N)
+*> \endverbatim
+*>
+*> \param[out] RESID
+*> \verbatim
+*>          RESID is REAL
+*>          If UPLO = 'L', norm(L*D*L' - A) / ( N * norm(A) * EPS )
+*>          If UPLO = 'U', norm(U*D*U' - A) / ( N * norm(A) * EPS )
+*> \endverbatim
+*
+*  Authors:
+*  ========
+*
+*> \author Univ. of Tennessee 
+*> \author Univ. of California Berkeley 
+*> \author Univ. of Colorado Denver 
+*> \author NAG Ltd. 
+*
+*> \date April 2013
+*
+*> \ingroup complex_lin
+*
+*  =====================================================================
+      SUBROUTINE CHET01_ROOK( UPLO, N, A, LDA, AFAC, LDAFAC, IPIV, C,
+     $                        LDC, RWORK, RESID )
+*
+*  -- LAPACK test 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            LDA, LDAFAC, LDC, N
+      REAL               RESID
+*     ..
+*     .. Array Arguments ..
+      INTEGER            IPIV( * )
+      REAL               RWORK( * )
+      COMPLEX            A( LDA, * ), AFAC( LDAFAC, * ), C( LDC, * )
+*     ..
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      REAL               ZERO, ONE
+      PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
+      COMPLEX            CZERO, CONE
+      PARAMETER          ( CZERO = ( 0.0E+0, 0.0E+0 ),
+     $                   CONE = ( 1.0E+0, 0.0E+0 ) )
+*     ..
+*     .. Local Scalars ..
+      INTEGER            I, INFO, J
+      REAL               ANORM, EPS
+*     ..
+*     .. External Functions ..
+      LOGICAL            LSAME
+      REAL               CLANHE, SLAMCH
+      EXTERNAL           LSAME, CLANHE, SLAMCH
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           CLASET, CLAVHE_ROOK
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          AIMAG, REAL
+*     ..
+*     .. Executable Statements ..
+*
+*     Quick exit if N = 0.
+*
+      IF( N.LE.0 ) THEN
+         RESID = ZERO
+         RETURN
+      END IF
+*
+*     Determine EPS and the norm of A.
+*
+      EPS = SLAMCH( 'Epsilon' )
+      ANORM = CLANHE( '1', UPLO, N, A, LDA, RWORK )
+*
+*     Check the imaginary parts of the diagonal elements and return with
+*     an error code if any are nonzero.
+*
+      DO 10 J = 1, N
+         IF( AIMAG( AFAC( J, J ) ).NE.ZERO ) THEN
+            RESID = ONE / EPS
+            RETURN
+         END IF
+   10 CONTINUE
+*
+*     Initialize C to the identity matrix.
+*
+      CALL CLASET( 'Full', N, N, CZERO, CONE, C, LDC )
+*
+*     Call CLAVHE_ROOK to form the product D * U' (or D * L' ).
+*
+      CALL CLAVHE_ROOK( UPLO, 'Conjugate', 'Non-unit', N, N, AFAC,
+     $                  LDAFAC, IPIV, C, LDC, INFO )
+*
+*     Call CLAVHE_ROOK again to multiply by U (or L ).
+*
+      CALL CLAVHE_ROOK( UPLO, 'No transpose', 'Unit', N, N, AFAC,
+     $                  LDAFAC, IPIV, C, LDC, INFO )
+*
+*     Compute the difference  C - A .
+*
+      IF( LSAME( UPLO, 'U' ) ) THEN
+         DO 30 J = 1, N
+            DO 20 I = 1, J - 1
+               C( I, J ) = C( I, J ) - A( I, J )
+   20       CONTINUE
+            C( J, J ) = C( J, J ) - REAL( A( J, J ) )
+   30    CONTINUE
+      ELSE
+         DO 50 J = 1, N
+            C( J, J ) = C( J, J ) - REAL( A( J, J ) )
+            DO 40 I = J + 1, N
+               C( I, J ) = C( I, J ) - A( I, J )
+   40       CONTINUE
+   50    CONTINUE
+      END IF
+*
+*     Compute norm( C - A ) / ( N * norm(A) * EPS )
+*
+      RESID = CLANHE( '1', UPLO, N, C, LDC, RWORK )
+*
+      IF( ANORM.LE.ZERO ) THEN
+         IF( RESID.NE.ZERO )
+     $      RESID = ONE / EPS
+      ELSE
+         RESID = ( ( RESID/REAL( N ) )/ANORM ) / EPS
+      END IF
+*
+      RETURN
+*
+*     End of CHET01_ROOK
+*
+      END
diff --git a/TESTING/LIN/clavhe_rook.f b/TESTING/LIN/clavhe_rook.f
new file mode 100644
index 00000000..4946a5e2
--- /dev/null
+++ b/TESTING/LIN/clavhe_rook.f
@@ -0,0 +1,614 @@
+*> \brief \b CLAVHE_ROOK
+*
+*  =========== DOCUMENTATION ===========
+*
+* Online html documentation available at 
+*            http://www.netlib.org/lapack/explore-html/ 
+*
+*  Definition:
+*  ===========
+*
+*       SUBROUTINE CLAVHE_ROOK( UPLO, TRANS, DIAG, N, NRHS, A, LDA, IPIV, B,
+*                               LDB, INFO )
+* 
+*       .. Scalar Arguments ..
+*       CHARACTER          DIAG, TRANS, UPLO
+*       INTEGER            INFO, LDA, LDB, N, NRHS
+*       ..
+*       .. Array Arguments ..
+*       INTEGER            IPIV( * )
+*       COMPLEX            A( LDA, * ), B( LDB, * )
+*       ..
+*  
+*
+*> \par Purpose:
+*  =============
+*>
+*> \verbatim
+*>
+*> CLAVHE_ROOK  performs one of the matrix-vector operations
+*>    x := A*x  or  x := A^H*x,
+*> where x is an N element vector and  A is one of the factors
+*> from the Hermitian factorization computed by CHETRF_ROOK.
+*>
+*> CHETRF_ROOK produces a factorization of the form
+*>      U * D * U^H     or     L * D * L^H,
+*> where U (or L) is a product of permutation and unit upper (lower)
+*> triangular matrices, U^H (or L^H) is the conjugate transpose of
+*> U (or L), and D is Hermitian and block diagonal with 1 x 1 and
+*> 2 x 2 diagonal blocks.  The multipliers for the transformations
+*> and the upper or lower triangular parts of the diagonal blocks
+*> are stored in the leading upper or lower triangle of the 2-D
+*> array A.
+*>
+*> If TRANS = 'N', multiplies by U  or U * D  (or L  or L * D)
+*> If TRANS = 'T', multiplies by U' or D * U' (or L' or D * L')
+*> If TRANS = 'C', multiplies by U' or D * U' (or L' or D * L')
+*> \endverbatim
+*
+*  Arguments:
+*  ==========
+*
+*> \param[in] UPLO
+*> \verbatim
+*>          UPLO is CHARACTER*1
+*>          Specifies whether the factor stored in A is upper or lower
+*>          triangular.
+*>          = 'U':  Upper triangular
+*>          = 'L':  Lower triangular
+*> \endverbatim
+*>
+*> \param[in] TRANS
+*> \verbatim
+*>          TRANS is CHARACTER*1
+*>          Specifies the operation to be performed:
+*>          = 'N':  x := A*x
+*>          = 'T':   x := A^H*x
+*>          = 'C':   x := A^H*x
+*> \endverbatim
+*>
+*> \param[in] DIAG
+*> \verbatim
+*>          DIAG is CHARACTER*1
+*>          Specifies whether or not the diagonal blocks are unit
+*>          matrices.  If the diagonal blocks are assumed to be unit,
+*>          then A = U or A = L, otherwise A = U*D or A = L*D.
+*>          = 'U':  Diagonal blocks are assumed to be unit matrices.
+*>          = 'N':  Diagonal blocks are assumed to be non-unit matrices.
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*>          N is INTEGER
+*>          The number of rows and columns 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 vectors
+*>          x to be multiplied by A.  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 CHETRF_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 CHETRF_ROOK.
+*>          If UPLO = 'U':
+*>             Only the last KB elements of IPIV are set.
+*>
+*>             If IPIV(k) > 0, then rows and columns k and IPIV(k) were
+*>             interchanged and D(k,k) is a 1-by-1 diagonal block.
+*>
+*>             If IPIV(k) < 0 and IPIV(k-1) < 0, then rows and
+*>             columns k and -IPIV(k) were interchanged and rows and
+*>             columns k-1 and -IPIV(k-1) were inerchaged,
+*>             D(k-1:k,k-1:k) is a 2-by-2 diagonal block.
+*>
+*>          If UPLO = 'L':
+*>             Only the first KB elements of IPIV are set.
+*>
+*>             If IPIV(k) > 0, then rows and columns k and IPIV(k)
+*>             were interchanged and D(k,k) is a 1-by-1 diagonal block.
+*>
+*>             If IPIV(k) < 0 and IPIV(k+1) < 0, then rows and
+*>             columns k and -IPIV(k) were interchanged and rows and
+*>             columns k+1 and -IPIV(k+1) were inerchaged,
+*>             D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
+*> \endverbatim
+*>
+*> \param[in,out] B
+*> \verbatim
+*>          B is COMPLEX array, dimension (LDB,NRHS)
+*>          On entry, B contains NRHS vectors of length N.
+*>          On exit, B is overwritten with the product A * B.
+*> \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 = -k, the k-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 April 2013
+*
+*> \ingroup complex_lin
+*
+*  =====================================================================
+      SUBROUTINE CLAVHE_ROOK( UPLO, TRANS, DIAG, N, NRHS, A, LDA, IPIV,
+     $                        B, LDB, INFO )
+*
+*  -- LAPACK test 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          DIAG, TRANS, 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            NOUNIT
+      INTEGER            J, K, KP
+      COMPLEX            D11, D12, D21, D22, T1, T2
+*     ..
+*     .. External Functions ..
+      LOGICAL            LSAME
+      EXTERNAL           LSAME
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           CGEMV, CGERU, CLACGV, CSCAL, CSWAP, XERBLA
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          ABS, CONJG, MAX
+*     ..
+*     .. Executable Statements ..
+*
+*     Test the input parameters.
+*
+      INFO = 0
+      IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
+         INFO = -1
+      ELSE IF( .NOT.LSAME( TRANS, 'N' ) .AND. .NOT.LSAME( TRANS, 'C' ) )
+     $          THEN
+         INFO = -2
+      ELSE IF( .NOT.LSAME( DIAG, 'U' ) .AND. .NOT.LSAME( DIAG, 'N' ) )
+     $          THEN
+         INFO = -3
+      ELSE IF( N.LT.0 ) THEN
+         INFO = -4
+      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
+         INFO = -6
+      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
+         INFO = -9
+      END IF
+      IF( INFO.NE.0 ) THEN
+         CALL XERBLA( 'CLAVHE_ROOK ', -INFO )
+         RETURN
+      END IF
+*
+*     Quick return if possible.
+*
+      IF( N.EQ.0 )
+     $   RETURN
+*
+      NOUNIT = LSAME( DIAG, 'N' )
+*------------------------------------------
+*
+*     Compute  B := A * B  (No transpose)
+*
+*------------------------------------------
+      IF( LSAME( TRANS, 'N' ) ) THEN
+*
+*        Compute  B := U*B
+*        where U = P(m)*inv(U(m))* ... *P(1)*inv(U(1))
+*
+         IF( LSAME( UPLO, 'U' ) ) THEN
+*
+*        Loop forward applying the transformations.
+*
+            K = 1
+   10       CONTINUE
+            IF( K.GT.N )
+     $         GO TO 30
+            IF( IPIV( K ).GT.0 ) THEN
+*
+*              1 x 1 pivot block
+*
+*              Multiply by the diagonal element if forming U * D.
+*
+               IF( NOUNIT )
+     $            CALL CSCAL( NRHS, A( K, K ), B( K, 1 ), LDB )
+*
+*              Multiply by  P(K) * inv(U(K))  if K > 1.
+*
+               IF( K.GT.1 ) THEN
+*
+*                 Apply the transformation.
+*
+                  CALL CGERU( K-1, NRHS, CONE, A( 1, K ), 1, B( K, 1 ),
+     $                        LDB, B( 1, 1 ), LDB )
+*
+*                 Interchange if P(K) != I.
+*
+                  KP = IPIV( K )
+                  IF( KP.NE.K )
+     $               CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+               END IF
+               K = K + 1
+            ELSE
+*
+*              2 x 2 pivot block
+*
+*              Multiply by the diagonal block if forming U * D.
+*
+               IF( NOUNIT ) THEN
+                  D11 = A( K, K )
+                  D22 = A( K+1, K+1 )
+                  D12 = A( K, K+1 )
+                  D21 = CONJG( D12 )
+                  DO 20 J = 1, NRHS
+                     T1 = B( K, J )
+                     T2 = B( K+1, J )
+                     B( K, J ) = D11*T1 + D12*T2
+                     B( K+1, J ) = D21*T1 + D22*T2
+   20             CONTINUE
+               END IF
+*
+*              Multiply by  P(K) * inv(U(K))  if K > 1.
+*
+               IF( K.GT.1 ) THEN
+*
+*                 Apply the transformations.
+*
+                  CALL CGERU( K-1, NRHS, CONE, A( 1, K ), 1, B( K, 1 ),
+     $                        LDB, B( 1, 1 ), LDB )
+                  CALL CGERU( K-1, NRHS, CONE, A( 1, K+1 ), 1,
+     $                        B( K+1, 1 ), LDB, B( 1, 1 ), LDB )
+*
+*                 Interchange if a permutation was applied at the
+*                 K-th step of the factorization.
+*
+*                 Swap the first of pair with IMAXth
+*
+                  KP = ABS( IPIV( K ) )
+                  IF( KP.NE.K )
+     $               CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+*
+*                 NOW swap the first of pair with Pth
+*
+                  KP = ABS( IPIV( K+1 ) )
+                  IF( KP.NE.K+1 )
+     $               CALL CSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ),
+     $                           LDB )
+               END IF
+               K = K + 2
+            END IF
+            GO TO 10
+   30       CONTINUE
+*
+*        Compute  B := L*B
+*        where L = P(1)*inv(L(1))* ... *P(m)*inv(L(m)) .
+*
+         ELSE
+*
+*           Loop backward applying the transformations to B.
+*
+            K = N
+   40       CONTINUE
+            IF( K.LT.1 )
+     $         GO TO 60
+*
+*           Test the pivot index.  If greater than zero, a 1 x 1
+*           pivot was used, otherwise a 2 x 2 pivot was used.
+*
+            IF( IPIV( K ).GT.0 ) THEN
+*
+*              1 x 1 pivot block:
+*
+*              Multiply by the diagonal element if forming L * D.
+*
+               IF( NOUNIT )
+     $            CALL CSCAL( NRHS, A( K, K ), B( K, 1 ), LDB )
+*
+*              Multiply by  P(K) * inv(L(K))  if K < N.
+*
+               IF( K.NE.N ) THEN
+                  KP = IPIV( K )
+*
+*                 Apply the transformation.
+*
+                  CALL CGERU( N-K, NRHS, CONE, A( K+1, K ), 1,
+     $                        B( K, 1 ), LDB, B( K+1, 1 ), LDB )
+*
+*                 Interchange if a permutation was applied at the
+*                 K-th step of the factorization.
+*
+                  IF( KP.NE.K )
+     $               CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+               END IF
+               K = K - 1
+*
+            ELSE
+*
+*              2 x 2 pivot block:
+*
+*              Multiply by the diagonal block if forming L * D.
+*
+               IF( NOUNIT ) THEN
+                  D11 = A( K-1, K-1 )
+                  D22 = A( K, K )
+                  D21 = A( K, K-1 )
+                  D12 = CONJG( D21 )
+                  DO 50 J = 1, NRHS
+                     T1 = B( K-1, J )
+                     T2 = B( K, J )
+                     B( K-1, J ) = D11*T1 + D12*T2
+                     B( K, J ) = D21*T1 + D22*T2
+   50             CONTINUE
+               END IF
+*
+*              Multiply by  P(K) * inv(L(K))  if K < N.
+*
+               IF( K.NE.N ) THEN
+*
+*                 Apply the transformation.
+*
+                  CALL CGERU( N-K, NRHS, CONE, A( K+1, K ), 1,
+     $                        B( K, 1 ), LDB, B( K+1, 1 ), LDB )
+                  CALL CGERU( N-K, NRHS, CONE, A( K+1, K-1 ), 1,
+     $                        B( K-1, 1 ), LDB, B( K+1, 1 ), LDB )
+*
+*                 Interchange if a permutation was applied at the
+*                 K-th step of the factorization.
+*
+*
+*                 Swap the second of pair with IMAXth
+*
+                  KP = ABS( IPIV( K ) )
+                  IF( KP.NE.K )
+     $               CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+*
+*                 NOW swap the first of pair with Pth
+*
+                  KP = ABS( IPIV( K-1 ) )
+                  IF( KP.NE.K-1 )
+     $               CALL CSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ),
+     $                           LDB )
+*
+               END IF
+               K = K - 2
+            END IF
+            GO TO 40
+   60       CONTINUE
+         END IF
+*--------------------------------------------------
+*
+*     Compute  B := A^H * B  (conjugate transpose)
+*
+*--------------------------------------------------
+      ELSE
+*
+*        Form  B := U^H*B
+*        where U  = P(m)*inv(U(m))* ... *P(1)*inv(U(1))
+*        and   U^H = inv(U^H(1))*P(1)* ... *inv(U^H(m))*P(m)
+*
+         IF( LSAME( UPLO, 'U' ) ) THEN
+*
+*           Loop backward applying the transformations.
+*
+            K = N
+   70       IF( K.LT.1 )
+     $         GO TO 90
+*
+*           1 x 1 pivot block.
+*
+            IF( IPIV( K ).GT.0 ) THEN
+               IF( K.GT.1 ) THEN
+*
+*                 Interchange if P(K) != I.
+*
+                  KP = IPIV( K )
+                  IF( KP.NE.K )
+     $               CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+*
+*                 Apply the transformation
+*                    y = y - B' conjg(x),
+*                 where x is a column of A and y is a row of B.
+*
+                  CALL CLACGV( NRHS, B( K, 1 ), LDB )
+                  CALL CGEMV( 'Conjugate', K-1, NRHS, CONE, B, LDB,
+     $                        A( 1, K ), 1, CONE, B( K, 1 ), LDB )
+                  CALL CLACGV( NRHS, B( K, 1 ), LDB )
+               END IF
+               IF( NOUNIT )
+     $            CALL CSCAL( NRHS, A( K, K ), B( K, 1 ), LDB )
+               K = K - 1
+*
+*           2 x 2 pivot block.
+*
+            ELSE
+               IF( K.GT.2 ) THEN
+*
+*                 Swap the second of pair with Pth
+*
+                  KP = ABS( IPIV( K ) )
+                  IF( KP.NE.K )
+     $               CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+*
+*                 Now swap the first of pair with IMAX(r)th
+*
+                  KP = ABS( IPIV( K-1 ) )
+                  IF( KP.NE.K-1 )
+     $               CALL CSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ),
+     $                           LDB )
+*
+*                 Apply the transformations
+*                    y = y - B' conjg(x),
+*                 where x is a block column of A and y is a block
+*                 row of B.
+*
+                  CALL CLACGV( NRHS, B( K, 1 ), LDB )
+                  CALL CGEMV( 'Conjugate', K-2, NRHS, CONE, B, LDB,
+     $                        A( 1, K ), 1, CONE, B( K, 1 ), LDB )
+                  CALL CLACGV( NRHS, B( K, 1 ), LDB )
+*
+                  CALL CLACGV( NRHS, B( K-1, 1 ), LDB )
+                  CALL CGEMV( 'Conjugate', K-2, NRHS, CONE, B, LDB,
+     $                        A( 1, K-1 ), 1, CONE, B( K-1, 1 ), LDB )
+                  CALL CLACGV( NRHS, B( K-1, 1 ), LDB )
+               END IF
+*
+*              Multiply by the diagonal block if non-unit.
+*
+               IF( NOUNIT ) THEN
+                  D11 = A( K-1, K-1 )
+                  D22 = A( K, K )
+                  D12 = A( K-1, K )
+                  D21 = CONJG( D12 )
+                  DO 80 J = 1, NRHS
+                     T1 = B( K-1, J )
+                     T2 = B( K, J )
+                     B( K-1, J ) = D11*T1 + D12*T2
+                     B( K, J ) = D21*T1 + D22*T2
+   80             CONTINUE
+               END IF
+               K = K - 2
+            END IF
+            GO TO 70
+   90       CONTINUE
+*
+*        Form  B := L^H*B
+*        where L  = P(1)*inv(L(1))* ... *P(m)*inv(L(m))
+*        and   L^H = inv(L^H(m))*P(m)* ... *inv(L^H(1))*P(1)
+*
+         ELSE
+*
+*           Loop forward applying the L-transformations.
+*
+            K = 1
+  100       CONTINUE
+            IF( K.GT.N )
+     $         GO TO 120
+*
+*           1 x 1 pivot block
+*
+            IF( IPIV( K ).GT.0 ) THEN
+               IF( K.LT.N ) THEN
+*
+*                 Interchange if P(K) != I.
+*
+                  KP = IPIV( K )
+                  IF( KP.NE.K )
+     $               CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+*
+*                 Apply the transformation
+*
+                  CALL CLACGV( NRHS, B( K, 1 ), LDB )
+                  CALL CGEMV( 'Conjugate', N-K, NRHS, CONE, B( K+1, 1 ),
+     $                       LDB, A( K+1, K ), 1, CONE, B( K, 1 ), LDB )
+                  CALL CLACGV( NRHS, B( K, 1 ), LDB )
+               END IF
+               IF( NOUNIT )
+     $            CALL CSCAL( NRHS, A( K, K ), B( K, 1 ), LDB )
+               K = K + 1
+*
+*           2 x 2 pivot block.
+*
+            ELSE
+               IF( K.LT.N-1 ) THEN
+*
+*                 Swap the first of pair with Pth
+*
+                  KP = ABS( IPIV( K ) )
+                  IF( KP.NE.K )
+     $               CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+*
+*                 Now swap the second of pair with IMAX(r)th
+*
+                  KP = ABS( IPIV( K+1 ) )
+                  IF( KP.NE.K+1 )
+     $               CALL CSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ),
+     $                           LDB )
+*
+*                 Apply the transformation
+*
+                  CALL CLACGV( NRHS, B( K+1, 1 ), LDB )
+                  CALL CGEMV( 'Conjugate', N-K-1, NRHS, CONE,
+     $                        B( K+2, 1 ), LDB, A( K+2, K+1 ), 1, CONE,
+     $                        B( K+1, 1 ), LDB )
+                  CALL CLACGV( NRHS, B( K+1, 1 ), LDB )
+*
+                  CALL CLACGV( NRHS, B( K, 1 ), LDB )
+                  CALL CGEMV( 'Conjugate', N-K-1, NRHS, CONE,
+     $                        B( K+2, 1 ), LDB, A( K+2, K ), 1, CONE,
+     $                        B( K, 1 ), LDB )
+                  CALL CLACGV( NRHS, B( K, 1 ), LDB )
+               END IF
+*
+*              Multiply by the diagonal block if non-unit.
+*
+               IF( NOUNIT ) THEN
+                  D11 = A( K, K )
+                  D22 = A( K+1, K+1 )
+                  D21 = A( K+1, K )
+                  D12 = CONJG( D21 )
+                  DO 110 J = 1, NRHS
+                     T1 = B( K, J )
+                     T2 = B( K+1, J )
+                     B( K, J ) = D11*T1 + D12*T2
+                     B( K+1, J ) = D21*T1 + D22*T2
+  110             CONTINUE
+               END IF
+               K = K + 2
+            END IF
+            GO TO 100
+  120       CONTINUE
+         END IF
+*
+      END IF
+      RETURN
+*
+*     End of CLAVHE_ROOK
+*
+      END
diff --git a/TESTING/LIN/zchkhe_rook.f b/TESTING/LIN/zchkhe_rook.f
new file mode 100644
index 00000000..19385ddb
--- /dev/null
+++ b/TESTING/LIN/zchkhe_rook.f
@@ -0,0 +1,850 @@
+*> \brief \b ZCHKHE_ROOK
+*
+*  =========== DOCUMENTATION ===========
+*
+* Online html documentation available at 
+*            http://www.netlib.org/lapack/explore-html/ 
+*
+*  Definition:
+*  ===========
+*
+*       SUBROUTINE ZCHKHE_ROOK( DOTYPE, NN, NVAL, NNB, NBVAL, NNS, NSVAL,
+*                               THRESH, TSTERR, NMAX, A, AFAC, AINV, B, X,
+*                               XACT, WORK, RWORK, IWORK, NOUT )
+* 
+*       .. Scalar Arguments ..
+*       LOGICAL            TSTERR
+*       INTEGER            NMAX, NN, NNB, NNS, NOUT
+*       DOUBLE PRECISION   THRESH
+*       ..
+*       .. Array Arguments ..
+*       LOGICAL            DOTYPE( * )
+*       INTEGER            IWORK( * ), NBVAL( * ), NSVAL( * ), NVAL( * )
+*       DOUBLE PRECISION   RWORK( * )
+*       COMPLEX*16         A( * ), AFAC( * ), AINV( * ), B( * ),
+*      $                   WORK( * ), X( * ), XACT( * )
+*       ..
+*  
+*
+*> \par Purpose:
+*  =============
+*>
+*> \verbatim
+*>
+*> ZCHKHE_ROOK tests ZHETRF_ROOK, -TRI_ROOK, -TRS_ROOK,
+*> and -CON_ROOK.
+*> \endverbatim
+*
+*  Arguments:
+*  ==========
+*
+*> \param[in] DOTYPE
+*> \verbatim
+*>          DOTYPE is LOGICAL array, dimension (NTYPES)
+*>          The matrix types to be used for testing.  Matrices of type j
+*>          (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) =
+*>          .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used.
+*> \endverbatim
+*>
+*> \param[in] NN
+*> \verbatim
+*>          NN is INTEGER
+*>          The number of values of N contained in the vector NVAL.
+*> \endverbatim
+*>
+*> \param[in] NVAL
+*> \verbatim
+*>          NVAL is INTEGER array, dimension (NN)
+*>          The values of the matrix dimension N.
+*> \endverbatim
+*>
+*> \param[in] NNB
+*> \verbatim
+*>          NNB is INTEGER
+*>          The number of values of NB contained in the vector NBVAL.
+*> \endverbatim
+*>
+*> \param[in] NBVAL
+*> \verbatim
+*>          NBVAL is INTEGER array, dimension (NBVAL)
+*>          The values of the blocksize NB.
+*> \endverbatim
+*>
+*> \param[in] NNS
+*> \verbatim
+*>          NNS is INTEGER
+*>          The number of values of NRHS contained in the vector NSVAL.
+*> \endverbatim
+*>
+*> \param[in] NSVAL
+*> \verbatim
+*>          NSVAL is INTEGER array, dimension (NNS)
+*>          The values of the number of right hand sides NRHS.
+*> \endverbatim
+*>
+*> \param[in] THRESH
+*> \verbatim
+*>          THRESH is DOUBLE PRECISION
+*>          The threshold value for the test ratios.  A result is
+*>          included in the output file if RESULT >= THRESH.  To have
+*>          every test ratio printed, use THRESH = 0.
+*> \endverbatim
+*>
+*> \param[in] TSTERR
+*> \verbatim
+*>          TSTERR is LOGICAL
+*>          Flag that indicates whether error exits are to be tested.
+*> \endverbatim
+*>
+*> \param[in] NMAX
+*> \verbatim
+*>          NMAX is INTEGER
+*>          The maximum value permitted for N, used in dimensioning the
+*>          work arrays.
+*> \endverbatim
+*>
+*> \param[out] A
+*> \verbatim
+*>          A is COMPLEX*16 array, dimension (NMAX*NMAX)
+*> \endverbatim
+*>
+*> \param[out] AFAC
+*> \verbatim
+*>          AFAC is COMPLEX*16 array, dimension (NMAX*NMAX)
+*> \endverbatim
+*>
+*> \param[out] AINV
+*> \verbatim
+*>          AINV is COMPLEX*16 array, dimension (NMAX*NMAX)
+*> \endverbatim
+*>
+*> \param[out] B
+*> \verbatim
+*>          B is COMPLEX*16 array, dimension (NMAX*NSMAX)
+*>          where NSMAX is the largest entry in NSVAL.
+*> \endverbatim
+*>
+*> \param[out] X
+*> \verbatim
+*>          X is COMPLEX*16 array, dimension (NMAX*NSMAX)
+*> \endverbatim
+*>
+*> \param[out] XACT
+*> \verbatim
+*>          XACT is COMPLEX*16 array, dimension (NMAX*NSMAX)
+*> \endverbatim
+*>
+*> \param[out] WORK
+*> \verbatim
+*>          WORK is COMPLEX*16 array, dimension
+*>                      (NMAX*max(3,NSMAX))
+*> \endverbatim
+*>
+*> \param[out] RWORK
+*> \verbatim
+*>          RWORK is DOUBLE PRECISION array, dimension
+*>                      (max(NMAX,2*NSMAX))
+*> \endverbatim
+*>
+*> \param[out] IWORK
+*> \verbatim
+*>          IWORK is INTEGER array, dimension (2*NMAX)
+*> \endverbatim
+*>
+*> \param[in] NOUT
+*> \verbatim
+*>          NOUT is INTEGER
+*>          The unit number for output.
+*> \endverbatim
+*
+*  Authors:
+*  ========
+*
+*> \author Univ. of Tennessee 
+*> \author Univ. of California Berkeley 
+*> \author Univ. of Colorado Denver 
+*> \author NAG Ltd. 
+*
+*> \date April 2013
+*
+*> \ingroup complex16_lin
+*
+*  =====================================================================
+      SUBROUTINE ZCHKHE_ROOK( DOTYPE, NN, NVAL, NNB, NBVAL, NNS, NSVAL,
+     $                        THRESH, TSTERR, NMAX, A, AFAC, AINV, B, X,
+     $                        XACT, WORK, RWORK, IWORK, NOUT )
+*
+*  -- LAPACK test routine (version 3.4.1) --
+*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*     April 2012
+*
+*     .. Scalar Arguments ..
+      LOGICAL            TSTERR
+      INTEGER            NMAX, NN, NNB, NNS, NOUT
+      DOUBLE PRECISION   THRESH
+*     ..
+*     .. Array Arguments ..
+      LOGICAL            DOTYPE( * )
+      INTEGER            IWORK( * ), NBVAL( * ), NSVAL( * ), NVAL( * )
+      DOUBLE PRECISION   RWORK( * )
+      COMPLEX*16         A( * ), AFAC( * ), AINV( * ), B( * ),
+     $                   WORK( * ), X( * ), XACT( * )
+*     ..
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      DOUBLE PRECISION   ZERO, ONE
+      PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
+      DOUBLE PRECISION   ONEHALF
+      PARAMETER          ( ONEHALF = 0.5E+0 )
+      DOUBLE PRECISION   EIGHT, SEVTEN
+      PARAMETER          ( EIGHT = 8.0E+0, SEVTEN = 17.0E+0 )
+      COMPLEX            CZERO
+      PARAMETER          ( CZERO = ( 0.0E+0, 0.0E+0 )  )
+      INTEGER            NTYPES
+      PARAMETER          ( NTYPES = 10 )
+      INTEGER            NTESTS
+      PARAMETER          ( NTESTS = 7 )
+*     ..
+*     .. Local Scalars ..
+      LOGICAL            TRFCON, ZEROT
+      CHARACTER          DIST, TYPE, UPLO, XTYPE
+      CHARACTER*3        PATH, MATPATH
+      INTEGER            I, I1, I2, IMAT, IN, INB, INFO, IOFF, IRHS,
+     $                   ITEMP, ITEMP2, IUPLO, IZERO, J, K, KL, KU, LDA,
+     $                   LWORK, MODE, N, NB, NERRS, NFAIL, NIMAT, NRHS,
+     $                   NRUN, NT
+      DOUBLE PRECISION   ALPHA, ANORM, CNDNUM, CONST, LAM_MAX, LAM_MIN,
+     $                   RCOND, RCONDC, STEMP
+*     ..
+*     .. Local Arrays ..
+      CHARACTER          UPLOS( 2 )
+      INTEGER            ISEED( 4 ), ISEEDY( 4 )
+      DOUBLE PRECISION   RESULT( NTESTS )
+      COMPLEX*16         BLOCK( 2, 2 ), CDUMMY( 1 )
+*     ..
+*     .. External Functions ..
+      DOUBLE PRECISION   ZLANGE, ZLANHE, DGET06
+      EXTERNAL           ZLANGE, ZLANHE, DGET06
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           ALAERH, ALAHD, ALASUM, ZERRHE, ZHEEVX, ZGET04,
+     $                   ZLACPY, ZLARHS, ZLATB4, ZLATMS, ZPOT02,
+     $                   ZPOT03, ZHECON_ROOK, ZHET01_ROOK, ZHETRF_ROOK,
+     $                   ZHETRI_ROOK, ZHETRS_ROOK, XLAENV
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          ABS, MAX, MIN, SQRT
+*     ..
+*     .. Scalars in Common ..
+      LOGICAL            LERR, OK
+      CHARACTER*32       SRNAMT
+      INTEGER            INFOT, NUNIT
+*     ..
+*     .. Common blocks ..
+      COMMON             / INFOC / INFOT, NUNIT, OK, LERR
+      COMMON             / SRNAMC / SRNAMT
+*     ..
+*     .. Data statements ..
+      DATA               ISEEDY / 1988, 1989, 1990, 1991 /
+      DATA               UPLOS / 'U', 'L' /
+*     ..
+*     .. Executable Statements ..
+*
+*     Initialize constants and the random number seed.
+*
+      ALPHA = ( ONE+SQRT( SEVTEN ) ) / EIGHT
+*
+*     Test path
+*
+      PATH( 1: 1 ) = 'Zomplex precision'
+      PATH( 2: 3 ) = 'HR'
+*
+*     Path to generate matrices
+*
+      MATPATH( 1: 1 ) = 'Zomplex precision'
+      MATPATH( 2: 3 ) = 'HE'
+*
+      NRUN = 0
+      NFAIL = 0
+      NERRS = 0
+      DO 10 I = 1, 4
+         ISEED( I ) = ISEEDY( I )
+   10 CONTINUE
+*
+*     Test the error exits
+*
+      IF( TSTERR )
+     $   CALL ZERRHE( PATH, NOUT )
+      INFOT = 0
+*
+*     Set the minimum block size for which the block routine should
+*     be used, which will be later returned by ILAENV
+*
+      CALL XLAENV( 2, 2 )
+*
+*     Do for each value of N in NVAL
+*
+      DO 270 IN = 1, NN
+         N = NVAL( IN )
+         LDA = MAX( N, 1 )
+         XTYPE = 'N'
+         NIMAT = NTYPES
+         IF( N.LE.0 )
+     $      NIMAT = 1
+*
+         IZERO = 0
+*
+*        Do for each value of matrix type IMAT
+*
+         DO 260 IMAT = 1, NIMAT
+*
+*           Do the tests only if DOTYPE( IMAT ) is true.
+*
+            IF( .NOT.DOTYPE( IMAT ) )
+     $         GO TO 260
+*
+*           Skip types 3, 4, 5, or 6 if the matrix size is too small.
+*
+            ZEROT = IMAT.GE.3 .AND. IMAT.LE.6
+            IF( ZEROT .AND. N.LT.IMAT-2 )
+     $         GO TO 260
+*
+*           Do first for UPLO = 'U', then for UPLO = 'L'
+*
+            DO 250 IUPLO = 1, 2
+               UPLO = UPLOS( IUPLO )
+*
+*                 Begin generate the test matrix A.
+*
+*                 Set up parameters with ZLATB4 for the matrix generator
+*                 based on the type of matrix to be generated.
+*
+                  CALL ZLATB4( MATPATH, IMAT, N, N, TYPE, KL, KU, ANORM,
+     $                         MODE, CNDNUM, DIST )
+*
+*                 Generate a matrix with ZLATMS.
+*
+                  SRNAMT = 'ZLATMS'
+                  CALL ZLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE,
+     $                         CNDNUM, ANORM, KL, KU, UPLO, A, LDA,
+     $                         WORK, INFO )
+*
+*                 Check error code from ZLATMS and handle error.
+*
+                  IF( INFO.NE.0 ) THEN
+                     CALL ALAERH( PATH, 'ZLATMS', INFO, 0, UPLO, N, N,
+     $                            -1, -1, -1, IMAT, NFAIL, NERRS, NOUT )
+*
+*                    Skip all tests for this generated matrix
+*
+                     GO TO 250
+                  END IF
+*
+*                 For matrix types 3-6, zero one or more rows and
+*                 columns of the matrix to test that INFO is returned
+*                 correctly.
+*
+                  IF( ZEROT ) THEN
+                     IF( IMAT.EQ.3 ) THEN
+                        IZERO = 1
+                     ELSE IF( IMAT.EQ.4 ) THEN
+                        IZERO = N
+                     ELSE
+                        IZERO = N / 2 + 1
+                     END IF
+*
+                     IF( IMAT.LT.6 ) THEN
+*
+*                       Set row and column IZERO to zero.
+*
+                        IF( IUPLO.EQ.1 ) THEN
+                           IOFF = ( IZERO-1 )*LDA
+                           DO 20 I = 1, IZERO - 1
+                              A( IOFF+I ) = CZERO
+   20                      CONTINUE
+                           IOFF = IOFF + IZERO
+                           DO 30 I = IZERO, N
+                              A( IOFF ) = CZERO
+                              IOFF = IOFF + LDA
+   30                      CONTINUE
+                        ELSE
+                           IOFF = IZERO
+                           DO 40 I = 1, IZERO - 1
+                              A( IOFF ) = CZERO
+                              IOFF = IOFF + LDA
+   40                      CONTINUE
+                           IOFF = IOFF - IZERO
+                           DO 50 I = IZERO, N
+                              A( IOFF+I ) = CZERO
+   50                      CONTINUE
+                        END IF
+                     ELSE
+                        IOFF = 0
+                        IF( IUPLO.EQ.1 ) THEN
+*
+*                          Set the first IZERO rows and columns to zero.
+*
+                           DO 70 J = 1, N
+                              I2 = MIN( J, IZERO )
+                              DO 60 I = 1, I2
+                                 A( IOFF+I ) = CZERO
+   60                         CONTINUE
+                              IOFF = IOFF + LDA
+   70                      CONTINUE
+                        ELSE
+*
+*                          Set the last IZERO rows and columns to zero.
+*
+                           DO 90 J = 1, N
+                              I1 = MAX( J, IZERO )
+                              DO 80 I = I1, N
+                                 A( IOFF+I ) = CZERO
+   80                         CONTINUE
+                              IOFF = IOFF + LDA
+   90                      CONTINUE
+                        END IF
+                     END IF
+                  ELSE
+                     IZERO = 0
+                  END IF
+*
+*                 End generate the test matrix A.
+*
+*
+*              Do for each value of NB in NBVAL
+*
+               DO 240 INB = 1, NNB
+*
+*                 Set the optimal blocksize, which will be later
+*                 returned by ILAENV.
+*
+                  NB = NBVAL( INB )
+                  CALL XLAENV( 1, NB )
+*
+*                 Copy the test matrix A into matrix AFAC which
+*                 will be factorized in place. This is needed to
+*                 preserve the test matrix A for subsequent tests.
+*
+                  CALL ZLACPY( UPLO, N, N, A, LDA, AFAC, LDA )
+*
+*                 Compute the L*D*L**T or U*D*U**T factorization of the
+*                 matrix. IWORK stores details of the interchanges and
+*                 the block structure of D. AINV is a work array for
+*                 block factorization, LWORK is the length of AINV.
+*
+                  LWORK = MAX( 2, NB )*LDA
+                  SRNAMT = 'ZHETRF_ROOK'
+                  CALL ZHETRF_ROOK( UPLO, N, AFAC, LDA, IWORK, AINV,
+     $                              LWORK, INFO )
+*
+*                 Adjust the expected value of INFO to account for
+*                 pivoting.
+*
+                  K = IZERO
+                  IF( K.GT.0 ) THEN
+  100                CONTINUE
+                     IF( IWORK( K ).LT.0 ) THEN
+                        IF( IWORK( K ).NE.-K ) THEN
+                           K = -IWORK( K )
+                           GO TO 100
+                        END IF
+                     ELSE IF( IWORK( K ).NE.K ) THEN
+                        K = IWORK( K )
+                        GO TO 100
+                     END IF
+                  END IF
+*
+*                 Check error code from ZHETRF_ROOK and handle error.
+*
+                  IF( INFO.NE.K)
+     $               CALL ALAERH( PATH, 'ZHETRF_ROOK', INFO, K,
+     $                            UPLO, N, N, -1, -1, NB, IMAT,
+     $                            NFAIL, NERRS, NOUT )
+*
+*                 Set the condition estimate flag if the INFO is not 0.
+*
+                  IF( INFO.NE.0 ) THEN
+                     TRFCON = .TRUE.
+                  ELSE
+                     TRFCON = .FALSE.
+                  END IF
+*
+*+    TEST 1
+*                 Reconstruct matrix from factors and compute residual.
+*
+                  CALL ZHET01_ROOK( UPLO, N, A, LDA, AFAC, LDA, IWORK,
+     $                              AINV, LDA, RWORK, RESULT( 1 ) )
+                  NT = 1
+*
+*+    TEST 2
+*                 Form the inverse and compute the residual,
+*                 if the factorization was competed without INFO > 0
+*                 (i.e. there is no zero rows and columns).
+*                 Do it only for the first block size.
+*
+                  IF( INB.EQ.1 .AND. .NOT.TRFCON ) THEN
+                     CALL ZLACPY( UPLO, N, N, AFAC, LDA, AINV, LDA )
+                     SRNAMT = 'ZHETRI_ROOK'
+                     CALL ZHETRI_ROOK( UPLO, N, AINV, LDA, IWORK, WORK,
+     $                                 INFO )
+*
+*                    Check error code from ZHETRI_ROOK and handle error.
+*
+                     IF( INFO.NE.0 )
+     $                  CALL ALAERH( PATH, 'ZHETRI_ROOK', INFO, -1,
+     $                               UPLO, N, N, -1, -1, -1, IMAT,
+     $                               NFAIL, NERRS, NOUT )
+*
+*                    Compute the residual for a Hermitian matrix times
+*                    its inverse.
+*
+                     CALL ZPOT03( UPLO, N, A, LDA, AINV, LDA, WORK, LDA,
+     $                            RWORK, RCONDC, RESULT( 2 ) )
+                     NT = 2
+                  END IF
+*
+*                 Print information about the tests that did not pass
+*                 the threshold.
+*
+                  DO 110 K = 1, NT
+                     IF( RESULT( K ).GE.THRESH ) THEN
+                        IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
+     $                     CALL ALAHD( NOUT, PATH )
+                        WRITE( NOUT, FMT = 9999 )UPLO, N, NB, IMAT, K,
+     $                     RESULT( K )
+                        NFAIL = NFAIL + 1
+                     END IF
+  110             CONTINUE
+                  NRUN = NRUN + NT
+*
+*+    TEST 3
+*                 Compute largest element in U or L
+*
+                  RESULT( 3 ) = ZERO
+                  STEMP = ZERO
+*
+                  CONST = ( ( ALPHA**2-ONE ) / ( ALPHA**2-ONEHALF ) ) /
+     $                    ( ONE-ALPHA )
+*
+                  IF( IUPLO.EQ.1 ) THEN
+*
+*                 Compute largest element in U
+*
+                     K = N
+  120                CONTINUE
+                     IF( K.LE.1 )
+     $                  GO TO 130
+*
+                     IF( IWORK( K ).GT.ZERO ) THEN
+*
+*                       Get max absolute value from elements
+*                       in column k in in U
+*
+                        STEMP = ZLANGE( 'M', K-1, 1,
+     $                          AFAC( ( K-1 )*LDA+1 ), LDA, RWORK )
+                     ELSE
+*
+*                       Get max absolute value from elements
+*                       in columns k and k-1 in U
+*
+                        STEMP = ZLANGE( 'M', K-2, 2,
+     $                          AFAC( ( K-2 )*LDA+1 ), LDA, RWORK )
+                        K = K - 1
+*
+                     END IF
+*
+*                    STEMP should be bounded by CONST
+*
+                     STEMP = STEMP - CONST + THRESH
+                     IF( STEMP.GT.RESULT( 3 ) )
+     $                  RESULT( 3 ) = STEMP
+*
+                     K = K - 1
+*
+                     GO TO 120
+  130                CONTINUE
+*
+                  ELSE
+*
+*                 Compute largest element in L
+*
+                     K = 1
+  140                CONTINUE
+                     IF( K.GE.N )
+     $                  GO TO 150
+*
+                     IF( IWORK( K ).GT.ZERO ) THEN
+*
+*                       Get max absolute value from elements
+*                       in column k in L
+*
+                        STEMP = ZLANGE( 'M', N-K, 1,
+     $                          AFAC( ( K-1 )*LDA+K+1 ), LDA, RWORK )
+                     ELSE
+*
+*                       Get max absolute value from elements
+*                       in columns k and k+1 in L
+*
+                        STEMP = ZLANGE( 'M', N-K-1, 2,
+     $                          AFAC( ( K-1 )*LDA+K+2 ), LDA, RWORK )
+                        K = K + 1
+*
+                     END IF
+*
+*                    STEMP should be bounded by CONST
+*
+                     STEMP = STEMP - CONST + THRESH
+                     IF( STEMP.GT.RESULT( 3 ) )
+     $                  RESULT( 3 ) = STEMP
+*
+                     K = K + 1
+*
+                     GO TO 140
+  150                CONTINUE
+                  END IF
+*
+*
+*+    TEST 4
+*                 Compute largest 2-Norm of 2-by-2 diag blocks
+*
+                  RESULT( 4 ) = ZERO
+                  STEMP = ZERO
+*
+                  CONST = ( ( ALPHA**2-ONE ) / ( ALPHA**2-ONEHALF ) )*
+     $                    ( ( ONE + ALPHA ) / ( ONE - ALPHA ) )
+*
+                  IF( IUPLO.EQ.1 ) THEN
+*
+*                    Loop backward for UPLO = 'U'
+*
+                     K = N
+  160                CONTINUE
+                     IF( K.LE.1 )
+     $                  GO TO 170
+*
+                     IF( IWORK( K ).LT.ZERO ) THEN
+*
+*                       Get the two eigenvalues of a 2-by-2 block,
+*                       store them in WORK array
+*
+                        BLOCK( 1, 1 ) = AFAC( ( K-2 )*LDA+K-1 )
+                        BLOCK( 2, 1 ) = AFAC( ( K-2 )*LDA+K )
+                        BLOCK( 1, 2 ) = BLOCK( 2, 1 )
+                        BLOCK( 2, 2 ) = AFAC( (K-1)*LDA+K )
+*
+                        CALL ZHEEVX( 'N', 'N', 'N', 'N', 2, BLOCK,
+     $                               2, WORK, CDUMMY, 1, CDUMMY, 1,
+     $                               ITEMP, ITEMP2, RWORK, STEMP,
+     $                               RWORK( 3 ), RWORK( 5 ), WORK( 3 ),
+     $                               4, RWORK( 7 ), INFO )
+*
+                        LAM_MAX = MAX( ABS( WORK( 1 ) ),
+     $                            ABS( WORK( 2 ) ) )
+                        LAM_MIN = MIN( ABS( WORK( 1 ) ),
+     $                            ABS( WORK( 2 ) ) )
+*
+                        STEMP = LAM_MAX / LAM_MIN
+*
+*                       STEMP should be bounded by CONST
+*
+                        STEMP = ABS( STEMP ) - CONST + THRESH
+                        IF( STEMP.GT.RESULT( 4 ) )
+     $                     RESULT( 4 ) = STEMP
+                        K = K - 1
+*
+                     END IF
+*
+                     K = K - 1
+*
+                     GO TO 160
+  170                CONTINUE
+*
+                  ELSE
+*
+*                    Loop forward for UPLO = 'L'
+*
+                     K = 1
+  180                CONTINUE
+                     IF( K.GE.N )
+     $                  GO TO 190
+*
+                     IF( IWORK( K ).LT.ZERO ) THEN
+*
+*                       Get the two eigenvalues of a 2-by-2 block,
+*                       store them in WORK array
+*
+                        BLOCK( 1, 1 ) = AFAC( ( K-1 )*LDA+K )
+                        BLOCK( 2, 1 ) = AFAC( ( K-1 )*LDA+K+1 )
+                        BLOCK( 1, 2 ) = BLOCK( 2, 1 )
+                        BLOCK( 2, 2 ) = AFAC( K*LDA+K+1 )
+*
+                        CALL ZHEEVX( 'N', 'N', 'N', 'N', 2, BLOCK,
+     $                               2, WORK, CDUMMY, 1, CDUMMY, 1,
+     $                               ITEMP, ITEMP2, RWORK, STEMP,
+     $                               RWORK( 3 ), RWORK( 5 ), WORK( 3 ),
+     $                               4, RWORK( 7 ), INFO )
+*
+                        LAM_MAX = MAX( ABS( WORK( 1 ) ),
+     $                            ABS( WORK( 2 ) ) )
+                        LAM_MIN = MIN( ABS( WORK( 1 ) ),
+     $                            ABS( WORK( 2 ) ) )
+*
+                        STEMP = LAM_MAX / LAM_MIN
+*
+*                       STEMP should be bounded by CONST
+*
+                        STEMP = ABS( STEMP ) - CONST + THRESH
+                        IF( STEMP.GT.RESULT( 4 ) )
+     $                     RESULT( 4 ) = STEMP
+                        K = K + 1
+*
+                     END IF
+*
+                     K = K + 1
+*
+                     GO TO 180
+  190                CONTINUE
+                  END IF
+*
+*                 Print information about the tests that did not pass
+*                 the threshold.
+*
+                  DO 200 K = 3, 4
+                     IF( RESULT( K ).GE.THRESH ) THEN
+                        IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
+     $                     CALL ALAHD( NOUT, PATH )
+                        WRITE( NOUT, FMT = 9999 )UPLO, N, NB, IMAT, K,
+     $                     RESULT( K )
+                        NFAIL = NFAIL + 1
+                     END IF
+  200             CONTINUE
+                  NRUN = NRUN + NT
+*
+*                 Skip the other tests if this is not the first block
+*                 size.
+*
+                  IF( INB.GT.1 )
+     $               GO TO 240
+*
+*                 Do only the condition estimate if INFO is not 0.
+*
+                  IF( TRFCON ) THEN
+                     RCONDC = ZERO
+                     GO TO 230
+                  END IF
+*
+                  DO 220 IRHS = 1, NNS
+                     NRHS = NSVAL( IRHS )
+*
+*                    Begin loop over NRHS values
+*
+*
+*+    TEST 5 ( Using TRS_ROOK)
+*                 Solve and compute residual for  A * X = B.
+*
+*                    Choose a set of NRHS random solution vectors
+*                    stored in XACT and set up the right hand side B
+*
+                     SRNAMT = 'ZLARHS'
+                     CALL ZLARHS( MATPATH, XTYPE, UPLO, ' ', N, N,
+     $                            KL, KU, NRHS, A, LDA, XACT, LDA,
+     $                            B, LDA, ISEED, INFO )
+                     CALL ZLACPY( 'Full', N, NRHS, B, LDA, X, LDA )
+*
+                     SRNAMT = 'ZHETRS_ROOK'
+                     CALL ZHETRS_ROOK( UPLO, N, NRHS, AFAC, LDA, IWORK,
+     $                                 X, LDA, INFO )
+*
+*                    Check error code from ZHETRS_ROOK and handle error.
+*
+                     IF( INFO.NE.0 )
+     $                  CALL ALAERH( PATH, 'ZHETRS_ROOK', INFO, 0,
+     $                               UPLO, N, N, -1, -1, NRHS, IMAT,
+     $                               NFAIL, NERRS, NOUT )
+*
+                     CALL ZLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA )
+*
+*                    Compute the residual for the solution
+*
+                     CALL ZPOT02( UPLO, N, NRHS, A, LDA, X, LDA, WORK,
+     $                            LDA, RWORK, RESULT( 5 ) )
+*
+*+    TEST 6
+*                 Check solution from generated exact solution.
+*
+                     CALL ZGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
+     $                            RESULT( 6 ) )
+*
+*                    Print information about the tests that did not pass
+*                    the threshold.
+*
+                     DO 210 K = 5, 6
+                        IF( RESULT( K ).GE.THRESH ) THEN
+                           IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
+     $                        CALL ALAHD( NOUT, PATH )
+                           WRITE( NOUT, FMT = 9998 )UPLO, N, NRHS,
+     $                        IMAT, K, RESULT( K )
+                           NFAIL = NFAIL + 1
+                        END IF
+  210                CONTINUE
+                     NRUN = NRUN + 2
+*
+*                    End loop over NRHS values
+*
+  220             CONTINUE
+*
+*+    TEST 7
+*                 Get an estimate of RCOND = 1/CNDNUM.
+*
+  230             CONTINUE
+                  ANORM = ZLANHE( '1', UPLO, N, A, LDA, RWORK )
+                  SRNAMT = 'ZHECON_ROOK'
+                  CALL ZHECON_ROOK( UPLO, N, AFAC, LDA, IWORK, ANORM,
+     $                              RCOND, WORK, INFO )
+*
+*                 Check error code from ZHECON_ROOK and handle error.
+*
+                  IF( INFO.NE.0 )
+     $               CALL ALAERH( PATH, 'ZHECON_ROOK', INFO, 0,
+     $                             UPLO, N, N, -1, -1, -1, IMAT,
+     $                             NFAIL, NERRS, NOUT )
+*
+*                 Compute the test ratio to compare to values of RCOND
+*
+                  RESULT( 7 ) = DGET06( RCOND, RCONDC )
+*
+*                 Print information about the tests that did not pass
+*                 the threshold.
+*
+                  IF( RESULT( 7 ).GE.THRESH ) THEN
+                     IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
+     $                  CALL ALAHD( NOUT, PATH )
+                     WRITE( NOUT, FMT = 9997 )UPLO, N, IMAT, 7,
+     $                  RESULT( 7 )
+                     NFAIL = NFAIL + 1
+                  END IF
+                  NRUN = NRUN + 1
+  240          CONTINUE
+*
+  250       CONTINUE
+  260    CONTINUE
+  270 CONTINUE
+*
+*     Print a summary of the results.
+*
+      CALL ALASUM( PATH, NOUT, NFAIL, NRUN, NERRS )
+*
+ 9999 FORMAT( ' UPLO = ''', A1, ''', N =', I5, ', NB =', I4, ', type ',
+     $      I2, ', test ', I2, ', ratio =', G12.5 )
+ 9998 FORMAT( ' UPLO = ''', A1, ''', N =', I5, ', NRHS=', I3, ', type ',
+     $      I2, ', test(', I2, ') =', G12.5 )
+ 9997 FORMAT( ' UPLO = ''', A1, ''', N =', I5, ',', 10X, ' type ', I2,
+     $      ', test(', I2, ') =', G12.5 )
+      RETURN
+*
+*     End of ZCHKHE_ROOK
+*
+      END
diff --git a/TESTING/LIN/zhet01_rook.f b/TESTING/LIN/zhet01_rook.f
new file mode 100644
index 00000000..36041ab0
--- /dev/null
+++ b/TESTING/LIN/zhet01_rook.f
@@ -0,0 +1,239 @@
+*> \brief \b ZHET01_ROOK
+*
+*  =========== DOCUMENTATION ===========
+*
+* Online html documentation available at 
+*            http://www.netlib.org/lapack/explore-html/ 
+*
+*  Definition:
+*  ===========
+*
+*       SUBROUTINE ZHET01_ROOK( UPLO, N, A, LDA, AFAC, LDAFAC, IPIV, C, LDC,
+*                               RWORK, RESID )
+* 
+*       .. Scalar Arguments ..
+*       CHARACTER          UPLO
+*       INTEGER            LDA, LDAFAC, LDC, N
+*       DOUBLE PRECISION   RESID
+*       ..
+*       .. Array Arguments ..
+*       INTEGER            IPIV( * )
+*       DOUBLE PRECISION   RWORK( * )
+*       COMPLEX*16         A( LDA, * ), AFAC( LDAFAC, * ), C( LDC, * )
+*       ..
+*  
+*
+*> \par Purpose:
+*  =============
+*>
+*> \verbatim
+*>
+*> ZHET01_ROOK reconstructs a complex Hermitian indefinite matrix A from its
+*> block L*D*L' or U*D*U' factorization and computes the residual
+*>    norm( C - A ) / ( N * norm(A) * EPS ),
+*> where C is the reconstructed matrix, EPS is the machine epsilon,
+*> L' is the transpose of L, and U' is the transpose of U.
+*> \endverbatim
+*
+*  Arguments:
+*  ==========
+*
+*> \param[in] UPLO
+*> \verbatim
+*>          UPLO is CHARACTER*1
+*>          Specifies whether the upper or lower triangular part of the
+*>          complex Hermitian matrix A is stored:
+*>          = 'U':  Upper triangular
+*>          = 'L':  Lower triangular
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*>          N is INTEGER
+*>          The number of rows and columns of the matrix A.  N >= 0.
+*> \endverbatim
+*>
+*> \param[in] A
+*> \verbatim
+*>          A is COMPLEX*16 array, dimension (LDA,N)
+*>          The original complex Hermitian matrix A.
+*> \endverbatim
+*>
+*> \param[in] LDA
+*> \verbatim
+*>          LDA is INTEGER
+*>          The leading dimension of the array A.  LDA >= max(1,N)
+*> \endverbatim
+*>
+*> \param[in] AFAC
+*> \verbatim
+*>          AFAC is COMPLEX*16 array, dimension (LDAFAC,N)
+*>          The factored form of the matrix A.  AFAC contains the block
+*>          diagonal matrix D and the multipliers used to obtain the
+*>          factor L or U from the block L*D*L' or U*D*U' factorization
+*>          as computed by CSYTRF_ROOK.
+*> \endverbatim
+*>
+*> \param[in] LDAFAC
+*> \verbatim
+*>          LDAFAC is INTEGER
+*>          The leading dimension of the array AFAC.  LDAFAC >= max(1,N).
+*> \endverbatim
+*>
+*> \param[in] IPIV
+*> \verbatim
+*>          IPIV is INTEGER array, dimension (N)
+*>          The pivot indices from CSYTRF_ROOK.
+*> \endverbatim
+*>
+*> \param[out] C
+*> \verbatim
+*>          C is COMPLEX*16 array, dimension (LDC,N)
+*> \endverbatim
+*>
+*> \param[in] LDC
+*> \verbatim
+*>          LDC is INTEGER
+*>          The leading dimension of the array C.  LDC >= max(1,N).
+*> \endverbatim
+*>
+*> \param[out] RWORK
+*> \verbatim
+*>          RWORK is DOUBLE PRECISION array, dimension (N)
+*> \endverbatim
+*>
+*> \param[out] RESID
+*> \verbatim
+*>          RESID is DOUBLE PRECISION
+*>          If UPLO = 'L', norm(L*D*L' - A) / ( N * norm(A) * EPS )
+*>          If UPLO = 'U', norm(U*D*U' - A) / ( N * norm(A) * EPS )
+*> \endverbatim
+*
+*  Authors:
+*  ========
+*
+*> \author Univ. of Tennessee 
+*> \author Univ. of California Berkeley 
+*> \author Univ. of Colorado Denver 
+*> \author NAG Ltd. 
+*
+*> \date April 2013
+*
+*> \ingroup complex16_lin
+*
+*  =====================================================================
+      SUBROUTINE ZHET01_ROOK( UPLO, N, A, LDA, AFAC, LDAFAC, IPIV, C,
+     $                        LDC, RWORK, RESID )
+*
+*  -- LAPACK test 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            LDA, LDAFAC, LDC, N
+      DOUBLE PRECISION   RESID
+*     ..
+*     .. Array Arguments ..
+      INTEGER            IPIV( * )
+      DOUBLE PRECISION   RWORK( * )
+      COMPLEX*16         A( LDA, * ), AFAC( LDAFAC, * ), C( LDC, * )
+*     ..
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      DOUBLE PRECISION   ZERO, ONE
+      PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
+      COMPLEX*16         CZERO, CONE
+      PARAMETER          ( CZERO = ( 0.0E+0, 0.0E+0 ),
+     $                   CONE = ( 1.0E+0, 0.0E+0 ) )
+*     ..
+*     .. Local Scalars ..
+      INTEGER            I, INFO, J
+      DOUBLE PRECISION   ANORM, EPS
+*     ..
+*     .. External Functions ..
+      LOGICAL            LSAME
+      DOUBLE PRECISION   ZLANHE, DLAMCH
+      EXTERNAL           LSAME, ZLANHE, DLAMCH
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           ZLASET, ZLAVHE_ROOK
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          DIMAG, DBLE
+*     ..
+*     .. Executable Statements ..
+*
+*     Quick exit if N = 0.
+*
+      IF( N.LE.0 ) THEN
+         RESID = ZERO
+         RETURN
+      END IF
+*
+*     Determine EPS and the norm of A.
+*
+      EPS = DLAMCH( 'Epsilon' )
+      ANORM = ZLANHE( '1', UPLO, N, A, LDA, RWORK )
+*
+*     Check the imaginary parts of the diagonal elements and return with
+*     an error code if any are nonzero.
+*
+      DO 10 J = 1, N
+         IF( DIMAG( AFAC( J, J ) ).NE.ZERO ) THEN
+            RESID = ONE / EPS
+            RETURN
+         END IF
+   10 CONTINUE
+*
+*     Initialize C to the identity matrix.
+*
+      CALL ZLASET( 'Full', N, N, CZERO, CONE, C, LDC )
+*
+*     Call ZLAVHE_ROOK to form the product D * U' (or D * L' ).
+*
+      CALL ZLAVHE_ROOK( UPLO, 'Conjugate', 'Non-unit', N, N, AFAC,
+     $                  LDAFAC, IPIV, C, LDC, INFO )
+*
+*     Call ZLAVHE_ROOK again to multiply by U (or L ).
+*
+      CALL ZLAVHE_ROOK( UPLO, 'No transpose', 'Unit', N, N, AFAC,
+     $                  LDAFAC, IPIV, C, LDC, INFO )
+*
+*     Compute the difference  C - A .
+*
+      IF( LSAME( UPLO, 'U' ) ) THEN
+         DO 30 J = 1, N
+            DO 20 I = 1, J - 1
+               C( I, J ) = C( I, J ) - A( I, J )
+   20       CONTINUE
+            C( J, J ) = C( J, J ) - DBLE( A( J, J ) )
+   30    CONTINUE
+      ELSE
+         DO 50 J = 1, N
+            C( J, J ) = C( J, J ) - DBLE( A( J, J ) )
+            DO 40 I = J + 1, N
+               C( I, J ) = C( I, J ) - A( I, J )
+   40       CONTINUE
+   50    CONTINUE
+      END IF
+*
+*     Compute norm( C - A ) / ( N * norm(A) * EPS )
+*
+      RESID = ZLANHE( '1', UPLO, N, C, LDC, RWORK )
+*
+      IF( ANORM.LE.ZERO ) THEN
+         IF( RESID.NE.ZERO )
+     $      RESID = ONE / EPS
+      ELSE
+         RESID = ( ( RESID/DBLE( N ) )/ANORM ) / EPS
+      END IF
+*
+      RETURN
+*
+*     End of ZHET01_ROOK
+*
+      END
diff --git a/TESTING/LIN/zlavhe_rook.f b/TESTING/LIN/zlavhe_rook.f
new file mode 100644
index 00000000..8470591a
--- /dev/null
+++ b/TESTING/LIN/zlavhe_rook.f
@@ -0,0 +1,613 @@
+*> \brief \b ZLAVHE_ROOK
+*
+*  =========== DOCUMENTATION ===========
+*
+* Online html documentation available at 
+*            http://www.netlib.org/lapack/explore-html/ 
+*
+*  Definition:
+*  ===========
+*
+*       SUBROUTINE ZLAVHE_ROOK( UPLO, TRANS, DIAG, N, NRHS, A, LDA, IPIV, B,
+*                               LDB, INFO )
+* 
+*       .. Scalar Arguments ..
+*       CHARACTER          DIAG, TRANS, UPLO
+*       INTEGER            INFO, LDA, LDB, N, NRHS
+*       ..
+*       .. Array Arguments ..
+*       INTEGER            IPIV( * )
+*       COMPLEX*16         A( LDA, * ), B( LDB, * )
+*       ..
+*  
+*
+*> \par Purpose:
+*  =============
+*>
+*> \verbatim
+*>
+*> ZLAVHE_ROOK  performs one of the matrix-vector operations
+*>    x := A*x  or  x := A^H*x,
+*> where x is an N element vector and  A is one of the factors
+*> from the Hermitian factorization computed by CHETRF_ROOK.
+*>
+*> ZHETRF_ROOK produces a factorization of the form
+*>      U * D * U^H     or     L * D * L^H,
+*> where U (or L) is a product of permutation and unit upper (lower)
+*> triangular matrices, U^H (or L^H) is the conjugate transpose of
+*> U (or L), and D is Hermitian and block diagonal with 1 x 1 and
+*> 2 x 2 diagonal blocks.  The multipliers for the transformations
+*> and the upper or lower triangular parts of the diagonal blocks
+*> are stored in the leading upper or lower triangle of the 2-D
+*> array A.
+*>
+*> If TRANS = 'N', multiplies by U  or U * D  (or L  or L * D)
+*> If TRANS = 'T', multiplies by U' or D * U' (or L' or D * L')
+*> If TRANS = 'C', multiplies by U' or D * U' (or L' or D * L')
+*> \endverbatim
+*
+*  Arguments:
+*  ==========
+*
+*> \param[in] UPLO
+*> \verbatim
+*>          UPLO is CHARACTER*1
+*>          Specifies whether the factor stored in A is upper or lower
+*>          triangular.
+*>          = 'U':  Upper triangular
+*>          = 'L':  Lower triangular
+*> \endverbatim
+*>
+*> \param[in] TRANS
+*> \verbatim
+*>          TRANS is CHARACTER*1
+*>          Specifies the operation to be performed:
+*>          = 'N':  x := A*x
+*>          = 'T':   x := A^H*x
+*> \endverbatim
+*>
+*> \param[in] DIAG
+*> \verbatim
+*>          DIAG is CHARACTER*1
+*>          Specifies whether or not the diagonal blocks are unit
+*>          matrices.  If the diagonal blocks are assumed to be unit,
+*>          then A = U or A = L, otherwise A = U*D or A = L*D.
+*>          = 'U':  Diagonal blocks are assumed to be unit matrices.
+*>          = 'N':  Diagonal blocks are assumed to be non-unit matrices.
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*>          N is INTEGER
+*>          The number of rows and columns 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 vectors
+*>          x to be multiplied by A.  NRHS >= 0.
+*> \endverbatim
+*>
+*> \param[in] A
+*> \verbatim
+*>          A is COMPLEX*16 array, dimension (LDA,N)
+*>          The block diagonal matrix D and the multipliers used to
+*>          obtain the factor U or L as computed by ZHETRF_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 ZHETRF_ROOK.
+*>          If UPLO = 'U':
+*>             Only the last KB elements of IPIV are set.
+*>
+*>             If IPIV(k) > 0, then rows and columns k and IPIV(k) were
+*>             interchanged and D(k,k) is a 1-by-1 diagonal block.
+*>
+*>             If IPIV(k) < 0 and IPIV(k-1) < 0, then rows and
+*>             columns k and -IPIV(k) were interchanged and rows and
+*>             columns k-1 and -IPIV(k-1) were inerchaged,
+*>             D(k-1:k,k-1:k) is a 2-by-2 diagonal block.
+*>
+*>          If UPLO = 'L':
+*>             Only the first KB elements of IPIV are set.
+*>
+*>             If IPIV(k) > 0, then rows and columns k and IPIV(k)
+*>             were interchanged and D(k,k) is a 1-by-1 diagonal block.
+*>
+*>             If IPIV(k) < 0 and IPIV(k+1) < 0, then rows and
+*>             columns k and -IPIV(k) were interchanged and rows and
+*>             columns k+1 and -IPIV(k+1) were inerchaged,
+*>             D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
+*> \endverbatim
+*>
+*> \param[in,out] B
+*> \verbatim
+*>          B is COMPLEX array, dimension (LDB,NRHS)
+*>          On entry, B contains NRHS vectors of length N.
+*>          On exit, B is overwritten with the product A * B.
+*> \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 = -k, the k-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 April 2013
+*
+*> \ingroup complex16_lin
+*
+*  =====================================================================
+      SUBROUTINE ZLAVHE_ROOK( UPLO, TRANS, DIAG, N, NRHS, A, LDA, IPIV,
+     $                        B, LDB, INFO )
+*
+*  -- LAPACK test 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          DIAG, TRANS, UPLO
+      INTEGER            INFO, LDA, LDB, N, NRHS
+*     ..
+*     .. Array Arguments ..
+      INTEGER            IPIV( * )
+      COMPLEX*16         A( LDA, * ), B( LDB, * )
+*     ..
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      COMPLEX*16         CONE
+      PARAMETER          ( CONE = ( 1.0E+0, 0.0E+0 ) )
+*     ..
+*     .. Local Scalars ..
+      LOGICAL            NOUNIT
+      INTEGER            J, K, KP
+      COMPLEX*16         D11, D12, D21, D22, T1, T2
+*     ..
+*     .. External Functions ..
+      LOGICAL            LSAME
+      EXTERNAL           LSAME
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           ZGEMV, ZGERU, ZLACGV, ZSCAL, ZSWAP, XERBLA
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          ABS, DCONJG, MAX
+*     ..
+*     .. Executable Statements ..
+*
+*     Test the input parameters.
+*
+      INFO = 0
+      IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
+         INFO = -1
+      ELSE IF( .NOT.LSAME( TRANS, 'N' ) .AND. .NOT.LSAME( TRANS, 'C' ) )
+     $          THEN
+         INFO = -2
+      ELSE IF( .NOT.LSAME( DIAG, 'U' ) .AND. .NOT.LSAME( DIAG, 'N' ) )
+     $          THEN
+         INFO = -3
+      ELSE IF( N.LT.0 ) THEN
+         INFO = -4
+      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
+         INFO = -6
+      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
+         INFO = -9
+      END IF
+      IF( INFO.NE.0 ) THEN
+         CALL XERBLA( 'ZLAVHE_ROOK ', -INFO )
+         RETURN
+      END IF
+*
+*     Quick return if possible.
+*
+      IF( N.EQ.0 )
+     $   RETURN
+*
+      NOUNIT = LSAME( DIAG, 'N' )
+*------------------------------------------
+*
+*     Compute  B := A * B  (No transpose)
+*
+*------------------------------------------
+      IF( LSAME( TRANS, 'N' ) ) THEN
+*
+*        Compute  B := U*B
+*        where U = P(m)*inv(U(m))* ... *P(1)*inv(U(1))
+*
+         IF( LSAME( UPLO, 'U' ) ) THEN
+*
+*        Loop forward applying the transformations.
+*
+            K = 1
+   10       CONTINUE
+            IF( K.GT.N )
+     $         GO TO 30
+            IF( IPIV( K ).GT.0 ) THEN
+*
+*              1 x 1 pivot block
+*
+*              Multiply by the diagonal element if forming U * D.
+*
+               IF( NOUNIT )
+     $            CALL ZSCAL( NRHS, A( K, K ), B( K, 1 ), LDB )
+*
+*              Multiply by  P(K) * inv(U(K))  if K > 1.
+*
+               IF( K.GT.1 ) THEN
+*
+*                 Apply the transformation.
+*
+                  CALL ZGERU( K-1, NRHS, CONE, A( 1, K ), 1, B( K, 1 ),
+     $                        LDB, B( 1, 1 ), LDB )
+*
+*                 Interchange if P(K) != I.
+*
+                  KP = IPIV( K )
+                  IF( KP.NE.K )
+     $               CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+               END IF
+               K = K + 1
+            ELSE
+*
+*              2 x 2 pivot block
+*
+*              Multiply by the diagonal block if forming U * D.
+*
+               IF( NOUNIT ) THEN
+                  D11 = A( K, K )
+                  D22 = A( K+1, K+1 )
+                  D12 = A( K, K+1 )
+                  D21 = DCONJG( D12 )
+                  DO 20 J = 1, NRHS
+                     T1 = B( K, J )
+                     T2 = B( K+1, J )
+                     B( K, J ) = D11*T1 + D12*T2
+                     B( K+1, J ) = D21*T1 + D22*T2
+   20             CONTINUE
+               END IF
+*
+*              Multiply by  P(K) * inv(U(K))  if K > 1.
+*
+               IF( K.GT.1 ) THEN
+*
+*                 Apply the transformations.
+*
+                  CALL ZGERU( K-1, NRHS, CONE, A( 1, K ), 1, B( K, 1 ),
+     $                        LDB, B( 1, 1 ), LDB )
+                  CALL ZGERU( K-1, NRHS, CONE, A( 1, K+1 ), 1,
+     $                        B( K+1, 1 ), LDB, B( 1, 1 ), LDB )
+*
+*                 Interchange if a permutation was applied at the
+*                 K-th step of the factorization.
+*
+*                 Swap the first of pair with IMAXth
+*
+                  KP = ABS( IPIV( K ) )
+                  IF( KP.NE.K )
+     $               CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+*
+*                 NOW swap the first of pair with Pth
+*
+                  KP = ABS( IPIV( K+1 ) )
+                  IF( KP.NE.K+1 )
+     $               CALL ZSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ),
+     $                           LDB )
+               END IF
+               K = K + 2
+            END IF
+            GO TO 10
+   30       CONTINUE
+*
+*        Compute  B := L*B
+*        where L = P(1)*inv(L(1))* ... *P(m)*inv(L(m)) .
+*
+         ELSE
+*
+*           Loop backward applying the transformations to B.
+*
+            K = N
+   40       CONTINUE
+            IF( K.LT.1 )
+     $         GO TO 60
+*
+*           Test the pivot index.  If greater than zero, a 1 x 1
+*           pivot was used, otherwise a 2 x 2 pivot was used.
+*
+            IF( IPIV( K ).GT.0 ) THEN
+*
+*              1 x 1 pivot block:
+*
+*              Multiply by the diagonal element if forming L * D.
+*
+               IF( NOUNIT )
+     $            CALL ZSCAL( NRHS, A( K, K ), B( K, 1 ), LDB )
+*
+*              Multiply by  P(K) * inv(L(K))  if K < N.
+*
+               IF( K.NE.N ) THEN
+                  KP = IPIV( K )
+*
+*                 Apply the transformation.
+*
+                  CALL ZGERU( N-K, NRHS, CONE, A( K+1, K ), 1,
+     $                        B( K, 1 ), LDB, B( K+1, 1 ), LDB )
+*
+*                 Interchange if a permutation was applied at the
+*                 K-th step of the factorization.
+*
+                  IF( KP.NE.K )
+     $               CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+               END IF
+               K = K - 1
+*
+            ELSE
+*
+*              2 x 2 pivot block:
+*
+*              Multiply by the diagonal block if forming L * D.
+*
+               IF( NOUNIT ) THEN
+                  D11 = A( K-1, K-1 )
+                  D22 = A( K, K )
+                  D21 = A( K, K-1 )
+                  D12 = DCONJG( D21 )
+                  DO 50 J = 1, NRHS
+                     T1 = B( K-1, J )
+                     T2 = B( K, J )
+                     B( K-1, J ) = D11*T1 + D12*T2
+                     B( K, J ) = D21*T1 + D22*T2
+   50             CONTINUE
+               END IF
+*
+*              Multiply by  P(K) * inv(L(K))  if K < N.
+*
+               IF( K.NE.N ) THEN
+*
+*                 Apply the transformation.
+*
+                  CALL ZGERU( N-K, NRHS, CONE, A( K+1, K ), 1,
+     $                        B( K, 1 ), LDB, B( K+1, 1 ), LDB )
+                  CALL ZGERU( N-K, NRHS, CONE, A( K+1, K-1 ), 1,
+     $                        B( K-1, 1 ), LDB, B( K+1, 1 ), LDB )
+*
+*                 Interchange if a permutation was applied at the
+*                 K-th step of the factorization.
+*
+*
+*                 Swap the second of pair with IMAXth
+*
+                  KP = ABS( IPIV( K ) )
+                  IF( KP.NE.K )
+     $               CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+*
+*                 NOW swap the first of pair with Pth
+*
+                  KP = ABS( IPIV( K-1 ) )
+                  IF( KP.NE.K-1 )
+     $               CALL ZSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ),
+     $                           LDB )
+*
+               END IF
+               K = K - 2
+            END IF
+            GO TO 40
+   60       CONTINUE
+         END IF
+*--------------------------------------------------
+*
+*     Compute  B := A^H * B  (conjugate transpose)
+*
+*--------------------------------------------------
+      ELSE
+*
+*        Form  B := U^H*B
+*        where U  = P(m)*inv(U(m))* ... *P(1)*inv(U(1))
+*        and   U^H = inv(U^H(1))*P(1)* ... *inv(U^H(m))*P(m)
+*
+         IF( LSAME( UPLO, 'U' ) ) THEN
+*
+*           Loop backward applying the transformations.
+*
+            K = N
+   70       IF( K.LT.1 )
+     $         GO TO 90
+*
+*           1 x 1 pivot block.
+*
+            IF( IPIV( K ).GT.0 ) THEN
+               IF( K.GT.1 ) THEN
+*
+*                 Interchange if P(K) != I.
+*
+                  KP = IPIV( K )
+                  IF( KP.NE.K )
+     $               CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+*
+*                 Apply the transformation
+*                    y = y - B' DCONJG(x),
+*                 where x is a column of A and y is a row of B.
+*
+                  CALL ZLACGV( NRHS, B( K, 1 ), LDB )
+                  CALL ZGEMV( 'Conjugate', K-1, NRHS, CONE, B, LDB,
+     $                        A( 1, K ), 1, CONE, B( K, 1 ), LDB )
+                  CALL ZLACGV( NRHS, B( K, 1 ), LDB )
+               END IF
+               IF( NOUNIT )
+     $            CALL ZSCAL( NRHS, A( K, K ), B( K, 1 ), LDB )
+               K = K - 1
+*
+*           2 x 2 pivot block.
+*
+            ELSE
+               IF( K.GT.2 ) THEN
+*
+*                 Swap the second of pair with Pth
+*
+                  KP = ABS( IPIV( K ) )
+                  IF( KP.NE.K )
+     $               CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+*
+*                 Now swap the first of pair with IMAX(r)th
+*
+                  KP = ABS( IPIV( K-1 ) )
+                  IF( KP.NE.K-1 )
+     $               CALL ZSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ),
+     $                           LDB )
+*
+*                 Apply the transformations
+*                    y = y - B' DCONJG(x),
+*                 where x is a block column of A and y is a block
+*                 row of B.
+*
+                  CALL ZLACGV( NRHS, B( K, 1 ), LDB )
+                  CALL ZGEMV( 'Conjugate', K-2, NRHS, CONE, B, LDB,
+     $                        A( 1, K ), 1, CONE, B( K, 1 ), LDB )
+                  CALL ZLACGV( NRHS, B( K, 1 ), LDB )
+*
+                  CALL ZLACGV( NRHS, B( K-1, 1 ), LDB )
+                  CALL ZGEMV( 'Conjugate', K-2, NRHS, CONE, B, LDB,
+     $                        A( 1, K-1 ), 1, CONE, B( K-1, 1 ), LDB )
+                  CALL ZLACGV( NRHS, B( K-1, 1 ), LDB )
+               END IF
+*
+*              Multiply by the diagonal block if non-unit.
+*
+               IF( NOUNIT ) THEN
+                  D11 = A( K-1, K-1 )
+                  D22 = A( K, K )
+                  D12 = A( K-1, K )
+                  D21 = DCONJG( D12 )
+                  DO 80 J = 1, NRHS
+                     T1 = B( K-1, J )
+                     T2 = B( K, J )
+                     B( K-1, J ) = D11*T1 + D12*T2
+                     B( K, J ) = D21*T1 + D22*T2
+   80             CONTINUE
+               END IF
+               K = K - 2
+            END IF
+            GO TO 70
+   90       CONTINUE
+*
+*        Form  B := L^H*B
+*        where L  = P(1)*inv(L(1))* ... *P(m)*inv(L(m))
+*        and   L^H = inv(L^H(m))*P(m)* ... *inv(L^H(1))*P(1)
+*
+         ELSE
+*
+*           Loop forward applying the L-transformations.
+*
+            K = 1
+  100       CONTINUE
+            IF( K.GT.N )
+     $         GO TO 120
+*
+*           1 x 1 pivot block
+*
+            IF( IPIV( K ).GT.0 ) THEN
+               IF( K.LT.N ) THEN
+*
+*                 Interchange if P(K) != I.
+*
+                  KP = IPIV( K )
+                  IF( KP.NE.K )
+     $               CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+*
+*                 Apply the transformation
+*
+                  CALL ZLACGV( NRHS, B( K, 1 ), LDB )
+                  CALL ZGEMV( 'Conjugate', N-K, NRHS, CONE, B( K+1, 1 ),
+     $                       LDB, A( K+1, K ), 1, CONE, B( K, 1 ), LDB )
+                  CALL ZLACGV( NRHS, B( K, 1 ), LDB )
+               END IF
+               IF( NOUNIT )
+     $            CALL ZSCAL( NRHS, A( K, K ), B( K, 1 ), LDB )
+               K = K + 1
+*
+*           2 x 2 pivot block.
+*
+            ELSE
+               IF( K.LT.N-1 ) THEN
+*
+*                 Swap the first of pair with Pth
+*
+                  KP = ABS( IPIV( K ) )
+                  IF( KP.NE.K )
+     $               CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+*
+*                 Now swap the second of pair with IMAX(r)th
+*
+                  KP = ABS( IPIV( K+1 ) )
+                  IF( KP.NE.K+1 )
+     $               CALL ZSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ),
+     $                           LDB )
+*
+*                 Apply the transformation
+*
+                  CALL ZLACGV( NRHS, B( K+1, 1 ), LDB )
+                  CALL ZGEMV( 'Conjugate', N-K-1, NRHS, CONE,
+     $                        B( K+2, 1 ), LDB, A( K+2, K+1 ), 1, CONE,
+     $                        B( K+1, 1 ), LDB )
+                  CALL ZLACGV( NRHS, B( K+1, 1 ), LDB )
+*
+                  CALL ZLACGV( NRHS, B( K, 1 ), LDB )
+                  CALL ZGEMV( 'Conjugate', N-K-1, NRHS, CONE,
+     $                        B( K+2, 1 ), LDB, A( K+2, K ), 1, CONE,
+     $                        B( K, 1 ), LDB )
+                  CALL ZLACGV( NRHS, B( K, 1 ), LDB )
+               END IF
+*
+*              Multiply by the diagonal block if non-unit.
+*
+               IF( NOUNIT ) THEN
+                  D11 = A( K, K )
+                  D22 = A( K+1, K+1 )
+                  D21 = A( K+1, K )
+                  D12 = DCONJG( D21 )
+                  DO 110 J = 1, NRHS
+                     T1 = B( K, J )
+                     T2 = B( K+1, J )
+                     B( K, J ) = D11*T1 + D12*T2
+                     B( K+1, J ) = D21*T1 + D22*T2
+  110             CONTINUE
+               END IF
+               K = K + 2
+            END IF
+            GO TO 100
+  120       CONTINUE
+         END IF
+*
+      END IF
+      RETURN
+*
+*     End of ZLAVHE_ROOK
+*
+      END