Sources Esope | Cast3M

itgl2l
C ITGL2L    SOURCE    CHAT      05/01/13    00:44:39     5004
CCC
C **********************************************************************
CCC
      SUBROUTINE INTEGRA_LS_2LEVELS(XTRI,XINI,NDIMX,LAM,NMODEL,TOLREL,
     .                        jmax,nescri,ues,nnumer,deltax,kerre,jtot)
 
      IMPLICIT INTEGER(I-N)
      integer nmodel,nitmax,nescri,nnumer,ues,j,kerre,
     .        ndimv,jmax,i,ii,k,n,ndimx,uesb,jl,jtot
      real*8  xtri(ndimx),xini(ndimx),lam,tolrel,deltax,
     .        x(5),dx(5),r(5),err,normainf,norm_l2_2,xnew(5),rnew(5),
     .        cotaere,lam0,lam1,lam2,f0,fp0,f1,f2,f3,s1,s2,aa,bb,cc,
     .        sig(4),vecm(4),err0,dla,dlb,ddll,beta
      logical escr,last
      integer             level,kmax
      common /linesearch/ level,kmax
      integer             augla
      real*8                    c
      common /auglagrang1/ augla
      common /auglagrang2/ c
 
c muestra tensiones
c      if (nescri.eq.1) then
c       call der_enerelas_dpral(xini,x,nmodel)
c       write(ues,*)' Tensiones iniciales'
c       write(ues,*)(x(i),i=1,3)
c       if (ndimx.eq.4) write(ues,*)' Variable interna inicial',x(4)
c      endif
 
c muestra iteraciones de line-search
      uesb=299
      escr=.false.
      if (nescri.eq.1) escr=.true.
 
c parametros generales
      kerre=0
      n=ndimx+1
      ndimv=1
      call equv(x,xini,ndimx)
        x(n)=lam
      beta=1.D-4
      j=0
 
      jtot=-1
      last=.false.
      jl=0
 
C CALCULOS PARA LAMBDA; SUPONEMOS SIG_INI=SIG_TRIAL
   2  continue
      if ((augla.eq.1).and.(jl.eq.0)) then
         jl=1
         goto 10
      endif
      jl=jl+1
C A; B
      call actualiza_lambda(x,ndimx,nmodel,lam,dla,dlb)
      ddll=-dla/dlb
      lam=lam+ddll
      if (lam.lt.0) then
        write(ues,*)' GAMMA negative !',lam,' => ',0.D0
        write(ues,*)'             dll:',ddll
        lam=0
      endif
C CHECK DL => GOTO END
      jtot=jtot+1
c      if (nescri.eq.1) write(ues,999)jtot,jl,lam,ddll,dla
      if (last) goto 100
      if (abs(ddll).lt.tolrel) last=.true.
      if (jl.gt.50+1) then
        write(ues,*)' ERROR: loop gamma no ha convergido en', 50
        kerre=1
        goto 100
      endif
 
C COMPUTE SIGMA_J FROM LAMBDA_J WITH LINE-SEARCH
  10  continue
      j=j+1
      call residuo_2levels(x,r,ndimx,ndimx,nmodel,xtri,sig,vecm,lam)
      f0=norm_l2_2(r,ndimx)
      fp0=-2.D0*f0
      call direction_2levels
     .     (x,r,dx,ndimx,ndimx,nmodel,sig,vecm,nnumer,deltax,lam)
      err0=normainf(dx,ndimx)
      lam0=1.D0
      call addvec(xnew,x,lam0,dx,ndimx)
      call residuo_2levels
     .     (xnew,rnew,ndimx,ndimx,nmodel,xtri,sig,vecm,lam)
      f1=norm_l2_2(rnew,ndimx)
      cotaere=f0+beta*lam0*fp0
      lam1=lam0
      if (escr) write(uesb,*)' more ',0,lam0
      if (f1.lt.cotaere) goto 80
      if ((f0.lt.tolrel).and.(f1.lt.tolrel)) goto 80
      if (kmax.le.0) goto 70
 
      k=1
      lam1=(-fp0)/2.D0/(f1-f0-fp0)
      if (lam1.lt.0.1D0) lam1=0.1D0
      call addvec(xnew,x,lam1,dx,ndimx)
      call residuo_2levels(xnew,rnew,n,ndimx,nmodel,xtri,sig,vecm,lam)
      f2=norm_l2_2(rnew,ndimx)
      cotaere=f0+beta*lam1*fp0
      if (escr) write(uesb,*)' more ',k,lam1
      if (f2.lt.cotaere) goto 80
      if (kmax.le.1) goto 70
 
  20  continue
      k=k+1
      lam2=(-fp0*lam1**2)/2.D0/(f2-f0-fp0*lam1)
      if (lam2.lt.(0.1D0*lam1)) lam2=0.1D0*lam1
      if (lam2.gt.(0.5D0*lam1)) lam2=0.5D0*lam1
      call addvec(xnew,x,lam2,dx,ndimx)
      call residuo_2levels(xnew,rnew,n,ndimx,nmodel,xtri,sig,vecm,lam)
      f3=norm_l2_2(rnew,ndimx)
      cotaere=f0+beta*lam2*fp0
      f1=f2
      f2=f3
      lam0=lam1
      lam1=lam2
      if (escr) write(uesb,*)' more ',k,lam2
      if (f3.lt.cotaere) goto 80
      if (k.lt.kmax) goto 20
  70  continue
      if (escr) write(uesb,*)' more line search',kmax
  80  continue
      err=err0*lam1/normainf(xnew,ndimx)
      call equv(x,xnew,ndimx)
c MUESTRA ITERACIONES EN SIGMA_J
c      if (nescri.eq.1) then
c         call Invari_t(x,3,aa)
c         call Invari_q(x,3,bb)
c         write(ues,999)jtot,j,(x(i),i=1,3),lam,err
c         write(ues,999)j,aa,bb,(x(i),i=4,n),err
c 999    format(2I4,'  ',6(E10.4,1x))
c        call der_enerelas_dpral(x,xnew,nmodel)
c        call Invari_t(xnew,3,cc)
c        call Invari_q(xnew,3,bb)
c        call Invari_p(xnew,3,aa)
c        write(ues,1999)j,aa,bb,cc*45.D0/atan(1.D0),err,rnew(n)
c1999    format(I3,' ',5(E14.8,1x))
c      endif
      if (err.lt.tolrel) goto 99
      if (j.lt.jmax) goto 10
      kerre=1
      write(ues,*)' ERROR: Loop interno sigmas no ha convergido'
  99  continue
c FIN ITERACIONES EN SIGMA_J
      jtot=jtot+j
      j=0
      goto 2
 
c RESULTADOS FINALES
 100  continue
      call equv(xini,x,ndimx)
 
c MUESTRA RESULTADOS
c      if (nescri.eq.1) then
c        write(ues,*)' Deform elas finales'
c        write(ues,*)(x(i),i=1,3)
c        write(ues,*)' Deform plas finales'
c        write(ues,*)(xtri(i)-x(i),i=1,3)
c        if (ndimx.eq.4) write(ues,*)' Variable interna final',x(4)
c        call der_enerelas_dpral(x,xini,nmodel)
c        write(ues,*)' Tensiones finales'
c        write(ues,*)(xini(i),i=1,3)
c        if (ndimx.eq.4) write(ues,*)' Variable interna inicial',x(4)
c        call equv(xini,x,ndimx)
c       endif
      end