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 / rsteispk.f < prev    next >
Text File  |  1996-09-28  |  13KB  |  334 lines
  1. c
  2. c     this driver tests  eispack  for the class of real symmetric tri-
  3. c     diagonal 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(rsteispk).
  9. c
  10. c     the dimension of  st  should be  nm  by  2.
  11. c     the dimension of  z  should be  nm  by  nm.
  12. c     the dimension of  w,norm,e,e2,ind,rv1,rv2,rv3,rv4,rv5,  and  rv6
  13. c     should be  nm.
  14. c     the dimension of  sthold  should be  nm  by  2.
  15. c     here nm = 20.
  16. c
  17.       double precision z(20,20),st(20,2),sthold(20,2),w(20),e(20),
  18.      x        e2(20),rv1( 20),rv2( 20),rv3( 20),rv4( 20),rv5( 20),  
  19.      x        rv6( 20),norm( 20),tcrit,machep,resdul,maxeig,
  20.      x        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.       mmb = 2
  34.    10 call rmatin(nm,mmb,n,mb,st,sthold,0)
  35.       write(iwrite,20) n
  36.    20 format(33h1the subdiagonal and diagonal in ,
  37.      x   49hthat order of the tridiagonal symmetric matrix  a,
  38.      x   10h  of order,i4,4h  is)
  39.       do  30  jj = 1,2
  40.          j = 3-jj
  41.    30    write(iwrite,40) (st(i,jj),i=j,n)
  42.    40    format(/(1x,1p5d24.16))
  43.       read(ireadc,50) mm,lb,ub,m11,no 
  44.    50 format(i4,2d24.16,2(4x,i4))
  45. c
  46. c     mm,lb,ub,m11,  and  no  are read from sysin after the matrix is
  47. c     generated.  mm,lb,  and  ub  specify to  bisect  the maximum
  48. c     number of eigenvalues and bounds for the interval which is to
  49. c     be searched. m11  and  no  specify to  tridib  the lower boundary
  50. c     index and the number of desired eigenvalues.
  51. c
  52.       do  230  icall = 1,10
  53.          if( icall .ne. 1 )  call  rmatin(nm,mmb,n,mb,st,sthold,1)
  54. c
  55. c
  56. c     if  tql1  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
  63.          go to  (70,75,80,85,89,90,95,100,110,115),  icall
  64. c
  65. c     rstwz  using  tql2
  66. c
  67.    70    write(iwrite,705)
  68.   705    format(42h1all of the eigenvalues and corresponding ,
  69.      x     53heigenvectors of the real symmetric tridiagonal matrix /
  70.      x     45h follow.  the path involving  tql2  was used.  )
  71.          do  72  i = 1,n
  72.             do  71  j = 1,n
  73.    71         z(i,j) = 0.0d0
  74.             z(i,i) = 1.0d0
  75.             e(i) = st(i,1)
  76.    72       w(i) = st(i,2)
  77.          call  tql2(nm,n,w,e,z,error)
  78.          write(iwrite,725) error
  79.   725    format(//29h *****error from tql2***** = ,i4)
  80.          do 73 i = 1,n
  81.             rv6(i) = w(i)
  82.    73    continue
  83.          m = n
  84.          if( error .ne. 0 ) m = error - 1
  85.          go to  130
  86. c
  87. c     rstwz  using  imtql2
  88. c     invoked from driver subroutine  rst.
  89. c
  90.    75    write(iwrite,755)
  91.   755    format(42h1all of the eigenvalues and corresponding ,
  92.      x     53heigenvectors of the real symmetric tridiagonal matrix /
  93.      x     47h follow.  the path involving  imtql2  was used.  )
  94.          do  77  i = 1,n
  95.             e(i) = st(i,1)
  96.    77       w(i) = st(i,2)
  97.          call  rst(nm,n,w,e,1,z,error)
  98.          write(iwrite,775) error
  99.   775    format(//31h *****error from imtql2***** = ,i4)
  100.          do 78 i = 1,n
  101.             rv5(i) = w(i)
  102.    78    continue
  103.          m = n
  104.          if( error .ne. 0 ) m = error - 1
  105.          go to  130
  106. c
  107. c     rstw  using  tql1
  108. c
  109.    80    write(iwrite,805)
  110.   805    format(24h1all of the eigenvalues ,
  111.      x     48hof the real symmetric tridiagonal matrix follow. /
  112.      x     36h the path involving  tql1  was used.  )
  113.          do  81  i = 1,n
  114.              e(i) = st(i,1)
  115.    81        w(i) = st(i,2)
  116.          call  tql1(n,w,e,error)
  117.          write(iwrite,815) error
  118.   815    format(//29h *****error from tql1***** = ,i4)
  119.          maxeig = 0.0d0
  120.          maxdif = 0.0d0
  121.          m = n
  122.          if( error .ne. 0 ) m = error - 1
  123.          if( m .eq. 0 ) go to 230
  124.          do 82 i = 1,m
  125.             if( dabs(w(i)) .gt. maxeig ) maxeig = dabs(w(i))
  126.             u = dabs(rv6(i) - w(i))
  127.             if( u .gt. maxdif ) maxdif = u
  128.    82    continue
  129.          if( maxeig .eq. 0.0d0 ) maxeig = 1.0d0
  130.          eigdif = maxdif/(maxeig*machep*dfloat(10*n))
  131.          write(iwrite,83) eigdif
  132.    83    format(//47h comparison of the eigenvalues from  tql1  with,
  133.      x          52h those from  tql2  gives the normalized difference  ,
  134.      x          1pd16.8)
  135.          go to 130
  136. c
  137. c     rstw  using  imtql1
  138. c     invoked from driver subroutine  rst.
  139. c
  140.    85    write(iwrite,855)
  141.   855    format(24h1all of the eigenvalues ,
  142.      x     48hof the real symmetric tridiagonal matrix follow. /
  143.      x     38h the path involving  imtql1  was used.  )
  144.          do  86  i = 1,n
  145.    86       w(i) = st(i,2)
  146.          call  rst(nm,n,w,st(1,1),0,z,error)
  147.          write(iwrite,865) error
  148.   865    format(//31h *****error from imtql1***** = ,i4)
  149.          maxeig = 0.0d0
  150.          maxdif = 0.0d0
  151.          m = n
  152.          if( error .ne. 0 ) m = error - 1
  153.          if( m .eq. 0 ) go to 230
  154.          do 87 i = 1,m
  155.             if( dabs(w(i)) .gt. maxeig ) maxeig = dabs(w(i))
  156.             u = dabs(rv5(i) - w(i))
  157.             if( u .gt. maxdif ) maxdif = u
  158.    87    continue
  159.          if( maxeig .eq. 0.0d0 ) maxeig = 1.0d0
  160.          eigdif = maxdif/(maxeig*machep*dfloat(10*n))
  161.          write(iwrite,88) eigdif
  162.    88    format(//49h comparison of the eigenvalues from  imtql1  with,
  163.      x          52h those from  imtql2  gives the normalized difference,
  164.      x          1pd16.8)
  165.          go to  130
  166. c
  167. c     rstw1z  ( usage here computes all the eigenvectors )
  168. c
  169.    89    write(iwrite,891)
  170.   891    format(43h1some of the eigenvalues and corresponding ,
  171.      x     53heigenvectors of the real symmetric tridiagonal matrix /
  172.      x     47h follow.  the path involving  imtqlv  was used.  )
  173.          do 892 i = 1,n
  174.   892       e2(i) = st(i,1) ** 2
  175.          call imtqlv(n,st(1,2),st(1,1),e2,w,ind,error,rv1)
  176.          write(iwrite,893) error
  177.   893    format(//31h *****error from imtqlv***** = ,i4)
  178.          m = n
  179.          if( error .ne. 0 ) m = error - 1
  180.          call tinvit(nm,n,st(1,2),st(1,1),e2,m,w,ind,z,error,rv1,rv2,
  181.      x               rv3,rv4,rv6)
  182.          write(iwrite,99) error
  183.          go to 130
  184. c
  185. c     rst1w1z  using  tsturm
  186. c
  187.    90    write(iwrite,91)
  188.    91    format(43h1some of the eigenvalues and corresponding ,
  189.      x     53heigenvectors of the real symmetric tridiagonal matrix /
  190.      x     47h follow.  the path involving  tsturm  was used.  )
  191.          eps1 = 0.0d0
  192.          do  92  i = 2,n
  193.    92       e2(i) = st(i,1)**2
  194.          call  tsturm(nm,n,eps1,st(1,2),st(1,1),e2,lb,ub,mm,m,w,z,
  195.      x                error,rv1,rv2,rv3,rv4,rv5,rv6)
  196.          write(iwrite,93) error
  197.    93    format(//31h *****error from tsturm***** = ,i4)
  198.          if( error .eq. 3*n + 1 ) go to 230
  199.          if( error .gt. 4*n ) m = error - 4*n - 1
  200.          go to  150
  201. c
  202. c     rst1w1z  using  bisect  and  tinvit
  203. c
  204.    95    write(iwrite,96)
  205.    96    format(43h1some of the eigenvalues and corresponding ,
  206.      x     53heigenvectors of the real symmetric tridiagonal matrix /
  207.      x     47h follow.  the path involving  bisect  was used.  )
  208.          eps1 = 0.0d0
  209.          do  97  i = 2,n
  210.    97       e2(i) = st(i,1)**2
  211.          call  bisect(n,eps1,st(1,2),st(1,1),e2,lb,ub,mm,m,w,ind,error,
  212.      x                rv4,rv5)
  213.          write(iwrite,98) error
  214.    98    format(//31h *****error from bisect***** = ,i4)
  215.          if( error .ne. 0 ) go to 230
  216.          call  tinvit(nm,n,st(1,2),st(1,1),e2,m,w,ind,z,
  217.      x                error,rv1,rv2,rv3,rv4,rv6)
  218.          write(iwrite,99) error
  219.    99    format(//31h *****error from tinvit***** = ,i4)
  220.          go to  150
  221. c
  222. c     rst1w  using  bisect
  223. c
  224.   100    write(iwrite,101)
  225.   101    format(25h1some of the eigenvalues ,
  226.      x     48hof the real symmetric tridiagonal matrix follow.)
  227.          eps1 = 0.0d0
  228.          do  102  i = 2,n
  229.   102       e2(i) = st(i,1)**2
  230.          call  bisect(n,eps1,st(1,2),st(1,1),e2,lb,ub,mm,m,w,ind,error,
  231.      x                rv4,rv5)
  232.          write(iwrite,103) error
  233.   103    format(//31h *****error from bisect***** = ,i4)
  234.          go to  150
  235. c
  236. c     rst1w1z  using  ratqr  and  tinvit
  237. c
  238.   110    write(iwrite,111)
  239.   111    format(43h1some of the eigenvalues and corresponding ,
  240.      x     53heigenvectors of the real symmetric tridiagonal matrix /
  241.      x     46h follow.  the path involving  ratqr  was used.  )
  242.          eps1 = 0.0d0
  243.          do 112  i = 2,n
  244.   112       e2(i) = st(i,1)**2
  245.          call  ratqr(n,eps1,st(1,2),st(1,1),e2,no+m11-1,w,
  246.      x               ind,rv1,.true.,0,error)
  247.          do 1125 i = 1,no
  248.             m = i + m11 - 1
  249.             w(i) = w(m)
  250.             ind(i) = ind(m)
  251.  1125    continue
  252.          write(iwrite,113) error
  253.   113    format(//30h *****error from ratqr***** = ,i4)
  254.          if( error .ne. 0 )  go to  230
  255.          m = no
  256.          call  tinvit(nm,n,st(1,2),st(1,1),e2,m,w,ind,z,error,
  257.      x                rv1,rv2,rv3,rv4,rv6)
  258.          write(iwrite,114) error
  259.   114    format(//31h *****error from tinvit***** = ,i4)
  260.          go to  150
  261. c
  262. c     rst1w  using  tridib
  263. c
  264.   115    write(iwrite,101)
  265.          eps1 = 0.0d0
  266.          do 117  i = 2,n
  267.   117       e2(i) = st(i,1)**2
  268.          call  tridib(n,eps1,st(1,2),st(1,1),e2,lb,ub,m11,no,w,ind,
  269.      x                error,rv4,rv5)
  270.          write(iwrite,118) error
  271.   118    format(//31h *****error from tridib***** = ,i4)
  272.          m = no
  273.          go to 150
  274. c
  275.   130    write(iwrite,140) (i,w(i),i=1,n)
  276.   140    format(///18h eigenvalues of  a/3(i4,2x,1pd24.16,3x))
  277.          if( icall .eq. 3 .or. icall .eq. 4 )  go to  230
  278.          go to  190
  279.   150    write(iwrite,160) m,lb,ub
  280.   160    format(///4h the ,i4,29h  eigenvalues of  a  between ,1pd24.16,
  281.      x        6h  and ,d24.16)
  282.          if( m .eq. 0 )  go to  230
  283.          write(iwrite,163) (i,w(i),i=1,m)
  284.   163    format(//3(i3,2x,1pd24.16,3x))
  285.          if( icall .eq. 8 .or. icall .eq. 10 ) go to 230
  286.   190    call  rsbwzr(nm,n,m,2,st,w,z,norm,resdul)
  287.          tcrit = resdul/(dfloat(10*n)*machep)
  288.          write(iwrite,200) tcrit,resdul 
  289.   200    format(///48h as a guide to evaluating the performance of the,
  290.      x     53h  eispack  codes, a number  x, related to the  1-norm /
  291.      x     56h of the residual matrix, may be normalized to the number,
  292.      x     56h  y  by dividing by the product of the machine precision/
  293.      x     57h machep  and  10*n  where  n  is the order of the matrix.,
  294.      x     36h  the  eispack  codes have performed /
  295.      x     56h (well,satisfactorily,poorly) according as the ratio  y ,
  296.      x     54h is (less than  1, less than  100, greater than  100)./
  297.      x     23h for this run,  y  is =,1pd14.6,11h  with  x =,d14.6)
  298.          do  220  kk = 1,m,3
  299.             kk1 = min0(kk+2,m)
  300.             irs = kk1-kk+1
  301.             go to  (204,208,212),  irs
  302.   204       write(iwrite,205) (w(k),k=kk,kk1),(norm(k),k=kk,kk1),
  303.      x           ((z(i,k),k=kk,kk1),i=1,n)
  304.   205    format(///1x,1(9x,15heigenvalue of a,8x)/1x,1(4x,1pd24.16,4x)//
  305.      x   1x,1(5x,23h1-norm of corresponding,4x)/
  306.      x   1x,1(9x,15hresidual vector,8x)/
  307.      x   1x,1(4x,24h!!ax-l*x!!/(!!x!!*!!a!!),4x)/
  308.      x   1x,1(8x,d16.8,8x)//
  309.      x   1x,1(4x,26hcorresponding eigenvector ,2x)/1(5x,d24.16,3x))
  310.             go to 220
  311.   208       write(iwrite,209) (w(k),k=kk,kk1),(norm(k),k=kk,kk1), 
  312.      x           ((z(i,k),k=kk,kk1),i=1,n)
  313.   209    format(///1x,2(9x,15heigenvalue of a,8x)/1x,2(4x,1pd24.16,4x)//
  314.      x   1x,2(5x,23h1-norm of corresponding,4x)/
  315.      x   1x,2(9x,15hresidual vector,8x)/
  316.      x   1x,2(4x,24h!!ax-l*x!!/(!!x!!*!!a!!),4x)/
  317.      x   1x,2(8x,d16.8,8x)//
  318.      x   1x,2(4x,26hcorresponding eigenvector ,2x)/(5x,d24.16,3x,
  319.      x   5x,d24.16,3x))
  320.             go to 220
  321.   212       write(iwrite,213) (w(k),k=kk,kk1),(norm(k),k=kk,kk1),
  322.      x           ((z(i,k),k=kk,kk1),i=1,n)
  323.   213    format(///1x,3(9x,15heigenvalue of a,8x)/1x,3(4x,1pd24.16,4x)//
  324.      x   1x,3(5x,23h1-norm of corresponding,4x)/
  325.      x   1x,3(9x,15hresidual vector,8x)/
  326.      x   1x,3(4x,24h!!ax-l*x!!/(!!x!!*!!a!!),4x)/
  327.      x   1x,3(8x,d16.8,8x)//
  328.      x   1x,3(4x,26hcorresponding eigenvector ,2x)/(5x,d24.16,3x,
  329.      x   5x,d24.16,3x,5x,d24.16,3x))
  330.   220    continue
  331.   230 continue
  332.       go to  10
  333.       end
  334.