--- rpl/lapack/lapack/zhetrs_aa.f	2018/05/29 07:18:21	1.4
+++ rpl/lapack/lapack/zhetrs_aa.f	2020/05/21 21:46:06	1.5
@@ -38,8 +38,8 @@
 *> \verbatim
 *>
 *> ZHETRS_AA solves a system of linear equations A*X = B with a complex
-*> hermitian matrix A using the factorization A = U*T*U**H or
-*> A = L*T*L**T computed by ZHETRF_AA.
+*> hermitian matrix A using the factorization A = U**H*T*U or
+*> A = L*T*L**H computed by ZHETRF_AA.
 *> \endverbatim
 *
 *  Arguments:
@@ -50,7 +50,7 @@
 *>          UPLO is CHARACTER*1
 *>          Specifies whether the details of the factorization are stored
 *>          as an upper or lower triangular matrix.
-*>          = 'U':  Upper triangular, form is A = U*T*U**H;
+*>          = 'U':  Upper triangular, form is A = U**H*T*U;
 *>          = 'L':  Lower triangular, form is A = L*T*L**H.
 *> \endverbatim
 *>
@@ -98,14 +98,16 @@
 *>          The leading dimension of the array B.  LDB >= max(1,N).
 *> \endverbatim
 *>
-*> \param[in] WORK
+*> \param[out] WORK
 *> \verbatim
-*>          WORK is DOUBLE array, dimension (MAX(1,LWORK))
+*>          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
 *> \endverbatim
 *>
 *> \param[in] LWORK
 *> \verbatim
-*>          LWORK is INTEGER, LWORK >= MAX(1,3*N-2).
+*>          LWORK is INTEGER
+*>          The dimension of the array WORK. LWORK >= max(1,3*N-2).
+*> \endverbatim
 *>
 *> \param[out] INFO
 *> \verbatim
@@ -199,61 +201,80 @@
 *
       IF( UPPER ) THEN
 *
-*        Solve A*X = B, where A = U*T*U**T.
+*        Solve A*X = B, where A = U**H*T*U.
 *
-*        Pivot, P**T * B
+*        1) Forward substitution with U**H
 *
-         DO K = 1, N
-            KP = IPIV( K )
-            IF( KP.NE.K )
-     $          CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
-         END DO
+         IF( N.GT.1 ) THEN
 *
-*        Compute (U \P**T * B) -> B    [ (U \P**T * B) ]
+*           Pivot, P**T * B -> B
 *
-         CALL ZTRSM('L', 'U', 'C', 'U', N-1, NRHS, ONE, A( 1, 2 ), LDA,
-     $               B( 2, 1 ), LDB)
+            DO K = 1, N
+               KP = IPIV( K )
+               IF( KP.NE.K )
+     $            CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+            END DO
 *
-*        Compute T \ B -> B   [ T \ (U \P**T * B) ]
+*           Compute U**H \ B -> B    [ (U**H \P**T * B) ]
 *
-         CALL ZLACPY( 'F', 1, N, A(1, 1), LDA+1, WORK(N), 1)
+            CALL ZTRSM( 'L', 'U', 'C', 'U', N-1, NRHS, ONE, A( 1, 2 ),
+     $                  LDA, B( 2, 1 ), LDB )
+         END IF
+*
+*        2) Solve with triangular matrix T
+*
+*        Compute T \ B -> B   [ T \ (U**H \P**T * B) ]
+*
+         CALL ZLACPY( 'F', 1, N, A(1, 1), LDA+1, WORK(N), 1 )
          IF( N.GT.1 ) THEN
              CALL ZLACPY( 'F', 1, N-1, A( 1, 2 ), LDA+1, WORK( 2*N ), 1)
-             CALL ZLACPY( 'F', 1, N-1, A( 1, 2 ), LDA+1, WORK( 1 ), 1)
+             CALL ZLACPY( 'F', 1, N-1, A( 1, 2 ), LDA+1, WORK( 1 ), 1 )
              CALL ZLACGV( N-1, WORK( 1 ), 1 )
          END IF
-         CALL ZGTSV(N, NRHS, WORK(1), WORK(N), WORK(2*N), B, LDB,
-     $              INFO)
+         CALL ZGTSV( N, NRHS, WORK(1), WORK(N), WORK(2*N), B, LDB,
+     $               INFO )
+*
+*        3) Backward substitution with U
+*
+         IF( N.GT.1 ) THEN
 *
-*        Compute (U**T \ B) -> B   [ U**T \ (T \ (U \P**T * B) ) ]
+*           Compute U \ B -> B   [ U \ (T \ (U**H \P**T * B) ) ]
 *
-         CALL ZTRSM( 'L', 'U', 'N', 'U', N-1, NRHS, ONE, A( 1, 2 ), LDA,
-     $               B(2, 1), LDB)
+            CALL ZTRSM( 'L', 'U', 'N', 'U', N-1, NRHS, ONE, A( 1, 2 ),
+     $                  LDA, B(2, 1), LDB)
 *
-*        Pivot, P * B  [ P * (U**T \ (T \ (U \P**T * B) )) ]
+*           Pivot, P * B  [ P * (U**H \ (T \ (U \P**T * B) )) ]
 *
-         DO K = N, 1, -1
-            KP = IPIV( K )
-            IF( KP.NE.K )
-     $         CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
-         END DO
+            DO K = N, 1, -1
+               KP = IPIV( K )
+               IF( KP.NE.K )
+     $            CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+            END DO
+         END IF
 *
       ELSE
 *
-*        Solve A*X = B, where A = L*T*L**T.
+*        Solve A*X = B, where A = L*T*L**H.
+*
+*        1) Forward substitution with L
 *
-*        Pivot, P**T * B
+         IF( N.GT.1 ) THEN
+*
+*           Pivot, P**T * B -> B
 *
-         DO K = 1, N
-            KP = IPIV( K )
-            IF( KP.NE.K )
-     $         CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
-         END DO
+            DO K = 1, N
+               KP = IPIV( K )
+               IF( KP.NE.K )
+     $            CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+            END DO
 *
-*        Compute (L \P**T * B) -> B    [ (L \P**T * B) ]
+*           Compute L \ B -> B    [ (L \P**T * B) ]
 *
-         CALL ZTRSM( 'L', 'L', 'N', 'U', N-1, NRHS, ONE, A( 2, 1 ), LDA,
-     $               B(2, 1), LDB)
+            CALL ZTRSM( 'L', 'L', 'N', 'U', N-1, NRHS, ONE, A( 2, 1 ),
+     $                  LDA, B(2, 1), LDB)
+         END IF
+*
+*        2) Solve with triangular matrix T
 *
 *        Compute T \ B -> B   [ T \ (L \P**T * B) ]
 *
@@ -266,18 +287,23 @@
          CALL ZGTSV(N, NRHS, WORK(1), WORK(N), WORK(2*N), B, LDB,
      $              INFO)
 *
-*        Compute (L**T \ B) -> B   [ L**T \ (T \ (L \P**T * B) ) ]
+*        3) Backward substitution with L**H
+*
+         IF( N.GT.1 ) THEN
 *
-         CALL ZTRSM( 'L', 'L', 'C', 'U', N-1, NRHS, ONE, A( 2, 1 ), LDA,
-     $              B( 2, 1 ), LDB)
+*           Compute L**H \ B -> B   [ L**H \ (T \ (L \P**T * B) ) ]
 *
-*        Pivot, P * B  [ P * (L**T \ (T \ (L \P**T * B) )) ]
+            CALL ZTRSM( 'L', 'L', 'C', 'U', N-1, NRHS, ONE, A( 2, 1 ),
+     $                  LDA, B( 2, 1 ), LDB)
 *
-         DO K = N, 1, -1
-            KP = IPIV( K )
-            IF( KP.NE.K )
-     $         CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
-         END DO
+*           Pivot, P * B  [ P * (L**H \ (T \ (L \P**T * B) )) ]
+*
+            DO K = N, 1, -1
+               KP = IPIV( K )
+               IF( KP.NE.K )
+     $            CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+            END DO
+         END IF
 *
       END IF
 *