--- rpl/lapack/lapack/dsytrf_aa_2stage.f	2018/05/29 14:49:25	1.1
+++ rpl/lapack/lapack/dsytrf_aa_2stage.f	2020/05/21 21:46:02	1.2
@@ -38,7 +38,7 @@
 *> DSYTRF_AA_2STAGE computes the factorization of a real symmetric matrix A
 *> using the Aasen's algorithm.  The form of the factorization is
 *>
-*>    A = U*T*U**T  or  A = L*T*L**T
+*>    A = U**T*T*U  or  A = L*T*L**T
 *>
 *> where U (or L) is a product of permutation and unit upper (lower)
 *> triangular matrices, and T is a symmetric band matrix with the
@@ -93,6 +93,7 @@
 *>
 *> \param[in] LTB
 *> \verbatim
+*>          LTB is INTEGER
 *>          The size of the array TB. LTB >= 4*N, internally
 *>          used to select NB such that LTB >= (3*NB+1)*N.
 *>
@@ -102,6 +103,22 @@
 *>          no error message related to LTB is issued by XERBLA.
 *> \endverbatim
 *>
+*> \param[out] IPIV
+*> \verbatim
+*>          IPIV is INTEGER array, dimension (N)
+*>          On exit, it contains the details of the interchanges, i.e.,
+*>          the row and column k of A were interchanged with the
+*>          row and column IPIV(k).
+*> \endverbatim
+*>
+*> \param[out] IPIV2
+*> \verbatim
+*>          IPIV2 is INTEGER array, dimension (N)
+*>          On exit, it contains the details of the interchanges, i.e.,
+*>          the row and column k of T were interchanged with the
+*>          row and column IPIV2(k).
+*> \endverbatim
+*>
 *> \param[out] WORK
 *> \verbatim
 *>          WORK is DOUBLE PRECISION workspace of size LWORK
@@ -109,6 +126,7 @@
 *>
 *> \param[in] LWORK
 *> \verbatim
+*>          LWORK is INTEGER
 *>          The size of WORK. LWORK >= N, internally used to select NB
 *>          such that LWORK >= N*NB.
 *>
@@ -118,22 +136,6 @@
 *>          no error message related to LWORK is issued by XERBLA.
 *> \endverbatim
 *>
-*> \param[out] IPIV
-*> \verbatim
-*>          IPIV is INTEGER array, dimension (N)
-*>          On exit, it contains the details of the interchanges, i.e.,
-*>          the row and column k of A were interchanged with the
-*>          row and column IPIV(k).
-*> \endverbatim
-*>
-*> \param[out] IPIV2
-*> \verbatim
-*>          IPIV is INTEGER array, dimension (N)
-*>          On exit, it contains the details of the interchanges, i.e.,
-*>          the row and column k of T were interchanged with the
-*>          row and column IPIV(k).
-*> \endverbatim
-*>
 *> \param[out] INFO
 *> \verbatim
 *>          INFO is INTEGER
@@ -191,7 +193,7 @@
       EXTERNAL           LSAME, ILAENV
 *     ..
 *     .. External Subroutines ..
-      EXTERNAL           XERBLA, DCOPY, DLACGV, DLACPY,
+      EXTERNAL           XERBLA, DCOPY, DLACPY,
      $                   DLASET, DGBTRF, DGEMM,  DGETRF, 
      $                   DSYGST, DSWAP, DTRSM 
 *     ..
@@ -273,7 +275,7 @@
       IF( UPPER ) THEN
 *
 *        .....................................................
-*        Factorize A as L*D*L**T using the upper triangle of A
+*        Factorize A as U**T*D*U using the upper triangle of A
 *        .....................................................
 *
          DO J = 0, NT-1
@@ -440,12 +442,14 @@ c               END IF
 *                    > Apply pivots to previous columns of L
                      CALL DSWAP( K-1, A( (J+1)*NB+1, I1 ), 1, 
      $                                A( (J+1)*NB+1, I2 ), 1 )
-*                    > Swap A(I1+1:M, I1) with A(I2, I1+1:M)               
-                     CALL DSWAP( I2-I1-1, A( I1, I1+1 ), LDA,
-     $                                    A( I1+1, I2 ), 1 )
+*                    > Swap A(I1+1:M, I1) with A(I2, I1+1:M)
+                     IF( I2.GT.(I1+1) )
+     $                  CALL DSWAP( I2-I1-1, A( I1, I1+1 ), LDA,
+     $                                       A( I1+1, I2 ), 1 )
 *                    > Swap A(I2+1:M, I1) with A(I2+1:M, I2)
-                     CALL DSWAP( N-I2, A( I1, I2+1 ), LDA,
-     $                                 A( I2, I2+1 ), LDA ) 
+                     IF( I2.LT.N )
+     $                  CALL DSWAP( N-I2, A( I1, I2+1 ), LDA,
+     $                                    A( I2, I2+1 ), LDA ) 
 *                    > Swap A(I1, I1) with A(I2, I2)
                      PIV = A( I1, I1 )
                      A( I1, I1 ) = A( I2, I2 )
@@ -614,11 +618,13 @@ c               END IF
                      CALL DSWAP( K-1, A( I1, (J+1)*NB+1 ), LDA, 
      $                                A( I2, (J+1)*NB+1 ), LDA )
 *                    > Swap A(I1+1:M, I1) with A(I2, I1+1:M)               
-                     CALL DSWAP( I2-I1-1, A( I1+1, I1 ), 1,
-     $                                    A( I2, I1+1 ), LDA )
+                     IF( I2.GT.(I1+1) )
+     $                  CALL DSWAP( I2-I1-1, A( I1+1, I1 ), 1,
+     $                                       A( I2, I1+1 ), LDA )
 *                    > Swap A(I2+1:M, I1) with A(I2+1:M, I2)
-                     CALL DSWAP( N-I2, A( I2+1, I1 ), 1,
-     $                                 A( I2+1, I2 ), 1 ) 
+                     IF( I2.LT.N )
+     $                  CALL DSWAP( N-I2, A( I2+1, I1 ), 1,
+     $                                    A( I2+1, I2 ), 1 ) 
 *                    > Swap A(I1, I1) with A(I2, I2)
                      PIV = A( I1, I1 )
                      A( I1, I1 ) = A( I2, I2 )
@@ -642,6 +648,8 @@ c     $                     (J+1)*NB+1,
 *     Factor the band matrix
       CALL DGBTRF( N, N, NB, NB, TB, LDTB, IPIV2, INFO )
 *
+      RETURN
+*
 *     End of DSYTRF_AA_2STAGE
 *
       END