Esope source files

   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
  21 :          call levma3(x,y,sig,ndata,a,ma,alpha,beta,nca,
  22 :      + ytest,dyda,chisq)
  23 :          deini = .false.
  24 :          do 113 j=1,ma
  25 :             aold(j)=a(j)
  26 :  113     continue
  27 :          goto 202
  28 :       endif
  29 : 
  30 :  201  ochisq = chisq
  31 : 
  32 : * calcul de chi2
  32 :       call levma3(x,y,sig,ndata,a,ma,covar,da,nca,
  33 :      + ytest,dyda,chisq)
  34 :  203  continue
  35 : 
  36 :       write(6,*) 'ochisq ',ochisq, ' chisq', chisq, 'alamda  ',alamda
  37 :       if (chisq.le.ochisq) then 
  38 : * si oui prend la nouvelle solution
  39 :          alamda = 0.1d0*alamda
  40 :          ochisq = chisq
  41 :          do 218 j=1,ma
  42 :            do 217 k=1,ma
  43 :              alpha(j,k)=covar(j,k)
  44 :  217       continue
  45 :            beta(j) = da(j)
  46 :            aold(j)=a(j)
  47 :  218     continue
  48 :       else
  49 : * si non augmente alamda et retour
  50 :         write(6,*) ' alamda augmente ' ,alamda
  51 :         alamda = 10.d0*alamda
  52 :         chisq = ochisq
  53 :         do 219 j= 1,ma
  54 :           a(j) = aold(j)
  55 :  219    continue
  56 :       endif
  57 : 
  58 : * propose nouveaux parametres 
  59 :  202  continue
  60 : * matrice d ajustment linearisee : augmente les elements diagonaux
  61 :       do 215 j=1,ma
  62 :          do 214 k=1,ma
  63 :             covar(j,k)=alpha(j,k)
  64 :  214     continue
  65 :          covar(j,j)=alpha(j,j)*(1. + alamda)
  66 :          da(j) = beta(j)
  67 :  215  continue
  68 : 
  69 : c      write(6,*) (covar(1,j),j=1,4) 
  70 : c      write(6,*) (covar(2,j),j=1,4) 
  71 : c      write(6,*) (covar(3,j),j=1,4) 
  72 : c      write(6,*) (covar(4,j),j=1,4)
  73 : c      write(6,*) 'da', (da(j),j=1,4) 
  73 :       call gaussk(covar,ma,nca,da,1,1)
  74 :       if (alamda.eq.0.) then
  74 :         call covsrt(covar,nca,ma)
  75 :         return
  76 :       endif
  77 : 
  78 : c      write(6,*) 'aold',(aold(lista(j)),j=1,4) 
  79 : c      write(6,*) 'a',(a(lista(j)),j=1,4) 
  80 : c      write(6,*) (lista(j),j=1,4)
  81 :       do 216 j =1,ma
  82 :         a(j) = a(j) + da(j)
  83 :  216   continue
  84 : c      write(6,*) 'da',(da(j),j=1,3) 
  85 :       return
  86 :       end
  87 :  
  88 :