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