Esope source files

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