Move LAPACK trunk into position.

1 files changed, 401 insertions, 0 deletions

diff --git a/SRC/zhptrs.f b/SRC/zhptrs.f
new file mode 100644
index 00000000..70719393
--- /dev/null
+++ b/SRC/zhptrs.f
@@ -0,0 +1,401 @@
+      SUBROUTINE ZHPTRS( UPLO, N, NRHS, AP, IPIV, B, LDB, INFO )
+*
+*  -- LAPACK routine (version 3.1) --
+*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
+*     November 2006
+*
+*     .. Scalar Arguments ..
+      CHARACTER          UPLO
+      INTEGER            INFO, LDB, N, NRHS
+*     ..
+*     .. Array Arguments ..
+      INTEGER            IPIV( * )
+      COMPLEX*16         AP( * ), B( LDB, * )
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  ZHPTRS solves a system of linear equations A*X = B with a complex
+*  Hermitian matrix A stored in packed format using the factorization
+*  A = U*D*U**H or A = L*D*L**H computed by ZHPTRF.
+*
+*  Arguments
+*  =========
+*
+*  UPLO    (input) CHARACTER*1
+*          Specifies whether the details of the factorization are stored
+*          as an upper or lower triangular matrix.
+*          = 'U':  Upper triangular, form is A = U*D*U**H;
+*          = 'L':  Lower triangular, form is A = L*D*L**H.
+*
+*  N       (input) INTEGER
+*          The order of the matrix A.  N >= 0.
+*
+*  NRHS    (input) INTEGER
+*          The number of right hand sides, i.e., the number of columns
+*          of the matrix B.  NRHS >= 0.
+*
+*  AP      (input) COMPLEX*16 array, dimension (N*(N+1)/2)
+*          The block diagonal matrix D and the multipliers used to
+*          obtain the factor U or L as computed by ZHPTRF, stored as a
+*          packed triangular matrix.
+*
+*  IPIV    (input) INTEGER array, dimension (N)
+*          Details of the interchanges and the block structure of D
+*          as determined by ZHPTRF.
+*
+*  B       (input/output) COMPLEX*16 array, dimension (LDB,NRHS)
+*          On entry, the right hand side matrix B.
+*          On exit, the solution matrix X.
+*
+*  LDB     (input) INTEGER
+*          The leading dimension of the array B.  LDB >= max(1,N).
+*
+*  INFO    (output) INTEGER
+*          = 0:  successful exit
+*          < 0: if INFO = -i, the i-th argument had an illegal value
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      COMPLEX*16         ONE
+      PARAMETER          ( ONE = ( 1.0D+0, 0.0D+0 ) )
+*     ..
+*     .. Local Scalars ..
+      LOGICAL            UPPER
+      INTEGER            J, K, KC, KP
+      DOUBLE PRECISION   S
+      COMPLEX*16         AK, AKM1, AKM1K, BK, BKM1, DENOM
+*     ..
+*     .. External Functions ..
+      LOGICAL            LSAME
+      EXTERNAL           LSAME
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           XERBLA, ZDSCAL, ZGEMV, ZGERU, ZLACGV, ZSWAP
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          DBLE, DCONJG, MAX
+*     ..
+*     .. Executable Statements ..
+*
+      INFO = 0
+      UPPER = LSAME( UPLO, 'U' )
+      IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
+         INFO = -1
+      ELSE IF( N.LT.0 ) THEN
+         INFO = -2
+      ELSE IF( NRHS.LT.0 ) THEN
+         INFO = -3
+      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
+         INFO = -7
+      END IF
+      IF( INFO.NE.0 ) THEN
+         CALL XERBLA( 'ZHPTRS', -INFO )
+         RETURN
+      END IF
+*
+*     Quick return if possible
+*
+      IF( N.EQ.0 .OR. NRHS.EQ.0 )
+     $   RETURN
+*
+      IF( UPPER ) THEN
+*
+*        Solve A*X = B, where A = U*D*U'.
+*
+*        First solve U*D*X = B, overwriting B with X.
+*
+*        K is the main loop index, decreasing from N to 1 in steps of
+*        1 or 2, depending on the size of the diagonal blocks.
+*
+         K = N
+         KC = N*( N+1 ) / 2 + 1
+   10    CONTINUE
+*
+*        If K < 1, exit from loop.
+*
+         IF( K.LT.1 )
+     $      GO TO 30
+*
+         KC = KC - K
+         IF( IPIV( K ).GT.0 ) THEN
+*
+*           1 x 1 diagonal block
+*
+*           Interchange rows K and IPIV(K).
+*
+            KP = IPIV( K )
+            IF( KP.NE.K )
+     $         CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+*
+*           Multiply by inv(U(K)), where U(K) is the transformation
+*           stored in column K of A.
+*
+            CALL ZGERU( K-1, NRHS, -ONE, AP( KC ), 1, B( K, 1 ), LDB,
+     $                  B( 1, 1 ), LDB )
+*
+*           Multiply by the inverse of the diagonal block.
+*
+            S = DBLE( ONE ) / DBLE( AP( KC+K-1 ) )
+            CALL ZDSCAL( NRHS, S, B( K, 1 ), LDB )
+            K = K - 1
+         ELSE
+*
+*           2 x 2 diagonal block
+*
+*           Interchange rows K-1 and -IPIV(K).
+*
+            KP = -IPIV( K )
+            IF( KP.NE.K-1 )
+     $         CALL ZSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ), LDB )
+*
+*           Multiply by inv(U(K)), where U(K) is the transformation
+*           stored in columns K-1 and K of A.
+*
+            CALL ZGERU( K-2, NRHS, -ONE, AP( KC ), 1, B( K, 1 ), LDB,
+     $                  B( 1, 1 ), LDB )
+            CALL ZGERU( K-2, NRHS, -ONE, AP( KC-( K-1 ) ), 1,
+     $                  B( K-1, 1 ), LDB, B( 1, 1 ), LDB )
+*
+*           Multiply by the inverse of the diagonal block.
+*
+            AKM1K = AP( KC+K-2 )
+            AKM1 = AP( KC-1 ) / AKM1K
+            AK = AP( KC+K-1 ) / DCONJG( AKM1K )
+            DENOM = AKM1*AK - ONE
+            DO 20 J = 1, NRHS
+               BKM1 = B( K-1, J ) / AKM1K
+               BK = B( K, J ) / DCONJG( AKM1K )
+               B( K-1, J ) = ( AK*BKM1-BK ) / DENOM
+               B( K, J ) = ( AKM1*BK-BKM1 ) / DENOM
+   20       CONTINUE
+            KC = KC - K + 1
+            K = K - 2
+         END IF
+*
+         GO TO 10
+   30    CONTINUE
+*
+*        Next solve U'*X = B, overwriting B with X.
+*
+*        K is the main loop index, increasing from 1 to N in steps of
+*        1 or 2, depending on the size of the diagonal blocks.
+*
+         K = 1
+         KC = 1
+   40    CONTINUE
+*
+*        If K > N, exit from loop.
+*
+         IF( K.GT.N )
+     $      GO TO 50
+*
+         IF( IPIV( K ).GT.0 ) THEN
+*
+*           1 x 1 diagonal block
+*
+*           Multiply by inv(U'(K)), where U(K) is the transformation
+*           stored in column K of A.
+*
+            IF( K.GT.1 ) THEN
+               CALL ZLACGV( NRHS, B( K, 1 ), LDB )
+               CALL ZGEMV( 'Conjugate transpose', K-1, NRHS, -ONE, B,
+     $                     LDB, AP( KC ), 1, ONE, B( K, 1 ), LDB )
+               CALL ZLACGV( NRHS, B( K, 1 ), LDB )
+            END IF
+*
+*           Interchange rows K and IPIV(K).
+*
+            KP = IPIV( K )
+            IF( KP.NE.K )
+     $         CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+            KC = KC + K
+            K = K + 1
+         ELSE
+*
+*           2 x 2 diagonal block
+*
+*           Multiply by inv(U'(K+1)), where U(K+1) is the transformation
+*           stored in columns K and K+1 of A.
+*
+            IF( K.GT.1 ) THEN
+               CALL ZLACGV( NRHS, B( K, 1 ), LDB )
+               CALL ZGEMV( 'Conjugate transpose', K-1, NRHS, -ONE, B,
+     $                     LDB, AP( KC ), 1, ONE, B( K, 1 ), LDB )
+               CALL ZLACGV( NRHS, B( K, 1 ), LDB )
+*
+               CALL ZLACGV( NRHS, B( K+1, 1 ), LDB )
+               CALL ZGEMV( 'Conjugate transpose', K-1, NRHS, -ONE, B,
+     $                     LDB, AP( KC+K ), 1, ONE, B( K+1, 1 ), LDB )
+               CALL ZLACGV( NRHS, B( K+1, 1 ), LDB )
+            END IF
+*
+*           Interchange rows K and -IPIV(K).
+*
+            KP = -IPIV( K )
+            IF( KP.NE.K )
+     $         CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+            KC = KC + 2*K + 1
+            K = K + 2
+         END IF
+*
+         GO TO 40
+   50    CONTINUE
+*
+      ELSE
+*
+*        Solve A*X = B, where A = L*D*L'.
+*
+*        First solve L*D*X = B, overwriting B with X.
+*
+*        K is the main loop index, increasing from 1 to N in steps of
+*        1 or 2, depending on the size of the diagonal blocks.
+*
+         K = 1
+         KC = 1
+   60    CONTINUE
+*
+*        If K > N, exit from loop.
+*
+         IF( K.GT.N )
+     $      GO TO 80
+*
+         IF( IPIV( K ).GT.0 ) THEN
+*
+*           1 x 1 diagonal block
+*
+*           Interchange rows K and IPIV(K).
+*
+            KP = IPIV( K )
+            IF( KP.NE.K )
+     $         CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+*
+*           Multiply by inv(L(K)), where L(K) is the transformation
+*           stored in column K of A.
+*
+            IF( K.LT.N )
+     $         CALL ZGERU( N-K, NRHS, -ONE, AP( KC+1 ), 1, B( K, 1 ),
+     $                     LDB, B( K+1, 1 ), LDB )
+*
+*           Multiply by the inverse of the diagonal block.
+*
+            S = DBLE( ONE ) / DBLE( AP( KC ) )
+            CALL ZDSCAL( NRHS, S, B( K, 1 ), LDB )
+            KC = KC + N - K + 1
+            K = K + 1
+         ELSE
+*
+*           2 x 2 diagonal block
+*
+*           Interchange rows K+1 and -IPIV(K).
+*
+            KP = -IPIV( K )
+            IF( KP.NE.K+1 )
+     $         CALL ZSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ), LDB )
+*
+*           Multiply by inv(L(K)), where L(K) is the transformation
+*           stored in columns K and K+1 of A.
+*
+            IF( K.LT.N-1 ) THEN
+               CALL ZGERU( N-K-1, NRHS, -ONE, AP( KC+2 ), 1, B( K, 1 ),
+     $                     LDB, B( K+2, 1 ), LDB )
+               CALL ZGERU( N-K-1, NRHS, -ONE, AP( KC+N-K+2 ), 1,
+     $                     B( K+1, 1 ), LDB, B( K+2, 1 ), LDB )
+            END IF
+*
+*           Multiply by the inverse of the diagonal block.
+*
+            AKM1K = AP( KC+1 )
+            AKM1 = AP( KC ) / DCONJG( AKM1K )
+            AK = AP( KC+N-K+1 ) / AKM1K
+            DENOM = AKM1*AK - ONE
+            DO 70 J = 1, NRHS
+               BKM1 = B( K, J ) / DCONJG( AKM1K )
+               BK = B( K+1, J ) / AKM1K
+               B( K, J ) = ( AK*BKM1-BK ) / DENOM
+               B( K+1, J ) = ( AKM1*BK-BKM1 ) / DENOM
+   70       CONTINUE
+            KC = KC + 2*( N-K ) + 1
+            K = K + 2
+         END IF
+*
+         GO TO 60
+   80    CONTINUE
+*
+*        Next solve L'*X = B, overwriting B with X.
+*
+*        K is the main loop index, decreasing from N to 1 in steps of
+*        1 or 2, depending on the size of the diagonal blocks.
+*
+         K = N
+         KC = N*( N+1 ) / 2 + 1
+   90    CONTINUE
+*
+*        If K < 1, exit from loop.
+*
+         IF( K.LT.1 )
+     $      GO TO 100
+*
+         KC = KC - ( N-K+1 )
+         IF( IPIV( K ).GT.0 ) THEN
+*
+*           1 x 1 diagonal block
+*
+*           Multiply by inv(L'(K)), where L(K) is the transformation
+*           stored in column K of A.
+*
+            IF( K.LT.N ) THEN
+               CALL ZLACGV( NRHS, B( K, 1 ), LDB )
+               CALL ZGEMV( 'Conjugate transpose', N-K, NRHS, -ONE,
+     $                     B( K+1, 1 ), LDB, AP( KC+1 ), 1, ONE,
+     $                     B( K, 1 ), LDB )
+               CALL ZLACGV( NRHS, B( K, 1 ), LDB )
+            END IF
+*
+*           Interchange rows K and IPIV(K).
+*
+            KP = IPIV( K )
+            IF( KP.NE.K )
+     $         CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+            K = K - 1
+         ELSE
+*
+*           2 x 2 diagonal block
+*
+*           Multiply by inv(L'(K-1)), where L(K-1) is the transformation
+*           stored in columns K-1 and K of A.
+*
+            IF( K.LT.N ) THEN
+               CALL ZLACGV( NRHS, B( K, 1 ), LDB )
+               CALL ZGEMV( 'Conjugate transpose', N-K, NRHS, -ONE,
+     $                     B( K+1, 1 ), LDB, AP( KC+1 ), 1, ONE,
+     $                     B( K, 1 ), LDB )
+               CALL ZLACGV( NRHS, B( K, 1 ), LDB )
+*
+               CALL ZLACGV( NRHS, B( K-1, 1 ), LDB )
+               CALL ZGEMV( 'Conjugate transpose', N-K, NRHS, -ONE,
+     $                     B( K+1, 1 ), LDB, AP( KC-( N-K ) ), 1, ONE,
+     $                     B( K-1, 1 ), LDB )
+               CALL ZLACGV( NRHS, B( K-1, 1 ), LDB )
+            END IF
+*
+*           Interchange rows K and -IPIV(K).
+*
+            KP = -IPIV( K )
+            IF( KP.NE.K )
+     $         CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+            KC = KC - ( N-K+2 )
+            K = K - 2
+         END IF
+*
+         GO TO 90
+  100    CONTINUE
+      END IF
+*
+      RETURN
+*
+*     End of ZHPTRS
+*
+      END