@process directive('" (') subroutine svd(nm,m,n,a,w,matu,u,matv,v,ierr,rv1) c integer i,j,k,l,m,n,i1,k1,l1,mn,nm,its,ierr real a(nm,n),w(n),u(nm,n),v(nm,n),rv1(n) real c,f,g,h,s,x,y,z,tst1,tst2,scale logical matu,matv c c this subroutine is a translation of the algol procedure svd, 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 the singular value decomposition c t c a=usv of a real m by n rectangular matrix. 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 u). c c n is the number of columns of a (and u) and the order of v. c c a contains the rectangular input matrix to be decomposed. c c matu should be set to .true. if the u matrix in the c decomposition is desired, and to .false. otherwise. c c matv should be set to .true. if the v matrix in the c decomposition is desired, and to .false. otherwise. c c on output c c a is unaltered (unless overwritten by u or v). 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 u contains the matrix u (orthogonal column vectors) of the c decomposition if matu has been set to .true. otherwise c u is used as a temporary array. u may coincide with a. c if an error exit is made, the columns of u corresponding c to indices of correct singular values should be correct. c c v contains the matrix v (orthogonal) of the decomposition if c matv has been set to .true. otherwise v is not referenced. c v may also coincide with a if u is not needed. if an error c exit is made, the columns of v corresponding to indices of c correct 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 Questions and comments should be directed to Alan K. Cline, c Pleasant Valley Software, 8603 Altus Cove, Austin, TX 78759. c Electronic mail to cline@cs.utexas.edu. c c this version dated january 1989. (for the IBM 3090vf) c c ------------------------------------------------------------------ c call xuflow(0) ierr = 0 c do 100 i = 1, m c do 100 j = 1, n u(i,j) = a(i,j) 100 continue c .......... householder reduction to bidiagonal form .......... g = 0.0e0 scale = 0.0e0 x = 0.0e0 c do 300 i = 1, n l = i + 1 rv1(i) = scale * g g = 0.0e0 s = 0.0e0 scale = 0.0e0 if (i .gt. m) go to 210 c do 120 k = i, m 120 scale = scale + abs(u(k,i)) c if (scale .eq. 0.0e0) go to 210 c do 130 k = i, m u(k,i) = u(k,i) / scale s = s + u(k,i)**2 130 continue c f = u(i,i) g = -sign(sqrt(s),f) h = f * g - s u(i,i) = f - g c c" ( ignore recrdeps do 150 j = l, n s = 0.0e0 c do 140 k = i, m 140 s = s + u(k,i) * u(k,j) c f = s / h c do 150 k = i, m u(k,j) = u(k,j) + f * u(k,i) 150 continue c do 200 k = i, m 200 u(k,i) = scale * u(k,i) c 210 w(i) = scale * g g = 0.0e0 s = 0.0e0 scale = 0.0e0 c if (i .gt. m .or. i .eq. n) go to 290 c" ( prefer vector do 220 k = l, n 220 scale = scale + abs(u(i,k)) c if (scale .eq. 0.0e0) go to 290 c do 230 k = l, n u(i,k) = u(i,k) / scale s = s + u(i,k)**2 230 continue c f = u(i,l) g = -sign(sqrt(s),f) h = f * g - s u(i,l) = f - g c do 240 k = l, n 240 rv1(k) = u(i,k) / h c c" ( ignore recrdeps do 260 j = l, m s = 0.0e0 c do 250 k = l, n 250 s = s + u(j,k) * u(i,k) c do 260 k = l, n u(j,k) = u(j,k) + s * rv1(k) 260 continue c c" ( prefer vector do 280 k = l, n 280 u(i,k) = scale * u(i,k) c 290 x = amax1(x,abs(w(i))+abs(rv1(i))) 300 continue c .......... accumulation of right-hand transformations .......... if (.not. matv) go to 410 c .......... for i=n step -1 until 1 do -- .......... do 400 i = n, 1, -1 if (i .eq. n) go to 390 if (g .eq. 0.0e0) go to 360 c do 320 j = l, n c .......... double division avoids possible underflow .......... 320 v(j,i) = (u(i,j) / u(i,l)) / g c c" ( ignore recrdeps do 350 j = l, n s = 0.0e0 c do 340 k = l, n 340 s = s + u(i,k) * v(k,j) c do 350 k = l, n v(k,j) = v(k,j) + s * v(k,i) 350 continue c 360 do 380 j = l, n v(i,j) = 0.0e0 v(j,i) = 0.0e0 380 continue c 390 v(i,i) = 1.0e0 g = rv1(i) l = i 400 continue c .......... accumulation of left-hand transformations .......... 410 if (.not. matu) go to 510 c ..........for i=min(m,n) step -1 until 1 do -- .......... mn = n if (m .lt. n) mn = m c do 500 i = mn, 1, -1 l = i + 1 g = w(i) if (i .eq. n) go to 430 c c" ( prefer vector do 420 j = l, n 420 u(i,j) = 0.0e0 c 430 if (g .eq. 0.0e0) go to 475 if (i .eq. mn) go to 460 c c" ( ignore recrdeps do 450 j = l, n s = 0.0e0 c do 440 k = l, m 440 s = s + u(k,i) * u(k,j) c .......... double division avoids possible underflow .......... f = (s / u(i,i)) / g c do 450 k = i, m u(k,j) = u(k,j) + f * u(k,i) 450 continue c 460 do 470 j = i, m 470 u(j,i) = u(j,i) / g c go to 490 c 475 do 480 j = i, m 480 u(j,i) = 0.0e0 c 490 u(i,i) = u(i,i) + 1.0e0 500 continue c .......... diagonalization of the bidiagonal form .......... 510 tst1 = x c .......... for k=n step -1 until 1 do -- .......... do 700 k = n, 1, -1 k1 = k - 1 its = 0 c .......... test for splitting. c for l=k step -1 until 1 do -- .......... 520 do 530 l = k, 1, -1 l1 = l - 1 tst2 = tst1 + abs(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 + abs(w(l1)) if (tst2 .eq. tst1) go to 540 530 continue c .......... cancellation of rv1(l) if l greater than 1 .......... 540 c = 0.0e0 s = 1.0e0 c do 560 i = l, k f = s * rv1(i) rv1(i) = c * rv1(i) tst2 = tst1 + abs(f) if (tst2 .eq. tst1) go to 565 g = w(i) c-- h = pythag(f,g) if (abs(f).gt.abs(g)) then h = abs(f)*sqrt(1e0+(g/f)**2) else if (g.ne.0e0) then h = abs(g)*sqrt((f/g)**2+1e0) else h = abs(f) endif w(i) = h c = g / h s = -f / h if (.not. matu) go to 560 c do 550 j = 1, m y = u(j,l1) z = u(j,i) u(j,l1) = y * c + z * s u(j,i) = -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.5e0 * (((g + z) / h) * ((g - z) / y) + y / h - h / y) c** g = pythag(f,1.0d0) if (abs(f).gt.1.0e0) then g = abs(f)*sqrt(1e0+(1.0e0/f)**2) else g = sqrt(f**2+1e0) endif f = x - (z / x) * z + (h / x) * (y / (f + sign(g,f)) - h) c .......... next qr transformation .......... c = 1.0e0 s = 1.0e0 c do 600 i1 = l, k1 i = i1 + 1 g = rv1(i) y = w(i) h = s * g g = c * g c-- z = pythag(f,h) if (abs(f).gt.abs(h)) then z = abs(f)*sqrt(1e0+(h/f)**2) else if (h.ne.0e0) then z = abs(h)*sqrt((f/h)**2+1e0) else z = abs(f) endif 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 if (.not. matv) go to 575 c do 570 j = 1, n x = v(j,i1) z = v(j,i) v(j,i1) = x * c + z * s v(j,i) = -x * s + z * c 570 continue c 575 continue c-- z = pythag(f,h) if (abs(f).gt.abs(h)) then z = abs(f)*sqrt(1e0+(h/f)**2) else if (h.ne.0e0) then z = abs(h)*sqrt((f/h)**2+1e0) else z = abs(f) endif w(i1) = z c .......... rotation can be arbitrary if z is zero .......... if (z .eq. 0.0e0) go to 580 c = f / z s = h / z 580 f = c * g + s * y x = -s * g + c * y if (.not. matu) go to 600 c do 590 j = 1, m y = u(j,i1) z = u(j,i) u(j,i1) = y * c + z * s u(j,i) = -y * s + z * c 590 continue c 600 continue c rv1(l) = 0.0e0 rv1(k) = f w(k) = x go to 520 c .......... convergence .......... 650 if (z .ge. 0.0e0) go to 700 c .......... w(k) is made non-negative .......... w(k) = -z if (.not. matv) go to 700 c do 690 j = 1, n 690 v(j,k) = -v(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