c
c     this driver tests  eispack  for the class of real generalized
c     matrix systems exhibiting the use of  eispack  to find all
c     the eigenvalues and eigenvectors or only all the eigenvalues for
c     the eigenproblem   a*x = (lambda)*b*x .
c
c     this driver is catalogued as  eispdrv4(rggeispk).
c
c     the dimension of  a ,  b  and  z  should be  nm  by  nm.
c     the dimension of  alfr,alfi,beta,  and  norm  should be  nm.
c     the dimension of  ahold  and  bhold  should be  nm  by  nm.
c     here nm = 20.
c
      double precision a( 20, 20),b( 20, 20),z( 20, 20),
     x        alfr( 20),alfi( 20),beta( 20),norm( 20),
     x        machep,tcrit,resdul,eps1
      double precision ahold( 20, 20),bhold( 20, 20)
      integer  error
      data iwrite/6/
c
c     .......... machep is a machine dependent parameter specifying
c                the relative precision of floating point arithmetic.
      machep = 1.0d0
    1 machep = 0.5d0*machep
      if (1.0d0 + machep .gt. 1.0d0 ) go to 1
      machep = 2.0d0*machep
c
c
      nm = 20
      open(unit=1,file='rggdata1')
      open(unit=2,file='rggdata2')
      rewind 1
      rewind 2
   10 call  rmatin(nm,n,a,b,ahold,bhold,0)
      write(iwrite,20) n
   20 format(30h1the full matrix  a   of order,
     x       i4,22h  is (printed by rows)  )
      do 30 i = 1,n
   30    write(iwrite,40) (a(i,j),j=1,n)
   40    format(/(1x,1p5d24.16))
      write(iwrite,41) n
   41 format(/////30h the full matrix  b   of order,
     x       i4,22h  is (printed by rows)  )
      do 42 i = 1,n
   42    write(iwrite,43) (b(i,j),j=1,n)
   43    format(/(1x,1p5d24.16))
      do  510  icall = 1,2
         if( icall  .ne. 1 )  call   rmatin(nm,n,a,b,ahold,bhold,1)
         go to  (70,75), icall
c
c     rggwz  using  qzit  and  qzvec
c     invoked from driver subroutine  rgg.
c
   70    write(iwrite,71)
   71    format(42h1all of the eigenvalues and corresponding ,
     x     40heigenvectors of the real general system /
     x     34h follow.  the path involving  qzit,
     x     12h  was used. )
         eps1 = 0.0d0
         call  rgg(nm,n,a,b,alfr,alfi,beta,1,z,error)
         write(iwrite,72) error
   72    format(//28h *****error from qzit***** = ,i4)
         go to  130
c
c     rggw  using  qzit
c     invoked from driver subroutine  rgg.
c
   75    write(iwrite,76)
   76    format(24h1all of the eigenvalues ,
     x     26hof the real general system /
     x     46h follow.  the path involving  qzit  was used.  )
         eps1 = 0.0d0
         call  rgg(nm,n,a,b,alfr,alfi,beta,0,a,error)
         write(iwrite,77) error
   77    format(//28h *****error from qzit***** = ,i4)
c
  130    write(iwrite,291)
  291    format(//50h each eigenvalue of  a*x-l*b*x  can be found from ,
     x   59h alfr,  alfi,  and  beta  by dividing  beta  into  alfr  to/
     x   59h produce the real part and dividing  beta  into  alfi  to p,
     x   26hroduce the imaginary part.)
         write(iwrite,292)
  292    format(///14x,7halfr(i),20x,7halfi(i),20x,7hbeta(i))
         do 295 i = 1,n
            write(iwrite,293) i,alfr(i),alfi(i),beta(i)
  293       format(i4,3(1pd24.16,3x))
  295    continue
         if( icall .gt. 1 ) go to 510
         call   rmatin(nm,n,a,b,ahold,bhold,1)
         call   rggwzr(nm,n,a,b,alfr,alfi,beta,z,norm,resdul)
         tcrit = resdul/(dfloat(10*n)*machep)
         write(iwrite,300) tcrit,resdul
  300    format(///48h as a guide to evaluating the performance of the,
     x     55h  eispack  codes, a number  x, related to the  1-norm  /
     x     56h of the residual matrix, may be normalized to the number,
     x     56h  y  by dividing by the product of the machine precision/
     x     57h machep  and  10*n  where  n  is the order of the matrix.,
     x     36h  the  eispack  codes have performed /
     x     56h (well,satisfactorily,poorly) according as the ratio  y ,
     x     54h is (less than  1, less than  100, greater than  100)./
     x     23h for this run,  y  is =,1pd16.8,11h  with  x =,d16.8)
         nm2 = n-2
         k = 1
  310    if( k .gt. n )  go to  510
         if( alfi(k) .eq. 0.0d0 )  go to  360
         if( k .gt. nm2 )  go to  340
         if( alfi(k+2) .eq. 0.0d0 )  go to  330
c
c        both eigenvectors are complex.
c
         write(iwrite,320) alfr(k),alfi(k),beta(k),alfr(k+2),alfi(k+2),
     x     beta(k+2),norm(k),norm(k+2),(z(i,k),z(i,k+1),z(i,k+2),
     x     z(i,k+3),i=1,n)
  320    format(///1x,2(7x,4halfr,14x,4halfi,14x,4hbeta,7x)
     x       /1x,6(1pd16.8,2x)
     x       //1x,2(8x,39h1-norm of corresponding residual vector,7x)//
     x       1x,2(7x,21h!!beta*a*x-alfa*b*x!!,26x)/
     x       2x,2(33h---------------------------------,1x,1h=,d16.8,3x)/
     x       1x,2(2x,31h!!x!!*(beta*!!a!!+!alfa!*!!b!!),21x)//
     x       1x,2(14x,25hcorresponding eigenvector,15x)/
     x       (2d24.16,6x,2d24.16,6x))
         k = k+4
         go to  310
c
c        one complex eigenvector and one real eigenvector in that
c        order.
c
  330    write(iwrite,320)alfr(k+2),alfi(k+2),beta(k+2),alfr(k),alfi(k),
     x       beta(k),norm(k+2),norm(k),(z(i,k+2),alfi(k+2),z(i,k),
     x       z(i,k+1),i=1,n )
         k = k+3
         go to  310
c
c        one complex eigenvector.
c
  340    write(iwrite,350) alfr(k),alfi(k),beta(k),norm(k),(z(i,k),
     x               z(i,k+1),i=1,n)
  350    format(///8x,4halfr,14x,4halfi,14x,4hbeta,
     x       23x/1x,3(1pd16.8,2x)//
     x       1x,8x,39h1-norm of corresponding residual vector,7x//
     x       1x,7x,21h!!beta*a*x-alfa*b*x!!/
     x       2x,33(1h-),1x,1h=,d16.8,8x/
     x       1x,2x,31h!!x!!*(beta*!!a!!+!alfa!*!!b!!)//
     x       1x,14x,25hcorresponding eigenvector,15x/
     x       (2d24.16,6x))
         go to  510
  360    if( k .eq. n )  go to  380
         if( alfi(k+1) .eq. 0.0d0 )  go to  370
c
c        one real eigenvector and one complex eigenvector in that
c        order.
c
         write(iwrite,320) alfr(k),alfi(k),beta(k),alfr(k+1),alfi(k+1),
     x                beta(k+1),norm(k),norm(k+1),
     x                (z(i,k),alfi(k),z(i,k+1),z(i,k+2),i=1,n)
         k = k+3
         go to  310
c
c        both eigenvectors are real.
c
  370    write(iwrite,320) alfr(k),alfi(k),beta(k),alfr(k+1),alfi(k+1),
     x          beta(k+1),norm(k),norm(k+1),
     x                (z(i,k),alfi(k),z(i,k+1),alfi(k+1),i=1,n)
         k = k+2
         go to  310
c
c        one real eigenvector.
c
  380    write(iwrite,350) alfr(k),alfi(k),beta(k),norm(k),(z(i,k),
     x          alfi(k),i=1,n)
         go to  510
  510 continue
      go to  10
      end