Télécharger invsvd.eso

Retour à la liste

Numérotation des lignes :

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

© Cast3M 2003 - Tous droits réservés.
Mentions légales