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 ssytrs2.o ssyconv.o stbcon.o \
+ ssytd2.o ssytf2.o ssytrd.o ssytrf.o ssytri.o ssytri2.o ssytri2x.o \
+ ssyswapr.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 \
crot.o cspcon.o cspmv.o cspr.o csprfs.o cspsv.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 csytrs2.o csyconv.o ctbcon.o ctbrfs.o ctbtrs.o ctgevc.o ctgex2.o \
+ csyr.o csyrfs.o csysv.o csysvx.o csytf2.o csytrf.o csytri.o csytri2.o csytri2x.o \
+ csyswapr.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 \
zrot.o zspcon.o zspmv.o zspr.o zsprfs.o zspsv.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 zsytrs2.o zsyconv.o ztbcon.o ztbrfs.o ztbtrs.o ztgevc.o ztgex2.o \
+ zsyr.o zsyrfs.o zsysv.o zsysvx.o zsytf2.o zsytrf.o zsytri.o zsytri2.o zsytri2x.o \
+ zsyswapr.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 CSYSWAPR( 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 ..
+ COMPLEX A(N,N)
+*
+* Purpose
+* =======
+*
+* CSYSWAPR 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) COMPLEX 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 CSYTRF.
+*
+* 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
+ COMPLEX TMP
+*
+* .. External Functions ..
+ LOGICAL LSAME
+ EXTERNAL LSAME
+* ..
+* .. External Subroutines ..
+ EXTERNAL CSWAP
+* ..
+* .. Executable Statements ..
+*
+ UPPER = LSAME( UPLO, 'U' )
+ IF (UPPER) THEN
+*
+* UPPER
+* first swap
+* - swap column I1 and I2 from I1 to I1-1
+ CALL CSWAP( 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 CSWAP ( 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 CSYSWAPR
+
--- /dev/null
+ SUBROUTINE CSYTRI2( 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( * )
+ COMPLEX A( LDA, * ), WORK( * )
+* ..
+*
+* Purpose
+* =======
+*
+* CSYTRI2 computes the inverse of a complex symmetric indefinite matrix
+* A using the factorization A = U*D*U**T or A = L*D*L**T computed by
+* CSYTRF. CSYTRI2 set the LEADING DIMENSION of the workspace
+* before calling CSYTRI2X 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) COMPLEX 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 CSYTRF.
+*
+* 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 CSYTRF.
+*
+* WORK (workspace) COMPLEX 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 CSYTRI2X
+* ..
+* .. Executable Statements ..
+*
+* Test the input parameters.
+*
+ INFO = 0
+ UPPER = LSAME( UPLO, 'U' )
+ LQUERY = ( LWORK.EQ.-1 )
+* Get blocksize
+ NBMAX = ILAENV( 1, 'CSYTRF', 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( 'CSYTRI2', -INFO )
+ RETURN
+ ELSE IF( LQUERY ) THEN
+ WORK(1)=(N+NBMAX+1)*(NBMAX+3)
+ RETURN
+ END IF
+ IF( N.EQ.0 )
+ $ RETURN
+
+ CALL CSYTRI2X( UPLO, N, A, LDA, IPIV, WORK, NBMAX, INFO )
+ RETURN
+*
+* End of CSYTRI2
+*
+ END
--- /dev/null
+ SUBROUTINE CSYTRI2X( 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, NB
+* ..
+* .. Array Arguments ..
+ INTEGER IPIV( * )
+ COMPLEX A( LDA, * ), WORK( N+NB+1,* )
+* ..
+*
+* Purpose
+* =======
+*
+* CSYTRI2X 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
+* CSYTRF.
+*
+* 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) COMPLEX 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 CSYTRF.
+*
+* 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 CSYTRF.
+*
+* WORK (workspace) COMPLEX 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 ..
+ COMPLEX ONE, ZERO
+ PARAMETER ( ONE = ( 1.0E+0, 0.0E+0 ),
+ $ ZERO = ( 0.0E+0, 0.0E+0 ) )
+* ..
+* .. Local Scalars ..
+ LOGICAL UPPER
+ INTEGER I, IINFO, IP, K, CUT, NNB
+ INTEGER COUNT
+ INTEGER J, U11, INVD
+
+ COMPLEX AK, AKKP1, AKP1, D, T
+ COMPLEX U01_I_J, U01_IP1_J
+ COMPLEX U11_I_J, U11_IP1_J
+* ..
+* .. External Functions ..
+ LOGICAL LSAME
+ EXTERNAL LSAME
+* ..
+* .. External Subroutines ..
+ EXTERNAL CSYSCONV, XERBLA, CTRTRI
+ EXTERNAL CGEMM, CTRMM, CSYSWAPR
+* ..
+* .. Intrinsic Functions ..
+ INTRINSIC 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( 'CSYTRI2X', -INFO )
+ RETURN
+ END IF
+ IF( N.EQ.0 )
+ $ RETURN
+*
+* Convert A
+* Workspace got Non-diag elements of D
+*
+ CALL CSYCONV( 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
+
+ IF( UPPER ) THEN
+*
+* invA = P * inv(U')*inv(D)*inv(U)*P'.
+*
+ CALL CTRTRI( 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) = ONE / A( K, K )
+ WORK(K,INVD+1) = 0
+ K=K+1
+ ELSE
+* 2 x 2 diagonal NNB
+ T = 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
+*
+ CALL CTRMM('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 CGEMM('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 CTRMM('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 CSYSWAPR( UPLO, N, A, I ,IP )
+ IF (I .GT. IP) CALL CSYSWAPR( UPLO, N, A, IP ,I )
+ ELSE
+ IP=-IPIV(I)
+ I=I+1
+ IF ( (I-1) .LT. IP)
+ $ CALL CSYSWAPR( UPLO, N, A, I-1 ,IP )
+ IF ( (I-1) .GT. IP)
+ $ CALL CSYSWAPR( 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 CTRTRI( 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) = ONE / A( K, K )
+ WORK(K,INVD+1) = 0
+ K=K-1
+ ELSE
+* 2 x 2 diagonal NNB
+ T = 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 CTRMM('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 CGEMM('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 CTRMM('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 CSYSWAPR( UPLO, N, A, I ,IP )
+ IF (I .GT. IP) CALL CSYSWAPR( UPLO, N, A, IP ,I )
+ ELSE
+ IP=-IPIV(I)
+ IF ( I .LT. IP) CALL CSYSWAPR( UPLO, N, A, I ,IP )
+ IF ( I .GT. IP) CALL CSYSWAPR( UPLO, N, A, IP ,I )
+ I=I-1
+ ENDIF
+ I=I-1
+ END DO
+ END IF
+*
+ RETURN
+*
+* End of CSYTRI2X
+*
+ END
+
* 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
*
* .. Scalar Arguments ..
CHARACTER UPLO
- INTEGER INFO, LDA, N
+ INTEGER INFO, LDA, N, NB
* ..
* .. Array Arguments ..
INTEGER IPIV( * )
EXTERNAL LSAME
* ..
* .. External Subroutines ..
- EXTERNAL DSCAL, DSYSCONV, DSWAP, XERBLA, DTRTRI
+ EXTERNAL DSYSCONV, XERBLA, DTRTRI
EXTERNAL DGEMM, DTRMM, DSYSWAPR
* ..
* .. Intrinsic Functions ..
- INTRINSIC ABS, MAX
+ INTRINSIC MAX
* ..
* .. Executable Statements ..
*
* 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
*
K=K+1
ELSE
* 2 x 2 diagonal NNB
- T = ABS( WORK(K+1,1) )
+ T = WORK(K+1,1)
AK = A( K, K ) / T
AKP1 = A( K+1, K+1 ) / T
AKKP1 = WORK(K+1,1) / T
*
* 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)
*
K=K-1
ELSE
* 2 x 2 diagonal NNB
- T = ABS(WORK(K-1,1))
+ T = WORK(K-1,1)
AK = A( K-1, K-1 ) / T
AKP1 = A( K, K ) / T
AKKP1 = WORK(K-1,1) / T
* 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
--- /dev/null
+ SUBROUTINE SSYSWAPR( 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 ..
+ REAL A(N,N)
+*
+* Purpose
+* =======
+*
+* SSYSWAPR 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) REAL 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 SSYTRF.
+*
+* 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
+ REAL TMP
+*
+* .. External Functions ..
+ LOGICAL LSAME
+ EXTERNAL LSAME
+* ..
+* .. External Subroutines ..
+ EXTERNAL SSWAP
+* ..
+* .. Executable Statements ..
+*
+ UPPER = LSAME( UPLO, 'U' )
+ IF (UPPER) THEN
+*
+* UPPER
+* first swap
+* - swap column I1 and I2 from I1 to I1-1
+ CALL SSWAP( 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 SSWAP( 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 SSYSWAPR
+
--- /dev/null
+ SUBROUTINE SSYTRI2( 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( * )
+ REAL A( LDA, * ), WORK( * )
+* ..
+*
+* Purpose
+* =======
+*
+* SSYTRI2 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
+* SSYTRF. SSYTRI2 set the LEADING DIMENSION of the workspace
+* before calling SSYTRI2X 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) REAL 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 SSYTRF.
+*
+* 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 SSYTRF.
+*
+* WORK (workspace) REAL 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 SSYTRI2X
+* ..
+* .. Executable Statements ..
+*
+* Test the input parameters.
+*
+ INFO = 0
+ UPPER = LSAME( UPLO, 'U' )
+ LQUERY = ( LWORK.EQ.-1 )
+* Get blocksize
+ NBMAX = ILAENV( 1, 'SSYTRF', 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( 'SSYTRI2', -INFO )
+ RETURN
+ ELSE IF( LQUERY ) THEN
+ WORK(1)=(N+NBMAX+1)*(NBMAX+3)
+ RETURN
+ END IF
+ IF( N.EQ.0 )
+ $ RETURN
+
+ CALL SSYTRI2X( UPLO, N, A, LDA, IPIV, WORK, NBMAX, INFO )
+ RETURN
+*
+* End of SSYTRI2
+*
+ END
--- /dev/null
+ SUBROUTINE SSYTRI2X( 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, NB
+* ..
+* .. Array Arguments ..
+ INTEGER IPIV( * )
+ REAL A( LDA, * ), WORK( N+NB+1,* )
+* ..
+*
+* Purpose
+* =======
+*
+* SSYTRI2X 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
+* SSYTRF.
+*
+* 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) REAL 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 SSYTRF.
+*
+* 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 SSYTRF.
+*
+* WORK (workspace) REAL 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 ..
+ REAL ONE, ZERO
+ PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 )
+* ..
+* .. Local Scalars ..
+ LOGICAL UPPER
+ INTEGER I, IINFO, IP, K, CUT, NNB
+ INTEGER COUNT
+ INTEGER J, U11, INVD
+
+ REAL AK, AKKP1, AKP1, D, T
+ REAL U01_I_J, U01_IP1_J
+ REAL U11_I_J, U11_IP1_J
+* ..
+* .. External Functions ..
+ LOGICAL LSAME
+ EXTERNAL LSAME
+* ..
+* .. External Subroutines ..
+ EXTERNAL SSYSCONV, XERBLA, STRTRI
+ EXTERNAL SGEMM, STRMM, SSYSWAPR
+* ..
+* .. Intrinsic Functions ..
+ INTRINSIC 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( 'SSYTRI2X', -INFO )
+ RETURN
+ END IF
+ IF( N.EQ.0 )
+ $ RETURN
+*
+* Convert A
+* Workspace got Non-diag elements of D
+*
+ CALL SSYCONV( 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
+
+ IF( UPPER ) THEN
+*
+* invA = P * inv(U')*inv(D)*inv(U)*P'.
+*
+ CALL STRTRI( 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 = 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
+*
+ CALL STRMM('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 SGEMM('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 STRMM('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 SSYSWAPR( UPLO, N, A, I ,IP )
+ IF (I .GT. IP) CALL SSYSWAPR( UPLO, N, A, IP ,I )
+ ELSE
+ IP=-IPIV(I)
+ I=I+1
+ IF ( (I-1) .LT. IP)
+ $ CALL SSYSWAPR( UPLO, N, A, I-1 ,IP )
+ IF ( (I-1) .GT. IP)
+ $ CALL SSYSWAPR( 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 STRTRI( 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 = 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 STRMM('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 SGEMM('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 STRMM('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 SSYSWAPR( UPLO, N, A, I ,IP )
+ IF (I .GT. IP) CALL SSYSWAPR( UPLO, N, A, IP ,I )
+ ELSE
+ IP=-IPIV(I)
+ IF ( I .LT. IP) CALL SSYSWAPR( UPLO, N, A, I ,IP )
+ IF ( I .GT. IP) CALL SSYSWAPR( UPLO, N, A, IP ,I )
+ I=I-1
+ ENDIF
+ I=I-1
+ END DO
+ END IF
+*
+ RETURN
+*
+* End of SSYTRI2X
+*
+ END
+
*
* .. Local Scalars ..
LOGICAL LQUERY
- INTEGER IINFO, LWKOPT, NB
+ INTEGER LWKOPT, NB
* ..
* .. External Functions ..
LOGICAL LSAME
* 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
--- /dev/null
+ SUBROUTINE ZSYSWAPR( 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 COMPLEX A(N,N)
+*
+* Purpose
+* =======
+*
+* ZSYSWAPR 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 COMPLEX 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 ZSYTRF.
+*
+* 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 COMPLEX TMP
+*
+* .. External Functions ..
+ LOGICAL LSAME
+ EXTERNAL LSAME
+* ..
+* .. External Subroutines ..
+ EXTERNAL ZSWAP
+* ..
+* .. Executable Statements ..
+*
+ UPPER = LSAME( UPLO, 'U' )
+ IF (UPPER) THEN
+*
+* UPPER
+* first swap
+* - swap column I1 and I2 from I1 to I1-1
+ CALL ZSWAP( 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 ZSWAP( 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 ZSYSWAPR
+
--- /dev/null
+ SUBROUTINE ZSYTRI2( 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 COMPLEX A( LDA, * ), WORK( * )
+* ..
+*
+* Purpose
+* =======
+*
+* ZSYTRI2 computes the inverse of a complex symmetric indefinite matrix
+* A using the factorization A = U*D*U**T or A = L*D*L**T computed by
+* ZSYTRF. ZSYTRI2 set the LEADING DIMENSION of the workspace
+* before calling ZSYTRI2X 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 COMPLEX 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 ZSYTRF.
+*
+* 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 ZSYTRF.
+*
+* WORK (workspace) DOUBLE COMPLEX 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 ZSYTRI2X
+* ..
+* .. Executable Statements ..
+*
+* Test the input parameters.
+*
+ INFO = 0
+ UPPER = LSAME( UPLO, 'U' )
+ LQUERY = ( LWORK.EQ.-1 )
+* Get blocksize
+ NBMAX = ILAENV( 1, 'ZSYTRF', 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( 'ZSYTRI2', -INFO )
+ RETURN
+ ELSE IF( LQUERY ) THEN
+ WORK(1)=(N+NBMAX+1)*(NBMAX+3)
+ RETURN
+ END IF
+ IF( N.EQ.0 )
+ $ RETURN
+
+ CALL ZSYTRI2X( UPLO, N, A, LDA, IPIV, WORK, NBMAX, INFO )
+ RETURN
+*
+* End of ZSYTRI2
+*
+ END
--- /dev/null
+ SUBROUTINE ZSYTRI2X( 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, NB
+* ..
+* .. Array Arguments ..
+ INTEGER IPIV( * )
+ DOUBLE COMPLEX A( LDA, * ), WORK( N+NB+1,* )
+* ..
+*
+* Purpose
+* =======
+*
+* ZSYTRI2X computes the inverse of a complex symmetric indefinite matrix
+* A using the factorization A = U*D*U**T or A = L*D*L**T computed by
+* ZSYTRF.
+*
+* 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 COMPLEX 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 ZSYTRF.
+*
+* 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 ZSYTRF.
+*
+* WORK (workspace) DOUBLE COMPLEX 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 COMPLEX ONE, ZERO
+ PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ),
+ $ ZERO = ( 0.0D+0, 0.0D+0 ) )
+* ..
+* .. Local Scalars ..
+ LOGICAL UPPER
+ INTEGER I, IINFO, IP, K, CUT, NNB
+ INTEGER COUNT
+ INTEGER J, U11, INVD
+
+ DOUBLE COMPLEX AK, AKKP1, AKP1, D, T
+ DOUBLE COMPLEX U01_I_J, U01_IP1_J
+ DOUBLE COMPLEX U11_I_J, U11_IP1_J
+* ..
+* .. External Functions ..
+ LOGICAL LSAME
+ EXTERNAL LSAME
+* ..
+* .. External Subroutines ..
+ EXTERNAL ZSYSCONV, XERBLA, ZTRTRI
+ EXTERNAL ZGEMM, ZTRMM, ZSYSWAPR
+* ..
+* .. Intrinsic Functions ..
+ INTRINSIC 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( 'ZSYTRI2X', -INFO )
+ RETURN
+ END IF
+ IF( N.EQ.0 )
+ $ RETURN
+*
+* Convert A
+* Workspace got Non-diag elements of D
+*
+ CALL ZSYCONV( 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
+
+ IF( UPPER ) THEN
+*
+* invA = P * inv(U')*inv(D)*inv(U)*P'.
+*
+ CALL ZTRTRI( 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 = 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
+*
+ CALL ZTRMM('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 ZGEMM('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 ZTRMM('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 ZSYSWAPR( UPLO, N, A, I ,IP )
+ IF (I .GT. IP) CALL ZSYSWAPR( UPLO, N, A, IP ,I )
+ ELSE
+ IP=-IPIV(I)
+ I=I+1
+ IF ( (I-1) .LT. IP)
+ $ CALL ZSYSWAPR( UPLO, N, A, I-1 ,IP )
+ IF ( (I-1) .GT. IP)
+ $ CALL ZSYSWAPR( 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 ZTRTRI( 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 = 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 ZTRMM('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 ZGEMM('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 ZTRMM('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 ZSYSWAPR( UPLO, N, A, I ,IP )
+ IF (I .GT. IP) CALL ZSYSWAPR( UPLO, N, A, IP ,I )
+ ELSE
+ IP=-IPIV(I)
+ IF ( I .LT. IP) CALL ZSYSWAPR( UPLO, N, A, I ,IP )
+ IF ( I .GT. IP) CALL ZSYSWAPR( UPLO, N, A, IP ,I )
+ I=I-1
+ ENDIF
+ I=I-1
+ END DO
+ END IF
+*
+ RETURN
+*
+* End of ZSYTRI2X
+*
+ END
+
* Purpose
* =======
*
-* CCHKSY tests CSYTRF, -TRI, -TRS, -TRS2, -RFS, and -CON.
+* CCHKSY tests CSYTRF, -TRI2, -TRS, -TRS2, -RFS, and -CON.
*
* Arguments
* =========
* .. External Subroutines ..
EXTERNAL ALAERH, ALAHD, ALASUM, CERRSY, CGET04, CLACPY,
$ CLARHS, CLATB4, CLATMS, CLATSY, CPOT05, CSYCON,
- $ CSYRFS, CSYT01, CSYT02, CSYT03, CSYTRF, CSYTRI,
+ $ CSYRFS, CSYT01, CSYT02, CSYT03, CSYTRF, CSYTRI2,
$ CSYTRS, XLAENV
* ..
* .. Intrinsic Functions ..
*
IF( INB.EQ.1 .AND. .NOT.TRFCON ) THEN
CALL CLACPY( UPLO, N, N, AFAC, LDA, AINV, LDA )
- SRNAMT = 'CSYTRI'
- CALL CSYTRI( UPLO, N, AINV, LDA, IWORK, WORK,
- $ INFO )
+ SRNAMT = 'CSYTRI2'
+ LWORK = (N+NB+1)*(NB+3)
+ CALL CSYTRI2( UPLO, N, AINV, LDA, IWORK, WORK,
+ $ LWORK, INFO )
*
* Check error code from CSYTRI.
*
* .. External Subroutines ..
EXTERNAL ALADHD, ALAERH, ALASVM, CERRVX, CGET04, CLACPY,
$ CLARHS, CLASET, CLATB4, CLATMS, CLATSY, CPOT05,
- $ CSYSV, CSYSVX, CSYT01, CSYT02, CSYTRF, CSYTRI,
+ $ CSYSV, CSYSVX, CSYT01, CSYT02, CSYTRF, CSYTRI2,
$ XLAENV
* ..
* .. Scalars in Common ..
* Compute inv(A) and take its norm.
*
CALL CLACPY( UPLO, N, N, AFAC, LDA, AINV, LDA )
- CALL CSYTRI( UPLO, N, AINV, LDA, IWORK, WORK,
- $ INFO )
+ LWORK = (N+NB+1)*(NB+3)
+ CALL CSYTRI2( UPLO, N, AINV, LDA, IWORK, WORK,
+ $ LWORK, INFO )
AINVNM = CLANSY( '1', UPLO, N, AINV, LDA, RWORK )
*
* Compute the 1-norm condition number of A.
* .. External Subroutines ..
EXTERNAL ALADHD, ALAERH, ALASVM, CERRVX, CGET04, CLACPY,
$ CLARHS, CLASET, CLATB4, CLATMS, CLATSY, CPOT05,
- $ CSYSV, CSYSVX, CSYT01, CSYT02, CSYTRF, CSYTRI,
+ $ CSYSV, CSYSVX, CSYT01, CSYT02, CSYTRF, CSYTRI2,
$ XLAENV, CSYSVXX
* ..
* .. Scalars in Common ..
* Compute inv(A) and take its norm.
*
CALL CLACPY( UPLO, N, N, AFAC, LDA, AINV, LDA )
- CALL CSYTRI( UPLO, N, AINV, LDA, IWORK, WORK,
- $ INFO )
+ LWORK = (N+NB+1)*(NB+3)
+ CALL CSYTRI2( UPLO, N, AINV, LDA, IWORK, WORK,
+ $ LWORK, INFO )
AINVNM = CLANSY( '1', UPLO, N, AINV, LDA, RWORK )
*
* Compute the 1-norm condition number of A.
* .. External Subroutines ..
EXTERNAL ALAESM, CHKXER, CSPCON, CSPRFS, CSPTRF, CSPTRI,
$ CSPTRS, CSYCON, CSYRFS, CSYTF2, CSYTRF, CSYTRI,
- $ CSYTRS
+ $ CSYTRI2, CSYTRS
* ..
* .. Scalars in Common ..
LOGICAL LERR, OK
CALL CSYTRI( 'U', 2, A, 1, IP, W, INFO )
CALL CHKXER( 'CSYTRI', INFOT, NOUT, LERR, OK )
*
+* CSYTRI2
+*
+ SRNAMT = 'CSYTRI2'
+ INFOT = 1
+ CALL CSYTRI2( '/', 0, A, 1, IP, W, 1, INFO )
+ CALL CHKXER( 'CSYTRI2', INFOT, NOUT, LERR, OK )
+ INFOT = 2
+ CALL CSYTRI2( 'U', -1, A, 1, IP, W, 1, INFO )
+ CALL CHKXER( 'CSYTRI2', INFOT, NOUT, LERR, OK )
+ INFOT = 4
+ CALL CSYTRI2( 'U', 2, A, 1, IP, W, 1, INFO )
+ CALL CHKXER( 'CSYTRI2', INFOT, NOUT, LERR, OK )
+*
* CSYTRS
*
SRNAMT = 'CSYTRS'
* .. External Subroutines ..
EXTERNAL ALAESM, CHKXER, CSPCON, CSPRFS, CSPTRF, CSPTRI,
$ CSPTRS, CSYCON, CSYRFS, CSYTF2, CSYTRF, CSYTRI,
- $ CSYTRS, CSYRFSX
+ $ CSYTRI2, CSYTRS, CSYRFSX
* ..
* .. Scalars in Common ..
LOGICAL LERR, OK
CALL CSYTRI( 'U', 2, A, 1, IP, W, INFO )
CALL CHKXER( 'CSYTRI', INFOT, NOUT, LERR, OK )
*
+* CSYTRI2
+*
+ SRNAMT = 'CSYTRI2'
+ INFOT = 1
+ CALL CSYTRI2( '/', 0, A, 1, IP, W, IW, INFO )
+ CALL CHKXER( 'CSYTRI2', INFOT, NOUT, LERR, OK )
+ INFOT = 2
+ CALL CSYTRI2( 'U', -1, A, 1, IP, W, IW, INFO )
+ CALL CHKXER( 'CSYTRI2', INFOT, NOUT, LERR, OK )
+ INFOT = 4
+ CALL CSYTRI2( 'U', 2, A, 1, IP, W, IW, INFO )
+ CALL CHKXER( 'CSYTRI2', INFOT, NOUT, LERR, OK )
+*
* CSYTRS
*
SRNAMT = 'CSYTRS'
* Purpose
* =======
*
-* SCHKSY tests SSYTRF, -TRI, -TRS, -TRS2, -RFS, and -CON.
+* SCHKSY tests SSYTRF, -TRI2, -TRS, -TRS2, -RFS, and -CON.
*
* Arguments
* =========
EXTERNAL ALAERH, ALAHD, ALASUM, SERRSY, SGET04, SLACPY,
$ SLARHS, SLATB4, SLATMS, SPOT02, SPOT03, SPOT05,
$ SSYCON, SSYCONV, SSYRFS, SSYT01, SSYTRF,
- $ SSYTRI, SSYTRS, SSYTRS2, XLAENV
+ $ SSYTRI2, SSYTRS, SSYTRS2, XLAENV
* ..
* .. Intrinsic Functions ..
INTRINSIC MAX, MIN
*
IF( INB.EQ.1 .AND. .NOT.TRFCON ) THEN
CALL SLACPY( UPLO, N, N, AFAC, LDA, AINV, LDA )
- SRNAMT = 'SSYTRI'
- CALL SSYTRI( UPLO, N, AINV, LDA, IWORK, WORK,
- $ INFO )
+ SRNAMT = 'SSYTRI2'
+ LWORK = (N+NB+1)*(NB+3)
+ CALL SSYTRI2( UPLO, N, AINV, LDA, IWORK, WORK,
+ $ LWORK, INFO )
*
-* Check error code from SSYTRI.
+* Check error code from SSYTRI2.
*
IF( INFO.NE.0 )
- $ CALL ALAERH( PATH, 'SSYTRI', INFO, -1, UPLO, N,
+ $ CALL ALAERH( PATH, 'SSYTRI2', INFO, -1, UPLO, N,
$ N, -1, -1, -1, IMAT, NFAIL, NERRS,
$ NOUT )
*
* .. External Subroutines ..
EXTERNAL ALADHD, ALAERH, ALASVM, SERRVX, SGET04, SLACPY,
$ SLARHS, SLASET, SLATB4, SLATMS, SPOT02, SPOT05,
- $ SSYSV, SSYSVX, SSYT01, SSYTRF, SSYTRI, XLAENV
+ $ SSYSV, SSYSVX, SSYT01, SSYTRF, SSYTRI2, XLAENV
* ..
* .. Scalars in Common ..
LOGICAL LERR, OK
* Compute inv(A) and take its norm.
*
CALL SLACPY( UPLO, N, N, AFAC, LDA, AINV, LDA )
- CALL SSYTRI( UPLO, N, AINV, LDA, IWORK, WORK,
- $ INFO )
+ LWORK = (N+NB+1)*(NB+3)
+ CALL SSYTRI2( UPLO, N, AINV, LDA, IWORK, WORK,
+ $ LWORK, INFO )
AINVNM = SLANSY( '1', UPLO, N, AINV, LDA, RWORK )
*
* Compute the 1-norm condition number of A.
* .. External Subroutines ..
EXTERNAL ALADHD, ALAERH, ALASVM, SERRVX, SGET04, SLACPY,
$ SLARHS, SLASET, SLATB4, SLATMS, SPOT02, SPOT05,
- $ SSYSV, SSYSVX, SSYT01, SSYTRF, SSYTRI, XLAENV,
+ $ SSYSV, SSYSVX, SSYT01, SSYTRF, SSYTRI2, XLAENV,
$ SSYSVXX
* ..
* .. Scalars in Common ..
* Compute inv(A) and take its norm.
*
CALL SLACPY( UPLO, N, N, AFAC, LDA, AINV, LDA )
- CALL SSYTRI( UPLO, N, AINV, LDA, IWORK, WORK,
- $ INFO )
+ LWORK = (N+NB+1)*(NB+3)
+ CALL SSYTRI2( UPLO, N, AINV, LDA, IWORK, WORK,
+ $ LWORK, INFO )
AINVNM = SLANSY( '1', UPLO, N, AINV, LDA, RWORK )
*
* Compute the 1-norm condition number of A.
* .. External Subroutines ..
EXTERNAL ALAESM, CHKXER, SSPCON, SSPRFS, SSPTRF, SSPTRI,
$ SSPTRS, SSYCON, SSYRFS, SSYTF2, SSYTRF, SSYTRI,
- $ SSYTRS
+ $ SSYTRI2, SSYTRS
* ..
* .. Scalars in Common ..
LOGICAL LERR, OK
CALL SSYTRI( 'U', 2, A, 1, IP, W, INFO )
CALL CHKXER( 'SSYTRI', INFOT, NOUT, LERR, OK )
*
+* SSYTRI2
+*
+ SRNAMT = 'SSYTRI2'
+ INFOT = 1
+ CALL SSYTRI2( '/', 0, A, 1, IP, W, IW, INFO )
+ CALL CHKXER( 'SSYTRI2', INFOT, NOUT, LERR, OK )
+ INFOT = 2
+ CALL SSYTRI2( 'U', -1, A, 1, IP, W, IW, INFO )
+ CALL CHKXER( 'SSYTRI2', INFOT, NOUT, LERR, OK )
+ INFOT = 4
+ CALL SSYTRI2( 'U', 2, A, 1, IP, W, IW, INFO )
+ CALL CHKXER( 'SSYTRI2', INFOT, NOUT, LERR, OK )
+*
* SSYTRS
*
SRNAMT = 'SSYTRS'
* .. External Subroutines ..
EXTERNAL ALAESM, CHKXER, SSPCON, SSPRFS, SSPTRF, SSPTRI,
$ SSPTRS, SSYCON, SSYRFS, SSYTF2, SSYTRF, SSYTRI,
- $ SSYTRS, SSYRFSX
+ $ SSYTRI2, SSYTRS, SSYRFSX
* ..
* .. Scalars in Common ..
LOGICAL LERR, OK
CALL SSYTRI( 'U', 2, A, 1, IP, W, INFO )
CALL CHKXER( 'SSYTRI', INFOT, NOUT, LERR, OK )
*
+* SSYTRI2
+*
+ SRNAMT = 'SSYTRI2'
+ INFOT = 1
+ CALL SSYTRI2( '/', 0, A, 1, IP, W, IW, INFO )
+ CALL CHKXER( 'SSYTRI', INFOT, NOUT, LERR, OK )
+ INFOT = 2
+ CALL SSYTRI2( 'U', -1, A, 1, IP, W, IW, INFO )
+ CALL CHKXER( 'SSYTRI', INFOT, NOUT, LERR, OK )
+ INFOT = 4
+ CALL SSYTRI2( 'U', 2, A, 1, IP, W, IW, INFO )
+ CALL CHKXER( 'SSYTRI', INFOT, NOUT, LERR, OK )
+*
* SSYTRS
*
SRNAMT = 'SSYTRS'
* Purpose
* =======
*
-* ZCHKSY tests ZSYTRF, -TRI, -TRS, -TRS2, -RFS, and -CON.
+* ZCHKSY tests ZSYTRF, -TRI2, -TRS, -TRS2, -RFS, and -CON.
*
* Arguments
* =========
EXTERNAL ALAERH, ALAHD, ALASUM, XLAENV, ZERRSY, ZGET04,
$ ZLACPY, ZLARHS, ZLATB4, ZLATMS, ZLATSY, ZPOT05,
$ ZSYCON, ZSYRFS, ZSYT01, ZSYT02, ZSYT03, ZSYTRF,
- $ ZSYTRI, ZSYTRS, ZSYTRS2
+ $ ZSYTRI2, ZSYTRS, ZSYTRS2
* ..
* .. Intrinsic Functions ..
INTRINSIC MAX, MIN
*
IF( INB.EQ.1 .AND. .NOT.TRFCON ) THEN
CALL ZLACPY( UPLO, N, N, AFAC, LDA, AINV, LDA )
- SRNAMT = 'ZSYTRI'
- CALL ZSYTRI( UPLO, N, AINV, LDA, IWORK, WORK,
- $ INFO )
+ SRNAMT = 'ZSYTRI2'
+ LWORK = (N+NB+1)*(NB+3)
+ CALL ZSYTRI2( UPLO, N, AINV, LDA, IWORK, WORK,
+ $ LWORK, INFO )
*
-* Check error code from ZSYTRI.
+* Check error code from ZSYTRI2.
*
IF( INFO.NE.0 )
- $ CALL ALAERH( PATH, 'ZSYTRI', INFO, 0, UPLO, N,
+ $ CALL ALAERH( PATH, 'ZSYTRI2', INFO, 0, UPLO, N,
$ N, -1, -1, -1, IMAT, NFAIL, NERRS,
$ NOUT )
*
EXTERNAL ALADHD, ALAERH, ALASVM, XLAENV, ZERRVX, ZGET04,
$ ZLACPY, ZLARHS, ZLASET, ZLATB4, ZLATMS, ZLATSY,
$ ZPOT05, ZSYSV, ZSYSVX, ZSYT01, ZSYT02, ZSYTRF,
- $ ZSYTRI
+ $ ZSYTRI2
* ..
* .. Scalars in Common ..
LOGICAL LERR, OK
* Compute inv(A) and take its norm.
*
CALL ZLACPY( UPLO, N, N, AFAC, LDA, AINV, LDA )
- CALL ZSYTRI( UPLO, N, AINV, LDA, IWORK, WORK,
- $ INFO )
+ LWORK = (N+NB+1)*(NB+3)
+ CALL ZSYTRI2( UPLO, N, AINV, LDA, IWORK, WORK,
+ $ LWORK, INFO )
AINVNM = ZLANSY( '1', UPLO, N, AINV, LDA, RWORK )
*
* Compute the 1-norm condition number of A.
EXTERNAL ALADHD, ALAERH, ALASVM, XLAENV, ZERRVX, ZGET04,
$ ZLACPY, ZLARHS, ZLASET, ZLATB4, ZLATMS, ZLATSY,
$ ZPOT05, ZSYSV, ZSYSVX, ZSYT01, ZSYT02, ZSYTRF,
- $ ZSYTRI, ZSYSVXX
+ $ ZSYTRI2, ZSYSVXX
* ..
* .. Scalars in Common ..
LOGICAL LERR, OK
* Compute inv(A) and take its norm.
*
CALL ZLACPY( UPLO, N, N, AFAC, LDA, AINV, LDA )
- CALL ZSYTRI( UPLO, N, AINV, LDA, IWORK, WORK,
- $ INFO )
+ LWORK = (N+NB+1)*(NB+3)
+ CALL ZSYTRI2( UPLO, N, AINV, LDA, IWORK, WORK,
+ $ LWORK, INFO )
AINVNM = ZLANSY( '1', UPLO, N, AINV, LDA, RWORK )
*
* Compute the 1-norm condition number of A.
* .. External Subroutines ..
EXTERNAL ALAESM, CHKXER, ZSPCON, ZSPRFS, ZSPTRF, ZSPTRI,
$ ZSPTRS, ZSYCON, ZSYRFS, ZSYTF2, ZSYTRF, ZSYTRI,
- $ ZSYTRS
+ $ ZSYTRI2, ZSYTRS
* ..
* .. Scalars in Common ..
LOGICAL LERR, OK
CALL ZSYTRI( 'U', 2, A, 1, IP, W, INFO )
CALL CHKXER( 'ZSYTRI', INFOT, NOUT, LERR, OK )
*
+* ZSYTRI2
+*
+ SRNAMT = 'ZSYTRI2'
+ INFOT = 1
+ CALL ZSYTRI2( '/', 0, A, 1, IP, W, 1, INFO )
+ CALL CHKXER( 'ZSYTRI2', INFOT, NOUT, LERR, OK )
+ INFOT = 2
+ CALL ZSYTRI2( 'U', -1, A, 1, IP, W, 1, INFO )
+ CALL CHKXER( 'ZSYTRI2', INFOT, NOUT, LERR, OK )
+ INFOT = 4
+ CALL ZSYTRI2( 'U', 2, A, 1, IP, W, 1, INFO )
+ CALL CHKXER( 'ZSYTRI2', INFOT, NOUT, LERR, OK )
+*
* ZSYTRS
*
SRNAMT = 'ZSYTRS'
* .. External Subroutines ..
EXTERNAL ALAESM, CHKXER, ZSPCON, ZSPRFS, ZSPTRF, ZSPTRI,
$ ZSPTRS, ZSYCON, ZSYRFS, ZSYTF2, ZSYTRF, ZSYTRI,
- $ ZSYTRS, ZSYRFSX
+ $ ZSYTRI2, ZSYTRS, ZSYRFSX
* ..
* .. Scalars in Common ..
LOGICAL LERR, OK
CALL ZSYTRI( 'U', 2, A, 1, IP, W, INFO )
CALL CHKXER( 'ZSYTRI', INFOT, NOUT, LERR, OK )
*
+* ZSYTRI2
+*
+ SRNAMT = 'ZSYTRI2'
+ INFOT = 1
+ CALL ZSYTRI2( '/', 0, A, 1, IP, W, 1, INFO )
+ CALL CHKXER( 'ZSYTRI2', INFOT, NOUT, LERR, OK )
+ INFOT = 2
+ CALL ZSYTRI2( 'U', -1, A, 1, IP, W, 1, INFO )
+ CALL CHKXER( 'ZSYTRI2', INFOT, NOUT, LERR, OK )
+ INFOT = 4
+ CALL ZSYTRI2( 'U', 2, A, 1, IP, W, 1, INFO )
+ CALL CHKXER( 'ZSYTRI2', INFOT, NOUT, LERR, OK )
+*
* ZSYTRS
*
SRNAMT = 'ZSYTRS'
projects / platform / upstream / lapack.git / commitdiff