c     program DRDSYMQL
c>> 1996-05-28 DRDSYMQL  Krogh Added external statement.
c>> 1994-10-19 DRDSYMQL  Krogh  Changes to use M77CON
c>> 1994-09-22 DRDSYMQL  CLL
c>> 1992-04-23 CLL
c>> 1992-03-04 DRDSYMQL  Krogh Initial version.
c     Demonstrate symmetric eigenvalue/eigenvector subroutine DSYMQL.
c     ------------------------------------------------------------------
c--D replaces "?": DR?SYMQL, ?SYMQL, ?VECP, ?MATP, ?DOT
c     ------------------------------------------------------------------
      integer I, IERR, J, LDA, N
      parameter (LDA = 4)
      double precision A(LDA,LDA), ASAV(LDA,LDA), ANORM, D(LDA,LDA)
      external DDOT
      double precision DDOT, EVAL(LDA), WORK(LDA)
      data A(1,1)          /  5.0d0 /
      data (A(2,J), J=1,2) /  4.0d0,  5.0d0 /
      data (A(3,J), J=1,3) /  1.0d0,  1.0d0, 4.0d0 /
      data (A(4,J), J=1,4) /  1.0d0,  1.0d0, 2.0d0, 4.0d0 /
      data ANORM / 11.0d0 /
      data N /LDA/
c     ------------------------------------------------------------------
      print*,'DRDSYMQL..  Demo driver for DSYMQL.'
c
c     First copy A() to ASAV() for later residual check.
c
      do 20 I = 1,N
         do 10 J = 1,I
            ASAV(I,J) = A(I,J)
            ASAV(J,I) = ASAV(I,J)
   10    continue
   20 continue
      call DSYMQL(A, LDA, N, EVAL, WORK, IERR)
      if (IERR .eq. 0) then
         call DVECP(EVAL, N, '0 Eigenvalues')
         call DMATP(A, LDA, N, N,  '0 Eigenvectors as column vectors')
c
c        As a check compute D = (ASAV*EVEC - EVEC*EVAL) / ANORM.
c        The EVEC's are in the array A().
c        Expect D to be close to machine precision.
c
         do 40 J = 1, N
            do 30 I = 1, N
               D(I, J) = (DDOT(N, ASAV(I,1), LDA, A(1,J), 1) -
     *                  A(I,J) * EVAL(J)) / ANORM
   30       continue
   40    continue
         call DMATP(D, LDA, N, N,
     *         '0 Residual matrix D = (A*EVEC - EVEC*EVAL) / ANORM')

      else
         print '(/a,i5)',' Convergence failure in DSYMQL, IERR =',IERR
      end if
      stop
      end