Sources Esope | Cast3M

levma2
C LEVMA2    SOURCE    CHAT      05/01/13    01:15:25     5004
       subroutine levma2(x,y,sig, ndata,a,ma,
     +ytest,dyda,chisq,alamda,alpha,beta,aold,covar,da)
c
c appele par mrqmin pour evaluer la matrice alpha et le vecteur beta
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      parameter (mmax=20)
      real*8 x(ndata),y(ndata),sig(ndata),alpha(ma,ma),beta(ma)
      real*8 covar(ma,ma),da(ma)
      real*8 dyda(ma*ndata),a(ma),ytest(ndata),aold(ma)
      logical deini
 
      nca = ma
      deini = .false.
      if (chisq.lt.0) deini = .true.
 
c       write(6,*) ytest(1),dyda(1),dyda(2),dyda(3),dyda(4)
* initialisation
      if (deini) then
         alamda = 0.001d0
         call levma3(x,y,sig,ndata,a,ma,alpha,beta,nca,
     + ytest,dyda,chisq)
         deini = .false.
         do 113 j=1,ma
            aold(j)=a(j)
 113     continue
         goto 202
      endif
 
 201  ochisq = chisq
 
* calcul de chi2
      call levma3(x,y,sig,ndata,a,ma,covar,da,nca,
     + ytest,dyda,chisq)
 203  continue
 
      write(6,*) 'ochisq ',ochisq, ' chisq', chisq, 'alamda  ',alamda
      if (chisq.le.ochisq) then
* si oui prend la nouvelle solution
         alamda = 0.1d0*alamda
         ochisq = chisq
         do 218 j=1,ma
           do 217 k=1,ma
             alpha(j,k)=covar(j,k)
 217       continue
           beta(j) = da(j)
           aold(j)=a(j)
 218     continue
      else
* si non augmente alamda et retour
        write(6,*) ' alamda augmente ' ,alamda
        alamda = 10.d0*alamda
        chisq = ochisq
        do 219 j= 1,ma
          a(j) = aold(j)
 219    continue
      endif
 
* propose nouveaux parametres
 202  continue
* matrice d ajustment linearisee : augmente les elements diagonaux
      do 215 j=1,ma
         do 214 k=1,ma
            covar(j,k)=alpha(j,k)
 214     continue
         covar(j,j)=alpha(j,j)*(1. + alamda)
         da(j) = beta(j)
 215  continue
 
c      write(6,*) (covar(1,j),j=1,4)
c      write(6,*) (covar(2,j),j=1,4)
c      write(6,*) (covar(3,j),j=1,4)
c      write(6,*) (covar(4,j),j=1,4)
c      write(6,*) 'da', (da(j),j=1,4)
      call gaussk(covar,ma,nca,da,1,1)
      if (alamda.eq.0.) then
        call covsrt(covar,nca,ma)
        return
      endif
 
c      write(6,*) 'aold',(aold(lista(j)),j=1,4)
c      write(6,*) 'a',(a(lista(j)),j=1,4)
c      write(6,*) (lista(j),j=1,4)
      do 216 j =1,ma
        a(j) = a(j) + da(j)
 216   continue
c      write(6,*) 'da',(da(j),j=1,3)
      return
      end