Esope source files

   1 : C JAC0J2    SOURCE    CHAT      05/01/13    00:47:53     5004           
   2 : CCC
   3 : C **********************************************************************
   4 : CCC
   5 :       SUBROUTINE JAC0J2 (KMAT,NDIMK,X,NDIMS,NDIMV,LAM,
   6 :      .                   nescri,ues,kerre)
   7 :       IMPLICIT INTEGER(I-N)
   8 :       integer ndims,nescri,ues,ndimv,i,j,k,kerre,ndimk
   9 :       real*8  x(ndims),lam,void,Kmat(ndimk,ndimk),g
  10 :       real*8  Amat(36),Gmat(36),vecn(6),vecm(6),Avecm(6),vecntA(6),
  11 :      .        vecEm(6),Hmat(36)
  12 :       kerre=0
  13 :       void =0.D0
  14 :         do i=1,6
  15 :          Avecm(i)=0.D0
  16 :          vecntA(i)=0.D0
  17 :          vecn(i)=0.D0
  18 :          vecm(i)=0.D0
  19 :          vecEm(i)=0.D0
  20 :         enddo
  21 :         do i=1,36
  22 :          amat(i)=0.D0
  23 :          gmat(i)=0.D0
  24 :          hmat(i)=0.D0
  25 :         enddo
  26 : c m
  26 :       call vflsigJ2(x,ndims,void,ndimv,vecm)
  27 : c n
  27 :       call vyisigJ2(x,ndims,void,ndimv,vecn)
  28 : c H
  28 :       call HessFlsigJ2(x,ndims,void,ndimv,Hmat,ndims)
  29 : c G=E
  29 :       call MatGenHook(Gmat,ndims,ndims)
  30 : c Em = E*m
  31 :       do i=1,ndims
  32 :        do j=1,ndims
  33 :          vecEm(i)=vecEm(i)+Gmat(ndims*(j-1)+i)*vecm(j)
  34 :        enddo
  35 :       enddo
  36 : c A = I+l*E*H
  37 :       do i=1,ndims
  38 :        do j=1,ndims
  39 :         do k=1,ndims
  40 :           Amat(ndims*(j-1)+i)=Amat(ndims*(j-1)+i)+
  41 :      .                        Gmat(ndims*(k-1)+i)*Hmat(ndims*(j-1)+k)
  42 :         enddo
  43 :         Amat(ndims*(j-1)+i)=lam*Amat(ndims*(j-1)+i)
  44 :        enddo
  45 :        Amat(ndims*(i-1)+i)=1.D0+Amat(ndims*(i-1)+i)
  46 :       enddo
  47 : c G = A-1
  47 :       call DescLU(Amat,ndims)
  47 :       call LUinv(Amat,Gmat,ndims)
  48 : c g = nt*A-1*Em
  49 :       g=0.D0
  50 :       do i=1,ndims
  51 :        do j=1,ndims
  52 :           g=g+vecn(i)*Gmat(i+ndims*(j-1))*vecEm(j)
  53 :        enddo
  54 :       enddo
  55 : c G*Em; nt*G
  56 :       do i=1,ndims
  57 :        do j=1,ndims
  58 :          Avecm(i)=Avecm(i)+Gmat(i+ndims*(j-1))*vecEm(j)
  59 :          vecntA(i)=vecntA(i)+vecn(j)*Gmat(j+ndims*(i-1))
  60 :        enddo
  61 :       enddo
  62 : c K = G-(G*m * nt*G)/g
  63 :       do i=1,ndims
  64 :        do j=1,ndims
  65 :           Kmat(i,j)=Gmat((j-1)*ndims+i)-Avecm(i)*vecntA(j)/g
  66 :        enddo
  67 :       enddo
  68 :       return
  69 :       end
  70 :  
  71 :