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 / rgeispak.f < prev    next >
Text File  |  1996-09-28  |  15KB  |  384 lines
  1. c
  2. c     this driver tests  eispack  for the class of real general matrices
  3. c     exhibiting the use of  eispack  to find all the eigenvalues and
  4. c     eigenvectors, only all the eigenvalues, or all the eigenvalues and
  5. c     some of the corresponding eigenvectors, using orthogonal or
  6. c     elementary similarity transformations to reduce the full matrix to
  7. c     hessenberg form, with and without balancing of the full matrix.
  8. c
  9. c     this driver is catalogued as  eispdrv4(rgeispak).
  10. c
  11. c     the dimension of  a,z,rm1, and  asave  should be  nm  by  nm.
  12. c     the dimension of  wr,wi,select,slhold,int,scale,ort,norm,rv1,
  13. c     and  rv2  should be  nm.
  14. c     the dimension of  ahold  should be  nm  by  nm.
  15. c     here nm = 20.
  16. c
  17.       double precision a(20,20),z(20,20),asave(20,20),rm1(20,20), 
  18.      x        ahold( 20, 20),ort( 20),wr( 20),wi(20),
  19.      x        scale( 20),rv1( 20),rv2( 20),norm( 20),
  20.      x        tcrit,machep,resdul
  21.       integer  int( 20),low,upp,error,cp,cp1,cp2
  22.       logical  select( 20),slhold( 20),cplx
  23.       data ireadc/5/,iwrite/6/
  24. c
  25. c     .......... machep is a machine dependent parameter specifying
  26. c                the relative precision of floating point arithmetic.
  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. c
  33.       nm = 20
  34.       open(unit=1,file='rgdatai')
  35.       open(unit=5,file='rgdatas')
  36.       rewind 1
  37.       rewind 5
  38.    10 call rmatin(nm,n,a,ahold,0)
  39.       write(iwrite,20) n
  40.    20 format(37h1the real general matrix  a  of order,
  41.      x       i4,23h  is (printed by rows)  )
  42.       do 30 i = 1,n
  43.    30    write(iwrite,40) (a(i,j),j=1,n)
  44.    40    format(/(1x,1p5d24.16))
  45.       read(ireadc,50) mm,(select(i),i=1,n)
  46.    50 format(i4/(72l1))
  47.       do  60  i = 1,n
  48.    60 slhold(i) = select(i)
  49. c
  50. c     mm  and  select  are read from sysin after the matrix is
  51. c     generated.  mm  specifies to  invit  the maximum number of
  52. c     eigenvectors that will be computed.  select  contains information
  53. c     about which eigenvectors are desired.
  54. c
  55.       do  510  icall = 1,12
  56.          icall2 = mod(icall,3)
  57.          if( icall2 .ne. 0 ) go to 80
  58.          do  70  i = 1,n
  59.    70    select(i) = slhold(i)
  60.    80    if( icall .ne. 1 )  call  rmatin(nm,n,a,ahold,1)
  61.          go to  (90,100,110,120,130,140,160,165,170,175,180,185), icall
  62. c
  63. c     rgewz
  64. c     invoked from driver subroutine  rg.
  65. c
  66.    90    write(iwrite,903)
  67.   903    format(41h1the eigenvalues and eigenvectors of the ,
  68.      x     50hreal general matrix follow.  this matrix has been /
  69.      x     56h balanced and has been reduced to hessenberg form using ,
  70.      x     27helementary transformations.  )
  71.          call  rg(nm,n,a,wr,wi,1,z,int,scale,error)
  72.          write(iwrite,905) error
  73.   905    format(//28h *****error from hqr2***** = ,i4)
  74.          go to  230
  75. c
  76. c     rgew
  77. c     invoked from driver subroutine  rg.
  78. c
  79.   100    write(iwrite,1003)
  80.  1003    format(24h1the eigenvalues of the ,
  81.      x     50hreal general matrix follow.  this matrix has been /
  82.      x     56h balanced and has been reduced to hessenberg form using ,
  83.      x     27helementary transformations.  )
  84.          call  rg(nm,n,a,wr,wi,0,a,int,scale,error)
  85.          write(iwrite,1005) error
  86.  1005    format(//27h *****error from hqr***** = ,i4)
  87.          go to  230
  88. c
  89. c     rgew1z
  90. c
  91.   110    write(iwrite,1103)
  92.  1103    format(50h1the eigenvalues and selected eigenvectors of the ,
  93.      x     50hreal general matrix follow.  this matrix has been /
  94.      x     56h balanced and has been reduced to hessenberg form using ,
  95.      x     27helementary transformations.  )
  96.          call  balanc(nm,n,a,low,upp,scale)
  97.          call  elmhes(nm,n,low,upp,a,int)
  98.          do  112  i = 1,n
  99.             do  112  j = 1,n
  100.   112         asave(i,j) = a(i,j)
  101.          call  hqr(nm,n,low,upp,a,wr,wi,error)
  102.          write(iwrite,113) error,low,upp
  103.   113    format(//35h *****error from hqr,low,upp***** = ,3i4)
  104.          call  invit(nm,n,asave,wr,wi,select,mm,m,z,error,rm1,rv1,rv2)
  105.          write(iwrite,114) error
  106.   114    format(//29h *****error from invit***** = ,i4)
  107.          call  elmbak(nm,low,upp,asave,int,m,z)
  108.          call  balbak(nm,n,low,upp,scale,m,z)
  109.          go to  230
  110. c
  111. c     rgnewz
  112. c
  113.   120    write(iwrite,1203)
  114.  1203    format(41h1the eigenvalues and eigenvectors of the ,
  115.      x     54hreal general matrix follow.  this matrix has not been /
  116.      x     56h balanced and has been reduced to hessenberg form using ,
  117.      x     27helementary transformations.  )
  118.          call  elmhes(nm,n,1,n,a,int)
  119.          call  eltran(nm,n,1,n,a,int,z)
  120.          call  hqr2(nm,n,1,n,a,wr,wi,z,error)
  121.          write(iwrite,1205) error
  122.  1205    format(//28h *****error from hqr2***** = ,i4)
  123.          go to  230
  124. c
  125. c     rgnew
  126. c
  127.   130    write(iwrite,1303)
  128.  1303    format(24h1the eigenvalues of the ,
  129.      x     54hreal general matrix follow.  this matrix has not been /
  130.      x     56h balanced and has been reduced to hessenberg form using ,
  131.      x     27helementary transformations.  )
  132.          call  elmhes(nm,n,1,n,a,int)
  133.          call  hqr(nm,n,1,n,a,wr,wi,error)
  134.          write(iwrite,1305) error
  135.  1305    format(//27h *****error from hqr***** = ,i4)
  136.          go to  230
  137. c
  138. c     rgnew1z
  139. c
  140.   140    write(iwrite,1403)
  141.  1403    format(50h1the eigenvalues and selected eigenvectors of the ,
  142.      x     54hreal general matrix follow.  this matrix has not been /
  143.      x     56h balanced and has been reduced to hessenberg form using ,
  144.      x     27helementary transformations.  )
  145.          call  elmhes(nm,n,1,n,a,int)
  146.          do  142  i = 1,n
  147.             do  142  j = 1,n
  148.   142         asave(i,j) = a(i,j)
  149.          call  hqr(nm,n,1,n,a,wr,wi,error)
  150.          write(iwrite,143) error
  151.   143    format(//27h *****error from hqr***** = ,i4)
  152.          call  invit(nm,n,asave,wr,wi,select,mm,m,z,error,rm1,rv1,rv2)
  153.          write(iwrite,144) error
  154.   144    format(//29h *****error from invit***** = ,i4)
  155.          call  elmbak(nm,1,n,asave,int,m,z)
  156.          go to 230
  157. c
  158. c     rgowz
  159. c
  160.   160    write(iwrite,1603)
  161.  1603    format(41h1the eigenvalues and eigenvectors of the ,
  162.      x     50hreal general matrix follow.  this matrix has been /
  163.      x     56h balanced and has been reduced to hessenberg form using ,
  164.      x     27horthogonal transformations.  )
  165.          call  balanc(nm,n,a,low,upp,scale)
  166.          call  orthes(nm,n,low,upp,a,ort)
  167.          call  ortran(nm,n,low,upp,a,ort,z)
  168.          call  hqr2(nm,n,low,upp,a,wr,wi,z,error)
  169.          write(iwrite,1605) error,low,upp
  170.  1605    format(//36h *****error from hqr2,low,upp***** = ,3i4)
  171.          call  balbak(nm,n,low,upp,scale,n,z)
  172.          go to  230
  173. c
  174. c     rgow
  175. c
  176.   165    write(iwrite,1653)
  177.  1653    format(24h1the eigenvalues of the ,
  178.      x     50hreal general matrix follow.  this matrix has been /
  179.      x     56h balanced and has been reduced to hessenberg form using ,
  180.      x     27horthogonal transformations.  )
  181.          call  balanc(nm,n,a,low,upp,scale)
  182.          call  orthes(nm,n,low,upp,a,ort)
  183.          call  hqr(nm,n,low,upp,a,wr,wi,error)
  184.          write(iwrite,1655) error,low,upp
  185.  1655    format(//35h *****error from hqr,low,upp***** = ,3i4)
  186.          go to  230
  187. c
  188. c     rgow1z
  189. c
  190.   170    write(iwrite,1703)
  191.  1703    format(50h1the eigenvalues and selected eigenvectors of the ,
  192.      x     50hreal general matrix follow.  this matrix has been /
  193.      x     56h balanced and has been reduced to hessenberg form using ,
  194.      x     27horthogonal transformations.  )
  195.          call  balanc(nm,n,a,low,upp,scale)
  196.          call  orthes(nm,n,low,upp,a,ort)
  197.          do  172  i = 1,n
  198.             do  172  j = 1,n
  199.   172         asave(i,j) = a(i,j)
  200.          call  hqr(nm,n,low,upp,a,wr,wi,error)
  201.          write(iwrite,173) error,low,upp
  202.   173    format(//35h *****error from hqr,low,upp***** = ,3i4)
  203.          call  invit(nm,n,asave,wr,wi,select,mm,m,z,error,rm1,rv1,rv2)
  204.          write(iwrite,174) error
  205.   174    format(//29h *****error from invit***** = ,i3)
  206.          call  ortbak(nm,low,upp,asave,ort,m,z)
  207.          call  balbak(nm,n,low,upp,scale,m,z)
  208.          go to  230
  209. c
  210. c     rgnowz
  211. c
  212.   175    write(iwrite,1753)
  213.  1753    format(41h1the eigenvalues and eigenvectors of the ,
  214.      x     54hreal general matrix follow.  this matrix has not been /
  215.      x     56h balanced and has been reduced to hessenberg form using ,
  216.      x     27horthogonal transformations.  )
  217.          call  orthes(nm,n,1,n,a,ort)
  218.          call  ortran(nm,n,1,n,a,ort,z)
  219.          call  hqr2(nm,n,1,n,a,wr,wi,z,error)
  220.          write(iwrite,1755) error
  221.  1755    format(//28h *****error from hqr2*****  ,i4)
  222.          go to  230
  223. c
  224. c     rgnow
  225. c
  226.   180    write(iwrite,1803)
  227.  1803    format(24h1the eigenvalues of the ,
  228.      x     54hreal general matrix follow.  this matrix has not been /
  229.      x     56h balanced and has been reduced to hessenberg form using ,
  230.      x     27horthogonal transformations.  )
  231.          call  orthes(nm,n,1,n,a,ort)
  232.          call  hqr(nm,n,1,n,a,wr,wi,error)
  233.         write(iwrite,1805) error
  234.  1805    format(//27h *****error from hqr***** = ,i3)
  235.          go to  230
  236. c
  237. c     rgnow1z
  238. c
  239.   185    write(iwrite,1853)
  240.  1853    format(50h1the eigenvalues and selected eigenvectors of the ,
  241.      x     54hreal general matrix follow.  this matrix has not been /
  242.      x     56h balanced and has been reduced to hessenberg form using ,
  243.      x     27horthogonal transformations.  )
  244.          call  orthes(nm,n,1,n,a,ort)
  245.          do  187  i = 1,n
  246.             do  187  j = 1,n
  247.   187         asave(i,j) = a(i,j)
  248.          call  hqr(nm,n,1,n,a,wr,wi,error)
  249.          write(iwrite,188) error
  250.   188    format(//27h *****error from hqr***** = ,i3)
  251.          call  invit(nm,n,asave,wr,wi,select,mm,m,z,error,rm1,rv1,rv2)
  252.          write(iwrite,189) error
  253.   189    format(//29h *****error from invit***** = ,i3)
  254.          call  ortbak(nm,1,n,asave,ort,m,z)
  255. c
  256.   230    if( icall2 .eq. 0 )  go to  250
  257.          write(iwrite,240) (i,wr(i),wi(i),i=1,n)
  258.   240    format(//29h all the eigenvalues computed/
  259.      x        2(i3,3x,1p2d24.16,2x))
  260.          go to  270
  261.   250    write(iwrite,260) (i,select(i),wr(i),wi(i),i=1,n)
  262.   260    format(//51h all the eigenvalues computed and those eigenvalues
  263.      x            ,9h selected/2(i3,1x,l1,1x,1p2d24.16,2x))
  264.   270    if( error .gt. 0 )  go to  510
  265.          if( icall2 .eq. 2 )  go to  510
  266.          call  rmatin(nm,n,a,ahold,1)
  267.          if( icall2 .ne. 0 )  go to  280
  268.          call  rgw1zr(nm,n,a,wr,wi,select,m,z,norm,resdul)
  269.          go to  290
  270.   280    call  rgwzr(nm,n,a,wr,wi,z,norm,resdul)
  271.   290    tcrit = resdul/(dfloat(10*n)*machep)
  272.          write(iwrite,300) tcrit,resdul
  273.   300    format(///48h as a guide to evaluating the performance of the,
  274.      x     53h  eispack  codes, a number  x, related to the  1-norm /
  275.      x     56h of the residual matrix, may be normalized to the number,
  276.      x     56h  y  by dividing by the product of the machine precision/
  277.      x     57h machep  and  10*n  where  n  is the order of the matrix.,
  278.      x     36h  the  eispack  codes have performed /
  279.      x     56h (well,satisfactorily,poorly) according as the ratio  y ,
  280.      x     54h is (less than  1, less than  100, greater than  100)./
  281.      x     23h for this run,  y  is =,1pd14.6,11h  with  x =,d14.6)
  282.          nm2 = n-2
  283.          k = 1
  284.          if( icall2 .eq. 0 )  go to  390
  285.   310      if( k .gt. n )  go to  510
  286.            if( wi(k) .eq. 0.0d0 )  go to  360
  287.            if( k .gt. nm2 )  go to  340
  288.           if( wi(k+2) .eq. 0.0d0 )  go to  330
  289. c          both digenvectors are complex.
  290.            write(iwrite,320)
  291.      x      wr(k),wi(k),wr(k+2),wi(k+2),norm(k),norm(k+2), 
  292.      x      (z(i,k),z(i,k+1),z(i,k+2),z(i,k+3),i=1,n)
  293.   320      format(///1x,2(19x,15heigenvalue of a,20x)/1x,2(1p2d24.16,6x)
  294.      x       //1x,2(8x,39h1-norm of corresponding residual vector,7x)/
  295.      x       1x,2(6x,26h!!ax-l*x!!/(!!x!!*!!a!!) =,d16.8,6x)//
  296.      x       1x,2(14x,25hcorresponding eigenvector,15x)/(1x,
  297.      x       2d24.16,6x,2d24.16,6x))
  298.            k = k+4
  299.            go to  310
  300. c          one complex eigenvector and one real eigenvector in that
  301. c          order.
  302.   330      write(iwrite,320)
  303.      x                wr(k+2),wi(k+2),wr(k),wi(k),norm(k+2),norm(k), 
  304.      x                (z(i,k+2),wi(k+2),z(i,k),z(i,k+1),i=1,n)
  305.            k = k+3
  306.            go to  310
  307. c          one complex eigenvector.
  308.   340      write(iwrite,350) wr(k),wi(k),norm(k),(z(i,k),z(i,k+1),i=1,n)
  309.   350      format(///1x,19x,15heigenvalue of a,20x/1x,1p2d24.16,6x//
  310.      x       1x,8x,39h1-norm of corresponding residual vector,7x/
  311.      x       1x,6x,26h!!ax-l*x!!/(!!x!!*!!a!!) =,d16.8,6x//
  312.      x       1x,14x,25hcorresponding eigenvector,15x/(1x,2d24.16,6x))
  313.            go to  510
  314.   360      if( k .eq. n )  go to  380
  315.            if( wi(k+1) .eq. 0.0d0 )  go to  370
  316. c          one real eigenvector and one complex eigenvector in that
  317. c          order.
  318.            write(iwrite,320)
  319.      x                wr(k),wi(k),wr(k+1),wi(k+1),norm(k),norm(k+1), 
  320.      x                (z(i,k),wi(k),z(i,k+1),z(i,k+2),i=1,n)
  321.            k = k+3
  322.            go to  310
  323. c          both eigenvectors are real.
  324.   370      write(iwrite,320)
  325.      x                wr(k),wi(k),wr(k+1),wi(k+1),norm(k),norm(k+1), 
  326.      x                (z(i,k),wi(k),z(i,k+1),wi(k+1),i=1,n)
  327.            k = k+2
  328.            go to  310
  329. c          one real eigenvector.
  330.   380      write(iwrite,350) wr(k),wi(k),norm(k),(z(i,k),wi(k),i=1,n)
  331.            go to  510
  332.   390    l = 1
  333.          cplx = .false.
  334.   400    cp1 = -1
  335.   410      if( k .gt. n ) if( cp1 )  510, 490, 500
  336.            if( .not. cplx )  go to  412
  337.              cplx = .false.
  338.              go to  415
  339.   412      if( wi(k) .ne. 0.0d0 )  cplx = .true.
  340.   415      kold = k
  341.            k = k+1
  342.            if( .not. select(kold) )  go to  410
  343.            lold = l
  344.            l = l+1
  345.            cp = 0
  346.            if( wi(kold) .eq. 0.0d0 )  go to  420
  347.              cp = 1
  348.              if( cplx )  k = k+1
  349.              l = l+1
  350.   420      if( cp1 )  430, 440, 440
  351.   430        lp1 = lold
  352.              kp1 = kold
  353.              cp1 = cp
  354.              go to  410
  355.   440        lp2 = lold
  356.              kp2 = kold
  357.              cp2 = cp
  358.            cp = cp1*2 + cp2 + 1
  359.            go to  (450,460,470,480),  cp
  360.   450      write(iwrite,320) wr(kp1),wi(kp1),wr(kp2),wi(kp2),norm(kp1),
  361.      x         norm(kp2),(z(i,lp1),wi(kp1),z(i,lp2),wi(kp2),i=1,n)
  362.            go to  400
  363.   460      write(iwrite,320) wr(kp1),wi(kp1),wr(kp2),wi(kp2),norm(kp1),
  364.      x        norm(kp2),(z(i,lp1),wi(kp1),z(i,lp2),z(i,lp2+1),i=1,n)
  365.            go to  400
  366.   470      kold = kp1
  367.            kp1 = kp2
  368.            kp2 = kold
  369.            lold = lp1
  370.            lp1 = lp2
  371.            lp2 = lold
  372.            go to  460
  373.   480      write(iwrite,320) wr(kp1),wi(kp1),wr(kp2),wi(kp2),norm(kp1),
  374.      x     norm(kp2),(z(i,lp1),z(i,lp1+1),z(i,lp2),z(i,lp2+1),i=1,n)
  375.            go to  400
  376.   490      write(iwrite,350) wr(kp1),wi(kp1),norm(kp1),
  377.      x                 (z(i,lp1),wi(kp1),i=1,n)
  378.            go to  510
  379.   500      write(iwrite,350) wr(kp1),wi(kp1),norm(kp1),
  380.      x                 (z(i,lp1),z(i,lp1+1),i=1,n)
  381.   510 continue
  382.       go to  10
  383.       end
  384.