Adding Assen for SV, TRF, SV from @iyamazak@icl.utk.edu

Contribution by Ichitaro Yamazaki, University of Tennessee

Note: need to add corresponding LAPACKE routines

25 files changed, 5732 insertions, 37 deletions

diff --git a/TESTING/LIN/CMakeLists.txt b/TESTING/LIN/CMakeLists.txt
index a924b12f..b8320977 100644
--- a/TESTING/LIN/CMakeLists.txt
+++ b/TESTING/LIN/CMakeLists.txt
@@ -10,10 +10,10 @@ set(SLINTST  schkaa.f
    schkeq.f schkgb.f schkge.f schkgt.f
    schklq.f schkpb.f schkpo.f schkps.f schkpp.f
    schkpt.f schkq3.f schkql.f schkqr.f schkrq.f
-   schksp.f schksy.f schksy_rook.f schktb.f schktp.f schktr.f
+   schksp.f schksy.f schksy_rook.f schksy_aasen.f schktb.f schktp.f schktr.f
    schktz.f
    sdrvgt.f sdrvls.f sdrvpb.f
-   sdrvpp.f sdrvpt.f sdrvsp.f sdrvsy.f sdrvsy_rook.f
+   sdrvpp.f sdrvpt.f sdrvsp.f sdrvsy.f sdrvsy_rook.f sdrvsy_aasen.f
    serrgt.f serrlq.f serrls.f
    serrpo.f serrps.f serrql.f serrqp.f serrqr.f
    serrrq.f serrsy.f serrtr.f serrtz.f serrvx.f
@@ -29,7 +29,7 @@ set(SLINTST  schkaa.f
    sqrt01.f sqrt01p.f sqrt02.f sqrt03.f sqrt11.f sqrt12.f
    sqrt13.f sqrt14.f sqrt15.f sqrt16.f sqrt17.f
    srqt01.f srqt02.f srqt03.f srzt01.f srzt02.f
-   sspt01.f ssyt01.f ssyt01_rook.f
+   sspt01.f ssyt01.f ssyt01_rook.f ssyt01_aasen.f 
    stbt02.f stbt03.f stbt05.f stbt06.f stpt01.f
    stpt02.f stpt03.f stpt05.f stpt06.f strt01.f
    strt02.f strt03.f strt05.f strt06.f
@@ -44,11 +44,11 @@ endif()
 
 set(CLINTST  cchkaa.f
    cchkeq.f cchkgb.f cchkge.f cchkgt.f
-   cchkhe.f cchkhe_rook.f cchkhp.f cchklq.f cchkpb.f
+   cchkhe.f cchkhe_rook.f cchkhe_aasen.f cchkhp.f cchklq.f cchkpb.f
    cchkpo.f cchkps.f cchkpp.f cchkpt.f cchkq3.f cchkql.f
    cchkqr.f cchkrq.f cchksp.f cchksy.f cchksy_rook.f cchktb.f
    cchktp.f cchktr.f cchktz.f
-   cdrvgt.f cdrvhe.f cdrvhe_rook.f cdrvhp.f
+   cdrvgt.f cdrvhe.f cdrvhe_rook.f cdrvhe_aasen.f cdrvhp.f
    cdrvls.f cdrvpb.f cdrvpp.f cdrvpt.f
    cdrvsp.f cdrvsy.f cdrvsy_rook.f
    cerrgt.f cerrhe.f cerrlq.f
@@ -58,7 +58,7 @@ set(CLINTST  cchkaa.f
    cgbt01.f cgbt02.f cgbt05.f cgelqs.f cgeqls.f cgeqrs.f
    cgerqs.f cget01.f cget02.f
    cget03.f cget04.f cget07.f cgtt01.f cgtt02.f
-   cgtt05.f chet01.f chet01_rook.f chpt01.f claipd.f claptm.f clarhs.f clatb4.f clatb5.f
+   cgtt05.f chet01.f chet01_rook.f chet01_aasen.f chpt01.f claipd.f claptm.f clarhs.f clatb4.f clatb5.f
    clatsp.f clatsy.f clattb.f clattp.f clattr.f
    clavhe.f clavhe_rook.f clavhp.f clavsp.f clavsy.f clavsy_rook.f clqt01.f
    clqt02.f clqt03.f cpbt01.f cpbt02.f cpbt05.f
@@ -87,10 +87,10 @@ set(DLINTST  dchkaa.f
    dchkeq.f dchkgb.f dchkge.f dchkgt.f
    dchklq.f dchkpb.f dchkpo.f dchkps.f dchkpp.f
    dchkpt.f dchkq3.f dchkql.f dchkqr.f dchkrq.f
-   dchksp.f dchksy.f dchksy_rook.f dchktb.f dchktp.f dchktr.f
+   dchksp.f dchksy.f dchksy_rook.f dchksy_aasen.f dchktb.f dchktp.f dchktr.f
    dchktz.f
    ddrvgt.f ddrvls.f ddrvpb.f
-   ddrvpp.f ddrvpt.f ddrvsp.f ddrvsy.f ddrvsy_rook.f
+   ddrvpp.f ddrvpt.f ddrvsp.f ddrvsy.f ddrvsy_rook.f ddrvsy_aasen.f
    derrgt.f derrlq.f derrls.f
    derrps.f derrql.f derrqp.f derrqr.f
    derrrq.f derrsy.f derrtr.f derrtz.f derrvx.f
@@ -106,7 +106,7 @@ set(DLINTST  dchkaa.f
    dqrt01.f dqrt01p.f dqrt02.f dqrt03.f dqrt11.f dqrt12.f
    dqrt13.f dqrt14.f dqrt15.f dqrt16.f dqrt17.f
    drqt01.f drqt02.f drqt03.f drzt01.f drzt02.f
-   dspt01.f dsyt01.f dsyt01_rook.f
+   dspt01.f dsyt01.f dsyt01_rook.f dsyt01_aasen.f
    dtbt02.f dtbt03.f dtbt05.f dtbt06.f dtpt01.f
    dtpt02.f dtpt03.f dtpt05.f dtpt06.f dtrt01.f
    dtrt02.f dtrt03.f dtrt05.f dtrt06.f
@@ -123,11 +123,11 @@ endif()
 
 set(ZLINTST  zchkaa.f
    zchkeq.f zchkgb.f zchkge.f zchkgt.f
-   zchkhe.f zchkhe_rook.f zchkhp.f zchklq.f zchkpb.f
+   zchkhe.f zchkhe_rook.f zchkhe_aasen.f zchkhp.f zchklq.f zchkpb.f
    zchkpo.f zchkps.f zchkpp.f zchkpt.f zchkq3.f zchkql.f
    zchkqr.f zchkrq.f zchksp.f zchksy.f zchksy_rook.f zchktb.f
    zchktp.f zchktr.f zchktz.f
-   zdrvgt.f zdrvhe.f zdrvhe_rook.f zdrvhp.f
+   zdrvgt.f zdrvhe.f zdrvhe_rook.f zdrvhe_aasen.f zdrvhp.f
    zdrvls.f zdrvpb.f zdrvpp.f zdrvpt.f
    zdrvsp.f zdrvsy.f zdrvsy_rook.f
    zerrgt.f zerrhe.f zerrlq.f
@@ -137,7 +137,7 @@ set(ZLINTST  zchkaa.f
    zgbt01.f zgbt02.f zgbt05.f zgelqs.f zgeqls.f zgeqrs.f
    zgerqs.f zget01.f zget02.f
    zget03.f zget04.f zget07.f zgtt01.f zgtt02.f
-   zgtt05.f zhet01.f zhet01.f zhet01_rook.f zhpt01.f zlaipd.f zlaptm.f zlarhs.f zlatb4.f zlatb5.f
+   zgtt05.f zhet01.f zhet01.f zhet01_rook.f zhet01_aasen.f zhpt01.f zlaipd.f zlaptm.f zlarhs.f zlatb4.f zlatb5.f
    zlatsp.f zlatsy.f zlattb.f zlattp.f zlattr.f
    zlavhe.f zlavhe_rook.f zlavhp.f zlavsp.f zlavsy.f zlavsy_rook.f zlqt01.f
    zlqt02.f zlqt03.f zpbt01.f zpbt02.f zpbt05.f
diff --git a/TESTING/LIN/Makefile b/TESTING/LIN/Makefile
index 57d4fe25..3232f0fd 100644
--- a/TESTING/LIN/Makefile
+++ b/TESTING/LIN/Makefile
@@ -51,10 +51,10 @@ SLINTST = schkaa.o \
    schkeq.o schkgb.o schkge.o schkgt.o \
    schklq.o schkpb.o schkpo.o schkps.o schkpp.o \
    schkpt.o schkq3.o schkql.o schkqr.o schkrq.o \
-   schksp.o schksy.o schksy_rook.o schktb.o schktp.o schktr.o \
+   schksp.o schksy.o schksy_rook.o schksy_aasen.o schktb.o schktp.o schktr.o \
    schktz.o \
    sdrvgt.o sdrvls.o sdrvpb.o \
-   sdrvpp.o sdrvpt.o sdrvsp.o  sdrvsy_rook.o\
+   sdrvpp.o sdrvpt.o sdrvsp.o  sdrvsy_rook.o sdrvsy_aasen.o\
    serrgt.o serrlq.o serrls.o \
    serrps.o serrql.o serrqp.o serrqr.o \
    serrrq.o serrtr.o serrtz.o \
@@ -70,7 +70,7 @@ SLINTST = schkaa.o \
    sqrt01.o sqrt01p.o sqrt02.o sqrt03.o sqrt11.o sqrt12.o \
    sqrt13.o sqrt14.o sqrt15.o sqrt16.o sqrt17.o \
    srqt01.o srqt02.o srqt03.o srzt01.o srzt02.o \
-   sspt01.o ssyt01.o ssyt01_rook.o \
+   sspt01.o ssyt01.o ssyt01_rook.o ssyt01_aasen.o\
    stbt02.o stbt03.o stbt05.o stbt06.o stpt01.o \
    stpt02.o stpt03.o stpt05.o stpt06.o strt01.o \
    strt02.o strt03.o strt05.o strt06.o \
@@ -86,11 +86,11 @@ endif
 
 CLINTST = cchkaa.o \
    cchkeq.o cchkgb.o cchkge.o cchkgt.o \
-   cchkhe.o cchkhe_rook.o cchkhp.o cchklq.o cchkpb.o \
+   cchkhe.o cchkhe_rook.o cchkhe_aasen.o cchkhp.o cchklq.o cchkpb.o \
    cchkpo.o cchkps.o cchkpp.o cchkpt.o cchkq3.o cchkql.o \
    cchkqr.o cchkrq.o cchksp.o cchksy.o cchksy_rook.o cchktb.o \
    cchktp.o cchktr.o cchktz.o \
-   cdrvgt.o cdrvhe_rook.o cdrvhp.o \
+   cdrvgt.o cdrvhe_rook.o cdrvhe_aasen.o cdrvhp.o \
    cdrvls.o cdrvpb.o cdrvpp.o cdrvpt.o \
    cdrvsp.o cdrvsy_rook.o \
    cerrgt.o cerrlq.o \
@@ -99,7 +99,7 @@ CLINTST = cchkaa.o \
    cgbt01.o cgbt02.o cgbt05.o cgelqs.o cgeqls.o cgeqrs.o \
    cgerqs.o cget01.o cget02.o \
    cget03.o cget04.o cget07.o cgtt01.o cgtt02.o \
-   cgtt05.o chet01.o chet01_rook.o chpt01.o claipd.o claptm.o clarhs.o clatb4.o clatb5.o \
+   cgtt05.o chet01.o chet01_rook.o chet01_aasen.o chpt01.o claipd.o claptm.o clarhs.o clatb4.o clatb5.o \
    clatsp.o clatsy.o clattb.o clattp.o clattr.o \
    clavhe.o clavhe_rook.o clavhp.o clavsp.o clavsy.o clavsy_rook.o clqt01.o \
    clqt02.o clqt03.o cpbt01.o cpbt02.o cpbt05.o \
@@ -129,10 +129,10 @@ DLINTST = dchkaa.o \
    dchkeq.o dchkgb.o dchkge.o dchkgt.o \
    dchklq.o dchkpb.o dchkpo.o dchkps.o dchkpp.o \
    dchkpt.o dchkq3.o dchkql.o dchkqr.o dchkrq.o \
-   dchksp.o dchksy.o dchksy_rook.o dchktb.o dchktp.o dchktr.o \
+   dchksp.o dchksy.o dchksy_rook.o dchksy_aasen.o dchktb.o dchktp.o dchktr.o \
    dchktz.o \
    ddrvgt.o ddrvls.o ddrvpb.o \
-   ddrvpp.o ddrvpt.o ddrvsp.o ddrvsy_rook.o \
+   ddrvpp.o ddrvpt.o ddrvsp.o ddrvsy_rook.o ddrvsy_aasen.o\
    derrgt.o derrlq.o derrls.o \
    derrps.o derrql.o derrqp.o derrqr.o \
    derrrq.o derrtr.o derrtz.o \
@@ -148,7 +148,7 @@ DLINTST = dchkaa.o \
    dqrt01.o dqrt01p.o  dqrt02.o dqrt03.o dqrt11.o dqrt12.o \
    dqrt13.o dqrt14.o dqrt15.o dqrt16.o dqrt17.o \
    drqt01.o drqt02.o drqt03.o drzt01.o drzt02.o \
-   dspt01.o dsyt01.o dsyt01_rook.o \
+   dspt01.o dsyt01.o dsyt01_rook.o dsyt01_aasen.o\
    dtbt02.o dtbt03.o dtbt05.o dtbt06.o dtpt01.o \
    dtpt02.o dtpt03.o dtpt05.o dtpt06.o dtrt01.o \
    dtrt02.o dtrt03.o dtrt05.o dtrt06.o \
@@ -165,11 +165,11 @@ endif
 
 ZLINTST = zchkaa.o \
    zchkeq.o zchkgb.o zchkge.o zchkgt.o \
-   zchkhe.o zchkhe_rook.o zchkhp.o zchklq.o zchkpb.o \
+   zchkhe.o zchkhe_rook.o zchkhe_aasen.o zchkhp.o zchklq.o zchkpb.o \
    zchkpo.o zchkps.o zchkpp.o zchkpt.o zchkq3.o zchkql.o \
    zchkqr.o zchkrq.o zchksp.o zchksy.o zchksy_rook.o zchktb.o \
    zchktp.o zchktr.o zchktz.o \
-   zdrvgt.o zdrvhe_rook.o zdrvhp.o \
+   zdrvgt.o zdrvhe_rook.o zdrvhe_aasen.o zdrvhp.o \
    zdrvls.o zdrvpb.o zdrvpp.o zdrvpt.o \
    zdrvsp.o zdrvsy_rook.o \
    zerrgt.o zerrlq.o \
@@ -178,7 +178,7 @@ ZLINTST = zchkaa.o \
    zgbt01.o zgbt02.o zgbt05.o zgelqs.o zgeqls.o zgeqrs.o \
    zgerqs.o zget01.o zget02.o \
    zget03.o zget04.o zget07.o zgtt01.o zgtt02.o \
-   zgtt05.o zhet01.o zhet01_rook.o zhpt01.o zlaipd.o zlaptm.o zlarhs.o zlatb4.o zlatb5.o \
+   zgtt05.o zhet01.o zhet01_rook.o zhet01_aasen.o zhpt01.o zlaipd.o zlaptm.o zlarhs.o zlatb4.o zlatb5.o \
    zlatsp.o zlatsy.o zlattb.o zlattp.o zlattr.o \
    zlavhe.o zlavhe_rook.o zlavhp.o zlavsp.o zlavsy.o zlavsy_rook.o zlqt01.o \
    zlqt02.o zlqt03.o zpbt01.o zpbt02.o zpbt05.o \
diff --git a/TESTING/LIN/aladhd.f b/TESTING/LIN/aladhd.f
index b2bb2e03..3a53e0bd 100644
--- a/TESTING/LIN/aladhd.f
+++ b/TESTING/LIN/aladhd.f
@@ -53,6 +53,8 @@
 *>                     with "rook" (bounded Bunch-Kaufman) pivoting
 *>             _SP:  Symmetric indefinite packed,
 *>                     with partial (Bunch-Kaufman) pivoting
+*>             _HA:  (complex) Hermitian ,
+*>                     Assen Algorithm
 *>             _HE:  (complex) Hermitian indefinite,
 *>                     with partial (Bunch-Kaufman) pivoting
 *>             _HR:  (complex) Hermitian indefinite,
@@ -275,7 +277,27 @@
          WRITE( IOUNIT, FMT = 9979 )3
          WRITE( IOUNIT, FMT = '( '' Messages:'' )' )
 *
-      ELSE IF( LSAMEN( 2, P2, 'HE' ) .OR. LSAMEN( 2, P2, 'HP' ) ) THEN
+      ELSE IF( LSAMEN( 2, P2, 'HA' ) ) THEN
+*
+*        HA: Hermitian 
+*            Aasen algorithm
+         WRITE( IOUNIT, FMT = 9971 )PATH, 'Hermitian'
+*
+         WRITE( IOUNIT, FMT = '( '' Matrix types:'' )' )
+         WRITE( IOUNIT, FMT = 9983 )
+*
+         WRITE( IOUNIT, FMT = '( '' Test ratios:'' )' )
+         WRITE( IOUNIT, FMT = 9974 )1
+         WRITE( IOUNIT, FMT = 9980 )2
+         WRITE( IOUNIT, FMT = 9979 )3
+         WRITE( IOUNIT, FMT = 9977 )4
+         WRITE( IOUNIT, FMT = 9978 )5
+         WRITE( IOUNIT, FMT = 9976 )6
+         WRITE( IOUNIT, FMT = '( '' Messages:'' )' )
+      
+      
+      ELSE IF( LSAMEN( 2, P2, 'HE' ) .OR. 
+     $         LSAMEN( 2, P2, 'HP' ) ) THEN
 *
 *        HE: Hermitian indefinite full
 *            with partial (Bunch-Kaufman) pivoting algorithm
@@ -336,6 +358,8 @@
      $      ' positive definite band matrices' )
  9993 FORMAT( / 1X, A3, ' drivers:  ', A9,
      $      ' positive definite tridiagonal' )
+ 9971 FORMAT( / 1X, A3, ' drivers:  ', A9, ' indefinite matrices',
+     $     ', "Aasen" Algorithm' )
  9992 FORMAT( / 1X, A3, ' drivers:  ', A9, ' indefinite matrices',
      $     ', "rook" (bounded Bunch-Kaufman) pivoting' )
  9991 FORMAT( / 1X, A3, ' drivers:  ', A9,
diff --git a/TESTING/LIN/alaerh.f b/TESTING/LIN/alaerh.f
index 2f45b1d2..4fec4522 100644
--- a/TESTING/LIN/alaerh.f
+++ b/TESTING/LIN/alaerh.f
@@ -490,6 +490,7 @@
       ELSE IF( LSAMEN( 2, P2, 'SY' )
      $         .OR. LSAMEN( 2, P2, 'SR' )
      $         .OR. LSAMEN( 2, P2, 'HE' )
+     $         .OR. LSAMEN( 2, P2, 'HA' )
      $         .OR. LSAMEN( 2, P2, 'HR' ) ) THEN
 *
 *        xSY: symmetric indefinite matrices
@@ -498,6 +499,8 @@
 *             with rook (bounded Bunch-Kaufman) pivoting;
 *        xHE: Hermitian indefinite matrices
 *             with partial (Bunch-Kaufman) pivoting.
+*        xHA: Hermitian matrices
+*             Aasen Algorithm
 *        xHR: Hermitian indefinite matrices
 *             with rook (bounded Bunch-Kaufman) pivoting;
 *
diff --git a/TESTING/LIN/alahd.f b/TESTING/LIN/alahd.f
index 8d56c798..995a5444 100644
--- a/TESTING/LIN/alahd.f
+++ b/TESTING/LIN/alahd.f
@@ -53,6 +53,8 @@
 *>                     with "rook" (bounded Bunch-Kaufman) pivoting
 *>             _SP:  Symmetric indefinite packed,
 *>                     with partial (Bunch-Kaufman) pivoting
+*>             _HA:  (complex) Hermitian ,
+*>                     with Aasen Algorithm
 *>             _HE:  (complex) Hermitian indefinite,
 *>                     with partial (Bunch-Kaufman) pivoting
 *>             _HR:  Symmetric indefinite,
@@ -355,6 +357,28 @@
          WRITE( IOUNIT, FMT = 9955 )8
          WRITE( IOUNIT, FMT = '( '' Messages:'' )' )
 *
+      ELSE IF( LSAMEN( 2, P2, 'HA' )  ) THEN
+*
+*        HA: Hermitian,
+*            with Assen Algorithm
+*
+         WRITE( IOUNIT, FMT = 9992 )PATH, 'Hermitian'
+*
+         WRITE( IOUNIT, FMT = '( '' Matrix types:'' )' )
+         WRITE( IOUNIT, FMT = 9972 )
+*
+         WRITE( IOUNIT, FMT = '( '' Test ratios:'' )' )
+         WRITE( IOUNIT, FMT = 9953 )1
+         WRITE( IOUNIT, FMT = 9961 )2
+         WRITE( IOUNIT, FMT = 9960 )3
+         WRITE( IOUNIT, FMT = 9960 )4
+         WRITE( IOUNIT, FMT = 9959 )5
+         WRITE( IOUNIT, FMT = 9958 )6
+         WRITE( IOUNIT, FMT = 9956 )7
+         WRITE( IOUNIT, FMT = 9957 )8
+         WRITE( IOUNIT, FMT = 9955 )9
+         WRITE( IOUNIT, FMT = '( '' Messages:'' )' )
+*
       ELSE IF( LSAMEN( 2, P2, 'HE' )  ) THEN
 *
 *        HE: Hermitian indefinite full,
diff --git a/TESTING/LIN/cchkaa.f b/TESTING/LIN/cchkaa.f
index 8a97fbb5..1f8b2c65 100644
--- a/TESTING/LIN/cchkaa.f
+++ b/TESTING/LIN/cchkaa.f
@@ -51,6 +51,7 @@
 *> CPT   12               List types on next line if 0 < NTYPES < 12
 *> CHE   10               List types on next line if 0 < NTYPES < 10
 *> CHR   10               List types on next line if 0 < NTYPES < 10
+*> CHA   10               List types on next line if 0 < NTYPES < 10
 *> CHP   10               List types on next line if 0 < NTYPES < 10
 *> CSY   11               List types on next line if 0 < NTYPES < 11
 *> CSR   11               List types on next line if 0 < NTYPES < 11
@@ -102,17 +103,17 @@
 *> \author Univ. of Colorado Denver
 *> \author NAG Ltd.
 *
-*> \date November 2015
+*> \date November 2016
 *
 *> \ingroup complex_lin
 *
 *  =====================================================================
       PROGRAM CCHKAA
 *
-*  -- LAPACK test routine (version 3.6.0) --
+*  -- LAPACK test routine (version 3.7.0) --
 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
-*     November 2015
+*     November 2016
 *
 *  =====================================================================
 *
@@ -641,6 +642,33 @@
             WRITE( NOUT, FMT = 9988 )PATH
          END IF
 *
+      ELSE IF( LSAMEN( 2, C2, 'HA' ) ) THEN
+*
+*        HA:  Hermitian matrices,
+*             Aasen Algorithm
+*
+         NTYPES = 10
+         CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT )
+*
+         IF( TSTCHK ) THEN
+            CALL CCHKHE_AASEN( DOTYPE, NN, NVAL, NNB2, NBVAL2, NNS,
+     $                         NSVAL, THRESH, TSTERR, LDA, 
+     $                         A( 1, 1 ), A( 1, 2 ), A( 1, 3 ),
+     $                         B( 1, 1 ), B( 1, 2 ), B( 1, 3 ),
+     $                         WORK, RWORK, IWORK, NOUT )
+         ELSE
+            WRITE( NOUT, FMT = 9989 )PATH
+         END IF
+*
+         IF( TSTDRV ) THEN
+            CALL CDRVHE_AASEN( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, 
+     $                         LDA, A( 1, 1 ), A( 1, 2 ), A( 1, 3 ), 
+     $                              B( 1, 1 ), B( 1, 2 ), B( 1, 3 ), 
+     $                         WORK, RWORK, IWORK, NOUT )
+         ELSE
+            WRITE( NOUT, FMT = 9988 )PATH
+         END IF
+*
       ELSE IF( LSAMEN( 2, C2, 'HR' ) ) THEN
 *
 *        HR:  Hermitian indefinite matrices,
@@ -867,7 +895,6 @@
          ELSE
             WRITE( NOUT, FMT = 9989 )PATH
          END IF
-
 *
       ELSE IF( LSAMEN( 2, C2, 'RQ' ) ) THEN
 *
@@ -926,7 +953,6 @@
          ELSE
             WRITE( NOUT, FMT = 9989 )PATH
          END IF
-
 *
       ELSE IF( LSAMEN( 2, C2, 'LS' ) ) THEN
 *
diff --git a/TESTING/LIN/cchkhe_aasen.f b/TESTING/LIN/cchkhe_aasen.f
new file mode 100644
index 00000000..355dd25b
--- /dev/null
+++ b/TESTING/LIN/cchkhe_aasen.f
@@ -0,0 +1,578 @@
+*> \brief \b CCHKHE_AASEN
+*
+*  =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+*            http://www.netlib.org/lapack/explore-html/
+*
+*  Definition:
+*  ===========
+*
+*       SUBROUTINE CCHKHE_AASEN( 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    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_AASEN tests CHETRF_AASEN, -TRS_AASEN.
+*> \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 (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 November 2016
+*
+*
+*> \ingroup complex_lin
+*
+*  =====================================================================
+      SUBROUTINE CCHKHE_AASEN( 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.7.0) --
+*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*     November 2016
+*
+      IMPLICIT NONE
+*
+*     .. 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
+      PARAMETER    ( ZERO = 0.0E+0 )
+      COMPLEX      CZERO
+      PARAMETER    ( CZERO = ( 0.0E+0, 0.0E+0 ) )
+      INTEGER      NTYPES
+      PARAMETER    ( NTYPES = 10 )
+      INTEGER      NTESTS
+      PARAMETER    ( NTESTS = 9 )
+*     ..
+*     .. Local Scalars ..
+      LOGICAL      TRFCON, ZEROT
+      CHARACTER    DIST, TYPE, UPLO, XTYPE
+      CHARACTER*3  PATH, MATPATH
+      INTEGER      I, I1, I2, IMAT, IN, INB, INFO, IOFF, IRHS,
+     $             IUPLO, IZERO, J, K, KL, KU, LDA, LWORK, MODE,
+     $             N, NB, NERRS, NFAIL, NIMAT, NRHS, NRUN, NT
+      REAL         ANORM, CNDNUM, RCOND, RCONDC
+*     ..
+*     .. Local Arrays ..
+      CHARACTER    UPLOS( 2 )
+      INTEGER      ISEED( 4 ), ISEEDY( 4 )
+      REAL         RESULT( NTESTS )
+*     ..
+*     .. External Functions ..
+      REAL         DGET06, CLANHE
+      EXTERNAL     DGET06, CLANHE
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL     ALAERH, ALAHD, ALASUM, XLAENV, CERRHE, CGET04,
+     $             ZHECON, CHERFS, CHET01, CHETRF_AASEN, ZHETRI2,
+     $             CHETRS_AASEN, CLACPY, CLAIPD, CLARHS, CLATB4, 
+     $             CLATMS, CPOT02, ZPOT03, ZPOT05
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC    REAL, IMAG, MAX, MIN
+*     ..
+*     .. 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.
+*
+*
+*     Test path
+*
+      PATH( 1: 1 ) = 'Complex precision'
+      PATH( 2: 3 ) = 'HA'
+*
+*     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 180 IN = 1, NN
+         N = NVAL( IN )
+         IF( N .GT. NMAX ) THEN
+            NFAIL = NFAIL + 1
+            WRITE(NOUT, 9995) 'M ', N, NMAX
+            GO TO 180
+         END IF
+         LDA = MAX( N, 1 )
+         XTYPE = 'N'
+         NIMAT = NTYPES
+         IF( N.LE.0 )
+     $      NIMAT = 1
+*
+         IZERO = 0
+         DO 170 IMAT = 1, NIMAT
+*
+*           Do the tests only if DOTYPE( IMAT ) is true.
+*
+            IF( .NOT.DOTYPE( IMAT ) )
+     $         GO TO 170
+*
+*           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 170
+*
+*           Do first for UPLO = 'U', then for UPLO = 'L'
+*
+            DO 160 IUPLO = 1, 2
+               UPLO = UPLOS( IUPLO )
+*
+*              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 160
+               END IF
+*
+*              For 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
+                     IF( IUPLO.EQ.1 ) THEN
+*
+*                       Set the first IZERO rows and columns to zero.
+*
+                        IOFF = 0
+                        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
+                        IZERO = 1
+                     ELSE
+*
+*                       Set the last IZERO rows and columns to zero.
+*
+                        IOFF = 0
+                        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 test matrix A.
+*
+*
+*              Set the imaginary part of the diagonals.
+*
+               CALL CLAIPD( N, A, LDA+1, 0 )
+*
+*              Do for each value of NB in NBVAL
+*
+               DO 150 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 = ( NB+1 )*LDA
+                  SRNAMT = 'CHETRF_AASEN'
+                  CALL CHETRF_AASEN( UPLO, N, AFAC, LDA, IWORK, AINV, 
+     $                               LWORK, INFO )
+*
+*                 Adjust the expected value of INFO to account for
+*                 pivoting.
+*
+                  IF( IZERO.GT.0 ) THEN
+                     J = 1
+                     K = IZERO
+  100                CONTINUE
+                     IF( J.EQ.K ) THEN
+                        K = IWORK( J )
+                     ELSE IF( IWORK( J ).EQ.K ) THEN
+                        K = J
+                     END IF
+                     IF( J.LT.K ) THEN
+                        J = J + 1
+                        GO TO 100
+                     END IF
+                  ELSE
+                     K = 0
+                  END IF
+*
+*                 Check error code from CHETRF and handle error.
+*
+                  IF( INFO.NE.K ) THEN
+                     CALL ALAERH( PATH, 'CHETRF_AASEN', INFO, K, UPLO, 
+     $                            N, N, -1, -1, NB, IMAT, NFAIL, NERRS, 
+     $                            NOUT )
+                  END IF
+*
+*                 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_AASEN( UPLO, N, A, LDA, AFAC, LDA, IWORK,
+     $                               AINV, LDA, RWORK, RESULT( 1 ) )
+                  NT = 1
+*
+*
+*                 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
+*
+*                 Do only the condition estimate if INFO is not 0.
+*
+                  IF( TRFCON ) THEN
+                     RCONDC = ZERO
+                     GO TO 140
+                  END IF
+*
+*                 Do for each value of NRHS in NSVAL.
+*
+                  DO 130 IRHS = 1, NNS
+                     NRHS = NSVAL( IRHS )
+*
+*+    TEST 3 (Using TRS)
+*                 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_AASEN'
+                     LWORK = 3*N-2
+                     CALL CHETRS_AASEN( UPLO, N, NRHS, AFAC, LDA, IWORK,
+     $                                  X, LDA, WORK, LWORK, INFO )
+*
+*                    Check error code from CHETRS and handle error.
+*
+                     IF( INFO.NE.0 ) THEN
+                        CALL ALAERH( PATH, 'CHETRS_AASEN', INFO, 0,
+     $                                UPLO, N, N, -1, -1, NRHS, IMAT,
+     $                                NFAIL, NERRS, NOUT )
+                     END IF
+*
+                     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( 2 ) )
+*
+*                    Print information about the tests that did not pass
+*                    the threshold.
+*
+                     DO 120 K = 2, 2
+                        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
+  120                CONTINUE
+                     NRUN = NRUN + 1
+*
+*                 End do for each value of NRHS in NSVAL.
+*
+  130             CONTINUE
+  140             CONTINUE
+  150          CONTINUE
+  160       CONTINUE
+  170    CONTINUE
+  180 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 )
+ 9995 FORMAT( ' Invalid input value: ', A4, '=', I6, '; must be <=',
+     $      I6 )
+      RETURN
+*
+*     End of CCHKHE_AASEN
+*
+      END
diff --git a/TESTING/LIN/cdrvhe_aasen.f b/TESTING/LIN/cdrvhe_aasen.f
new file mode 100644
index 00000000..617f307a
--- /dev/null
+++ b/TESTING/LIN/cdrvhe_aasen.f
@@ -0,0 +1,529 @@
+*> \brief \b CDRVHE_AASEN
+*
+*  =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+*            http://www.netlib.org/lapack/explore-html/
+*
+*  Definition:
+*  ===========
+*
+*       SUBROUTINE CDRVHE_AASEN( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, NMAX,
+*                                A, AFAC, AINV, B, X, XACT, WORK, RWORK, IWORK,
+*                                NOUT )
+*
+*       .. Scalar Arguments ..
+*       LOGICAL            TSTERR
+*       INTEGER            NMAX, NN, NOUT, NRHS
+*       REAL               THRESH
+*       ..
+*       .. Array Arguments ..
+*       LOGICAL            DOTYPE( * )
+*       INTEGER            IWORK( * ), NVAL( * )
+*       REAL               RWORK( * )
+*       COMPLEX            A( * ), AFAC( * ), AINV( * ), B( * ),
+*      $                   WORK( * ), X( * ), XACT( * )
+*       ..
+*
+*
+*> \par Purpose:
+*  =============
+*>
+*> \verbatim
+*>
+*> CDRVHE_AASEN tests the driver routine CHESV_AASEN.
+*> \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] NRHS
+*> \verbatim
+*>          NRHS is INTEGER
+*>          The number of right hand side vectors to be generated for
+*>          each linear system.
+*> \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*NRHS)
+*> \endverbatim
+*>
+*> \param[out] X
+*> \verbatim
+*>          X is COMPLEX array, dimension (NMAX*NRHS)
+*> \endverbatim
+*>
+*> \param[out] XACT
+*> \verbatim
+*>          XACT is COMPLEX array, dimension (NMAX*NRHS)
+*> \endverbatim
+*>
+*> \param[out] WORK
+*> \verbatim
+*>          WORK is COMPLEX array, dimension (NMAX*max(2,NRHS))
+*> \endverbatim
+*>
+*> \param[out] RWORK
+*> \verbatim
+*>          RWORK is REAL array, dimension (NMAX+2*NRHS)
+*> \endverbatim
+*>
+*> \param[out] IWORK
+*> \verbatim
+*>          IWORK is INTEGER array, dimension (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 November 2016
+*
+*> \ingroup complex_lin
+*
+*  =====================================================================
+      SUBROUTINE CDRVHE_AASEN( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR,
+     $                         NMAX, A, AFAC, AINV, B, X, XACT, WORK,
+     $                         RWORK, IWORK, NOUT )
+*
+*  -- LAPACK test routine (version 3.7.0) --
+*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*     November 2016
+*
+*     .. Scalar Arguments ..
+      LOGICAL            TSTERR
+      INTEGER            NMAX, NN, NOUT, NRHS
+      REAL               THRESH
+*     ..
+*     .. Array Arguments ..
+      LOGICAL            DOTYPE( * )
+      INTEGER            IWORK( * ), NVAL( * )
+      REAL               RWORK( * )
+      COMPLEX            A( * ), AFAC( * ), AINV( * ), B( * ),
+     $                   WORK( * ), X( * ), XACT( * )
+*     ..
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      REAL               ONE, ZERO
+      PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
+      INTEGER            NTYPES, NTESTS
+      PARAMETER          ( NTYPES = 10, NTESTS = 3 )
+      INTEGER            NFACT
+      PARAMETER          ( NFACT = 2 )
+*     ..
+*     .. Local Scalars ..
+      LOGICAL            ZEROT
+      CHARACTER          DIST, FACT, TYPE, UPLO, XTYPE
+      CHARACTER*3        MATPATH, PATH
+      INTEGER            I, I1, I2, IFACT, IMAT, IN, INFO, IOFF, IUPLO,
+     $                   IZERO, J, K, K1, KL, KU, LDA, LWORK, MODE, N,
+     $                   NB, NBMIN, NERRS, NFAIL, NIMAT, NRUN, NT
+      REAL               AINVNM, ANORM, CNDNUM, RCOND, RCONDC
+*     ..
+*     .. Local Arrays ..
+      CHARACTER          FACTS( NFACT ), UPLOS( 2 )
+      INTEGER            ISEED( 4 ), ISEEDY( 4 )
+      REAL               RESULT( NTESTS )
+*     ..
+*     .. External Functions ..
+      REAL               CLANHE, SGET06
+      EXTERNAL           CLANHE, SGET06
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           ALADHD, ALAERH, ALASVM, XLAENV, CERRVX,
+     $                   CGET04, CLACPY, CLARHS, CLATB4, CLATMS,
+     $                   CHESV_AASEN, CHET01_AASEN, CPOT02,
+     $                   CHETRF_AASEN
+*     ..
+*     .. Scalars in Common ..
+      LOGICAL            LERR, OK
+      CHARACTER*32       SRNAMT
+      INTEGER            INFOT, NUNIT
+*     ..
+*     .. Common blocks ..
+      COMMON             / INFOC / INFOT, NUNIT, OK, LERR
+      COMMON             / SRNAMC / SRNAMT
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          CMPLX, MAX, MIN
+*     ..
+*     .. Data statements ..
+      DATA               ISEEDY / 1988, 1989, 1990, 1991 /
+      DATA               UPLOS / 'U', 'L' / , FACTS / 'F', 'N' /
+*     ..
+*     .. Executable Statements ..
+*
+*     Initialize constants and the random number seed.
+*
+*     Test path
+*
+      PATH( 1: 1 ) = 'Complex precision'
+      PATH( 2: 3 ) = 'HA'
+*
+*     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
+      LWORK = MAX( 2*NMAX, NMAX*NRHS )
+*
+*     Test the error exits
+*
+      IF( TSTERR )
+     $   CALL CERRVX( PATH, NOUT )
+      INFOT = 0
+*
+*     Set the block size and minimum block size for testing.
+*
+      NB = 1
+      NBMIN = 2
+      CALL XLAENV( 1, NB )
+      CALL XLAENV( 2, NBMIN )
+*
+*     Do for each value of N in NVAL
+*
+      DO 180 IN = 1, NN
+         N = NVAL( IN )
+         LDA = MAX( N, 1 )
+         XTYPE = 'N'
+         NIMAT = NTYPES
+         IF( N.LE.0 )
+     $      NIMAT = 1
+*
+         DO 170 IMAT = 1, NIMAT
+*
+*           Do the tests only if DOTYPE( IMAT ) is true.
+*
+            IF( .NOT.DOTYPE( IMAT ) )
+     $         GO TO 170
+*
+*           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 170
+*
+*           Do first for UPLO = 'U', then for UPLO = 'L'
+*
+            DO 160 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 )
+                     GO TO 160
+                  END IF
+*
+*                 For 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 ) = ZERO
+   20                      CONTINUE
+                           IOFF = IOFF + IZERO
+                           DO 30 I = IZERO, N
+                              A( IOFF ) = ZERO
+                              IOFF = IOFF + LDA
+   30                      CONTINUE
+                        ELSE
+                           IOFF = IZERO
+                           DO 40 I = 1, IZERO - 1
+                              A( IOFF ) = ZERO
+                              IOFF = IOFF + LDA
+   40                      CONTINUE
+                           IOFF = IOFF - IZERO
+                           DO 50 I = IZERO, N
+                              A( IOFF+I ) = ZERO
+   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 ) = ZERO
+   60                         CONTINUE
+                              IOFF = IOFF + LDA
+   70                      CONTINUE
+                           IZERO = 1
+                        ELSE
+*
+*                       Set the first IZERO rows and columns to zero.
+*
+                           IOFF = 0
+                           DO 90 J = 1, N
+                              I1 = MAX( J, IZERO )
+                              DO 80 I = I1, N
+                                 A( IOFF+I ) = ZERO
+   80                         CONTINUE
+                              IOFF = IOFF + LDA
+   90                      CONTINUE
+                        END IF
+                     END IF
+                  ELSE
+                     IZERO = 0
+                  END IF
+*
+*                 End generate the test matrix A.
+*
+*
+               DO 150 IFACT = 1, NFACT
+*
+*                 Do first for FACT = 'F', then for other values.
+*
+                  FACT = FACTS( IFACT )
+*
+*                 Compute the condition number for comparison with
+*                 the value returned by CHESVX.
+*
+                  IF( ZEROT ) THEN
+                     IF( IFACT.EQ.1 )
+     $                  GO TO 150
+                     RCONDC = ZERO
+*
+                  ELSE IF( IFACT.EQ.1 ) THEN
+*
+*                    Compute the 1-norm of A.
+*
+                     ANORM = CLANHE( '1', UPLO, N, A, LDA, RWORK )
+*
+*                    Factor the matrix A.
+*
+c                     CALL CLACPY( UPLO, N, N, A, LDA, AFAC, LDA )
+c                     SRNAMT = 'CHETRF_AASEN'
+c                     CALL CHETRF_AASEN( UPLO, N, AFAC, LDA, IWORK,
+c     $                            WORK, LWORK, INFO )
+*
+*                    Compute inv(A) and take its norm.
+*
+c                     CALL CLACPY( UPLO, N, N, AFAC, LDA, AINV, LDA )
+c                     LWORK = (N+NB+1)*(NB+3)
+c                     SRNAMT = 'CHETRI2'
+c                     CALL CHETRI2( UPLO, N, AINV, LDA, IWORK, WORK,
+c     $                            LWORK, INFO )
+c                     AINVNM = CLANHE( '1', UPLO, N, AINV, LDA, RWORK )
+*
+*                    Compute the 1-norm condition number of A.
+*
+c                     IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN
+c                        RCONDC = ONE
+c                     ELSE
+c                        RCONDC = ( ONE / ANORM ) / AINVNM
+c                     END IF
+                  END IF
+*
+*                 Form an exact solution and set the right hand side.
+*
+                  SRNAMT = 'CLARHS'
+                  CALL CLARHS( MATPATH, XTYPE, UPLO, ' ', N, N, KL, KU,
+     $                         NRHS, A, LDA, XACT, LDA, B, LDA, ISEED,
+     $                         INFO )
+                  XTYPE = 'C'
+*
+*                 --- Test CHESV_AASEN  ---
+*
+                  IF( IFACT.EQ.2 ) THEN
+                     CALL CLACPY( UPLO, N, N, A, LDA, AFAC, LDA )
+                     CALL CLACPY( 'Full', N, NRHS, B, LDA, X, LDA )
+*
+*                    Factor the matrix and solve the system using CHESV_AASEN.
+*
+                     SRNAMT = 'CHESV_AASEN '
+                     CALL CHESV_AASEN( UPLO, N, NRHS, AFAC, LDA, IWORK,
+     $                                 X, LDA, WORK, LWORK, INFO )
+*
+*                    Adjust the expected value of INFO to account for
+*                    pivoting.
+*
+                     IF( IZERO.GT.0 ) THEN
+                        J = 1
+                        K = IZERO
+  100                   CONTINUE
+                        IF( J.EQ.K ) THEN
+                           K = IWORK( J )
+                        ELSE IF( IWORK( J ).EQ.K ) THEN
+                           K = J
+                        END IF
+                        IF( J.LT.K ) THEN
+                           J = J + 1
+                           GO TO 100
+                        END IF
+                     ELSE
+                        K = 0
+                     END IF
+*
+*                    Check error code from CHESV_AASEN .
+*					 
+                     IF( INFO.NE.K ) THEN
+                        CALL ALAERH( PATH, 'CHESV_AASEN', INFO, K, 
+     $                               UPLO, N, N, -1, -1, NRHS, 
+     $                               IMAT, NFAIL, NERRS, NOUT )
+                        GO TO 120
+                     ELSE IF( INFO.NE.0 ) THEN
+                        GO TO 120
+                     END IF
+*
+*                    Reconstruct matrix from factors and compute
+*                    residual.
+*
+                     CALL CHET01_AASEN( UPLO, N, A, LDA, AFAC, LDA,
+     $                                  IWORK, AINV, LDA, RWORK, 
+     $                                  RESULT( 1 ) )
+*
+*                    Compute residual of the computed solution.
+*
+                     CALL CLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA )
+                     CALL CPOT02( UPLO, N, NRHS, A, LDA, X, LDA, WORK,
+     $                            LDA, RWORK, RESULT( 2 ) )
+*
+*                    Check solution from generated exact solution.
+*
+                     CALL CGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
+     $                            RESULT( 3 ) )
+                     NT = 3
+*
+*                    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 ALADHD( NOUT, PATH )
+                           WRITE( NOUT, FMT = 9999 )'CHESV_AASEN ',
+     $                         UPLO, N, IMAT, K, RESULT( K )
+                           NFAIL = NFAIL + 1
+                        END IF
+  110                CONTINUE
+                     NRUN = NRUN + NT
+  120                CONTINUE
+                  END IF
+*
+  150          CONTINUE
+*
+  160       CONTINUE
+  170    CONTINUE
+  180 CONTINUE
+*
+*     Print a summary of the results.
+*
+      CALL ALASVM( PATH, NOUT, NFAIL, NRUN, NERRS )
+*
+ 9999 FORMAT( 1X, A, ', UPLO=''', A1, ''', N =', I5, ', type ', I2,
+     $      ', test ', I2, ', ratio =', G12.5 )
+      RETURN
+*
+*     End of CDRVHE_AASEN
+*
+      END
diff --git a/TESTING/LIN/cerrvx.f b/TESTING/LIN/cerrvx.f
index 52ca890d..380dc699 100644
--- a/TESTING/LIN/cerrvx.f
+++ b/TESTING/LIN/cerrvx.f
@@ -93,7 +93,7 @@
      $                   CHESV, CHESV_ROOK, CHESVX, CHKXER, CHPSV,
      $                   CHPSVX, CPBSV, CPBSVX, CPOSV, CPOSVX, CPPSV,
      $                   CPPSVX, CPTSV, CPTSVX, CSPSV, CSPSVX, CSYSV,
-     $                   CSYSV_ROOK, CSYSVX
+     $                   CSYSV_AASEN, CSYSV_ROOK, CSYSVX
 *     ..
 *     .. Scalars in Common ..
       LOGICAL            LERR, OK
@@ -632,6 +632,25 @@
      $                RCOND, R1, R2, W, 3, RW, INFO )
          CALL CHKXER( 'CHESVX', INFOT, NOUT, LERR, OK )
 *
+      ELSE IF( LSAMEN( 2, C2, 'HA' ) ) THEN
+*
+*        CHESV_AASEN
+*
+        SRNAMT = 'CHESV_AASEN'
+        INFOT = 1
+        CALL CHESV_AASEN( '/', 0, 0, A, 1, IP, B, 1, W, 1, INFO )
+        CALL CHKXER( 'CHESV_AASEN', INFOT, NOUT, LERR, OK )
+        INFOT = 2
+        CALL CHESV_AASEN( 'U', -1, 0, A, 1, IP, B, 1, W, 1, INFO )
+        CALL CHKXER( 'CHESV_AASEN', INFOT, NOUT, LERR, OK )
+        INFOT = 3
+        CALL CHESV_AASEN( 'U', 0, -1, A, 1, IP, B, 1, W, 1, INFO )
+        CALL CHKXER( 'CHESV_AASEN', INFOT, NOUT, LERR, OK )
+        INFOT = 8
+        CALL CHESV_AASEN( 'U', 2, 0, A, 2, IP, B, 1, W, 1, INFO )
+        CALL CHKXER( 'CHESV_AASEN', INFOT, NOUT, LERR, OK )
+*
+
       ELSE IF( LSAMEN( 2, C2, 'HR' ) ) THEN
 *
 *        CHESV_ROOK
diff --git a/TESTING/LIN/chet01_aasen.f b/TESTING/LIN/chet01_aasen.f
new file mode 100644
index 00000000..d87a6102
--- /dev/null
+++ b/TESTING/LIN/chet01_aasen.f
@@ -0,0 +1,267 @@
+*> \brief \b CHET01_AASEN
+*
+*  =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+*            http://www.netlib.org/lapack/explore-html/
+*
+*  Definition:
+*  ===========
+*
+*       SUBROUTINE CHET01_AASEN( UPLO, N, A, LDA, AFAC, LDAFAC, IPIV, 
+*                          C, LDC, RWORK, RESID )
+*
+*       .. Scalar Arguments ..
+*       CHARACTER          UPLO
+*       INTEGER            LDA, LDAFAC, LDC, N
+*       COMPLEX            RESID
+*       ..
+*       .. Array Arguments ..
+*       INTEGER            IPIV( * )
+*       COMPLEX            A( LDA, * ), AFAC( LDAFAC, * ), C( LDC, * ),
+*      $                   RWORK( * )
+*       ..
+*
+*
+*> \par Purpose:
+*  =============
+*>
+*> \verbatim
+*>
+*> CHET01_AASEN reconstructs a 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 and EPS is the machine epsilon.
+*> \endverbatim
+*
+*  Arguments:
+*  ==========
+*
+*> \param[in] UPLO
+*> \verbatim
+*>          UPLO is CHARACTER*1
+*>          Specifies whether the upper or lower triangular part of the
+*>          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 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 CHETRF.
+*> \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 CHETRF.
+*> \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 COMPLEX array, dimension (N)
+*> \endverbatim
+*>
+*> \param[out] RESID
+*> \verbatim
+*>          RESID is COMPLEX
+*>          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 November 2016
+*
+*
+*> \ingroup complex_lin
+*
+*  =====================================================================
+      SUBROUTINE CHET01_AASEN( UPLO, N, A, LDA, AFAC, LDAFAC, IPIV, C,
+     $                         LDC, RWORK, RESID )
+*
+*  -- LAPACK test routine (version 3.7.0) --
+*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*     November 2016
+*
+*     .. Scalar Arguments ..
+      CHARACTER          UPLO
+      INTEGER            LDA, LDAFAC, LDC, N
+      REAL               RESID
+*     ..
+*     .. Array Arguments ..
+      INTEGER            IPIV( * )
+      COMPLEX            A( LDA, * ), AFAC( LDAFAC, * ), C( LDC, * ),
+     $                   RWORK( * )
+*     ..
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      COMPLEX         CZERO, CONE
+      PARAMETER          ( CZERO = ( 0.0E+0, 0.0E+0 ), 
+     $                     CONE  = ( 1.0E+0, 0.0E+0 ) )
+      REAL               ZERO, ONE
+      PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
+*     ..
+*     .. Local Scalars ..
+      INTEGER            I, J
+      REAL               ANORM, EPS
+*     ..
+*     .. External Functions ..
+      LOGICAL            LSAME
+      REAL               SLAMCH, CLANHE
+      EXTERNAL           LSAME, SLAMCH, CLANHE
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           CLASET, CLAVHE
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          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 = SLAMCH( 'Epsilon' )
+      ANORM = CLANHE( '1', UPLO, N, A, LDA, RWORK )
+*
+*     Initialize C to the tridiagonal matrix T.
+*
+      CALL CLASET( 'Full', N, N, CZERO, CZERO, C, LDC )
+      CALL CLACPY( 'F', 1, N, AFAC( 1, 1 ), LDAFAC+1, C( 1, 1 ), LDC+1 )
+      IF( N.GT.1 ) THEN
+         IF( LSAME( UPLO, 'U' ) ) THEN
+            CALL CLACPY( 'F', 1, N-1, AFAC( 1, 2 ), LDAFAC+1, C( 1, 2 ),
+     $                   LDC+1 )
+            CALL CLACPY( 'F', 1, N-1, AFAC( 1, 2 ), LDAFAC+1, C( 2, 1 ),
+     $                   LDC+1 )
+            CALL CLACGV( N-1, C( 2, 1 ), LDC+1 )
+         ELSE
+            CALL CLACPY( 'F', 1, N-1, AFAC( 2, 1 ), LDAFAC+1, C( 1, 2 ),
+     $                   LDC+1 )
+            CALL CLACPY( 'F', 1, N-1, AFAC( 2, 1 ), LDAFAC+1, C( 2, 1 ),
+     $                   LDC+1 )
+            CALL CLACGV( N-1, C( 1, 2 ), LDC+1 )
+         ENDIF
+      ENDIF
+*
+*     Call CTRMM to form the product U' * D (or L * D ).
+*
+      IF( LSAME( UPLO, 'U' ) ) THEN
+         CALL CTRMM( 'Left', UPLO, 'Conjugate transpose', 'Unit', N-1,
+     $               N, CONE, AFAC( 1, 2 ), LDAFAC, C( 2, 1 ), LDC )
+      ELSE
+         CALL CTRMM( 'Left', UPLO, 'No transpose', 'Unit', N-1, N,
+     $               CONE, AFAC( 2, 1 ), LDAFAC, C( 2, 1 ), LDC )
+      END IF
+*
+*     Call CTRMM again to multiply by U (or L ).
+*
+      IF( LSAME( UPLO, 'U' ) ) THEN
+         CALL CTRMM( 'Right', UPLO, 'No transpose', 'Unit', N, N-1,
+     $               CONE, AFAC( 1, 2 ), LDAFAC, C( 1, 2 ), LDC )
+      ELSE
+         CALL CTRMM( 'Right', UPLO, 'Conjugate transpose', 'Unit', N,
+     $               N-1, CONE, AFAC( 2, 1 ), LDAFAC, C( 1, 2 ), LDC )
+      END IF
+*
+*     Apply hermitian pivots
+*
+      DO J = N, 1, -1
+         I = IPIV( J )
+         IF( I.NE.J )
+     $      CALL CSWAP( N, C( J, 1 ), LDC, C( I, 1 ), LDC )
+      END DO
+      DO J = N, 1, -1
+         I = IPIV( J )
+         IF( I.NE.J )
+     $      CALL CSWAP( N, C( 1, J ), 1, C( 1, I ), 1 )
+      END DO
+*
+*
+*     Compute the difference  C - A .
+*
+      IF( LSAME( UPLO, 'U' ) ) THEN
+         DO J = 1, N
+            DO I = 1, J
+               C( I, J ) = C( I, J ) - A( I, J )
+            END DO
+         END DO
+      ELSE
+         DO J = 1, N
+            DO I = J, N
+               C( I, J ) = C( I, J ) - A( I, J )
+            END DO
+         END DO
+      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 / DBLE( N ) ) / ANORM ) / EPS
+      END IF
+*
+      RETURN
+*
+*     End of CHET01_AASEN
+*
+      END
diff --git a/TESTING/LIN/dchkaa.f b/TESTING/LIN/dchkaa.f
index 49478c93..9e7a14ab 100644
--- a/TESTING/LIN/dchkaa.f
+++ b/TESTING/LIN/dchkaa.f
@@ -49,6 +49,7 @@
 *> DPP    9               List types on next line if 0 < NTYPES <  9
 *> DPB    8               List types on next line if 0 < NTYPES <  8
 *> DPT   12               List types on next line if 0 < NTYPES < 12
+*> DSA   10               List types on next line if 0 < NTYPES < 10
 *> DSY   10               List types on next line if 0 < NTYPES < 10
 *> DSR   10               List types on next line if 0 < NTYPES < 10
 *> DSP   10               List types on next line if 0 < NTYPES < 10
@@ -664,6 +665,34 @@
             WRITE( NOUT, FMT = 9988 )PATH
          END IF
 *
+      ELSE IF( LSAMEN( 2, C2, 'SA' ) ) THEN
+*
+*        SY:  symmetric indefinite matrices,
+*             with partial (Aasen's) pivoting algorithm
+*
+         NTYPES = 10
+         CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT )
+*
+         IF( TSTCHK ) THEN
+            CALL DCHKSY_AASEN( DOTYPE, NN, NVAL, NNB2, NBVAL2, NNS, 
+     $                         NSVAL, THRESH, TSTERR, LDA,
+     $                         A( 1, 1 ), A( 1, 2 ), A( 1, 3 ),  
+     $                         B( 1, 1 ), B( 1, 2 ), B( 1, 3 ),
+     $                         WORK, RWORK, IWORK, NOUT )
+         ELSE
+            WRITE( NOUT, FMT = 9989 )PATH
+         END IF
+*
+         IF( TSTDRV ) THEN
+            CALL DDRVSY_AASEN( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR,
+     $                         LDA, A( 1, 1 ), A( 1, 2 ), A( 1, 3 ),
+     $                         B( 1, 1 ), B( 1, 2 ), B( 1, 3 ),
+     $                         WORK, RWORK, IWORK, NOUT )
+         ELSE
+            WRITE( NOUT, FMT = 9988 )PATH
+         END IF
+*
+*
       ELSE IF( LSAMEN( 2, C2, 'SP' ) ) THEN
 *
 *        SP:  symmetric indefinite packed matrices,
diff --git a/TESTING/LIN/dchksy_aasen.f b/TESTING/LIN/dchksy_aasen.f
new file mode 100644
index 00000000..041ef754
--- /dev/null
+++ b/TESTING/LIN/dchksy_aasen.f
@@ -0,0 +1,577 @@
+*> \brief \b DCHKSY_AASEN
+*
+*  =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+*            http://www.netlib.org/lapack/explore-html/
+*
+*  Definition:
+*  ===========
+*
+*       SUBROUTINE DCHKSY_AASEN( 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   A( * ), AFAC( * ), AINV( * ), B( * ),
+*      $                   RWORK( * ), WORK( * ), X( * ), XACT( * )
+*       ..
+*
+*
+*> \par Purpose:
+*  =============
+*>
+*> \verbatim
+*>
+*> DCHKSY_AASEN tests DSYTRF_AASEN, -TRS_AASEN.
+*> \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 DOUBLE PRECISION array, dimension (NMAX*NMAX)
+*> \endverbatim
+*>
+*> \param[out] AFAC
+*> \verbatim
+*>          AFAC is DOUBLE PRECISION array, dimension (NMAX*NMAX)
+*> \endverbatim
+*>
+*> \param[out] AINV
+*> \verbatim
+*>          AINV is DOUBLE PRECISION array, dimension (NMAX*NMAX)
+*> \endverbatim
+*>
+*> \param[out] B
+*> \verbatim
+*>          B is DOUBLE PRECISION array, dimension (NMAX*NSMAX)
+*>          where NSMAX is the largest entry in NSVAL.
+*> \endverbatim
+*>
+*> \param[out] X
+*> \verbatim
+*>          X is DOUBLE PRECISION array, dimension (NMAX*NSMAX)
+*> \endverbatim
+*>
+*> \param[out] XACT
+*> \verbatim
+*>          XACT is DOUBLE PRECISION array, dimension (NMAX*NSMAX)
+*> \endverbatim
+*>
+*> \param[out] WORK
+*> \verbatim
+*>          WORK is DOUBLE PRECISION 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 November 2016
+*
+*
+*> \ingroup double_lin
+*
+*  =====================================================================
+      SUBROUTINE DCHKSY_AASEN( 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.7.0) --
+*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*     November 2016
+*
+      IMPLICIT NONE
+*
+*     .. Scalar Arguments ..
+      LOGICAL            TSTERR
+      INTEGER            NN, NNB, NNS, NMAX, NOUT
+      DOUBLE PRECISION   THRESH
+*     ..
+*     .. Array Arguments ..
+      LOGICAL            DOTYPE( * )
+      INTEGER            IWORK( * ), NBVAL( * ), NSVAL( * ), NVAL( * )
+      DOUBLE PRECISION   A( * ), AFAC( * ), AINV( * ), B( * ),
+     $                   RWORK( * ), WORK( * ), X( * ), XACT( * )
+*     ..
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      DOUBLE PRECISION   ZERO, ONE
+      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
+      INTEGER            NTYPES
+      PARAMETER          ( NTYPES = 10 )
+      INTEGER            NTESTS
+      PARAMETER          ( NTESTS = 9 )
+*     ..
+*     .. Local Scalars ..
+      LOGICAL            TRFCON, ZEROT
+      CHARACTER          DIST, TYPE, UPLO, XTYPE
+      CHARACTER*3        PATH, MATPATH
+      INTEGER            I, I1, I2, IMAT, IN, INB, INFO, IOFF, IRHS,
+     $                   IUPLO, IZERO, J, K, KL, KU, LDA, LWORK, MODE,
+     $                   N, NB, NERRS, NFAIL, NIMAT, NRHS, NRUN, NT
+      DOUBLE PRECISION   ANORM, CNDNUM, RCONDC
+*     ..
+*     .. Local Arrays ..
+      CHARACTER          UPLOS( 2 )
+      INTEGER            ISEED( 4 ), ISEEDY( 4 )
+      DOUBLE PRECISION   RESULT( NTESTS )
+*     ..
+*     .. External Functions ..
+      DOUBLE PRECISION   DGET06, DLANSY
+      EXTERNAL           DGET06, DLANSY
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           ALAERH, ALAHD, ALASUM, DERRSY, DGET04, DLACPY,
+     $                   DLARHS, DLATB4, DLATMS, DPOT02, DPOT03, DPOT05,
+     $                   DSYCON, DSYRFS, DSYT01, DSYTRF_AASEN,
+     $                   DSYTRI2, DSYTRS_AASEN, XLAENV
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          MAX, MIN
+*     ..
+*     .. 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.
+*
+*     Test path
+*
+      PATH( 1: 1 ) = 'Double precision'
+      PATH( 2: 3 ) = 'SA'
+*
+*     Path to generate matrices
+*
+      MATPATH( 1: 1 ) = 'Double precision'
+      MATPATH( 2: 3 ) = 'SY'
+      NRUN = 0
+      NFAIL = 0
+      NERRS = 0
+      DO 10 I = 1, 4
+         ISEED( I ) = ISEEDY( I )
+   10 CONTINUE
+*
+*     Test the error exits
+*
+      IF( TSTERR )
+     $   CALL DERRSY( 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 180 IN = 1, NN
+         N = NVAL( IN )
+         IF( N .GT. NMAX ) THEN
+            NFAIL = NFAIL + 1
+            WRITE(NOUT, 9995) 'M ', N, NMAX
+            GO TO 180
+         END IF
+         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 170 IMAT = 1, NIMAT
+*
+*           Do the tests only if DOTYPE( IMAT ) is true.
+*
+            IF( .NOT.DOTYPE( IMAT ) )
+     $         GO TO 170
+*
+*           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 170
+*
+*           Do first for UPLO = 'U', then for UPLO = 'L'
+*
+            DO 160 IUPLO = 1, 2
+               UPLO = UPLOS( IUPLO )
+*
+*              Begin generate the test matrix A.
+*
+*
+*              Set up parameters with DLATB4 for the matrix generator
+*              based on the type of matrix to be generated.
+*
+               CALL DLATB4( MATPATH, IMAT, N, N, TYPE, KL, KU,
+     $                      ANORM, MODE, CNDNUM, DIST )
+*
+*              Generate a matrix with DLATMS.
+*
+               SRNAMT = 'DLATMS'
+               CALL DLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE,
+     $                      CNDNUM, ANORM, KL, KU, UPLO, A, LDA, WORK,
+     $                      INFO )
+*
+*              Check error code from DLATMS and handle error.
+*
+               IF( INFO.NE.0 ) THEN
+                  CALL ALAERH( PATH, 'DLATMS', INFO, 0, UPLO, N, N, -1,
+     $                         -1, -1, IMAT, NFAIL, NERRS, NOUT )
+*
+*                    Skip all tests for this generated matrix
+*
+                  GO TO 160
+               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 ) = ZERO
+   20                   CONTINUE
+                        IOFF = IOFF + IZERO
+                        DO 30 I = IZERO, N
+                           A( IOFF ) = ZERO
+                           IOFF = IOFF + LDA
+   30                   CONTINUE
+                     ELSE
+                        IOFF = IZERO
+                        DO 40 I = 1, IZERO - 1
+                           A( IOFF ) = ZERO
+                           IOFF = IOFF + LDA
+   40                   CONTINUE
+                        IOFF = IOFF - IZERO
+                        DO 50 I = IZERO, N
+                           A( IOFF+I ) = ZERO
+   50                   CONTINUE
+                     END IF
+                  ELSE
+                     IF( IUPLO.EQ.1 ) THEN
+*
+*                       Set the first IZERO rows and columns to zero.
+*
+                        IOFF = 0
+                        DO 70 J = 1, N
+                           I2 = MIN( J, IZERO )
+                           DO 60 I = 1, I2
+                              A( IOFF+I ) = ZERO
+   60                      CONTINUE
+                           IOFF = IOFF + LDA
+   70                   CONTINUE
+                        IZERO = 1
+                     ELSE
+*
+*                       Set the last IZERO rows and columns to zero.
+*
+                        IOFF = 0
+                        DO 90 J = 1, N
+                           I1 = MAX( J, IZERO )
+                           DO 80 I = I1, N
+                              A( IOFF+I ) = ZERO
+   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 150 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 DLACPY( 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.
+*
+                  SRNAMT = 'DSYTRF_AASEN'
+                  LWORK = N*NB + N
+                  CALL DSYTRF_AASEN( UPLO, N, AFAC, LDA, IWORK, AINV, 
+     $                               LWORK, INFO )
+*
+*                 Adjust the expected value of INFO to account for
+*                 pivoting.
+*
+                  IF( IZERO.GT.0 ) THEN
+                     J = 1
+                     K = IZERO
+  100                CONTINUE
+                     IF( J.EQ.K ) THEN
+                        K = IWORK( J )
+                     ELSE IF( IWORK( J ).EQ.K ) THEN
+                        K = J
+                     END IF
+                     IF( J.LT.K ) THEN
+                        J = J + 1
+                        GO TO 100
+                     END IF
+                  ELSE
+                     K = 0
+                  END IF
+*
+*                 Check error code from DSYTRF and handle error.
+*
+                  IF( INFO.NE.K ) THEN
+                     CALL ALAERH( PATH, 'DSYTRF_AASEN', INFO, K, UPLO, 
+     $                            N, N, -1, -1, NB, IMAT, NFAIL, NERRS, 
+     $                            NOUT )
+                  END IF
+*
+*                 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 DSYT01_AASEN( UPLO, N, A, LDA, AFAC, LDA, IWORK,
+     $                               AINV, LDA, RWORK, RESULT( 1 ) )
+                  NT = 1
+*
+*
+*                 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
+*
+*                 Do only the condition estimate if INFO is not 0.
+*
+                  IF( TRFCON ) THEN
+                     RCONDC = ZERO
+                     GO TO 140
+                  END IF
+*
+*                 Do for each value of NRHS in NSVAL.
+*
+                  DO 130 IRHS = 1, NNS
+                     NRHS = NSVAL( IRHS )
+*
+*+    TEST 3 ( Using TRS)
+*                 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 = 'DLARHS'
+                     CALL DLARHS( MATPATH, XTYPE, UPLO, ' ', N, N,
+     $                            KL, KU, NRHS, A, LDA, XACT, LDA,
+     $                            B, LDA, ISEED, INFO )
+                     CALL DLACPY( 'Full', N, NRHS, B, LDA, X, LDA )
+*
+                     SRNAMT = 'DSYTRS_AASEN'
+                     LWORK = 3*N-2
+                     CALL DSYTRS_AASEN( UPLO, N, NRHS, AFAC, LDA, 
+     $                                  IWORK, X, LDA, WORK, LWORK,
+     $                                  INFO )
+*
+*                    Check error code from DSYTRS and handle error.
+*
+                     IF( INFO.NE.0 ) THEN
+                        CALL ALAERH( PATH, 'DSYTRS_AASEN', INFO, 0,
+     $                               UPLO, N, N, -1, -1, NRHS, IMAT,
+     $                               NFAIL, NERRS, NOUT )
+                     END IF
+*
+                     CALL DLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA )
+*
+*                    Compute the residual for the solution
+*
+                     CALL DPOT02( UPLO, N, NRHS, A, LDA, X, LDA, WORK,
+     $                            LDA, RWORK, RESULT( 2 ) )
+*
+*
+*                    Print information about the tests that did not pass
+*                    the threshold.
+*
+                     DO 120 K = 2, 2
+                        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
+  120                CONTINUE
+                     NRUN = NRUN + 1
+*
+*                 End do for each value of NRHS in NSVAL.
+*
+  130             CONTINUE
+  140             CONTINUE
+  150          CONTINUE
+  160       CONTINUE
+  170    CONTINUE
+  180 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 )
+ 9995 FORMAT( ' Invalid input value: ', A4, '=', I6, '; must be <=',
+     $      I6 )
+      RETURN
+*
+*     End of DCHKSY_AASEN
+*
+      END
diff --git a/TESTING/LIN/ddrvsy_aasen.f b/TESTING/LIN/ddrvsy_aasen.f
new file mode 100644
index 00000000..a3520cb1
--- /dev/null
+++ b/TESTING/LIN/ddrvsy_aasen.f
@@ -0,0 +1,517 @@
+*> \brief \b DDRVSY_AASEN
+*
+*  =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+*            http://www.netlib.org/lapack/explore-html/
+*
+*  Definition:
+*  ===========
+*
+*       SUBROUTINE DDRVSY_AASEN( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, NMAX,
+*                          A, AFAC, AINV, B, X, XACT, WORK, RWORK, IWORK,
+*                          NOUT )
+*
+*       .. Scalar Arguments ..
+*       LOGICAL            TSTERR
+*       INTEGER            NMAX, NN, NOUT, NRHS
+*       DOUBLE PRECISION   THRESH
+*       ..
+*       .. Array Arguments ..
+*       LOGICAL            DOTYPE( * )
+*       INTEGER            IWORK( * ), NVAL( * )
+*       DOUBLE PRECISION   A( * ), AFAC( * ), AINV( * ), B( * ),
+*      $                   RWORK( * ), WORK( * ), X( * ), XACT( * )
+*       ..
+*
+*
+*> \par Purpose:
+*  =============
+*>
+*> \verbatim
+*>
+*> DDRVSY_AASEN tests the driver routine DSYSV_AASEN.
+*> \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] NRHS
+*> \verbatim
+*>          NRHS is INTEGER
+*>          The number of right hand side vectors to be generated for
+*>          each linear system.
+*> \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 DOUBLE PRECISION array, dimension (NMAX*NMAX)
+*> \endverbatim
+*>
+*> \param[out] AFAC
+*> \verbatim
+*>          AFAC is DOUBLE PRECISION array, dimension (NMAX*NMAX)
+*> \endverbatim
+*>
+*> \param[out] AINV
+*> \verbatim
+*>          AINV is DOUBLE PRECISION array, dimension (NMAX*NMAX)
+*> \endverbatim
+*>
+*> \param[out] B
+*> \verbatim
+*>          B is DOUBLE PRECISION array, dimension (NMAX*NRHS)
+*> \endverbatim
+*>
+*> \param[out] X
+*> \verbatim
+*>          X is DOUBLE PRECISION array, dimension (NMAX*NRHS)
+*> \endverbatim
+*>
+*> \param[out] XACT
+*> \verbatim
+*>          XACT is DOUBLE PRECISION array, dimension (NMAX*NRHS)
+*> \endverbatim
+*>
+*> \param[out] WORK
+*> \verbatim
+*>          WORK is DOUBLE PRECISION array, dimension (NMAX*max(2,NRHS))
+*> \endverbatim
+*>
+*> \param[out] RWORK
+*> \verbatim
+*>          RWORK is DOUBLE PRECISION array, dimension (NMAX+2*NRHS)
+*> \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 November 2016
+*
+*> \ingroup double_lin
+*
+*  =====================================================================
+      SUBROUTINE DDRVSY_AASEN( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR,
+     $                         NMAX, A, AFAC, AINV, B, X, XACT, WORK,
+     $                         RWORK, IWORK, NOUT )
+*
+*  -- LAPACK test routine (version 3.7.0) --
+*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*     November 2016
+*
+*     .. Scalar Arguments ..
+      LOGICAL            TSTERR
+      INTEGER            NMAX, NN, NOUT, NRHS
+      DOUBLE PRECISION   THRESH
+*     ..
+*     .. Array Arguments ..
+      LOGICAL            DOTYPE( * )
+      INTEGER            IWORK( * ), NVAL( * )
+      DOUBLE PRECISION   A( * ), AFAC( * ), AINV( * ), B( * ),
+     $                   RWORK( * ), WORK( * ), X( * ), XACT( * )
+*     ..
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      DOUBLE PRECISION   ONE, ZERO
+      PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
+      INTEGER            NTYPES, NTESTS
+      PARAMETER          ( NTYPES = 10, NTESTS = 3 )
+      INTEGER            NFACT
+      PARAMETER          ( NFACT = 2 )
+*     ..
+*     .. Local Scalars ..
+      LOGICAL            ZEROT
+      CHARACTER          DIST, FACT, TYPE, UPLO, XTYPE
+      CHARACTER*3        MATPATH, PATH
+      INTEGER            I, I1, I2, IFACT, IMAT, IN, INFO, IOFF, IUPLO,
+     $                   IZERO, J, K, K1, KL, KU, LDA, LWORK, MODE, N,
+     $                   NB, NBMIN, NERRS, NFAIL, NIMAT, NRUN, NT
+      DOUBLE PRECISION   AINVNM, ANORM, CNDNUM, RCOND, RCONDC
+*     ..
+*     .. Local Arrays ..
+      CHARACTER          FACTS( NFACT ), UPLOS( 2 )
+      INTEGER            ISEED( 4 ), ISEEDY( 4 )
+      DOUBLE PRECISION   RESULT( NTESTS )
+*     ..
+*     .. External Functions ..
+      DOUBLE PRECISION   DGET06, DLANSY
+      EXTERNAL           DGET06, DLANSY
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           ALADHD, ALAERH, ALASVM, DERRVX, DGET04, DLACPY,
+     $                   DLARHS, DLASET, DLATB4, DLATMS, DPOT02, DPOT05,
+     $                   DSYSV_AASEN, DSYT01_AASEN, DSYTRF_AASEN, XLAENV
+*     ..
+*     .. Scalars in Common ..
+      LOGICAL            LERR, OK
+      CHARACTER*32       SRNAMT
+      INTEGER            INFOT, NUNIT
+*     ..
+*     .. Common blocks ..
+      COMMON             / INFOC / INFOT, NUNIT, OK, LERR
+      COMMON             / SRNAMC / SRNAMT
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          MAX, MIN
+*     ..
+*     .. Data statements ..
+      DATA               ISEEDY / 1988, 1989, 1990, 1991 /
+      DATA               UPLOS / 'U', 'L' / , FACTS / 'F', 'N' /
+*     ..
+*     .. Executable Statements ..
+*
+*     Initialize constants and the random number seed.
+*
+*     Test path
+*
+      PATH( 1: 1 ) = 'Double precision'
+      PATH( 2: 3 ) = 'SA'
+*
+*     Path to generate matrices
+*
+      MATPATH( 1: 1 ) = 'Double precision'
+      MATPATH( 2: 3 ) = 'SY'
+*
+      NRUN = 0
+      NFAIL = 0
+      NERRS = 0
+      DO 10 I = 1, 4
+         ISEED( I ) = ISEEDY( I )
+   10 CONTINUE
+      LWORK = MAX( 2*NMAX, NMAX*NRHS )
+*
+*     Test the error exits
+*
+      IF( TSTERR )
+     $   CALL DERRVX( PATH, NOUT )
+      INFOT = 0
+*
+*     Set the block size and minimum block size for testing.
+*
+      NB = 1
+      NBMIN = 2
+      CALL XLAENV( 1, NB )
+      CALL XLAENV( 2, NBMIN )
+*
+*     Do for each value of N in NVAL
+*
+      DO 180 IN = 1, NN
+         N = NVAL( IN )
+         LDA = MAX( N, 1 )
+         XTYPE = 'N'
+         NIMAT = NTYPES
+         IF( N.LE.0 )
+     $      NIMAT = 1
+*
+         DO 170 IMAT = 1, NIMAT
+*
+*           Do the tests only if DOTYPE( IMAT ) is true.
+*
+            IF( .NOT.DOTYPE( IMAT ) )
+     $         GO TO 170
+*
+*           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 170
+*
+*           Do first for UPLO = 'U', then for UPLO = 'L'
+*
+            DO 160 IUPLO = 1, 2
+               UPLO = UPLOS( IUPLO )
+*
+*              Set up parameters with DLATB4 and generate a test matrix
+*              with DLATMS.
+*
+               CALL DLATB4( MATPATH, IMAT, N, N, TYPE, KL, KU, ANORM,
+     $                      MODE, CNDNUM, DIST )
+*
+               SRNAMT = 'DLATMS'
+               CALL DLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE,
+     $                      CNDNUM, ANORM, KL, KU, UPLO, A, LDA, WORK,
+     $                      INFO )
+*
+*              Check error code from DLATMS.
+*
+               IF( INFO.NE.0 ) THEN
+                  CALL ALAERH( PATH, 'DLATMS', INFO, 0, UPLO, N, N, -1,
+     $                         -1, -1, IMAT, NFAIL, NERRS, NOUT )
+                  GO TO 160
+               END IF
+*
+*              For 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 ) = ZERO
+   20                   CONTINUE
+                        IOFF = IOFF + IZERO
+                        DO 30 I = IZERO, N
+                           A( IOFF ) = ZERO
+                           IOFF = IOFF + LDA
+   30                   CONTINUE
+                     ELSE
+                        IOFF = IZERO
+                        DO 40 I = 1, IZERO - 1
+                           A( IOFF ) = ZERO
+                           IOFF = IOFF + LDA
+   40                   CONTINUE
+                        IOFF = IOFF - IZERO
+                        DO 50 I = IZERO, N
+                           A( IOFF+I ) = ZERO
+   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 ) = ZERO
+   60                      CONTINUE
+                           IOFF = IOFF + LDA
+   70                   CONTINUE
+                        IZERO = 1
+                     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 ) = ZERO
+   80                      CONTINUE
+                           IOFF = IOFF + LDA
+   90                   CONTINUE
+                     END IF
+                  END IF
+               ELSE
+                  IZERO = 0
+               END IF
+*
+               DO 150 IFACT = 1, NFACT
+*
+*                 Do first for FACT = 'F', then for other values.
+*
+                  FACT = FACTS( IFACT )
+*
+*                 Compute the condition number for comparison with
+*                 the value returned by DSYSVX.
+*
+                  IF( ZEROT ) THEN
+                     IF( IFACT.EQ.1 )
+     $                  GO TO 150
+                     RCONDC = ZERO
+*
+                  ELSE IF( IFACT.EQ.1 ) THEN
+*
+*                    Compute the 1-norm of A.
+*
+                     ANORM = DLANSY( '1', UPLO, N, A, LDA, RWORK )
+*
+*                    Factor the matrix A.
+*
+c                     CALL DLACPY( UPLO, N, N, A, LDA, AFAC, LDA )
+c                     CALL DSYTRF( UPLO, N, AFAC, LDA, IWORK, WORK,
+c     $                            LWORK, INFO )
+*
+*                    Compute inv(A) and take its norm.
+*
+c                     CALL DLACPY( UPLO, N, N, AFAC, LDA, AINV, LDA )
+c                     LWORK = (N+NB+1)*(NB+3)
+c                     SRNAMT = 'DSYTRI2'
+c                     CALL DSYTRI2( UPLO, N, AINV, LDA, IWORK, WORK,
+c     $                            LWORK, INFO )
+c                     AINVNM = DLANSY( '1', UPLO, N, AINV, LDA, RWORK )
+*
+*                    Compute the 1-norm condition number of A.
+*
+c                     IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN
+c                        RCONDC = ONE
+c                     ELSE
+c                        RCONDC = ( ONE / ANORM ) / AINVNM
+c                     END IF
+                  END IF
+*
+*                 Form an exact solution and set the right hand side.
+*
+                  SRNAMT = 'DLARHS'
+                  CALL DLARHS( MATPATH, XTYPE, UPLO, ' ', N, N, KL, KU,
+     $                         NRHS, A, LDA, XACT, LDA, B, LDA, ISEED,
+     $                         INFO )
+                  XTYPE = 'C'
+*
+*                 --- Test DSYSV_AASEN  ---
+*
+                  IF( IFACT.EQ.2 ) THEN
+                     CALL DLACPY( UPLO, N, N, A, LDA, AFAC, LDA )
+                     CALL DLACPY( 'Full', N, NRHS, B, LDA, X, LDA )
+*
+*                    Factor the matrix and solve the system using DSYSV_AASEN.
+*
+                     SRNAMT = 'DSYSV_AASEN'
+                     CALL DSYSV_AASEN( UPLO, N, NRHS, AFAC, LDA, IWORK,
+     $                                 X, LDA, WORK, LWORK, INFO )
+*
+*                    Adjust the expected value of INFO to account for
+*                    pivoting.
+*
+                     IF( IZERO.GT.0 ) THEN
+                        J = 1
+                        K = IZERO
+  100                   CONTINUE
+                        IF( J.EQ.K ) THEN
+                           K = IWORK( J )
+                        ELSE IF( IWORK( J ).EQ.K ) THEN
+                           K = J
+                        END IF
+                        IF( J.LT.K ) THEN
+                           J = J + 1
+                           GO TO 100
+                        END IF
+                     ELSE
+                        K = 0
+                     END IF
+*
+*                    Check error code from DSYSV_AASEN .
+*
+                     IF( INFO.NE.K ) THEN
+                        CALL ALAERH( PATH, 'DSYSV_AASEN ', INFO, K,
+     $                               UPLO, N, N, -1, -1, NRHS,
+     $                               IMAT, NFAIL, NERRS, NOUT )
+                        GO TO 120
+                     ELSE IF( INFO.NE.0 ) THEN
+                        GO TO 120
+                     END IF
+*
+*                    Reconstruct matrix from factors and compute
+*                    residual.
+*
+                     CALL DSYT01_AASEN( UPLO, N, A, LDA, AFAC, LDA,
+     $                                  IWORK, AINV, LDA, RWORK, 
+     $                                  RESULT( 1 ) )
+*
+*                    Compute residual of the computed solution.
+*
+                     CALL DLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA )
+                     CALL DPOT02( UPLO, N, NRHS, A, LDA, X, LDA, WORK,
+     $                            LDA, RWORK, RESULT( 2 ) )
+*
+*                    Check solution from generated exact solution.
+*
+                     CALL DGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
+     $                            RESULT( 3 ) )
+                     NT = 3
+*
+*                    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 ALADHD( NOUT, PATH )
+                           WRITE( NOUT, FMT = 9999 )'DSYSV_AASEN ',
+     $                        UPLO, N, IMAT, K, RESULT( K )
+                           NFAIL = NFAIL + 1
+                        END IF
+  110                CONTINUE
+                     NRUN = NRUN + NT
+  120                CONTINUE
+                  END IF
+*
+  150          CONTINUE
+*
+  160       CONTINUE
+  170    CONTINUE
+  180 CONTINUE
+*
+*     Print a summary of the results.
+*
+      CALL ALASVM( PATH, NOUT, NFAIL, NRUN, NERRS )
+*
+ 9999 FORMAT( 1X, A, ', UPLO=''', A1, ''', N =', I5, ', type ', I2,
+     $      ', test ', I2, ', ratio =', G12.5 )
+      RETURN
+*
+*     End of DDRVSY_AASEN
+*
+      END
diff --git a/TESTING/LIN/derrvx.f b/TESTING/LIN/derrvx.f
index 36eb8d6a..2affbe2f 100644
--- a/TESTING/LIN/derrvx.f
+++ b/TESTING/LIN/derrvx.f
@@ -91,7 +91,7 @@
       EXTERNAL           CHKXER, DGBSV, DGBSVX, DGESV, DGESVX, DGTSV,
      $                   DGTSVX, DPBSV, DPBSVX, DPOSV, DPOSVX, DPPSV,
      $                   DPPSVX, DPTSV, DPTSVX, DSPSV, DSPSVX, DSYSV,
-     $                   DSYSV_ROOK, DSYSVX
+     $                   DSYSV_AASEN, DSYSV_ROOK, DSYSVX
 *     ..
 *     .. Scalars in Common ..
       LOGICAL            LERR, OK
@@ -627,6 +627,25 @@
      $                RCOND, R1, R2, W, 3, IW, INFO )
          CALL CHKXER( 'DSYSVX', INFOT, NOUT, LERR, OK )
 *
+      ELSE IF( LSAMEN( 2, C2, 'SA' ) ) THEN
+*
+*        DSYSV_AASEN
+*
+        SRNAMT = 'DSYSV_AASEN'
+        INFOT = 1
+        CALL DSYSV_AASEN( '/', 0, 0, A, 1, IP, B, 1, W, 1, INFO )
+        CALL CHKXER( 'DSYSV_AASEN', INFOT, NOUT, LERR, OK )
+        INFOT = 2
+        CALL DSYSV_AASEN( 'U', -1, 0, A, 1, IP, B, 1, W, 1, INFO )
+        CALL CHKXER( 'DSYSV_AASEN', INFOT, NOUT, LERR, OK )
+        INFOT = 3
+        CALL DSYSV_AASEN( 'U', 0, -1, A, 1, IP, B, 1, W, 1, INFO )
+        CALL CHKXER( 'DSYSV_AASEN', INFOT, NOUT, LERR, OK )
+        INFOT = 8
+        CALL DSYSV_AASEN( 'U', 2, 0, A, 2, IP, B, 1, W, 1, INFO )
+        CALL CHKXER( 'DSYSV_AASEN', INFOT, NOUT, LERR, OK )
+*
+
       ELSE IF( LSAMEN( 2, C2, 'SR' ) ) THEN
 *
 *        DSYSV_ROOK
diff --git a/TESTING/LIN/dsyt01_aasen.f b/TESTING/LIN/dsyt01_aasen.f
new file mode 100644
index 00000000..f3756771
--- /dev/null
+++ b/TESTING/LIN/dsyt01_aasen.f
@@ -0,0 +1,263 @@
+*> \brief \b DSYT01
+*
+*  =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+*            http://www.netlib.org/lapack/explore-html/
+*
+*  Definition:
+*  ===========
+*
+*       SUBROUTINE DSYT01( 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   A( LDA, * ), AFAC( LDAFAC, * ), C( LDC, * ),
+*      $                   RWORK( * )
+*       ..
+*
+*
+*> \par Purpose:
+*  =============
+*>
+*> \verbatim
+*>
+*> DSYT01 reconstructs a symmetric 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 and EPS is the machine epsilon.
+*> \endverbatim
+*
+*  Arguments:
+*  ==========
+*
+*> \param[in] UPLO
+*> \verbatim
+*>          UPLO is CHARACTER*1
+*>          Specifies whether the upper or lower triangular part of the
+*>          symmetric matrix A is stored:
+*>          = 'U':  Upper triangular
+*>          = 'L':  Lower triangular
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*>          N is INTEGER
+*>          The number of rows and columns of the matrix A.  N >= 0.
+*> \endverbatim
+*>
+*> \param[in] A
+*> \verbatim
+*>          A is DOUBLE PRECISION array, dimension (LDA,N)
+*>          The original symmetric 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 DOUBLE PRECISION 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 DSYTRF.
+*> \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 DSYTRF.
+*> \endverbatim
+*>
+*> \param[out] C
+*> \verbatim
+*>          C is DOUBLE PRECISION 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 November 2016
+*
+*  @precisions fortran d -> s
+*
+*> \ingroup double_lin
+*
+*  =====================================================================
+      SUBROUTINE DSYT01_AASEN( UPLO, N, A, LDA, AFAC, LDAFAC, IPIV, C,
+     $                         LDC, RWORK, RESID )
+*
+*  -- LAPACK test routine (version 3.5.0) --
+*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*     November 2016
+*
+*     .. Scalar Arguments ..
+      CHARACTER          UPLO
+      INTEGER            LDA, LDAFAC, LDC, N
+      DOUBLE PRECISION   RESID
+*     ..
+*     .. Array Arguments ..
+      INTEGER            IPIV( * )
+      DOUBLE PRECISION   A( LDA, * ), AFAC( LDAFAC, * ), C( LDC, * ),
+     $                   RWORK( * )
+*     ..
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      DOUBLE PRECISION   ZERO, ONE
+      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
+*     ..
+*     .. Local Scalars ..
+      INTEGER            I, J
+      DOUBLE PRECISION   ANORM, EPS
+*     ..
+*     .. External Functions ..
+      LOGICAL            LSAME
+      DOUBLE PRECISION   DLAMCH, DLANSY
+      EXTERNAL           LSAME, DLAMCH, DLANSY
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           DLASET, DLAVSY
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          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 = DLANSY( '1', UPLO, N, A, LDA, RWORK )
+*
+*     Initialize C to the tridiagonal matrix T.
+*
+      CALL DLASET( 'Full', N, N, ZERO, ZERO, C, LDC )
+      CALL DLACPY( 'F', 1, N, AFAC( 1, 1 ), LDAFAC+1, C( 1, 1 ), LDC+1 )
+      IF( N.GT.1 ) THEN
+         IF( LSAME( UPLO, 'U' ) ) THEN
+            CALL DLACPY( 'F', 1, N-1, AFAC( 1, 2 ), LDAFAC+1, C( 1, 2 ),
+     $                   LDC+1 )
+            CALL DLACPY( 'F', 1, N-1, AFAC( 1, 2 ), LDAFAC+1, C( 2, 1 ),
+     $                   LDC+1 )
+         ELSE
+            CALL DLACPY( 'F', 1, N-1, AFAC( 2, 1 ), LDAFAC+1, C( 1, 2 ),
+     $                   LDC+1 )
+            CALL DLACPY( 'F', 1, N-1, AFAC( 2, 1 ), LDAFAC+1, C( 2, 1 ),
+     $                   LDC+1 )
+         ENDIF
+      ENDIF
+*
+*     Call DTRMM to form the product U' * D (or L * D ).
+*
+      IF( LSAME( UPLO, 'U' ) ) THEN
+         CALL DTRMM( 'Left', UPLO, 'Transpose', 'Unit', N-1, N,
+     $               ONE, AFAC( 1, 2 ), LDAFAC, C( 2, 1 ), LDC )
+      ELSE
+         CALL DTRMM( 'Left', UPLO, 'No transpose', 'Unit', N-1, N,
+     $               ONE, AFAC( 2, 1 ), LDAFAC, C( 2, 1 ), LDC )
+      END IF
+*
+*     Call DTRMM again to multiply by U (or L ).
+*
+      IF( LSAME( UPLO, 'U' ) ) THEN
+         CALL DTRMM( 'Right', UPLO, 'No transpose', 'Unit', N, N-1,
+     $               ONE, AFAC( 1, 2 ), LDAFAC, C( 1, 2 ), LDC )
+      ELSE
+         CALL DTRMM( 'Right', UPLO, 'Transpose', 'Unit', N, N-1,
+     $               ONE, AFAC( 2, 1 ), LDAFAC, C( 1, 2 ), LDC )
+      END IF
+*
+*     Apply symmetric pivots
+*
+      DO J = N, 1, -1
+         I = IPIV( J )
+         IF( I.NE.J )
+     $      CALL DSWAP( N, C( J, 1 ), LDC, C( I, 1 ), LDC )
+      END DO
+      DO J = N, 1, -1
+         I = IPIV( J )
+         IF( I.NE.J )
+     $      CALL DSWAP( N, C( 1, J ), 1, C( 1, I ), 1 )
+      END DO
+*
+*
+*     Compute the difference  C - A .
+*
+      IF( LSAME( UPLO, 'U' ) ) THEN
+         DO J = 1, N
+            DO I = 1, J
+               C( I, J ) = C( I, J ) - A( I, J )
+            END DO
+         END DO
+      ELSE
+         DO J = 1, N
+            DO I = J, N
+               C( I, J ) = C( I, J ) - A( I, J )
+            END DO
+         END DO
+      END IF
+*
+*     Compute norm( C - A ) / ( N * norm(A) * EPS )
+*
+      RESID = DLANSY( '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 DSYT01
+*
+      END
diff --git a/TESTING/LIN/schkaa.f b/TESTING/LIN/schkaa.f
index 821710ce..a6bfa101 100644
--- a/TESTING/LIN/schkaa.f
+++ b/TESTING/LIN/schkaa.f
@@ -664,6 +664,33 @@
             WRITE( NOUT, FMT = 9988 )PATH
          END IF
 *
+      ELSE IF( LSAMEN( 2, C2, 'SA' ) ) THEN
+*
+*        SY:  symmetric indefinite matrices,
+*             with partial (Aasen's) pivoting algorithm
+*
+         NTYPES = 10
+         CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT )
+*
+         IF( TSTCHK ) THEN
+            CALL SCHKSY_AASEN( DOTYPE, NN, NVAL, NNB2, NBVAL2, NNS,
+     $                         NSVAL, THRESH, TSTERR, LDA,
+     $                         A( 1, 1 ), A( 1, 2 ), A( 1, 3 ),
+     $                         B( 1, 1 ), B( 1, 2 ), B( 1, 3 ),
+     $                         WORK, RWORK, IWORK, NOUT )
+         ELSE
+            WRITE( NOUT, FMT = 9989 )PATH
+         END IF
+*
+         IF( TSTDRV ) THEN
+            CALL SDRVSY_AASEN( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR,
+     $                         LDA, A( 1, 1 ), A( 1, 2 ), A( 1, 3 ),
+     $                         B( 1, 1 ), B( 1, 2 ), B( 1, 3 ),
+     $                         WORK, RWORK, IWORK, NOUT )
+        ELSE
+           WRITE( NOUT, FMT = 9988 )PATH
+        END IF
+*
       ELSE IF( LSAMEN( 2, C2, 'SP' ) ) THEN
 *
 *        SP:  symmetric indefinite packed matrices,
diff --git a/TESTING/LIN/schksy_aasen.f b/TESTING/LIN/schksy_aasen.f
new file mode 100644
index 00000000..659e3fd3
--- /dev/null
+++ b/TESTING/LIN/schksy_aasen.f
@@ -0,0 +1,578 @@
+*> \brief \b SCHKSY_AASEN
+*
+*  =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+*            http://www.netlib.org/lapack/explore-html/
+*
+*  Definition:
+*  ===========
+*
+*       SUBROUTINE SCHKSY_AAEN( 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       A( * ), AFAC( * ), AINV( * ), B( * ),
+*      $           RWORK( * ), WORK( * ), X( * ), XACT( * )
+*       ..
+*
+*
+*> \par Purpose:
+*  =============
+*>
+*> \verbatim
+*>
+*> SCHKSY_AASEN tests SSYTRF_AASEN, -TRS_AASEN.
+*> \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 REAL array, dimension (NMAX*NMAX)
+*> \endverbatim
+*>
+*> \param[out] AFAC
+*> \verbatim
+*>          AFAC is REAL array, dimension (NMAX*NMAX)
+*> \endverbatim
+*>
+*> \param[out] AINV
+*> \verbatim
+*>          AINV is REAL array, dimension (NMAX*NMAX)
+*> \endverbatim
+*>
+*> \param[out] B
+*> \verbatim
+*>          B is REAL array, dimension (NMAX*NSMAX)
+*>          where NSMAX is the largest entry in NSVAL.
+*> \endverbatim
+*>
+*> \param[out] X
+*> \verbatim
+*>          X is REAL array, dimension (NMAX*NSMAX)
+*> \endverbatim
+*>
+*> \param[out] XACT
+*> \verbatim
+*>          XACT is REAL array, dimension (NMAX*NSMAX)
+*> \endverbatim
+*>
+*> \param[out] WORK
+*> \verbatim
+*>          WORK is REAL 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 November 2016
+*
+*
+*> \ingroup real_lin
+*
+*  =====================================================================
+      SUBROUTINE SCHKSY_AASEN( 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.7.0) --
+*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*     November 2016
+*
+      IMPLICIT NONE
+*
+*     .. Scalar Arguments ..
+      LOGICAL      TSTERR
+      INTEGER      NN, NNB, NNS, NMAX, NOUT
+      REAL         THRESH
+*     ..
+*     .. Array Arguments ..
+      LOGICAL      DOTYPE( * )
+      INTEGER      IWORK( * ), NBVAL( * ), NSVAL( * ), NVAL( * )
+      REAL         A( * ), AFAC( * ), AINV( * ), B( * ),
+     $             RWORK( * ), WORK( * ), X( * ), XACT( * )
+*     ..
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      REAL         ZERO
+      PARAMETER    ( ZERO = 0.0E+0 )
+      INTEGER      NTYPES
+      PARAMETER    ( NTYPES = 10 )
+      INTEGER      NTESTS
+      PARAMETER    ( NTESTS = 9 )
+*     ..
+*     .. Local Scalars ..
+      LOGICAL      TRFCON, ZEROT
+      CHARACTER    DIST, TYPE, UPLO, XTYPE
+      CHARACTER*3  PATH, MATPATH
+      INTEGER      I, I1, I2, IMAT, IN, INB, INFO, IOFF, IRHS,
+     $             IUPLO, IZERO, J, K, KL, KU, LDA, LWORK, MODE,
+     $             N, NB, NERRS, NFAIL, NIMAT, NRHS, NRUN, NT
+      REAL         ANORM, CNDNUM, RCONDC
+*     ..
+*     .. Local Arrays ..
+      CHARACTER    UPLOS( 2 )
+      INTEGER      ISEED( 4 ), ISEEDY( 4 )
+      REAL         RESULT( NTESTS )
+*     ..
+*     .. External Functions ..
+      REAL         DGET06, SLANSY
+      EXTERNAL     DGET06, SLANSY
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL     ALAERH, ALAHD, ALASUM, SERRSY, SGET04, SLACPY,
+     $             SLARHS, SLATB4, SLATMS, SPOT02, DPOT03, DPOT05,
+     $             DSYCON, SSYRFS, SSYT01_AASEN, SSYTRF_AASEN,
+     $             DSYTRI2, SSYTRS_AASEN, XLAENV
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC    MAX, MIN
+*     ..
+*     .. 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.
+*
+*
+*     Test path
+*
+      PATH( 1: 1 ) = 'Single precision'
+      PATH( 2: 3 ) = 'SA'
+*
+*     Path to generate matrices
+*
+      MATPATH( 1: 1 ) = 'Single precision'
+      MATPATH( 2: 3 ) = 'SY'
+      NRUN = 0
+      NFAIL = 0
+      NERRS = 0
+      DO 10 I = 1, 4
+         ISEED( I ) = ISEEDY( I )
+   10 CONTINUE
+*
+*     Test the error exits
+*
+      IF( TSTERR )
+     $   CALL SERRSY( 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 180 IN = 1, NN
+         N = NVAL( IN )
+         IF( N .GT. NMAX ) THEN
+            NFAIL = NFAIL + 1
+            WRITE(NOUT, 9995) 'M ', N, NMAX
+            GO TO 180
+         END IF
+         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 170 IMAT = 1, NIMAT
+*
+*           Do the tests only if DOTYPE( IMAT ) is true.
+*
+            IF( .NOT.DOTYPE( IMAT ) )
+     $         GO TO 170
+*
+*           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 170
+*
+*           Do first for UPLO = 'U', then for UPLO = 'L'
+*
+            DO 160 IUPLO = 1, 2
+               UPLO = UPLOS( IUPLO )
+*
+*              Begin generate the test matrix A.
+*
+*
+*              Set up parameters with SLATB4 for the matrix generator
+*              based on the type of matrix to be generated.
+*
+               CALL SLATB4( MATPATH, IMAT, N, N, TYPE, KL, KU,
+     $                      ANORM, MODE, CNDNUM, DIST )
+*
+*              Generate a matrix with SLATMS.
+*
+               SRNAMT = 'SLATMS'
+               CALL SLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE,
+     $                      CNDNUM, ANORM, KL, KU, UPLO, A, LDA, WORK,
+     $                      INFO )
+*
+*              Check error code from SLATMS and handle error.
+*
+               IF( INFO.NE.0 ) THEN
+                  CALL ALAERH( PATH, 'SLATMS', INFO, 0, UPLO, N, N, -1,
+     $                         -1, -1, IMAT, NFAIL, NERRS, NOUT )
+*
+*                    Skip all tests for this generated matrix
+*
+                  GO TO 160
+               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 ) = ZERO
+   20                   CONTINUE
+                        IOFF = IOFF + IZERO
+                        DO 30 I = IZERO, N
+                           A( IOFF ) = ZERO
+                           IOFF = IOFF + LDA
+   30                   CONTINUE
+                     ELSE
+                        IOFF = IZERO
+                        DO 40 I = 1, IZERO - 1
+                           A( IOFF ) = ZERO
+                           IOFF = IOFF + LDA
+   40                   CONTINUE
+                        IOFF = IOFF - IZERO
+                        DO 50 I = IZERO, N
+                           A( IOFF+I ) = ZERO
+   50                   CONTINUE
+                     END IF
+                  ELSE
+                     IF( IUPLO.EQ.1 ) THEN
+*
+*                       Set the first IZERO rows and columns to zero.
+*
+                        IOFF = 0
+                        DO 70 J = 1, N
+                           I2 = MIN( J, IZERO )
+                           DO 60 I = 1, I2
+                              A( IOFF+I ) = ZERO
+   60                      CONTINUE
+                           IOFF = IOFF + LDA
+   70                   CONTINUE
+                        IZERO = 1
+                     ELSE
+*
+*                       Set the last IZERO rows and columns to zero.
+*
+                        IOFF = 0
+                        DO 90 J = 1, N
+                           I1 = MAX( J, IZERO )
+                           DO 80 I = I1, N
+                              A( IOFF+I ) = ZERO
+   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 150 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 SLACPY( 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.
+*
+                  SRNAMT = 'SSYTRF_AASEN'
+                  LWORK = N*NB + N
+                  CALL SSYTRF_AASEN( UPLO, N, AFAC, LDA, IWORK, AINV, 
+     $                               LWORK, INFO )
+*
+*                 Adjust the expected value of INFO to account for
+*                 pivoting.
+*
+                  IF( IZERO.GT.0 ) THEN
+                     J = 1
+                     K = IZERO
+  100                CONTINUE
+                     IF( J.EQ.K ) THEN
+                        K = IWORK( J )
+                     ELSE IF( IWORK( J ).EQ.K ) THEN
+                        K = J
+                     END IF
+                     IF( J.LT.K ) THEN
+                        J = J + 1
+                        GO TO 100
+                     END IF
+                  ELSE
+                     K = 0
+                  END IF
+*
+*                 Check error code from SSYTRF and handle error.
+*
+                  IF( INFO.NE.K ) THEN
+                     CALL ALAERH( PATH, 'SSYTRF_AASEN', INFO, K, UPLO, 
+     $                            N, N, -1, -1, NB, IMAT, NFAIL, NERRS, 
+     $                            NOUT )
+                  END IF
+*
+*                 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 SSYT01_AASEN( UPLO, N, A, LDA, AFAC, LDA, IWORK,
+     $                               AINV, LDA, RWORK, RESULT( 1 ) )
+                  NT = 1
+*
+*
+*                 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
+*
+*                 Do only the condition estimate if INFO is not 0.
+*
+                  IF( TRFCON ) THEN
+                     RCONDC = ZERO
+                     GO TO 140
+                  END IF
+*
+*                 Do for each value of NRHS in NSVAL.
+*
+                  DO 130 IRHS = 1, NNS
+                     NRHS = NSVAL( IRHS )
+*
+*+    TEST 3 ( Using TRS)
+*                 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 = 'SLARHS'
+                     CALL SLARHS( MATPATH, XTYPE, UPLO, ' ', N, N,
+     $                            KL, KU, NRHS, A, LDA, XACT, LDA,
+     $                            B, LDA, ISEED, INFO )
+                     CALL SLACPY( 'Full', N, NRHS, B, LDA, X, LDA )
+*
+                     SRNAMT = 'SSYTRS_AASEN'
+                     LWORK = 3*N-2
+                     CALL SSYTRS_AASEN( UPLO, N, NRHS, AFAC, LDA, 
+     $                                  IWORK, X, LDA, WORK, LWORK,
+     $                                  INFO )
+*
+*                    Check error code from SSYTRS and handle error.
+*
+                     IF( INFO.NE.0 ) THEN
+                        CALL ALAERH( PATH, 'SSYTRS_AASEN', INFO, 0,
+     $                               UPLO, N, N, -1, -1, NRHS, IMAT, 
+     $                               NFAIL, NERRS, NOUT )
+                     END IF
+*
+                     CALL SLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA )
+*
+*                    Compute the residual for the solution
+*
+                     CALL SPOT02( UPLO, N, NRHS, A, LDA, X, LDA, WORK,
+     $                            LDA, RWORK, RESULT( 2 ) )
+*
+*
+*                    Print information about the tests that did not pass
+*                    the threshold.
+*
+                     DO 120 K = 2, 2
+                        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
+  120                CONTINUE
+                     NRUN = NRUN + 1
+*
+*                 End do for each value of NRHS in NSVAL.
+*
+  130             CONTINUE
+  140             CONTINUE
+  150          CONTINUE
+  160       CONTINUE
+  170    CONTINUE
+  180 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 )
+ 9995 FORMAT( ' Invalid input value: ', A4, '=', I6, '; must be <=',
+     $      I6 )
+      RETURN
+*
+*     End of SCHKSY_AASEN
+*
+      END
diff --git a/TESTING/LIN/sdrvsy_aasen.f b/TESTING/LIN/sdrvsy_aasen.f
new file mode 100644
index 00000000..afa5dc49
--- /dev/null
+++ b/TESTING/LIN/sdrvsy_aasen.f
@@ -0,0 +1,517 @@
+*> \brief \b SDRVSY_AASEN
+*
+*  =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+*            http://www.netlib.org/lapack/explore-html/
+*
+*  Definition:
+*  ===========
+*
+*       SUBROUTINE SDRVSY_AASEN( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, NMAX,
+*                          A, AFAC, AINV, B, X, XACT, WORK, RWORK, IWORK,
+*                          NOUT )
+*
+*       .. Scalar Arguments ..
+*       LOGICAL            TSTERR
+*       INTEGER            NMAX, NN, NOUT, NRHS
+*       REAL               THRESH
+*       ..
+*       .. Array Arguments ..
+*       LOGICAL            DOTYPE( * )
+*       INTEGER            IWORK( * ), NVAL( * )
+*       REAL               A( * ), AFAC( * ), AINV( * ), B( * ),
+*      $                   RWORK( * ), WORK( * ), X( * ), XACT( * )
+*       ..
+*
+*
+*> \par Purpose:
+*  =============
+*>
+*> \verbatim
+*>
+*> SDRVSY_AASEN tests the driver routine SSYSV_AASEN.
+*> \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] NRHS
+*> \verbatim
+*>          NRHS is INTEGER
+*>          The number of right hand side vectors to be generated for
+*>          each linear system.
+*> \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 REAL array, dimension (NMAX*NMAX)
+*> \endverbatim
+*>
+*> \param[out] AFAC
+*> \verbatim
+*>          AFAC is REAL array, dimension (NMAX*NMAX)
+*> \endverbatim
+*>
+*> \param[out] AINV
+*> \verbatim
+*>          AINV is REAL array, dimension (NMAX*NMAX)
+*> \endverbatim
+*>
+*> \param[out] B
+*> \verbatim
+*>          B is REAL array, dimension (NMAX*NRHS)
+*> \endverbatim
+*>
+*> \param[out] X
+*> \verbatim
+*>          X is REAL array, dimension (NMAX*NRHS)
+*> \endverbatim
+*>
+*> \param[out] XACT
+*> \verbatim
+*>          XACT is REAL array, dimension (NMAX*NRHS)
+*> \endverbatim
+*>
+*> \param[out] WORK
+*> \verbatim
+*>          WORK is REAL array, dimension (NMAX*max(2,NRHS))
+*> \endverbatim
+*>
+*> \param[out] RWORK
+*> \verbatim
+*>          RWORK is REAL array, dimension (NMAX+2*NRHS)
+*> \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 November 2016
+*
+*> \ingroup real_lin
+*
+*  =====================================================================
+      SUBROUTINE SDRVSY_AASEN( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR,
+     $                         NMAX, A, AFAC, AINV, B, X, XACT, WORK,
+     $                         RWORK, IWORK, NOUT )
+*
+*  -- LAPACK test routine (version 3.7.0) --
+*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*     November 2016
+*
+*     .. Scalar Arguments ..
+      LOGICAL            TSTERR
+      INTEGER            NMAX, NN, NOUT, NRHS
+      REAL               THRESH
+*     ..
+*     .. Array Arguments ..
+      LOGICAL            DOTYPE( * )
+      INTEGER            IWORK( * ), NVAL( * )
+      REAL               A( * ), AFAC( * ), AINV( * ), B( * ),
+     $                   RWORK( * ), WORK( * ), X( * ), XACT( * )
+*     ..
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      REAL               ONE, ZERO
+      PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
+      INTEGER            NTYPES, NTESTS
+      PARAMETER          ( NTYPES = 10, NTESTS = 3 )
+      INTEGER            NFACT
+      PARAMETER          ( NFACT = 2 )
+*     ..
+*     .. Local Scalars ..
+      LOGICAL            ZEROT
+      CHARACTER          DIST, FACT, TYPE, UPLO, XTYPE
+      CHARACTER*3        MATPATH, PATH
+      INTEGER            I, I1, I2, IFACT, IMAT, IN, INFO, IOFF, IUPLO,
+     $                   IZERO, J, K, K1, KL, KU, LDA, LWORK, MODE, N,
+     $                   NB, NBMIN, NERRS, NFAIL, NIMAT, NRUN, NT
+      REAL               AINVNM, ANORM, CNDNUM, RCOND, RCONDC
+*     ..
+*     .. Local Arrays ..
+      CHARACTER          FACTS( NFACT ), UPLOS( 2 )
+      INTEGER            ISEED( 4 ), ISEEDY( 4 )
+      REAL               RESULT( NTESTS )
+*     ..
+*     .. External Functions ..
+      REAL               DGET06, SLANSY
+      EXTERNAL           DGET06, SLANSY
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           ALADHD, ALAERH, ALASVM, SERRVX, SGET04, SLACPY,
+     $                   SLARHS, SLASET, SLATB4, SLATMS, SPOT02, DPOT05,
+     $                   SSYSV_AASEN, SSYT01_AASEN, SSYTRF_AASEN, XLAENV
+*     ..
+*     .. Scalars in Common ..
+      LOGICAL            LERR, OK
+      CHARACTER*32       SRNAMT
+      INTEGER            INFOT, NUNIT
+*     ..
+*     .. Common blocks ..
+      COMMON             / INFOC / INFOT, NUNIT, OK, LERR
+      COMMON             / SRNAMC / SRNAMT
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          MAX, MIN
+*     ..
+*     .. Data statements ..
+      DATA               ISEEDY / 1988, 1989, 1990, 1991 /
+      DATA               UPLOS / 'U', 'L' / , FACTS / 'F', 'N' /
+*     ..
+*     .. Executable Statements ..
+*
+*     Initialize constants and the random number seed.
+*
+*     Test path
+*
+      PATH( 1: 1 ) = 'Single precision'
+      PATH( 2: 3 ) = 'SA'
+*
+*     Path to generate matrices
+*
+      MATPATH( 1: 1 ) = 'Single precision'
+      MATPATH( 2: 3 ) = 'SY'
+*
+      NRUN = 0
+      NFAIL = 0
+      NERRS = 0
+      DO 10 I = 1, 4
+         ISEED( I ) = ISEEDY( I )
+   10 CONTINUE
+      LWORK = MAX( 2*NMAX, NMAX*NRHS )
+*
+*     Test the error exits
+*
+      IF( TSTERR )
+     $   CALL SERRVX( PATH, NOUT )
+      INFOT = 0
+*
+*     Set the block size and minimum block size for testing.
+*
+      NB = 1
+      NBMIN = 2
+      CALL XLAENV( 1, NB )
+      CALL XLAENV( 2, NBMIN )
+*
+*     Do for each value of N in NVAL
+*
+      DO 180 IN = 1, NN
+         N = NVAL( IN )
+         LDA = MAX( N, 1 )
+         XTYPE = 'N'
+         NIMAT = NTYPES
+         IF( N.LE.0 )
+     $      NIMAT = 1
+*
+         DO 170 IMAT = 1, NIMAT
+*
+*           Do the tests only if DOTYPE( IMAT ) is true.
+*
+            IF( .NOT.DOTYPE( IMAT ) )
+     $         GO TO 170
+*
+*           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 170
+*
+*           Do first for UPLO = 'U', then for UPLO = 'L'
+*
+            DO 160 IUPLO = 1, 2
+               UPLO = UPLOS( IUPLO )
+*
+*              Set up parameters with SLATB4 and generate a test matrix
+*              with SLATMS.
+*
+               CALL SLATB4( MATPATH, IMAT, N, N, TYPE, KL, KU, ANORM,
+     $                      MODE, CNDNUM, DIST )
+*
+               SRNAMT = 'SLATMS'
+               CALL SLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE,
+     $                      CNDNUM, ANORM, KL, KU, UPLO, A, LDA, WORK,
+     $                      INFO )
+*
+*              Check error code from SLATMS.
+*
+               IF( INFO.NE.0 ) THEN
+                  CALL ALAERH( PATH, 'SLATMS', INFO, 0, UPLO, N, N, -1,
+     $                         -1, -1, IMAT, NFAIL, NERRS, NOUT )
+                  GO TO 160
+               END IF
+*
+*              For 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 ) = ZERO
+   20                   CONTINUE
+                        IOFF = IOFF + IZERO
+                        DO 30 I = IZERO, N
+                           A( IOFF ) = ZERO
+                           IOFF = IOFF + LDA
+   30                   CONTINUE
+                     ELSE
+                        IOFF = IZERO
+                        DO 40 I = 1, IZERO - 1
+                           A( IOFF ) = ZERO
+                           IOFF = IOFF + LDA
+   40                   CONTINUE
+                        IOFF = IOFF - IZERO
+                        DO 50 I = IZERO, N
+                           A( IOFF+I ) = ZERO
+   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 ) = ZERO
+   60                      CONTINUE
+                           IOFF = IOFF + LDA
+   70                   CONTINUE
+                        IZERO = 1
+                     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 ) = ZERO
+   80                      CONTINUE
+                           IOFF = IOFF + LDA
+   90                   CONTINUE
+                     END IF
+                  END IF
+               ELSE
+                  IZERO = 0
+               END IF
+*
+               DO 150 IFACT = 1, NFACT
+*
+*                 Do first for FACT = 'F', then for other values.
+*
+                  FACT = FACTS( IFACT )
+*
+*                 Compute the condition number for comparison with
+*                 the value returned by SSYSVX.
+*
+                  IF( ZEROT ) THEN
+                     IF( IFACT.EQ.1 )
+     $                  GO TO 150
+                     RCONDC = ZERO
+*
+                  ELSE IF( IFACT.EQ.1 ) THEN
+*
+*                    Compute the 1-norm of A.
+*
+                     ANORM = SLANSY( '1', UPLO, N, A, LDA, RWORK )
+*
+*                    Factor the matrix A.
+*
+c                     CALL SLACPY( UPLO, N, N, A, LDA, AFAC, LDA )
+c                     CALL SSYTRF( UPLO, N, AFAC, LDA, IWORK, WORK,
+c     $                            LWORK, INFO )
+*
+*                    Compute inv(A) and take its norm.
+*
+c                     CALL SLACPY( UPLO, N, N, AFAC, LDA, AINV, LDA )
+c                     LWORK = (N+NB+1)*(NB+3)
+c                     SRNAMT = 'DSYTRI2'
+c                     CALL DSYTRI2( UPLO, N, AINV, LDA, IWORK, WORK,
+c     $                            LWORK, INFO )
+c                     AINVNM = SLANSY( '1', UPLO, N, AINV, LDA, RWORK )
+*
+*                    Compute the 1-norm condition number of A.
+*
+c                     IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN
+c                        RCONDC = ONE
+c                     ELSE
+c                        RCONDC = ( ONE / ANORM ) / AINVNM
+c                     END IF
+                  END IF
+*
+*                 Form an exact solution and set the right hand side.
+*
+                  SRNAMT = 'SLARHS'
+                  CALL SLARHS( MATPATH, XTYPE, UPLO, ' ', N, N, KL, KU,
+     $                         NRHS, A, LDA, XACT, LDA, B, LDA, ISEED,
+     $                         INFO )
+                  XTYPE = 'C'
+*
+*                 --- Test SSYSV_AASEN  ---
+*
+                  IF( IFACT.EQ.2 ) THEN
+                     CALL SLACPY( UPLO, N, N, A, LDA, AFAC, LDA )
+                     CALL SLACPY( 'Full', N, NRHS, B, LDA, X, LDA )
+*
+*                    Factor the matrix and solve the system using SSYSV_AASEN.
+*
+                     SRNAMT = 'SSYSV_AASEN'
+                     CALL SSYSV_AASEN( UPLO, N, NRHS, AFAC, LDA, IWORK,
+     $                                 X, LDA, WORK, LWORK, INFO )
+*
+*                    Adjust the expected value of INFO to account for
+*                    pivoting.
+*
+                     IF( IZERO.GT.0 ) THEN
+                        J = 1
+                        K = IZERO
+  100                   CONTINUE
+                        IF( J.EQ.K ) THEN
+                           K = IWORK( J )
+                        ELSE IF( IWORK( J ).EQ.K ) THEN
+                           K = J
+                        END IF
+                        IF( J.LT.K ) THEN
+                           J = J + 1
+                           GO TO 100
+                        END IF
+                     ELSE
+                        K = 0
+                     END IF
+*
+*                    Check error code from SSYSV_AASEN .
+*
+                     IF( INFO.NE.K ) THEN
+                        CALL ALAERH( PATH, 'SSYSV_AASEN ', INFO, K, 
+     $                               UPLO, N, N, -1, -1, NRHS, 
+     $                               IMAT, NFAIL, NERRS, NOUT )
+                        GO TO 120
+                     ELSE IF( INFO.NE.0 ) THEN
+                        GO TO 120
+                     END IF
+*
+*                    Reconstruct matrix from factors and compute
+*                    residual.
+*
+                     CALL SSYT01_AASEN( UPLO, N, A, LDA, AFAC, LDA,
+     $                                  IWORK, AINV, LDA, RWORK, 
+     $                                  RESULT( 1 ) )
+*
+*                    Compute residual of the computed solution.
+*
+                     CALL SLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA )
+                     CALL SPOT02( UPLO, N, NRHS, A, LDA, X, LDA, WORK,
+     $                            LDA, RWORK, RESULT( 2 ) )
+*
+*                    Check solution from generated exact solution.
+*
+                     CALL SGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
+     $                            RESULT( 3 ) )
+                     NT = 3
+*
+*                    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 ALADHD( NOUT, PATH )
+                           WRITE( NOUT, FMT = 9999 )'SSYSV_AASEN ',
+     $                        UPLO, N, IMAT, K, RESULT( K )
+                           NFAIL = NFAIL + 1
+                        END IF
+  110                CONTINUE
+                     NRUN = NRUN + NT
+  120                CONTINUE
+                  END IF
+*
+  150          CONTINUE
+*
+  160       CONTINUE
+  170    CONTINUE
+  180 CONTINUE
+*
+*     Print a summary of the results.
+*
+      CALL ALASVM( PATH, NOUT, NFAIL, NRUN, NERRS )
+*
+ 9999 FORMAT( 1X, A, ', UPLO=''', A1, ''', N =', I5, ', type ', I2,
+     $      ', test ', I2, ', ratio =', G12.5 )
+      RETURN
+*
+*     End of SDRVSY_AASEN
+*
+      END
diff --git a/TESTING/LIN/serrvx.f b/TESTING/LIN/serrvx.f
index 7d999c4a..41334f5c 100644
--- a/TESTING/LIN/serrvx.f
+++ b/TESTING/LIN/serrvx.f
@@ -91,7 +91,7 @@
       EXTERNAL           CHKXER, SGBSV, SGBSVX, SGESV, SGESVX, SGTSV,
      $                   SGTSVX, SPBSV, SPBSVX, SPOSV, SPOSVX, SPPSV,
      $                   SPPSVX, SPTSV, SPTSVX, SSPSV, SSPSVX, SSYSV,
-     $                   SSYSV_ROOK, SSYSVX
+     $                   SSYSV_AASEN, SSYSV_ROOK, SSYSVX
 *     ..
 *     .. Scalars in Common ..
       LOGICAL            LERR, OK
@@ -627,6 +627,24 @@
      $                RCOND, R1, R2, W, 3, IW, INFO )
          CALL CHKXER( 'SSYSVX', INFOT, NOUT, LERR, OK )
 *
+      ELSE IF( LSAMEN( 2, C2, 'SA' ) ) THEN
+*
+*        SSYSV_AASEN
+*
+        SRNAMT = 'SSYSV_AASEN'
+        INFOT = 1
+        CALL SSYSV_AASEN( '/', 0, 0, A, 1, IP, B, 1, W, 1, INFO )
+        CALL CHKXER( 'SSYSV_AASEN', INFOT, NOUT, LERR, OK )
+        INFOT = 2
+        CALL SSYSV_AASEN( 'U', -1, 0, A, 1, IP, B, 1, W, 1, INFO )
+        CALL CHKXER( 'SSYSV_AASEN', INFOT, NOUT, LERR, OK )
+        INFOT = 3
+        CALL SSYSV_AASEN( 'U', 0, -1, A, 1, IP, B, 1, W, 1, INFO )
+        CALL CHKXER( 'SSYSV_AASEN', INFOT, NOUT, LERR, OK )
+        INFOT = 8
+        CALL SSYSV_AASEN( 'U', 2, 0, A, 2, IP, B, 1, W, 1, INFO )
+        CALL CHKXER( 'SSYSV_AASEN', INFOT, NOUT, LERR, OK )
+*
       ELSE IF( LSAMEN( 2, C2, 'SR' ) ) THEN
 *
 *        SSYSV_ROOK
diff --git a/TESTING/LIN/ssyt01_aasen.f b/TESTING/LIN/ssyt01_aasen.f
new file mode 100644
index 00000000..01163938
--- /dev/null
+++ b/TESTING/LIN/ssyt01_aasen.f
@@ -0,0 +1,262 @@
+*> \brief \b SSYT01_AASEN
+*
+*  =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+*            http://www.netlib.org/lapack/explore-html/
+*
+*  Definition:
+*  ===========
+*
+*       SUBROUTINE SSYT01_AASEN( 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               A( LDA, * ), AFAC( LDAFAC, * ), C( LDC, * ),
+*      $                   RWORK( * )
+*       ..
+*
+*
+*> \par Purpose:
+*  =============
+*>
+*> \verbatim
+*>
+*> SSYT01_AASEN reconstructs a symmetric 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 and EPS is the machine epsilon.
+*> \endverbatim
+*
+*  Arguments:
+*  ==========
+*
+*> \param[in] UPLO
+*> \verbatim
+*>          UPLO is CHARACTER*1
+*>          Specifies whether the upper or lower triangular part of the
+*>          symmetric matrix A is stored:
+*>          = 'U':  Upper triangular
+*>          = 'L':  Lower triangular
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*>          N is INTEGER
+*>          The number of rows and columns of the matrix A.  N >= 0.
+*> \endverbatim
+*>
+*> \param[in] A
+*> \verbatim
+*>          A is REAL array, dimension (LDA,N)
+*>          The original symmetric 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 REAL 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 SSYTRF.
+*> \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 SSYTRF.
+*> \endverbatim
+*>
+*> \param[out] C
+*> \verbatim
+*>          C is REAL 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 November 2016
+*
+*
+*> \ingroup real_lin
+*
+*  =====================================================================
+      SUBROUTINE SSYT01_AASEN( UPLO, N, A, LDA, AFAC, LDAFAC, IPIV, C,
+     $                         LDC, RWORK, RESID )
+*
+*  -- LAPACK test routine (version 3.7.0) --
+*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*     November 2016
+*
+*     .. Scalar Arguments ..
+      CHARACTER          UPLO
+      INTEGER            LDA, LDAFAC, LDC, N
+      REAL               RESID
+*     ..
+*     .. Array Arguments ..
+      INTEGER            IPIV( * )
+      REAL               A( LDA, * ), AFAC( LDAFAC, * ), C( LDC, * ),
+     $                   RWORK( * )
+*     ..
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      REAL   ZERO, ONE
+      PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
+*     ..
+*     .. Local Scalars ..
+      INTEGER            I, J
+      REAL               ANORM, EPS
+*     ..
+*     .. External Functions ..
+      LOGICAL            LSAME
+      REAL               SLAMCH, SLANSY
+      EXTERNAL           LSAME, SLAMCH, SLANSY
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           SLASET, SLAVSY
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          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 = SLAMCH( 'Epsilon' )
+      ANORM = SLANSY( '1', UPLO, N, A, LDA, RWORK )
+*
+*     Initialize C to the tridiagonal matrix T.
+*
+      CALL SLASET( 'Full', N, N, ZERO, ZERO, C, LDC )
+      CALL SLACPY( 'F', 1, N, AFAC( 1, 1 ), LDAFAC+1, C( 1, 1 ), LDC+1 )
+      IF( N.GT.1 ) THEN
+         IF( LSAME( UPLO, 'U' ) ) THEN
+            CALL SLACPY( 'F', 1, N-1, AFAC( 1, 2 ), LDAFAC+1, C( 1, 2 ),
+     $                   LDC+1 )
+            CALL SLACPY( 'F', 1, N-1, AFAC( 1, 2 ), LDAFAC+1, C( 2, 1 ),
+     $                   LDC+1 )
+         ELSE
+            CALL SLACPY( 'F', 1, N-1, AFAC( 2, 1 ), LDAFAC+1, C( 1, 2 ),
+     $                   LDC+1 )
+            CALL SLACPY( 'F', 1, N-1, AFAC( 2, 1 ), LDAFAC+1, C( 2, 1 ),
+     $                   LDC+1 )
+         ENDIF
+      ENDIF
+*
+*     Call STRMM to form the product U' * D (or L * D ).
+*
+      IF( LSAME( UPLO, 'U' ) ) THEN
+         CALL STRMM( 'Left', UPLO, 'Transpose', 'Unit', N-1, N,
+     $               ONE, AFAC( 1, 2 ), LDAFAC, C( 2, 1 ), LDC )
+      ELSE
+         CALL STRMM( 'Left', UPLO, 'No transpose', 'Unit', N-1, N,
+     $               ONE, AFAC( 2, 1 ), LDAFAC, C( 2, 1 ), LDC )
+      END IF
+*
+*     Call STRMM again to multiply by U (or L ).
+*
+      IF( LSAME( UPLO, 'U' ) ) THEN
+         CALL STRMM( 'Right', UPLO, 'No transpose', 'Unit', N, N-1,
+     $               ONE, AFAC( 1, 2 ), LDAFAC, C( 1, 2 ), LDC )
+      ELSE
+         CALL STRMM( 'Right', UPLO, 'Transpose', 'Unit', N, N-1,
+     $               ONE, AFAC( 2, 1 ), LDAFAC, C( 1, 2 ), LDC )
+      END IF
+*
+*     Apply symmetric pivots
+*
+      DO J = N, 1, -1
+         I = IPIV( J )
+         IF( I.NE.J )
+     $      CALL SSWAP( N, C( J, 1 ), LDC, C( I, 1 ), LDC )
+      END DO
+      DO J = N, 1, -1
+         I = IPIV( J )
+         IF( I.NE.J )
+     $      CALL SSWAP( N, C( 1, J ), 1, C( 1, I ), 1 )
+      END DO
+*
+*
+*     Compute the difference  C - A .
+*
+      IF( LSAME( UPLO, 'U' ) ) THEN
+         DO J = 1, N
+            DO I = 1, J
+               C( I, J ) = C( I, J ) - A( I, J )
+            END DO
+         END DO
+      ELSE
+         DO J = 1, N
+            DO I = J, N
+               C( I, J ) = C( I, J ) - A( I, J )
+            END DO
+         END DO
+      END IF
+*
+*     Compute norm( C - A ) / ( N * norm(A) * EPS )
+*
+      RESID = SLANSY( '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 SSYT01_AASEN
+*
+      END
diff --git a/TESTING/LIN/zchkaa.f b/TESTING/LIN/zchkaa.f
index 9212afc9..90b98a2e 100644
--- a/TESTING/LIN/zchkaa.f
+++ b/TESTING/LIN/zchkaa.f
@@ -50,6 +50,7 @@
 *> ZPB    8               List types on next line if 0 < NTYPES <  8
 *> ZPT   12               List types on next line if 0 < NTYPES < 12
 *> ZHE   10               List types on next line if 0 < NTYPES < 10
+*> ZHA   10               List types on next line if 0 < NTYPES < 10
 *> ZHR   10               List types on next line if 0 < NTYPES < 10
 *> ZHP   10               List types on next line if 0 < NTYPES < 10
 *> ZSY   11               List types on next line if 0 < NTYPES < 11
@@ -102,17 +103,17 @@
 *> \author Univ. of Colorado Denver
 *> \author NAG Ltd.
 *
-*> \date November 2015
+*> \date November 2016
 *
 *> \ingroup complex16_lin
 *
 *  =====================================================================
       PROGRAM ZCHKAA
 *
-*  -- LAPACK test routine (version 3.6.0) --
+*  -- LAPACK test routine (version 3.7.0) --
 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
-*     November 2015
+*     November 2016
 *
 *  =====================================================================
 *
@@ -639,6 +640,33 @@
             WRITE( NOUT, FMT = 9988 )PATH
          END IF
 *
+      ELSE IF( LSAMEN( 2, C2, 'HA' ) ) THEN
+*
+*        HA:  Hermitian indefinite matrices,
+*             with partial (Aasen's) pivoting algorithm
+*
+         NTYPES = 10
+         CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT )
+*
+         IF( TSTCHK ) THEN
+            CALL ZCHKHE_AASEN( DOTYPE, NN, NVAL, NNB2, NBVAL2, NNS,
+     $                         NSVAL, THRESH, TSTERR, LDA, 
+     $                         A( 1, 1 ), A( 1, 2 ), A( 1, 3 ),
+     $                         B( 1, 1 ), B( 1, 2 ), B( 1, 3 ),
+     $                         WORK, RWORK, IWORK, NOUT )
+         ELSE
+            WRITE( NOUT, FMT = 9989 )PATH
+         END IF
+*
+         IF( TSTDRV ) THEN
+            CALL ZDRVHE_AASEN( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR,
+     $                         LDA, A( 1, 1 ), A( 1, 2 ), A( 1, 3 ), 
+     $                              B( 1, 1 ), B( 1, 2 ), B( 1, 3 ), 
+     $                         WORK, RWORK, IWORK, NOUT )
+         ELSE
+            WRITE( NOUT, FMT = 9988 )PATH
+         END IF
+*
       ELSE IF( LSAMEN( 2, C2, 'HR' ) ) THEN
 *
 *        HR:  Hermitian indefinite matrices,
diff --git a/TESTING/LIN/zchkhe_aasen.f b/TESTING/LIN/zchkhe_aasen.f
new file mode 100644
index 00000000..1ebe5ef4
--- /dev/null
+++ b/TESTING/LIN/zchkhe_aasen.f
@@ -0,0 +1,579 @@
+*> \brief \b ZCHKHE_AASEN
+*
+*  =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+*            http://www.netlib.org/lapack/explore-html/
+*
+*  Definition:
+*  ===========
+*
+*       SUBROUTINE ZCHKHE_AASEN( 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            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_AASEN tests ZHETRF_AASEN, -TRS_AASEN.
+*> \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 (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 November 2016
+*
+*
+*> \ingroup complex16_lin
+*
+*  =====================================================================
+      SUBROUTINE ZCHKHE_AASEN( 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.7.0) --
+*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*     November 2016
+*
+      IMPLICIT NONE
+*
+*     .. 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
+      PARAMETER          ( ZERO = 0.0D+0 )
+      COMPLEX*16         CZERO
+      PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ) )
+      INTEGER            NTYPES
+      PARAMETER          ( NTYPES = 10 )
+      INTEGER            NTESTS
+      PARAMETER          ( NTESTS = 9 )
+*     ..
+*     .. Local Scalars ..
+      LOGICAL            TRFCON, ZEROT
+      CHARACTER          DIST, TYPE, UPLO, XTYPE
+      CHARACTER*3        PATH, MATPATH
+      INTEGER            I, I1, I2, IMAT, IN, INB, INFO, IOFF, IRHS,
+     $                   IUPLO, IZERO, J, K, KL, KU, LDA, LWORK, MODE,
+     $                   N, NB, NERRS, NFAIL, NIMAT, NRHS, NRUN, NT
+      DOUBLE PRECISION   ANORM, CNDNUM, RCOND, RCONDC
+*     ..
+*     .. Local Arrays ..
+      CHARACTER          UPLOS( 2 )
+      INTEGER            ISEED( 4 ), ISEEDY( 4 )
+      DOUBLE PRECISION   RESULT( NTESTS )
+*     ..
+*     .. External Functions ..
+      DOUBLE PRECISION   DGET06, ZLANHE
+      EXTERNAL           DGET06, ZLANHE
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           ALAERH, ALAHD, ALASUM, XLAENV, ZERRHE, ZGET04,
+     $                   ZHECON, ZHERFS, ZHET01, ZHETRF_AASEN, ZHETRI2,
+     $                   ZHETRS_AASEN, ZLACPY, ZLAIPD, ZLARHS, ZLATB4, 
+     $                   ZLATMS, ZPOT02, ZPOT03, ZPOT05
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          MAX, MIN
+*     ..
+*     .. 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.
+*
+*     Test path
+*
+      PATH( 1: 1 ) = 'Zomplex precision'
+      PATH( 2: 3 ) = 'HA'
+*
+*     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 180 IN = 1, NN
+         N = NVAL( IN )
+         IF( N .GT. NMAX ) THEN
+            NFAIL = NFAIL + 1
+            WRITE(NOUT, 9995) 'M ', N, NMAX
+            GO TO 180
+         END IF
+         LDA = MAX( N, 1 )
+         XTYPE = 'N'
+         NIMAT = NTYPES
+         IF( N.LE.0 )
+     $      NIMAT = 1
+*
+         IZERO = 0
+         DO 170 IMAT = 1, NIMAT
+*
+*           Do the tests only if DOTYPE( IMAT ) is true.
+*
+            IF( .NOT.DOTYPE( IMAT ) )
+     $         GO TO 170
+*
+*           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 170
+*
+*           Do first for UPLO = 'U', then for UPLO = 'L'
+*
+            DO 160 IUPLO = 1, 2
+               UPLO = UPLOS( IUPLO )
+*
+*              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 160
+               END IF
+*
+*              For 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
+                     IF( IUPLO.EQ.1 ) THEN
+*
+*                       Set the first IZERO rows and columns to zero.
+*
+                        IOFF = 0
+                        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
+                        IZERO = 1
+                     ELSE
+*
+*                       Set the last IZERO rows and columns to zero.
+*
+                        IOFF = 0
+                        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 test matrix A.
+*
+*
+*              Set the imaginary part of the diagonals.
+*
+               CALL ZLAIPD( N, A, LDA+1, 0 )
+*
+*              Do for each value of NB in NBVAL
+*
+               DO 150 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 = ( NB+1 )*LDA
+                  SRNAMT = 'ZHETRF_AASEN'
+                  CALL ZHETRF_AASEN( UPLO, N, AFAC, LDA, IWORK, AINV, 
+     $                               LWORK, INFO )
+*
+*                 Adjust the expected value of INFO to account for
+*                 pivoting.
+*
+                  IF( IZERO.GT.0 ) THEN
+                     J = 1
+                     K = IZERO
+  100                CONTINUE
+                     IF( J.EQ.K ) THEN
+                        K = IWORK( J )
+                     ELSE IF( IWORK( J ).EQ.K ) THEN
+                        K = J
+                     END IF
+                     IF( J.LT.K ) THEN
+                        J = J + 1
+                        GO TO 100
+                     END IF
+                  ELSE
+                     K = 0
+                  END IF
+*
+*                 Check error code from ZHETRF and handle error.
+*
+                  IF( INFO.NE.K ) THEN
+                     CALL ALAERH( PATH, 'ZHETRF_AASEN', INFO, K, UPLO, 
+     $                            N, N, -1, -1, NB, IMAT, NFAIL, NERRS, 
+     $                            NOUT )
+                  END IF
+*
+*                 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_AASEN( UPLO, N, A, LDA, AFAC, LDA, IWORK,
+     $                               AINV, LDA, RWORK, RESULT( 1 ) )
+                  NT = 1
+*
+*
+*                 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
+*
+*                 Do only the condition estimate if INFO is not 0.
+*
+                  IF( TRFCON ) THEN
+                     RCONDC = ZERO
+                     GO TO 140
+                  END IF
+*
+*                 Do for each value of NRHS in NSVAL.
+*
+                  DO 130 IRHS = 1, NNS
+                     NRHS = NSVAL( IRHS )
+*
+*+    TEST 3 (Using TRS)
+*                 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_AASEN'
+                     LWORK = 3*N-2
+                     CALL ZHETRS_AASEN( UPLO, N, NRHS, AFAC, LDA, IWORK,
+     $                                  X, LDA, WORK, LWORK, INFO )
+*
+*                    Check error code from ZHETRS and handle error.
+*
+                     IF( INFO.NE.0 ) THEN
+                        CALL ALAERH( PATH, 'ZHETRS', INFO, 0, UPLO, N,
+     $                                N, -1, -1, NRHS, IMAT, NFAIL,
+     $                                NERRS, NOUT )
+                     END IF
+*
+                     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( 2 ) )
+*
+*                    Print information about the tests that did not pass
+*                    the threshold.
+*
+                     DO 120 K = 2, 2
+                        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
+  120                CONTINUE
+                     NRUN = NRUN + 1
+*
+*                 End do for each value of NRHS in NSVAL.
+*
+  130             CONTINUE
+  140             CONTINUE
+  150          CONTINUE
+  160       CONTINUE
+  170    CONTINUE
+  180 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 )
+c 9997 FORMAT( ' UPLO = ''', A1, ''', N =', I5, ',', 10X, ' type ', I2,
+c     $      ', test(', I2, ') =', G12.5 )
+ 9995 FORMAT( ' Invalid input value: ', A4, '=', I6, '; must be <=',
+     $      I6 )
+      RETURN
+*
+*     End of ZCHKHE_AASEN
+*
+      END
diff --git a/TESTING/LIN/zdrvhe_aasen.f b/TESTING/LIN/zdrvhe_aasen.f
new file mode 100644
index 00000000..0170873e
--- /dev/null
+++ b/TESTING/LIN/zdrvhe_aasen.f
@@ -0,0 +1,525 @@
+*> \brief \b ZDRVHE_AASEN
+*
+*  =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+*            http://www.netlib.org/lapack/explore-html/
+*
+*  Definition:
+*  ===========
+*
+*       SUBROUTINE ZDRVHE_AASEN( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, NMAX,
+*                                A, AFAC, AINV, B, X, XACT, WORK, RWORK, IWORK,
+*                                NOUT )
+*
+*       .. Scalar Arguments ..
+*       LOGICAL            TSTERR
+*       INTEGER            NMAX, NN, NOUT, NRHS
+*       DOUBLE PRECISION   THRESH
+*       ..
+*       .. Array Arguments ..
+*       LOGICAL            DOTYPE( * )
+*       INTEGER            IWORK( * ), NVAL( * )
+*       DOUBLE PRECISION   RWORK( * )
+*       COMPLEX*16         A( * ), AFAC( * ), AINV( * ), B( * ),
+*      $                   WORK( * ), X( * ), XACT( * )
+*       ..
+*
+*
+*> \par Purpose:
+*  =============
+*>
+*> \verbatim
+*>
+*> ZDRVHE_AASEN tests the driver routine ZHESV_AASEN.
+*> \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] NRHS
+*> \verbatim
+*>          NRHS is INTEGER
+*>          The number of right hand side vectors to be generated for
+*>          each linear system.
+*> \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*NRHS)
+*> \endverbatim
+*>
+*> \param[out] X
+*> \verbatim
+*>          X is COMPLEX*16 array, dimension (NMAX*NRHS)
+*> \endverbatim
+*>
+*> \param[out] XACT
+*> \verbatim
+*>          XACT is COMPLEX*16 array, dimension (NMAX*NRHS)
+*> \endverbatim
+*>
+*> \param[out] WORK
+*> \verbatim
+*>          WORK is COMPLEX*16 array, dimension (NMAX*max(2,NRHS))
+*> \endverbatim
+*>
+*> \param[out] RWORK
+*> \verbatim
+*>          RWORK is DOUBLE PRECISION array, dimension (NMAX+2*NRHS)
+*> \endverbatim
+*>
+*> \param[out] IWORK
+*> \verbatim
+*>          IWORK is INTEGER array, dimension (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 November 2016
+*
+*> \ingroup complex16_lin
+*
+*  =====================================================================
+      SUBROUTINE ZDRVHE_AASEN( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR,
+     $                         NMAX, A, AFAC, AINV, B, X, XACT, WORK, 
+     $                         RWORK, IWORK, NOUT )
+*
+*  -- LAPACK test routine (version 3.7.0) --
+*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*     November 2016
+*
+*     .. Scalar Arguments ..
+      LOGICAL            TSTERR
+      INTEGER            NMAX, NN, NOUT, NRHS
+      DOUBLE PRECISION   THRESH
+*     ..
+*     .. Array Arguments ..
+      LOGICAL            DOTYPE( * )
+      INTEGER            IWORK( * ), NVAL( * )
+      DOUBLE PRECISION   RWORK( * )
+      COMPLEX*16         A( * ), AFAC( * ), AINV( * ), B( * ),
+     $                   WORK( * ), X( * ), XACT( * )
+*     ..
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      DOUBLE PRECISION   ONE, ZERO
+      PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
+      INTEGER            NTYPES, NTESTS
+      PARAMETER          ( NTYPES = 10, NTESTS = 3 )
+      INTEGER            NFACT
+      PARAMETER          ( NFACT = 2 )
+*     ..
+*     .. Local Scalars ..
+      LOGICAL            ZEROT
+      CHARACTER          DIST, FACT, TYPE, UPLO, XTYPE
+      CHARACTER*3        MATPATH, PATH
+      INTEGER            I, I1, I2, IFACT, IMAT, IN, INFO, IOFF, IUPLO,
+     $                   IZERO, J, K, K1, KL, KU, LDA, LWORK, MODE, N,
+     $                   NB, NBMIN, NERRS, NFAIL, NIMAT, NRUN, NT
+      DOUBLE PRECISION   AINVNM, ANORM, CNDNUM, RCOND, RCONDC
+*     ..
+*     .. Local Arrays ..
+      CHARACTER          FACTS( NFACT ), UPLOS( 2 )
+      INTEGER            ISEED( 4 ), ISEEDY( 4 )
+      DOUBLE PRECISION   RESULT( NTESTS )
+*     ..
+*     .. External Functions ..
+      DOUBLE PRECISION   DGET06, ZLANHE
+      EXTERNAL           DGET06, ZLANHE
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           ALADHD, ALAERH, ALASVM, XLAENV, ZERRVX, ZGET04,
+     $                   ZHESV_AASEN, ZHET01_AASEN, ZHETRF_AASEN, 
+     $                   ZHETRI2, ZLACPY, ZLAIPD, ZLARHS, ZLATB4, ZLATMS,
+     $                   ZPOT02
+*     ..
+*     .. Scalars in Common ..
+      LOGICAL            LERR, OK
+      CHARACTER*32       SRNAMT
+      INTEGER            INFOT, NUNIT
+*     ..
+*     .. Common blocks ..
+      COMMON             / INFOC / INFOT, NUNIT, OK, LERR
+      COMMON             / SRNAMC / SRNAMT
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          DCMPLX, MAX, MIN
+*     ..
+*     .. Data statements ..
+      DATA               ISEEDY / 1988, 1989, 1990, 1991 /
+      DATA               UPLOS / 'U', 'L' / , FACTS / 'F', 'N' /
+*     ..
+*     .. Executable Statements ..
+*
+*     Initialize constants and the random number seed.
+*
+*     Test path
+*
+      PATH( 1: 1 ) = 'Zomplex precision'
+      PATH( 2: 3 ) = 'HA'
+*
+*     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
+      LWORK = MAX( 2*NMAX, NMAX*NRHS )
+*
+*     Test the error exits
+*
+      IF( TSTERR )
+     $   CALL ZERRVX( PATH, NOUT )
+      INFOT = 0
+*
+*     Set the block size and minimum block size for testing.
+*
+      NB = 1
+      NBMIN = 2
+      CALL XLAENV( 1, NB )
+      CALL XLAENV( 2, NBMIN )
+*
+*     Do for each value of N in NVAL
+*
+      DO 180 IN = 1, NN
+         N = NVAL( IN )
+         LDA = MAX( N, 1 )
+         XTYPE = 'N'
+         NIMAT = NTYPES
+         IF( N.LE.0 )
+     $      NIMAT = 1
+*
+         DO 170 IMAT = 1, NIMAT
+*
+*           Do the tests only if DOTYPE( IMAT ) is true.
+*
+            IF( .NOT.DOTYPE( IMAT ) )
+     $         GO TO 170
+*
+*           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 170
+*
+*           Do first for UPLO = 'U', then for UPLO = 'L'
+*
+            DO 160 IUPLO = 1, 2
+               UPLO = UPLOS( IUPLO )
+*
+*              Begin generate the test matrix A.
+*
+*              Set up parameters with ZLATB4 and generate a test matrix
+*              with ZLATMS.
+*
+               CALL ZLATB4( MATPATH, IMAT, N, N, TYPE, KL, KU, ANORM,
+     $                       MODE, CNDNUM, DIST )
+*
+               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.
+*
+               IF( INFO.NE.0 ) THEN
+                  CALL ALAERH( PATH, 'ZLATMS', INFO, 0, UPLO, N, N, -1,
+     $                         -1, -1, IMAT, NFAIL, NERRS, NOUT )
+                  GO TO 160
+               END IF
+*
+*              For 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 ) = ZERO
+   20                   CONTINUE
+                        IOFF = IOFF + IZERO
+                        DO 30 I = IZERO, N
+                           A( IOFF ) = ZERO
+                           IOFF = IOFF + LDA
+   30                   CONTINUE
+                     ELSE
+                        IOFF = IZERO
+                        DO 40 I = 1, IZERO - 1
+                           A( IOFF ) = ZERO
+                           IOFF = IOFF + LDA
+   40                   CONTINUE
+                        IOFF = IOFF - IZERO
+                        DO 50 I = IZERO, N
+                           A( IOFF+I ) = ZERO
+   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 ) = ZERO
+   60                      CONTINUE
+                           IOFF = IOFF + LDA
+   70                   CONTINUE
+                        IZERO = 1
+                     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 ) = ZERO
+   80                      CONTINUE
+                           IOFF = IOFF + LDA
+   90                   CONTINUE
+                     END IF
+                  END IF
+               ELSE
+                  IZERO = 0
+               END IF
+*
+*              Set the imaginary part of the diagonals.
+*
+               CALL ZLAIPD( N, A, LDA+1, 0 )
+*
+               DO 150 IFACT = 1, NFACT
+*
+*                 Do first for FACT = 'F', then for other values.
+*
+                  FACT = FACTS( IFACT )
+*
+*                 Compute the condition number for comparison with
+*                 the value returned by ZHESVX.
+*
+                  IF( ZEROT ) THEN
+                     IF( IFACT.EQ.1 )
+     $                  GO TO 150
+                     RCONDC = ZERO
+*
+                  ELSE IF( IFACT.EQ.1 ) THEN
+*
+*                    Compute the 1-norm of A.
+*
+                     ANORM = ZLANHE( '1', UPLO, N, A, LDA, RWORK )
+*
+*                    Factor the matrix A.
+*
+c                     CALL ZLACPY( UPLO, N, N, A, LDA, AFAC, LDA )
+c                     CALL ZHETRF( UPLO, N, AFAC, LDA, IWORK, WORK,
+c     $                            LWORK, INFO )
+*
+*                    Compute inv(A) and take its norm.
+*
+c                     CALL ZLACPY( UPLO, N, N, AFAC, LDA, AINV, LDA )
+c                     LWORK = (N+NB+1)*(NB+3)
+c                     CALL ZHETRI2( UPLO, N, AINV, LDA, IWORK, WORK,
+c     $                            LWORK, INFO )
+c                     AINVNM = ZLANHE( '1', UPLO, N, AINV, LDA, RWORK )
+*
+*                    Compute the 1-norm condition number of A.
+*
+c                     IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN
+c                        RCONDC = ONE
+c                     ELSE
+c                        RCONDC = ( ONE / ANORM ) / AINVNM
+c                     END IF
+                  END IF
+*
+*                 Form an exact solution and set the right hand side.
+*
+                  SRNAMT = 'ZLARHS'
+                  CALL ZLARHS( MATPATH, XTYPE, UPLO, ' ', N, N, KL, KU,
+     $                         NRHS, A, LDA, XACT, LDA, B, LDA, ISEED,
+     $                         INFO )
+                  XTYPE = 'C'
+*
+*                 --- Test ZHESV_AASEN  ---
+*
+                  IF( IFACT.EQ.2 ) THEN
+                     CALL ZLACPY( UPLO, N, N, A, LDA, AFAC, LDA )
+                     CALL ZLACPY( 'Full', N, NRHS, B, LDA, X, LDA )
+*
+*                    Factor the matrix and solve the system using ZHESV.
+*
+                     SRNAMT = 'ZHESV_AASEN '
+                     CALL ZHESV_AASEN( UPLO, N, NRHS, AFAC, LDA, IWORK,
+     $                                 X, LDA, WORK, LWORK, INFO )
+*
+*                    Adjust the expected value of INFO to account for
+*                    pivoting.
+*
+                     IF( IZERO.GT.0 ) THEN
+                        J = 1
+                        K = IZERO
+  100                   CONTINUE
+                        IF( J.EQ.K ) THEN
+                           K = IWORK( J )
+                        ELSE IF( IWORK( J ).EQ.K ) THEN
+                           K = J
+                        END IF
+                        IF( J.LT.K ) THEN
+                           J = J + 1
+                           GO TO 100
+                        END IF
+                     ELSE
+                        K = 0
+                     END IF
+*
+*                    Check error code from ZHESV .
+*
+                     IF( INFO.NE.K ) THEN
+                        CALL ALAERH( PATH, 'ZHESV_AASEN', INFO, K, UPLO, N,
+     $                               N, -1, -1, NRHS, IMAT, NFAIL,
+     $                               NERRS, NOUT )
+                        GO TO 120
+                     ELSE IF( INFO.NE.0 ) THEN
+                        GO TO 120
+                     END IF
+*
+*                    Reconstruct matrix from factors and compute
+*                    residual.
+*
+                     CALL ZHET01_AASEN( UPLO, N, A, LDA, AFAC, LDA,
+     $                                  IWORK, AINV, LDA, RWORK, 
+     $                                  RESULT( 1 ) )
+*
+*                    Compute residual of the computed solution.
+*
+                     CALL ZLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA )
+                     CALL ZPOT02( UPLO, N, NRHS, A, LDA, X, LDA, WORK,
+     $                            LDA, RWORK, RESULT( 2 ) )
+*
+*                    Check solution from generated exact solution.
+*
+                     CALL ZGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
+     $                            RESULT( 3 ) )
+                     NT = 3
+*
+*                    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 ALADHD( NOUT, PATH )
+                           WRITE( NOUT, FMT = 9999 )'ZHESV_AASEN', UPLO, N,
+     $                        IMAT, K, RESULT( K )
+                           NFAIL = NFAIL + 1
+                        END IF
+  110                CONTINUE
+                     NRUN = NRUN + NT
+  120                CONTINUE
+                  END IF
+*
+  150          CONTINUE
+*
+  160       CONTINUE
+  170    CONTINUE
+  180 CONTINUE
+*
+*     Print a summary of the results.
+*
+      CALL ALASVM( PATH, NOUT, NFAIL, NRUN, NERRS )
+*
+ 9999 FORMAT( 1X, A, ', UPLO=''', A1, ''', N =', I5, ', type ', I2,
+     $      ', test ', I2, ', ratio =', G12.5 )
+       RETURN
+*
+*     End of ZDRVHE_AASEN
+*
+      END
diff --git a/TESTING/LIN/zerrvx.f b/TESTING/LIN/zerrvx.f
index 11a35aa4..1bb2d222 100644
--- a/TESTING/LIN/zerrvx.f
+++ b/TESTING/LIN/zerrvx.f
@@ -93,7 +93,7 @@
      $                   ZGTSVX, ZHESV, ZHESV_ROOK, ZHESVX, ZHPSV,
      $                   ZHPSVX, ZPBSV, ZPBSVX, ZPOSV, ZPOSVX, ZPPSV,
      $                   ZPPSVX, ZPTSV, ZPTSVX, ZSPSV, ZSPSVX, ZSYSV,
-     $                   ZSYSV_ROOK, ZSYSVX
+     $                   ZSYSV_AASEN, ZSYSV_ROOK, ZSYSVX
 *     ..
 *     .. Scalars in Common ..
       LOGICAL            LERR, OK
@@ -634,6 +634,25 @@
      $                RCOND, R1, R2, W, 3, RW, INFO )
          CALL CHKXER( 'ZHESVX', INFOT, NOUT, LERR, OK )
 *
+      ELSE IF( LSAMEN( 2, C2, 'HA' ) ) THEN
+*
+*        ZHESV_AASEN
+*
+        SRNAMT = 'ZHESV_AASEN'
+        INFOT = 1
+        CALL ZHESV_AASEN( '/', 0, 0, A, 1, IP, B, 1, W, 1, INFO )
+        CALL CHKXER( 'ZHESV_AASEN', INFOT, NOUT, LERR, OK )
+        INFOT = 2
+        CALL ZHESV_AASEN( 'U', -1, 0, A, 1, IP, B, 1, W, 1, INFO )
+        CALL CHKXER( 'ZHESV_AASEN', INFOT, NOUT, LERR, OK )
+        INFOT = 3
+        CALL ZHESV_AASEN( 'U', 0, -1, A, 1, IP, B, 1, W, 1, INFO )
+        CALL CHKXER( 'ZHESV_AASEN', INFOT, NOUT, LERR, OK )
+        INFOT = 8
+        CALL ZHESV_AASEN( 'U', 2, 0, A, 2, IP, B, 1, W, 1, INFO )
+        CALL CHKXER( 'ZHESV_AASEN', INFOT, NOUT, LERR, OK )
+*
+
       ELSE IF( LSAMEN( 2, C2, 'HR' ) ) THEN
 *
 *        ZHESV_ROOK
diff --git a/TESTING/LIN/zhet01_aasen.f b/TESTING/LIN/zhet01_aasen.f
new file mode 100644
index 00000000..89b87990
--- /dev/null
+++ b/TESTING/LIN/zhet01_aasen.f
@@ -0,0 +1,267 @@
+*> \brief \b ZHET01_AASEN
+*
+*  =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+*            http://www.netlib.org/lapack/explore-html/
+*
+*  Definition:
+*  ===========
+*
+*       SUBROUTINE ZHET01_AASEN( UPLO, N, A, LDA, AFAC, LDAFAC, IPIV,
+*                          C, LDC, RWORK, RESID )
+*
+*       .. Scalar Arguments ..
+*       CHARACTER          UPLO
+*       INTEGER            LDA, LDAFAC, LDC, N
+*       COMPLEX*16   RESID
+*       ..
+*       .. Array Arguments ..
+*       INTEGER            IPIV( * )
+*       COMPLEX*16   A( LDA, * ), AFAC( LDAFAC, * ), C( LDC, * ),
+*      $                   RWORK( * )
+*       ..
+*
+*
+*> \par Purpose:
+*  =============
+*>
+*> \verbatim
+*>
+*> ZHET01_AASEN reconstructs a 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 and EPS is the machine epsilon.
+*> \endverbatim
+*
+*  Arguments:
+*  ==========
+*
+*> \param[in] UPLO
+*> \verbatim
+*>          UPLO is CHARACTER*1
+*>          Specifies whether the upper or lower triangular part of the
+*>          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 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 ZHETRF.
+*> \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 ZHETRF.
+*> \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 COMPLEX*16 array, dimension (N)
+*> \endverbatim
+*>
+*> \param[out] RESID
+*> \verbatim
+*>          RESID is COMPLEX*16
+*>          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 November 2016
+*
+*
+*> \ingroup complex16_lin
+*
+*  =====================================================================
+      SUBROUTINE ZHET01_AASEN( UPLO, N, A, LDA, AFAC, LDAFAC, IPIV, C,
+     $                         LDC, RWORK, RESID )
+*
+*  -- LAPACK test routine (version 3.7.0) --
+*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*     November 2016
+*
+*     .. Scalar Arguments ..
+      CHARACTER          UPLO
+      INTEGER            LDA, LDAFAC, LDC, N
+      DOUBLE PRECISION   RESID
+*     ..
+*     .. Array Arguments ..
+      INTEGER            IPIV( * )
+      COMPLEX*16         A( LDA, * ), AFAC( LDAFAC, * ), C( LDC, * ),
+     $                   RWORK( * )
+*     ..
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      COMPLEX*16         CZERO, CONE
+      PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ), 
+     $                     CONE  = ( 1.0D+0, 0.0D+0 ) )
+      DOUBLE PRECISION   ZERO, ONE
+      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
+*     ..
+*     .. Local Scalars ..
+      INTEGER            I, J
+      DOUBLE PRECISION   ANORM, EPS
+*     ..
+*     .. External Functions ..
+      LOGICAL            LSAME
+      DOUBLE PRECISION   DLAMCH, ZLANHE
+      EXTERNAL           LSAME, DLAMCH, ZLANHE
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           ZLASET, ZLAVHE
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          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 )
+*
+*     Initialize C to the tridiagonal matrix T.
+*
+      CALL ZLASET( 'Full', N, N, CZERO, CZERO, C, LDC )
+      CALL ZLACPY( 'F', 1, N, AFAC( 1, 1 ), LDAFAC+1, C( 1, 1 ), LDC+1 )
+      IF( N.GT.1 ) THEN
+         IF( LSAME( UPLO, 'U' ) ) THEN
+            CALL ZLACPY( 'F', 1, N-1, AFAC( 1, 2 ), LDAFAC+1, C( 1, 2 ),
+     $                   LDC+1 )
+            CALL ZLACPY( 'F', 1, N-1, AFAC( 1, 2 ), LDAFAC+1, C( 2, 1 ),
+     $                   LDC+1 )
+            CALL ZLACGV( N-1, C( 2, 1 ), LDC+1 )
+         ELSE
+            CALL ZLACPY( 'F', 1, N-1, AFAC( 2, 1 ), LDAFAC+1, C( 1, 2 ),
+     $                   LDC+1 )
+            CALL ZLACPY( 'F', 1, N-1, AFAC( 2, 1 ), LDAFAC+1, C( 2, 1 ),
+     $                   LDC+1 )
+            CALL ZLACGV( N-1, C( 1, 2 ), LDC+1 )
+         ENDIF
+      ENDIF
+*
+*     Call ZTRMM to form the product U' * D (or L * D ).
+*
+      IF( LSAME( UPLO, 'U' ) ) THEN
+         CALL ZTRMM( 'Left', UPLO, 'Conjugate transpose', 'Unit', N-1,
+     $               N, CONE, AFAC( 1, 2 ), LDAFAC, C( 2, 1 ), LDC )
+      ELSE
+         CALL ZTRMM( 'Left', UPLO, 'No transpose', 'Unit', N-1, N,
+     $               CONE, AFAC( 2, 1 ), LDAFAC, C( 2, 1 ), LDC )
+      END IF
+*
+*     Call ZTRMM again to multiply by U (or L ).
+*
+      IF( LSAME( UPLO, 'U' ) ) THEN
+         CALL ZTRMM( 'Right', UPLO, 'No transpose', 'Unit', N, N-1,
+     $               CONE, AFAC( 1, 2 ), LDAFAC, C( 1, 2 ), LDC )
+      ELSE
+         CALL ZTRMM( 'Right', UPLO, 'Conjugate transpose', 'Unit', N,
+     $               N-1, CONE, AFAC( 2, 1 ), LDAFAC, C( 1, 2 ), LDC )
+      END IF
+*
+*     Apply hermitian pivots
+*
+      DO J = N, 1, -1
+         I = IPIV( J )
+         IF( I.NE.J )
+     $      CALL ZSWAP( N, C( J, 1 ), LDC, C( I, 1 ), LDC )
+      END DO
+      DO J = N, 1, -1
+         I = IPIV( J )
+         IF( I.NE.J )
+     $      CALL ZSWAP( N, C( 1, J ), 1, C( 1, I ), 1 )
+      END DO
+*
+*
+*     Compute the difference  C - A .
+*
+      IF( LSAME( UPLO, 'U' ) ) THEN
+         DO J = 1, N
+            DO I = 1, J
+               C( I, J ) = C( I, J ) - A( I, J )
+            END DO
+         END DO
+      ELSE
+         DO J = 1, N
+            DO I = J, N
+               C( I, J ) = C( I, J ) - A( I, J )
+            END DO
+         END DO
+      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_AASEN
+*
+      END