Esope source files

   1 : C JAC21     SOURCE    CHAT      05/01/13    00:48:02     5004           
   2 : CCC
   3 : C **********************************************************************
   4 : CCC
   5 :       SUBROUTINE JAC21 (KMAT,NDIMK,X,NDIMX,LAM,NDIMS,
   6 :      .                nmodel,nescri,ues,nnumer,deltax,kerre)
   7 :       IMPLICIT INTEGER(I-N)
   8 :       integer ndims,nmodel,nescri,nnumer,ues,ndimv,ndimx
   9 :       integer i,j,k,kerre,ndimk
  10 :       real*8    x(ndimx),lam,void,Kmat(ndimk,ndimk),g,deltax
  11 :       real*8    Amat(64),Gmat(64),vecn(8),vecm(8),Avecm(8),vecntA(8),
  12 :      .          vecEm(8),Hmat(64)
  13 :       kerre=0
  14 :       ndimv=2
  15 :       void =0.D0
  16 :       do i=1,8
  17 :           vecntA(i)=0.D0
  18 :           Avecm(i)=0.D0
  19 :           vecn(i)=0.D0
  20 :           vecm(i)=0.D0
  21 :           vecEm(i)=0.D0
  22 :         enddo
  23 :       do i=1,64
  24 :           Amat(i)=0.D0
  25 :           Gmat(i)=0.D0
  26 :           Hmat(i)=0.D0
  27 :         enddo
  28 : c m
  28 :       call VecFlMAC(x,ndimx,vecm,nmodel)
  29 : c n
  29 :       call VecYiMAC(x,ndimx,vecn,nmodel)
  30 : c H
  30 :       call HessMAC(x,ndimx,Hmat,ndimx,nnumer,deltax,nmodel)
  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 :