subroutine minfit(nm,m,n,a,w,ip,b,ierr,rv1) c integer i,j,k,l,m,n,ii,ip,i1,kk,k1,ll,l1,m1,nm,its,ierr double precision a(nm,n),w(n),b(nm,ip),rv1(n) double precision c,f,g,h,s,x,y,z,tst1,tst2,scale,pythag c c this subroutine is a translation of the algol procedure minfit, c num. math. 14, 403-420(1970) by golub and reinsch. c handbook for auto. comp., vol ii-linear algebra, 134-151(1971). c c this subroutine determines, towards the solution of the linear c t c system ax=b, the singular value decomposition a=usv of a real c t c m by n rectangular matrix, forming u b rather than u. householder c bidiagonalization and a variant of the qr algorithm are used. c c on input c c nm must be set to the row dimension of two-dimensional c array parameters as declared in the calling program c dimension statement. note that nm must be at least c as large as the maximum of m and n. c c m is the number of rows of a and b. c c n is the number of columns of a and the order of v. c c a contains the rectangular coefficient matrix of the system. c c ip is the number of columns of b. ip can be zero. c c b contains the constant column matrix of the system c if ip is not zero. otherwise b is not referenced. c c on output c c a has been overwritten by the matrix v (orthogonal) of the c decomposition in its first n rows and columns. if an c error exit is made, the columns of v corresponding to c indices of correct singular values should be correct. c c w contains the n (non-negative) singular values of a (the c diagonal elements of s). they are unordered. if an c error exit is made, the singular values should be correct c for indices ierr+1,ierr+2,...,n. c c t c b has been overwritten by u b. if an error exit is made, c t c the rows of u b corresponding to indices of correct c singular values should be correct. c c ierr is set to c zero for normal return, c k if the k-th singular value has not been c determined after 30 iterations. c c rv1 is a temporary storage array. c c calls pythag for dsqrt(a*a + b*b) . c c questions and comments should be directed to burton s. garbow, c mathematics and computer science div, argonne national laboratory c c this version dated august 1983. c c ------------------------------------------------------------------ c ierr = 0 c .......... householder reduction to bidiagonal form .......... g = 0.0d0 scale = 0.0d0 x = 0.0d0 c do 300 i = 1, n l = i + 1 rv1(i) = scale * g g = 0.0d0 s = 0.0d0 scale = 0.0d0 if (i .gt. m) go to 210 c do 120 k = i, m 120 scale = scale + dabs(a(k,i)) c if (scale .eq. 0.0d0) go to 210 c do 130 k = i, m a(k,i) = a(k,i) / scale s = s + a(k,i)**2 130 continue c f = a(i,i) g = -dsign(dsqrt(s),f) h = f * g - s a(i,i) = f - g if (i .eq. n) go to 160 c do 150 j = l, n s = 0.0d0 c do 140 k = i, m 140 s = s + a(k,i) * a(k,j) c f = s / h c do 150 k = i, m a(k,j) = a(k,j) + f * a(k,i) 150 continue c 160 if (ip .eq. 0) go to 190 c do 180 j = 1, ip s = 0.0d0 c do 170 k = i, m 170 s = s + a(k,i) * b(k,j) c f = s / h c do 180 k = i, m b(k,j) = b(k,j) + f * a(k,i) 180 continue c 190 do 200 k = i, m 200 a(k,i) = scale * a(k,i) c 210 w(i) = scale * g g = 0.0d0 s = 0.0d0 scale = 0.0d0 if (i .gt. m .or. i .eq. n) go to 290 c do 220 k = l, n 220 scale = scale + dabs(a(i,k)) c if (scale .eq. 0.0d0) go to 290 c do 230 k = l, n a(i,k) = a(i,k) / scale s = s + a(i,k)**2 230 continue c f = a(i,l) g = -dsign(dsqrt(s),f) h = f * g - s a(i,l) = f - g c do 240 k = l, n 240 rv1(k) = a(i,k) / h c if (i .eq. m) go to 270 c do 260 j = l, m s = 0.0d0 c do 250 k = l, n 250 s = s + a(j,k) * a(i,k) c do 260 k = l, n a(j,k) = a(j,k) + s * rv1(k) 260 continue c 270 do 280 k = l, n 280 a(i,k) = scale * a(i,k) c 290 x = dmax1(x,dabs(w(i))+dabs(rv1(i))) 300 continue c .......... accumulation of right-hand transformations. c for i=n step -1 until 1 do -- .......... do 400 ii = 1, n i = n + 1 - ii if (i .eq. n) go to 390 if (g .eq. 0.0d0) go to 360 c do 320 j = l, n c .......... double division avoids possible underflow .......... 320 a(j,i) = (a(i,j) / a(i,l)) / g c do 350 j = l, n s = 0.0d0 c do 340 k = l, n 340 s = s + a(i,k) * a(k,j) c do 350 k = l, n a(k,j) = a(k,j) + s * a(k,i) 350 continue c 360 do 380 j = l, n a(i,j) = 0.0d0 a(j,i) = 0.0d0 380 continue c 390 a(i,i) = 1.0d0 g = rv1(i) l = i 400 continue c if (m .ge. n .or. ip .eq. 0) go to 510 m1 = m + 1 c do 500 i = m1, n c do 500 j = 1, ip b(i,j) = 0.0d0 500 continue c .......... diagonalization of the bidiagonal form .......... 510 tst1 = x c .......... for k=n step -1 until 1 do -- .......... do 700 kk = 1, n k1 = n - kk k = k1 + 1 its = 0 c .......... test for splitting. c for l=k step -1 until 1 do -- .......... 520 do 530 ll = 1, k l1 = k - ll l = l1 + 1 tst2 = tst1 + dabs(rv1(l)) if (tst2 .eq. tst1) go to 565 c .......... rv1(1) is always zero, so there is no exit c through the bottom of the loop .......... tst2 = tst1 + dabs(w(l1)) if (tst2 .eq. tst1) go to 540 530 continue c .......... cancellation of rv1(l) if l greater than 1 .......... 540 c = 0.0d0 s = 1.0d0 c do 560 i = l, k f = s * rv1(i) rv1(i) = c * rv1(i) tst2 = tst1 + dabs(f) if (tst2 .eq. tst1) go to 565 g = w(i) h = pythag(f,g) w(i) = h c = g / h s = -f / h if (ip .eq. 0) go to 560 c do 550 j = 1, ip y = b(l1,j) z = b(i,j) b(l1,j) = y * c + z * s b(i,j) = -y * s + z * c 550 continue c 560 continue c .......... test for convergence .......... 565 z = w(k) if (l .eq. k) go to 650 c .......... shift from bottom 2 by 2 minor .......... if (its .eq. 30) go to 1000 its = its + 1 x = w(l) y = w(k1) g = rv1(k1) h = rv1(k) f = 0.5d0 * (((g + z) / h) * ((g - z) / y) + y / h - h / y) g = pythag(f,1.0d0) f = x - (z / x) * z + (h / x) * (y / (f + dsign(g,f)) - h) c .......... next qr transformation .......... c = 1.0d0 s = 1.0d0 c do 600 i1 = l, k1 i = i1 + 1 g = rv1(i) y = w(i) h = s * g g = c * g z = pythag(f,h) rv1(i1) = z c = f / z s = h / z f = x * c + g * s g = -x * s + g * c h = y * s y = y * c c do 570 j = 1, n x = a(j,i1) z = a(j,i) a(j,i1) = x * c + z * s a(j,i) = -x * s + z * c 570 continue c z = pythag(f,h) w(i1) = z c .......... rotation can be arbitrary if z is zero .......... if (z .eq. 0.0d0) go to 580 c = f / z s = h / z 580 f = c * g + s * y x = -s * g + c * y if (ip .eq. 0) go to 600 c do 590 j = 1, ip y = b(i1,j) z = b(i,j) b(i1,j) = y * c + z * s b(i,j) = -y * s + z * c 590 continue c 600 continue c rv1(l) = 0.0d0 rv1(k) = f w(k) = x go to 520 c .......... convergence .......... 650 if (z .ge. 0.0d0) go to 700 c .......... w(k) is made non-negative .......... w(k) = -z c do 690 j = 1, n 690 a(j,k) = -a(j,k) c 700 continue c go to 1001 c .......... set error -- no convergence to a c singular value after 30 iterations .......... 1000 ierr = k 1001 return end