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 / rleispak.f < prev    next >
Text File  |  1996-09-28  |  4KB  |  117 lines
  1. c
  2. c     this driver tests  eispack  for the class of real matrices
  3. c     exhibiting the use of  eispack  to find the singular values
  4. c     and the solution to the equation  a*x = b .
  5. c
  6. c     this driver is catalogued as  eispdrv4(rleispak).
  7. c
  8. c     the dimension of  a,u,  and  v  should be  nm  by  nm.
  9. c     the dimension of  sigma  and  rv1  should be  nm.
  10. c     the dimension of  ahold  should be  nm  by  nm.
  11. c     here  nm = 20.
  12. c
  13.       double precision a( 20, 20),u( 20, 20),v( 20, 20),
  14.      x        sigma( 20),rv1( 20),ahold( 20, 20),machep,
  15.      x        resdul,tcrit
  16.       integer error
  17.       data iwrite/6/
  18. c
  19. c     .......... machep is a machine dependent parameter specifying
  20. c                the relative precision of floating point arithmetic.
  21.       machep = 1.0d0
  22.     1 machep = 0.5d0*machep
  23.       if (1.0d0 + machep .gt. 1.0d0 ) go to 1
  24.       machep = 2.0d0*machep
  25. c
  26. c
  27.       nm = 20
  28.    10 call  rmatin(nm,nm,m,n,a,ahold,ahold,0,0)
  29.       ip = m
  30.       write(iwrite,20) m,n
  31.    20 format(30h1the matrix  a , row dimension  ,i4,
  32.      x         22h  and column dimension  ,i4,
  33.      x         18h, printed by rows.)
  34.       do 40 i = 1,m
  35.          write(iwrite,30) (a(i,j),j=1,n) 
  36.    30    format(/(1x,1p5d24.16))
  37.    40 continue
  38.       write(iwrite,50) m,ip
  39.    50 format(///30h the matrix  b , row dimension  ,i4,
  40.      x       22h  and column dimension  ,i4,
  41.      x       25h, is the identity matrix.)
  42. c
  43.       do 230 icall = 1,2
  44.          if( icall .eq. 1 ) go to 90
  45.          call  rmatin(nm,nm,m,n,a,ahold,ahold,0,1)
  46.          go to 100
  47. c
  48. c     rl  using  svd
  49. c
  50.    90    write(iwrite,91)
  51.    91    format(//43h1all the singular values of the real matrix,
  52.      x          44h follow.  the path involving  svd  was used. )
  53.          call  svd(nm,m,n,a,sigma,.true.,u,.true.,v,error,rv1)
  54.          write(iwrite,92) error
  55.    92    format(//28h *****error from svd***** = ,i4)
  56.          go to 200
  57. c
  58. c     rl  using  minfit
  59. c
  60.   100    if( ip .eq. 0 ) go to 230
  61.          write(iwrite,101)
  62.   101    format(44h1all the singular values of the real matrix ,
  63.      x          46hfollow.  the path involving  minfit  was used.)
  64.          do 103 i = 1,m
  65.             do 102 j = 1,ip
  66.                v(i,j) = 0.0d0
  67.   102       continue
  68.             v(i,i) = 1.0d0
  69.   103    continue
  70.          call  minfit(nm,m,n,a,sigma,ip,v,error,rv1)
  71.          write(iwrite,104) error
  72.   104    format(//31h *****error from minfit***** = ,i4)
  73.          do 107 i = 1,m
  74.             do 106 j = 1,n
  75.                u(i,j) = v(j,i)
  76.   106       continue
  77.   107    continue
  78.          do 109 i = 1,n
  79.             do 108 j = 1,n
  80.                v(i,j) = a(i,j)
  81.   108       continue
  82.   109    continue
  83.   200    write(iwrite,201) (sigma(i),i=1,n)
  84.   201    format(//26h the singular values of a /3(1pd24.16,5x))
  85.          call  rmatin(nm,nm,m,n,a,ahold,ahold,0,1)
  86.          call  rlres(nm,m,n,a,u,sigma,v,rv1,resdul)
  87.          tcrit = resdul/(machep*dfloat(10*max0(m,n)))
  88.          write(iwrite,202) tcrit,resdul
  89.   202    format(///49h as a guide to evaluating the performance of the ,
  90.      x      47h eispack  codes  svd  and  minfit, a number  x, /
  91.      x      47h related to the  1-norm  of the residual matrix,
  92.      x      44h for the singular value decomposition of  a, /
  93.      x      48h may be normalized to the number  y  by dividing,
  94.      x      48h by the product of the machine precision  machep /
  95.      x      41h and the larger of the matrix dimensions.,
  96.      x      36h  the  eispack  codes have performed /
  97.      x      56h (well,satisfactorily,poorly) according as the ratio  y ,
  98.      x      54h is (less than  1, less than  100, greater than  100)./
  99.      x      23h for this run,  y  is =,1pd14.6,11h  with  x =,d14.6)
  100.          write(iwrite,203) m,n
  101.   203    format(///30h the matrix  u , row dimension ,i4,
  102.      x          21h and column dimension ,i4,
  103.      x          10h, follows. )
  104.          do 205 i = 1,m
  105.             write(iwrite,204) (u(i,j),j=1,n)
  106.   204       format(/(1x,1p5d24.16))
  107.   205    continue
  108.          write(iwrite,206) n
  109.   206    format(///25h the matrix  v , of order ,i4,
  110.      x          10h, follows. )
  111.          do 207 i = 1,n
  112.             write(iwrite,204) (v(i,j),j=1,n)  
  113.   207    continue
  114.   230 continue
  115.       go to 10
  116.       end
  117.