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