Add SYCONV routine: convert back and forth the factorization returned by SYTRF to be able to call SYTRS2.
Modify SYSV that now is calling SYTRS2 instead of SYTRS (and also SYCONV to convert and revert the factorization returned by SYTRF).
Modify testing to have TRS but also TRS2 tested in the LIN testing for SY.
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 \
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 \
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 \
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 \
--- /dev/null
+ 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
EXTERNAL ILAENV, LSAME
* ..
* .. External Subroutines ..
- EXTERNAL CSYTRF, CSYTRS, XERBLA
+ EXTERNAL CSYCONV, CSYTRF, CSYTRS2, XERBLA
* ..
* .. Intrinsic Functions ..
INTRINSIC MAX
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
*
--- /dev/null
+ 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
--- /dev/null
+ 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
* -- 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
*
* .. Local Scalars ..
LOGICAL LQUERY
- INTEGER LWKOPT, NB
+ INTEGER LWKOPT, NB, IINFO
* ..
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME, ILAENV
* ..
* .. External Subroutines ..
- EXTERNAL DSYTRF, DSYTRS, XERBLA
+ EXTERNAL DSYCONV, DSYTRF, DSYTRS2, XERBLA
* ..
* .. Intrinsic Functions ..
INTRINSIC MAX
* 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
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
*
--- /dev/null
+ 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
--- /dev/null
+ 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
* -- 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
*
* .. Local Scalars ..
LOGICAL LQUERY
- INTEGER LWKOPT, NB
+ INTEGER LWKOPT, NB, IINFO
* ..
* .. External Functions ..
LOGICAL LSAME
EXTERNAL ILAENV, LSAME
* ..
* .. External Subroutines ..
- EXTERNAL SSYTRF, SSYTRS, XERBLA
+ EXTERNAL SSYCONV, SSYTRF, SSYTRS2, XERBLA
* ..
* .. Intrinsic Functions ..
INTRINSIC MAX
* 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
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
*
--- /dev/null
+ 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
--- /dev/null
+ 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
*
* .. Local Scalars ..
LOGICAL LQUERY
- INTEGER LWKOPT, NB
+ INTEGER LWKOPT, NB, IINFO
* ..
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME, ILAENV
* ..
* .. External Subroutines ..
- EXTERNAL XERBLA, ZSYTRF, ZSYTRS
+ EXTERNAL ZSYCONV, XERBLA, ZSYTRF, ZSYTRS2
* ..
* .. Intrinsic Functions ..
INTRINSIC MAX
* 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
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
*
--- /dev/null
+ 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
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
* Purpose
* =======
*
-* CCHKSY tests CSYTRF, -TRI, -TRS, -RFS, and -CON.
+* CCHKSY tests CSYTRF, -TRI, -TRS, -TRS2, -RFS, and -CON.
*
* Arguments
* =========
INTEGER NTYPES
PARAMETER ( NTYPES = 11 )
INTEGER NTESTS
- PARAMETER ( NTESTS = 8 )
+ PARAMETER ( NTESTS = 9 )
* ..
* .. Local Scalars ..
LOGICAL TRFCON, ZEROT
$ LDA, RWORK, RESULT( 3 ) )
*
*+ TEST 4
+* Solve and compute residual for A * X = B.
+*
+ SRNAMT = 'CLARHS'
+ CALL CLARHS( PATH, XTYPE, UPLO, ' ', N, N, KL, KU,
+ $ NRHS, A, LDA, XACT, LDA, B, LDA,
+ $ ISEED, INFO )
+ CALL CLACPY( 'Full', N, NRHS, B, LDA, X, LDA )
+*
+ SRNAMT = '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'
$ NERRS, NOUT )
*
CALL CGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
- $ RESULT( 5 ) )
+ $ RESULT( 6 ) )
CALL CPOT05( UPLO, N, NRHS, A, LDA, B, LDA, X, LDA,
$ XACT, LDA, RWORK, RWORK( NRHS+1 ),
- $ RESULT( 6 ) )
+ $ RESULT( 7 ) )
*
* Print information about the tests that did not pass
* the threshold.
*
- DO 120 K = 3, 7
+ DO 120 K = 3, 8
IF( RESULT( K ).GE.THRESH ) THEN
IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
$ CALL ALAHD( NOUT, PATH )
NRUN = NRUN + 5
130 CONTINUE
*
-*+ TEST 8
+*+ TEST 9
* Get an estimate of RCOND = 1/CNDNUM.
*
140 CONTINUE
$ CALL ALAERH( PATH, '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
* Purpose
* =======
*
-* DCHKSY tests DSYTRF, -TRI, -TRS, -RFS, and -CON.
+* DCHKSY tests DSYTRF, -TRI, -TRS, -TRS2, -RFS, and -CON.
*
* Arguments
* =========
INTEGER NTYPES
PARAMETER ( NTYPES = 10 )
INTEGER NTESTS
- PARAMETER ( NTESTS = 8 )
+ PARAMETER ( NTESTS = 9 )
* ..
* .. Local Scalars ..
LOGICAL TRFCON, ZEROT
* .. 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
DO 130 IRHS = 1, NNS
NRHS = NSVAL( IRHS )
*
-*+ TEST 3
+*+ TEST 3 ( Using TRS)
* Solve and compute residual for A * X = B.
*
SRNAMT = 'DLARHS'
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'
$ 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 )
NRUN = NRUN + 5
130 CONTINUE
*
-*+ TEST 8
+*+ TEST 9
* Get an estimate of RCOND = 1/CNDNUM.
*
140 CONTINUE
$ 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
* Purpose
* =======
*
-* SCHKSY tests SSYTRF, -TRI, -TRS, -RFS, and -CON.
+* SCHKSY tests SSYTRF, -TRI, -TRS, -TRS2, -RFS, and -CON.
*
* Arguments
* =========
INTEGER NTYPES
PARAMETER ( NTYPES = 10 )
INTEGER NTESTS
- PARAMETER ( NTESTS = 8 )
+ PARAMETER ( NTESTS = 9 )
* ..
* .. Local Scalars ..
LOGICAL TRFCON, ZEROT
* .. 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
DO 130 IRHS = 1, NNS
NRHS = NSVAL( IRHS )
*
-*+ TEST 3
+*+ TEST 3 (Using DSYTRS)
* Solve and compute residual for A * X = B.
*
SRNAMT = 'SLARHS'
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'
$ 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 )
NRUN = NRUN + 5
130 CONTINUE
*
-*+ TEST 8
+*+ TEST 9
* Get an estimate of RCOND = 1/CNDNUM.
*
140 CONTINUE
$ 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
* Purpose
* =======
*
-* ZCHKSY tests ZSYTRF, -TRI, -TRS, -RFS, and -CON.
+* ZCHKSY tests ZSYTRF, -TRI, -TRS, -TRS2, -RFS, and -CON.
*
* Arguments
* =========
INTEGER NTYPES
PARAMETER ( NTYPES = 11 )
INTEGER NTESTS
- PARAMETER ( NTESTS = 8 )
+ PARAMETER ( NTESTS = 9 )
* ..
* .. Local Scalars ..
LOGICAL TRFCON, ZEROT
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
DO 130 IRHS = 1, NNS
NRHS = NSVAL( IRHS )
*
-*+ TEST 3
+*+ TEST 3 (Using ZSYTRS)
* Solve and compute residual for A * X = B.
*
SRNAMT = 'ZLARHS'
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'
$ NERRS, NOUT )
*
CALL ZGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
- $ RESULT( 5 ) )
+ $ RESULT( 6 ) )
CALL ZPOT05( UPLO, N, NRHS, A, LDA, B, LDA, X, LDA,
$ XACT, LDA, RWORK, RWORK( NRHS+1 ),
- $ RESULT( 6 ) )
+ $ RESULT( 7 ) )
*
* Print information about the tests that did not pass
* the threshold.
*
- DO 120 K = 3, 7
+ DO 120 K = 3, 8
IF( RESULT( K ).GE.THRESH ) THEN
IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
$ CALL ALAHD( NOUT, PATH )
NRUN = NRUN + 5
130 CONTINUE
*
-*+ TEST 8
+*+ TEST 9
* Get an estimate of RCOND = 1/CNDNUM.
*
140 CONTINUE
$ CALL ALAERH( PATH, '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
projects / platform / upstream / lapack.git / commitdiff