Télécharger levma2.eso

Retour à la liste

Numérotation des lignes :

  1. C LEVMA2 SOURCE CHAT 05/01/13 01:15:25 5004
  2. subroutine levma2(x,y,sig, ndata,a,ma,
  3. +ytest,dyda,chisq,alamda,alpha,beta,aold,covar,da)
  4. c
  5. c appele par mrqmin pour evaluer la matrice alpha et le vecteur beta
  6. IMPLICIT INTEGER(I-N)
  7. IMPLICIT REAL*8 (A-H,O-Z)
  8. parameter (mmax=20)
  9. real*8 x(ndata),y(ndata),sig(ndata),alpha(ma,ma),beta(ma)
  10. real*8 covar(ma,ma),da(ma)
  11. real*8 dyda(ma*ndata),a(ma),ytest(ndata),aold(ma)
  12. logical deini
  13.  
  14. nca = ma
  15. deini = .false.
  16. if (chisq.lt.0) deini = .true.
  17.  
  18. c write(6,*) ytest(1),dyda(1),dyda(2),dyda(3),dyda(4)
  19. * initialisation
  20. if (deini) then
  21. alamda = 0.001d0
  22. call levma3(x,y,sig,ndata,a,ma,alpha,beta,nca,
  23. + ytest,dyda,chisq)
  24. deini = .false.
  25. do 113 j=1,ma
  26. aold(j)=a(j)
  27. 113 continue
  28. goto 202
  29. endif
  30.  
  31. 201 ochisq = chisq
  32.  
  33. * calcul de chi2
  34. call levma3(x,y,sig,ndata,a,ma,covar,da,nca,
  35. + ytest,dyda,chisq)
  36. 203 continue
  37.  
  38. write(6,*) 'ochisq ',ochisq, ' chisq', chisq, 'alamda ',alamda
  39. if (chisq.le.ochisq) then
  40. * si oui prend la nouvelle solution
  41. alamda = 0.1d0*alamda
  42. ochisq = chisq
  43. do 218 j=1,ma
  44. do 217 k=1,ma
  45. alpha(j,k)=covar(j,k)
  46. 217 continue
  47. beta(j) = da(j)
  48. aold(j)=a(j)
  49. 218 continue
  50. else
  51. * si non augmente alamda et retour
  52. write(6,*) ' alamda augmente ' ,alamda
  53. alamda = 10.d0*alamda
  54. chisq = ochisq
  55. do 219 j= 1,ma
  56. a(j) = aold(j)
  57. 219 continue
  58. endif
  59.  
  60. * propose nouveaux parametres
  61. 202 continue
  62. * matrice d ajustment linearisee : augmente les elements diagonaux
  63. do 215 j=1,ma
  64. do 214 k=1,ma
  65. covar(j,k)=alpha(j,k)
  66. 214 continue
  67. covar(j,j)=alpha(j,j)*(1. + alamda)
  68. da(j) = beta(j)
  69. 215 continue
  70.  
  71. c write(6,*) (covar(1,j),j=1,4)
  72. c write(6,*) (covar(2,j),j=1,4)
  73. c write(6,*) (covar(3,j),j=1,4)
  74. c write(6,*) (covar(4,j),j=1,4)
  75. c write(6,*) 'da', (da(j),j=1,4)
  76. call gaussk(covar,ma,nca,da,1,1)
  77. if (alamda.eq.0.) then
  78. call covsrt(covar,nca,ma)
  79. return
  80. endif
  81.  
  82. c write(6,*) 'aold',(aold(lista(j)),j=1,4)
  83. c write(6,*) 'a',(a(lista(j)),j=1,4)
  84. c write(6,*) (lista(j),j=1,4)
  85. do 216 j =1,ma
  86. a(j) = a(j) + da(j)
  87. 216 continue
  88. c write(6,*) 'da',(da(j),j=1,3)
  89. return
  90. end
  91.  
  92.  
  93.  

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