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 / rbleispk.f < prev    next >
Text File  |  1996-09-28  |  4KB  |  114 lines
  1. c
  2. c     this driver tests  eispack  for the class of real
  3. c     band matrices exhibiting the use of  eispack  to find the
  4. c     solution to the equation  a*x = b .
  5. c
  6. c     this driver is catalogued as  eispdrv4(rbleispk).
  7. c
  8. c     the dimension of  a  and  ahold  should be  nm  by  mbb.
  9. c     the dimension of  x,b,  and  bhold  should be  nm  by  nm.
  10. c     the dimension of  resdul  should be  nm.
  11. c     the dimension of  rv  should be  nm*mbb.
  12. c     the dimension of  rv1  and  rv6  should be  nm.
  13. c     here  nm = 20  and  mbb = 39.
  14. c
  15. c
  16.       double precision a( 20, 39),x( 20, 20),b( 20, 20),rv( 780),
  17.      x        resdul( 20),rv1( 20),rv6( 20),machep,det,tcrit,
  18.      x        resmax,e21
  19.       double precision ahold( 20, 39),bhold( 20, 20)
  20.       integer error
  21.       data ireadc/5/,iwrite/6/
  22. c
  23. c     .......... machep is a machine dependent parameter specifying
  24. c                the relative precision of floating point arithmetic.
  25.       machep = 1.0d0
  26.     1 machep = 0.5d0*machep
  27.       if (1.0d0 + machep .gt. 1.0d0 ) go to 1
  28.       machep = 2.0d0*machep
  29. c
  30. c
  31.       nm = 20
  32.       mbb = 39
  33.    10 call  rmatin(nm,mbb,n,mbw,a,ahold,bhold,0,0)
  34.       call  rmatin(nm,nm,n,ip,x,ahold,bhold,1,0)
  35.       read(ireadc,15) ie21
  36.    15 format(i6)
  37.       e21 = dfloat(ie21)
  38. c
  39.       if( e21 .eq. 1.0d0 ) go to 19
  40.       write(iwrite,17) n,mbw
  41.    17 format(53h1the highest superdiagonal to the lowest subdiagonal ,
  42.      x       31hof the band matrix  a  of order ,i4,2x,
  43.      x       19hand full band width,i4,4h  is )
  44.       mbh = (mbw + 1)/2 - 1
  45.       do 18 i = 1,mbh
  46.          ibw = mbw - i + 1
  47.          k = n - (mbh - i + 1)
  48.          write(iwrite,30)(a(j,ibw),j=1,k)
  49.    18 continue
  50.       mbh = mbh + 1
  51.       go to 23
  52.    19 write(iwrite,20) n,mbw
  53.    20 format(40h1the diagonal to the lowest subdiagonal ,
  54.      x       41hof the band symmetric matrix  a  of order ,i4,2x,
  55.      x       19hand half band width,i4,4h  is )
  56.    22 mbh = mbw
  57.    23 do 40 jj = 1,mbh
  58.          j = mbh -jj + 1
  59.          write(iwrite,30)(a(i,j),i=jj,n)
  60.    30    format(/(1x,1p5d24.16))
  61.    40 continue
  62.       write(iwrite,50) n,ip
  63.    50 format(///35h the matrix  b  with row dimension  ,i4,
  64.      x       23h  and column dimension  ,i4,18h  printed by rows.)
  65.       do 70 i = 1,n
  66.          write(iwrite,60) (x(i,j),j=1,ip)
  67.    60    format(/(1x,1p5d24.16))
  68.    70 continue
  69. c
  70. c     rbl  using  bandv
  71. c
  72.       do 905 i = 1,n
  73.          rv1(i) = 0.0d0
  74.   905 continue
  75.       call  bandv(nm,n,mbw,a,e21,ip,rv1,x,error,780,rv,rv6)
  76.       write(iwrite,91) error
  77.    91 format(//30h *****error from bandv***** = ,i4)
  78.       det = 1.0d0
  79.       do 92 i = 1,n
  80.          det = dabs(rv(i))*det
  81.    92 continue
  82.       write(iwrite,93) det
  83.    93 format(//46h the unsigned determinant of the matrix  a  = ,
  84.      x       1pd24.16)
  85.       call  rmatin(nm,mbb,n,mbw,a,ahold,bhold,0,1)
  86.       call  rmatin(nm,nm,n,ip,b,ahold,bhold,1,1)
  87.       call  rblres(nm,n,mbw,ip,a,b,x,resdul,e21)
  88.       resmax = 0.0d0
  89.       do 935 i = 1,ip
  90.          if( resdul(i) .gt. resmax ) resmax = resdul(i)
  91.   935 continue
  92.       tcrit = resmax/(machep*dfloat(n))
  93.       write(iwrite,94) tcrit,resmax
  94.    94 format(///49h as a guide to evaluating the performance of the ,
  95.      x   47h eispack  code  bandv  for the solution of band /
  96.      x   46h linear equations, a number  x, defined as the,
  97.      x   47h maximum of a set of numbers  r  related to the /
  98.      x   43h 1-norms  of the separate residual vectors, ,
  99.      x   51h may be normalized to the number  y  by dividing by /
  100.      x   46h the product of the machine precision  machep ,
  101.      x   51h and the order of the matrix.  bandv  has performed /
  102.      x   56h (well,satisfactorily,poorly) according as the ratio  y ,
  103.      x   54h is (less than  1, less than  100, greater than  100)./
  104.      x   23h for this run,  y  is =,1pd14.6,11h  with  x =,d14.6)
  105.       do 96 i = 1,ip
  106.          write(iwrite,945) i,(x(j,i),j=1,n) 
  107.   945    format(///44h the solution vector corresponding to column ,
  108.      x          i3,15h of  b  follows /(1x,1p5d24.16))
  109.          write(iwrite,95) resdul(i)
  110.    95    format(//21h for this vector r = ,1pd14.6)
  111.    96 continue
  112.       go to 10
  113.       end
  114.