diff options
| -rw-r--r-- | SRC/Makefile | 8 | ||||
| -rw-r--r-- | SRC/csyconv.f | 302 | ||||
| -rw-r--r-- | SRC/csysv.f | 12 | ||||
| -rw-r--r-- | SRC/csytrs2.f | 272 | ||||
| -rw-r--r-- | SRC/dsyconv.f | 302 | ||||
| -rw-r--r-- | SRC/dsysv.f | 17 | ||||
| -rw-r--r-- | SRC/dsytrs2.f | 272 | ||||
| -rw-r--r-- | SRC/ssyconv.f | 302 | ||||
| -rw-r--r-- | SRC/ssysv.f | 17 | ||||
| -rw-r--r-- | SRC/ssytrs2.f | 272 | ||||
| -rw-r--r-- | SRC/zsyconv.f | 302 | ||||
| -rw-r--r-- | SRC/zsysv.f | 15 | ||||
| -rw-r--r-- | SRC/zsytrs2.f | 272 | ||||
| -rw-r--r-- | TESTING/LIN/alahd.f | 28 | ||||
| -rw-r--r-- | TESTING/LIN/cchksy.f | 52 | ||||
| -rw-r--r-- | TESTING/LIN/dchksy.f | 61 | ||||
| -rw-r--r-- | TESTING/LIN/schksy.f | 60 | ||||
| -rw-r--r-- | TESTING/LIN/zchksy.f | 59 |
18 files changed, 2548 insertions, 77 deletions
diff --git a/SRC/Makefile b/SRC/Makefile index ce43806d..67315be6 100644 --- a/SRC/Makefile +++ b/SRC/Makefile @@ -131,7 +131,7 @@ SLASRC = \ ssptrf.o ssptri.o ssptrs.o sstegr.o sstein.o sstev.o sstevd.o sstevr.o \ sstevx.o ssycon.o ssyev.o ssyevd.o ssyevr.o ssyevx.o ssygs2.o \ ssygst.o ssygv.o ssygvd.o ssygvx.o ssyrfs.o ssysv.o ssysvx.o \ - ssytd2.o ssytf2.o ssytrd.o ssytrf.o ssytri.o ssytrs.o stbcon.o \ + ssytd2.o ssytf2.o ssytrd.o ssytrf.o ssytri.o ssytrs.o ssytrs2.o ssyconv.o stbcon.o \ stbrfs.o stbtrs.o stgevc.o stgex2.o stgexc.o stgsen.o \ stgsja.o stgsna.o stgsy2.o stgsyl.o stpcon.o stprfs.o stptri.o \ stptrs.o \ @@ -195,7 +195,7 @@ CLASRC = \ cspsvx.o csptrf.o csptri.o csptrs.o csrscl.o cstedc.o \ cstegr.o cstein.o csteqr.o csycon.o csymv.o \ csyr.o csyrfs.o csysv.o csysvx.o csytf2.o csytrf.o csytri.o \ - csytrs.o ctbcon.o ctbrfs.o ctbtrs.o ctgevc.o ctgex2.o \ + csytrs.o csytrs2.o csyconv.o ctbcon.o ctbrfs.o ctbtrs.o ctgevc.o ctgex2.o \ ctgexc.o ctgsen.o ctgsja.o ctgsna.o ctgsy2.o ctgsyl.o ctpcon.o \ ctprfs.o ctptri.o \ ctptrs.o ctrcon.o ctrevc.o ctrexc.o ctrrfs.o ctrsen.o ctrsna.o \ @@ -264,7 +264,7 @@ DLASRC = \ dstevx.o dsycon.o dsyev.o dsyevd.o dsyevr.o \ dsyevx.o dsygs2.o dsygst.o dsygv.o dsygvd.o dsygvx.o dsyrfs.o \ dsysv.o dsysvx.o \ - dsytd2.o dsytf2.o dsytrd.o dsytrf.o dsytri.o dsytrs.o dtbcon.o \ + dsytd2.o dsytf2.o dsytrd.o dsytrf.o dsytri.o dsytrs.o dsytrs2.o dsyconv.o dtbcon.o \ dtbrfs.o dtbtrs.o dtgevc.o dtgex2.o dtgexc.o dtgsen.o \ dtgsja.o dtgsna.o dtgsy2.o dtgsyl.o dtpcon.o dtprfs.o dtptri.o \ dtptrs.o \ @@ -332,7 +332,7 @@ ZLASRC = \ zspsvx.o zsptrf.o zsptri.o zsptrs.o zdrscl.o zstedc.o \ zstegr.o zstein.o zsteqr.o zsycon.o zsymv.o \ zsyr.o zsyrfs.o zsysv.o zsysvx.o zsytf2.o zsytrf.o zsytri.o \ - zsytrs.o ztbcon.o ztbrfs.o ztbtrs.o ztgevc.o ztgex2.o \ + zsytrs.o zsytrs2.o zsyconv.o ztbcon.o ztbrfs.o ztbtrs.o ztgevc.o ztgex2.o \ ztgexc.o ztgsen.o ztgsja.o ztgsna.o ztgsy2.o ztgsyl.o ztpcon.o \ ztprfs.o ztptri.o \ ztptrs.o ztrcon.o ztrevc.o ztrexc.o ztrrfs.o ztrsen.o ztrsna.o \ diff --git a/SRC/csyconv.f b/SRC/csyconv.f new file mode 100644 index 00000000..74306022 --- /dev/null +++ b/SRC/csyconv.f @@ -0,0 +1,302 @@ + SUBROUTINE CSYCONV( UPLO, WAY, N, A, LDA, IPIV, WORK, INFO ) +* +* -- LAPACK PROTOTYPE routine (version 3.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, WAY + INTEGER INFO, LDA, N +* .. +* .. Array Arguments .. + INTEGER IPIV( * ) + COMPLEX A( LDA, * ), WORK( * ) +* .. +* +* Purpose +* ======= +* +* CSYCONV convert A given by TRF into L and D and vice-versa. +* Get Non-diag elements of D (returned in workspace) and +* apply or reverse permutation done in TRF. +* +* 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**T; +* = 'L': Lower triangular, form is A = L*D*L**T. +* +* WAY (input) CHARACTER*1 +* = 'C': Convert +* = 'R': Revert +* +* N (input) INTEGER +* The order of the matrix A. N >= 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 CSYTRF. +* +* 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 CSYTRF. +* +* WORK (workspace) COMPLEX array, dimension (N) +* +* LWORK (input) INTEGER +* The length of WORK. LWORK >=1. +* LWORK = N +* +* If LWORK = -1, then a workspace query is assumed; the routine +* only calculates the optimal size of the WORK array, returns +* this value as the first entry of the WORK array, and no error +* message related to LWORK is issued by XERBLA. +* +* INFO (output) INTEGER +* = 0: successful exit +* < 0: if INFO = -i, the i-th argument had an illegal value +* +* ===================================================================== +* +* .. Parameters .. + COMPLEX ZERO + PARAMETER ( ZERO = (0.0E+0,0.0E+0) ) +* .. +* .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME +* +* .. External Subroutines .. + EXTERNAL XERBLA +* .. Local Scalars .. + LOGICAL UPPER, CONVERT + INTEGER I, IP + COMPLEX TEMP +* .. +* .. Executable Statements .. +* + INFO = 0 + UPPER = LSAME( UPLO, 'U' ) + CONVERT = LSAME( WAY, 'C' ) + IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN + INFO = -1 + ELSE IF( .NOT.CONVERT .AND. .NOT.LSAME( WAY, 'R' ) ) THEN + INFO = -2 + ELSE IF( N.LT.0 ) THEN + INFO = -3 + ELSE IF( LDA.LT.MAX( 1, N ) ) THEN + INFO = -5 + + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'CSYCONV', -INFO ) + RETURN + END IF +* +* Quick return if possible +* + IF( N.EQ.0 ) + $ RETURN +* + IF( UPPER ) THEN +* +* A is UPPER +* +* Convert A (A is upper) +* +* Convert VALUE +* + IF ( CONVERT ) THEN + I=N + WORK(1)=ZERO + DO WHILE ( I .GT. 1 ) + IF( IPIV(I) .LT. 0 ) THEN + WORK(I)=A(I-1,I) + A(I-1,I)=ZERO + I=I-1 + ELSE + WORK(I)=ZERO + ENDIF + I=I-1 + END DO +* +* Convert PERMUTATIONS +* + I=N + DO WHILE ( I .GE. 1 ) + IF( IPIV(I) .GT. 0) THEN + IP=IPIV(I) + IF( I .LT. N) THEN + DO 12 J= I+1,N + TEMP=A(IP,J) + A(IP,J)=A(I,J) + A(I,J)=TEMP + 12 CONTINUE + ENDIF + ELSE + IP=-IPIV(I) + IF( I .LT. N) THEN + DO 13 J= I+1,N + TEMP=A(IP,J) + A(IP,J)=A(I-1,J) + A(I-1,J)=TEMP + 13 CONTINUE + ENDIF + I=I-1 + ENDIF + I=I-1 + END DO + + ELSE +* +* Revert A (A is upper) +* +* +* Revert PERMUTATIONS +* + I=1 + DO WHILE ( I .LE. N ) + IF( IPIV(I) .GT. 0 ) THEN + IP=IPIV(I) + IF( I .LT. N) THEN + DO J= I+1,N + TEMP=A(IP,J) + A(IP,J)=A(I,J) + A(I,J)=TEMP + END DO + ENDIF + ELSE + IP=-IPIV(I) + I=I+1 + IF( I .LT. N) THEN + DO J= I+1,N + TEMP=A(IP,J) + A(IP,J)=A(I-1,J) + A(I-1,J)=TEMP + END DO + ENDIF + ENDIF + I=I+1 + END DO +* +* Revert VALUE +* + I=N + DO WHILE ( I .GT. 1 ) + IF( IPIV(I) .LT. 0 ) THEN + A(I-1,I)=WORK(I) + I=I-1 + ENDIF + I=I-1 + END DO + END IF + ELSE +* +* A is LOWER +* + IF ( CONVERT ) THEN +* +* Convert A (A is lower) +* +* +* Convert VALUE +* + I=1 + WORK(N)=ZERO + DO WHILE ( I .LE. N ) + IF( I.LT.N .AND. IPIV(I) .LT. 0 ) THEN + WORK(I)=A(I+1,I) + A(I+1,I)=ZERO + I=I+1 + ELSE + WORK(I)=ZERO + ENDIF + I=I+1 + END DO +* +* Convert PERMUTATIONS +* + I=1 + DO WHILE ( I .LE. N ) + IF( IPIV(I) .GT. 0 ) THEN + IP=IPIV(I) + IF (I .GT. 1) THEN + DO 22 J= 1,I-1 + TEMP=A(IP,J) + A(IP,J)=A(I,J) + A(I,J)=TEMP + 22 CONTINUE + ENDIF + ELSE + IP=-IPIV(I) + IF ( I .GT. 1 ) THEN + DO 23 J= 1,I-1 + TEMP=A(IP,J) + A(IP,J)=A(I+1,J) + A(I+1,J)=TEMP + 23 CONTINUE + ENDIF + I=I+1 + ENDIF + I=I+1 + END DO + ELSE +* +* Revert A (A is lower) +* +* +* Revert PERMUTATIONS +* + I=N + DO WHILE ( I .GE. 1 ) + IF( IPIV(I) .GT. 0 ) THEN + IP=IPIV(I) + IF (I .GT. 1) THEN + DO J= 1,I-1 + TEMP=A(I,J) + A(I,J)=A(IP,J) + A(IP,J)=TEMP + END DO + ENDIF + ELSE + IP=-IPIV(I) + I=I-1 + IF (I .GT. 1) THEN + DO J= 1,I-1 + TEMP=A(I+1,J) + A(I+1,J)=A(IP,J) + A(IP,J)=TEMP + END DO + ENDIF + ENDIF + I=I-1 + END DO +* +* Revert VALUE +* + I=1 + DO WHILE ( I .LE. N-1 ) + IF( IPIV(I) .LT. 0 ) THEN + A(I+1,I)=WORK(I) + I=I+1 + ENDIF + I=I+1 + END DO + END IF + END IF + + RETURN +* +* End of CSYCONV +* + END diff --git a/SRC/csysv.f b/SRC/csysv.f index 855405f2..a9832e53 100644 --- a/SRC/csysv.f +++ b/SRC/csysv.f @@ -113,7 +113,7 @@ EXTERNAL ILAENV, LSAME * .. * .. External Subroutines .. - EXTERNAL CSYTRF, CSYTRS, XERBLA + EXTERNAL CSYCONV, CSYTRF, CSYTRS2, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX @@ -160,9 +160,17 @@ CALL CSYTRF( 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 CSYTRS( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, INFO ) + CALL CSYTRS2( 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/csytrs2.f b/SRC/csytrs2.f new file mode 100644 index 00000000..bd670c1b --- /dev/null +++ b/SRC/csytrs2.f @@ -0,0 +1,272 @@ + SUBROUTINE CSYTRS2( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, + $ WORK, INFO ) +* +* -- LAPACK PROTOTYPE routine (version 3.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 +* ======= +* +* CSYTRS2 solves a system of linear equations A*X = B with a COMPLEX +* symmetric 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**T; +* = 'L': Lower triangular, form is A = L*D*L**T. +* +* 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 CSYTRF. +* +* 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 CSYTRF. +* +* 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 + COMPLEX AK, AKM1, AKM1K, BK, BKM1, DENOM +* .. +* .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME +* .. +* .. External Subroutines .. + EXTERNAL CSCAL, CSWAP, CTRSM, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC 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( 'CSYTRS2', -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 + CALL CSCAL( NRHS, ONE / A( I, I ), B( I, 1 ), N ) + 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 ) / AKM1K + DENOM = AKM1*AK - ONE + DO 15 J = 1, NRHS + BKM1 = B( I-1, J ) / AKM1K + BK = B( I, J ) / 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','T','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 + CALL CSCAL( NRHS, ONE / A( I, I ), B( I, 1 ), N ) + ELSE + AKM1K = WORK(I) + AKM1 = A( I, I ) / AKM1K + AK = A( I+1, I+1 ) / AKM1K + DENOM = AKM1*AK - ONE + DO 25 J = 1, NRHS + BKM1 = B( I, J ) / 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','T','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 CSYTRS2 +* + END diff --git a/SRC/dsyconv.f b/SRC/dsyconv.f new file mode 100644 index 00000000..0666098a --- /dev/null +++ b/SRC/dsyconv.f @@ -0,0 +1,302 @@ + SUBROUTINE DSYCONV( UPLO, WAY, N, A, LDA, IPIV, WORK, INFO ) +* +* -- LAPACK PROTOTYPE routine (version 3.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, WAY + INTEGER INFO, LDA, N +* .. +* .. Array Arguments .. + INTEGER IPIV( * ) + DOUBLE PRECISION A( LDA, * ), WORK( * ) +* .. +* +* Purpose +* ======= +* +* DSYCONV convert A given by TRF into L and D and vice-versa. +* Get Non-diag elements of D (returned in workspace) and +* apply or reverse permutation done in TRF. +* +* 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**T; +* = 'L': Lower triangular, form is A = L*D*L**T. +* +* WAY (input) CHARACTER*1 +* = 'C': Convert +* = 'R': Revert +* +* N (input) INTEGER +* The order of the matrix A. N >= 0. +* +* A (input) DOUBLE PRECISION array, dimension (LDA,N) +* The block diagonal matrix D and the multipliers used to +* obtain the factor U or L as computed by DSYTRF. +* +* 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 DSYTRF. +* +* WORK (workspace) DOUBLE PRECISION array, dimension (N) +* +* LWORK (input) INTEGER +* The length of WORK. LWORK >=1. +* LWORK = N +* +* If LWORK = -1, then a workspace query is assumed; the routine +* only calculates the optimal size of the WORK array, returns +* this value as the first entry of the WORK array, and no error +* message related to LWORK is issued by XERBLA. +* +* INFO (output) INTEGER +* = 0: successful exit +* < 0: if INFO = -i, the i-th argument had an illegal value +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ZERO + PARAMETER ( ZERO = 0.0D+0 ) +* .. +* .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME +* +* .. External Subroutines .. + EXTERNAL XERBLA +* .. Local Scalars .. + LOGICAL UPPER, CONVERT + INTEGER I, IP + DOUBLE PRECISION TEMP +* .. +* .. Executable Statements .. +* + INFO = 0 + UPPER = LSAME( UPLO, 'U' ) + CONVERT = LSAME( WAY, 'C' ) + IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN + INFO = -1 + ELSE IF( .NOT.CONVERT .AND. .NOT.LSAME( WAY, 'R' ) ) THEN + INFO = -2 + ELSE IF( N.LT.0 ) THEN + INFO = -3 + ELSE IF( LDA.LT.MAX( 1, N ) ) THEN + INFO = -5 + + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'DSYCONV', -INFO ) + RETURN + END IF +* +* Quick return if possible +* + IF( N.EQ.0 ) + $ RETURN +* + IF( UPPER ) THEN +* +* A is UPPER +* +* Convert A (A is upper) +* +* Convert VALUE +* + IF ( CONVERT ) THEN + I=N + WORK(1)=ZERO + DO WHILE ( I .GT. 1 ) + IF( IPIV(I) .LT. 0 ) THEN + WORK(I)=A(I-1,I) + A(I-1,I)=ZERO + I=I-1 + ELSE + WORK(I)=ZERO + ENDIF + I=I-1 + END DO +* +* Convert PERMUTATIONS +* + I=N + DO WHILE ( I .GE. 1 ) + IF( IPIV(I) .GT. 0) THEN + IP=IPIV(I) + IF( I .LT. N) THEN + DO 12 J= I+1,N + TEMP=A(IP,J) + A(IP,J)=A(I,J) + A(I,J)=TEMP + 12 CONTINUE + ENDIF + ELSE + IP=-IPIV(I) + IF( I .LT. N) THEN + DO 13 J= I+1,N + TEMP=A(IP,J) + A(IP,J)=A(I-1,J) + A(I-1,J)=TEMP + 13 CONTINUE + ENDIF + I=I-1 + ENDIF + I=I-1 + END DO + + ELSE +* +* Revert A (A is upper) +* +* +* Revert PERMUTATIONS +* + I=1 + DO WHILE ( I .LE. N ) + IF( IPIV(I) .GT. 0 ) THEN + IP=IPIV(I) + IF( I .LT. N) THEN + DO J= I+1,N + TEMP=A(IP,J) + A(IP,J)=A(I,J) + A(I,J)=TEMP + END DO + ENDIF + ELSE + IP=-IPIV(I) + I=I+1 + IF( I .LT. N) THEN + DO J= I+1,N + TEMP=A(IP,J) + A(IP,J)=A(I-1,J) + A(I-1,J)=TEMP + END DO + ENDIF + ENDIF + I=I+1 + END DO +* +* Revert VALUE +* + I=N + DO WHILE ( I .GT. 1 ) + IF( IPIV(I) .LT. 0 ) THEN + A(I-1,I)=WORK(I) + I=I-1 + ENDIF + I=I-1 + END DO + END IF + ELSE +* +* A is LOWER +* + IF ( CONVERT ) THEN +* +* Convert A (A is lower) +* +* +* Convert VALUE +* + I=1 + WORK(N)=ZERO + DO WHILE ( I .LE. N ) + IF( I.LT.N .AND. IPIV(I) .LT. 0 ) THEN + WORK(I)=A(I+1,I) + A(I+1,I)=ZERO + I=I+1 + ELSE + WORK(I)=ZERO + ENDIF + I=I+1 + END DO +* +* Convert PERMUTATIONS +* + I=1 + DO WHILE ( I .LE. N ) + IF( IPIV(I) .GT. 0 ) THEN + IP=IPIV(I) + IF (I .GT. 1) THEN + DO 22 J= 1,I-1 + TEMP=A(IP,J) + A(IP,J)=A(I,J) + A(I,J)=TEMP + 22 CONTINUE + ENDIF + ELSE + IP=-IPIV(I) + IF (I .GT. 1) THEN + DO 23 J= 1,I-1 + TEMP=A(IP,J) + A(IP,J)=A(I+1,J) + A(I+1,J)=TEMP + 23 CONTINUE + ENDIF + I=I+1 + ENDIF + I=I+1 + END DO + ELSE +* +* Revert A (A is lower) +* +* +* Revert PERMUTATIONS +* + I=N + DO WHILE ( I .GE. 1 ) + IF( IPIV(I) .GT. 0 ) THEN + IP=IPIV(I) + IF (I .GT. 1) THEN + DO J= 1,I-1 + TEMP=A(I,J) + A(I,J)=A(IP,J) + A(IP,J)=TEMP + END DO + ENDIF + ELSE + IP=-IPIV(I) + I=I-1 + IF (I .GT. 1) THEN + DO J= 1,I-1 + TEMP=A(I+1,J) + A(I+1,J)=A(IP,J) + A(IP,J)=TEMP + END DO + ENDIF + ENDIF + I=I-1 + END DO +* +* Revert VALUE +* + I=1 + DO WHILE ( I .LE. N-1 ) + IF( IPIV(I) .LT. ZERO ) THEN + A(I+1,I)=WORK(I) + I=I+1 + ENDIF + I=I+1 + END DO + END IF + END IF + + RETURN +* +* End of DSYCONV +* + END diff --git a/SRC/dsysv.f b/SRC/dsysv.f index 3abe723e..f8dbac5b 100644 --- a/SRC/dsysv.f +++ b/SRC/dsysv.f @@ -4,7 +4,7 @@ * -- LAPACK driver routine (version 3.2) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* November 2006 +* May 2010 * * .. Scalar Arguments .. CHARACTER UPLO @@ -105,7 +105,7 @@ * * .. Local Scalars .. LOGICAL LQUERY - INTEGER LWKOPT, NB + INTEGER LWKOPT, NB, IINFO * .. * .. External Functions .. LOGICAL LSAME @@ -113,7 +113,7 @@ EXTERNAL LSAME, ILAENV * .. * .. External Subroutines .. - EXTERNAL DSYTRF, DSYTRS, XERBLA + EXTERNAL DSYCONV, DSYTRF, DSYTRS2, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX @@ -123,6 +123,7 @@ * Test the input parameters. * INFO = 0 + IINFO = 0 LQUERY = ( LWORK.EQ.-1 ) IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN INFO = -1 @@ -160,9 +161,17 @@ CALL DSYTRF( UPLO, N, A, LDA, IPIV, WORK, LWORK, INFO ) IF( INFO.EQ.0 ) THEN * +* Convert A +* + CALL DSYCONV( UPLO, 'C', N, A, LDA, IPIV, WORK, IINFO ) +* * Solve the system A*X = B, overwriting B with X. * - CALL DSYTRS( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, INFO ) + CALL DSYTRS2( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, WORK, INFO ) +* +* Revert A +* + CALL DSYCONV( UPLO, 'R', N, A, LDA, IPIV, WORK, IINFO ) * END IF * diff --git a/SRC/dsytrs2.f b/SRC/dsytrs2.f new file mode 100644 index 00000000..658aed24 --- /dev/null +++ b/SRC/dsytrs2.f @@ -0,0 +1,272 @@ + SUBROUTINE DSYTRS2( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, + $ WORK, INFO ) +* +* -- LAPACK PROTOTYPE routine (version 3.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 PRECISION A( LDA, * ), B( LDB, * ), WORK( * ) +* .. +* +* Purpose +* ======= +* +* DSYTRS2 solves a system of linear equations A*X = B with a real +* symmetric matrix A using the factorization A = U*D*U**T or +* A = L*D*L**T computed by DSYTRF and converted by DSYCONV. +* +* 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**T; +* = 'L': Lower triangular, form is A = L*D*L**T. +* +* 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 PRECISION array, dimension (LDA,N) +* The block diagonal matrix D and the multipliers used to +* obtain the factor U or L as computed by DSYTRF. +* +* 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 DSYTRF. +* +* B (input/output) DOUBLE PRECISION 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 PRECISION ONE + PARAMETER ( ONE = 1.0D+0 ) +* .. +* .. Local Scalars .. + LOGICAL UPPER + INTEGER I, J, K, KP + DOUBLE PRECISION AK, AKM1, AKM1K, BK, BKM1, DENOM +* .. +* .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME +* .. +* .. External Subroutines .. + EXTERNAL DSCAL, DSWAP, DTRSM, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC 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( 'DSYTRS2', -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 DSWAP( 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 DSWAP( 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 DTRSM('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 + CALL DSCAL( NRHS, ONE / A( I, I ), B( I, 1 ), N ) + 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 ) / AKM1K + DENOM = AKM1*AK - ONE + DO 15 J = 1, NRHS + BKM1 = B( I-1, J ) / AKM1K + BK = B( I, J ) / 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 DTRSM('L','U','T','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 DSWAP( 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 DSWAP( 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 DSWAP( 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 DSWAP( 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 DTRSM('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 + CALL DSCAL( NRHS, ONE / A( I, I ), B( I, 1 ), N ) + ELSE + AKM1K = WORK(I) + AKM1 = A( I, I ) / AKM1K + AK = A( I+1, I+1 ) / AKM1K + DENOM = AKM1*AK - ONE + DO 25 J = 1, NRHS + BKM1 = B( I, J ) / 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 DTRSM('L','L','T','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 DSWAP( 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 DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) + K=K-2 + ENDIF + END DO +* + END IF +* + RETURN +* +* End of DSYTRS2 +* + END diff --git a/SRC/ssyconv.f b/SRC/ssyconv.f new file mode 100644 index 00000000..5a966225 --- /dev/null +++ b/SRC/ssyconv.f @@ -0,0 +1,302 @@ + SUBROUTINE SSYCONV( UPLO, WAY, N, A, LDA, IPIV, WORK, INFO ) +* +* -- LAPACK PROTOTYPE routine (version 3.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, WAY + INTEGER INFO, LDA, N +* .. +* .. Array Arguments .. + INTEGER IPIV( * ) + REAL A( LDA, * ), WORK( * ) +* .. +* +* Purpose +* ======= +* +* SSYCONV convert A given by TRF into L and D and vice-versa. +* Get Non-diag elements of D (returned in workspace) and +* apply or reverse permutation done in TRF. +* +* 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**T; +* = 'L': Lower triangular, form is A = L*D*L**T. +* +* WAY (input) CHARACTER*1 +* = 'C': Convert +* = 'R': Revert +* +* N (input) INTEGER +* The order of the matrix A. N >= 0. +* +* A (input) REAL array, dimension (LDA,N) +* The block diagonal matrix D and the multipliers used to +* obtain the factor U or L as computed by SSYTRF. +* +* 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 SSYTRF. +* +* WORK (workspace) REAL array, dimension (N) +* +* LWORK (input) INTEGER +* The length of WORK. LWORK >=1. +* LWORK = N +* +* If LWORK = -1, then a workspace query is assumed; the routine +* only calculates the optimal size of the WORK array, returns +* this value as the first entry of the WORK array, and no error +* message related to LWORK is issued by XERBLA. +* +* INFO (output) INTEGER +* = 0: successful exit +* < 0: if INFO = -i, the i-th argument had an illegal value +* +* ===================================================================== +* +* .. Parameters .. + REAL ZERO + PARAMETER ( ZERO = 0.0E+0 ) +* .. +* .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME +* +* .. External Subroutines .. + EXTERNAL XERBLA +* .. Local Scalars .. + LOGICAL UPPER, CONVERT + INTEGER I, IP + REAL TEMP +* .. +* .. Executable Statements .. +* + INFO = 0 + UPPER = LSAME( UPLO, 'U' ) + CONVERT = LSAME( WAY, 'C' ) + IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN + INFO = -1 + ELSE IF( .NOT.CONVERT .AND. .NOT.LSAME( WAY, 'R' ) ) THEN + INFO = -2 + ELSE IF( N.LT.0 ) THEN + INFO = -3 + ELSE IF( LDA.LT.MAX( 1, N ) ) THEN + INFO = -5 + + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'SSYCONV', -INFO ) + RETURN + END IF +* +* Quick return if possible +* + IF( N.EQ.0 ) + $ RETURN +* + IF( UPPER ) THEN +* +* A is UPPER +* +* Convert A (A is upper) +* +* Convert VALUE +* + IF ( CONVERT ) THEN + I=N + WORK(1)=ZERO + DO WHILE ( I .GT. 1 ) + IF( IPIV(I) .LT. 0 ) THEN + WORK(I)=A(I-1,I) + A(I-1,I)=ZERO + I=I-1 + ELSE + WORK(I)=ZERO + ENDIF + I=I-1 + END DO +* +* Convert PERMUTATIONS +* + I=N + DO WHILE ( I .GE. 1 ) + IF( IPIV(I) .GT. 0) THEN + IP=IPIV(I) + IF( I .LT. N) THEN + DO 12 J= I+1,N + TEMP=A(IP,J) + A(IP,J)=A(I,J) + A(I,J)=TEMP + 12 CONTINUE + ENDIF + ELSE + IP=-IPIV(I) + IF( I .LT. N) THEN + DO 13 J= I+1,N + TEMP=A(IP,J) + A(IP,J)=A(I-1,J) + A(I-1,J)=TEMP + 13 CONTINUE + ENDIF + I=I-1 + ENDIF + I=I-1 + END DO + + ELSE +* +* Revert A (A is upper) +* +* +* Revert PERMUTATIONS +* + I=1 + DO WHILE ( I .LE. N ) + IF( IPIV(I) .GT. 0 ) THEN + IP=IPIV(I) + IF( I .LT. N) THEN + DO J= I+1,N + TEMP=A(IP,J) + A(IP,J)=A(I,J) + A(I,J)=TEMP + END DO + ENDIF + ELSE + IP=-IPIV(I) + I=I+1 + IF( I .LT. N) THEN + DO J= I+1,N + TEMP=A(IP,J) + A(IP,J)=A(I-1,J) + A(I-1,J)=TEMP + END DO + ENDIF + ENDIF + I=I+1 + END DO +* +* Revert VALUE +* + I=N + DO WHILE ( I .GT. 1 ) + IF( IPIV(I) .LT. 0 ) THEN + A(I-1,I)=WORK(I) + I=I-1 + ENDIF + I=I-1 + END DO + END IF + ELSE +* +* A is LOWER +* + IF ( CONVERT ) THEN +* +* Convert A (A is lower) +* +* +* Convert VALUE +* + I=1 + WORK(N)=ZERO + DO WHILE ( I .LE. N ) + IF( I.LT.N .AND. IPIV(I) .LT. 0 ) THEN + WORK(I)=A(I+1,I) + A(I+1,I)=ZERO + I=I+1 + ELSE + WORK(I)=ZERO + ENDIF + I=I+1 + END DO +* +* Convert PERMUTATIONS +* + I=1 + DO WHILE ( I .LE. N ) + IF( IPIV(I) .GT. 0 ) THEN + IP=IPIV(I) + IF (I .GT. 1) THEN + DO 22 J= 1,I-1 + TEMP=A(IP,J) + A(IP,J)=A(I,J) + A(I,J)=TEMP + 22 CONTINUE + ENDIF + ELSE + IP=-IPIV(I) + IF (I .GT. 1) THEN + DO 23 J= 1,I-1 + TEMP=A(IP,J) + A(IP,J)=A(I+1,J) + A(I+1,J)=TEMP + 23 CONTINUE + ENDIF + I=I+1 + ENDIF + I=I+1 + END DO + ELSE +* +* Revert A (A is lower) +* +* +* Revert PERMUTATIONS +* + I=N + DO WHILE ( I .GE. 1 ) + IF( IPIV(I) .GT. 0 ) THEN + IP=IPIV(I) + IF (I .GT. 1) THEN + DO J= 1,I-1 + TEMP=A(I,J) + A(I,J)=A(IP,J) + A(IP,J)=TEMP + END DO + ENDIF + ELSE + IP=-IPIV(I) + I=I-1 + IF (I .GT. 1) THEN + DO J= 1,I-1 + TEMP=A(I+1,J) + A(I+1,J)=A(IP,J) + A(IP,J)=TEMP + END DO + ENDIF + ENDIF + I=I-1 + END DO +* +* Revert VALUE +* + I=1 + DO WHILE ( I .LE. N-1 ) + IF( IPIV(I) .LT. 0 ) THEN + A(I+1,I)=WORK(I) + I=I+1 + ENDIF + I=I+1 + END DO + END IF + END IF + + RETURN +* +* End of SSYCONV +* + END diff --git a/SRC/ssysv.f b/SRC/ssysv.f index bcefac89..a7262b65 100644 --- a/SRC/ssysv.f +++ b/SRC/ssysv.f @@ -4,7 +4,7 @@ * -- LAPACK driver routine (version 3.2) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* November 2006 +* May 2010 * * .. Scalar Arguments .. CHARACTER UPLO @@ -105,7 +105,7 @@ * * .. Local Scalars .. LOGICAL LQUERY - INTEGER LWKOPT, NB + INTEGER LWKOPT, NB, IINFO * .. * .. External Functions .. LOGICAL LSAME @@ -113,7 +113,7 @@ EXTERNAL ILAENV, LSAME * .. * .. External Subroutines .. - EXTERNAL SSYTRF, SSYTRS, XERBLA + EXTERNAL SSYCONV, SSYTRF, SSYTRS2, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX @@ -123,6 +123,7 @@ * Test the input parameters. * INFO = 0 + IINFO = 0 LQUERY = ( LWORK.EQ.-1 ) IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN INFO = -1 @@ -160,9 +161,17 @@ CALL SSYTRF( UPLO, N, A, LDA, IPIV, WORK, LWORK, INFO ) IF( INFO.EQ.0 ) THEN * +* Convert A +* + CALL SSYCONV( UPLO, 'C', N, A, LDA, IPIV, WORK, IINFO ) +* * Solve the system A*X = B, overwriting B with X. * - CALL SSYTRS( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, INFO ) + CALL SSYTRS2( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, WORK, INFO ) +* +* Revert A +* + CALL SSYCONV( UPLO, 'R', N, A, LDA, IPIV, WORK, IINFO ) * END IF * diff --git a/SRC/ssytrs2.f b/SRC/ssytrs2.f new file mode 100644 index 00000000..cb98bbd9 --- /dev/null +++ b/SRC/ssytrs2.f @@ -0,0 +1,272 @@ + SUBROUTINE SSYTRS2( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, + $ WORK, INFO ) +* +* -- LAPACK PROTOTYPE routine (version 3.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( * ) + REAL A( LDA, * ), B( LDB, * ), WORK( * ) +* .. +* +* Purpose +* ======= +* +* SSYTRS2 solves a system of linear equations A*X = B with a real +* symmetric matrix A using the factorization A = U*D*U**T or +* A = L*D*L**T computed by SSYTRF and converted by SSYCONV. +* +* 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**T; +* = 'L': Lower triangular, form is A = L*D*L**T. +* +* 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) REAL array, dimension (LDA,N) +* The block diagonal matrix D and the multipliers used to +* obtain the factor U or L as computed by SSYTRF. +* +* 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 SSYTRF. +* +* B (input/output) REAL 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 .. + REAL ONE + PARAMETER ( ONE = 1.0E+0 ) +* .. +* .. Local Scalars .. + LOGICAL UPPER + INTEGER I, J, K, KP + REAL AK, AKM1, AKM1K, BK, BKM1, DENOM +* .. +* .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME +* .. +* .. External Subroutines .. + EXTERNAL SSCAL, SSWAP, STRSM, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC 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( 'SSYTRS2', -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 SSWAP( 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 SSWAP( 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 STRSM('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 + CALL SSCAL( NRHS, ONE / A( I, I ), B( I, 1 ), N ) + 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 ) / AKM1K + DENOM = AKM1*AK - ONE + DO 15 J = 1, NRHS + BKM1 = B( I-1, J ) / AKM1K + BK = B( I, J ) / 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 STRSM('L','U','T','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 SSWAP( 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 SSWAP( 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 SSWAP( 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 SSWAP( 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 STRSM('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 + CALL SSCAL( NRHS, ONE / A( I, I ), B( I, 1 ), N ) + ELSE + AKM1K = WORK(I) + AKM1 = A( I, I ) / AKM1K + AK = A( I+1, I+1 ) / AKM1K + DENOM = AKM1*AK - ONE + DO 25 J = 1, NRHS + BKM1 = B( I, J ) / 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 STRSM('L','L','T','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 SSWAP( 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 SSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) + K=K-2 + ENDIF + END DO +* + END IF +* + RETURN +* +* End of SSYTRS2 +* + END diff --git a/SRC/zsyconv.f b/SRC/zsyconv.f new file mode 100644 index 00000000..205c650e --- /dev/null +++ b/SRC/zsyconv.f @@ -0,0 +1,302 @@ + SUBROUTINE ZSYCONV( UPLO, WAY, N, A, LDA, IPIV, WORK, INFO ) +* +* -- LAPACK PROTOTYPE routine (version 3.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, WAY + INTEGER INFO, LDA, N +* .. +* .. Array Arguments .. + INTEGER IPIV( * ) + DOUBLE COMPLEX A( LDA, * ), WORK( * ) +* .. +* +* Purpose +* ======= +* +* ZSYCONV convert A given by TRF into L and D and vice-versa. +* Get Non-diag elements of D (returned in workspace) and +* apply or reverse permutation done in TRF. +* +* 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**T; +* = 'L': Lower triangular, form is A = L*D*L**T. +* +* WAY (input) CHARACTER*1 +* = 'C': Convert +* = 'R': Revert +* +* N (input) INTEGER +* The order of the matrix A. N >= 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 ZSYTRF. +* +* 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 ZSYTRF. +* +* WORK (workspace) DOUBLE COMPLEX array, dimension (N) +* +* LWORK (input) INTEGER +* The length of WORK. LWORK >=1. +* LWORK = N +* +* If LWORK = -1, then a workspace query is assumed; the routine +* only calculates the optimal size of the WORK array, returns +* this value as the first entry of the WORK array, and no error +* message related to LWORK is issued by XERBLA. +* +* INFO (output) INTEGER +* = 0: successful exit +* < 0: if INFO = -i, the i-th argument had an illegal value +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE COMPLEX ZERO + PARAMETER ( ZERO = (0.0D+0,0.0D+0) ) +* .. +* .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME +* +* .. External Subroutines .. + EXTERNAL XERBLA +* .. Local Scalars .. + LOGICAL UPPER, CONVERT + INTEGER I, IP + DOUBLE COMPLEX TEMP +* .. +* .. Executable Statements .. +* + INFO = 0 + UPPER = LSAME( UPLO, 'U' ) + CONVERT = LSAME( WAY, 'C' ) + IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN + INFO = -1 + ELSE IF( .NOT.CONVERT .AND. .NOT.LSAME( WAY, 'R' ) ) THEN + INFO = -2 + ELSE IF( N.LT.0 ) THEN + INFO = -3 + ELSE IF( LDA.LT.MAX( 1, N ) ) THEN + INFO = -5 + + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'ZSYCONV', -INFO ) + RETURN + END IF +* +* Quick return if possible +* + IF( N.EQ.0 ) + $ RETURN +* + IF( UPPER ) THEN +* +* A is UPPER +* +* Convert A (A is upper) +* +* Convert VALUE +* + IF ( CONVERT ) THEN + I=N + WORK(1)=ZERO + DO WHILE ( I .GT. 1 ) + IF( IPIV(I) .LT. 0 ) THEN + WORK(I)=A(I-1,I) + A(I-1,I)=ZERO + I=I-1 + ELSE + WORK(I)=ZERO + ENDIF + I=I-1 + END DO +* +* Convert PERMUTATIONS +* + I=N + DO WHILE ( I .GE. 1 ) + IF( IPIV(I) .GT. 0) THEN + IP=IPIV(I) + IF( I .LT. N) THEN + DO 12 J= I+1,N + TEMP=A(IP,J) + A(IP,J)=A(I,J) + A(I,J)=TEMP + 12 CONTINUE + ENDIF + ELSE + IP=-IPIV(I) + IF( I .LT. N) THEN + DO 13 J= I+1,N + TEMP=A(IP,J) + A(IP,J)=A(I-1,J) + A(I-1,J)=TEMP + 13 CONTINUE + ENDIF + I=I-1 + ENDIF + I=I-1 + END DO + + ELSE +* +* Revert A (A is upper) +* +* +* Revert PERMUTATIONS +* + I=1 + DO WHILE ( I .LE. N ) + IF( IPIV(I) .GT. 0 ) THEN + IP=IPIV(I) + IF( I .LT. N) THEN + DO J= I+1,N + TEMP=A(IP,J) + A(IP,J)=A(I,J) + A(I,J)=TEMP + END DO + ENDIF + ELSE + IP=-IPIV(I) + I=I+1 + IF( I .LT. N) THEN + DO J= I+1,N + TEMP=A(IP,J) + A(IP,J)=A(I-1,J) + A(I-1,J)=TEMP + END DO + ENDIF + ENDIF + I=I+1 + END DO +* +* Revert VALUE +* + I=N + DO WHILE ( I .GT. 1 ) + IF( IPIV(I) .LT. 0 ) THEN + A(I-1,I)=WORK(I) + I=I-1 + ENDIF + I=I-1 + END DO + END IF + ELSE +* +* A is LOWER +* + IF ( CONVERT ) THEN +* +* Convert A (A is lower) +* +* +* Convert VALUE +* + I=1 + WORK(N)=ZERO + DO WHILE ( I .LE. N ) + IF( I.LT.N .AND. IPIV(I) .LT. 0 ) THEN + WORK(I)=A(I+1,I) + A(I+1,I)=ZERO + I=I+1 + ELSE + WORK(I)=ZERO + ENDIF + I=I+1 + END DO +* +* Convert PERMUTATIONS +* + I=1 + DO WHILE ( I .LE. N ) + IF( IPIV(I) .GT. 0 ) THEN + IP=IPIV(I) + IF (I .GT. 1) THEN + DO 22 J= 1,I-1 + TEMP=A(IP,J) + A(IP,J)=A(I,J) + A(I,J)=TEMP + 22 CONTINUE + ENDIF + ELSE + IP=-IPIV(I) + IF (I .GT. 1) THEN + DO 23 J= 1,I-1 + TEMP=A(IP,J) + A(IP,J)=A(I+1,J) + A(I+1,J)=TEMP + 23 CONTINUE + ENDIF + I=I+1 + ENDIF + I=I+1 + END DO + ELSE +* +* Revert A (A is lower) +* +* +* Revert PERMUTATIONS +* + I=N + DO WHILE ( I .GE. 1 ) + IF( IPIV(I) .GT. 0 ) THEN + IP=IPIV(I) + IF (I .GT. 1) THEN + DO J= 1,I-1 + TEMP=A(I,J) + A(I,J)=A(IP,J) + A(IP,J)=TEMP + END DO + ENDIF + ELSE + IP=-IPIV(I) + I=I-1 + IF (I .GT. 1) THEN + DO J= 1,I-1 + TEMP=A(I+1,J) + A(I+1,J)=A(IP,J) + A(IP,J)=TEMP + END DO + ENDIF + ENDIF + I=I-1 + END DO +* +* Revert VALUE +* + I=1 + DO WHILE ( I .LE. N-1 ) + IF( IPIV(I) .LT. 0 ) THEN + A(I+1,I)=WORK(I) + I=I+1 + ENDIF + I=I+1 + END DO + END IF + END IF + + RETURN +* +* End of ZSYCONV +* + END diff --git a/SRC/zsysv.f b/SRC/zsysv.f index 2520e165..bf7baab7 100644 --- a/SRC/zsysv.f +++ b/SRC/zsysv.f @@ -105,7 +105,7 @@ * * .. Local Scalars .. LOGICAL LQUERY - INTEGER LWKOPT, NB + INTEGER LWKOPT, NB, IINFO * .. * .. External Functions .. LOGICAL LSAME @@ -113,7 +113,7 @@ EXTERNAL LSAME, ILAENV * .. * .. External Subroutines .. - EXTERNAL XERBLA, ZSYTRF, ZSYTRS + EXTERNAL ZSYCONV, XERBLA, ZSYTRF, ZSYTRS2 * .. * .. Intrinsic Functions .. INTRINSIC MAX @@ -123,6 +123,7 @@ * Test the input parameters. * INFO = 0 + IINFO = 0 LQUERY = ( LWORK.EQ.-1 ) IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN INFO = -1 @@ -160,9 +161,17 @@ CALL ZSYTRF( 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 ZSYTRS( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, INFO ) + CALL ZSYTRS2( 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/zsytrs2.f b/SRC/zsytrs2.f new file mode 100644 index 00000000..88d54408 --- /dev/null +++ b/SRC/zsytrs2.f @@ -0,0 +1,272 @@ + SUBROUTINE ZSYTRS2( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, + $ WORK, INFO ) +* +* -- LAPACK PROTOTYPE routine (version 3.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 +* ======= +* +* ZSYTRS2 solves a system of linear equations A*X = B with a real +* symmetric 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**T; +* = 'L': Lower triangular, form is A = L*D*L**T. +* +* 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 ZSYTRF. +* +* 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 ZSYTRF. +* +* 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 COMPLEX AK, AKM1, AKM1K, BK, BKM1, DENOM +* .. +* .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME +* .. +* .. External Subroutines .. + EXTERNAL ZSCAL, ZSWAP, ZTRSM, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC 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( 'ZSYTRS2', -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 + CALL ZSCAL( NRHS, ONE / A( I, I ), B( I, 1 ), N ) + 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 ) / AKM1K + DENOM = AKM1*AK - ONE + DO 15 J = 1, NRHS + BKM1 = B( I-1, J ) / AKM1K + BK = B( I, J ) / 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','T','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 + CALL ZSCAL( NRHS, ONE / A( I, I ), B( I, 1 ), N ) + ELSE + AKM1K = WORK(I) + AKM1 = A( I, I ) / AKM1K + AK = A( I+1, I+1 ) / AKM1K + DENOM = AKM1*AK - ONE + DO 25 J = 1, NRHS + BKM1 = B( I, J ) / 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','T','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 ZSYTRS2 +* + END diff --git a/TESTING/LIN/alahd.f b/TESTING/LIN/alahd.f index db2ad217..d147bc20 100644 --- a/TESTING/LIN/alahd.f +++ b/TESTING/LIN/alahd.f @@ -223,9 +223,35 @@ WRITE( IOUNIT, FMT = 9955 )7 WRITE( IOUNIT, FMT = '( '' Messages:'' )' ) * - ELSE IF( LSAMEN( 2, P2, 'SY' ) .OR. LSAMEN( 2, P2, 'SP' ) ) THEN + ELSE IF( LSAMEN( 2, P2, 'SY' ) ) THEN * * SY: Symmetric indefinite full +* + IF( LSAME( C3, 'Y' ) ) THEN + WRITE( IOUNIT, FMT = 9992 )PATH, 'Symmetric' + ELSE + WRITE( IOUNIT, FMT = 9991 )PATH, 'Symmetric' + 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, 'SP' ) ) THEN +* * SP: Symmetric indefinite packed * IF( LSAME( C3, 'Y' ) ) THEN diff --git a/TESTING/LIN/cchksy.f b/TESTING/LIN/cchksy.f index f0f90506..31a6ba21 100644 --- a/TESTING/LIN/cchksy.f +++ b/TESTING/LIN/cchksy.f @@ -22,7 +22,7 @@ * Purpose * ======= * -* CCHKSY tests CSYTRF, -TRI, -TRS, -RFS, and -CON. +* CCHKSY tests CSYTRF, -TRI, -TRS, -TRS2, -RFS, and -CON. * * Arguments * ========= @@ -94,7 +94,7 @@ INTEGER NTYPES PARAMETER ( NTYPES = 11 ) INTEGER NTESTS - PARAMETER ( NTESTS = 8 ) + PARAMETER ( NTESTS = 9 ) * .. * .. Local Scalars .. LOGICAL TRFCON, ZEROT @@ -406,12 +406,40 @@ $ 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 = 'CSYTRS2' + CALL CSYCONV( UPLO, 'C', N, AFAC, LDA, IWORK, WORK, + $ INFO ) + CALL CSYTRS2( UPLO, N, NRHS, AFAC, LDA, IWORK, X, + $ LDA, WORK, INFO ) + CALL CSYCONV( UPLO, 'R', N, AFAC, LDA, IWORK, WORK, + $ INFO ) +* +* Check error code from CSYTRS2. +* + IF( INFO.NE.0 ) + $ CALL ALAERH( PATH, 'CSYTRS2', INFO, 0, UPLO, N, + $ N, -1, -1, NRHS, IMAT, NFAIL, + $ NERRS, NOUT ) +* + CALL CLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA ) + CALL CSYT02( 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 = 'CSYRFS' @@ -428,15 +456,15 @@ $ 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 ) @@ -448,7 +476,7 @@ NRUN = NRUN + 5 130 CONTINUE * -*+ TEST 8 +*+ TEST 9 * Get an estimate of RCOND = 1/CNDNUM. * 140 CONTINUE @@ -463,16 +491,16 @@ $ CALL ALAERH( PATH, 'CSYCON', 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 ) + WRITE( NOUT, FMT = 9997 )UPLO, N, IMAT, 9, + $ RESULT( 9 ) NFAIL = NFAIL + 1 END IF NRUN = NRUN + 1 diff --git a/TESTING/LIN/dchksy.f b/TESTING/LIN/dchksy.f index 97e040e7..4bedee89 100644 --- a/TESTING/LIN/dchksy.f +++ b/TESTING/LIN/dchksy.f @@ -21,7 +21,7 @@ * Purpose * ======= * -* DCHKSY tests DSYTRF, -TRI, -TRS, -RFS, and -CON. +* DCHKSY tests DSYTRF, -TRI, -TRS, -TRS2, -RFS, and -CON. * * Arguments * ========= @@ -93,7 +93,7 @@ INTEGER NTYPES PARAMETER ( NTYPES = 10 ) INTEGER NTESTS - PARAMETER ( NTESTS = 8 ) + PARAMETER ( NTESTS = 9 ) * .. * .. Local Scalars .. LOGICAL TRFCON, ZEROT @@ -116,8 +116,8 @@ * .. External Subroutines .. EXTERNAL ALAERH, ALAHD, ALASUM, DERRSY, DGET04, DLACPY, $ DLARHS, DLATB4, DLATMS, DPOT02, DPOT03, DPOT05, - $ DSYCON, DSYRFS, DSYT01, DSYTRF, DSYTRI, DSYTRS, - $ XLAENV + $ DSYCON, DSYCONV, DSYRFS, DSYT01, DSYTRF, DSYTRI, + $ DSYTRS, DSYTRS2, XLAENV * .. * .. Intrinsic Functions .. INTRINSIC MAX, MIN @@ -371,7 +371,7 @@ DO 130 IRHS = 1, NNS NRHS = NSVAL( IRHS ) * -*+ TEST 3 +*+ TEST 3 ( Using TRS) * Solve and compute residual for A * X = B. * SRNAMT = 'DLARHS' @@ -395,13 +395,42 @@ CALL DPOT02( UPLO, N, NRHS, A, LDA, X, LDA, WORK, $ LDA, RWORK, RESULT( 3 ) ) * -*+ TEST 4 +*+ TEST 4 (Using TRS2) +* +* Solve and compute residual for A * X = B. +* + SRNAMT = 'DLARHS' + CALL DLARHS( PATH, XTYPE, UPLO, ' ', N, N, KL, KU, + $ NRHS, A, LDA, XACT, LDA, B, LDA, + $ ISEED, INFO ) + CALL DLACPY( 'Full', N, NRHS, B, LDA, X, LDA ) +* + SRNAMT = 'DSYTRS2' + CALL DSYCONV( UPLO, 'C', N, AFAC, LDA, IWORK, WORK, + $ INFO ) + CALL DSYTRS2( UPLO, N, NRHS, AFAC, LDA, IWORK, X, + $ LDA, WORK, INFO ) + CALL DSYCONV( UPLO, 'R', N, AFAC, LDA, IWORK, WORK, + $ INFO ) +* +* Check error code from DSYTRS2. +* + IF( INFO.NE.0 ) + $ CALL ALAERH( PATH, 'DSYTRS2', INFO, 0, UPLO, N, + $ N, -1, -1, NRHS, IMAT, NFAIL, + $ NERRS, NOUT ) +* + CALL DLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA ) + CALL DPOT02( UPLO, N, NRHS, A, LDA, X, LDA, WORK, + $ LDA, RWORK, RESULT( 4 ) ) +* +*+ TEST 5 * Check solution from generated exact solution. * CALL DGET04( 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 = 'DSYRFS' @@ -418,15 +447,15 @@ $ NERRS, NOUT ) * CALL DGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC, - $ RESULT( 5 ) ) + $ RESULT( 6 ) ) CALL DPOT05( 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 ) @@ -438,7 +467,7 @@ NRUN = NRUN + 5 130 CONTINUE * -*+ TEST 8 +*+ TEST 9 * Get an estimate of RCOND = 1/CNDNUM. * 140 CONTINUE @@ -453,16 +482,16 @@ $ CALL ALAERH( PATH, 'DSYCON', 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 diff --git a/TESTING/LIN/schksy.f b/TESTING/LIN/schksy.f index 2065292a..32739b68 100644 --- a/TESTING/LIN/schksy.f +++ b/TESTING/LIN/schksy.f @@ -21,7 +21,7 @@ * Purpose * ======= * -* SCHKSY tests SSYTRF, -TRI, -TRS, -RFS, and -CON. +* SCHKSY tests SSYTRF, -TRI, -TRS, -TRS2, -RFS, and -CON. * * Arguments * ========= @@ -93,7 +93,7 @@ INTEGER NTYPES PARAMETER ( NTYPES = 10 ) INTEGER NTESTS - PARAMETER ( NTESTS = 8 ) + PARAMETER ( NTESTS = 9 ) * .. * .. Local Scalars .. LOGICAL TRFCON, ZEROT @@ -116,8 +116,8 @@ * .. External Subroutines .. EXTERNAL ALAERH, ALAHD, ALASUM, SERRSY, SGET04, SLACPY, $ SLARHS, SLATB4, SLATMS, SPOT02, SPOT03, SPOT05, - $ SSYCON, SSYRFS, SSYT01, SSYTRF, SSYTRI, SSYTRS, - $ XLAENV + $ SSYCON, SSYCONV, SSYRFS, SSYT01, SSYTRF, SSYTRI, + $ SSYTRS, SSYTRS2, XLAENV * .. * .. Intrinsic Functions .. INTRINSIC MAX, MIN @@ -371,7 +371,7 @@ DO 130 IRHS = 1, NNS NRHS = NSVAL( IRHS ) * -*+ TEST 3 +*+ TEST 3 (Using DSYTRS) * Solve and compute residual for A * X = B. * SRNAMT = 'SLARHS' @@ -395,13 +395,41 @@ CALL SPOT02( UPLO, N, NRHS, A, LDA, X, LDA, WORK, $ LDA, RWORK, RESULT( 3 ) ) * -*+ TEST 4 +*+ TEST 4 (Using DSYTRS2) +* Solve and compute residual for A * X = B. +* + SRNAMT = 'SLARHS' + CALL SLARHS( PATH, XTYPE, UPLO, ' ', N, N, KL, KU, + $ NRHS, A, LDA, XACT, LDA, B, LDA, + $ ISEED, INFO ) + CALL SLACPY( 'Full', N, NRHS, B, LDA, X, LDA ) +* + SRNAMT = 'DSYTRS2' + CALL SSYCONV( UPLO, 'C', N, AFAC, LDA, IWORK, WORK, + $ INFO ) + CALL SSYTRS2( UPLO, N, NRHS, AFAC, LDA, IWORK, X, + $ LDA, WORK, INFO ) + CALL SSYCONV( UPLO, 'R', N, AFAC, LDA, IWORK, WORK, + $ INFO ) +* +* Check error code from SSYTRS2. +* + IF( INFO.NE.0 ) + $ CALL ALAERH( PATH, 'SSYTRS2', INFO, 0, UPLO, N, + $ N, -1, -1, NRHS, IMAT, NFAIL, + $ NERRS, NOUT ) +* + CALL SLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA ) + CALL SPOT02( UPLO, N, NRHS, A, LDA, X, LDA, WORK, + $ LDA, RWORK, RESULT( 4 ) ) +* +*+ TEST 5 * Check solution from generated exact solution. * CALL SGET04( 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 = 'SSYRFS' @@ -418,15 +446,15 @@ $ NERRS, NOUT ) * CALL SGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC, - $ RESULT( 5 ) ) + $ RESULT( 6 ) ) CALL SPOT05( 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 ) @@ -438,7 +466,7 @@ NRUN = NRUN + 5 130 CONTINUE * -*+ TEST 8 +*+ TEST 9 * Get an estimate of RCOND = 1/CNDNUM. * 140 CONTINUE @@ -453,16 +481,16 @@ $ CALL ALAERH( PATH, 'SSYCON', 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 ) + WRITE( NOUT, FMT = 9997 )UPLO, N, IMAT, 9, + $ RESULT( 9 ) NFAIL = NFAIL + 1 END IF NRUN = NRUN + 1 diff --git a/TESTING/LIN/zchksy.f b/TESTING/LIN/zchksy.f index fcc794ae..3995c08e 100644 --- a/TESTING/LIN/zchksy.f +++ b/TESTING/LIN/zchksy.f @@ -22,7 +22,7 @@ * Purpose * ======= * -* ZCHKSY tests ZSYTRF, -TRI, -TRS, -RFS, and -CON. +* ZCHKSY tests ZSYTRF, -TRI, -TRS, -TRS2, -RFS, and -CON. * * Arguments * ========= @@ -94,7 +94,7 @@ INTEGER NTYPES PARAMETER ( NTYPES = 11 ) INTEGER NTESTS - PARAMETER ( NTESTS = 8 ) + PARAMETER ( NTESTS = 9 ) * .. * .. Local Scalars .. LOGICAL TRFCON, ZEROT @@ -118,7 +118,7 @@ EXTERNAL ALAERH, ALAHD, ALASUM, XLAENV, ZERRSY, ZGET04, $ ZLACPY, ZLARHS, ZLATB4, ZLATMS, ZLATSY, ZPOT05, $ ZSYCON, ZSYRFS, ZSYT01, ZSYT02, ZSYT03, ZSYTRF, - $ ZSYTRI, ZSYTRS + $ ZSYTRI, ZSYTRS, ZSYTRS2 * .. * .. Intrinsic Functions .. INTRINSIC MAX, MIN @@ -381,7 +381,7 @@ DO 130 IRHS = 1, NNS NRHS = NSVAL( IRHS ) * -*+ TEST 3 +*+ TEST 3 (Using ZSYTRS) * Solve and compute residual for A * X = B. * SRNAMT = 'ZLARHS' @@ -405,13 +405,42 @@ CALL ZSYT02( UPLO, N, NRHS, A, LDA, X, LDA, WORK, $ LDA, RWORK, RESULT( 3 ) ) * -*+ TEST 4 +*+ TEST 4 (Using ZSYTRS2) +* 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 = 'ZSYTRS2' + CALL ZSYCONV( UPLO, 'C', N, AFAC, LDA, IWORK, WORK, + $ INFO ) + CALL ZSYTRS2( UPLO, N, NRHS, AFAC, LDA, IWORK, X, + $ LDA, WORK, INFO ) + CALL ZSYCONV( UPLO, 'R', N, AFAC, LDA, IWORK, WORK, + $ INFO ) +* +* Check error code from ZSYTRS. +* + IF( INFO.NE.0 ) + $ CALL ALAERH( PATH, 'ZSYTRS', INFO, 0, UPLO, N, + $ N, -1, -1, NRHS, IMAT, NFAIL, + $ NERRS, NOUT ) +* + CALL ZLACPY( 'Full', N, NRHS, B, LDA, WORK, LDA ) + CALL ZSYT02( 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 = 'ZSYRFS' @@ -428,15 +457,15 @@ $ 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 ) @@ -448,7 +477,7 @@ NRUN = NRUN + 5 130 CONTINUE * -*+ TEST 8 +*+ TEST 9 * Get an estimate of RCOND = 1/CNDNUM. * 140 CONTINUE @@ -463,16 +492,16 @@ $ CALL ALAERH( PATH, 'ZSYCON', 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 |