************************************************************************ * Program Name: TESTEVX * ************************************************************************ * * * Description: * * * * This program demonstrates a problem with DSYEVX where the * * eigenvectors returned do not pass the math verification tests. * * * * 1) Call CM415 to create the test matrix. The problem occurs * * when the matrix is of order greater than 6. * * * * 2) When N=7, the following matrix is created: * * * * 1.0 0.0 0.0 0.0 0.0 0.0 1.0 * * 0.0 1.0 0.0 0.0 0.0 0.0 2.0 * * 0.0 0.0 1.0 0.0 0.0 0.0 3.0 * * 0.0 0.0 0.0 1.0 0.0 0.0 4.0 * * 0.0 0.0 0.0 0.0 1.0 0.0 5.0 * * 0.0 0.0 0.0 0.0 0.0 1.0 6.0 * * 1.0 2.0 3.0 4.0 5.0 6.0 7.0 * * * * 3) When DSYEVX is called with JOBZ=V and RANGE=A, then * * following eigenvalues are returned: * * * * (1) -6 * * (2) 1 * * (3) 1 * * (4) 1 * * (5) 1 * * (6) 1 * * (7) 14 * * * * 4) The eigenvectors returned with these eigenvalues pass the * * DSYT21 math checks. * * * * 5) A failure occurs when DSYEVX is called with JOBZ=V and * * RANGE=V with VL=-1 and VU=11. With these values of VL * * and VU only the 5 eigenvalues with a value of 1 will be * * returned. * * The correct eigenvalues are returned but the eigenvectors * * returned cause the DSYT22 math checks to fail. * * * ************************************************************************ PROGRAM TESTEVX IMPLICIT NONE INTEGER MAXN, NB, LWORK PARAMETER (MAXN=7,NB=64,LWORK=(NB+3)*MAXN) * -------------------------------------------------------------------- * The essential variables for calling dsyevx. * -------------------------------------------------------------------- INTEGER N,IL,IU,LDA,LDZ,M CHARACTER*1 JOBZ, RANGE, UPLO PARAMETER (LDA=MAXN,LDZ=MAXN) DOUBLE PRECISION A(LDA,MAXN),AA(LDA,MAXN) DOUBLE PRECISION Z(LDZ,MAXN) DOUBLE PRECISION W(MAXN) DOUBLE PRECISION VL,VU,ABSTOL DOUBLE PRECISION WORK(LWORK),IWORK(5*MAXN) INTEGER INFO INTEGER IFAIL(MAXN) * -------------------------------------------------------------------- * Additional variables. * -------------------------------------------------------------------- INTEGER I,J, IUNIT PARAMETER (IUNIT = 4) ! output device for printing REAL*8 THRESH PARAMETER (THRESH=100.0D0) REAL*8 RSLT21(2),RSLT22(2) REAL*8 WORKT21(2*MAXN**2), WORKT22(2*MAXN**2) * -------------------------------------------------------------------- * .. * .. Executable Statements .. N = 7 CALL CM415 (AA,LDA,N) JOBZ='V' RANGE='V' VL = -1.0D0 VU = 11.0D0 UPLO = 'L' ABSTOL = 0.0D0 CALL DLACPY(' ',N,N,AA,LDA,A,LDA) CALL DSYEVX( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL, IU, $ ABSTOL, M, W, Z, LDZ, WORK, LWORK, IWORK, $ IFAIL, INFO ) *----------------------------------------------------------------------* * Run math checks to verify the eigenvectors work and are * * orthogonal. * * If the number of eigenvalues returned is equal to N then call * * DSYT21. Otherwise, call DSYT22. * * * * DSYT21 Performs the following math checks: * * * * RESULT(1) = | A - U S U' | / ( |A| n ulp ) *and* * * RESULT(2) = | I - UU' | / ( n ulp ) * * * * DSYT22 Performs the following math checks: * * * * RESULT(1) = | U' A U - S | / ( |A| m ulp ) *and* * * RESULT(2) = | I - U'U | / ( m ulp ) * * * *----------------------------------------------------------------------* CALL DLACPY(' ',N,N,AA,LDA,A,LDA) IF (M .EQ. N) THEN CALL DSYT21(1,UPLO,N,0,A,LDA,W(1),0,Z(1,1),LDZ,0,0,0, + WORKT21,RSLT21) IF ((RSLT21(1) .GE. THRESH) .OR. (RSLT21(2) .GE. THRESH)) + THEN WRITE(IUNIT,2580) 'DSYT21',N,M,THRESH WRITE(IUNIT,2585) RSLT21(1),RSLT21(2) ENDIF ELSE CALL DSYT22(1,UPLO,N,M,0,A,LDA,W(1),0,Z(1,1),LDZ, + 0,0,0,WORKT22,RSLT22) IF ((RSLT22(1) .GE. THRESH) .OR. (RSLT22(2) .GE. THRESH)) + THEN WRITE(IUNIT,2580) 'DSYT22',N,M,THRESH WRITE(IUNIT,2595) RSLT22(1),RSLT22(2) ENDIF ENDIF 2580 FORMAT(/,6X,'*** Math check results for ',A7, + ' with N=',I4,' and M=',I4,/, + 6x,'*** exceeded the threshold value of ',D24.8) 2585 FORMAT(/,6X,'*** Math check: | A - U S U\' | / ( |A| n ulp )', + /,6X,'*** result was ',D24.8, + /,6X,'*** Math check: | I - UU\' | / ( n ulp )', + /,6X,'*** result was ',D24.8) 2595 FORMAT(/,6X,'*** Math check: | U\' A U - S | / ( |A| m ulp )', + /,6X,'*** result was ',D24.8, + /,6X,'*** Math check: | I - U\'U | / ( m ulp )', + /,6X,'*** result was ',D24.8) END ************************************************************************ * Subroutine Name: CM415 * *----------------------------------------------------------------------* * * * Description: Create test matrix using algorithm * * * *----------------------------------------------------------------------* * * * Parameter definitions: * * * * Input: * * LDA I = The leading dimension of matrix A * * N I = The order of matrix A * * * * Output: * * A R*8 = Real sysmetric matrix of order N * * * ************************************************************************ SUBROUTINE CM415 (A,LDA,N ) INTEGER I, J, LDA, N REAL*8 A(LDA,*) DO 100, I = 1, N DO 110, J = 1, N IF (I .EQ. N) THEN A(I,J) = J ELSEIF (J .EQ. N) THEN A(I,J) = I ELSEIF (I .EQ. J) THEN A(I,J) = 1 ELSE A(I,J) = 0 ENDIF 110 CONTINUE 100 CONTINUE RETURN END