Merge pull request #314 from iyamazaki/aasen

fix corner-case in Aasen (and typos in documentations)

Author: langou

SRC/chesv_aa.f

@@ -42,7 +42,7 @@
 *> matrices.
 *>
 *> Aasen's algorithm is used to factor A as
-*>    A = U * T * U**H,  if UPLO = 'U', or
+*>    A = U**H * T * U,  if UPLO = 'U', or
 *>    A = L * T * L**H,  if UPLO = 'L',
 *> where U (or L) is a product of permutation and unit upper (lower)
 *> triangular matrices, and T is Hermitian and tridiagonal. The factored form
@@ -86,7 +86,7 @@
 *>
 *>          On exit, if INFO = 0, the tridiagonal matrix T and the
 *>          multipliers used to obtain the factor U or L from the
-*>          factorization A = U*T*U**H or A = L*T*L**H as computed by
+*>          factorization A = U**H*T*U or A = L*T*L**H as computed by
 *>          CHETRF_AA.
 *> \endverbatim
 *>
@@ -230,7 +230,7 @@ SUBROUTINE CHESV_AA( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, WORK,
          RETURN
       END IF
 *
-*     Compute the factorization A = U*T*U**H or A = L*T*L**H.
+*     Compute the factorization A = U**H*T*U or A = L*T*L**H.
 *
       CALL CHETRF_AA( UPLO, N, A, LDA, IPIV, WORK, LWORK, INFO )
       IF( INFO.EQ.0 ) THEN

SRC/chesv_aa_2stage.f

@@ -43,7 +43,7 @@
 *> matrices.
 *>
 *> Aasen's 2-stage algorithm is used to factor A as
-*>    A = U * T * U**H,  if UPLO = 'U', or
+*>    A = U**H * T * U,  if UPLO = 'U', or
 *>    A = L * T * L**H,  if UPLO = 'L',
 *> where U (or L) is a product of permutation and unit upper (lower)
 *> triangular matrices, and T is Hermitian and band. The matrix T is
@@ -257,7 +257,7 @@ SUBROUTINE CHESV_AA_2STAGE( UPLO, N, NRHS, A, LDA, TB, LTB,
       END IF
 *
 *
-*     Compute the factorization A = U*T*U**H or A = L*T*L**H.
+*     Compute the factorization A = U**H*T*U or A = L*T*L**H.
 *
       CALL CHETRF_AA_2STAGE( UPLO, N, A, LDA, TB, LTB, IPIV, IPIV2,
      $                       WORK, LWORK, INFO )

SRC/chetrf_aa.f

@@ -37,7 +37,7 @@
 *> CHETRF_AA computes the factorization of a complex hermitian matrix A
 *> using the Aasen's algorithm.  The form of the factorization is
 *>
-*>    A = U*T*U**H  or  A = L*T*L**H
+*>    A = U**H*T*U  or  A = L*T*L**H
 *>
 *> where U (or L) is a product of permutation and unit upper (lower)
 *> triangular matrices, and T is a hermitian tridiagonal matrix.
@@ -223,7 +223,7 @@ SUBROUTINE CHETRF_AA( UPLO, N, A, LDA, IPIV, WORK, LWORK, INFO)
       IF( UPPER ) THEN
 *
 *        .....................................................
-*        Factorize A as L*D*L**H using the upper triangle of A
+*        Factorize A as U**H*D*U using the upper triangle of A
 *        .....................................................
 *
 *        copy first row A(1, 1:N) into H(1:n) (stored in WORK(1:N))

SRC/chetrf_aa_2stage.f

@@ -38,7 +38,7 @@
 *> CHETRF_AA_2STAGE computes the factorization of a real hermitian 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 hermitian band matrix with the
@@ -277,7 +277,7 @@ SUBROUTINE CHETRF_AA_2STAGE( UPLO, N, A, LDA, TB, LTB, IPIV,
       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

SRC/chetrs_aa.f

@@ -37,7 +37,7 @@
 *> \verbatim
 *>
 *> CHETRS_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
+*> hermitian matrix A using the factorization A = U**H*T*U or
 *> A = L*T*L**H computed by CHETRF_AA.
 *> \endverbatim
 *
@@ -49,7 +49,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
 *>
@@ -200,24 +200,31 @@ SUBROUTINE CHETRS_AA( UPLO, N, NRHS, A, LDA, IPIV, B, LDB,
 *
       IF( UPPER ) THEN
 *
-*        Solve A*X = B, where A = U*T*U**T.
+*        Solve A*X = B, where A = U**H*T*U.
 *
-*        P**T * B
+*        1) Forward substitution with U**H
 *
-         K = 1
-         DO WHILE ( K.LE.N )
-            KP = IPIV( K )
-            IF( KP.NE.K )
-     $          CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
-            K = K + 1
-         END DO
+         IF( N.GT.1 ) THEN
+*
+*           Pivot, P**T * B -> B
+*
+            K = 1
+            DO WHILE ( K.LE.N )
+               KP = IPIV( K )
+               IF( KP.NE.K )
+     $            CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+               K = K + 1
+            END DO
 *
-*        Compute (U \P**T * B) -> B    [ (U \P**T * B) ]
+*           Compute U**H \ B -> B    [ (U**H \P**T * B) ]
+*
+            CALL CTRSM( 'L', 'U', 'C', 'U', N-1, NRHS, ONE, A( 1, 2 ),
+     $                  LDA, B( 2, 1 ), LDB)
+         END IF
 *
-         CALL CTRSM('L', 'U', 'C', 'U', N-1, NRHS, ONE, A( 1, 2 ), LDA,
-     $               B( 2, 1 ), LDB)
+*        2) Solve with triangular matrix T
 *
-*        Compute T \ B -> B   [ T \ (U \P**T * B) ]
+*        Compute T \ B -> B   [ T \ (U**H \P**T * B) ]
 *
          CALL CLACPY( 'F', 1, N, A(1, 1), LDA+1, WORK(N), 1)
          IF( N.GT.1 ) THEN
@@ -228,65 +235,82 @@ SUBROUTINE CHETRS_AA( UPLO, N, NRHS, A, LDA, IPIV, B, LDB,
          CALL CGTSV(N, NRHS, WORK(1), WORK(N), WORK(2*N), B, LDB,
      $              INFO)
 *
-*        Compute (U**T \ B) -> B   [ U**T \ (T \ (U \P**T * B) ) ]
+*        3) Backward substitution with U
 *
-         CALL CTRSM( 'L', 'U', 'N', 'U', N-1, NRHS, ONE, A( 1, 2 ), LDA,
-     $               B(2, 1), LDB)
+         IF( N.GT.1 ) THEN
+*
+*           Compute U \ B -> B   [ U \ (T \ (U**H \P**T * B) ) ]
 *
-*        Pivot, P * B  [ P * (U**T \ (T \ (U \P**T * B) )) ]
+            CALL CTRSM( 'L', 'U', 'N', 'U', N-1, NRHS, ONE, A( 1, 2 ),
+     $                  LDA, B(2, 1), LDB)
 *
-         K = N
-         DO WHILE ( K.GE.1 )
-            KP = IPIV( K )
-            IF( KP.NE.K )
-     $         CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
-            K = K - 1
-         END DO
+*           Pivot, P * B  -> B [ P * (U \ (T \ (U**H \P**T * B) )) ]
+*
+            K = N
+            DO WHILE ( K.GE.1 )
+               KP = IPIV( K )
+               IF( KP.NE.K )
+     $            CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+               K = K - 1
+            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.
 *
-*        Pivot, P**T * B
+*        1) Forward substitution with L
 *
-         K = 1
-         DO WHILE ( K.LE.N )
-            KP = IPIV( K )
-            IF( KP.NE.K )
-     $         CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
-            K = K + 1
-         END DO
+         IF( N.GT.1 ) THEN
+*
+*           Pivot, P**T * B -> B
+*
+            K = 1
+            DO WHILE ( K.LE.N )
+               KP = IPIV( K )
+               IF( KP.NE.K )
+     $            CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+               K = K + 1
+            END DO
 *
-*        Compute (L \P**T * B) -> B    [ (L \P**T * B) ]
+*           Compute L \ B -> B    [ (L \P**T * B) ]
+*
+            CALL CTRSM( 'L', 'L', 'N', 'U', N-1, NRHS, ONE, A( 2, 1),
+     $                  LDA, B(2, 1), LDB )
+         END IF
 *
-         CALL CTRSM( 'L', 'L', 'N', 'U', N-1, NRHS, ONE, A( 2, 1), LDA,
-     $               B(2, 1), LDB)
+*        2) Solve with triangular matrix T
 *
 *        Compute T \ B -> B   [ T \ (L \P**T * B) ]
 *
          CALL CLACPY( 'F', 1, N, A(1, 1), LDA+1, WORK(N), 1)
          IF( N.GT.1 ) THEN
-             CALL CLACPY( 'F', 1, N-1, A( 2, 1 ), LDA+1, WORK( 1 ), 1)
+             CALL CLACPY( 'F', 1, N-1, A( 2, 1 ), LDA+1, WORK( 1 ), 1 )
              CALL CLACPY( 'F', 1, N-1, A( 2, 1 ), LDA+1, WORK( 2*N ), 1)
              CALL CLACGV( N-1, WORK( 2*N ), 1 )
          END IF
          CALL CGTSV(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
 *
-         CALL CTRSM( 'L', 'L', 'C', 'U', N-1, NRHS, ONE, A( 2, 1 ), LDA,
-     $              B( 2, 1 ), LDB)
+         IF( N.GT.1 ) THEN
+*
+*           Compute (L**H \ B) -> B   [ L**H \ (T \ (L \P**T * B) ) ]
 *
-*        Pivot, P * B  [ P * (L**T \ (T \ (L \P**T * B) )) ]
+            CALL CTRSM( 'L', 'L', 'C', 'U', N-1, NRHS, ONE, A( 2, 1 ),
+     $                  LDA, B( 2, 1 ), LDB )
 *
-         K = N
-         DO WHILE ( K.GE.1 )
-            KP = IPIV( K )
-            IF( KP.NE.K )
-     $         CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
-            K = K - 1
-         END DO
+*           Pivot, P * B -> B  [ P * (L**H \ (T \ (L \P**T * B) )) ]
+*
+            K = N
+            DO WHILE ( K.GE.1 )
+               KP = IPIV( K )
+               IF( KP.NE.K )
+     $            CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+               K = K - 1
+            END DO
+         END IF
 *
       END IF
 *

SRC/chetrs_aa_2stage.f

@@ -38,7 +38,7 @@
 *> \verbatim
 *>
 *> CHETRS_AA_2STAGE solves a system of linear equations A*X = B with a real
-*> hermitian matrix A using the factorization A = U*T*U**T or
+*> hermitian matrix A using the factorization A = U**T*T*U or
 *> A = L*T*L**T computed by CHETRF_AA_2STAGE.
 *> \endverbatim
 *
@@ -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**T;
+*>          = 'U':  Upper triangular, form is A = U**T*T*U;
 *>          = 'L':  Lower triangular, form is A = L*T*L**T.
 *> \endverbatim
 *>
@@ -210,15 +210,15 @@ SUBROUTINE CHETRS_AA_2STAGE( UPLO, N, NRHS, A, LDA, TB, LTB,
 *
       IF( UPPER ) THEN
 *
-*        Solve A*X = B, where A = U*T*U**T.
+*        Solve A*X = B, where A = U**T*T*U.
 *
          IF( N.GT.NB ) THEN
 *
-*           Pivot, P**T * B
+*           Pivot, P**T * B -> B
 *
             CALL CLASWP( NRHS, B, LDB, NB+1, N, IPIV, 1 )
 *
-*           Compute (U**T \P**T * B) -> B    [ (U**T \P**T * B) ]
+*           Compute (U**T \ B) -> B    [ (U**T \P**T * B) ]
 *
             CALL CTRSM( 'L', 'U', 'C', 'U', N-NB, NRHS, ONE, A(1, NB+1),
      $                 LDA, B(NB+1, 1), LDB)

SRC/csysv_aa.f

@@ -42,7 +42,7 @@
 *> matrices.
 *>
 *> Aasen's algorithm is used to factor A as
-*>    A = U * T * U**T,  if UPLO = 'U', or
+*>    A = U**T * T * U,  if UPLO = 'U', or
 *>    A = L * T * L**T,  if UPLO = 'L',
 *> where U (or L) is a product of permutation and unit upper (lower)
 *> triangular matrices, and T is symmetric tridiagonal. The factored
@@ -86,7 +86,7 @@
 *>
 *>          On exit, if INFO = 0, the tridiagonal matrix T and the
 *>          multipliers used to obtain the factor U or L from the
-*>          factorization A = U*T*U**T or A = L*T*L**T as computed by
+*>          factorization A = U**T*T*U or A = L*T*L**T as computed by
 *>          CSYTRF.
 *> \endverbatim
 *>
@@ -230,7 +230,7 @@ SUBROUTINE CSYSV_AA( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, WORK,
          RETURN
       END IF
 *
-*     Compute the factorization A = U*T*U**T or A = L*T*L**T.
+*     Compute the factorization A = U**T*T*U or A = L*T*L**T.
 *
       CALL CSYTRF_AA( UPLO, N, A, LDA, IPIV, WORK, LWORK, INFO )
       IF( INFO.EQ.0 ) THEN

SRC/csysv_aa_2stage.f

@@ -43,8 +43,8 @@
 *> matrices.
 *>
 *> Aasen's 2-stage algorithm is used to factor A as
-*>    A = U * T * U**H,  if UPLO = 'U', or
-*>    A = L * T * L**H,  if UPLO = 'L',
+*>    A = U**T * T * U,  if UPLO = 'U', or
+*>    A = L * T * L**T,  if UPLO = 'L',
 *> where U (or L) is a product of permutation and unit upper (lower)
 *> triangular matrices, and T is symmetric and band. The matrix T is
 *> then LU-factored with partial pivoting. The factored form of A
@@ -257,7 +257,7 @@ SUBROUTINE CSYSV_AA_2STAGE( UPLO, N, NRHS, A, LDA, TB, LTB,
       END IF
 *
 *
-*     Compute the factorization A = U*T*U**H or A = L*T*L**H.
+*     Compute the factorization A = U**T*T*U or A = L*T*L**T.
 *
       CALL CSYTRF_AA_2STAGE( UPLO, N, A, LDA, TB, LTB, IPIV, IPIV2,
      $                       WORK, LWORK, INFO )

SRC/csytrf_aa.f

@@ -37,7 +37,7 @@
 *> CSYTRF_AA computes the factorization of a complex 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 complex symmetric tridiagonal matrix.
@@ -223,7 +223,7 @@ SUBROUTINE CSYTRF_AA( UPLO, N, A, LDA, IPIV, WORK, LWORK, INFO)
       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
 *        .....................................................
 *
 *        Copy first row A(1, 1:N) into H(1:n) (stored in WORK(1:N))

SRC/csytrf_aa_2stage.f

@@ -38,7 +38,7 @@
 *> CSYTRF_AA_2STAGE computes the factorization of a complex 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 complex symmetric band matrix with the
@@ -275,7 +275,7 @@ SUBROUTINE CSYTRF_AA_2STAGE( UPLO, N, A, LDA, TB, LTB, IPIV,
       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