From 9205713fbc07fa5bcca7d17e74394430762d9aad Mon Sep 17 00:00:00 2001 From: julie Date: Tue, 2 Nov 2010 18:53:38 +0000 Subject: [PATCH] [xSYTRS/xSYSV] Hide the call to syconv inside trs2 to avoid changing the SYTRS interface. Update the testing accordingly [DSYTRI2] Comit dsytri2 to get some feedback Update the testing accordingly DSYTRI2 is the Level 3 blas Version of DSYTRI The actual routine that does the work is DSYTRI2X (name can be changed) DSYTRI2 is just a wrapper to allow to hide the 2D Workspace required by the routine. The interface had to be changed to integrate the possibility of doing a workspace query. DSYTRI2x implementation will be documented in a LAWN. This algorithm was inspired by the following paper: "Families of Algorithms Related to the Inversion of a Symmetric Positive Definite Matrix" PAOLO BIENTINESI Duke University and BRIAN GUNTER Delft University of Technology and ROBERT A. VAN DE GEIJN The University of Texas at Austin --- SRC/Makefile | 3 +- SRC/chesv.f | 12 +- SRC/chetrs2.f | 12 +- SRC/csysv.f | 12 +- SRC/csytrs2.f | 12 +- SRC/dsyconv.f | 9 +- SRC/dsysv.f | 16 +- SRC/dsyswapr.f | 125 +++++++++++++ SRC/dsytri2.f | 127 +++++++++++++ SRC/dsytri2x.f | 499 ++++++++++++++++++++++++++++++++++++++++++++++++++ SRC/dsytrs2.f | 21 ++- SRC/ssysv.f | 12 +- SRC/ssytrs2.f | 12 +- SRC/zhesv.f | 10 +- SRC/zhetrs2.f | 12 +- SRC/zsysv.f | 12 +- SRC/zsytrs2.f | 12 +- TESTING/LIN/cchkhe.f | 4 - TESTING/LIN/cchksy.f | 4 - TESTING/LIN/dchksy.f | 22 +-- TESTING/LIN/ddrvsy.f | 7 +- TESTING/LIN/ddrvsyx.f | 7 +- TESTING/LIN/derrsy.f | 15 +- TESTING/LIN/derrsyx.f | 15 +- TESTING/LIN/schksy.f | 4 - TESTING/LIN/zchkhe.f | 4 - TESTING/LIN/zchksy.f | 4 - 27 files changed, 880 insertions(+), 124 deletions(-) create mode 100644 SRC/dsyswapr.f create mode 100644 SRC/dsytri2.f create mode 100644 SRC/dsytri2x.f diff --git a/SRC/Makefile b/SRC/Makefile index 52db7fc..2ac3888 100644 --- a/SRC/Makefile +++ b/SRC/Makefile @@ -268,7 +268,8 @@ 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 dsytrs2.o dsyconv.o dtbcon.o \ + dsytd2.o dsytf2.o dsytrd.o dsytrf.o dsytri.o dsytri2.o dsytri2x.o \ + dsyswapr.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 \ diff --git a/SRC/chesv.f b/SRC/chesv.f index 47d5f68..c73e162 100644 --- a/SRC/chesv.f +++ b/SRC/chesv.f @@ -105,7 +105,7 @@ * * .. Local Scalars .. LOGICAL LQUERY - INTEGER IINFO, LWKOPT, NB + INTEGER LWKOPT, NB * .. * .. External Functions .. LOGICAL LSAME @@ -113,7 +113,7 @@ EXTERNAL ILAENV, LSAME * .. * .. External Subroutines .. - EXTERNAL CHETRF, CHETRS2, CSYCONV, XERBLA + EXTERNAL CHETRF, CHETRS2, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX @@ -160,18 +160,10 @@ 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 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 * WORK( 1 ) = LWKOPT diff --git a/SRC/chetrs2.f b/SRC/chetrs2.f index 1204a59..266cde6 100644 --- a/SRC/chetrs2.f +++ b/SRC/chetrs2.f @@ -73,7 +73,7 @@ * .. * .. Local Scalars .. LOGICAL UPPER - INTEGER I, J, K, KP + INTEGER I, IINFO, J, K, KP REAL S COMPLEX AK, AKM1, AKM1K, BK, BKM1, DENOM * .. @@ -82,7 +82,7 @@ EXTERNAL LSAME * .. * .. External Subroutines .. - EXTERNAL CSCAL, CSWAP, CTRSM, XERBLA + EXTERNAL CSCAL, CSYCONV, CSWAP, CTRSM, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC CONJG, MAX, REAL @@ -112,6 +112,10 @@ IF( N.EQ.0 .OR. NRHS.EQ.0 ) $ RETURN * +* Convert A +* + CALL CSYCONV( UPLO, 'C', N, A, LDA, IPIV, WORK, IINFO ) +* IF( UPPER ) THEN * * Solve A*X = B, where A = U*D*U'. @@ -268,6 +272,10 @@ * END IF * +* Revert A +* + CALL CSYCONV( UPLO, 'R', N, A, LDA, IPIV, WORK, IINFO ) +* RETURN * * End of CHETRS2 diff --git a/SRC/csysv.f b/SRC/csysv.f index 0654af7..1c97a79 100644 --- a/SRC/csysv.f +++ b/SRC/csysv.f @@ -105,7 +105,7 @@ * * .. Local Scalars .. LOGICAL LQUERY - INTEGER IINFO, LWKOPT, NB + INTEGER LWKOPT, NB * .. * .. External Functions .. LOGICAL LSAME @@ -113,7 +113,7 @@ EXTERNAL ILAENV, LSAME * .. * .. External Subroutines .. - EXTERNAL CSYCONV, CSYTRF, CSYTRS2, XERBLA + EXTERNAL CSYTRF, CSYTRS2, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX @@ -160,18 +160,10 @@ 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 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 * WORK( 1 ) = LWKOPT diff --git a/SRC/csytrs2.f b/SRC/csytrs2.f index bcdf52a..d3bb8ca 100644 --- a/SRC/csytrs2.f +++ b/SRC/csytrs2.f @@ -73,7 +73,7 @@ * .. * .. Local Scalars .. LOGICAL UPPER - INTEGER I, J, K, KP + INTEGER I, IINFO, J, K, KP COMPLEX AK, AKM1, AKM1K, BK, BKM1, DENOM * .. * .. External Functions .. @@ -81,7 +81,7 @@ EXTERNAL LSAME * .. * .. External Subroutines .. - EXTERNAL CSCAL, CSWAP, CTRSM, XERBLA + EXTERNAL CSCAL, CSYCONV, CSWAP, CTRSM, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX @@ -111,6 +111,10 @@ IF( N.EQ.0 .OR. NRHS.EQ.0 ) $ RETURN * +* Convert A +* + CALL CSYCONV( UPLO, 'C', N, A, LDA, IPIV, WORK, IINFO ) +* IF( UPPER ) THEN * * Solve A*X = B, where A = U*D*U'. @@ -265,6 +269,10 @@ * END IF * +* Revert A +* + CALL CSYCONV( UPLO, 'R', N, A, LDA, IPIV, WORK, IINFO ) +* RETURN * * End of CSYTRS2 diff --git a/SRC/dsyconv.f b/SRC/dsyconv.f index 15e9a73..3c74183 100644 --- a/SRC/dsyconv.f +++ b/SRC/dsyconv.f @@ -1,12 +1,11 @@ SUBROUTINE DSYCONV( UPLO, WAY, N, A, LDA, IPIV, WORK, INFO ) * -* -- LAPACK PROTOTYPE routine (version 3.2.2) -- -* -* -- Written by Julie Langou of the Univ. of TN -- -* May 2010 -* +* -- LAPACK PROTOTYPE routine (version 3.3.0) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2010 +* +* -- Written by Julie Langou of the Univ. of TN -- * * .. Scalar Arguments .. CHARACTER UPLO, WAY diff --git a/SRC/dsysv.f b/SRC/dsysv.f index 0501858..1bcd9aa 100644 --- a/SRC/dsysv.f +++ b/SRC/dsysv.f @@ -105,7 +105,7 @@ * * .. Local Scalars .. LOGICAL LQUERY - INTEGER IINFO, LWKOPT, NB + INTEGER LWKOPT, NB * .. * .. External Functions .. LOGICAL LSAME @@ -113,7 +113,7 @@ EXTERNAL LSAME, ILAENV * .. * .. External Subroutines .. - EXTERNAL DSYCONV, DSYTRF, DSYTRS2, XERBLA + EXTERNAL DSYTRF, DSYTRS2, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX @@ -161,17 +161,9 @@ CALL DSYTRF( UPLO, N, A, LDA, IPIV, WORK, LWORK, INFO ) IF( INFO.EQ.0 ) THEN * -* Convert A +* Solve the system A*X = B, overwriting B with X. * - CALL DSYCONV( UPLO, 'C', N, A, LDA, IPIV, WORK, IINFO ) -* -* Solve the system A*X = B, overwriting B with X. -* - CALL DSYTRS2( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, WORK, INFO ) -* -* Revert A -* - CALL DSYCONV( UPLO, 'R', N, A, LDA, IPIV, WORK, IINFO ) + CALL DSYTRS2( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, WORK, INFO ) * END IF * diff --git a/SRC/dsyswapr.f b/SRC/dsyswapr.f new file mode 100644 index 0000000..3ea542f --- /dev/null +++ b/SRC/dsyswapr.f @@ -0,0 +1,125 @@ + SUBROUTINE DSYSWAPR( UPLO, N, A, I1, I2) +* +* -- LAPACK routine (version 3.3.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2010 +* +* .. Scalar Arguments .. + CHARACTER UPLO + INTEGER I1, I2, N +* .. +* .. Array Arguments .. + DOUBLE PRECISION A(N,N) +* +* Purpose +* ======= +* +* DSYSWAPR swaps two rows of a lower or upper matrix +* +* 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. +* +* A (input/output) DOUBLE PRECISION array, dimension (LDA,N) +* On entry, the NB diagonal matrix D and the multipliers +* used to obtain the factor U or L as computed by DSYTRF. +* +* On exit, if INFO = 0, the (symmetric) inverse of the original +* matrix. If UPLO = 'U', the upper triangular part of the +* inverse is formed and the part of A below the diagonal is not +* referenced; if UPLO = 'L' the lower triangular part of the +* inverse is formed and the part of A above the diagonal is +* not referenced. +* +* I1 (input) INTEGER +* Index of the first row to swap +* +* I2 (input) INTEGER +* Index of the second row to swap +* +* ===================================================================== +* +* .. +* .. Local Scalars .. + LOGICAL UPPER + INTEGER I + DOUBLE PRECISION TMP +* +* .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME +* .. +* .. External Subroutines .. + EXTERNAL DSWAP +* .. +* .. Executable Statements .. +* + UPPER = LSAME( UPLO, 'U' ) + IF (UPPER) THEN +* +* UPPER +* first swap +* - swap column I1 and I2 from I1 to I1-1 + CALL DSWAP( I1-1, A(1,I1), 1, A(1,I2), 1 ) +* +* second swap : +* - swap A(I1,I1) and A(I2,I2) +* - swap row I1 from I1+1 to I2-1 with col I2 from I1+1 to I2-1 + TMP=A(I1,I1) + A(I1,I1)=A(I2,I2) + A(I2,I2)=TMP +* + DO I=1,I2-I1-1 + TMP=A(I1,I1+I) + A(I1,I1+I)=A(I1+I,I2) + A(I1+I,I2)=TMP + END DO +* +* third swap +* - swap row I1 and I2 from I2+1 to N + DO I=I2+1,N + TMP=A(I1,I) + A(I1,I)=A(I2,I) + A(I2,I)=TMP + END DO +* + ELSE +* +* LOWER +* first swap +* - swap row I1 and I2 from I1 to I1-1 + CALL DSWAP( I1-1, A(I1,1), N, A(I2,1), N ) +* +* second swap : +* - swap A(I1,I1) and A(I2,I2) +* - swap col I1 from I1+1 to I2-1 with row I2 from I1+1 to I2-1 + TMP=A(I1,I1) + A(I1,I1)=A(I2,I2) + A(I2,I2)=TMP +* + DO I=1,I2-I1-1 + TMP=A(I1+I,I1) + A(I1+I,I1)=A(I2,I1+I) + A(I2,I1+I)=TMP + END DO +* +* third swap +* - swap col I1 and I2 from I2+1 to N + DO I=I2+1,N + TMP=A(I,I1) + A(I,I1)=A(I,I2) + A(I,I2)=TMP + END DO +* + ENDIF + END SUBROUTINE DSYSWAPR + diff --git a/SRC/dsytri2.f b/SRC/dsytri2.f new file mode 100644 index 0000000..c04a0bd --- /dev/null +++ b/SRC/dsytri2.f @@ -0,0 +1,127 @@ + SUBROUTINE DSYTRI2( UPLO, N, A, LDA, IPIV, WORK, LWORK, INFO ) +* +* -- LAPACK routine (version 3.3.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2010 +* +* -- Written by Julie Langou of the Univ. of TN -- +* +* .. Scalar Arguments .. + CHARACTER UPLO + INTEGER INFO, LDA, LWORK, N +* .. +* .. Array Arguments .. + INTEGER IPIV( * ) + DOUBLE PRECISION A( LDA, * ), WORK( * ) +* .. +* +* Purpose +* ======= +* +* DSYTRI2 computes the inverse of a real symmetric indefinite matrix +* A using the factorization A = U*D*U**T or A = L*D*L**T computed by +* DSYTRF. DSYTRI2 set the LEADING DIMENSION of the workspace +* before calling DSYTRI2X that actually compute the inverse. +* +* 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. +* +* A (input/output) DOUBLE PRECISION array, dimension (LDA,N) +* On entry, the NB diagonal matrix D and the multipliers +* used to obtain the factor U or L as computed by DSYTRF. +* +* On exit, if INFO = 0, the (symmetric) inverse of the original +* matrix. If UPLO = 'U', the upper triangular part of the +* inverse is formed and the part of A below the diagonal is not +* referenced; if UPLO = 'L' the lower triangular part of the +* inverse is formed and the part of A above the diagonal is +* not referenced. +* +* 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 NB structure of D +* as determined by DSYTRF. +* +* WORK (workspace) DOUBLE PRECISION array, dimension (N+NB+1)*(NB+3) +* +* LWORK (input) INTEGER +* The dimension of the array WORK. +* WORK is size >= (N+NB+1)*(NB+3) +* If LDWORK = -1, then a workspace query is assumed; the routine +* 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 LDWORK is issued by XERBLA. +* +* INFO (output) INTEGER +* = 0: successful exit +* < 0: if INFO = -i, the i-th argument had an illegal value +* > 0: if INFO = i, D(i,i) = 0; the matrix is singular and its +* inverse could not be computed. +* +* ===================================================================== +* +* .. Local Scalars .. + LOGICAL UPPER, LQUERY + INTEGER MINSIZE, NBMAX +* .. +* .. External Functions .. + LOGICAL LSAME + INTEGER ILAENV + EXTERNAL LSAME, ILAENV +* .. +* .. External Subroutines .. + EXTERNAL DSYTRI2X +* .. +* .. Executable Statements .. +* +* Test the input parameters. +* + INFO = 0 + UPPER = LSAME( UPLO, 'U' ) + LQUERY = ( LWORK.EQ.-1 ) +* Get blocksize + NBMAX = ILAENV( 1, 'DSYTRF', UPLO, N, -1, -1, -1 ) + MINSIZE = (N+NBMAX+1)*(NBMAX+3) +* + IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN + INFO = -1 + ELSE IF( N.LT.0 ) THEN + INFO = -2 + ELSE IF( LDA.LT.MAX( 1, N ) ) THEN + INFO = -4 + ELSE IF (LWORK .LT. MINSIZE .AND. .NOT.LQUERY ) THEN + INFO = -7 + END IF +* +* Quick return if possible +* +* + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'DSYTRI2', -INFO ) + RETURN + ELSE IF( LQUERY ) THEN + WORK(1)=(N+NBMAX+1)*(NBMAX+3) + RETURN + END IF + IF( N.EQ.0 ) + $ RETURN + + CALL DSYTRI2X( UPLO, N, A, LDA, IPIV, WORK, NBMAX, INFO ) + RETURN +* +* End of DSYTRI2 +* + END diff --git a/SRC/dsytri2x.f b/SRC/dsytri2x.f new file mode 100644 index 0000000..d8ca317 --- /dev/null +++ b/SRC/dsytri2x.f @@ -0,0 +1,499 @@ + SUBROUTINE DSYTRI2X( UPLO, N, A, LDA, IPIV, WORK, NB, INFO ) +* +* -- LAPACK routine (version 3.3.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2010 +* +* -- Written by Julie Langou of the Univ. of TN -- +* +* .. Scalar Arguments .. + CHARACTER UPLO + INTEGER INFO, LDA, N +* .. +* .. Array Arguments .. + INTEGER IPIV( * ) + DOUBLE PRECISION A( LDA, * ), WORK( N+NB+1,* ) +* .. +* +* Purpose +* ======= +* +* DSYTRI2X computes the inverse of a real symmetric indefinite matrix +* A using the factorization A = U*D*U**T or A = L*D*L**T computed by +* DSYTRF. +* +* 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. +* +* A (input/output) DOUBLE PRECISION array, dimension (LDA,N) +* On entry, the NNB diagonal matrix D and the multipliers +* used to obtain the factor U or L as computed by DSYTRF. +* +* On exit, if INFO = 0, the (symmetric) inverse of the original +* matrix. If UPLO = 'U', the upper triangular part of the +* inverse is formed and the part of A below the diagonal is not +* referenced; if UPLO = 'L' the lower triangular part of the +* inverse is formed and the part of A above the diagonal is +* not referenced. +* +* 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 NNB structure of D +* as determined by DSYTRF. +* +* WORK (workspace) DOUBLE PRECISION array, dimension (N+NNB+1,NNB+3) +* +* NB (input) INTEGER +* Block size +* +* INFO (output) INTEGER +* = 0: successful exit +* < 0: if INFO = -i, the i-th argument had an illegal value +* > 0: if INFO = i, D(i,i) = 0; the matrix is singular and its +* inverse could not be computed. +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ONE, ZERO + PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) +* .. +* .. Local Scalars .. + LOGICAL UPPER + INTEGER I, IINFO, IP, K, CUT, NNB + INTEGER COUNT + INTEGER J, U11, INVD + + DOUBLE PRECISION AK, AKKP1, AKP1, D, T + DOUBLE PRECISION U01_I_J, U01_IP1_J + DOUBLE PRECISION U11_I_J, U11_IP1_J +* .. +* .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME +* .. +* .. External Subroutines .. + EXTERNAL DSCAL, DSYSCONV, DSWAP, XERBLA, DTRTRI + EXTERNAL DGEMM, DTRMM, DSYSWAPR +* .. +* .. Intrinsic Functions .. + INTRINSIC ABS, MAX +* .. +* .. Executable Statements .. +* +* Test the input parameters. +* + 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( LDA.LT.MAX( 1, N ) ) THEN + INFO = -4 + END IF +* +* Quick return if possible +* +* + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'DSYTRI2X', -INFO ) + RETURN + END IF + IF( N.EQ.0 ) + $ RETURN +* +* Convert A +* Workspace got Non-diag elements of D +* + CALL DSYCONV( UPLO, 'C', N, A, LDA, IPIV, WORK, IINFO ) +* +* Check that the diagonal matrix D is nonsingular. +* + IF( UPPER ) THEN +* +* Upper triangular storage: examine D from bottom to top +* + DO INFO = N, 1, -1 + IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO ) + $ RETURN + END DO + ELSE +* +* Lower triangular storage: examine D from top to bottom. +* + DO INFO = 1, N + IF( IPIV( INFO ).GT.0 .AND. A( INFO, INFO ).EQ.ZERO ) + $ RETURN + END DO + END IF + INFO = 0 +* +* Splitting Workspace +* U01 is a block (N,NB+1) +* The first element of U01 is in WORK(1,1) +* U11 is a block (NB+1,NB+1) +* The first element of U11 is in WORK(N+1,1) + U11 = N +* INVD is a block (N,2) +* The first element of INVD is in WORK(1,INVD) + INVD = NB+2 + COORD_J=NB+3 + + IF( UPPER ) THEN +* +* invA = P * inv(U')*inv(D)*inv(U)*P'. +* + CALL DTRTRI( UPLO, 'U', N, A, LDA, INFO ) +* +* inv(D) and inv(D)*inv(U) +* + K=1 + DO WHILE ( K .LE. N ) + IF( IPIV( K ).GT.0 ) THEN +* 1 x 1 diagonal NNB + WORK(K,INVD) = 1/ A( K, K ) + WORK(K,INVD+1) = 0 + K=K+1 + ELSE +* 2 x 2 diagonal NNB + T = ABS( WORK(K+1,1) ) + AK = A( K, K ) / T + AKP1 = A( K+1, K+1 ) / T + AKKP1 = WORK(K+1,1) / T + D = T*( AK*AKP1-ONE ) + WORK(K,INVD) = AKP1 / D + WORK(K+1,INVD+1) = AK / D + WORK(K,INVD+1) = -AKKP1 / D + WORK(K+1,INVD) = -AKKP1 / D + K=K+2 + END IF + END DO +* +* inv(U') = (inv(U))' +* +* inv(U')*inv(D)*inv(U) +* + CUT=N + DO WHILE (CUT .GT. 0) + NNB=NB + IF (CUT .LE. NNB) THEN + NNB=CUT + ELSE + COUNT = 0 +* count negative elements, + DO I=CUT+1-NNB,CUT + IF (IPIV(I) .LT. 0) COUNT=COUNT+1 + END DO +* need a even number for a clear cut + IF (MOD(COUNT,2) .EQ. 1) NNB=NNB+1 + END IF + + CUT=CUT-NNB +* +* U01 Block +* + DO I=1,CUT + DO J=1,NNB + WORK(I,J)=A(I,CUT+J) + END DO + END DO +* +* U11 Block +* + DO I=1,NNB + WORK(U11+I,I)=ONE + DO J=1,I-1 + WORK(U11+I,J)=ZERO + END DO + DO J=I+1,NNB + WORK(U11+I,J)=A(CUT+I,CUT+J) + END DO + END DO +* +* invD*U01 +* + I=1 + DO WHILE (I .LE. CUT) + IF (IPIV(I) > 0) THEN + DO J=1,NNB + WORK(I,J)=WORK(I,INVD)*WORK(I,J) + END DO + I=I+1 + ELSE + DO J=1,NNB + U01_I_J = WORK(I,J) + U01_IP1_J = WORK(I+1,J) + WORK(I,J)=WORK(I,INVD)*U01_I_J+ + $ WORK(I,INVD+1)*U01_IP1_J + WORK(I+1,J)=WORK(I+1,INVD)*U01_I_J+ + $ WORK(I+1,INVD+1)*U01_IP1_J + END DO + I=I+2 + END IF + END DO +* +* invD1*U11 +* + I=1 + DO WHILE (I .LE. NNB) + IF (IPIV(CUT+I) > 0) THEN + DO J=I,NNB + WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J) + END DO + I=I+1 + ELSE + DO J=I,NNB + U11_I_J = WORK(U11+I,J) + U11_IP1_J = WORK(U11+I+1,J) + WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J) + + $ WORK(CUT+I,INVD+1)*WORK(U11+I+1,J) + WORK(U11+I+1,J)=WORK(CUT+I+1,INVD)*U11_I_J+ + $ WORK(CUT+I+1,INVD+1)*U11_IP1_J + END DO + I=I+2 + END IF + END DO +* +* U11T*invD1*U11->U11 +* +* SUBROUTINE DTRMM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB) + + CALL DTRMM('L','U','T','U',NNB, NNB, + $ ONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1) + +* +* U01'invD*U01->A(CUT+I,CUT+J) +* + CALL DGEMM('T','N',NNB,NNB,CUT,ONE,A(1,CUT+1),LDA, + $ WORK,N+NB+1, ZERO, A(CUT+1,CUT+1), LDA) +* +* U11 = U11T*invD1*U11 + U01'invD*U01 (Prem + Deus) +* + DO I=1,NNB + DO J=I,NNB + A(CUT+I,CUT+J)=A(CUT+I,CUT+J)+WORK(U11+I,J) + END DO + END DO +* +* U01 = U00T*invD0*U01 +* + CALL DTRMM('L',UPLO,'T','U',CUT, NNB, + $ ONE,A,LDA,WORK,N+NB+1) + +* +* Update U01 +* + DO I=1,CUT + DO J=1,NNB + A(I,CUT+J)=WORK(I,J) + END DO + END DO +* Next Block + END DO +* +* Apply PERMUTATIONS P and P': P * inv(U')*inv(D)*inv(U) *P' +* + I=1 + DO WHILE ( I .LE. N ) + IF( IPIV(I) .GT. 0 ) THEN + IP=IPIV(I) + IF (I .LT. IP) CALL DSYSWAPR( UPLO, N, A, I ,IP ) + IF (I .GT. IP) CALL DSYSWAPR( UPLO, N, A, IP ,I ) + ELSE + IP=-IPIV(I) + I=I+1 + IF ( (I-1) .LT. IP) + $ CALL DSYSWAPR( UPLO, N, A, I-1 ,IP ) + IF ( (I-1) .GT. IP) + $ CALL DSYSWAPR( UPLO, N, A, IP ,I-1 ) + ENDIF + I=I+1 + END DO + + DO I=1,N + DO J=I+1,N + A(J,I)=A(I,J) + END DO + END DO + ELSE +* +* LOWER... +* +* invA = P * inv(U')*inv(D)*inv(U)*P'. +* + CALL DTRTRI( UPLO, 'U', N, A, LDA, INFO ) +* +* inv(D) and inv(D)*inv(U) +* + K=N + DO WHILE ( K .GE. 1 ) + IF( IPIV( K ).GT.0 ) THEN +* 1 x 1 diagonal NNB + WORK(K,INVD) = 1/ A( K, K ) + WORK(K,INVD+1) = 0 + K=K-1 + ELSE +* 2 x 2 diagonal NNB + T = ABS(WORK(K-1,1)) + AK = A( K-1, K-1 ) / T + AKP1 = A( K, K ) / T + AKKP1 = WORK(K-1,1) / T + D = T*( AK*AKP1-ONE ) + WORK(K-1,INVD) = AKP1 / D + WORK(K,INVD) = AK / D + WORK(K,INVD+1) = -AKKP1 / D + WORK(K-1,INVD+1) = -AKKP1 / D + K=K-2 + END IF + END DO +* +* inv(U') = (inv(U))' +* +* inv(U')*inv(D)*inv(U) +* + CUT=0 + DO WHILE (CUT .LT. N) + NNB=NB + IF (CUT + NNB .GE. N) THEN + NNB=N-CUT + ELSE + COUNT = 0 +* count negative elements, + DO I=CUT+1,CUT+NNB + IF (IPIV(I) .LT. 0) COUNT=COUNT+1 + END DO +* need a even number for a clear cut + IF (MOD(COUNT,2) .EQ. 1) NNB=NNB+1 + END IF +* L21 Block + DO I=1,N-CUT-NNB + DO J=1,NNB + WORK(I,J)=A(CUT+NNB+I,CUT+J) + END DO + END DO +* L11 Block + DO I=1,NNB + WORK(U11+I,I)=ONE + DO J=I+1,NNB + WORK(U11+I,J)=ZERO + END DO + DO J=1,I-1 + WORK(U11+I,J)=A(CUT+I,CUT+J) + END DO + END DO +* +* invD*L21 +* + I=N-CUT-NNB + DO WHILE (I .GE. 1) + IF (IPIV(CUT+NNB+I) > 0) THEN + DO J=1,NNB + WORK(I,J)=WORK(CUT+NNB+I,INVD)*WORK(I,J) + END DO + I=I-1 + ELSE + DO J=1,NNB + U01_I_J = WORK(I,J) + U01_IP1_J = WORK(I-1,J) + WORK(I,J)=WORK(CUT+NNB+I,INVD)*U01_I_J+ + $ WORK(CUT+NNB+I,INVD+1)*U01_IP1_J + WORK(I-1,J)=WORK(CUT+NNB+I-1,INVD+1)*U01_I_J+ + $ WORK(CUT+NNB+I-1,INVD)*U01_IP1_J + END DO + I=I-2 + END IF + END DO +* +* invD1*L11 +* + I=NNB + DO WHILE (I .GE. 1) + IF (IPIV(CUT+I) > 0) THEN + DO J=1,NNB + WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J) + END DO + I=I-1 + ELSE + DO J=1,NNB + U11_I_J = WORK(U11+I,J) + U11_IP1_J = WORK(U11+I-1,J) + WORK(U11+I,J)=WORK(CUT+I,INVD)*WORK(U11+I,J) + + $ WORK(CUT+I,INVD+1)*U11_IP1_J + WORK(U11+I-1,J)=WORK(CUT+I-1,INVD+1)*U11_I_J+ + $ WORK(CUT+I-1,INVD)*U11_IP1_J + END DO + I=I-2 + END IF + END DO +* +* U11T*invD1*U11->U11 +* + CALL DTRMM('L',UPLO,'T','U',NNB, NNB, + $ ONE,A(CUT+1,CUT+1),LDA,WORK(U11+1,1),N+NB+1) +* +* L21T*invD2*L21->A(CUT+I,CUT+J) +* + CALL DGEMM('T','N',NNB,NNB,N-NNB-CUT,ONE,A(CUT+NNB+1,CUT+1) + $ ,LDA,WORK,N+NB+1, ZERO, A(CUT+1,CUT+1), LDA) + +* +* L11 = L11T*invD1*L11 + U01'invD*U01 (Prem + Deus) +* + DO I=1,NNB + DO J=1,I + A(CUT+I,CUT+J)=A(CUT+I,CUT+J)+WORK(U11+I,J) + END DO + END DO +* +* U01 = L22T*invD2*L21 +* + CALL DTRMM('L',UPLO,'T','U', N-NNB-CUT, NNB, + $ ONE,A(CUT+NNB+1,CUT+NNB+1),LDA,WORK,N+NB+1) + +* Update L21 + DO I=1,N-CUT-NNB + DO J=1,NNB + A(CUT+NNB+I,CUT+J)=WORK(I,J) + END DO + END DO +* Next Block + CUT=CUT+NNB + END DO +* +* Apply PERMUTATIONS P and P': P * inv(U')*inv(D)*inv(U) *P' +* + I=N + DO WHILE ( I .GE. 1 ) + IF( IPIV(I) .GT. 0 ) THEN + IP=IPIV(I) + IF (I .LT. IP) CALL DSYSWAPR( UPLO, N, A, I ,IP ) + IF (I .GT. IP) CALL DSYSWAPR( UPLO, N, A, IP ,I ) + ELSE + IP=-IPIV(I) + IF ( I .LT. IP) CALL DSYSWAPR( UPLO, N, A, I ,IP ) + IF ( I .GT. IP) CALL DSYSWAPR( UPLO, N, A, IP ,I ) + I=I-1 + ENDIF + I=I-1 + END DO + END IF +* + RETURN +* +* End of DSYTRI2X +* + END + diff --git a/SRC/dsytrs2.f b/SRC/dsytrs2.f index 120f2ff..d00e815 100644 --- a/SRC/dsytrs2.f +++ b/SRC/dsytrs2.f @@ -1,13 +1,12 @@ SUBROUTINE DSYTRS2( 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 PROTOTYPE routine (version 3.3.0) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2010 +* +* -- Written by Julie Langou of the Univ. of TN -- * * .. Scalar Arguments .. CHARACTER UPLO @@ -73,7 +72,7 @@ * .. * .. Local Scalars .. LOGICAL UPPER - INTEGER I, J, K, KP + INTEGER I, IINFO, J, K, KP DOUBLE PRECISION AK, AKM1, AKM1K, BK, BKM1, DENOM * .. * .. External Functions .. @@ -81,7 +80,7 @@ EXTERNAL LSAME * .. * .. External Subroutines .. - EXTERNAL DSCAL, DSWAP, DTRSM, XERBLA + EXTERNAL DSCAL, DSYCONV, DSWAP, DTRSM, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX @@ -111,6 +110,10 @@ IF( N.EQ.0 .OR. NRHS.EQ.0 ) $ RETURN * +* Convert A +* + CALL DSYCONV( UPLO, 'C', N, A, LDA, IPIV, WORK, IINFO ) +* IF( UPPER ) THEN * * Solve A*X = B, where A = U*D*U'. @@ -265,6 +268,10 @@ * END IF * +* Revert A +* + CALL DSYCONV( UPLO, 'R', N, A, LDA, IPIV, WORK, IINFO ) +* RETURN * * End of DSYTRS2 diff --git a/SRC/ssysv.f b/SRC/ssysv.f index ce2ddd7..f1c5cf2 100644 --- a/SRC/ssysv.f +++ b/SRC/ssysv.f @@ -105,7 +105,7 @@ * * .. Local Scalars .. LOGICAL LQUERY - INTEGER IINFO, LWKOPT, NB + INTEGER LWKOPT, NB * .. * .. External Functions .. LOGICAL LSAME @@ -113,7 +113,7 @@ EXTERNAL ILAENV, LSAME * .. * .. External Subroutines .. - EXTERNAL SSYCONV, SSYTRF, SSYTRS2, XERBLA + EXTERNAL SSYTRF, SSYTRS2, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX @@ -161,18 +161,10 @@ 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 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 * WORK( 1 ) = LWKOPT diff --git a/SRC/ssytrs2.f b/SRC/ssytrs2.f index 93c009c..885373d 100644 --- a/SRC/ssytrs2.f +++ b/SRC/ssytrs2.f @@ -73,7 +73,7 @@ * .. * .. Local Scalars .. LOGICAL UPPER - INTEGER I, J, K, KP + INTEGER I, IINFO, J, K, KP REAL AK, AKM1, AKM1K, BK, BKM1, DENOM * .. * .. External Functions .. @@ -81,7 +81,7 @@ EXTERNAL LSAME * .. * .. External Subroutines .. - EXTERNAL SSCAL, SSWAP, STRSM, XERBLA + EXTERNAL SSCAL, SSYCONV, SSWAP, STRSM, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX @@ -111,6 +111,10 @@ IF( N.EQ.0 .OR. NRHS.EQ.0 ) $ RETURN * +* Convert A +* + CALL SSYCONV( UPLO, 'C', N, A, LDA, IPIV, WORK, IINFO ) +* IF( UPPER ) THEN * * Solve A*X = B, where A = U*D*U'. @@ -265,6 +269,10 @@ * END IF * +* Revert A +* + CALL SSYCONV( UPLO, 'R', N, A, LDA, IPIV, WORK, IINFO ) +* RETURN * * End of SSYTRS2 diff --git a/SRC/zhesv.f b/SRC/zhesv.f index f539a27..cedbfa0 100644 --- a/SRC/zhesv.f +++ b/SRC/zhesv.f @@ -113,7 +113,7 @@ EXTERNAL LSAME, ILAENV * .. * .. External Subroutines .. - EXTERNAL XERBLA, ZHETRF, ZHETRS2, ZSYCONV + EXTERNAL XERBLA, ZHETRF, ZHETRS2 * .. * .. Intrinsic Functions .. INTRINSIC MAX @@ -160,18 +160,10 @@ 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 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 * WORK( 1 ) = LWKOPT diff --git a/SRC/zhetrs2.f b/SRC/zhetrs2.f index 91a53cd..8ce2cc8 100644 --- a/SRC/zhetrs2.f +++ b/SRC/zhetrs2.f @@ -73,7 +73,7 @@ * .. * .. Local Scalars .. LOGICAL UPPER - INTEGER I, J, K, KP + INTEGER I, IINFO, J, K, KP DOUBLE PRECISION S DOUBLE COMPLEX AK, AKM1, AKM1K, BK, BKM1, DENOM * .. @@ -82,7 +82,7 @@ EXTERNAL LSAME * .. * .. External Subroutines .. - EXTERNAL ZLACGV, ZSCAL, ZSWAP, ZTRSM, XERBLA + EXTERNAL ZLACGV, ZSCAL, ZSYCONV, ZSWAP, ZTRSM, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC DBLE, DCONJG, MAX @@ -112,6 +112,10 @@ IF( N.EQ.0 .OR. NRHS.EQ.0 ) $ RETURN * +* Convert A +* + CALL ZSYCONV( UPLO, 'C', N, A, LDA, IPIV, WORK, IINFO ) +* IF( UPPER ) THEN * * Solve A*X = B, where A = U*D*U'. @@ -268,6 +272,10 @@ * END IF * +* Revert A +* + CALL ZSYCONV( UPLO, 'R', N, A, LDA, IPIV, WORK, IINFO ) +* RETURN * * End of ZHETRS2 diff --git a/SRC/zsysv.f b/SRC/zsysv.f index 19ef672..8332abe 100644 --- a/SRC/zsysv.f +++ b/SRC/zsysv.f @@ -105,7 +105,7 @@ * * .. Local Scalars .. LOGICAL LQUERY - INTEGER IINFO, LWKOPT, NB + INTEGER LWKOPT, NB * .. * .. External Functions .. LOGICAL LSAME @@ -113,7 +113,7 @@ EXTERNAL LSAME, ILAENV * .. * .. External Subroutines .. - EXTERNAL ZSYCONV, XERBLA, ZSYTRF, ZSYTRS2 + EXTERNAL XERBLA, ZSYTRF, ZSYTRS2 * .. * .. Intrinsic Functions .. INTRINSIC MAX @@ -161,18 +161,10 @@ 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 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 * WORK( 1 ) = LWKOPT diff --git a/SRC/zsytrs2.f b/SRC/zsytrs2.f index 4656923..2c53092 100644 --- a/SRC/zsytrs2.f +++ b/SRC/zsytrs2.f @@ -73,7 +73,7 @@ * .. * .. Local Scalars .. LOGICAL UPPER - INTEGER I, J, K, KP + INTEGER I, IINFO, J, K, KP DOUBLE COMPLEX AK, AKM1, AKM1K, BK, BKM1, DENOM * .. * .. External Functions .. @@ -81,7 +81,7 @@ EXTERNAL LSAME * .. * .. External Subroutines .. - EXTERNAL ZSCAL, ZSWAP, ZTRSM, XERBLA + EXTERNAL ZSCAL, ZSYCONV, ZSWAP, ZTRSM, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX @@ -111,6 +111,10 @@ IF( N.EQ.0 .OR. NRHS.EQ.0 ) $ RETURN * +* Convert A +* + CALL ZSYCONV( UPLO, 'C', N, A, LDA, IPIV, WORK, IINFO ) +* IF( UPPER ) THEN * * Solve A*X = B, where A = U*D*U'. @@ -265,6 +269,10 @@ * END IF * +* Revert A +* + CALL ZSYCONV( UPLO, 'R', N, A, LDA, IPIV, WORK, IINFO ) +* RETURN * * End of ZSYTRS2 diff --git a/TESTING/LIN/cchkhe.f b/TESTING/LIN/cchkhe.f index 77b7e8f..534d70e 100644 --- a/TESTING/LIN/cchkhe.f +++ b/TESTING/LIN/cchkhe.f @@ -409,12 +409,8 @@ 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. * diff --git a/TESTING/LIN/cchksy.f b/TESTING/LIN/cchksy.f index bec003d..927d1c0 100644 --- a/TESTING/LIN/cchksy.f +++ b/TESTING/LIN/cchksy.f @@ -415,12 +415,8 @@ 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. * diff --git a/TESTING/LIN/dchksy.f b/TESTING/LIN/dchksy.f index 25ee05e..6f49079 100644 --- a/TESTING/LIN/dchksy.f +++ b/TESTING/LIN/dchksy.f @@ -21,7 +21,7 @@ * Purpose * ======= * -* DCHKSY tests DSYTRF, -TRI, -TRS, -TRS2, -RFS, and -CON. +* DCHKSY tests DSYTRF, -TRI2, -TRS, -TRS2, -RFS, and -CON. * * Arguments * ========= @@ -108,6 +108,7 @@ CHARACTER UPLOS( 2 ) INTEGER ISEED( 4 ), ISEEDY( 4 ) DOUBLE PRECISION RESULT( NTESTS ) + DOUBLE PRECISION MYWORK( NTESTS ) * .. * .. External Functions .. DOUBLE PRECISION DGET06, DLANSY @@ -116,8 +117,8 @@ * .. External Subroutines .. EXTERNAL ALAERH, ALAHD, ALASUM, DERRSY, DGET04, DLACPY, $ DLARHS, DLATB4, DLATMS, DPOT02, DPOT03, DPOT05, - $ DSYCON, DSYCONV, DSYRFS, DSYT01, DSYTRF, - $ DSYTRI, DSYTRS, DSYTRS2, XLAENV + $ DSYCON, DSYRFS, DSYT01, DSYTRF, + $ DSYTRI2, DSYTRS, DSYTRS2, XLAENV * .. * .. Intrinsic Functions .. INTRINSIC MAX, MIN @@ -325,14 +326,15 @@ * IF( INB.EQ.1 .AND. .NOT.TRFCON ) THEN CALL DLACPY( UPLO, N, N, AFAC, LDA, AINV, LDA ) - SRNAMT = 'DSYTRI' - CALL DSYTRI( UPLO, N, AINV, LDA, IWORK, WORK, - $ INFO ) + SRNAMT = 'DSYTRI2' + LWORK = (N+NB+1)*(NB+3) + CALL DSYTRI2( UPLO, N, AINV, LDA, IWORK, WORK, + $ LWORK, INFO ) * -* Check error code from DSYTRI. +* Check error code from DSYTRI2. * IF( INFO.NE.0 ) - $ CALL ALAERH( PATH, 'DSYTRI', INFO, -1, UPLO, N, + $ CALL ALAERH( PATH, 'DSYTRI2', INFO, -1, UPLO, N, $ N, -1, -1, -1, IMAT, NFAIL, NERRS, $ NOUT ) * @@ -406,12 +408,8 @@ 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. * diff --git a/TESTING/LIN/ddrvsy.f b/TESTING/LIN/ddrvsy.f index 33bf412..40ffad4 100644 --- a/TESTING/LIN/ddrvsy.f +++ b/TESTING/LIN/ddrvsy.f @@ -106,7 +106,7 @@ * .. External Subroutines .. EXTERNAL ALADHD, ALAERH, ALASVM, DERRVX, DGET04, DLACPY, $ DLARHS, DLASET, DLATB4, DLATMS, DPOT02, DPOT05, - $ DSYSV, DSYSVX, DSYT01, DSYTRF, DSYTRI, XLAENV + $ DSYSV, DSYSVX, DSYT01, DSYTRF, DSYTRI2, XLAENV * .. * .. Scalars in Common .. LOGICAL LERR, OK @@ -294,8 +294,9 @@ * Compute inv(A) and take its norm. * CALL DLACPY( UPLO, N, N, AFAC, LDA, AINV, LDA ) - CALL DSYTRI( UPLO, N, AINV, LDA, IWORK, WORK, - $ INFO ) + LWORK = (N+NB+1)*(NB+3) + CALL DSYTRI2( UPLO, N, AINV, LDA, IWORK, WORK, + $ LWORK, INFO ) AINVNM = DLANSY( '1', UPLO, N, AINV, LDA, RWORK ) * * Compute the 1-norm condition number of A. diff --git a/TESTING/LIN/ddrvsyx.f b/TESTING/LIN/ddrvsyx.f index 11b1116..fca0f37 100644 --- a/TESTING/LIN/ddrvsyx.f +++ b/TESTING/LIN/ddrvsyx.f @@ -112,7 +112,7 @@ * .. External Subroutines .. EXTERNAL ALADHD, ALAERH, ALASVM, DERRVX, DGET04, DLACPY, $ DLARHS, DLASET, DLATB4, DLATMS, DPOT02, DPOT05, - $ DSYSV, DSYSVX, DSYT01, DSYTRF, DSYTRI, XLAENV, + $ DSYSV, DSYSVX, DSYT01, DSYTRF, DSYTRI2, XLAENV, $ DSYSVXX * .. * .. Scalars in Common .. @@ -301,8 +301,9 @@ * Compute inv(A) and take its norm. * CALL DLACPY( UPLO, N, N, AFAC, LDA, AINV, LDA ) - CALL DSYTRI( UPLO, N, AINV, LDA, IWORK, WORK, - $ INFO ) + LWORK = (N+NB+1)*(NB+3) + CALL DSYTRI2( UPLO, N, AINV, LDA, IWORK, WORK, + $ LWORK, INFO ) AINVNM = DLANSY( '1', UPLO, N, AINV, LDA, RWORK ) * * Compute the 1-norm condition number of A. diff --git a/TESTING/LIN/derrsy.f b/TESTING/LIN/derrsy.f index 27098be..3618681 100644 --- a/TESTING/LIN/derrsy.f +++ b/TESTING/LIN/derrsy.f @@ -47,7 +47,7 @@ * .. External Subroutines .. EXTERNAL ALAESM, CHKXER, DSPCON, DSPRFS, DSPTRF, DSPTRI, $ DSPTRS, DSYCON, DSYRFS, DSYTF2, DSYTRF, DSYTRI, - $ DSYTRS + $ DSYTRI2, DSYTRS * .. * .. Scalars in Common .. LOGICAL LERR, OK @@ -130,6 +130,19 @@ CALL DSYTRI( 'U', 2, A, 1, IP, W, INFO ) CALL CHKXER( 'DSYTRI', INFOT, NOUT, LERR, OK ) * +* DSYTRI2 +* + SRNAMT = 'DSYTRI2' + INFOT = 1 + CALL DSYTRI2( '/', 0, A, 1, IP, W, IW, INFO ) + CALL CHKXER( 'DSYTRI2', INFOT, NOUT, LERR, OK ) + INFOT = 2 + CALL DSYTRI2( 'U', -1, A, 1, IP, W, IW, INFO ) + CALL CHKXER( 'DSYTRI2', INFOT, NOUT, LERR, OK ) + INFOT = 4 + CALL DSYTRI2( 'U', 2, A, 1, IP, W, IW, INFO ) + CALL CHKXER( 'DSYTRI2', INFOT, NOUT, LERR, OK ) +* * DSYTRS * SRNAMT = 'DSYTRS' diff --git a/TESTING/LIN/derrsyx.f b/TESTING/LIN/derrsyx.f index 3c65ff4..086a29b 100644 --- a/TESTING/LIN/derrsyx.f +++ b/TESTING/LIN/derrsyx.f @@ -53,7 +53,7 @@ * .. External Subroutines .. EXTERNAL ALAESM, CHKXER, DSPCON, DSPRFS, DSPTRF, DSPTRI, $ DSPTRS, DSYCON, DSYRFS, DSYTF2, DSYTRF, DSYTRI, - $ DSYTRS, DSYRFSX + $ DSYTRI2, DSYTRS, DSYRFSX * .. * .. Scalars in Common .. LOGICAL LERR, OK @@ -137,6 +137,19 @@ CALL DSYTRI( 'U', 2, A, 1, IP, W, INFO ) CALL CHKXER( 'DSYTRI', INFOT, NOUT, LERR, OK ) * +* DSYTRI2 +* + SRNAMT = 'DSYTRI2' + INFOT = 1 + CALL DSYTRI2( '/', 0, A, 1, IP, W, IW, INFO ) + CALL CHKXER( 'DSYTRI2', INFOT, NOUT, LERR, OK ) + INFOT = 2 + CALL DSYTRI2( 'U', -1, A, 1, IP, W, IW, INFO ) + CALL CHKXER( 'DSYTRI2', INFOT, NOUT, LERR, OK ) + INFOT = 4 + CALL DSYTRI2( 'U', 2, A, 1, IP, W, IW, INFO ) + CALL CHKXER( 'DSYTRI2', INFOT, NOUT, LERR, OK ) +* * DSYTRS * SRNAMT = 'DSYTRS' diff --git a/TESTING/LIN/schksy.f b/TESTING/LIN/schksy.f index 0d78f87..681aca6 100644 --- a/TESTING/LIN/schksy.f +++ b/TESTING/LIN/schksy.f @@ -405,12 +405,8 @@ 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. * diff --git a/TESTING/LIN/zchkhe.f b/TESTING/LIN/zchkhe.f index e6524ca..a030a1d 100644 --- a/TESTING/LIN/zchkhe.f +++ b/TESTING/LIN/zchkhe.f @@ -409,12 +409,8 @@ 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. * diff --git a/TESTING/LIN/zchksy.f b/TESTING/LIN/zchksy.f index 58911a9..a82c454 100644 --- a/TESTING/LIN/zchksy.f +++ b/TESTING/LIN/zchksy.f @@ -415,12 +415,8 @@ 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. * -- 2.7.4