home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / eispack-1.0-src.tgz / tar.out / contrib / eispack / ex / rggeispk.f < prev    next >
Text File  |  1996-09-28  |  7KB  |  175 lines
  1. c
  2. c     this driver tests  eispack  for the class of real generalized
  3. c     matrix systems exhibiting the use of  eispack  to find all
  4. c     the eigenvalues and eigenvectors or only all the eigenvalues for
  5. c     the eigenproblem   a*x = (lambda)*b*x .
  6. c
  7. c     this driver is catalogued as  eispdrv4(rggeispk).
  8. c
  9. c     the dimension of  a ,  b  and  z  should be  nm  by  nm.
  10. c     the dimension of  alfr,alfi,beta,  and  norm  should be  nm.
  11. c     the dimension of  ahold  and  bhold  should be  nm  by  nm.
  12. c     here nm = 20.
  13. c
  14.       double precision a( 20, 20),b( 20, 20),z( 20, 20),
  15.      x        alfr( 20),alfi( 20),beta( 20),norm( 20),
  16.      x        machep,tcrit,resdul,eps1
  17.       double precision ahold( 20, 20),bhold( 20, 20)
  18.       integer  error
  19.       data iwrite/6/
  20. c
  21. c     .......... machep is a machine dependent parameter specifying
  22. c                the relative precision of floating point arithmetic.
  23.       machep = 1.0d0
  24.     1 machep = 0.5d0*machep
  25.       if (1.0d0 + machep .gt. 1.0d0 ) go to 1
  26.       machep = 2.0d0*machep
  27. c
  28. c
  29.       nm = 20
  30.       open(unit=1,file='rggdata1')
  31.       open(unit=2,file='rggdata2')
  32.       rewind 1
  33.       rewind 2
  34.    10 call  rmatin(nm,n,a,b,ahold,bhold,0)
  35.       write(iwrite,20) n
  36.    20 format(30h1the full matrix  a   of order,
  37.      x       i4,22h  is (printed by rows)  )
  38.       do 30 i = 1,n
  39.    30    write(iwrite,40) (a(i,j),j=1,n)
  40.    40    format(/(1x,1p5d24.16))
  41.       write(iwrite,41) n 
  42.    41 format(/////30h the full matrix  b   of order,
  43.      x       i4,22h  is (printed by rows)  )
  44.       do 42 i = 1,n
  45.    42    write(iwrite,43) (b(i,j),j=1,n)
  46.    43    format(/(1x,1p5d24.16))
  47.       do  510  icall = 1,2
  48.          if( icall  .ne. 1 )  call   rmatin(nm,n,a,b,ahold,bhold,1)
  49.          go to  (70,75), icall
  50. c
  51. c     rggwz  using  qzit  and  qzvec
  52. c     invoked from driver subroutine  rgg.
  53. c
  54.    70    write(iwrite,71)
  55.    71    format(42h1all of the eigenvalues and corresponding ,
  56.      x     40heigenvectors of the real general system /
  57.      x     34h follow.  the path involving  qzit,
  58.      x     12h  was used. )
  59.          eps1 = 0.0d0
  60.          call  rgg(nm,n,a,b,alfr,alfi,beta,1,z,error)
  61.          write(iwrite,72) error
  62.    72    format(//28h *****error from qzit***** = ,i4)
  63.          go to  130
  64. c
  65. c     rggw  using  qzit
  66. c     invoked from driver subroutine  rgg.
  67. c
  68.    75    write(iwrite,76)
  69.    76    format(24h1all of the eigenvalues ,
  70.      x     26hof the real general system /
  71.      x     46h follow.  the path involving  qzit  was used.  )
  72.          eps1 = 0.0d0
  73.          call  rgg(nm,n,a,b,alfr,alfi,beta,0,a,error)
  74.          write(iwrite,77) error
  75.    77    format(//28h *****error from qzit***** = ,i4)
  76. c
  77.   130    write(iwrite,291)
  78.   291    format(//50h each eigenvalue of  a*x-l*b*x  can be found from ,
  79.      x   59h alfr,  alfi,  and  beta  by dividing  beta  into  alfr  to/
  80.      x   59h produce the real part and dividing  beta  into  alfi  to p,
  81.      x   26hroduce the imaginary part.)
  82.          write(iwrite,292)
  83.   292    format(///14x,7halfr(i),20x,7halfi(i),20x,7hbeta(i))
  84.          do 295 i = 1,n
  85.             write(iwrite,293) i,alfr(i),alfi(i),beta(i)
  86.   293       format(i4,3(1pd24.16,3x))
  87.   295    continue
  88.          if( icall .gt. 1 ) go to 510
  89.          call   rmatin(nm,n,a,b,ahold,bhold,1)
  90.          call   rggwzr(nm,n,a,b,alfr,alfi,beta,z,norm,resdul)
  91.          tcrit = resdul/(dfloat(10*n)*machep)
  92.          write(iwrite,300) tcrit,resdul
  93.   300    format(///48h as a guide to evaluating the performance of the,
  94.      x     55h  eispack  codes, a number  x, related to the  1-norm  /
  95.      x     56h of the residual matrix, may be normalized to the number,
  96.      x     56h  y  by dividing by the product of the machine precision/
  97.      x     57h machep  and  10*n  where  n  is the order of the matrix.,
  98.      x     36h  the  eispack  codes have performed /
  99.      x     56h (well,satisfactorily,poorly) according as the ratio  y ,
  100.      x     54h is (less than  1, less than  100, greater than  100)./
  101.      x     23h for this run,  y  is =,1pd16.8,11h  with  x =,d16.8)
  102.          nm2 = n-2
  103.          k = 1
  104.   310    if( k .gt. n )  go to  510
  105.          if( alfi(k) .eq. 0.0d0 )  go to  360
  106.          if( k .gt. nm2 )  go to  340
  107.          if( alfi(k+2) .eq. 0.0d0 )  go to  330
  108. c
  109. c        both eigenvectors are complex.
  110. c
  111.          write(iwrite,320) alfr(k),alfi(k),beta(k),alfr(k+2),alfi(k+2),
  112.      x     beta(k+2),norm(k),norm(k+2),(z(i,k),z(i,k+1),z(i,k+2),
  113.      x     z(i,k+3),i=1,n)
  114.   320    format(///1x,2(7x,4halfr,14x,4halfi,14x,4hbeta,7x)
  115.      x       /1x,6(1pd16.8,2x)
  116.      x       //1x,2(8x,39h1-norm of corresponding residual vector,7x)//
  117.      x       1x,2(7x,21h!!beta*a*x-alfa*b*x!!,26x)/
  118.      x       2x,2(33h---------------------------------,1x,1h=,d16.8,3x)/
  119.      x       1x,2(2x,31h!!x!!*(beta*!!a!!+!alfa!*!!b!!),21x)//
  120.      x       1x,2(14x,25hcorresponding eigenvector,15x)/
  121.      x       (2d24.16,6x,2d24.16,6x))
  122.          k = k+4
  123.          go to  310
  124. c
  125. c        one complex eigenvector and one real eigenvector in that
  126. c        order.
  127. c
  128.   330    write(iwrite,320)alfr(k+2),alfi(k+2),beta(k+2),alfr(k),alfi(k),
  129.      x       beta(k),norm(k+2),norm(k),(z(i,k+2),alfi(k+2),z(i,k),
  130.      x       z(i,k+1),i=1,n )
  131.          k = k+3
  132.          go to  310
  133. c
  134. c        one complex eigenvector.
  135. c
  136.   340    write(iwrite,350) alfr(k),alfi(k),beta(k),norm(k),(z(i,k),
  137.      x               z(i,k+1),i=1,n)
  138.   350    format(///8x,4halfr,14x,4halfi,14x,4hbeta,
  139.      x       23x/1x,3(1pd16.8,2x)//
  140.      x       1x,8x,39h1-norm of corresponding residual vector,7x//
  141.      x       1x,7x,21h!!beta*a*x-alfa*b*x!!/
  142.      x       2x,33(1h-),1x,1h=,d16.8,8x/
  143.      x       1x,2x,31h!!x!!*(beta*!!a!!+!alfa!*!!b!!)//
  144.      x       1x,14x,25hcorresponding eigenvector,15x/
  145.      x       (2d24.16,6x))
  146.          go to  510
  147.   360    if( k .eq. n )  go to  380
  148.          if( alfi(k+1) .eq. 0.0d0 )  go to  370
  149. c
  150. c        one real eigenvector and one complex eigenvector in that
  151. c        order.
  152. c
  153.          write(iwrite,320) alfr(k),alfi(k),beta(k),alfr(k+1),alfi(k+1),
  154.      x                beta(k+1),norm(k),norm(k+1),
  155.      x                (z(i,k),alfi(k),z(i,k+1),z(i,k+2),i=1,n)
  156.          k = k+3
  157.          go to  310
  158. c
  159. c        both eigenvectors are real.
  160. c
  161.   370    write(iwrite,320) alfr(k),alfi(k),beta(k),alfr(k+1),alfi(k+1),
  162.      x          beta(k+1),norm(k),norm(k+1),
  163.      x                (z(i,k),alfi(k),z(i,k+1),alfi(k+1),i=1,n)
  164.          k = k+2
  165.          go to  310
  166. c
  167. c        one real eigenvector.
  168. c
  169.   380    write(iwrite,350) alfr(k),alfi(k),beta(k),norm(k),(z(i,k),
  170.      x          alfi(k),i=1,n)
  171.          go to  510
  172.   510 continue
  173.       go to  10
  174.       end
  175.