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 / cgeispak.f < prev    next >
Text File  |  1996-09-28  |  13KB  |  322 lines
  1. c
  2. c     this driver tests  eispack  for the class of complex general ma-
  3. c     trices exhibiting the use of  eispack  to find all the eigenvalues
  4. c     and eigenvectors, only all the eigenvalues, or all the eigen-
  5. c     values and some of the corresponding eigenvectors, using unitary
  6. c     or elementary similarity transformations to reduce the full matrix
  7. c     to hessenberg form, with and without balancing of the full matrix.
  8. c
  9. c     this driver is catalogued as  eispdrv4(cgeispak).
  10. c
  11. c     the dimension of  ar,ai,zr,zi,asaver,asavei,rm1, and  rm2
  12. c     should be  nm  by  nm.
  13. c     the dimension of  wr,wi,select,slhold,int,scale,ortr,orti,norm,
  14. c     rv1,  and  rv2  should be  nm.
  15. c     the dimension of  arhold  and  aihold  should be  nm  by  nm.
  16. c     here nm = 20.
  17. c
  18.       double precision ar( 20, 20),ai( 20, 20),zr( 20, 20),zi( 20, 20),
  19.      x        asaver( 20, 20),asavei( 20, 20),rm1( 20, 20),rm2( 20,20),
  20.      x        ortr( 20),orti( 20),wr( 20),wi( 20),
  21.      x        scale( 20),rv1( 20),rv2( 20),norm( 20),
  22.      x        tcrit,machep,resdul
  23.       double precision arhold( 20, 20),aihold( 20, 20)      
  24.       integer  int( 20),upp,error
  25.       logical  select( 20),slhold( 20)
  26.       data ireadc/5/,iwrite/6/
  27. c
  28. c     .......... machep is a machine dependent parameter specifying
  29. c                the relative precision of floating point arithmetic.
  30. c
  31.       machep = 1.0d0
  32.     1 machep = 0.5d0*machep
  33.       if (1.0d0 + machep .gt. 1.0d0 ) go to 1
  34.       machep = 2.0d0*machep
  35. c
  36. c
  37.       nm = 20
  38.    10 call cmatin(nm,n,ar,ai,arhold,aihold,0)
  39.       write(iwrite,20) n 
  40.    20 format(50h1the complex general matrix  a = (ar,ai)  of order,
  41.      x       i4,22h  is (printed by rows)  )
  42.       do  30  i = 1,n
  43.    30    write(iwrite,40) (ar(i,j),ai(i,j),j=1,n)
  44.    40    format(/2(1p2d24.16,4x))
  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  300  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  cmatin(nm,n,ar,ai,arhold,aihold,1)
  61.          go to  (90,100,110,120,130,140,160,165,170,175,180,185), icall
  62. c
  63. c     cgewz
  64. c
  65.    90    write(iwrite,903)
  66.   903    format(41h1the eigenvalues and eigenvectors of the ,
  67.      x     53hcomplex general matrix follow.  this matrix has been /
  68.      x     56h balanced and has been reduced to hessenberg form using ,
  69.      x     27helementary transformations.  )
  70.          call  cbal(nm,n,ar,ai,low,upp,scale)
  71.          call  comhes(nm,n,low,upp,ar,ai,int)
  72.          call  comlr2(nm,n,low,upp,int,ar,ai,wr,wi,zr,zi,error)
  73.          write(iwrite,905) error,low,upp 
  74.   905    format(//38h *****error from comlr2,low,upp***** =,3i4)
  75.          call  cbabk2(nm,n,low,upp,scale,n,zr,zi)
  76.          go to  190
  77. c
  78. c     cgew
  79. c
  80.   100    write(iwrite,1003)
  81.  1003    format(24h1the eigenvalues of the ,
  82.      x     53hcomplex general matrix follow.  this matrix has been /
  83.      x     56h balanced and has been reduced to hessenberg form using ,
  84.      x     27helementary transformations.  )
  85.          call  cbal(nm,n,ar,ai,low,upp,scale)
  86.          call  comhes(nm,n,low,upp,ar,ai,int)
  87.          call  comlr(nm,n,low,upp,ar,ai,wr,wi,error)
  88.          write(iwrite,1005) error,low,upp  
  89.  1005    format(//37h *****error from comlr,low,upp***** = ,3i4)
  90.          go to  190
  91. c
  92. c     cgew1z
  93. c
  94.   110    write(iwrite,1103)
  95.  1103    format(50h1the eigenvalues and selected eigenvectors of the ,
  96.      x     53hcomplex general matrix follow.  this matrix has been /
  97.      x     56h balanced and has been reduced to hessenberg form using ,
  98.      x     27helementary transformations.  )
  99.          call  cbal(nm,n,ar,ai,low,upp,scale)
  100.          call  comhes(nm,n,low,upp,ar,ai,int)
  101.          do  112  i = 1,n
  102.             do  112  j = 1,n
  103.               asaver(i,j) = ar(i,j)
  104.   112         asavei(i,j) = ai(i,j)
  105.          call  comlr(nm,n,low,upp,ar,ai,wr,wi,error)
  106.          write(iwrite,113) error,low,upp  
  107.   113    format(//37h *****error from comlr,low,upp***** = ,3i4)
  108.          call  cinvit(nm,n,asaver,asavei,wr,wi,select,mm,m,zr,zi,
  109.      x                error,rm1,rm2,rv1,rv2)
  110.          write(iwrite,114) error
  111.   114    format(//30h *****error from cinvit***** = ,i4)
  112.          call  combak(nm,low,upp,asaver,asavei,int,m,zr,zi)
  113.          call  cbabk2(nm,n,low,upp,scale,m,zr,zi)
  114.          go to  190
  115. c
  116. c     cgnewz
  117. c
  118.   120    write(iwrite,1203)
  119.  1203    format(41h1the eigenvalues and eigenvectors of the ,
  120.      x     57hcomplex general matrix follow.  this matrix has not been /
  121.      x     56h balanced and has been reduced to hessenberg form using ,
  122.      x     27helementary transformations.  )
  123.          call  comhes(nm,n,1,n,ar,ai,int)
  124.          call  comlr2(nm,n,1,n,int,ar,ai,wr,wi,zr,zi,error)
  125.          write(iwrite,1205) error
  126.  1205    format(//30h *****error from comlr2***** = ,i4)
  127.          go to  190
  128. c
  129. c     cgnew
  130. c
  131.   130    write(iwrite,1303)
  132.  1303    format(24h1the eigenvalues of the ,
  133.      x     57hcomplex general matrix follow.  this matrix has not been /
  134.      x     56h balanced and has been reduced to hessenberg form using ,
  135.      x     27helementary transformations.  )
  136.          call  comhes(nm,n,1,n,ar,ai,int)
  137.          call  comlr(nm,n,1,n,ar,ai,wr,wi,error)
  138.          write(iwrite,1305) error
  139.  1305    format(//29h *****error from comlr***** = ,i4)
  140.          go to  190
  141. c
  142. c     cgnew1z
  143. c
  144.   140    write(iwrite,1403)
  145.  1403    format(50h1the eigenvalues and selected eigenvectors of the ,
  146.      x     57hcomplex general matrix follow.  this matrix has not been /
  147.      x     56h balanced and has been reduced to hessenberg form using ,
  148.      x     27helementary transformations.  )
  149.          call  comhes(nm,n,1,n,ar,ai,int)
  150.          do  142  i = 1,n
  151.             do  142  j = 1,n
  152.               asaver(i,j) = ar(i,j)
  153.   142         asavei(i,j) = ai(i,j)
  154.          call  comlr(nm,n,1,n,ar,ai,wr,wi,error)
  155.          write(iwrite,143) error
  156.   143    format(//29h *****error from comlr***** = ,i4)
  157.          call  cinvit(nm,n,asaver,asavei,wr,wi,select,mm,m,zr,zi,
  158.      x                error,rm1,rm2,rv1,rv2)
  159.          write(iwrite,144) error
  160.   144    format(//30h *****error from cinvit***** = ,i4)
  161.          call  combak(nm,1,n,asaver,asavei,int,m,zr,zi)
  162.          go to  190
  163. c
  164. c     cguwz
  165. c     invoked from driver subroutine  cg.
  166. c
  167.   160    write(iwrite,1603)
  168.  1603    format(41h1the eigenvalues and eigenvectors of the ,
  169.      x     53hcomplex general matrix follow.  this matrix has been /
  170.      x     56h balanced and has been reduced to hessenberg form using ,
  171.      x     24hunitary transformations.  )
  172.          call  cg(nm,n,ar,ai,wr,wi,1,zr,zi,scale,ortr,orti,error)
  173.          write(iwrite,1605) error
  174.  1605    format(//30h *****error from comqr2***** = ,i4)
  175.          go to  190
  176. c
  177. c     cguw
  178. c     invoked from driver subroutine  cg.
  179. c
  180.   165    write(iwrite,1653)
  181.  1653    format(24h1the eigenvalues of the ,
  182.      x     53hcomplex general matrix follow.  this matrix has been /
  183.      x     56h balanced and has been reduced to hessenberg form using ,
  184.      x     24hunitary transformations.  )
  185.          call  cg(nm,n,ar,ai,wr,wi,0,ar,ai,scale,ortr,orti,error)
  186.          write(iwrite,1655) error
  187.  1655    format(//29h *****error from comqr***** = ,i4)
  188.          go to  190
  189. c
  190. c     cguw1z
  191. c
  192.   170    write(iwrite,1703)
  193.  1703    format(50h1the eigenvalues and selected eigenvectors of the ,
  194.      x     53hcomplex general matrix follow.  this matrix has been /
  195.      x     56h balanced and has been reduced to hessenberg form using ,
  196.      x     24hunitary transformations.  )
  197.          call  cbal(nm,n,ar,ai,low,upp,scale)
  198.          call  corth(nm,n,low,upp,ar,ai,ortr,orti)
  199.          do  172  i = 1,n
  200.             do  172  j = 1,n
  201.               asaver(i,j) = ar(i,j)
  202.   172         asavei(i,j) = ai(i,j)
  203.          call  comqr(nm,n,low,upp,ar,ai,wr,wi,error)
  204.          write(iwrite,173) error,low,upp
  205.   173    format(//37h *****error from comqr,low,upp***** = ,3i4)
  206.          call  cinvit(nm,n,asaver,asavei,wr,wi,select,mm,m,zr,zi,
  207.      x                error,rm1,rm2,rv1,rv2)
  208.          write(iwrite,174) error
  209.   174    format(//30h *****error from cinvit***** = ,i3)
  210.          call  cortb(nm,low,upp,asaver,asavei,ortr,orti,m,zr,zi)
  211.          call  cbabk2(nm,n,low,upp,scale,m,zr,zi)
  212.          go to  190
  213. c
  214. c     cgnuwz
  215. c
  216.   175    write(iwrite,1753)
  217.  1753    format(41h1the eigenvalues and eigenvectors of the ,
  218.      x     57hcomplex general matrix follow.  this matrix has not been /
  219.      x     56h balanced and has been reduced to hessenberg form using ,
  220.      x     24hunitary transformations.  )
  221.          call  corth(nm,n,1,n,ar,ai,ortr,orti)
  222.          call  comqr2(nm,n,1,n,ortr,orti,ar,ai,wr,wi,zr,zi,error)
  223.          write(iwrite,1755) error
  224.  1755    format(//30h *****error from comqr2***** = ,i4)
  225.          go to  190
  226. c
  227. c     cgnuw
  228. c
  229.   180    write(iwrite,1803)
  230.  1803 format(24h1the eigenvalues of the ,
  231.      x  57hcomplex general matrix follow.  this matrix has not been /
  232.      x  56h balanced and has been reduced to hessenberg form using ,
  233.      x  24hunitary transformations.  )
  234.          call  corth(nm,n,1,n,ar,ai,ortr,orti)
  235.          call  comqr(nm,n,1,n,ar,ai,wr,wi,error)
  236.          write(iwrite,1805) error
  237.  1805    format(//29h *****error from comqr***** = ,i3)
  238.          go to  190
  239. c
  240. c     cgnuw1z
  241. c
  242.   185    write(iwrite,1853)
  243.  1853    format(50h1the eigenvalues and selected eigenvectors of the ,
  244.      x     57hcomplex general matrix follow.  this matrix has not been /
  245.      x     56h balanced and has been reduced to hessenberg form using ,
  246.      x     24hunitary transformations.  )
  247.          call  corth(nm,n,1,n,ar,ai,ortr,orti)
  248.          do  187  i = 1,n
  249.             do  187  j = 1,n
  250.               asaver(i,j) = ar(i,j)
  251.   187         asavei(i,j) = ai(i,j)
  252.          call  comqr(nm,n,1,n,ar,ai,wr,wi,error)
  253.          write(iwrite,188) error
  254.   188    format(//29h *****error from comqr***** = ,i3)
  255.          call  cinvit(nm,n,asaver,asavei,wr,wi,select,mm,m,zr,zi,
  256.      x                error,rm1,rm2,rv1,rv2)
  257.          write(iwrite,189) error
  258.   189    format(//30h *****error from cinvit***** = ,i3)
  259.          call  cortb(nm,1,n,asaver,asavei,ortr,orti,m,zr,zi)
  260. c
  261.   190    if( icall2 .eq. 0 ) write(iwrite,195)
  262.      x                                   (i,select(i),wr(i),wi(i),i=1,n)
  263.   195    format(//51h all the eigenvalues computed and those eigenvalues
  264.      x          ,9h selected/2(i3,1x,l1,1x,1p2d24.16,2x))
  265.          if( icall2 .ne. 0 ) write(iwrite,200)(i,wr(i),wi(i),i=1,n)
  266.   200    format(//29h all the eigenvalues computed/
  267.      x          2(i3,3x,1p2d24.16,2x))
  268.          if( error .gt. 0 )  go to  300
  269.          if( icall2 .eq. 2 )  go to  300
  270.          call  cmatin(nm,n,ar,ai,arhold,aihold,1)
  271.          if( icall2 .eq. 0 ) call  cgw1zr(nm,n,ar,ai,wr,wi,select,
  272.      x                                    m,zr,zi,norm,resdul)
  273.          if( icall2 .eq. 1 ) call  cgwzr(nm,n,ar,ai,wr,wi,zr,zi,
  274.      x                                   norm,resdul)
  275.          tcrit = resdul/(dfloat(10*n)*machep)
  276.          write(iwrite,220) tcrit,resdul
  277.   220    format(///48h as a guide to evaluating the performance of the,
  278.      x     53h  eispack  codes, a number  x, related to the  1-norm /
  279.      x     56h of the residual matrix, may be normalized to the number,
  280.      x     56h  y  by dividing by the product of the machine precision/
  281.      x     57h machep  and  10*n  where  n  is the order of the matrix.,
  282.      x     36h  the  eispack  codes have performed /
  283.      x     56h (well,satisfactorily,poorly) according as the ratio  y ,
  284.      x     54h is (less than  1, less than  100, greater than  100)./
  285.      x     23h for this run,  y  is =,1pd14.6,11h  with  x =,d14.6)
  286.          if( icall2 .eq. 0 )  go to  270
  287.          do  260  k = 1,n,2
  288.             if( k .eq. n ) go to  240
  289.             write(iwrite,230)
  290.      x                    wr(k),wi(k),wr(k+1),wi(k+1),norm(k),norm(k+1),
  291.      x            (zr(i,k),zi(i,k),zr(i,k+1),zi(i,k+1),i=1,n)
  292.   230      format(///1x,2(19x,15heigenvalue of a,20x)/1x,2(1p2d24.16,6x)
  293.      x       //1x,2(8x,39h1-norm of corresponding residual vector,7x)/
  294.      x       1x,2(6x,26h!!ax-l*x!!/(!!x!!*!!a!!) =,d16.8,6x)//
  295.      x       1x,2(14x,25hcorresponding eigenvector,15x)/( 1x,
  296.      x       2d24.16,6x,2d24.16,6x))
  297.             go to  260
  298.   240       write(iwrite,250)wr(k),wi(k),norm(k),(zr(i,k),zi(i,k),i=1,n)
  299.   250      format(///1x,19x,15heigenvalue of a,20x/1x,1p2d24.16,6x//
  300.      x       1x,8x,39h1-norm of corresponding residual vector,7x/
  301.      x       1x,6x,26h!!ax-l*x!!/(!!x!!*!!a!!) =,d16.8,6x//
  302.      x       1x,14x,25hcorresponding eigenvector,15x/(1x,2d24.16,6x))
  303.   260       continue
  304.          go to  300
  305.   270    l = 1
  306.          k = 0
  307.          do  290  kk = 1,n
  308.            if( .not. select(kk) )  go to  290
  309.            if( k .ne. 0 )  go to  280
  310.            k = kk
  311.            go to  290
  312.   280      write(iwrite,230) wr(k),wi(k),wr(kk),wi(kk),norm(k),norm(kk),
  313.      x                   (zr(i,l),zi(i,l),zr(i,l+1),zi(i,l+1),i=1,n)
  314.            k = 0
  315.            l = l+2
  316.   290      continue
  317.          if( k .ne. 0 )  write(iwrite,250) wr(k),wi(k),norm(k),
  318.      x                       (zr(i,l),zi(i,l),i=1,n)
  319.   300 continue
  320.       go to  10
  321.       end
  322.