home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / eispack-1.0-src.tgz / tar.out / contrib / eispack / svd.f < prev    next >
Text File  |  1996-09-28  |  10KB  |  359 lines
  1.       subroutine svd(nm,m,n,a,w,matu,u,matv,v,ierr,rv1)
  2. c
  3.       integer i,j,k,l,m,n,ii,i1,kk,k1,ll,l1,mn,nm,its,ierr
  4.       double precision a(nm,n),w(n),u(nm,n),v(nm,n),rv1(n)
  5.       double precision c,f,g,h,s,x,y,z,tst1,tst2,scale,pythag
  6.       logical matu,matv
  7. c
  8. c     this subroutine is a translation of the algol procedure svd,
  9. c     num. math. 14, 403-420(1970) by golub and reinsch.
  10. c     handbook for auto. comp., vol ii-linear algebra, 134-151(1971).
  11. c
  12. c     this subroutine determines the singular value decomposition
  13. c          t
  14. c     a=usv  of a real m by n rectangular matrix.  householder
  15. c     bidiagonalization and a variant of the qr algorithm are used.
  16. c
  17. c     on input
  18. c
  19. c        nm must be set to the row dimension of two-dimensional
  20. c          array parameters as declared in the calling program
  21. c          dimension statement.  note that nm must be at least
  22. c          as large as the maximum of m and n.
  23. c
  24. c        m is the number of rows of a (and u).
  25. c
  26. c        n is the number of columns of a (and u) and the order of v.
  27. c
  28. c        a contains the rectangular input matrix to be decomposed.
  29. c
  30. c        matu should be set to .true. if the u matrix in the
  31. c          decomposition is desired, and to .false. otherwise.
  32. c
  33. c        matv should be set to .true. if the v matrix in the
  34. c          decomposition is desired, and to .false. otherwise.
  35. c
  36. c     on output
  37. c
  38. c        a is unaltered (unless overwritten by u or v).
  39. c
  40. c        w contains the n (non-negative) singular values of a (the
  41. c          diagonal elements of s).  they are unordered.  if an
  42. c          error exit is made, the singular values should be correct
  43. c          for indices ierr+1,ierr+2,...,n.
  44. c
  45. c        u contains the matrix u (orthogonal column vectors) of the
  46. c          decomposition if matu has been set to .true.  otherwise
  47. c          u is used as a temporary array.  u may coincide with a.
  48. c          if an error exit is made, the columns of u corresponding
  49. c          to indices of correct singular values should be correct.
  50. c
  51. c        v contains the matrix v (orthogonal) of the decomposition if
  52. c          matv has been set to .true.  otherwise v is not referenced.
  53. c          v may also coincide with a if u is not needed.  if an error
  54. c          exit is made, the columns of v corresponding to indices of
  55. c          correct singular values should be correct.
  56. c
  57. c        ierr is set to
  58. c          zero       for normal return,
  59. c          k          if the k-th singular value has not been
  60. c                     determined after 30 iterations.
  61. c
  62. c        rv1 is a temporary storage array.
  63. c
  64. c     calls pythag for  dsqrt(a*a + b*b) .
  65. c
  66. c     questions and comments should be directed to burton s. garbow,
  67. c     mathematics and computer science div, argonne national laboratory
  68. c
  69. c     this version dated august 1983.
  70. c
  71. c     ------------------------------------------------------------------
  72. c
  73.       ierr = 0
  74. c
  75.       do 100 i = 1, m
  76. c
  77.          do 100 j = 1, n
  78.             u(i,j) = a(i,j)
  79.   100 continue
  80. c     .......... householder reduction to bidiagonal form ..........
  81.       g = 0.0d0
  82.       scale = 0.0d0
  83.       x = 0.0d0
  84. c
  85.       do 300 i = 1, n
  86.          l = i + 1
  87.          rv1(i) = scale * g
  88.          g = 0.0d0
  89.          s = 0.0d0
  90.          scale = 0.0d0
  91.          if (i .gt. m) go to 210
  92. c
  93.          do 120 k = i, m
  94.   120    scale = scale + dabs(u(k,i))
  95. c
  96.          if (scale .eq. 0.0d0) go to 210
  97. c
  98.          do 130 k = i, m
  99.             u(k,i) = u(k,i) / scale
  100.             s = s + u(k,i)**2
  101.   130    continue
  102. c
  103.          f = u(i,i)
  104.          g = -dsign(dsqrt(s),f)
  105.          h = f * g - s
  106.          u(i,i) = f - g
  107.          if (i .eq. n) go to 190
  108. c
  109.          do 150 j = l, n
  110.             s = 0.0d0
  111. c
  112.             do 140 k = i, m
  113.   140       s = s + u(k,i) * u(k,j)
  114. c
  115.             f = s / h
  116. c
  117.             do 150 k = i, m
  118.                u(k,j) = u(k,j) + f * u(k,i)
  119.   150    continue
  120. c
  121.   190    do 200 k = i, m
  122.   200    u(k,i) = scale * u(k,i)
  123. c
  124.   210    w(i) = scale * g
  125.          g = 0.0d0
  126.          s = 0.0d0
  127.          scale = 0.0d0
  128.          if (i .gt. m .or. i .eq. n) go to 290
  129. c
  130.          do 220 k = l, n
  131.   220    scale = scale + dabs(u(i,k))
  132. c
  133.          if (scale .eq. 0.0d0) go to 290
  134. c
  135.          do 230 k = l, n
  136.             u(i,k) = u(i,k) / scale
  137.             s = s + u(i,k)**2
  138.   230    continue
  139. c
  140.          f = u(i,l)
  141.          g = -dsign(dsqrt(s),f)
  142.          h = f * g - s
  143.          u(i,l) = f - g
  144. c
  145.          do 240 k = l, n
  146.   240    rv1(k) = u(i,k) / h
  147. c
  148.          if (i .eq. m) go to 270
  149. c
  150.          do 260 j = l, m
  151.             s = 0.0d0
  152. c
  153.             do 250 k = l, n
  154.   250       s = s + u(j,k) * u(i,k)
  155. c
  156.             do 260 k = l, n
  157.                u(j,k) = u(j,k) + s * rv1(k)
  158.   260    continue
  159. c
  160.   270    do 280 k = l, n
  161.   280    u(i,k) = scale * u(i,k)
  162. c
  163.   290    x = dmax1(x,dabs(w(i))+dabs(rv1(i)))
  164.   300 continue
  165. c     .......... accumulation of right-hand transformations ..........
  166.       if (.not. matv) go to 410
  167. c     .......... for i=n step -1 until 1 do -- ..........
  168.       do 400 ii = 1, n
  169.          i = n + 1 - ii
  170.          if (i .eq. n) go to 390
  171.          if (g .eq. 0.0d0) go to 360
  172. c
  173.          do 320 j = l, n
  174. c     .......... double division avoids possible underflow ..........
  175.   320    v(j,i) = (u(i,j) / u(i,l)) / g
  176. c
  177.          do 350 j = l, n
  178.             s = 0.0d0
  179. c
  180.             do 340 k = l, n
  181.   340       s = s + u(i,k) * v(k,j)
  182. c
  183.             do 350 k = l, n
  184.                v(k,j) = v(k,j) + s * v(k,i)
  185.   350    continue
  186. c
  187.   360    do 380 j = l, n
  188.             v(i,j) = 0.0d0
  189.             v(j,i) = 0.0d0
  190.   380    continue
  191. c
  192.   390    v(i,i) = 1.0d0
  193.          g = rv1(i)
  194.          l = i
  195.   400 continue
  196. c     .......... accumulation of left-hand transformations ..........
  197.   410 if (.not. matu) go to 510
  198. c     ..........for i=min(m,n) step -1 until 1 do -- ..........
  199.       mn = n
  200.       if (m .lt. n) mn = m
  201. c
  202.       do 500 ii = 1, mn
  203.          i = mn + 1 - ii
  204.          l = i + 1
  205.          g = w(i)
  206.          if (i .eq. n) go to 430
  207. c
  208.          do 420 j = l, n
  209.   420    u(i,j) = 0.0d0
  210. c
  211.   430    if (g .eq. 0.0d0) go to 475
  212.          if (i .eq. mn) go to 460
  213. c
  214.          do 450 j = l, n
  215.             s = 0.0d0
  216. c
  217.             do 440 k = l, m
  218.   440       s = s + u(k,i) * u(k,j)
  219. c     .......... double division avoids possible underflow ..........
  220.             f = (s / u(i,i)) / g
  221. c
  222.             do 450 k = i, m
  223.                u(k,j) = u(k,j) + f * u(k,i)
  224.   450    continue
  225. c
  226.   460    do 470 j = i, m
  227.   470    u(j,i) = u(j,i) / g
  228. c
  229.          go to 490
  230. c
  231.   475    do 480 j = i, m
  232.   480    u(j,i) = 0.0d0
  233. c
  234.   490    u(i,i) = u(i,i) + 1.0d0
  235.   500 continue
  236. c     .......... diagonalization of the bidiagonal form ..........
  237.   510 tst1 = x
  238. c     .......... for k=n step -1 until 1 do -- ..........
  239.       do 700 kk = 1, n
  240.          k1 = n - kk
  241.          k = k1 + 1
  242.          its = 0
  243. c     .......... test for splitting.
  244. c                for l=k step -1 until 1 do -- ..........
  245.   520    do 530 ll = 1, k
  246.             l1 = k - ll
  247.             l = l1 + 1
  248.             tst2 = tst1 + dabs(rv1(l))
  249.             if (tst2 .eq. tst1) go to 565
  250. c     .......... rv1(1) is always zero, so there is no exit
  251. c                through the bottom of the loop ..........
  252.             tst2 = tst1 + dabs(w(l1))
  253.             if (tst2 .eq. tst1) go to 540
  254.   530    continue
  255. c     .......... cancellation of rv1(l) if l greater than 1 ..........
  256.   540    c = 0.0d0
  257.          s = 1.0d0
  258. c
  259.          do 560 i = l, k
  260.             f = s * rv1(i)
  261.             rv1(i) = c * rv1(i)
  262.             tst2 = tst1 + dabs(f)
  263.             if (tst2 .eq. tst1) go to 565
  264.             g = w(i)
  265.             h = pythag(f,g)
  266.             w(i) = h
  267.             c = g / h
  268.             s = -f / h
  269.             if (.not. matu) go to 560
  270. c
  271.             do 550 j = 1, m
  272.                y = u(j,l1)
  273.                z = u(j,i)
  274.                u(j,l1) = y * c + z * s
  275.                u(j,i) = -y * s + z * c
  276.   550       continue
  277. c
  278.   560    continue
  279. c     .......... test for convergence ..........
  280.   565    z = w(k)
  281.          if (l .eq. k) go to 650
  282. c     .......... shift from bottom 2 by 2 minor ..........
  283.          if (its .eq. 30) go to 1000
  284.          its = its + 1
  285.          x = w(l)
  286.          y = w(k1)
  287.          g = rv1(k1)
  288.          h = rv1(k)
  289.          f = 0.5d0 * (((g + z) / h) * ((g - z) / y) + y / h - h / y)
  290.          g = pythag(f,1.0d0)
  291.          f = x - (z / x) * z + (h / x) * (y / (f + dsign(g,f)) - h)
  292. c     .......... next qr transformation ..........
  293.          c = 1.0d0
  294.          s = 1.0d0
  295. c
  296.          do 600 i1 = l, k1
  297.             i = i1 + 1
  298.             g = rv1(i)
  299.             y = w(i)
  300.             h = s * g
  301.             g = c * g
  302.             z = pythag(f,h)
  303.             rv1(i1) = z
  304.             c = f / z
  305.             s = h / z
  306.             f = x * c + g * s
  307.             g = -x * s + g * c
  308.             h = y * s
  309.             y = y * c
  310.             if (.not. matv) go to 575
  311. c
  312.             do 570 j = 1, n
  313.                x = v(j,i1)
  314.                z = v(j,i)
  315.                v(j,i1) = x * c + z * s
  316.                v(j,i) = -x * s + z * c
  317.   570       continue
  318. c
  319.   575       z = pythag(f,h)
  320.             w(i1) = z
  321. c     .......... rotation can be arbitrary if z is zero ..........
  322.             if (z .eq. 0.0d0) go to 580
  323.             c = f / z
  324.             s = h / z
  325.   580       f = c * g + s * y
  326.             x = -s * g + c * y
  327.             if (.not. matu) go to 600
  328. c
  329.             do 590 j = 1, m
  330.                y = u(j,i1)
  331.                z = u(j,i)
  332.                u(j,i1) = y * c + z * s
  333.                u(j,i) = -y * s + z * c
  334.   590       continue
  335. c
  336.   600    continue
  337. c
  338.          rv1(l) = 0.0d0
  339.          rv1(k) = f
  340.          w(k) = x
  341.          go to 520
  342. c     .......... convergence ..........
  343.   650    if (z .ge. 0.0d0) go to 700
  344. c     .......... w(k) is made non-negative ..........
  345.          w(k) = -z
  346.          if (.not. matv) go to 700
  347. c
  348.          do 690 j = 1, n
  349.   690    v(j,k) = -v(j,k)
  350. c
  351.   700 continue
  352. c
  353.       go to 1001
  354. c     .......... set error -- no convergence to a
  355. c                singular value after 30 iterations ..........
  356.  1000 ierr = k
  357.  1001 return
  358.       end
  359.