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