Adding new hetri routines

18 files changed, 1640 insertions, 36 deletions

diff --git a/SRC/CMakeLists.txt b/SRC/CMakeLists.txt
index 4f9f8a33..4cbfff1a 100644
--- a/SRC/CMakeLists.txt
+++ b/SRC/CMakeLists.txt
@@ -172,7 +172,8 @@ set(CLASRC
    checon.f cheev.f  cheevd.f cheevr.f cheevx.f chegs2.f chegst.f 
    chegv.f  chegvd.f chegvx.f cherfs.f chesv.f  chesvx.f chetd2.f 
    chetf2.f chetrd.f 
-   chetrf.f chetri.f chetrs.f chetrs2.f chgeqz.f chpcon.f chpev.f  chpevd.f 
+   chetrf.f chetri.f chetri2.f  chetri2x.f cheswapr.f
+   chetrs.f chetrs2.f chgeqz.f chpcon.f chpev.f  chpevd.f 
    chpevx.f chpgst.f chpgv.f  chpgvd.f chpgvx.f chprfs.f chpsv.f  
    chpsvx.f 
    chptrd.f chptrf.f chptri.f chptrs.f chsein.f chseqr.f clabrd.f 
@@ -312,7 +313,8 @@ set(ZLASRC
    zhecon.f zheev.f  zheevd.f zheevr.f zheevx.f zhegs2.f zhegst.f 
    zhegv.f  zhegvd.f zhegvx.f zherfs.f zhesv.f  zhesvx.f zhetd2.f 
    zhetf2.f zhetrd.f 
-   zhetrf.f zhetri.f zhetrs.f zhetrs2.f zhgeqz.f zhpcon.f zhpev.f  zhpevd.f 
+   zhetrf.f zhetri.f zhetri2.f zhetri2x.f zheswapr.f
+   zhetrs.f zhetrs2.f zhgeqz.f zhpcon.f zhpev.f  zhpevd.f 
    zhpevx.f zhpgst.f zhpgv.f  zhpgvd.f zhpgvx.f zhprfs.f zhpsv.f  
    zhpsvx.f 
    zhptrd.f zhptrf.f zhptri.f zhptrs.f zhsein.f zhseqr.f zlabrd.f 
diff --git a/SRC/Makefile b/SRC/Makefile
index 9fe35696..59756832 100644
--- a/SRC/Makefile
+++ b/SRC/Makefile
@@ -173,7 +173,8 @@ CLASRC = \
    checon.o cheev.o  cheevd.o cheevr.o cheevx.o chegs2.o chegst.o \
    chegv.o  chegvd.o chegvx.o cherfs.o chesv.o  chesvx.o chetd2.o \
    chetf2.o chetrd.o \
-   chetrf.o chetri.o chetrs.o chetrs2.o chgeqz.o chpcon.o chpev.o  chpevd.o \
+   chetrf.o chetri.o chetri2.o chetri2x.o cheswapr.o \
+   chetrs.o chetrs2.o chgeqz.o chpcon.o chpev.o  chpevd.o \
    chpevx.o chpgst.o chpgv.o  chpgvd.o chpgvx.o chprfs.o chpsv.o  \
    chpsvx.o \
    chptrd.o chptrf.o chptri.o chptrs.o chsein.o chseqr.o clabrd.o \
@@ -312,7 +313,8 @@ ZLASRC = \
    zhecon.o zheev.o  zheevd.o zheevr.o zheevx.o zhegs2.o zhegst.o \
    zhegv.o  zhegvd.o zhegvx.o zherfs.o zhesv.o  zhesvx.o zhetd2.o \
    zhetf2.o zhetrd.o \
-   zhetrf.o zhetri.o zhetrs.o zhetrs2.o zhgeqz.o zhpcon.o zhpev.o  zhpevd.o \
+   zhetrf.o zhetri.o zhetri2.o zhetri2x.o zheswapr.o \
+   zhetrs.o zhetrs2.o zhgeqz.o zhpcon.o zhpev.o  zhpevd.o \
    zhpevx.o zhpgst.o zhpgv.o  zhpgvd.o zhpgvx.o zhprfs.o zhpsv.o  \
    zhpsvx.o \
    zhptrd.o zhptrf.o zhptri.o zhptrs.o zhsein.o zhseqr.o zlabrd.o \
diff --git a/SRC/cheswapr.f b/SRC/cheswapr.f
new file mode 100644
index 00000000..3bb2c226
--- /dev/null
+++ b/SRC/cheswapr.f
@@ -0,0 +1,136 @@
+      SUBROUTINE CHESWAPR( UPLO, N, A, I1, I2)
+*
+*  -- LAPACK auxiliary routine (version 3.3.0) --
+*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*     November 2010
+*
+*     .. Scalar Arguments ..
+      CHARACTER        UPLO
+      INTEGER          I1, I2, N
+*     ..
+*     .. Array Arguments ..
+      COMPLEX          A(N,N)
+*
+*  Purpose
+*  =======
+*
+*  CHESWAPR applies an elementary permutation on the rows and the columns of
+*  a hermitian matrix.
+*
+*  Arguments
+*  =========
+*
+*  UPLO    (input) CHARACTER*1
+*          Specifies whether the details of the factorization are stored
+*          as an upper or lower triangular matrix.
+*          = 'U':  Upper triangular, form is A = U*D*U**T;
+*          = 'L':  Lower triangular, form is A = L*D*L**T.
+*
+*  N       (input) INTEGER
+*          The order of the matrix A.  N >= 0.
+*
+*  A       (input/output) COMPLEX array, dimension (LDA,N)
+*          On entry, the NB diagonal matrix D and the multipliers
+*          used to obtain the factor U or L as computed by CSYTRF.
+*
+*          On exit, if INFO = 0, the (symmetric) inverse of the original
+*          matrix.  If UPLO = 'U', the upper triangular part of the
+*          inverse is formed and the part of A below the diagonal is not
+*          referenced; if UPLO = 'L' the lower triangular part of the
+*          inverse is formed and the part of A above the diagonal is
+*          not referenced.
+*
+*  I1      (input) INTEGER
+*          Index of the first row to swap
+*
+*  I2      (input) INTEGER
+*          Index of the second row to swap
+*
+*  =====================================================================
+*
+*     ..
+*     .. Local Scalars ..
+      LOGICAL            UPPER
+      INTEGER            I
+      COMPLEX            TMP
+*
+*     .. External Functions ..
+      LOGICAL            LSAME
+      EXTERNAL           LSAME
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           CSWAP
+*     ..
+*     .. Executable Statements ..
+*
+      UPPER = LSAME( UPLO, 'U' )
+      IF (UPPER) THEN
+*
+*         UPPER
+*         first swap
+*          - swap column I1 and I2 from I1 to I1-1 
+         CALL CSWAP( I1-1, A(1,I1), 1, A(1,I2), 1 )
+*
+*          second swap :
+*          - swap A(I1,I1) and A(I2,I2)
+*          - swap row I1 from I1+1 to I2-1 with col I2 from I1+1 to I2-1
+*          - swap A(I2,I1) and A(I1,I2)
+     
+         TMP=A(I1,I1)
+         A(I1,I1)=A(I2,I2)
+         A(I2,I2)=TMP
+*
+         DO I=1,I2-I1-1
+            TMP=A(I1,I1+I)
+            A(I1,I1+I)=CONJG(A(I1+I,I2))
+            A(I1+I,I2)=CONJG(TMP)
+         END DO
+*
+          A(I1,I2)=CONJG(A(I1,I2))
+
+*
+*          third swap
+*          - swap row I1 and I2 from I2+1 to N
+         DO I=I2+1,N
+            TMP=A(I1,I)
+            A(I1,I)=A(I2,I)
+            A(I2,I)=TMP
+         END DO
+*
+        ELSE
+*
+*         LOWER
+*         first swap
+*          - swap row I1 and I2 from 1 to I1-1 
+         CALL CSWAP ( I1-1, A(I1,1), N, A(I2,1), N )
+*
+*         second swap :
+*          - swap A(I1,I1) and A(I2,I2)
+*          - swap col I1 from I1+1 to I2-1 with row I2 from I1+1 to I2-1     
+*          - swap A(I2,I1) and A(I1,I2)
+
+          TMP=A(I1,I1)
+          A(I1,I1)=A(I2,I2)
+          A(I2,I2)=TMP
+*
+          DO I=1,I2-I1-1
+             TMP=A(I1+I,I1)
+             A(I1+I,I1)=CONJG(A(I2,I1+I))
+             A(I2,I1+I)=CONJG(TMP)
+          END DO
+*
+          A(I2,I1)=CONJG(A(I2,I1))
+*
+*         third swap
+*          - swap col I1 and I2 from I2+1 to N
+          DO I=I2+1,N
+             TMP=A(I,I1)
+             A(I,I1)=A(I,I2)
+             A(I,I2)=TMP
+          END DO
+*
+      ENDIF
+      
+      END SUBROUTINE CHESWAPR
+
diff --git a/SRC/chetri2.f b/SRC/chetri2.f
new file mode 100644
index 00000000..f8c55ee9
--- /dev/null
+++ b/SRC/chetri2.f
@@ -0,0 +1,127 @@
+      SUBROUTINE CHETRI2( UPLO, N, A, LDA, IPIV, WORK, LWORK, INFO )
+*
+*  -- LAPACK routine (version 3.3.0) --
+*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*     November 2010
+*
+*  -- Written by Julie Langou of the Univ. of TN    --
+*
+*     .. Scalar Arguments ..
+      CHARACTER          UPLO
+      INTEGER            INFO, LDA, LWORK, N
+*     ..
+*     .. Array Arguments ..
+      INTEGER            IPIV( * )
+      COMPLEX            A( LDA, * ), WORK( * )
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  CHETRI2 computes the inverse of a complex hermitian indefinite matrix
+*  A using the factorization A = U*D*U**T or A = L*D*L**T computed by
+*  CHETRF. CHETRI2 set the LEADING DIMENSION of the workspace
+*  before calling CHETRI2X that actually compute the inverse.
+*
+*  Arguments
+*  =========
+*
+*  UPLO    (input) CHARACTER*1
+*          Specifies whether the details of the factorization are stored
+*          as an upper or lower triangular matrix.
+*          = 'U':  Upper triangular, form is A = U*D*U**T;
+*          = 'L':  Lower triangular, form is A = L*D*L**T.
+*
+*  N       (input) INTEGER
+*          The order of the matrix A.  N >= 0.
+*
+*  A       (input/output) COMPLEX array, dimension (LDA,N)
+*          On entry, the NB diagonal matrix D and the multipliers
+*          used to obtain the factor U or L as computed by CHETRF.
+*
+*          On exit, if INFO = 0, the (symmetric) inverse of the original
+*          matrix.  If UPLO = 'U', the upper triangular part of the
+*          inverse is formed and the part of A below the diagonal is not
+*          referenced; if UPLO = 'L' the lower triangular part of the
+*          inverse is formed and the part of A above the diagonal is
+*          not referenced.
+*
+*  LDA     (input) INTEGER
+*          The leading dimension of the array A.  LDA >= max(1,N).
+*
+*  IPIV    (input) INTEGER array, dimension (N)
+*          Details of the interchanges and the NB structure of D
+*          as determined by CHETRF.
+*
+*  WORK    (workspace) COMPLEX array, dimension (N+NB+1)*(NB+3)
+*
+*  LWORK   (input) INTEGER
+*          The dimension of the array WORK.
+*          WORK is size >= (N+NB+1)*(NB+3)
+*          If LDWORK = -1, then a workspace query is assumed; the routine
+*           calculates:
+*              - the optimal size of the WORK array, returns
+*          this value as the first entry of the WORK array,
+*              - and no error message related to LDWORK is issued by XERBLA.
+*
+*  INFO    (output) INTEGER
+*          = 0: successful exit
+*          < 0: if INFO = -i, the i-th argument had an illegal value
+*          > 0: if INFO = i, D(i,i) = 0; the matrix is singular and its
+*               inverse could not be computed.
+*
+*  =====================================================================
+*
+*     .. Local Scalars ..
+      LOGICAL            UPPER, LQUERY
+      INTEGER            MINSIZE, NBMAX
+*     ..
+*     .. External Functions ..
+      LOGICAL            LSAME
+      INTEGER            ILAENV
+      EXTERNAL           LSAME, ILAENV
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           CHETRI2X
+*     ..
+*     .. Executable Statements ..
+*
+*     Test the input parameters.
+*
+      INFO = 0
+      UPPER = LSAME( UPLO, 'U' )
+      LQUERY = ( LWORK.EQ.-1 )
+*     Get blocksize
+      NBMAX = ILAENV( 1, 'CHETRF', UPLO, N, -1, -1, -1 )
+      MINSIZE = (N+NBMAX+1)*(NBMAX+3)
+*
+      IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
+         INFO = -1
+      ELSE IF( N.LT.0 ) THEN
+         INFO = -2
+      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
+         INFO = -4
+      ELSE IF (LWORK .LT. MINSIZE .AND. .NOT.LQUERY ) THEN
+         INFO = -7
+      END IF
+*
+*     Quick return if possible
+*
+*
+      IF( INFO.NE.0 ) THEN
+         CALL XERBLA( 'CHETRI2', -INFO )
+         RETURN
+      ELSE IF( LQUERY ) THEN
+         WORK(1)=(N+NBMAX+1)*(NBMAX+3)
+         RETURN
+      END IF
+      IF( N.EQ.0 )
+     $   RETURN
+      
+      CALL CHETRI2X( UPLO, N, A, LDA, IPIV, WORK, NBMAX, INFO )
+      RETURN
+*
+*     End of CHETRI2
+*
+      END
diff --git a/SRC/chetri2x.f b/SRC/chetri2x.f
new file mode 100644
index 00000000..91699060
--- /dev/null
+++ b/SRC/chetri2x.f
@@ -0,0 +1,508 @@
+      SUBROUTINE CHETRI2X( UPLO, N, A, LDA, IPIV, WORK, NB, INFO )
+*
+*  -- LAPACK routine (version 3.3.0) --
+*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*     November 2010
+*
+*  -- Written by Julie Langou of the Univ. of TN    --
+*
+*     .. Scalar Arguments ..
+      CHARACTER          UPLO
+      INTEGER            INFO, LDA, N, NB
+*     ..
+*     .. Array Arguments ..
+      INTEGER            IPIV( * )
+      COMPLEX            A( LDA, * ), WORK( N+NB+1,* )
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  CHETRI2X computes the inverse of a complex Hermitian indefinite matrix
+*  A using the factorization A = U*D*U**H or A = L*D*L**H computed by
+*  CHETRF.
+*
+*  Arguments
+*  =========
+*
+*  UPLO    (input) CHARACTER*1
+*          Specifies whether the details of the factorization are stored
+*          as an upper or lower triangular matrix.
+*          = 'U':  Upper triangular, form is A = U*D*U**H;
+*          = 'L':  Lower triangular, form is A = L*D*L**H.
+*
+*  N       (input) INTEGER
+*          The order of the matrix A.  N >= 0.
+*
+*  A       (input/output) COMPLEX array, dimension (LDA,N)
+*          On entry, the NNB diagonal matrix D and the multipliers
+*          used to obtain the factor U or L as computed by CHETRF.
+*
+*          On exit, if INFO = 0, the (symmetric) inverse of the original
+*          matrix.  If UPLO = 'U', the upper triangular part of the
+*          inverse is formed and the part of A below the diagonal is not
+*          referenced; if UPLO = 'L' the lower triangular part of the
+*          inverse is formed and the part of A above the diagonal is
+*          not referenced.
+*
+*  LDA     (input) INTEGER
+*          The leading dimension of the array A.  LDA >= max(1,N).
+*
+*  IPIV    (input) INTEGER array, dimension (N)
+*          Details of the interchanges and the NNB structure of D
+*          as determined by CHETRF.
+*
+*  WORK    (workspace) COMPLEX array, dimension (N+NNB+1,NNB+3)
+*
+*  NB      (input) INTEGER
+*          Block size
+*
+*  INFO    (output) INTEGER
+*          = 0: successful exit
+*          < 0: if INFO = -i, the i-th argument had an illegal value
+*          > 0: if INFO = i, D(i,i) = 0; the matrix is singular and its
+*               inverse could not be computed.
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      REAL               ONE
+      COMPLEX            CONE, ZERO
+      PARAMETER          ( ONE = 1.0E+0,
+     $                   CONE = ( 1.0E+0, 0.0E+0 ),
+     $                   ZERO = ( 0.0E+0, 0.0E+0 ) )
+*     ..
+*     .. Local Scalars ..
+      LOGICAL            UPPER
+      INTEGER            I, IINFO, IP, K, CUT, NNB
+      INTEGER            COUNT
+      INTEGER            J, U11, INVD
+
+      COMPLEX   AK, AKKP1, AKP1, D, T
+      COMPLEX   U01_I_J, U01_IP1_J
+      COMPLEX   U11_I_J, U11_IP1_J
+*     ..
+*     .. External Functions ..
+      LOGICAL            LSAME
+      EXTERNAL           LSAME
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           CSYCONV, XERBLA, CTRTRI
+      EXTERNAL           CGEMM, CTRMM, CHESWAPR
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          MAX
+*     ..
+*     .. Executable Statements ..
+*
+*     Test the input parameters.
+*
+      INFO = 0
+      UPPER = LSAME( UPLO, 'U' )
+      IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
+         INFO = -1
+      ELSE IF( N.LT.0 ) THEN
+         INFO = -2
+      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
+         INFO = -4
+      END IF
+*
+*     Quick return if possible
+*
+*
+      IF( INFO.NE.0 ) THEN
+         CALL XERBLA( 'CHETRI2X', -INFO )
+         RETURN
+      END IF
+      IF( N.EQ.0 )
+     $   RETURN
+*
+*     Convert A
+*     Workspace got Non-diag elements of D
+*
+      CALL CSYCONV( UPLO, 'C', N, A, LDA, IPIV, WORK, IINFO )
+*
+*     Check that the diagonal matrix D is nonsingular.
+*
+      IF( UPPER ) THEN
+*
+*        Upper triangular storage: examine D from bottom to top
+*
+         DO INFO = N, 1, -1
+            IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO )
+     $         RETURN
+         END DO
+      ELSE
+*
+*        Lower triangular storage: examine D from top to bottom.
+*
+         DO INFO = 1, N
+            IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO )
+     $         RETURN
+         END DO
+      END IF
+      INFO = 0
+*
+*  Splitting Workspace
+*     U01 is a block (N,NB+1) 
+*     The first element of U01 is in WORK(1,1)
+*     U11 is a block (NB+1,NB+1)
+*     The first element of U11 is in WORK(N+1,1)
+      U11 = N
+*     INVD is a block (N,2)
+*     The first element of INVD is in WORK(1,INVD)
+      INVD = NB+2
+
+      IF( UPPER ) THEN
+*
+*        invA = P * inv(U')*inv(D)*inv(U)*P'.
+*
+        CALL CTRTRI( UPLO, 'U', N, A, LDA, INFO )
+*
+*       inv(D) and inv(D)*inv(U)
+* 
+        K=1
+        DO WHILE ( K .LE. N )
+         IF( IPIV( K ).GT.0 ) THEN
+*           1 x 1 diagonal NNB
+             WORK(K,INVD) = ONE / REAL ( A( K, K ) )
+             WORK(K,INVD+1) = 0
+            K=K+1
+         ELSE
+*           2 x 2 diagonal NNB
+             T = ABS ( WORK(K+1,1) )
+             AK = REAL ( A( K, K ) ) / T
+             AKP1 = REAL ( A( K+1, K+1 ) ) / T
+             AKKP1 = WORK(K+1,1)  / T
+             D = T*( AK*AKP1-ONE )
+             WORK(K,INVD) = AKP1 / D
+             WORK(K+1,INVD+1) = AK / D
+             WORK(K,INVD+1) = -AKKP1 / D  
+             WORK(K+1,INVD) = CONJG (WORK(K,INVD+1) )
+            K=K+2
+         END IF
+        END DO
+*
+*       inv(U') = (inv(U))'
+*
+*       inv(U')*inv(D)*inv(U)
+*
+        CUT=N
+        DO WHILE (CUT .GT. 0)
+           NNB=NB
+           IF (CUT .LE. NNB) THEN
+              NNB=CUT
+           ELSE
+              COUNT = 0
+*             count negative elements, 
+              DO I=CUT+1-NNB,CUT
+                  IF (IPIV(I) .LT. 0) COUNT=COUNT+1
+              END DO
+*             need a even number for a clear cut
+              IF (MOD(COUNT,2) .EQ. 1) NNB=NNB+1
+           END IF
+
+           CUT=CUT-NNB
+*
+*          U01 Block 
+*
+           DO I=1,CUT
+             DO J=1,NNB
+              WORK(I,J)=A(I,CUT+J)
+             END DO
+           END DO
+*
+*          U11 Block
+*
+           DO I=1,NNB
+             WORK(U11+I,I)=CONE
+             DO J=1,I-1
+                WORK(U11+I,J)=ZERO
+             END DO
+             DO J=I+1,NNB
+                WORK(U11+I,J)=A(CUT+I,CUT+J)
+             END DO
+           END DO
+*
+*          invD*U01
+*
+           I=1
+           DO WHILE (I .LE. CUT)
+             IF (IPIV(I) > 0) THEN
+                DO J=1,NNB
+                    WORK(I,J)=WORK(I,INVD)*WORK(I,J)
+                END DO
+                I=I+1
+             ELSE
+                DO J=1,NNB
+                   U01_I_J = WORK(I,J)
+                   U01_IP1_J = WORK(I+1,J)
+                   WORK(I,J)=WORK(I,INVD)*U01_I_J+
+     $                      WORK(I,INVD+1)*U01_IP1_J
+                   WORK(I+1,J)=WORK(I+1,INVD)*U01_I_J+
+     $                      WORK(I+1,INVD+1)*U01_IP1_J
+                END DO
+                I=I+2
+             END IF
+           END DO
+*
+*        invD1*U11
+*
+           I=1
+           DO WHILE (I .LE. NNB)
+             IF (IPIV(CUT+I) > 0) THEN
+                DO J=I,NNB
+                    WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J)
+                END DO
+                I=I+1
+             ELSE
+                DO J=I,NNB
+                   U11_I_J = WORK(U11+I,J)
+                   U11_IP1_J = WORK(U11+I+1,J)
+                WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J) +
+     $                      WORK(CUT+I,INVD+1)*WORK(U11+I+1,J)
+                WORK(U11+I+1,J)=WORK(CUT+I+1,INVD)*U11_I_J+
+     $                      WORK(CUT+I+1,INVD+1)*U11_IP1_J
+                END DO
+                I=I+2
+             END IF
+           END DO
+*    
+*       U11T*invD1*U11->U11
+*
+        CALL CTRMM('L','U','C','U',NNB, NNB,
+     $             CONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1)
+*
+*          U01'invD*U01->A(CUT+I,CUT+J)
+*
+         CALL CGEMM('C','N',NNB,NNB,CUT,CONE,A(1,CUT+1),LDA,
+     $              WORK,N+NB+1, ZERO, A(CUT+1,CUT+1), LDA)
+*
+*        U11 =  U11T*invD1*U11 + U01'invD*U01
+*
+         DO I=1,NNB
+            DO J=I,NNB
+              A(CUT+I,CUT+J)=A(CUT+I,CUT+J)+WORK(U11+I,J)
+            END DO
+         END DO
+*
+*        U01 =  U00T*invD0*U01
+*
+         CALL CTRMM('L',UPLO,'C','U',CUT, NNB,
+     $             CONE,A,LDA,WORK,N+NB+1)
+
+*
+*        Update U01
+*
+         DO I=1,CUT
+           DO J=1,NNB
+            A(I,CUT+J)=WORK(I,J)
+           END DO
+         END DO
+*
+*      Next Block
+*
+       END DO
+*
+*        Apply PERMUTATIONS P and P': P * inv(U')*inv(D)*inv(U) *P'
+*  
+            I=1
+            DO WHILE ( I .LE. N )
+               IF( IPIV(I) .GT. 0 ) THEN
+                  IP=IPIV(I)
+                 IF (I .LT. IP) CALL CHESWAPR( UPLO, N, A, I ,IP )
+                 IF (I .GT. IP) CALL CHESWAPR( UPLO, N, A, IP ,I )
+               ELSE
+                 IP=-IPIV(I)
+                 I=I+1
+                 IF ( (I-1) .LT. IP) 
+     $                  CALL CHESWAPR( UPLO, N, A, I-1 ,IP )
+                 IF ( (I-1) .GT. IP) 
+     $                  CALL CHESWAPR( UPLO, N, A, IP ,I-1 )
+              ENDIF
+               I=I+1
+            END DO
+      ELSE
+*
+*        LOWER...
+*
+*        invA = P * inv(U')*inv(D)*inv(U)*P'.
+*
+         CALL CTRTRI( UPLO, 'U', N, A, LDA, INFO )
+*
+*       inv(D) and inv(D)*inv(U)
+* 
+        K=N
+        DO WHILE ( K .GE. 1 )
+         IF( IPIV( K ).GT.0 ) THEN
+*           1 x 1 diagonal NNB
+             WORK(K,INVD) = ONE / REAL ( A( K, K ) )
+             WORK(K,INVD+1) = 0
+            K=K-1
+         ELSE
+*           2 x 2 diagonal NNB
+             T = ABS ( WORK(K-1,1) )
+             AK = REAL ( A( K-1, K-1 ) ) / T
+             AKP1 = REAL ( A( K, K ) ) / T
+             AKKP1 = WORK(K-1,1) / T
+             D = T*( AK*AKP1-ONE )
+             WORK(K-1,INVD) = AKP1 / D
+             WORK(K,INVD) = AK / D
+             WORK(K,INVD+1) = -AKKP1 / D  
+             WORK(K-1,INVD+1) = CONJG (WORK(K,INVD+1) )
+            K=K-2
+         END IF
+        END DO
+*
+*       inv(U') = (inv(U))'
+*
+*       inv(U')*inv(D)*inv(U)
+*
+        CUT=0
+        DO WHILE (CUT .LT. N)
+           NNB=NB
+           IF (CUT + NNB .GE. N) THEN
+              NNB=N-CUT
+           ELSE
+              COUNT = 0
+*             count negative elements, 
+              DO I=CUT+1,CUT+NNB
+                  IF (IPIV(I) .LT. 0) COUNT=COUNT+1
+              END DO
+*             need a even number for a clear cut
+              IF (MOD(COUNT,2) .EQ. 1) NNB=NNB+1
+           END IF
+*      L21 Block
+           DO I=1,N-CUT-NNB
+             DO J=1,NNB
+              WORK(I,J)=A(CUT+NNB+I,CUT+J)
+             END DO
+           END DO
+*     L11 Block
+           DO I=1,NNB
+             WORK(U11+I,I)=CONE
+             DO J=I+1,NNB
+                WORK(U11+I,J)=ZERO
+             END DO
+             DO J=1,I-1
+                WORK(U11+I,J)=A(CUT+I,CUT+J)
+             END DO
+           END DO
+*
+*          invD*L21
+*
+           I=N-CUT-NNB
+           DO WHILE (I .GE. 1)
+             IF (IPIV(CUT+NNB+I) > 0) THEN
+                DO J=1,NNB
+                    WORK(I,J)=WORK(CUT+NNB+I,INVD)*WORK(I,J)
+                END DO
+                I=I-1
+             ELSE
+                DO J=1,NNB
+                   U01_I_J = WORK(I,J)
+                   U01_IP1_J = WORK(I-1,J)
+                   WORK(I,J)=WORK(CUT+NNB+I,INVD)*U01_I_J+
+     $                        WORK(CUT+NNB+I,INVD+1)*U01_IP1_J
+                   WORK(I-1,J)=WORK(CUT+NNB+I-1,INVD+1)*U01_I_J+
+     $                        WORK(CUT+NNB+I-1,INVD)*U01_IP1_J
+                END DO
+                I=I-2
+             END IF
+           END DO
+*
+*        invD1*L11
+*
+           I=NNB
+           DO WHILE (I .GE. 1)
+             IF (IPIV(CUT+I) > 0) THEN
+                DO J=1,NNB
+                    WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J)
+                END DO
+                I=I-1
+             ELSE
+                DO J=1,NNB
+                   U11_I_J = WORK(U11+I,J)
+                   U11_IP1_J = WORK(U11+I-1,J)
+                WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J) +
+     $                      WORK(CUT+I,INVD+1)*U11_IP1_J
+                WORK(U11+I-1,J)=WORK(CUT+I-1,INVD+1)*U11_I_J+
+     $                      WORK(CUT+I-1,INVD)*U11_IP1_J
+                END DO
+                I=I-2
+             END IF
+           END DO
+*    
+*       L11T*invD1*L11->L11
+*
+        CALL CTRMM('L',UPLO,'C','U',NNB, NNB,
+     $             CONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1)
+
+        IF ( (CUT+NNB) .LT. N ) THEN
+*
+*          L21T*invD2*L21->A(CUT+I,CUT+J)
+*
+         CALL CGEMM('C','N',NNB,NNB,N-NNB-CUT,CONE,A(CUT+NNB+1,CUT+1)
+     $             ,LDA,WORK,N+NB+1, ZERO, A(CUT+1,CUT+1), LDA)
+       
+*
+*        L11 =  L11T*invD1*L11 + U01'invD*U01
+*
+         DO I=1,NNB
+            DO J=1,I
+              A(CUT+I,CUT+J)=A(CUT+I,CUT+J)+WORK(U11+I,J)
+            END DO
+         END DO
+*
+*        L01 =  L22T*invD2*L21
+*
+         CALL CTRMM('L',UPLO,'C','U', N-NNB-CUT, NNB,
+     $             CONE,A(CUT+NNB+1,CUT+NNB+1),LDA,WORK,N+NB+1)
+
+*      Update L21
+         DO I=1,N-CUT-NNB
+           DO J=1,NNB
+              A(CUT+NNB+I,CUT+J)=WORK(I,J)
+           END DO
+         END DO
+       ELSE
+*
+*        L11 =  L11T*invD1*L11
+*
+         DO I=1,NNB
+            DO J=1,I
+              A(CUT+I,CUT+J)=WORK(U11+I,J)
+            END DO
+         END DO
+       END IF
+*
+*      Next Block
+*
+           CUT=CUT+NNB
+       END DO
+*
+*        Apply PERMUTATIONS P and P': P * inv(U')*inv(D)*inv(U) *P'
+* 
+            I=N
+            DO WHILE ( I .GE. 1 )
+               IF( IPIV(I) .GT. 0 ) THEN
+                  IP=IPIV(I)
+                 IF (I .LT. IP) CALL CHESWAPR( UPLO, N, A, I ,IP  )
+                 IF (I .GT. IP) CALL CHESWAPR( UPLO, N, A, IP ,I )
+               ELSE
+                 IP=-IPIV(I)
+                 IF ( I .LT. IP) CALL CHESWAPR( UPLO, N, A, I ,IP )
+                 IF ( I .GT. IP) CALL CHESWAPR(  UPLO, N, A, IP ,I )
+                 I=I-1
+               ENDIF
+               I=I-1
+            END DO
+      END IF
+*
+      RETURN
+*
+*     End of CHETRI2X
+*
+      END
+
diff --git a/SRC/zheswapr.f b/SRC/zheswapr.f
new file mode 100644
index 00000000..0cead401
--- /dev/null
+++ b/SRC/zheswapr.f
@@ -0,0 +1,136 @@
+      SUBROUTINE ZHESWAPR( UPLO, N, A, I1, I2)
+*
+*  -- LAPACK auxiliary routine (version 3.3.0) --
+*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*     November 2010
+*
+*     .. Scalar Arguments ..
+      CHARACTER        UPLO
+      INTEGER          I1, I2, N
+*     ..
+*     .. Array Arguments ..
+      COMPLEX*16          A(N,N)
+*
+*  Purpose
+*  =======
+*
+*  ZHESWAPR applies an elementary permutation on the rows and the columns of
+*  a hermitian matrix.
+*
+*  Arguments
+*  =========
+*
+*  UPLO    (input) CHARACTER*1
+*          Specifies whether the details of the factorization are stored
+*          as an upper or lower triangular matrix.
+*          = 'U':  Upper triangular, form is A = U*D*U**T;
+*          = 'L':  Lower triangular, form is A = L*D*L**T.
+*
+*  N       (input) INTEGER
+*          The order of the matrix A.  N >= 0.
+*
+*  A       (input/output) COMPLEX*16 array, dimension (LDA,N)
+*          On entry, the NB diagonal matrix D and the multipliers
+*          used to obtain the factor U or L as computed by CSYTRF.
+*
+*          On exit, if INFO = 0, the (symmetric) inverse of the original
+*          matrix.  If UPLO = 'U', the upper triangular part of the
+*          inverse is formed and the part of A below the diagonal is not
+*          referenced; if UPLO = 'L' the lower triangular part of the
+*          inverse is formed and the part of A above the diagonal is
+*          not referenced.
+*
+*  I1      (input) INTEGER
+*          Index of the first row to swap
+*
+*  I2      (input) INTEGER
+*          Index of the second row to swap
+*
+*  =====================================================================
+*
+*     ..
+*     .. Local Scalars ..
+      LOGICAL            UPPER
+      INTEGER            I
+      COMPLEX*16            TMP
+*
+*     .. External Functions ..
+      LOGICAL            LSAME
+      EXTERNAL           LSAME
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           ZSWAP
+*     ..
+*     .. Executable Statements ..
+*
+      UPPER = LSAME( UPLO, 'U' )
+      IF (UPPER) THEN
+*
+*         UPPER
+*         first swap
+*          - swap column I1 and I2 from I1 to I1-1 
+         CALL ZSWAP( I1-1, A(1,I1), 1, A(1,I2), 1 )
+*
+*          second swap :
+*          - swap A(I1,I1) and A(I2,I2)
+*          - swap row I1 from I1+1 to I2-1 with col I2 from I1+1 to I2-1
+*          - swap A(I2,I1) and A(I1,I2)
+     
+         TMP=A(I1,I1)
+         A(I1,I1)=A(I2,I2)
+         A(I2,I2)=TMP
+*
+         DO I=1,I2-I1-1
+            TMP=A(I1,I1+I)
+            A(I1,I1+I)=DCONJG(A(I1+I,I2))
+            A(I1+I,I2)=DCONJG(TMP)
+         END DO
+*
+          A(I1,I2)=DCONJG(A(I1,I2))
+
+*
+*          third swap
+*          - swap row I1 and I2 from I2+1 to N
+         DO I=I2+1,N
+            TMP=A(I1,I)
+            A(I1,I)=A(I2,I)
+            A(I2,I)=TMP
+         END DO
+*
+        ELSE
+*
+*         LOWER
+*         first swap
+*          - swap row I1 and I2 from 1 to I1-1 
+         CALL ZSWAP ( I1-1, A(I1,1), N, A(I2,1), N )
+*
+*         second swap :
+*          - swap A(I1,I1) and A(I2,I2)
+*          - swap col I1 from I1+1 to I2-1 with row I2 from I1+1 to I2-1     
+*          - swap A(I2,I1) and A(I1,I2)
+
+          TMP=A(I1,I1)
+          A(I1,I1)=A(I2,I2)
+          A(I2,I2)=TMP
+*
+          DO I=1,I2-I1-1
+             TMP=A(I1+I,I1)
+             A(I1+I,I1)=DCONJG(A(I2,I1+I))
+             A(I2,I1+I)=DCONJG(TMP)
+          END DO
+*
+          A(I2,I1)=DCONJG(A(I2,I1))
+*
+*         third swap
+*          - swap col I1 and I2 from I2+1 to N
+          DO I=I2+1,N
+             TMP=A(I,I1)
+             A(I,I1)=A(I,I2)
+             A(I,I2)=TMP
+          END DO
+*
+      ENDIF
+      
+      END SUBROUTINE ZHESWAPR
+
diff --git a/SRC/zhetri2.f b/SRC/zhetri2.f
new file mode 100644
index 00000000..7938b898
--- /dev/null
+++ b/SRC/zhetri2.f
@@ -0,0 +1,127 @@
+      SUBROUTINE ZHETRI2( UPLO, N, A, LDA, IPIV, WORK, LWORK, INFO )
+*
+*  -- LAPACK routine (version 3.3.0) --
+*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*     November 2010
+*
+*  -- Written by Julie Langou of the Univ. of TN    --
+*
+*     .. Scalar Arguments ..
+      CHARACTER          UPLO
+      INTEGER            INFO, LDA, LWORK, N
+*     ..
+*     .. Array Arguments ..
+      INTEGER            IPIV( * )
+      COMPLEX*16         A( LDA, * ), WORK( * )
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  ZHETRI2 computes the inverse of a COMPLEX*16 hermitian indefinite matrix
+*  A using the factorization A = U*D*U**T or A = L*D*L**T computed by
+*  ZHETRF. ZHETRI2 set the LEADING DIMENSION of the workspace
+*  before calling ZHETRI2X that actually compute the inverse.
+*
+*  Arguments
+*  =========
+*
+*  UPLO    (input) CHARACTER*1
+*          Specifies whether the details of the factorization are stored
+*          as an upper or lower triangular matrix.
+*          = 'U':  Upper triangular, form is A = U*D*U**T;
+*          = 'L':  Lower triangular, form is A = L*D*L**T.
+*
+*  N       (input) INTEGER
+*          The order of the matrix A.  N >= 0.
+*
+*  A       (input/output) COMPLEX*16 array, dimension (LDA,N)
+*          On entry, the NB diagonal matrix D and the multipliers
+*          used to obtain the factor U or L as computed by ZHETRF.
+*
+*          On exit, if INFO = 0, the (symmetric) inverse of the original
+*          matrix.  If UPLO = 'U', the upper triangular part of the
+*          inverse is formed and the part of A below the diagonal is not
+*          referenced; if UPLO = 'L' the lower triangular part of the
+*          inverse is formed and the part of A above the diagonal is
+*          not referenced.
+*
+*  LDA     (input) INTEGER
+*          The leading dimension of the array A.  LDA >= max(1,N).
+*
+*  IPIV    (input) INTEGER array, dimension (N)
+*          Details of the interchanges and the NB structure of D
+*          as determined by ZHETRF.
+*
+*  WORK    (workspace) COMPLEX*16 array, dimension (N+NB+1)*(NB+3)
+*
+*  LWORK   (input) INTEGER
+*          The dimension of the array WORK.
+*          WORK is size >= (N+NB+1)*(NB+3)
+*          If LDWORK = -1, then a workspace query is assumed; the routine
+*           calculates:
+*              - the optimal size of the WORK array, returns
+*          this value as the first entry of the WORK array,
+*              - and no error message related to LDWORK is issued by XERBLA.
+*
+*  INFO    (output) INTEGER
+*          = 0: successful exit
+*          < 0: if INFO = -i, the i-th argument had an illegal value
+*          > 0: if INFO = i, D(i,i) = 0; the matrix is singular and its
+*               inverse could not be computed.
+*
+*  =====================================================================
+*
+*     .. Local Scalars ..
+      LOGICAL            UPPER, LQUERY
+      INTEGER            MINSIZE, NBMAX
+*     ..
+*     .. External Functions ..
+      LOGICAL            LSAME
+      INTEGER            ILAENV
+      EXTERNAL           LSAME, ILAENV
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           ZHETRI2X
+*     ..
+*     .. Executable Statements ..
+*
+*     Test the input parameters.
+*
+      INFO = 0
+      UPPER = LSAME( UPLO, 'U' )
+      LQUERY = ( LWORK.EQ.-1 )
+*     Get blocksize
+      NBMAX = ILAENV( 1, 'ZHETRF', UPLO, N, -1, -1, -1 )
+      MINSIZE = (N+NBMAX+1)*(NBMAX+3)
+*
+      IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
+         INFO = -1
+      ELSE IF( N.LT.0 ) THEN
+         INFO = -2
+      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
+         INFO = -4
+      ELSE IF (LWORK .LT. MINSIZE .AND. .NOT.LQUERY ) THEN
+         INFO = -7
+      END IF
+*
+*     Quick return if possible
+*
+*
+      IF( INFO.NE.0 ) THEN
+         CALL XERBLA( 'ZHETRI2', -INFO )
+         RETURN
+      ELSE IF( LQUERY ) THEN
+         WORK(1)=(N+NBMAX+1)*(NBMAX+3)
+         RETURN
+      END IF
+      IF( N.EQ.0 )
+     $   RETURN
+      
+      CALL ZHETRI2X( UPLO, N, A, LDA, IPIV, WORK, NBMAX, INFO )
+      RETURN
+*
+*     End of ZHETRI2
+*
+      END
diff --git a/SRC/zhetri2x.f b/SRC/zhetri2x.f
new file mode 100644
index 00000000..0d3f489e
--- /dev/null
+++ b/SRC/zhetri2x.f
@@ -0,0 +1,508 @@
+      SUBROUTINE ZHETRI2X( UPLO, N, A, LDA, IPIV, WORK, NB, INFO )
+*
+*  -- LAPACK routine (version 3.3.0) --
+*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*     November 2010
+*
+*  -- Written by Julie Langou of the Univ. of TN    --
+*
+*     .. Scalar Arguments ..
+      CHARACTER          UPLO
+      INTEGER            INFO, LDA, N, NB
+*     ..
+*     .. Array Arguments ..
+      INTEGER            IPIV( * )
+      COMPLEX*16            A( LDA, * ), WORK( N+NB+1,* )
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  ZHETRI2X computes the inverse of a COMPLEX*16 Hermitian indefinite matrix
+*  A using the factorization A = U*D*U**H or A = L*D*L**H computed by
+*  ZHETRF.
+*
+*  Arguments
+*  =========
+*
+*  UPLO    (input) CHARACTER*1
+*          Specifies whether the details of the factorization are stored
+*          as an upper or lower triangular matrix.
+*          = 'U':  Upper triangular, form is A = U*D*U**H;
+*          = 'L':  Lower triangular, form is A = L*D*L**H.
+*
+*  N       (input) INTEGER
+*          The order of the matrix A.  N >= 0.
+*
+*  A       (input/output) COMPLEX*16 array, dimension (LDA,N)
+*          On entry, the NNB diagonal matrix D and the multipliers
+*          used to obtain the factor U or L as computed by ZHETRF.
+*
+*          On exit, if INFO = 0, the (symmetric) inverse of the original
+*          matrix.  If UPLO = 'U', the upper triangular part of the
+*          inverse is formed and the part of A below the diagonal is not
+*          referenced; if UPLO = 'L' the lower triangular part of the
+*          inverse is formed and the part of A above the diagonal is
+*          not referenced.
+*
+*  LDA     (input) INTEGER
+*          The leading dimension of the array A.  LDA >= max(1,N).
+*
+*  IPIV    (input) INTEGER array, dimension (N)
+*          Details of the interchanges and the NNB structure of D
+*          as determined by ZHETRF.
+*
+*  WORK    (workspace) COMPLEX*16 array, dimension (N+NNB+1,NNB+3)
+*
+*  NB      (input) INTEGER
+*          Block size
+*
+*  INFO    (output) INTEGER
+*          = 0: successful exit
+*          < 0: if INFO = -i, the i-th argument had an illegal value
+*          > 0: if INFO = i, D(i,i) = 0; the matrix is singular and its
+*               inverse could not be computed.
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      REAL               ONE
+      COMPLEX*16            CONE, ZERO
+      PARAMETER          ( ONE = 1.0D+0,
+     $                   CONE = ( 1.0D+0, 0.0D+0 ),
+     $                   ZERO = ( 0.0D+0, 0.0D+0 ) )
+*     ..
+*     .. Local Scalars ..
+      LOGICAL            UPPER
+      INTEGER            I, IINFO, IP, K, CUT, NNB
+      INTEGER            COUNT
+      INTEGER            J, U11, INVD
+
+      COMPLEX*16   AK, AKKP1, AKP1, D, T
+      COMPLEX*16   U01_I_J, U01_IP1_J
+      COMPLEX*16   U11_I_J, U11_IP1_J
+*     ..
+*     .. External Functions ..
+      LOGICAL            LSAME
+      EXTERNAL           LSAME
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           ZSYCONV, XERBLA, ZTRTRI
+      EXTERNAL           ZGEMM, ZTRMM, ZHESWAPR
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          MAX
+*     ..
+*     .. Executable Statements ..
+*
+*     Test the input parameters.
+*
+      INFO = 0
+      UPPER = LSAME( UPLO, 'U' )
+      IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
+         INFO = -1
+      ELSE IF( N.LT.0 ) THEN
+         INFO = -2
+      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
+         INFO = -4
+      END IF
+*
+*     Quick return if possible
+*
+*
+      IF( INFO.NE.0 ) THEN
+         CALL XERBLA( 'ZHETRI2X', -INFO )
+         RETURN
+      END IF
+      IF( N.EQ.0 )
+     $   RETURN
+*
+*     Convert A
+*     Workspace got Non-diag elements of D
+*
+      CALL ZSYCONV( UPLO, 'C', N, A, LDA, IPIV, WORK, IINFO )
+*
+*     Check that the diagonal matrix D is nonsingular.
+*
+      IF( UPPER ) THEN
+*
+*        Upper triangular storage: examine D from bottom to top
+*
+         DO INFO = N, 1, -1
+            IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO )
+     $         RETURN
+         END DO
+      ELSE
+*
+*        Lower triangular storage: examine D from top to bottom.
+*
+         DO INFO = 1, N
+            IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO )
+     $         RETURN
+         END DO
+      END IF
+      INFO = 0
+*
+*  Splitting Workspace
+*     U01 is a block (N,NB+1) 
+*     The first element of U01 is in WORK(1,1)
+*     U11 is a block (NB+1,NB+1)
+*     The first element of U11 is in WORK(N+1,1)
+      U11 = N
+*     INVD is a block (N,2)
+*     The first element of INVD is in WORK(1,INVD)
+      INVD = NB+2
+
+      IF( UPPER ) THEN
+*
+*        invA = P * inv(U')*inv(D)*inv(U)*P'.
+*
+        CALL ZTRTRI( UPLO, 'U', N, A, LDA, INFO )
+*
+*       inv(D) and inv(D)*inv(U)
+* 
+        K=1
+        DO WHILE ( K .LE. N )
+         IF( IPIV( K ).GT.0 ) THEN
+*           1 x 1 diagonal NNB
+             WORK(K,INVD) = ONE / REAL ( A( K, K ) )
+             WORK(K,INVD+1) = 0
+            K=K+1
+         ELSE
+*           2 x 2 diagonal NNB
+             T = ABS ( WORK(K+1,1) )
+             AK = REAL ( A( K, K ) ) / T
+             AKP1 = REAL ( A( K+1, K+1 ) ) / T
+             AKKP1 = WORK(K+1,1)  / T
+             D = T*( AK*AKP1-ONE )
+             WORK(K,INVD) = AKP1 / D
+             WORK(K+1,INVD+1) = AK / D
+             WORK(K,INVD+1) = -AKKP1 / D  
+             WORK(K+1,INVD) = DCONJG (WORK(K,INVD+1) )
+            K=K+2
+         END IF
+        END DO
+*
+*       inv(U') = (inv(U))'
+*
+*       inv(U')*inv(D)*inv(U)
+*
+        CUT=N
+        DO WHILE (CUT .GT. 0)
+           NNB=NB
+           IF (CUT .LE. NNB) THEN
+              NNB=CUT
+           ELSE
+              COUNT = 0
+*             count negative elements, 
+              DO I=CUT+1-NNB,CUT
+                  IF (IPIV(I) .LT. 0) COUNT=COUNT+1
+              END DO
+*             need a even number for a clear cut
+              IF (MOD(COUNT,2) .EQ. 1) NNB=NNB+1
+           END IF
+
+           CUT=CUT-NNB
+*
+*          U01 Block 
+*
+           DO I=1,CUT
+             DO J=1,NNB
+              WORK(I,J)=A(I,CUT+J)
+             END DO
+           END DO
+*
+*          U11 Block
+*
+           DO I=1,NNB
+             WORK(U11+I,I)=CONE
+             DO J=1,I-1
+                WORK(U11+I,J)=ZERO
+             END DO
+             DO J=I+1,NNB
+                WORK(U11+I,J)=A(CUT+I,CUT+J)
+             END DO
+           END DO
+*
+*          invD*U01
+*
+           I=1
+           DO WHILE (I .LE. CUT)
+             IF (IPIV(I) > 0) THEN
+                DO J=1,NNB
+                    WORK(I,J)=WORK(I,INVD)*WORK(I,J)
+                END DO
+                I=I+1
+             ELSE
+                DO J=1,NNB
+                   U01_I_J = WORK(I,J)
+                   U01_IP1_J = WORK(I+1,J)
+                   WORK(I,J)=WORK(I,INVD)*U01_I_J+
+     $                      WORK(I,INVD+1)*U01_IP1_J
+                   WORK(I+1,J)=WORK(I+1,INVD)*U01_I_J+
+     $                      WORK(I+1,INVD+1)*U01_IP1_J
+                END DO
+                I=I+2
+             END IF
+           END DO
+*
+*        invD1*U11
+*
+           I=1
+           DO WHILE (I .LE. NNB)
+             IF (IPIV(CUT+I) > 0) THEN
+                DO J=I,NNB
+                    WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J)
+                END DO
+                I=I+1
+             ELSE
+                DO J=I,NNB
+                   U11_I_J = WORK(U11+I,J)
+                   U11_IP1_J = WORK(U11+I+1,J)
+                WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J) +
+     $                      WORK(CUT+I,INVD+1)*WORK(U11+I+1,J)
+                WORK(U11+I+1,J)=WORK(CUT+I+1,INVD)*U11_I_J+
+     $                      WORK(CUT+I+1,INVD+1)*U11_IP1_J
+                END DO
+                I=I+2
+             END IF
+           END DO
+*    
+*       U11T*invD1*U11->U11
+*
+        CALL ZTRMM('L','U','C','U',NNB, NNB,
+     $             CONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1)
+*
+*          U01'invD*U01->A(CUT+I,CUT+J)
+*
+         CALL ZGEMM('C','N',NNB,NNB,CUT,CONE,A(1,CUT+1),LDA,
+     $              WORK,N+NB+1, ZERO, A(CUT+1,CUT+1), LDA)
+*
+*        U11 =  U11T*invD1*U11 + U01'invD*U01
+*
+         DO I=1,NNB
+            DO J=I,NNB
+              A(CUT+I,CUT+J)=A(CUT+I,CUT+J)+WORK(U11+I,J)
+            END DO
+         END DO
+*
+*        U01 =  U00T*invD0*U01
+*
+         CALL ZTRMM('L',UPLO,'C','U',CUT, NNB,
+     $             CONE,A,LDA,WORK,N+NB+1)
+
+*
+*        Update U01
+*
+         DO I=1,CUT
+           DO J=1,NNB
+            A(I,CUT+J)=WORK(I,J)
+           END DO
+         END DO
+*
+*      Next Block
+*
+       END DO
+*
+*        Apply PERMUTATIONS P and P': P * inv(U')*inv(D)*inv(U) *P'
+*  
+            I=1
+            DO WHILE ( I .LE. N )
+               IF( IPIV(I) .GT. 0 ) THEN
+                  IP=IPIV(I)
+                 IF (I .LT. IP) CALL ZHESWAPR( UPLO, N, A, I ,IP )
+                 IF (I .GT. IP) CALL ZHESWAPR( UPLO, N, A, IP ,I )
+               ELSE
+                 IP=-IPIV(I)
+                 I=I+1
+                 IF ( (I-1) .LT. IP) 
+     $                  CALL ZHESWAPR( UPLO, N, A, I-1 ,IP )
+                 IF ( (I-1) .GT. IP) 
+     $                  CALL ZHESWAPR( UPLO, N, A, IP ,I-1 )
+              ENDIF
+               I=I+1
+            END DO
+      ELSE
+*
+*        LOWER...
+*
+*        invA = P * inv(U')*inv(D)*inv(U)*P'.
+*
+         CALL ZTRTRI( UPLO, 'U', N, A, LDA, INFO )
+*
+*       inv(D) and inv(D)*inv(U)
+* 
+        K=N
+        DO WHILE ( K .GE. 1 )
+         IF( IPIV( K ).GT.0 ) THEN
+*           1 x 1 diagonal NNB
+             WORK(K,INVD) = ONE / REAL ( A( K, K ) )
+             WORK(K,INVD+1) = 0
+            K=K-1
+         ELSE
+*           2 x 2 diagonal NNB
+             T = ABS ( WORK(K-1,1) )
+             AK = REAL ( A( K-1, K-1 ) ) / T
+             AKP1 = REAL ( A( K, K ) ) / T
+             AKKP1 = WORK(K-1,1) / T
+             D = T*( AK*AKP1-ONE )
+             WORK(K-1,INVD) = AKP1 / D
+             WORK(K,INVD) = AK / D
+             WORK(K,INVD+1) = -AKKP1 / D  
+             WORK(K-1,INVD+1) = DCONJG (WORK(K,INVD+1) )
+            K=K-2
+         END IF
+        END DO
+*
+*       inv(U') = (inv(U))'
+*
+*       inv(U')*inv(D)*inv(U)
+*
+        CUT=0
+        DO WHILE (CUT .LT. N)
+           NNB=NB
+           IF (CUT + NNB .GE. N) THEN
+              NNB=N-CUT
+           ELSE
+              COUNT = 0
+*             count negative elements, 
+              DO I=CUT+1,CUT+NNB
+                  IF (IPIV(I) .LT. 0) COUNT=COUNT+1
+              END DO
+*             need a even number for a clear cut
+              IF (MOD(COUNT,2) .EQ. 1) NNB=NNB+1
+           END IF
+*      L21 Block
+           DO I=1,N-CUT-NNB
+             DO J=1,NNB
+              WORK(I,J)=A(CUT+NNB+I,CUT+J)
+             END DO
+           END DO
+*     L11 Block
+           DO I=1,NNB
+             WORK(U11+I,I)=CONE
+             DO J=I+1,NNB
+                WORK(U11+I,J)=ZERO
+             END DO
+             DO J=1,I-1
+                WORK(U11+I,J)=A(CUT+I,CUT+J)
+             END DO
+           END DO
+*
+*          invD*L21
+*
+           I=N-CUT-NNB
+           DO WHILE (I .GE. 1)
+             IF (IPIV(CUT+NNB+I) > 0) THEN
+                DO J=1,NNB
+                    WORK(I,J)=WORK(CUT+NNB+I,INVD)*WORK(I,J)
+                END DO
+                I=I-1
+             ELSE
+                DO J=1,NNB
+                   U01_I_J = WORK(I,J)
+                   U01_IP1_J = WORK(I-1,J)
+                   WORK(I,J)=WORK(CUT+NNB+I,INVD)*U01_I_J+
+     $                        WORK(CUT+NNB+I,INVD+1)*U01_IP1_J
+                   WORK(I-1,J)=WORK(CUT+NNB+I-1,INVD+1)*U01_I_J+
+     $                        WORK(CUT+NNB+I-1,INVD)*U01_IP1_J
+                END DO
+                I=I-2
+             END IF
+           END DO
+*
+*        invD1*L11
+*
+           I=NNB
+           DO WHILE (I .GE. 1)
+             IF (IPIV(CUT+I) > 0) THEN
+                DO J=1,NNB
+                    WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J)
+                END DO
+                I=I-1
+             ELSE
+                DO J=1,NNB
+                   U11_I_J = WORK(U11+I,J)
+                   U11_IP1_J = WORK(U11+I-1,J)
+                WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J) +
+     $                      WORK(CUT+I,INVD+1)*U11_IP1_J
+                WORK(U11+I-1,J)=WORK(CUT+I-1,INVD+1)*U11_I_J+
+     $                      WORK(CUT+I-1,INVD)*U11_IP1_J
+                END DO
+                I=I-2
+             END IF
+           END DO
+*    
+*       L11T*invD1*L11->L11
+*
+        CALL ZTRMM('L',UPLO,'C','U',NNB, NNB,
+     $             CONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1)
+
+        IF ( (CUT+NNB) .LT. N ) THEN
+*
+*          L21T*invD2*L21->A(CUT+I,CUT+J)
+*
+         CALL ZGEMM('C','N',NNB,NNB,N-NNB-CUT,CONE,A(CUT+NNB+1,CUT+1)
+     $             ,LDA,WORK,N+NB+1, ZERO, A(CUT+1,CUT+1), LDA)
+       
+*
+*        L11 =  L11T*invD1*L11 + U01'invD*U01
+*
+         DO I=1,NNB
+            DO J=1,I
+              A(CUT+I,CUT+J)=A(CUT+I,CUT+J)+WORK(U11+I,J)
+            END DO
+         END DO
+*
+*        L01 =  L22T*invD2*L21
+*
+         CALL ZTRMM('L',UPLO,'C','U', N-NNB-CUT, NNB,
+     $             CONE,A(CUT+NNB+1,CUT+NNB+1),LDA,WORK,N+NB+1)
+
+*      Update L21
+         DO I=1,N-CUT-NNB
+           DO J=1,NNB
+              A(CUT+NNB+I,CUT+J)=WORK(I,J)
+           END DO
+         END DO
+       ELSE
+*
+*        L11 =  L11T*invD1*L11
+*
+         DO I=1,NNB
+            DO J=1,I
+              A(CUT+I,CUT+J)=WORK(U11+I,J)
+            END DO
+         END DO
+       END IF
+*
+*      Next Block
+*
+           CUT=CUT+NNB
+       END DO
+*
+*        Apply PERMUTATIONS P and P': P * inv(U')*inv(D)*inv(U) *P'
+* 
+            I=N
+            DO WHILE ( I .GE. 1 )
+               IF( IPIV(I) .GT. 0 ) THEN
+                  IP=IPIV(I)
+                 IF (I .LT. IP) CALL ZHESWAPR( UPLO, N, A, I ,IP  )
+                 IF (I .GT. IP) CALL ZHESWAPR( UPLO, N, A, IP ,I )
+               ELSE
+                 IP=-IPIV(I)
+                 IF ( I .LT. IP) CALL ZHESWAPR( UPLO, N, A, I ,IP )
+                 IF ( I .GT. IP) CALL ZHESWAPR(  UPLO, N, A, IP ,I )
+                 I=I-1
+               ENDIF
+               I=I-1
+            END DO
+      END IF
+*
+      RETURN
+*
+*     End of ZHETRI2X
+*
+      END
+
diff --git a/TESTING/LIN/cchkhe.f b/TESTING/LIN/cchkhe.f
index 1f0b6db5..aff1883d 100644
--- a/TESTING/LIN/cchkhe.f
+++ b/TESTING/LIN/cchkhe.f
@@ -22,7 +22,7 @@
 *  Purpose
 *  =======
 *
-*  CCHKHE tests CHETRF, -TRI, -TRS, -TRS2, -RFS, and -CON.
+*  CCHKHE tests CHETRF, -TRI2, -TRS, -TRS2, -RFS, and -CON.
 *
 *  Arguments
 *  =========
@@ -116,9 +116,9 @@
 *     ..
 *     .. External Subroutines ..
       EXTERNAL           ALAERH, ALAHD, ALASUM, CERRHE, CGET04, CHECON,
-     $                   CHERFS, CHET01, CHETRF, CHETRI, CHETRS, CLACPY,
-     $                   CLAIPD, CLARHS, CLATB4, CLATMS, CPOT02, CPOT03,
-     $                   CPOT05, XLAENV
+     $                   CHERFS, CHET01, CHETRF, CHETRI2, CHETRS, 
+     $                   CLACPY, CLAIPD, CLARHS, CLATB4, CLATMS, CPOT02,
+     $                   CPOT03, CPOT05, XLAENV
 *     ..
 *     .. Intrinsic Functions ..
       INTRINSIC          MAX, MIN
@@ -329,9 +329,10 @@
 *
                   IF( INB.EQ.1 .AND. .NOT.TRFCON ) THEN
                      CALL CLACPY( UPLO, N, N, AFAC, LDA, AINV, LDA )
-                     SRNAMT = 'CHETRI'
-                     CALL CHETRI( UPLO, N, AINV, LDA, IWORK, WORK,
-     $                            INFO )
+                     SRNAMT = 'CHETRI2'
+                     LWORK = (N+NB+1)*(NB+3)
+                     CALL CHETRI2( UPLO, N, AINV, LDA, IWORK, WORK,
+     $                            LWORK, INFO )
 *
 *                 Check error code from CHETRI.
 *
diff --git a/TESTING/LIN/cdrvhe.f b/TESTING/LIN/cdrvhe.f
index cba64094..a70a526a 100644
--- a/TESTING/LIN/cdrvhe.f
+++ b/TESTING/LIN/cdrvhe.f
@@ -106,7 +106,7 @@
 *     ..
 *     .. External Subroutines ..
       EXTERNAL           ALADHD, ALAERH, ALASVM, CERRVX, CGET04, CHESV,
-     $                   CHESVX, CHET01, CHETRF, CHETRI, CLACPY, CLAIPD,
+     $                   CHESVX, CHET01, CHETRF, CHETRI2, CLACPY, CLAIPD,
      $                   CLARHS, CLASET, CLATB4, CLATMS, CPOT02, CPOT05,
      $                   XLAENV
 *     ..
@@ -300,8 +300,9 @@
 *                    Compute inv(A) and take its norm.
 *
                      CALL CLACPY( UPLO, N, N, AFAC, LDA, AINV, LDA )
-                     CALL CHETRI( UPLO, N, AINV, LDA, IWORK, WORK,
-     $                            INFO )
+                     LWORK = (N+NB+1)*(NB+3)
+                     CALL CHETRI2( UPLO, N, AINV, LDA, IWORK, WORK,
+     $                            LWORK, INFO )
                      AINVNM = CLANHE( '1', UPLO, N, AINV, LDA, RWORK )
 *
 *                    Compute the 1-norm condition number of A.
diff --git a/TESTING/LIN/cdrvhex.f b/TESTING/LIN/cdrvhex.f
index bf16a294..56aebba8 100644
--- a/TESTING/LIN/cdrvhex.f
+++ b/TESTING/LIN/cdrvhex.f
@@ -112,7 +112,7 @@
 *     ..
 *     .. External Subroutines ..
       EXTERNAL           ALADHD, ALAERH, ALASVM, CERRVX, CGET04, CHESV,
-     $                   CHESVX, CHET01, CHETRF, CHETRI, CLACPY, CLAIPD,
+     $                   CHESVX, CHET01, CHETRF, CHETRI2, CLACPY, CLAIPD,
      $                   CLARHS, CLASET, CLATB4, CLATMS, CPOT02, CPOT05,
      $                   XLAENV, CHESVXX
 *     ..
@@ -306,8 +306,9 @@
 *                    Compute inv(A) and take its norm.
 *
                      CALL CLACPY( UPLO, N, N, AFAC, LDA, AINV, LDA )
-                     CALL CHETRI( UPLO, N, AINV, LDA, IWORK, WORK,
-     $                            INFO )
+                     LWORK = (N+NB+1)*(NB+3)
+                     CALL CHETRI2( UPLO, N, AINV, LDA, IWORK, WORK,
+     $                            LWORK, INFO )
                      AINVNM = CLANHE( '1', UPLO, N, AINV, LDA, RWORK )
 *
 *                    Compute the 1-norm condition number of A.
diff --git a/TESTING/LIN/cerrhe.f b/TESTING/LIN/cerrhe.f
index c58cd36e..a9b616f0 100644
--- a/TESTING/LIN/cerrhe.f
+++ b/TESTING/LIN/cerrhe.f
@@ -48,8 +48,8 @@
 *     ..
 *     .. External Subroutines ..
       EXTERNAL           ALAESM, CHECON, CHERFS, CHETF2, CHETRF, CHETRI,
-     $                   CHETRS, CHKXER, CHPCON, CHPRFS, CHPTRF, CHPTRI,
-     $                   CHPTRS
+     $                   CHETRI2, CHETRS, CHKXER, CHPCON, CHPRFS, CHPTRF,
+     $                   CHPTRI, CHPTRS
 *     ..
 *     .. Scalars in Common ..
       LOGICAL            LERR, OK
@@ -130,6 +130,19 @@
          CALL CHETRI( 'U', 2, A, 1, IP, W, INFO )
          CALL CHKXER( 'CHETRI', INFOT, NOUT, LERR, OK )
 *
+*        CHETRI2
+*
+         SRNAMT = 'CHETRI2'
+         INFOT = 1
+         CALL CHETRI2( '/', 0, A, 1, IP, W, 1, INFO )
+         CALL CHKXER( 'CHETRI2', INFOT, NOUT, LERR, OK )
+         INFOT = 2
+         CALL CHETRI2( 'U', -1, A, 1, IP, W, 1, INFO )
+         CALL CHKXER( 'CHETRI2', INFOT, NOUT, LERR, OK )
+         INFOT = 4
+         CALL CHETRI2( 'U', 2, A, 1, IP, W, 1, INFO )
+         CALL CHKXER( 'CHETRI2', INFOT, NOUT, LERR, OK )
+*
 *        CHETRS
 *
          SRNAMT = 'CHETRS'
diff --git a/TESTING/LIN/cerrhex.f b/TESTING/LIN/cerrhex.f
index e6bd10e5..a159abd8 100644
--- a/TESTING/LIN/cerrhex.f
+++ b/TESTING/LIN/cerrhex.f
@@ -54,8 +54,8 @@
 *     ..
 *     .. External Subroutines ..
       EXTERNAL           ALAESM, CHECON, CHERFS, CHETF2, CHETRF, CHETRI,
-     $                   CHETRS, CHKXER, CHPCON, CHPRFS, CHPTRF, CHPTRI,
-     $                   CHPTRS, CHERFSX
+     $                   CHETRI2, CHETRS, CHKXER, CHPCON, CHPRFS, CHPTRF,
+     $                   CHPTRI, CHPTRS, CHERFSX
 *     ..
 *     .. Scalars in Common ..
       LOGICAL            LERR, OK
@@ -137,6 +137,19 @@
          CALL CHETRI( 'U', 2, A, 1, IP, W, INFO )
          CALL CHKXER( 'CHETRI', INFOT, NOUT, LERR, OK )
 *
+*        CHETRI2
+*
+         SRNAMT = 'CHETRI2'
+         INFOT = 1
+         CALL CHETRI2( '/', 0, A, 1, IP, W, 1, INFO )
+         CALL CHKXER( 'CHETRI2', INFOT, NOUT, LERR, OK )
+         INFOT = 2
+         CALL CHETRI2( 'U', -1, A, 1, IP, W, 1, INFO )
+         CALL CHKXER( 'CHETRI2', INFOT, NOUT, LERR, OK )
+         INFOT = 4
+         CALL CHETRI2( 'U', 2, A, 1, IP, W, 1, INFO )
+         CALL CHKXER( 'CHETRI2', INFOT, NOUT, LERR, OK )
+*
 *        CHETRS
 *
          SRNAMT = 'CHETRS'
diff --git a/TESTING/LIN/zchkhe.f b/TESTING/LIN/zchkhe.f
index 1d04bbd4..9dd87c70 100644
--- a/TESTING/LIN/zchkhe.f
+++ b/TESTING/LIN/zchkhe.f
@@ -22,7 +22,7 @@
 *  Purpose
 *  =======
 *
-*  ZCHKHE tests ZHETRF, -TRI, -TRS, -TRS2, -RFS, and -CON.
+*  ZCHKHE tests ZHETRF, -TRI2, -TRS, -TRS2, -RFS, and -CON.
 *
 *  Arguments
 *  =========
@@ -116,7 +116,7 @@
 *     ..
 *     .. External Subroutines ..
       EXTERNAL           ALAERH, ALAHD, ALASUM, XLAENV, ZERRHE, ZGET04,
-     $                   ZHECON, ZHERFS, ZHET01, ZHETRF, ZHETRI, ZHETRS,
+     $                   ZHECON, ZHERFS, ZHET01, ZHETRF, ZHETRI2, ZHETRS,
      $                   ZLACPY, ZLAIPD, ZLARHS, ZLATB4, ZLATMS, ZPOT02,
      $                   ZPOT03, ZPOT05
 *     ..
@@ -329,9 +329,10 @@
 *
                   IF( INB.EQ.1 .AND. .NOT.TRFCON ) THEN
                      CALL ZLACPY( UPLO, N, N, AFAC, LDA, AINV, LDA )
-                     SRNAMT = 'ZHETRI'
-                     CALL ZHETRI( UPLO, N, AINV, LDA, IWORK, WORK,
-     $                            INFO )
+                     SRNAMT = 'ZHETRI2'
+                     LWORK = (N+NB+1)*(NB+3)
+                     CALL ZHETRI2( UPLO, N, AINV, LDA, IWORK, WORK,
+     $                            LWORK, INFO )
 *
 *                 Check error code from ZHETRI.
 *
diff --git a/TESTING/LIN/zdrvhe.f b/TESTING/LIN/zdrvhe.f
index fdf70f6d..bf428507 100644
--- a/TESTING/LIN/zdrvhe.f
+++ b/TESTING/LIN/zdrvhe.f
@@ -106,7 +106,7 @@
 *     ..
 *     .. External Subroutines ..
       EXTERNAL           ALADHD, ALAERH, ALASVM, XLAENV, ZERRVX, ZGET04,
-     $                   ZHESV, ZHESVX, ZHET01, ZHETRF, ZHETRI, ZLACPY,
+     $                   ZHESV, ZHESVX, ZHET01, ZHETRF, ZHETRI2, ZLACPY,
      $                   ZLAIPD, ZLARHS, ZLASET, ZLATB4, ZLATMS, ZPOT02,
      $                   ZPOT05
 *     ..
@@ -300,8 +300,9 @@
 *                    Compute inv(A) and take its norm.
 *
                      CALL ZLACPY( UPLO, N, N, AFAC, LDA, AINV, LDA )
-                     CALL ZHETRI( UPLO, N, AINV, LDA, IWORK, WORK,
-     $                            INFO )
+                     LWORK = (N+NB+1)*(NB+3)
+                     CALL ZHETRI2( UPLO, N, AINV, LDA, IWORK, WORK,
+     $                            LWORK, INFO )
                      AINVNM = ZLANHE( '1', UPLO, N, AINV, LDA, RWORK )
 *
 *                    Compute the 1-norm condition number of A.
diff --git a/TESTING/LIN/zdrvhex.f b/TESTING/LIN/zdrvhex.f
index d0461376..4ad4738f 100644
--- a/TESTING/LIN/zdrvhex.f
+++ b/TESTING/LIN/zdrvhex.f
@@ -112,7 +112,7 @@
 *     ..
 *     .. External Subroutines ..
       EXTERNAL           ALADHD, ALAERH, ALASVM, XLAENV, ZERRVX, ZGET04,
-     $                   ZHESV, ZHESVX, ZHET01, ZHETRF, ZHETRI, ZLACPY,
+     $                   ZHESV, ZHESVX, ZHET01, ZHETRF, ZHETRI2, ZLACPY,
      $                   ZLAIPD, ZLARHS, ZLASET, ZLATB4, ZLATMS, ZPOT02,
      $                   ZPOT05, ZHESVXX
 *     ..
@@ -306,8 +306,9 @@
 *                    Compute inv(A) and take its norm.
 *
                      CALL ZLACPY( UPLO, N, N, AFAC, LDA, AINV, LDA )
-                     CALL ZHETRI( UPLO, N, AINV, LDA, IWORK, WORK,
-     $                            INFO )
+                     LWORK = (N+NB+1)*(NB+3)
+                     CALL ZHETRI2( UPLO, N, AINV, LDA, IWORK, WORK,
+     $                            LWORK, INFO )
                      AINVNM = ZLANHE( '1', UPLO, N, AINV, LDA, RWORK )
 *
 *                    Compute the 1-norm condition number of A.
diff --git a/TESTING/LIN/zerrhe.f b/TESTING/LIN/zerrhe.f
index 60980a90..2961640d 100644
--- a/TESTING/LIN/zerrhe.f
+++ b/TESTING/LIN/zerrhe.f
@@ -48,8 +48,8 @@
 *     ..
 *     .. External Subroutines ..
       EXTERNAL           ALAESM, CHKXER, ZHECON, ZHERFS, ZHETF2, ZHETRF,
-     $                   ZHETRI, ZHETRS, ZHPCON, ZHPRFS, ZHPTRF, ZHPTRI,
-     $                   ZHPTRS
+     $                   ZHETRI, ZHETRI2, ZHETRS, ZHPCON, ZHPRFS, ZHPTRF,
+     $                   ZHPTRI, ZHPTRS
 *     ..
 *     .. Scalars in Common ..
       LOGICAL            LERR, OK
@@ -132,6 +132,19 @@
          CALL ZHETRI( 'U', 2, A, 1, IP, W, INFO )
          CALL CHKXER( 'ZHETRI', INFOT, NOUT, LERR, OK )
 *
+*        ZHETRI2
+*
+         SRNAMT = 'ZHETRI2'
+         INFOT = 1
+         CALL ZHETRI2( '/', 0, A, 1, IP, W, 1, INFO )
+         CALL CHKXER( 'ZHETRI2', INFOT, NOUT, LERR, OK )
+         INFOT = 2
+         CALL ZHETRI2( 'U', -1, A, 1, IP, W, 1, INFO )
+         CALL CHKXER( 'ZHETRI2', INFOT, NOUT, LERR, OK )
+         INFOT = 4
+         CALL ZHETRI2( 'U', 2, A, 1, IP, W, 1, INFO )
+         CALL CHKXER( 'ZHETRI2', INFOT, NOUT, LERR, OK )
+*
 *        ZHETRS
 *
          SRNAMT = 'ZHETRS'
diff --git a/TESTING/LIN/zerrhex.f b/TESTING/LIN/zerrhex.f
index 0777423d..674295ec 100644
--- a/TESTING/LIN/zerrhex.f
+++ b/TESTING/LIN/zerrhex.f
@@ -54,8 +54,8 @@
 *     ..
 *     .. External Subroutines ..
       EXTERNAL           ALAESM, CHKXER, ZHECON, ZHERFS, ZHETF2, ZHETRF,
-     $                   ZHETRI, ZHETRS, ZHPCON, ZHPRFS, ZHPTRF, ZHPTRI,
-     $                   ZHPTRS, ZHERFSX
+     $                   ZHETRI, ZHETRI2, ZHETRS, ZHPCON, ZHPRFS, ZHPTRF,
+     $                   ZHPTRI, ZHPTRS, ZHERFSX
 *     ..
 *     .. Scalars in Common ..
       LOGICAL            LERR, OK
@@ -139,6 +139,19 @@
          CALL ZHETRI( 'U', 2, A, 1, IP, W, INFO )
          CALL CHKXER( 'ZHETRI', INFOT, NOUT, LERR, OK )
 *
+*        ZHETRI2
+*
+         SRNAMT = 'ZHETRI2'
+         INFOT = 1
+         CALL ZHETRI2( '/', 0, A, 1, IP, W, 1, INFO )
+         CALL CHKXER( 'ZHETRI2', INFOT, NOUT, LERR, OK )
+         INFOT = 2
+         CALL ZHETRI2( 'U', -1, A, 1, IP, W, 1, INFO )
+         CALL CHKXER( 'ZHETRI2', INFOT, NOUT, LERR, OK )
+         INFOT = 4
+         CALL ZHETRI2( 'U', 2, A, 1, IP, W, 1, INFO )
+         CALL CHKXER( 'ZHETRI2', INFOT, NOUT, LERR, OK )
+*
 *        ZHETRS
 *
          SRNAMT = 'ZHETRS'