c     program DRDEVVUN
c>> 1996-05-28 DRDEVVUN  Krogh Added external statement.
c>> 1994-10-19 DRDEVVUN  Krogh  Changes to use M77CON
c>> 1994-09-22 DRDEVVUN  CLL
c>> 1992-04-23 CLL
c>> 1992-03-04 DRDEVVUN  Krogh Initial version.
c     Demonstrate unsymmetric eigenvalue/eigenvector subroutine DEVVUN.
c     ------------------------------------------------------------------
c--D replaces "?": DR?EVVUN, ?EVVUN, ?VECP, ?MATP, ?DOT
c     ------------------------------------------------------------------
      integer I, J, LDA, N
      parameter (LDA = 3)
      integer IFLAG(LDA)
      double precision A(LDA,LDA), VR(LDA), VI(LDA), VEC(LDA, LDA)
      double precision WORK(LDA), ASAV(LDA, LDA), D(LDA, LDA), ANORM
      external DDOT
      double precision DDOT
      data (A(1,I), I=1,3) /  8.0d0, -1.0d0, -5.0d0 /
      data (A(2,I), I=1,3) / -4.0d0,  4.0d0, -2.0d0 /
      data (A(3,I), I=1,3) / 18.0d0, -5.0d0, -7.0d0 /
      data ANORM / 30.0d0 /
      data N / LDA /
c     ------------------------------------------------------------------
      print*,'DRDEVVUN..  Demo driver for DEVVUN.'
c
c     First copy A() to ASAV() for later residual check.
c
         do 20 J = 1, N
            do 10 I = 1, N
               ASAV(I, J) = A(I, J)
   10       continue
   20    continue

         call DEVVUN(A(1, 1), LDA, N, VR, VI, VEC, IFLAG, WORK)

         print '(a, I2)',' IFLAG(1) =',IFLAG(1)
         if (IFLAG(1) .le. 2) then
            call DVECP(VR, N, ' Real part of the eigenvalues')
            call DVECP(VI, N, ' Imaginary part of the eigenvalues')
            call DMATP(VEC, LDA, N, N,
     *         '0 Eigenvectors as columns or pairs of columns')
c
c           As a check compute D = (ASAV*VEC - VEC*EVAL) / ANORM.
c           Expect D to be close to the machine precision.
c
            do 50 J = 1, N
               if (VI(J) .eq. 0.0d0) then
c                 Compute residual for a real eigenvalue and eigenvector
                  do 30 I = 1, N
                     D(I, J) = (DDOT(N, ASAV(I,1), LDA, VEC(1,J), 1) -
     *                  VEC(I,J) * VR(J)) / ANORM
   30             continue
               else if (VI(J) .gt. 0.0d0) then
                  do 40 I = 1, N
                     D(I, J) = (DDOT(N,ASAV(I,1),LDA,VEC(1,J),1) -
     *                  VEC(I,J)*VR(J)+VEC(I,J+1)*VI(J))/ANORM
                     D(I, J+1) = (DDOT(N,ASAV(I,1),LDA,VEC(1,J+1),1) -
     *                  VEC(I,J)*VI(J)-VEC(I,J+1)*VR(J))/ANORM
   40             continue
               end if
   50       continue
            call DMATP(D, LDA, N, N,
     *      '0 Packed residual matrix D = (A*EVEC - EVEC*EVAL) / ANORM')
         else
            print'(/a)', ' Failure in DEVVUN. '
         end if
      stop
      end