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 / rlres.f < prev    next >
Text File  |  1996-09-28  |  2KB  |  91 lines
  1.       subroutine rlres(nm,m,n,a,u,s,v,norm,resdul)
  2. c
  3.       double precision a(nm,n),u(nm,n),v(nm,n),s(n),norm(n)
  4.       double precision sum1,sum2,x,resdul,norma,suma,sumuv
  5. c
  6. c     this subroutine performs a residual test for the decomposition
  7. c                              t
  8. c     of a matrix  a  into  usv .
  9. c
  10. c     this subroutine is catalogued as eispdrv4(rlres).
  11. c
  12. c     input.
  13. c
  14. c       nm  is the row dimension of two-dimensional array parameters
  15. c         as declared in the calling program dimension statement;
  16. c
  17. c       m  is the row dimension of the matrices  a  and  u.
  18. c
  19. c       n  is the column dimension of the matrices  a  and u , and
  20. c         the order of  v;
  21. c
  22. c       a  is a real general matrix of dimemsion  m  by  n;
  23. c
  24. c       u  is a matrix with orthogonal columns of dimension  m  by  n;
  25. c
  26. c       s  is a diagonal matrix stored as a vector of dimension  n;
  27. c
  28. c       v  is an orthogonal matrix of dimension  n  by  n.
  29. c
  30. c     output.
  31. c
  32. c       norm  is an array such that for each k,
  33. c                                              t
  34. c         norm(k) = !!a*v(k)-s(k)*u(k)!!  + !!a *u(k)-s(k)*v(k)!!
  35. c                   ---------------------------------------------
  36. c                           !!a!!*(!!u(k)!! + !!v(k)!!)
  37. c
  38. c       resdul  is the largest element of  norm.
  39. c
  40. c     ------------------------------------------------------------------
  41. c
  42.       resdul = 0.0d0
  43.       norma = 0.0d0
  44. c
  45.       do 20 i = 1,m
  46.          suma = 0.0d0
  47. c
  48.          do 10 j = 1,n
  49.             suma = suma + dabs(a(i,j))
  50.    10    continue
  51. c
  52.          norma = dmax1(norma,suma)
  53.    20 continue
  54. c
  55.       if( norma .eq. 0.0d0 ) norma = 1.0d0
  56. c
  57.       do 70 k = 1,n
  58.          sumuv = 0.0d0
  59.          sum1 = 0.0d0
  60. c
  61.          do 40 i = 1,m
  62.             sumuv = sumuv + dabs(u(i,k))
  63.             x = -s(k)*u(i,k)
  64. c
  65.             do 30 j = 1,n
  66.                x = x + a(i,j)*v(j,k)
  67.    30       continue
  68. c
  69.             sum1 = sum1 + dabs(x)
  70.    40    continue
  71. c
  72.          sum2 = 0.0d0
  73. c
  74.          do 60 i = 1,n
  75.             sumuv = sumuv + dabs(v(i,k))
  76.             x = -s(k)*v(i,k)
  77. c
  78.             do 50 j = 1,m
  79.                x = x + a(j,i)*u(j,k)
  80.    50       continue
  81. c
  82.             sum2 = sum2 + dabs(x)
  83.    60    continue
  84. c
  85.          norm(k) = (sum1 + sum2)/(sumuv*norma)
  86.          resdul = dmax1(norm(k),resdul)
  87.    70 continue
  88. c
  89.       return
  90.       end
  91.