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