Télécharger levmar.eso

Retour à la liste

Numérotation des lignes :

  1. C LEVMAR SOURCE KICH 16/07/12 21:15:04 9027
  2. SUBROUTINE LEVMAR
  3. *
  4. * moteur pour algorithme levenberg-marquardt
  5. *
  6. IMPLICIT INTEGER(I-N)
  7. IMPLICIT REAL*8(A-H,O-Z)
  8.  
  9. -INC CCOPTIO
  10. -INC CCREEL
  11. -INC SMLENTI
  12. -INC SMLREEL
  13. SEGMENT MLREE4.MLREEL,MLREE5.MLREEL,MLREE6.MLREEL,MLREE8.MLREEL,
  14. + MLREE9.MLREEL,MLRE11.MLREEL,MLRE12.MLREEL,MLRE13.MLREEL
  15.  
  16. PARAMETER (NTYP=11,NPRO = 10)
  17. CHARACTER*4 LISTYP(NTYP)
  18. integer iprocm(NPRO,2)
  19. real*8 xmordo(NPRO)
  20. logical dpro
  21. DATA LISTYP / 'ABSC', 'ORDO', 'PARA', 'TORD', 'DYDA',
  22. & 'SIGM', 'PROC', 'ALPH', 'BETA', 'CHI2',
  23. & 'PARO' /
  24.  
  25. data iprocm/0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/
  26. SAVE iprocm
  27. data xmordo/1.D0,1.D0,1.D0,1.D0,1.D0,1.D0,1.D0,1.D0,1.D0,1.D0/
  28. SAVE xmordo
  29. c
  30. dpro =.true.
  31. lapro = 0
  32. c
  33. jg = NTYP
  34. segini mlenti
  35. c
  36. call LIROBJ('LISTREEL',ip1,0,iret)
  37. IF (IERR.NE.0) GOTO 900
  38. if(iret.eq.1) then
  39. lect(4) = ip1
  40. call LIROBJ('LISTREEL',ip1,1,iret)
  41. IF (IERR.NE.0) GOTO 900
  42. lect(5) = ip1
  43. endif
  44. c
  45. idon = 0
  46.  
  47. 10 continue
  48. CALL LIRMOT(LISTYP,NTYP,iplac,0)
  49. IF (IERR.NE.0) GOTO 900
  50. if (iplac.gt.0) then
  51. if (iplac.ne.7) then
  52. call LIROBJ('LISTREEL',ip1,1,iret)
  53. else if (iplac.eq.7) then
  54. call LIROBJ('PROCEDUR',ip1,1,iret)
  55. do kpro = 1,NPRO
  56. if (iprocm(kpro,1).eq.ip1) then
  57. iprocc = iprocm(kpro,2)
  58. lapro = kpro
  59. dpro = .false.
  60. endif
  61. if (iprocm(kpro,1).eq.0.and.dpro) then
  62. iprocm(kpro,1) = ip1
  63. iprocc = 0
  64. lapro = kpro
  65. dpro = .false.
  66. endif
  67. enddo
  68. endif
  69. IF (IERR.NE.0) GOTO 900
  70. lect(iplac) = ip1
  71. idon = idon + 1
  72. goto 10
  73. endif
  74.  
  75.  
  76.  
  77. if (iprocc.eq.0) then
  78. if (lect(7).le.0) then
  79. moterr(1:8) ='PROCEDUR'
  80. call erreur(37)
  81. GOTO 900
  82. endif
  83. mlree1 = lect(1)
  84. mlree2 = lect(2)
  85. mlree6 = lect(6)
  86. mlree3 = lect(3)
  87. if (mlree1.le.0 .or. mlree2.le.0 .or. mlree6.le.0 .or.
  88. & mlree3.le.0) then
  89. call erreur(201)
  90. GOTO 900
  91. endif
  92. segact,mlree1,mlree2,mlree6
  93. n1 = mlree1.prog(/1)
  94. n2 = mlree2.prog(/1)
  95. n6 = mlree6.prog(/1)
  96. if (n1.LE.0 .or. n1.NE.n2 .or. n1.NE.n6 .or. n2.NE.n6) then
  97. segdes,mlree1,mlree2,mlree6
  98. call ERREUR(217)
  99. GOTO 900
  100. endif
  101. segact,mlree3
  102. ma = mlree3.prog(/1)
  103. segdes,mlree3
  104. if (ma.LE.0) then
  105. segdes,mlree1,mlree2,mlree6
  106. call ERREUR(21)
  107. GOTO 900
  108. endif
  109. segini,mlre11=mlree3
  110. lect(3) = mlre11
  111. jg = ma
  112. segini mlre12
  113. lect(11) = mlre12
  114. jg = ma*ma
  115. segini mlree8
  116. lect(8) = mlree8
  117. jg = ma
  118. segini mlree9
  119. lect(9) = mlree9
  120. jg = 2
  121. segini mlreel
  122. prog(1) = -1.
  123. prog(2) = -1.
  124. lect(10) = mlreel
  125. c
  126. * dimensionner critere
  127. ip2 = mlree2
  128. call ecrobj('LISTREEL',ip2)
  129. call opobje(14)
  130. call maximu(1)
  131. call lirree(umordo,1,iretou)
  132. if (umordo.GT.XSPETI.and.umordo.lt.1.d0) then
  133. xmordo(lapro) = umordo * umordo
  134. else
  135. xmordo(lapro) = 1.D0
  136. endif
  137. c
  138. else
  139. c* else if (iprocc.gt.0) then
  140. do idon = 1, NTYP
  141. if (lect(idon).LE.0) then
  142. * write(ioimp,*) '!! IPROCC =',iprocc
  143. * iprocc = 0
  144. iprocm(lapro,1) = 0
  145. iprocm(lapro,2) = 0
  146. if (idon.eq.7) then
  147. moterr(1:8) ='PROCEDUR'
  148. call erreur(37)
  149. else
  150. call erreur(201)
  151. endif
  152. GOTO 900
  153. endif
  154. enddo
  155.  
  156. call levma1(mlenti)
  157. endif
  158.  
  159. mlreel = lect(10)
  160. chi2 = prog(1)
  161. alamda = prog(2)
  162. *
  163. * write(ioimp,*)'levm',xmordo(lapro),chi2,(abs(chi2)/xmordo(lapro))
  164. if(((abs(chi2)/xmordo(lapro)).lt.1.e-6.and.chi2.gt.0)
  165. &.or.iprocc.gt.100.or.alamda.gt.1.e11) goto 80
  166. *
  167. iprocc = iprocc + 1
  168. iprocm(lapro,2) = iprocc
  169. ip1 = lect(11)
  170. call nomobj('LISTREEL','WLEVPARO',ip1)
  171. call ecrobj('LISTREEL',ip1)
  172. call ECRCHA('PARO')
  173. ip1 = lect(9)
  174. call nomobj('LISTREEL','WLEVBETA',ip1)
  175. call ecrobj('LISTREEL',ip1)
  176. call ECRCHA('BETA')
  177. ip1 = lect(8)
  178. call nomobj('LISTREEL','WLEVALPH',ip1)
  179. call ecrobj('LISTREEL',ip1)
  180. call ECRCHA('ALPH')
  181. ip1 = lect(7)
  182. call ecrobj('PROCEDUR',ip1)
  183. call ECRCHA('PROC')
  184. ip1 = lect(10)
  185. call nomobj('LISTREEL','WLEVCHI2',ip1)
  186. call ecrobj('LISTREEL',ip1)
  187. call ECRCHA('CHI2')
  188. ip1 = lect(6)
  189. call ecrobj('LISTREEL',ip1)
  190. call ECRCHA('SIGM')
  191. ip1 = lect(3)
  192. call nomobj('LISTREEL','WLEVPARA',ip1)
  193. call ecrobj('LISTREEL',ip1)
  194. call ECRCHA('PARA')
  195. ip1 = lect(2)
  196. call ecrobj('LISTREEL',ip1)
  197. call ECRCHA('ORDO')
  198. ip1 = lect(1)
  199. call ecrobj('LISTREEL',ip1)
  200. call ECRCHA('ABSC')
  201. call ECRCHA('LEVM')
  202. ip1 = lect(3)
  203. call ecrobj('LISTREEL',ip1)
  204. ip1 = lect(1)
  205. call ecrobj('LISTREEL',ip1)
  206. ip1 = lect(7)
  207. call ecrobj('PROCEDUR',ip1)
  208. GOTO 900
  209.  
  210. 80 continue
  211. if (alamda.gt.1.e11) write(ioimp,*) 'je ne sais pas faire mieux'
  212. * oubli des objets intermediares
  213. ip1 = lect(8)
  214. CALL NOMOBJ('ANNULE ','WLEVALPH',ip1)
  215. ip1 = lect(9)
  216. CALL NOMOBJ('ANNULE ','WLEVBETA',ip1)
  217. ip1 = lect(11)
  218. CALL NOMOBJ('ANNULE ','WLEVPARO',ip1)
  219.  
  220. mlree1 = lect(1)
  221. mlree2 = lect(2)
  222. mlree3 = lect(6)
  223. segdes mlree1,mlree2,mlree3
  224. mlree1 = lect(3)
  225. segini,mlree2=mlree1
  226. ip1 = mlree1
  227. CALL NOMOBJ('ANNULE ','WLEVPARA',ip1)
  228. ipout = mlree2
  229. mlreel = lect(10)
  230. chi2 = prog(1)
  231. ip1 = mlreel
  232. CALL NOMOBJ('ANNULE ','WLEVPARA',ip1)
  233. call ECRREE(chi2)
  234. call ECROBJ('LISTREEL', ipout)
  235. iprocm(lapro,1) = 0
  236. iprocm(lapro,2) = 0
  237.  
  238. 900 CONTINUE
  239. segsup mlenti
  240.  
  241. RETURN
  242. END
  243.  
  244.  
  245.  
  246.  
  247.  

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