Esope source files

   1 : C JAC0RH    SOURCE    CHAT      05/01/13    00:47:56     5004           
   2 : CCC
   3 : C **********************************************************************
   4 : CCC
   5 :       SUBROUTINE JAC0RHMC (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 : *      write(6,*) ' dans jac0rhmc ',ndimk,x,ndims,ndimv,lam,nescri,ues,
  13 : *     >  kerre
  14 :       kerre=0
  15 :       void =0.D0
  16 :         do i=1,6
  17 :          Avecm(i)=0.D0
  18 :          vecntA(i)=0.D0
  19 :          vecn(i)=0.D0
  20 :          vecm(i)=0.D0
  21 :          vecEm(i)=0.D0
  22 :         enddo
  23 :         do i=1,36
  24 :          amat(i)=0.D0
  25 :          gmat(i)=0.D0
  26 :          hmat(i)=0.D0
  27 :         enddo
  28 : c m
  28 :       call vflsigRHMC(x,ndims,void,ndimv,vecm)
  29 : c n
  29 :       call vyisigRHMC(x,ndims,void,ndimv,vecn)
  30 : c H
  30 :       call HessFlsigRHMC(x,ndims,void,ndimv,Hmat,ndims)
  31 : c G=E
  31 :       call MatGenHook(Gmat,ndims,ndims)
  32 : *      write (6,*) ' jac0 1'
  33 : c Em = E*m
  34 :       do i=1,ndims
  35 :        do j=1,ndims
  36 :          vecEm(i)=vecEm(i)+Gmat(ndims*(j-1)+i)*vecm(j)
  37 :        enddo
  38 :       enddo
  39 : *      write (6,*) ' jac0 1a'
  40 : c A = I+l*E*H
  41 :       do i=1,ndims
  42 :        do j=1,ndims
  43 :           i1=ndims*(j-1)+i
  44 : 
  45 :         do k=1,ndims 
  46 :           Amat(i1)=Amat(i1)+
  47 :      .                        Gmat(ndims*(k-1)+i)*Hmat(ndims*(j-1)+k)
  48 :         enddo
  49 :         Amat(i1)=lam*Amat(i1)
  50 :        enddo
  51 :        Amat(ndims*(i-1)+i)=1.D0+Amat(ndims*(i-1)+i)
  52 :       enddo
  53 : c G = A-1
  54 : *      write (6,*) ' jac0 2'
  54 :        call DescLU(Amat,ndims)
  54 :       call LUinv(Amat,Gmat,ndims)
  55 : c g = nt*A-1*Em
  56 :       g=0.D0
  57 :       do i=1,ndims
  58 :        do j=1,ndims
  59 :           g=g+vecn(i)*Gmat(i+ndims*(j-1))*vecEm(j)
  60 :        enddo
  61 :       enddo
  62 : c G*Em; nt*G
  63 :       do i=1,ndims
  64 :        do j=1,ndims
  65 :          Avecm(i)=Avecm(i)+Gmat(i+ndims*(j-1))*vecEm(j)
  66 :          vecntA(i)=vecntA(i)+vecn(j)*Gmat(j+ndims*(i-1))
  67 :        enddo
  68 :       enddo
  69 : *           write (6,*) ' jac0 3'
  70 : c K = G-(G*m * nt*G)/g
  71 :       do i=1,ndims
  72 :        do j=1,ndims
  73 :           Kmat(i,j)=Gmat((j-1)*ndims+i)-Avecm(i)*vecntA(j)/g
  74 :        enddo
  75 :       enddo
  76 :       return
  77 :       end
  78 :  
  79 :  
  80 :