Sources Esope | Cast3M

itgldp
C ITGLDP    SOURCE    CB215821  16/04/21    21:17:18     8920
C  LINESEARCH for FEFP models
C **********************************************************************
CCC
cc Trabaja en direcciones principales y las incognitas son DEFO ELAS!
 
      SUBROUTINE INTEGRA_LS_DPRAL(XTRI,XINI,NDIMX,LAM,NMODEL,TOLREL,
     .                           jmax,nescri,ues,nnumer,deltax,kerre,j)
 
      IMPLICIT INTEGER(I-N)
      integer nmodel,nitmax,nescri,nnumer,ues,j,kerre,
     .        ndimv,jmax,i,ii,k,n,ndimx,uesb
      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,aux,beta
      logical escr
c kmax=0 => no line search
      integer             level,kmax
      common /linesearch/ level,kmax
 
c muestra tensiones
c      if (nescri.eq.1) then
c       write(ues,*)' Variables iniciales'
c       write(ues,*)(xini(i),i=1,ndimx)
c       write(ues,*)' Variables trial'
c       write(ues,*)(xtri(i),i=1,ndimx)
c       write(ues,*)' Lambda inicial'
c       write(ues,*)lam
c       write(ues,*)' Variables varias'
c       write(ues,*)level,kmax,nmodel,nnumer,deltax
c       call der_enerelas_dpral(xini,x,nmodel)
c       write(ues,*)' Tensiones iniciales'
c       write(ues,*)(x(i),i=1,3)
c       call Invari_t(x,3,cc)
c       call Invari_J2(x,3,bb)
c       call Invari_I1(x,3,aa)
c       write(ues,999)0,aa,bb,cc*45.D0/atan(1.D0)
c       if (ndimx.eq.4) write(ues,*)' Variable interna inicial',x(4)
c      endif
 
c muestra iteraciones de line-search
      uesb=199
      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
 
  10  continue
      j=j+1
      call residuo_dpral(x,r,n,ndimx,nmodel,xtri,sig,vecm)
      if (nescri.eq.2) then
       write(ues,*)' TRAS residuo_dpral'
       write(ues,*)(x(i),i=1,n)
       write(ues,*)(r(i),i=1,n)
       write(ues,*)(vecm(i),i=1,ndimx)
      endif
 
c control constrain gamma positiva
      if (x(n).le.0.D0) then
       aux=0.D0
       do i=1,ndimx
        aux=aux+vecm(i)*(x(i)-xtri(i))
       enddo
       if (aux.gt.0.D0)
     .    write(ues,*)' ERROR: Active constrain gam>0 !',x(n),aux
c        if (x(n).lt.0.D0)
c     .    write(ues,*)' GAMMA negative !',x(n),' => ',0.D0
       x(n)=0.D0
      endif
 
c line-search cero
      f0=norm_l2_2(r,n)
      fp0=-2.D0*f0
      call direction_dpral(x,r,dx,n,ndimx,nmodel,sig,vecm,nnumer,deltax)
c      if (nescri.eq.2) then
c       write(ues,*)' TRAS direction_dpral'
c       write(ues,*)(x(i),i=1,n)
c       write(ues,*)(r(i),i=1,n)
c       write(ues,*)(dx(i),i=1,n)
c      endif
      err0=normainf(dx,ndimx+1)
      lam0=1.D0
      call addvec(xnew,x,lam0,dx,n)
      call residuo_dpral(xnew,rnew,n,ndimx,nmodel,xtri,sig,vecm)
      f1=norm_l2_2(rnew,n)
      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
c line-search uno
      k=1
      lam1=(-fp0)/2.D0/(f1-f0-fp0)
      if (lam1.lt.0.1D0) lam1=0.1D0
      call addvec(xnew,x,lam1,dx,n)
      call residuo_dpral(xnew,rnew,n,ndimx,nmodel,xtri,sig,vecm)
      f2=norm_l2_2(rnew,n)
      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
c posteriores line-search
      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,n)
      call residuo_dpral(xnew,rnew,n,ndimx,nmodel,xtri,sig,vecm)
      f3=norm_l2_2(rnew,n)
      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
c      if (escr)
      if (escr) write(uesb,*)' more line search',kmax
  80  continue
c final line-search
      err=err0*lam1/max(normainf(xnew,ndimx),ABS(xnew(ndimx+1)))
      call equv(x,xnew,n)
c muestra resultados iteraciones
c      if (nescri.eq.1) then
c        call der_enerelas_dpral(x,xnew,nmodel)
c        call Invari_t(xnew,3,cc)
c        call Invari_J2(xnew,3,bb)
c        call Invari_I1(xnew,3,aa)
c        write(ues,999)j,(x(i),i=1,n),err
c        write(ues,999)j,aa,bb,cc*45.D0/atan(1.D0),err,rnew(n)
999     format(I3,' ',5(E14.8,1x))
c      endif
c control de convergencia
      if (err.lt.tolrel) goto 99
      if (j.lt.jmax) goto 10
      kerre=1
  99  continue
      call equv(xini,x,ndimx)
      lam=x(n)
 
c muestra resultados finales
c      if (nescri.eq.2) 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