Add hetrs2 for c and z, a Level BLAS 3 version of hetrs (same than rev 751)
authorjulie <julielangou@users.noreply.github.com>
Thu, 30 Sep 2010 06:58:55 +0000 (06:58 +0000)
committerjulie <julielangou@users.noreply.github.com>
Thu, 30 Sep 2010 06:58:55 +0000 (06:58 +0000)
SRC/Makefile
SRC/chesv.f
SRC/chetrs2.f [new file with mode: 0644]
SRC/zhesv.f
SRC/zhetrs2.f [new file with mode: 0644]
TESTING/LIN/alahd.f
TESTING/LIN/cchkhe.f
TESTING/LIN/zchkhe.f

index 2cf9323..52db7fc 100644 (file)
@@ -169,7 +169,7 @@ CLASRC = \
    checon.o cheev.o  cheevd.o cheevr.o cheevx.o chegs2.o chegst.o \
    chegv.o  chegvd.o chegvx.o cherfs.o chesv.o  chesvx.o chetd2.o \
    chetf2.o chetrd.o \
-   chetrf.o chetri.o chetrs.o chgeqz.o chpcon.o chpev.o  chpevd.o \
+   chetrf.o chetri.o chetrs.o chetrs2.o chgeqz.o chpcon.o chpev.o  chpevd.o \
    chpevx.o chpgst.o chpgv.o  chpgvd.o chpgvx.o chprfs.o chpsv.o  \
    chpsvx.o \
    chptrd.o chptrf.o chptri.o chptrs.o chsein.o chseqr.o clabrd.o \
@@ -305,7 +305,7 @@ ZLASRC = \
    zhecon.o zheev.o  zheevd.o zheevr.o zheevx.o zhegs2.o zhegst.o \
    zhegv.o  zhegvd.o zhegvx.o zherfs.o zhesv.o  zhesvx.o zhetd2.o \
    zhetf2.o zhetrd.o \
-   zhetrf.o zhetri.o zhetrs.o zhgeqz.o zhpcon.o zhpev.o  zhpevd.o \
+   zhetrf.o zhetri.o zhetrs.o zhetrs2.o zhgeqz.o zhpcon.o zhpev.o  zhpevd.o \
    zhpevx.o zhpgst.o zhpgv.o  zhpgvd.o zhpgvx.o zhprfs.o zhpsv.o  \
    zhpsvx.o \
    zhptrd.o zhptrf.o zhptri.o zhptrs.o zhsein.o zhseqr.o zlabrd.o \
index c2c1975..47d5f68 100644 (file)
 *
 *     .. Local Scalars ..
       LOGICAL            LQUERY
-      INTEGER            LWKOPT, NB
+      INTEGER            IINFO, LWKOPT, NB
 *     ..
 *     .. External Functions ..
       LOGICAL            LSAME
       EXTERNAL           ILAENV, LSAME
 *     ..
 *     .. External Subroutines ..
-      EXTERNAL           CHETRF, CHETRS, XERBLA
+      EXTERNAL           CHETRF, CHETRS2, CSYCONV, XERBLA
 *     ..
 *     .. Intrinsic Functions ..
       INTRINSIC          MAX
       CALL CHETRF( UPLO, N, A, LDA, IPIV, WORK, LWORK, INFO )
       IF( INFO.EQ.0 ) THEN
 *
+*        Convert A
+*
+         CALL CSYCONV( UPLO, 'C', N, A, LDA, IPIV, WORK, IINFO )
+*
 *        Solve the system A*X = B, overwriting B with X.
 *
-         CALL CHETRS( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, INFO )
+         CALL CHETRS2( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, WORK, INFO )
+*
+*        Revert A
+*
+         CALL CSYCONV( UPLO, 'R', N, A, LDA, IPIV, WORK, IINFO )
 *
       END IF
 *
diff --git a/SRC/chetrs2.f b/SRC/chetrs2.f
new file mode 100644 (file)
index 0000000..1204a59
--- /dev/null
@@ -0,0 +1,275 @@
+      SUBROUTINE CHETRS2( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, 
+     $                    WORK, INFO )
+*
+*  -- LAPACK PROTOTYPE routine (version 3.2.2) --
+*
+*  -- Written by Julie Langou of the Univ. of TN    --
+*     May 2010
+*
+*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*
+*     .. Scalar Arguments ..
+      CHARACTER          UPLO
+      INTEGER            INFO, LDA, LDB, N, NRHS
+*     ..
+*     .. Array Arguments ..
+      INTEGER            IPIV( * )
+      COMPLEX            A( LDA, * ), B( LDB, * ), WORK( * )
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  CHETRS2 solves a system of linear equations A*X = B with a COMPLEX
+*  Hermitian matrix A using the factorization A = U*D*U**T or
+*  A = L*D*L**T computed by CSYTRF and converted by CSYCONV.
+*
+*  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.
+*
+*  A       (input) COMPLEX array, dimension (LDA,N)
+*          The block diagonal matrix D and the multipliers used to
+*          obtain the factor U or L as computed by CHETRF.
+*
+*  LDA     (input) INTEGER
+*          The leading dimension of the array A.  LDA >= max(1,N).
+*
+*  IPIV    (input) INTEGER array, dimension (N)
+*          Details of the interchanges and the block structure of D
+*          as determined by CHETRF.
+*
+*  B       (input/output) COMPLEX 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).
+*
+*  WORK    (workspace) COMPLEX array, dimension (N)
+*
+*  INFO    (output) INTEGER
+*          = 0:  successful exit
+*          < 0:  if INFO = -i, the i-th argument had an illegal value
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      COMPLEX            ONE
+      PARAMETER          ( ONE = (1.0E+0,0.0E+0) )
+*     ..
+*     .. Local Scalars ..
+      LOGICAL            UPPER
+      INTEGER            I, J, K, KP
+      REAL               S
+      COMPLEX            AK, AKM1, AKM1K, BK, BKM1, DENOM
+*     ..
+*     .. External Functions ..
+      LOGICAL            LSAME
+      EXTERNAL           LSAME
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           CSCAL, CSWAP, CTRSM, XERBLA
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          CONJG, MAX, REAL
+*     ..
+*     .. 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( LDA.LT.MAX( 1, N ) ) THEN
+         INFO = -5
+      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
+         INFO = -8
+      END IF
+      IF( INFO.NE.0 ) THEN
+         CALL XERBLA( 'CHETRS2', -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'.
+*
+*       P' * B  
+        K=N
+        DO WHILE ( K .GE. 1 )
+         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 CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+            K=K-1
+         ELSE
+*           2 x 2 diagonal block
+*           Interchange rows K-1 and -IPIV(K).
+            KP = -IPIV( K )
+            IF( KP.EQ.-IPIV( K-1 ) )
+     $         CALL CSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ), LDB )
+            K=K-2
+         END IF
+        END DO
+*
+*  Compute (U \P' * B) -> B    [ (U \P' * B) ]
+*
+        CALL CTRSM('L','U','N','U',N,NRHS,ONE,A,N,B,N)
+*
+*  Compute D \ B -> B   [ D \ (U \P' * B) ]
+*       
+         I=N
+         DO WHILE ( I .GE. 1 )
+            IF( IPIV(I) .GT. 0 ) THEN
+              S = REAL( ONE ) / REAL( A( I, I ) )
+              CALL CSSCAL( NRHS, S, B( I, 1 ), LDB )
+            ELSEIF ( I .GT. 1) THEN
+               IF ( IPIV(I-1) .EQ. IPIV(I) ) THEN
+                  AKM1K = WORK(I)
+                  AKM1 = A( I-1, I-1 ) / AKM1K
+                  AK = A( I, I ) / CONJG( AKM1K )
+                  DENOM = AKM1*AK - ONE
+                  DO 15 J = 1, NRHS
+                     BKM1 = B( I-1, J ) / AKM1K
+                     BK = B( I, J ) / CONJG( AKM1K )
+                     B( I-1, J ) = ( AK*BKM1-BK ) / DENOM
+                     B( I, J ) = ( AKM1*BK-BKM1 ) / DENOM
+ 15              CONTINUE
+               I = I - 1
+               ENDIF
+            ENDIF
+            I = I - 1
+         END DO
+*
+*      Compute (U' \ B) -> B   [ U' \ (D \ (U \P' * B) ) ]
+*
+         CALL CTRSM('L','U','C','U',N,NRHS,ONE,A,N,B,N)
+*
+*       P * B  [ P * (U' \ (D \ (U \P' * B) )) ]
+*
+        K=1
+        DO WHILE ( K .LE. N )
+         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 CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+            K=K+1
+         ELSE
+*           2 x 2 diagonal block
+*           Interchange rows K-1 and -IPIV(K).
+            KP = -IPIV( K )
+            IF( K .LT. N .AND. KP.EQ.-IPIV( K+1 ) )
+     $         CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+            K=K+2
+         ENDIF
+        END DO
+*
+      ELSE
+*
+*        Solve A*X = B, where A = L*D*L'.
+*
+*       P' * B  
+        K=1
+        DO WHILE ( K .LE. N )
+         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 CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+            K=K+1
+         ELSE
+*           2 x 2 diagonal block
+*           Interchange rows K and -IPIV(K+1).
+            KP = -IPIV( K+1 )
+            IF( KP.EQ.-IPIV( K ) )
+     $         CALL CSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ), LDB )
+            K=K+2
+         ENDIF
+        END DO
+*
+*  Compute (L \P' * B) -> B    [ (L \P' * B) ]
+*
+        CALL CTRSM('L','L','N','U',N,NRHS,ONE,A,N,B,N)
+*
+*  Compute D \ B -> B   [ D \ (L \P' * B) ]
+*       
+         I=1
+         DO WHILE ( I .LE. N )
+            IF( IPIV(I) .GT. 0 ) THEN
+              S = REAL( ONE ) / REAL( A( I, I ) )
+              CALL CSSCAL( NRHS, S, B( I, 1 ), LDB )
+            ELSE
+                  AKM1K = WORK(I)
+                  AKM1 = A( I, I ) / CONJG( AKM1K )
+                  AK = A( I+1, I+1 ) / AKM1K
+                  DENOM = AKM1*AK - ONE
+                  DO 25 J = 1, NRHS
+                     BKM1 = B( I, J ) / CONJG( AKM1K )
+                     BK = B( I+1, J ) / AKM1K
+                     B( I, J ) = ( AK*BKM1-BK ) / DENOM
+                     B( I+1, J ) = ( AKM1*BK-BKM1 ) / DENOM
+ 25              CONTINUE
+                  I = I + 1
+            ENDIF
+            I = I + 1
+         END DO
+*
+*  Compute (L' \ B) -> B   [ L' \ (D \ (L \P' * B) ) ]
+* 
+        CALL CTRSM('L','L','C','U',N,NRHS,ONE,A,N,B,N)
+*
+*       P * B  [ P * (L' \ (D \ (L \P' * B) )) ]
+*
+        K=N
+        DO WHILE ( K .GE. 1 )
+         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 CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+            K=K-1
+         ELSE
+*           2 x 2 diagonal block
+*           Interchange rows K-1 and -IPIV(K).
+            KP = -IPIV( K )
+            IF( K.GT.1 .AND. KP.EQ.-IPIV( K-1 ) )
+     $         CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+            K=K-2
+         ENDIF
+        END DO
+*
+      END IF
+*
+      RETURN
+*
+*     End of CHETRS2
+*
+      END
index e457d98..f539a27 100644 (file)
 *
 *     .. Local Scalars ..
       LOGICAL            LQUERY
-      INTEGER            LWKOPT, NB
+      INTEGER            IINFO, LWKOPT, NB
 *     ..
 *     .. External Functions ..
       LOGICAL            LSAME
       EXTERNAL           LSAME, ILAENV
 *     ..
 *     .. External Subroutines ..
-      EXTERNAL           XERBLA, ZHETRF, ZHETRS
+      EXTERNAL           XERBLA, ZHETRF, ZHETRS2, ZSYCONV
 *     ..
 *     .. Intrinsic Functions ..
       INTRINSIC          MAX
       CALL ZHETRF( UPLO, N, A, LDA, IPIV, WORK, LWORK, INFO )
       IF( INFO.EQ.0 ) THEN
 *
+*        Convert A
+*
+         CALL ZSYCONV( UPLO, 'C', N, A, LDA, IPIV, WORK, IINFO )
+*
 *        Solve the system A*X = B, overwriting B with X.
 *
-         CALL ZHETRS( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, INFO )
+         CALL ZHETRS2( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, WORK, INFO )
+*
+*        Revert A
+*
+         CALL ZSYCONV( UPLO, 'R', N, A, LDA, IPIV, WORK, IINFO )
 *
       END IF
 *
diff --git a/SRC/zhetrs2.f b/SRC/zhetrs2.f
new file mode 100644 (file)
index 0000000..91a53cd
--- /dev/null
@@ -0,0 +1,275 @@
+      SUBROUTINE ZHETRS2( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, 
+     $                    WORK, INFO )
+*
+*  -- LAPACK PROTOTYPE routine (version 3.2.2) --
+*
+*  -- Written by Julie Langou of the Univ. of TN    --
+*     May 2010
+*
+*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*
+*     .. Scalar Arguments ..
+      CHARACTER          UPLO
+      INTEGER            INFO, LDA, LDB, N, NRHS
+*     ..
+*     .. Array Arguments ..
+      INTEGER            IPIV( * )
+      DOUBLE COMPLEX   A( LDA, * ), B( LDB, * ), WORK( * )
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  ZHETRS2 solves a system of linear equations A*X = B with a real
+*  Hermitian matrix A using the factorization A = U*D*U**T or
+*  A = L*D*L**T computed by ZSYTRF and converted by ZSYCONV.
+*
+*  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.
+*
+*  A       (input) DOUBLE COMPLEX array, dimension (LDA,N)
+*          The block diagonal matrix D and the multipliers used to
+*          obtain the factor U or L as computed by ZHETRF.
+*
+*  LDA     (input) INTEGER
+*          The leading dimension of the array A.  LDA >= max(1,N).
+*
+*  IPIV    (input) INTEGER array, dimension (N)
+*          Details of the interchanges and the block structure of D
+*          as determined by ZHETRF.
+*
+*  B       (input/output) DOUBLE COMPLEX 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).
+*
+*  WORK    (workspace) REAL array, dimension (N)
+*
+*  INFO    (output) INTEGER
+*          = 0:  successful exit
+*          < 0:  if INFO = -i, the i-th argument had an illegal value
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      DOUBLE COMPLEX     ONE
+      PARAMETER          ( ONE = (1.0D+0,0.0D+0) )
+*     ..
+*     .. Local Scalars ..
+      LOGICAL            UPPER
+      INTEGER            I, J, K, KP
+      DOUBLE PRECISION   S
+      DOUBLE COMPLEX     AK, AKM1, AKM1K, BK, BKM1, DENOM
+*     ..
+*     .. External Functions ..
+      LOGICAL            LSAME
+      EXTERNAL           LSAME
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           ZLACGV, ZSCAL, ZSWAP, ZTRSM, XERBLA
+*     ..
+*     .. 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( LDA.LT.MAX( 1, N ) ) THEN
+         INFO = -5
+      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
+         INFO = -8
+      END IF
+      IF( INFO.NE.0 ) THEN
+         CALL XERBLA( 'ZHETRS2', -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'.
+*
+*       P' * B  
+        K=N
+        DO WHILE ( K .GE. 1 )
+         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 )
+            K=K-1
+         ELSE
+*           2 x 2 diagonal block
+*           Interchange rows K-1 and -IPIV(K).
+            KP = -IPIV( K )
+            IF( KP.EQ.-IPIV( K-1 ) )
+     $         CALL ZSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ), LDB )
+            K=K-2
+         END IF
+        END DO
+*
+*  Compute (U \P' * B) -> B    [ (U \P' * B) ]
+*
+        CALL ZTRSM('L','U','N','U',N,NRHS,ONE,A,N,B,N)
+*
+*  Compute D \ B -> B   [ D \ (U \P' * B) ]
+*       
+         I=N
+         DO WHILE ( I .GE. 1 )
+            IF( IPIV(I) .GT. 0 ) THEN
+              S = DBLE( ONE ) / DBLE( A( I, I ) )
+              CALL ZDSCAL( NRHS, S, B( I, 1 ), LDB )
+            ELSEIF ( I .GT. 1) THEN
+               IF ( IPIV(I-1) .EQ. IPIV(I) ) THEN
+                  AKM1K = WORK(I)
+                  AKM1 = A( I-1, I-1 ) / AKM1K
+                  AK = A( I, I ) / DCONJG( AKM1K )
+                  DENOM = AKM1*AK - ONE
+                  DO 15 J = 1, NRHS
+                     BKM1 = B( I-1, J ) / AKM1K
+                     BK = B( I, J ) / DCONJG( AKM1K )
+                     B( I-1, J ) = ( AK*BKM1-BK ) / DENOM
+                     B( I, J ) = ( AKM1*BK-BKM1 ) / DENOM
+ 15              CONTINUE
+               I = I - 1
+               ENDIF
+            ENDIF
+            I = I - 1
+         END DO
+*
+*      Compute (U' \ B) -> B   [ U' \ (D \ (U \P' * B) ) ]
+*
+         CALL ZTRSM('L','U','C','U',N,NRHS,ONE,A,N,B,N)
+*
+*       P * B  [ P * (U' \ (D \ (U \P' * B) )) ]
+*
+        K=1
+        DO WHILE ( K .LE. N )
+         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 )
+            K=K+1
+         ELSE
+*           2 x 2 diagonal block
+*           Interchange rows K-1 and -IPIV(K).
+            KP = -IPIV( K )
+            IF( K .LT. N .AND. KP.EQ.-IPIV( K+1 ) )
+     $         CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+            K=K+2
+         ENDIF
+        END DO
+*
+      ELSE
+*
+*        Solve A*X = B, where A = L*D*L'.
+*
+*       P' * B  
+        K=1
+        DO WHILE ( K .LE. N )
+         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 )
+            K=K+1
+         ELSE
+*           2 x 2 diagonal block
+*           Interchange rows K and -IPIV(K+1).
+            KP = -IPIV( K+1 )
+            IF( KP.EQ.-IPIV( K ) )
+     $         CALL ZSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ), LDB )
+            K=K+2
+         ENDIF
+        END DO
+*
+*  Compute (L \P' * B) -> B    [ (L \P' * B) ]
+*
+        CALL ZTRSM('L','L','N','U',N,NRHS,ONE,A,N,B,N)
+*
+*  Compute D \ B -> B   [ D \ (L \P' * B) ]
+*       
+         I=1
+         DO WHILE ( I .LE. N )
+            IF( IPIV(I) .GT. 0 ) THEN
+              S = DBLE( ONE ) / DBLE( A( I, I ) )
+              CALL ZDSCAL( NRHS, S, B( I, 1 ), LDB )
+            ELSE
+                  AKM1K = WORK(I)
+                  AKM1 = A( I, I ) / DCONJG( AKM1K )
+                  AK = A( I+1, I+1 ) / AKM1K
+                  DENOM = AKM1*AK - ONE
+                  DO 25 J = 1, NRHS
+                     BKM1 = B( I, J ) / DCONJG( AKM1K )
+                     BK = B( I+1, J ) / AKM1K
+                     B( I, J ) = ( AK*BKM1-BK ) / DENOM
+                     B( I+1, J ) = ( AKM1*BK-BKM1 ) / DENOM
+ 25              CONTINUE
+                  I = I + 1
+            ENDIF
+            I = I + 1
+         END DO
+*
+*  Compute (L' \ B) -> B   [ L' \ (D \ (L \P' * B) ) ]
+* 
+        CALL ZTRSM('L','L','C','U',N,NRHS,ONE,A,N,B,N)
+*
+*       P * B  [ P * (L' \ (D \ (L \P' * B) )) ]
+*
+        K=N
+        DO WHILE ( K .GE. 1 )
+         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 )
+            K=K-1
+         ELSE
+*           2 x 2 diagonal block
+*           Interchange rows K-1 and -IPIV(K).
+            KP = -IPIV( K )
+            IF( K.GT.1 .AND. KP.EQ.-IPIV( K-1 ) )
+     $         CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+            K=K-2
+         ENDIF
+        END DO
+*
+      END IF
+*
+      RETURN
+*
+*     End of ZHETRS2
+*
+      END
index 6046d45..56c893f 100644 (file)
          WRITE( IOUNIT, FMT = 9955 )8
          WRITE( IOUNIT, FMT = '( '' Messages:'' )' )
 *
-      ELSE IF( LSAMEN( 2, P2, 'HE' ) .OR. LSAMEN( 2, P2, 'HP' ) ) THEN
+      ELSE IF( LSAMEN( 2, P2, 'HE' )  ) THEN
 *
 *        HE: Hermitian indefinite full
+*
+         IF( LSAME( C3, 'E' ) ) THEN
+            WRITE( IOUNIT, FMT = 9992 )PATH, 'Hermitian'
+         ELSE
+            WRITE( IOUNIT, FMT = 9991 )PATH, 'Hermitian'
+         END IF
+         WRITE( IOUNIT, FMT = '( '' Matrix types:'' )' )
+         IF( SORD ) THEN
+            WRITE( IOUNIT, FMT = 9972 )
+         ELSE
+            WRITE( IOUNIT, FMT = 9971 )
+         END IF
+         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, 'HP' ) ) THEN
+*
 *        HP: Hermitian indefinite packed
 *
          IF( LSAME( C3, 'E' ) ) THEN
index 99ad04e..77b7e8f 100644 (file)
@@ -22,7 +22,7 @@
 *  Purpose
 *  =======
 *
-*  CCHKHE tests CHETRF, -TRI, -TRS, -RFS, and -CON.
+*  CCHKHE tests CHETRF, -TRI, -TRS, -TRS2, -RFS, and -CON.
 *
 *  Arguments
 *  =========
@@ -94,7 +94,7 @@
       INTEGER            NTYPES
       PARAMETER          ( NTYPES = 10 )
       INTEGER            NTESTS
-      PARAMETER          ( NTESTS = 8 )
+      PARAMETER          ( NTESTS = 9 )
 *     ..
 *     .. Local Scalars ..
       LOGICAL            TRFCON, ZEROT
      $                            LDA, RWORK, RESULT( 3 ) )
 *
 *+    TEST 4
+*                 Solve and compute residual for  A * X = B.
+*
+                     SRNAMT = 'CLARHS'
+                     CALL CLARHS( PATH, 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 = 'CHETRS2'
+                     CALL CSYCONV( UPLO, 'C', N, AFAC, LDA, IWORK, WORK, 
+     $                            INFO )
+                     CALL CHETRS2( UPLO, N, NRHS, AFAC, LDA, IWORK, X,
+     $                            LDA, WORK, INFO )
+                     CALL CSYCONV( UPLO, 'R', N, AFAC, LDA, IWORK, WORK,
+     $                            INFO )
+*
+*                 Check error code from CHETRS2.
+*
+                     IF( INFO.NE.0 )
+     $                  CALL ALAERH( PATH, 'CHETRS2', INFO, 0, UPLO, N,
+     $                               N, -1, -1, NRHS, IMAT, NFAIL,
+     $                               NERRS, NOUT )
+*
+                     CALL CLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA )
+                     CALL CPOT02( UPLO, N, NRHS, A, LDA, X, LDA, WORK,
+     $                            LDA, RWORK, RESULT( 4 ) )
+*
+*+    TEST 5
 *                 Check solution from generated exact solution.
 *
                      CALL CGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
-     $                            RESULT( 4 ) )
+     $                            RESULT( 5 ) )
 *
-*+    TESTS 5, 6, and 7
+*+    TESTS 6, 7, and 8
 *                 Use iterative refinement to improve the solution.
 *
                      SRNAMT = 'CHERFS'
      $                               NERRS, NOUT )
 *
                      CALL CGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
-     $                            RESULT( 5 ) )
+     $                            RESULT( 6 ) )
                      CALL CPOT05( UPLO, N, NRHS, A, LDA, B, LDA, X, LDA,
      $                            XACT, LDA, RWORK, RWORK( NRHS+1 ),
-     $                            RESULT( 6 ) )
+     $                            RESULT( 7 ) )
 *
 *                    Print information about the tests that did not pass
 *                    the threshold.
 *
-                     DO 120 K = 3, 7
+                     DO 120 K = 3, 8
                         IF( RESULT( K ).GE.THRESH ) THEN
                            IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
      $                        CALL ALAHD( NOUT, PATH )
                      NRUN = NRUN + 5
   130             CONTINUE
 *
-*+    TEST 8
+*+    TEST 9
 *                 Get an estimate of RCOND = 1/CNDNUM.
 *
   140             CONTINUE
      $               CALL ALAERH( PATH, 'CHECON', INFO, 0, UPLO, N, N,
      $                            -1, -1, -1, IMAT, NFAIL, NERRS, NOUT )
 *
-                  RESULT( 8 ) = SGET06( RCOND, RCONDC )
+                  RESULT( 9 ) = SGET06( RCOND, RCONDC )
 *
 *                 Print information about the tests that did not pass
 *                 the threshold.
 *
-                  IF( RESULT( 8 ).GE.THRESH ) THEN
+                  IF( RESULT( 9 ).GE.THRESH ) THEN
                      IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
      $                  CALL ALAHD( NOUT, PATH )
                      WRITE( NOUT, FMT = 9997 )UPLO, N, IMAT, 8,
-     $                  RESULT( 8 )
+     $                  RESULT( 9 )
                      NFAIL = NFAIL + 1
                   END IF
                   NRUN = NRUN + 1
index 151b4f7..e6524ca 100644 (file)
@@ -22,7 +22,7 @@
 *  Purpose
 *  =======
 *
-*  ZCHKHE tests ZHETRF, -TRI, -TRS, -RFS, and -CON.
+*  ZCHKHE tests ZHETRF, -TRI, -TRS, -TRS2, -RFS, and -CON.
 *
 *  Arguments
 *  =========
@@ -94,7 +94,7 @@
       INTEGER            NTYPES
       PARAMETER          ( NTYPES = 10 )
       INTEGER            NTESTS
-      PARAMETER          ( NTESTS = 8 )
+      PARAMETER          ( NTESTS = 9 )
 *     ..
 *     .. Local Scalars ..
       LOGICAL            TRFCON, ZEROT
      $                            LDA, RWORK, RESULT( 3 ) )
 *
 *+    TEST 4
+*                 Solve and compute residual for  A * X = B.
+*
+                     SRNAMT = 'ZLARHS'
+                     CALL ZLARHS( PATH, 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 = 'ZHETRS2'
+                     CALL ZSYCONV( UPLO, 'C', N, AFAC, LDA, IWORK, WORK, 
+     $                            INFO )
+                     CALL ZHETRS2( UPLO, N, NRHS, AFAC, LDA, IWORK, X,
+     $                            LDA, WORK, INFO )
+                     CALL ZSYCONV( UPLO, 'R', N, AFAC, LDA, IWORK, WORK,
+     $                            INFO )
+*
+*                 Check error code from ZSYTRS2.
+*
+                     IF( INFO.NE.0 )
+     $                  CALL ALAERH( PATH, 'ZHETRS2', INFO, 0, UPLO, N,
+     $                               N, -1, -1, NRHS, IMAT, NFAIL,
+     $                               NERRS, NOUT )
+*
+                     CALL ZLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA )
+                     CALL ZPOT02( UPLO, N, NRHS, A, LDA, X, LDA, WORK,
+     $                            LDA, RWORK, RESULT( 4 ) )
+*
+*+    TEST 5
 *                 Check solution from generated exact solution.
 *
                      CALL ZGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
-     $                            RESULT( 4 ) )
+     $                            RESULT( 5 ) )
 *
-*+    TESTS 5, 6, and 7
+*+    TESTS 6, 7, and 8
 *                 Use iterative refinement to improve the solution.
 *
                      SRNAMT = 'ZHERFS'
      $                               NERRS, NOUT )
 *
                      CALL ZGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
-     $                            RESULT( 5 ) )
+     $                            RESULT( 6 ) )
                      CALL ZPOT05( UPLO, N, NRHS, A, LDA, B, LDA, X, LDA,
      $                            XACT, LDA, RWORK, RWORK( NRHS+1 ),
-     $                            RESULT( 6 ) )
+     $                            RESULT( 7 ) )
 *
 *                    Print information about the tests that did not pass
 *                    the threshold.
 *
-                     DO 120 K = 3, 7
+                     DO 120 K = 3, 8
                         IF( RESULT( K ).GE.THRESH ) THEN
                            IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
      $                        CALL ALAHD( NOUT, PATH )
                      NRUN = NRUN + 5
   130             CONTINUE
 *
-*+    TEST 8
+*+    TEST 9
 *                 Get an estimate of RCOND = 1/CNDNUM.
 *
   140             CONTINUE
      $               CALL ALAERH( PATH, 'ZHECON', INFO, 0, UPLO, N, N,
      $                            -1, -1, -1, IMAT, NFAIL, NERRS, NOUT )
 *
-                  RESULT( 8 ) = DGET06( RCOND, RCONDC )
+                  RESULT( 9 ) = DGET06( RCOND, RCONDC )
 *
 *                 Print information about the tests that did not pass
 *                 the threshold.
 *
-                  IF( RESULT( 8 ).GE.THRESH ) THEN
+                  IF( RESULT( 9 ).GE.THRESH ) THEN
                      IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
      $                  CALL ALAHD( NOUT, PATH )
-                     WRITE( NOUT, FMT = 9997 )UPLO, N, IMAT, 8,
-     $                  RESULT( 8 )
+                     WRITE( NOUT, FMT = 9997 )UPLO, N, IMAT, 9,
+     $                  RESULT( 9 )
                      NFAIL = NFAIL + 1
                   END IF
                   NRUN = NRUN + 1