Esope source files

   1 : C ITG0J2    SOURCE    CHAT      05/01/13    00:44:19     5004           
   2 : CCC
   3 : C **********************************************************************
   4 : CCC
   5 :       SUBROUTINE INTEGRA0J2 (XTRI,X,NDIMX,LAM,DDEFPL,NDIMS,
   6 :      .                       tolrel,nitmax,nescri,ues,kerre,iiter)
   7 :       IMPLICIT INTEGER(I-N)
   8 :       integer ndims,nitmax,nescri,ndimx,ues
   9 :       real*8  xtri(ndimx),x(ndimx),ddefpl(ndims),lam,tolrel
  10 :       integer i,j,iiter,kerre,ndimv
  11 :       real*8  dl,lres,siginv(3),void(1),auxr1,vtLUiw,auxmax1,auxmax2
  12 :       real*8  xres(6),dx(6),Amat(36),Gmat(36),vecn(6),vecm(6)
  13 :       ndimv=1
  14 :       kerre=0
  15 :       void(1)=0.D0
  16 :         do i=1,6
  17 :          xres(i)=0.D0
  18 :          dx(i)=0.D0
  19 :          vecn(i)=0.D0
  20 :          vecm(i)=0.D0
  21 :         enddo
  22 :         do i=1,36
  23 :          amat(i)=0.D0
  24 :          gmat(i)=0.D0
  25 :         enddo
  26 :       iiter=-1      
  26 :       call MatGenHookinv(Gmat,ndimx,ndims)
  27 :         dl=1.D0
  28 : 10    continue
  29 :       iiter=iiter+1
  29 :       call yielddJ2(x,ndims,void,ndimv,lres)
  29 :       call vflsigJ2(x,ndims,void,ndimv,vecm)
  30 :       do i=1,ndimx
  31 :          xres(i)=lam*vecm(i)
  32 :          do j=1,ndimx
  33 :             xres(i)=xres(i)+Gmat((i-1)*ndimx+j)*(x(j)-xtri(j))
  34 :          enddo
  35 :       enddo
  36 :       auxmax1=dabs(dx(1))
  37 :       auxmax2=dabs(x(1))
  38 :       do i=2,ndimx
  39 :        if (dabs(dx(i)).gt.auxmax1) auxmax1=dabs(dx(i))
  40 :        if (dabs(x(i)).gt.auxmax2) auxmax2=dabs(x(i))
  41 :       enddo
  42 :       auxr1=max(auxmax1,dabs(dl))/max(auxmax2,1.D0)
  43 :       if (nescri.eq.1) then
  44 :         write(ues,'(I5,3X,E12.6)')iiter,auxr1
  45 :       endif
  46 :       if (auxr1.lt.tolrel) then
  47 : c CONVERGIDO
  48 :          do i=1,ndims
  49 :             ddefpl(i)=0.D0
  50 :             do j=1,ndims
  51 :                ddefpl(i)=ddefpl(i)+Gmat((i-1)*ndimx+j)*(xtri(j)-x(j))
  52 :             enddo
  53 :          enddo
  54 :          return
  55 :       endif
  56 : c NO CONVERGIDO
  57 :       if (iiter.eq.nitmax) then
  58 :          kerre=1
  59 :          return
  60 :       endif
  60 :       call vyisigJ2(x,ndims,void,ndimv,vecn)
  60 :       call HessFlsigJ2(x,ndimx,void,ndimv,Amat,ndimx)
  61 :       do i=1,ndimx*ndimx
  62 :          Amat(i)=Gmat(i)+lam*Amat(i)
  63 :       enddo
  63 :       call DescLU(Amat,ndimx)
  64 :       dl=(lres-vtLUiw(vecn,Amat,xres,ndimx))
  65 :      .         /vtLUiw(vecn,Amat,vecm,ndimx)
  66 :       do i=1,ndimx
  67 :          xres(i)=-xres(i)-dl*vecm(i)
  68 :       enddo
  68 :       call LUiw(Amat,xres,dx,ndimx)
  69 :       lam=lam+dl
  70 :       do i=1,ndimx
  71 :          x(i)=x(i)+dx(i)
  72 :       enddo
  73 :       go to 10
  74 :       end
  75 :  
  76 :  
  77 :