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