1 subroutine dsvdc(x,ldx,n,p,s,e,u,ldu,v,ldv,work,job,info) 2 integer ldx,n,p,ldu,ldv,job,info 3 double precision x(ldx,1),s(1),e(1),u(ldu,1),v(ldv,1),work(1) 4c 5c 6c dsvdc is a subroutine to reduce a double precision nxp matrix x 7c by orthogonal transformations u and v to diagonal form. the 8c diagonal elements s(i) are the singular values of x. the 9c columns of u are the corresponding left singular vectors, 10c and the columns of v the right singular vectors. 11c 12c on entry 13c 14c x double precision(ldx,p), where ldx.ge.n. 15c x contains the matrix whose singular value 16c decomposition is to be computed. x is 17c destroyed by dsvdc. 18c 19c ldx integer. 20c ldx is the leading dimension of the array x. 21c 22c n integer. 23c n is the number of rows of the matrix x. 24c 25c p integer. 26c p is the number of columns of the matrix x. 27c 28c ldu integer. 29c ldu is the leading dimension of the array u. 30c (see below). 31c 32c ldv integer. 33c ldv is the leading dimension of the array v. 34c (see below). 35c 36c work double precision(n). 37c work is a scratch array. 38c 39c job integer. 40c job controls the computation of the singular 41c vectors. it has the decimal expansion ab 42c with the following meaning 43c 44c a.eq.0 do not compute the left singular 45c vectors. 46c a.eq.1 return the n left singular vectors 47c in u. 48c a.ge.2 return the first min(n,p) singular 49c vectors in u. 50c b.eq.0 do not compute the right singular 51c vectors. 52c b.eq.1 return the right singular vectors 53c in v. 54c 55c on return 56c 57c s double precision(mm), where mm=min(n+1,p). 58c the first min(n,p) entries of s contain the 59c singular values of x arranged in descending 60c order of magnitude. 61c 62c e double precision(p), 63c e ordinarily contains zeros. however see the 64c discussion of info for exceptions. 65c 66c u double precision(ldu,k), where ldu.ge.n. if 67c joba.eq.1 then k.eq.n, if joba.ge.2 68c then k.eq.min(n,p). 69c u contains the matrix of left singular vectors. 70c u is not referenced if joba.eq.0. if n.le.p 71c or if joba.eq.2, then u may be identified with x 72c in the subroutine call. 73c 74c v double precision(ldv,p), where ldv.ge.p. 75c v contains the matrix of right singular vectors. 76c v is not referenced if job.eq.0. if p.le.n, 77c then v may be identified with x in the 78c subroutine call. 79c 80c info integer. 81c the singular values (and their corresponding 82c singular vectors) s(info+1),s(info+2),...,s(m) 83c are correct (here m=min(n,p)). thus if 84c info.eq.0, all the singular values and their 85c vectors are correct. in any event, the matrix 86c b = trans(u)*x*v is the bidiagonal matrix 87c with the elements of s on its diagonal and the 88c elements of e on its super-diagonal (trans(u) 89c is the transpose of u). thus the singular 90c values of x and b are the same. 91c 92c linpack. this version dated 08/14/78 . 93c correction made to shift 2/84. 94c g.w. stewart, university of maryland, argonne national lab. 95c 96c dsvdc uses the following functions and subprograms. 97c 98c external drot 99c blas daxpy,ddot,dscal,dswap,dnrm2,drotg 100c fortran dabs,dmax1,max0,min0,mod,dsqrt 101c 102c internal variables 103c 104 integer i,iter,j,jobu,k,kase,kk,l,ll,lls,lm1,lp1,ls,lu,m,maxit, 105 * mm,mm1,mp1,nct,nctp1,ncu,nrt,nrtp1 106 double precision ddot,t,r 107 double precision b,c,cs,el,emm1,f,g,dnrm2,scale,shift,sl,sm,sn, 108 * smm1,t1,test,ztest 109 logical wantu,wantv 110c 111c 112c set the maximum number of iterations. 113c 114 maxit = 1000 115c 116c determine what is to be computed. 117c 118 wantu = .false. 119 wantv = .false. 120 jobu = mod(job,100)/10 121 ncu = n 122 if (jobu .gt. 1) ncu = min0(n,p) 123 if (jobu .ne. 0) wantu = .true. 124 if (mod(job,10) .ne. 0) wantv = .true. 125c 126c reduce x to bidiagonal form, storing the diagonal elements 127c in s and the super-diagonal elements in e. 128c 129 info = 0 130 nct = min0(n-1,p) 131 nrt = max0(0,min0(p-2,n)) 132 lu = max0(nct,nrt) 133 if (lu .lt. 1) go to 170 134 do 160 l = 1, lu 135 lp1 = l + 1 136 if (l .gt. nct) go to 20 137c 138c compute the transformation for the l-th column and 139c place the l-th diagonal in s(l). 140c 141 s(l) = dnrm2(n-l+1,x(l,l),1) 142 if (s(l) .eq. 0.0d0) go to 10 143 if (x(l,l) .ne. 0.0d0) s(l) = dsign(s(l),x(l,l)) 144 call dscal(n-l+1,1.0d0/s(l),x(l,l),1) 145 x(l,l) = 1.0d0 + x(l,l) 146 10 continue 147 s(l) = -s(l) 148 20 continue 149 if (p .lt. lp1) go to 50 150 do 40 j = lp1, p 151 if (l .gt. nct) go to 30 152 if (s(l) .eq. 0.0d0) go to 30 153c 154c apply the transformation. 155c 156 t = -ddot(n-l+1,x(l,l),1,x(l,j),1)/x(l,l) 157 call daxpy(n-l+1,t,x(l,l),1,x(l,j),1) 158 30 continue 159c 160c place the l-th row of x into e for the 161c subsequent calculation of the row transformation. 162c 163 e(j) = x(l,j) 164 40 continue 165 50 continue 166 if (.not.wantu .or. l .gt. nct) go to 70 167c 168c place the transformation in u for subsequent back 169c multiplication. 170c 171 do 60 i = l, n 172 u(i,l) = x(i,l) 173 60 continue 174 70 continue 175 if (l .gt. nrt) go to 150 176c 177c compute the l-th row transformation and place the 178c l-th super-diagonal in e(l). 179c 180 e(l) = dnrm2(p-l,e(lp1),1) 181 if (e(l) .eq. 0.0d0) go to 80 182 if (e(lp1) .ne. 0.0d0) e(l) = dsign(e(l),e(lp1)) 183 call dscal(p-l,1.0d0/e(l),e(lp1),1) 184 e(lp1) = 1.0d0 + e(lp1) 185 80 continue 186 e(l) = -e(l) 187 if (lp1 .gt. n .or. e(l) .eq. 0.0d0) go to 120 188c 189c apply the transformation. 190c 191 do 90 i = lp1, n 192 work(i) = 0.0d0 193 90 continue 194 do 100 j = lp1, p 195 call daxpy(n-l,e(j),x(lp1,j),1,work(lp1),1) 196 100 continue 197 do 110 j = lp1, p 198 call daxpy(n-l,-e(j)/e(lp1),work(lp1),1,x(lp1,j),1) 199 110 continue 200 120 continue 201 if (.not.wantv) go to 140 202c 203c place the transformation in v for subsequent 204c back multiplication. 205c 206 do 130 i = lp1, p 207 v(i,l) = e(i) 208 130 continue 209 140 continue 210 150 continue 211 160 continue 212 170 continue 213c 214c set up the final bidiagonal matrix or order m. 215c 216 m = min0(p,n+1) 217 nctp1 = nct + 1 218 nrtp1 = nrt + 1 219 if (nct .lt. p) s(nctp1) = x(nctp1,nctp1) 220 if (n .lt. m) s(m) = 0.0d0 221 if (nrtp1 .lt. m) e(nrtp1) = x(nrtp1,m) 222 e(m) = 0.0d0 223c 224c if required, generate u. 225c 226 if (.not.wantu) go to 300 227 if (ncu .lt. nctp1) go to 200 228 do 190 j = nctp1, ncu 229 do 180 i = 1, n 230 u(i,j) = 0.0d0 231 180 continue 232 u(j,j) = 1.0d0 233 190 continue 234 200 continue 235 if (nct .lt. 1) go to 290 236 do 280 ll = 1, nct 237 l = nct - ll + 1 238 if (s(l) .eq. 0.0d0) go to 250 239 lp1 = l + 1 240 if (ncu .lt. lp1) go to 220 241 do 210 j = lp1, ncu 242 t = -ddot(n-l+1,u(l,l),1,u(l,j),1)/u(l,l) 243 call daxpy(n-l+1,t,u(l,l),1,u(l,j),1) 244 210 continue 245 220 continue 246 call dscal(n-l+1,-1.0d0,u(l,l),1) 247 u(l,l) = 1.0d0 + u(l,l) 248 lm1 = l - 1 249 if (lm1 .lt. 1) go to 240 250 do 230 i = 1, lm1 251 u(i,l) = 0.0d0 252 230 continue 253 240 continue 254 go to 270 255 250 continue 256 do 260 i = 1, n 257 u(i,l) = 0.0d0 258 260 continue 259 u(l,l) = 1.0d0 260 270 continue 261 280 continue 262 290 continue 263 300 continue 264c 265c if it is required, generate v. 266c 267 if (.not.wantv) go to 350 268 do 340 ll = 1, p 269 l = p - ll + 1 270 lp1 = l + 1 271 if (l .gt. nrt) go to 320 272 if (e(l) .eq. 0.0d0) go to 320 273 do 310 j = lp1, p 274 t = -ddot(p-l,v(lp1,l),1,v(lp1,j),1)/v(lp1,l) 275 call daxpy(p-l,t,v(lp1,l),1,v(lp1,j),1) 276 310 continue 277 320 continue 278 do 330 i = 1, p 279 v(i,l) = 0.0d0 280 330 continue 281 v(l,l) = 1.0d0 282 340 continue 283 350 continue 284c 285c main iteration loop for the singular values. 286c 287 mm = m 288 iter = 0 289 360 continue 290c 291c quit if all the singular values have been found. 292c 293c ...exit 294 if (m .eq. 0) go to 620 295c 296c if too many iterations have been performed, set 297c flag and return. 298c 299 if (iter .lt. maxit) go to 370 300 info = m 301c ......exit 302 go to 620 303 370 continue 304c 305c this section of the program inspects for 306c negligible elements in the s and e arrays. on 307c completion the variables kase and l are set as follows. 308c 309c kase = 1 if s(m) and e(l-1) are negligible and l.lt.m 310c kase = 2 if s(l) is negligible and l.lt.m 311c kase = 3 if e(l-1) is negligible, l.lt.m, and 312c s(l), ..., s(m) are not negligible (qr step). 313c kase = 4 if e(m-1) is negligible (convergence). 314c 315 do 390 ll = 1, m 316 l = m - ll 317c ...exit 318 if (l .eq. 0) go to 400 319 test = dabs(s(l)) + dabs(s(l+1)) 320 ztest = test + dabs(e(l)) 321 if (ztest .ne. test) go to 380 322 e(l) = 0.0d0 323c ......exit 324 go to 400 325 380 continue 326 390 continue 327 400 continue 328 if (l .ne. m - 1) go to 410 329 kase = 4 330 go to 480 331 410 continue 332 lp1 = l + 1 333 mp1 = m + 1 334 do 430 lls = lp1, mp1 335 ls = m - lls + lp1 336c ...exit 337 if (ls .eq. l) go to 440 338 test = 0.0d0 339 if (ls .ne. m) test = test + dabs(e(ls)) 340 if (ls .ne. l + 1) test = test + dabs(e(ls-1)) 341 ztest = test + dabs(s(ls)) 342 if (ztest .ne. test) go to 420 343 s(ls) = 0.0d0 344c ......exit 345 go to 440 346 420 continue 347 430 continue 348 440 continue 349 if (ls .ne. l) go to 450 350 kase = 3 351 go to 470 352 450 continue 353 if (ls .ne. m) go to 460 354 kase = 1 355 go to 470 356 460 continue 357 kase = 2 358 l = ls 359 470 continue 360 480 continue 361 l = l + 1 362c 363c perform the task indicated by kase. 364c 365 go to (490,520,540,570), kase 366c 367c deflate negligible s(m). 368c 369 490 continue 370 mm1 = m - 1 371 f = e(m-1) 372 e(m-1) = 0.0d0 373 do 510 kk = l, mm1 374 k = mm1 - kk + l 375 t1 = s(k) 376 call drotg(t1,f,cs,sn) 377 s(k) = t1 378 if (k .eq. l) go to 500 379 f = -sn*e(k-1) 380 e(k-1) = cs*e(k-1) 381 500 continue 382 if (wantv) call drot(p,v(1,k),1,v(1,m),1,cs,sn) 383 510 continue 384 go to 610 385c 386c split at negligible s(l). 387c 388 520 continue 389 f = e(l-1) 390 e(l-1) = 0.0d0 391 do 530 k = l, m 392 t1 = s(k) 393 call drotg(t1,f,cs,sn) 394 s(k) = t1 395 f = -sn*e(k) 396 e(k) = cs*e(k) 397 if (wantu) call drot(n,u(1,k),1,u(1,l-1),1,cs,sn) 398 530 continue 399 go to 610 400c 401c perform one qr step. 402c 403 540 continue 404c 405c calculate the shift. 406c 407 scale = dmax1(dabs(s(m)),dabs(s(m-1)),dabs(e(m-1)), 408 * dabs(s(l)),dabs(e(l))) 409 sm = s(m)/scale 410 smm1 = s(m-1)/scale 411 emm1 = e(m-1)/scale 412 sl = s(l)/scale 413 el = e(l)/scale 414 b = ((smm1 + sm)*(smm1 - sm) + emm1**2)/2.0d0 415 c = (sm*emm1)**2 416 shift = 0.0d0 417 if (b .eq. 0.0d0 .and. c .eq. 0.0d0) go to 550 418 shift = dsqrt(b**2+c) 419 if (b .lt. 0.0d0) shift = -shift 420 shift = c/(b + shift) 421 550 continue 422 f = (sl + sm)*(sl - sm) + shift 423 g = sl*el 424c 425c chase zeros. 426c 427 mm1 = m - 1 428 do 560 k = l, mm1 429 call drotg(f,g,cs,sn) 430 if (k .ne. l) e(k-1) = f 431 f = cs*s(k) + sn*e(k) 432 e(k) = cs*e(k) - sn*s(k) 433 g = sn*s(k+1) 434 s(k+1) = cs*s(k+1) 435 if (wantv) call drot(p,v(1,k),1,v(1,k+1),1,cs,sn) 436 call drotg(f,g,cs,sn) 437 s(k) = f 438 f = cs*e(k) + sn*s(k+1) 439 s(k+1) = -sn*e(k) + cs*s(k+1) 440 g = sn*e(k+1) 441 e(k+1) = cs*e(k+1) 442 if (wantu .and. k .lt. n) 443 * call drot(n,u(1,k),1,u(1,k+1),1,cs,sn) 444 560 continue 445 e(m-1) = f 446 iter = iter + 1 447 go to 610 448c 449c convergence. 450c 451 570 continue 452c 453c make the singular value positive. 454c 455 if (s(l) .ge. 0.0d0) go to 580 456 s(l) = -s(l) 457 if (wantv) call dscal(p,-1.0d0,v(1,l),1) 458 580 continue 459c 460c order the singular value. 461c 462 590 if (l .eq. mm) go to 600 463c ...exit 464 if (s(l) .ge. s(l+1)) go to 600 465 t = s(l) 466 s(l) = s(l+1) 467 s(l+1) = t 468 if (wantv .and. l .lt. p) 469 * call dswap(p,v(1,l),1,v(1,l+1),1) 470 if (wantu .and. l .lt. n) 471 * call dswap(n,u(1,l),1,u(1,l+1),1) 472 l = l + 1 473 go to 590 474 600 continue 475 iter = 0 476 m = m - 1 477 610 continue 478 go to 360 479 620 continue 480 return 481 end 482