home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / eispack-1.0-src.tgz / tar.out / contrib / eispack / tests / rspeispk.f < prev    next >
Text File  |  1996-09-28  |  12KB  |  321 lines
  1. c
  2. c     this driver tests  eispack  for the class of real symmetric
  3. c     packed matrices exhibiting the use of  eispack  to find all the
  4. c     eigenvalues and eigenvectors, only all the eigenvalues, some of
  5. c     the eigenvalues and the corresponding eigenvectors, or only some
  6. c     of the eigenvalues.
  7. c
  8. c     this driver is catalogued as  eispdrv4(rspeispk).
  9. c
  10. c     the dimension of  a  should be  nnn  and the dimension
  11. c     of  z  should be  nm  by  nm.
  12. c     the dimension of  w,norm,d,e,e2,ind,rv1,rv2,rv3,rv4,rv5,  and
  13. c     rv6  should be  nm.
  14. c     the dimension of  ahold  should be  nnn.
  15. c     here nm = 20 and  nnn = nm*(nm+1)/2.
  16. c
  17.       double precision a(210),z(20,20),ahold(210),w(20),d(20),e(20),
  18.      x        e2( 20),rv1( 20),rv2( 20),rv3( 20),rv4( 20),
  19.      x        rv5( 20),rv6( 20),norm( 20),tcrit,machep,resdul,
  20.      x        maxeig,maxdif,eigdif,u,lb,ub,eps1
  21.       integer  ind( 20),error
  22.       data ireadc/5/,iwrite/6/
  23. c
  24. c     .......... machep is a machine dependent parameter specifying
  25. c                the relative precision of floating point arithmetic.
  26. c
  27.       machep = 1.0d0
  28.     1 machep = 0.5d0*machep
  29.       if (1.0d0 + machep .gt. 1.0d0 ) go to 1
  30.       machep = 2.0d0*machep
  31. c
  32.       nm = 20
  33.       nnn = (nm*(nm+1))/2
  34.    10 call rmatin(nnn,n,a,ahold,0)
  35.       write(iwrite,20) n
  36.    20 format(50h1the symmetric matrix  a   in packed form of order,
  37.      x       i4,22h  is (printed by rows)  )
  38.       k=1
  39.       do 35 l=1,n
  40.          k=k+l-1
  41.          m=k+l-1
  42.    35    write(iwrite,40) (a(i),i=k,m) 
  43.    40    format(/(1x,1p5d24.16))
  44.       read(ireadc,50) mm,lb,ub,m11,no 
  45.    50 format(i4,2d24.16,2(4x,i4))
  46. c
  47. c     mm,lb,ub,m11,  and  no  are read from sysin after the matrix is
  48. c     generated.  mm,lb,  and  ub  specify to  bisect  the maximum
  49. c     number of eigenvalues and the bounds for the interval which is
  50. c     to be searched.  m11  and  no  specify to  tridib  the lower
  51. c     boundary index and the number of desired eigenvalues.
  52. c
  53.       do  230  icall = 1,10
  54.          if( icall .ne. 1 )  call  rmatin(nnn,n,a,ahold,1)
  55. c
  56. c     if  tqlrat  path (label 80) is taken then  tql2  path (label 70)
  57. c     must also be taken in order that the measure of performance be
  58. c     meaningful.
  59. c     if  imtql1  path (label 85) is taken then  imtql2  path (label 75)
  60. c     must also be taken in order that the measure of performance be
  61. c     meaningful.
  62. c     if  tql2  (imtql2)  path fails, then  tqlrat  (imtql1)  path is
  63. c     omitted and printout flagged with  -1.0.
  64. c
  65.          go to  (70,75,80,85,89,90,95,100,110,115),  icall
  66. c
  67. c     rspwz  using  tql2
  68. c     invoked from driver subroutine  rsp.
  69. c
  70.    70    write(iwrite,71)
  71.    71    format(42h1all of the eigenvalues and corresponding ,
  72.      x     49heigenvectors of the real symmetric matrix follow./
  73.      x     36h the path involving  tql2  was used.  )
  74.          call  rsp(nm,n,nnn,a,w,1,z,e,e,error)
  75.          write(iwrite,72) error
  76.    72    format(//29h *****error from tql2***** = ,i4)
  77.          do 73 i = 1,n
  78.             rv6(i) = w(i)
  79.    73    continue
  80.          m = n
  81.          if( error.ne.0) m = error - 1
  82.          go to  130
  83. c
  84. c     rspwz  using  imtql2
  85. c
  86.    75    write(iwrite,76)
  87.    76    format(42h1all of the eigenvalues and corresponding ,
  88.      x     49heigenvectors of the real symmetric matrix follow./
  89.      x     38h the path involving  imtql2  was used.  )
  90.          do 767 i=1,n
  91.            do 765 j=1,n
  92.   765        z(i,j)=0.0d0
  93.   767      z(i,i)=1.0d0
  94.          call  tred3(n,nnn,a,w,e,e)
  95.          call  imtql2(nm,n,w,e,z,error)
  96.          write(iwrite,77) error
  97.    77    format(//31h *****error from imtql2***** = ,i4)
  98.          do 78 i = 1,n
  99.             rv5(i) = w(i)
  100.    78    continue
  101.          m = n
  102.          if( error .ne. 0 ) m = error - 1
  103.          call  trbak3(nm,n,nnn,a,m,z)
  104.          go to  130
  105. c
  106. c     rspw  using  tqlrat
  107. c     invoked from driver subroutine  rsp.
  108. c
  109.    80    write(iwrite,805)
  110.   805    format(24h1all of the eigenvalues ,
  111.      x     36hof the real symmetric matrix follow./
  112.      x     38h the path involving  tqlrat  was used.  )
  113.          call  rsp(nm,n,nnn,a,w,0,z,e,e2,error)
  114.          write(iwrite,807) error
  115.   807    format(//31h *****error from tqlrat***** = ,i4)
  116.          maxeig = 0.0d0
  117.          maxdif = 0.0d0
  118.          m = n
  119.          if( error .ne. 0 ) m = error - 1
  120.          if( m .eq. 0 ) go to 230
  121.          do 81 i = 1,n
  122.             if( dabs(w(i)) .gt. maxeig ) maxeig = dabs(w(i))
  123.             u = dabs(rv6(i) - w(i))
  124.             if( u .gt. maxdif ) maxdif = u
  125.    81    continue
  126.          if( maxeig .eq. 0.0d0 ) maxeig = 1.0d0
  127.          eigdif = maxdif/(maxeig*machep*dfloat(10*n))
  128.          write(iwrite,82)eigdif
  129.    82    format(//49h comparison of the eigenvalues from  tqlrat  with,
  130.      x          52h those from  tql2  gives the normalized difference  ,
  131.      x          1pd16.8)
  132.          go to 130
  133. c
  134. c     rspw  using  imtql1
  135. c
  136.    85    write(iwrite,855)
  137.   855    format(24h1all of the eigenvalues ,
  138.      x     36hof the real symmetric matrix follow./
  139.      x     38h the path involving  imtql1  was used.  )
  140.          call  tred3(n,nnn,a,w,e,e)
  141.          call  imtql1(n,w,e,error)
  142.          write(iwrite,857) error
  143.   857    format(//31h *****error from imtql1***** = ,i4)
  144.          maxeig = 0.0d0
  145.          maxdif = 0.0d0
  146.          m = n
  147.          if( error .ne. 0 ) m = error - 1
  148.          if( m .eq. 0 ) go to 230
  149.          do 86 i = 1,m
  150.             if( dabs(w(i)) .gt. maxeig ) maxeig = dabs(w(i))
  151.             u = dabs(rv5(i) - w(i))
  152.             if( u .gt. maxdif ) maxdif = u
  153.    86    continue
  154.          if( maxeig .eq. 0.0d0 ) maxeig = 1.0d0
  155.          eigdif = maxdif/(maxeig*machep*dfloat(10*n))
  156.          write(iwrite,87) eigdif
  157.    87    format(//49h comparison of the eigenvalues from  imtql1  with,
  158.      x          52h those from  imtql2  gives the normalized difference,
  159.      x          1pd16.8)
  160.          go to  130
  161. c
  162. c     rspw1z  ( usage here computes all the eigenvectors )
  163. c
  164.   890    format(43h1some of the eigenvalues and corresponding ,
  165.      x     49heigenvectors of the real symmetric matrix follow./
  166.      x     38h the path involving  imtqlv  was used.  )
  167.    89    write(iwrite,890)
  168.          call  tred3(n,nnn,a,d,e,e2)
  169.          call imtqlv(n,d,e,e2,w,ind,error,rv1)
  170.          write(iwrite,895) error
  171.   895    format(//31h *****error from imtqlv***** = ,i4)
  172.          m = n
  173.          if( error .ne. 0 ) m = error - 1
  174.          call tinvit(nm,n,d,e,e2,m,w,ind,z,error,rv1,rv2,rv3,rv4,rv6)
  175.          write(iwrite,98) error
  176.          call  trbak3(nm,n,nnn,a,m,z)
  177.          go to 130
  178. c
  179. c     rsp1w1z  using  tsturm
  180. c
  181.    90    write(iwrite,91)
  182.    91    format(43h1some of the eigenvalues and corresponding ,
  183.      x     49heigenvectors of the real symmetric matrix follow./
  184.      x     38h the path involving  tsturm  was used.  )
  185.          eps1 = 0.0d0
  186.          call  tred3(n,nnn,a,d,e,e2)
  187.          call  tsturm(nm,n,eps1,d,e,e2,lb,ub,mm,m,w,z,error,rv1,rv2,
  188.      x                rv3,rv4,rv5,rv6)
  189.          write(iwrite,92) error
  190.    92    format(//31h *****error from tsturm***** = ,i4)
  191.          if( error .eq. 3*n + 1 ) go to 230
  192.          if( error .ne. 0 ) m = error - 4*n + 1
  193.          call  trbak3(nm,n,nnn,a,m,z)
  194.          go to  150
  195. c
  196. c     rsp1w1z  using  bisect  and  tinvit
  197. c
  198.    95    write(iwrite,96)
  199.    96    format(43h1some of the eigenvalues and corresponding ,
  200.      x     49heigenvectors of the real symmetric matrix follow./
  201.      x     38h the path involving  bisect  was used.  )
  202.          eps1 = 0.0d0
  203.          call  tred3(n,nnn,a,d,e,e2)
  204.          call  bisect(n,eps1,d,e,e2,lb,ub,mm,m,w,ind,error,rv4,rv5)
  205.          write(iwrite,97) error
  206.    97    format(//31h *****error from bisect***** = ,i4)
  207.          if( error .ne. 0 ) go to 230
  208.          call  tinvit(nm,n,d,e,e2,m,w,ind,z,error,rv1,rv2,rv3,rv4,rv6)
  209.          write(iwrite,98) error
  210.    98    format(//31h *****error from tinvit***** = ,i4)
  211.          call  trbak3(nm,n,nnn,a,m,z)
  212.          go to  150
  213. c
  214. c     rsp1w  using  bisect
  215. c
  216.   100    write(iwrite,101)
  217.   101    format(25h1some of the eigenvalues ,
  218.      x     36hof the real symmetric matrix follow. )
  219.          eps1 = 0.0d0
  220.          call  tred3(n,nnn,a,d,e,e2)
  221.          call  bisect(n,eps1,d,e,e2,lb,ub,mm,m,w,ind,error,rv4,rv5)
  222.          write(iwrite,102) error
  223.   102    format(//31h *****error from bisect***** = ,i4)
  224.          go to  150
  225. c
  226. c     rsp1w1z  using  tridib  and  tinvit
  227. c
  228.   110    write(iwrite,111)
  229.   111    format(43h1some of the eigenvalues and corresponding ,
  230.      x     49heigenvectors of the real symmetrix matrix follow./
  231.      x     38h the path involving  tridib  was used.  )
  232.          eps1 = 0.0d0
  233.          call  tred3(n,nnn,a,d,e,e2)
  234.          call  tridib(n,eps1,d,e,e2,lb,ub,m11,no,w,ind,error,rv4,rv5)
  235.          write(iwrite,112) error
  236.   112    format(//31h *****error from tridib***** = ,i4)
  237.          if( error .ne. 0 )  go to  230
  238.          m = no
  239.          write(iwrite,113) lb,ub
  240.   113    format(//49h the eigenvalues as determined by  tridib  are in
  241.      x          , 13h the interval,1pd16.8,5h  to ,d16.8)
  242.          call  tinvit(nm,n,d,e,e2,m,w,ind,z,error,rv1,rv2,rv3,rv4,rv6)
  243.          write(iwrite,114) error 
  244.   114    format(//31h *****error from  tinvit***** = ,i4)
  245.          call  trbak3(nm,n,nnn,a,m,z)
  246.          go to 150
  247. c
  248. c     rsp1w  using  tridib
  249. c
  250.   115    write(iwrite,101)
  251.          eps1 = 0.0d0
  252.          call  tred3(n,nnn,a,d,e,e2)
  253.          call  tridib(n,eps1,d,e,e2,lb,ub,m11,no,w,ind,error,rv4,rv5)
  254.          write(iwrite,117) error
  255.   117    format(//31h *****error from tridib***** = ,i4)
  256.          m = no
  257.          write(iwrite,113) lb,ub
  258.          if( error .ne. 0 )  go to  230
  259.          go to 150
  260. c
  261.   130    write(iwrite,140) (i,w(i),i=1,n)
  262.   140    format(///18h eigenvalues of  a/3(i3,2x,1pd24.16,3x))
  263.          if( icall .eq. 3 .or. icall .eq. 4 )  go to  230
  264.          go to  190
  265.   150    write(iwrite,160) m,lb,ub 
  266.   160    format(///4h the ,i4,29h  eigenvalues of  a  between ,1pd24.16,
  267.      x        6h  and ,d24.16)
  268.          if( m .eq. 0 )  go to  230
  269.          write(iwrite,163) (i,w(i),i=1,m)  
  270.   163    format(//3(i3,2x,1pd24.16,3x))
  271.          if( icall .eq. 8 .or. icall .eq. 10 )  go to  230
  272.   190    call  rmatin(nnn,n,a,ahold,1)
  273.          call  rspwzr(nm,n,nnn,m,a,w,z,norm,resdul)
  274.          tcrit = resdul/(dfloat(10*n)*machep)
  275.          write(iwrite,200) tcrit,resdul 
  276.   200    format(///48h as a guide to evaluating the performance of the,
  277.      x     53h  eispack  codes, a number  x, related to the  1-norm /
  278.      x     56h of the residual matrix, may be normalized to the number,
  279.      x     56h  y  by dividing by the product of the machine precision/
  280.      x     57h machep  and  10*n  where  n  is the order of the matrix.,
  281.      x     36h  the  eispack  codes have performed /
  282.      x     56h (well,satisfactorily,poorly) according as the ratio  y ,
  283.      x     54h is (less than  1, less than  100, greater than  100)./
  284.      x     23h for this run,  y  is =,1pd14.6,11h  with  x =,d14.6)
  285.          do  220  kk = 1,m,3
  286.             kk1 = min0(kk+2,m)
  287.             irs = kk1-kk+1
  288.             go to  (204,208,212),  irs
  289.   204       write(iwrite,205) (w(k),k=kk,kk1),(norm(k),k=kk,kk1), 
  290.      x           ((z(i,k),k=kk,kk1),i=1,n)
  291.   205    format(///1x,1(9x,15heigenvalue of a,8x)/1x,1(4x,1pd24.16,4x)//
  292.      x   1x,1(5x,23h1-norm of corresponding,4x)/
  293.      x   1x,1(9x,15hresidual vector,8x)/
  294.      x   1x,1(4x,24h!!ax-l*x!!/(!!x!!*!!a!!),4x)/
  295.      x   1x,1(8x,d16.8,8x)//
  296.      x   1x,1(4x,26hcorresponding eigenvector ,2x)/1(5x,d24.16,3x))
  297.             go to 220
  298.   208       write(iwrite,209) (w(k),k=kk,kk1),(norm(k),k=kk,kk1),
  299.      x           ((z(i,k),k=kk,kk1),i=1,n)
  300.   209    format(///1x,2(9x,15heigenvalue of a,8x)/1x,2(4x,1pd24.16,4x)//
  301.      x   1x,2(5x,23h1-norm of corresponding,4x)/
  302.      x   1x,2(9x,15hresidual vector,8x)/
  303.      x   1x,2(4x,24h!!ax-l*x!!/(!!x!!*!!a!!),4x)/
  304.      x   1x,2(8x,d16.8,8x)//
  305.      x   1x,2(4x,26hcorresponding eigenvector ,2x)/(5x,d24.16,3x,
  306.      x   5x,d24.16,3x))
  307.             go to 220
  308.   212       write(iwrite,213) (w(k),k=kk,kk1),(norm(k),k=kk,kk1),
  309.      x           ((z(i,k),k=kk,kk1),i=1,n)
  310.   213    format(///1x,3(9x,15heigenvalue of a,8x)/1x,3(4x,1pd24.16,4x)//
  311.      x   1x,3(5x,23h1-norm of corresponding,4x)/
  312.      x   1x,3(9x,15hresidual vector,8x)/
  313.      x   1x,3(4x,24h!!ax-l*x!!/(!!x!!*!!a!!),4x)/
  314.      x   1x,3(8x,d16.8,8x)//
  315.      x   1x,3(4x,26hcorresponding eigenvector ,2x)/(5x,d24.16,3x,
  316.      x   5x,d24.16,3x,5x,d24.16,3x))
  317.   220    continue
  318.   230 continue
  319.       go to  10
  320.       end
  321.