C MTCLDD    SOURCE    PV        07/11/23    21:17:54     5978           
CCC
C **********************************************************************
CCC
      SUBROUTINE MTC_LS_DPRAL_DENSI(EKMAT,NDIMK,X,NDIMX,LAM,
     .                              xdensi,xmat,
     .                              nmodel,nescri,ues,nnumer,deltax)

      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (a-h,o-z)
      integer ndimk,ndimx,nmodel,nescri,ues,nnumer
      integer i,j,k,ndims,ndimv
      real*8  EKmat(ndimk,ndimk),x(ndimx),lam,deltax
      real*8  amat(16),aamat(16),eemat(16),vecm(4),vecnaux(4),vecn(4),
     .        gmat(16),Avecm(4),vecntA(4),kmat(4,4),g,void,sig(3)
      real*8  xdensi,dxdensi,xdensi2,xmat(*),resu,resuaux
      real*8  Mmat(4,3),kmat2(4,4),kmat3(4,4),ekmat2(3,3),ekmat3(3,3),
     .        vecmaux(4)

      call zzero(amat ,16)
      call zzero(aamat,16)
      call zzero(eemat,16)
      call zzero(gmat ,16)
      call zzero(kmat ,16)
      call zzero(vecm,4)
      call zzero(vecn,4)
      call zzero(vecnaux,4)
      call zzero(vecn,4)
      call zzero(Avecm,4)
      call zzero(vecntA,4)
      call zzero(sig,3)

      call zzero(mmat,12)
      call zzero(kmat2,16)
      call zzero(kmat3,16)
      call zzero(ekmat2,9)
      call zzero(ekmat3,9)
      call zzero(vecmaux,4)

      void=0.D0
      ndimv=1
      ndims=3
      call der_enerelas_dpral(x,sig,nmodel)
c m,n,A
      if (ndimx.eq.3) then
       call yieldd(sig,ndims,void,ndimv,resu,nmodel)
       call vflsig(sig,ndims,void,ndimv,vecm,nmodel)
       call vyisig(sig,ndims,void,ndimv,vecnaux,nmodel)
       call HessFlsig(sig,ndims,void,ndimv,amat,ndimx,nmodel)
      else if (ndimx.eq.4) then 
       call yieldd(sig,ndims,x(ndimx),ndimv,resu,nmodel)
       call vflsig(sig,ndims,x(ndimx),ndimv,vecm,nmodel)
       call vflvar(sig,ndims,x(ndimx),ndimv,vecm(ndimx),nmodel)
       call vyisig(sig,ndims,x(ndimx),ndimv,vecnaux,nmodel)
       call vyivar(sig,ndims,x(ndimx),ndimv,vecnaux(ndimx),nmodel)
       call HessFlsig(sig,ndims,x(ndimx),ndimv,amat,ndimx,nmodel)
      endif
c E=d2_ener (ampliada de 3 a ndimx con 1 en la diagonal)
      call der2_enerelas_dpral(x,eemat,ndimx,nmodel)
c n^T=n^T*E
c AA=A*E
      do i=1,ndimx
       vecn(i)=0.D0
       do j=1,ndimx
         vecn(i)=vecn(i)+vecnaux(j)*eemat(j+(i-1)*ndimx)
         aamat(i+(j-1)*ndimx)=0.D0
         do k=1,ndimx
           aamat(i+(j-1)*ndimx)=aamat(i+(j-1)*ndimx)+
     .                         amat(i+(k-1)*ndimx)*eemat(k+(j-1)*ndimx)
         enddo 
       enddo
      enddo
      call zzero(amat,16 )
c A=I+l*AA
      do i=1,ndimx
       amat(i+(i-1)*ndimx)=1.D0
       do j=1,ndimx
         amat(i+(j-1)*ndimx)=amat(i+(j-1)*ndimx)+
     .                       lam*aamat(i+(j-1)*ndimx)
       enddo
      enddo
c G = A-1
      call DescLU(Amat,ndimx)
      call LUinv(Amat,Gmat,ndimx)
c g = nt*A-1*m
      g=0.D0
      do i=1,ndimx
       do j=1,ndimx
          g=g+vecn(i)*Gmat(i+ndimx*(j-1))*vecm(j)
       enddo
      enddo
c G*m; nt*G
      do i=1,ndimx
       do j=1,ndimx
         Avecm(i)=Avecm(i)+Gmat(i+ndimx*(j-1))*vecm(j)
         vecntA(i)=vecntA(i)+vecn(j)*Gmat(j+ndimx*(i-1))
       enddo
      enddo
c K = G-(G*m * nt*G)/g
      do i=1,ndimx
       do j=1,ndimx
          Kmat(i,j)=Gmat((j-1)*ndimx+i)-Avecm(i)*vecntA(j)/g
       enddo
      enddo
CCCC 
C densi
CCCC
c Vecm y f perturbados

c      write(399,*)
c      dxrel = 0.1
c      do i = 1,10
c      dx = xdensi * (dxrel**i)
c      xdensi2 = xdensi + dx
c      call carac_mate_rhmc_densi(XMAT,xdensi2)
c      call yieldd(sig,ndims,void,ndimv,resuaux,nmodel)
c      call vflsig(sig,ndims,void,ndimv,vecmaux,nmodel)
c      write(399,*)dxrel**i,(resuaux-resu)/dx
c      write(399,*)'  ',((vecmaux(j)-vecm(j))/dx,j=1,3)
c      call carac_mate_rhmc_densi(XMAT,xdensi)
c      enddo
c      do i = 1,10
c      dx = - xdensi * (dxrel**i)
c      xdensi2 = xdensi + dx
c      call carac_mate_rhmc_densi(XMAT,xdensi2)
c      call yieldd(sig,ndims,void,ndimv,resuaux,nmodel)
c      call vflsig(sig,ndims,void,ndimv,vecmaux,nmodel)
c      write(399,*)-dxrel**i,(resuaux-resu)/dx
c      write(399,*)'  ',((vecmaux(j)-vecm(j))/dx,j=1,3)
c      call carac_mate_rhmc_densi(XMAT,xdensi)
c      enddo
CC
      dxrel = 0.1D0
      dx = (1.D0 + xdensi) * (dxrel**6)
      xdensi2 = xdensi + dx
c      call carac_mate_rhmc_densi(XMAT,xdensi2)
      call carac_mate_densi(XMAT,xdensi2,nmodel)
      if (ndimx.eq.3) then
       call yieldd(sig,ndims,void,ndimv,resuaux,nmodel)
       call vflsig(sig,ndims,void,ndimv,vecmaux,nmodel)
      else if (ndimx.eq.4) then 
       call yieldd(sig,ndims,x(ndimx),ndimv,resuaux,nmodel)
       call vflsig(sig,ndims,x(ndimx),ndimv,vecmaux,nmodel)
       call vflvar(sig,ndims,x(ndimx),ndimv,vecmaux(ndimx),nmodel)
      endif
      xdensi2 = xdensi - dx
c      call carac_mate_rhmc_densi(XMAT,xdensi2)
      call carac_mate_densi(XMAT,xdensi2,nmodel)
      if (ndimx.eq.3) then
       call yieldd(sig,ndims,void,ndimv,resu,nmodel)
       call vflsig(sig,ndims,void,ndimv,vecm,nmodel) 
      else if (ndimx.eq.4) then 
       call yieldd(sig,ndims,x(ndimx),ndimv,resu,nmodel)
       call vflsig(sig,ndims,x(ndimx),ndimv,vecm,nmodel)
       call vflvar(sig,ndims,x(ndimx),ndimv,vecm(ndimx),nmodel)
      endif
      dx = dx * 2.D0
c      call carac_mate_rhmc_densi(XMAT,xdensi)
      call carac_mate_densi(XMAT,xdensi,nmodel)
c Mmat = derivada de vecm (3 veces)
      do i=1,ndimx
          Mmat(i,1)=(vecmaux(i)-vecm(i))/dx
          Mmat(i,2)=Mmat(i,1)
          Mmat(i,3)=Mmat(i,1)
      enddo
c Kmat2 = Kmat * Mmat
c Kmat3 = Avecm/g (3 veces)
      do i=1,ndimx
       do j=1,3
         do k=1,ndimx
          Kmat2(i,j)=Kmat2(i,j)+Kmat(i,k)*Mmat(k,j)
         enddo
         Kmat3(i,j)=Avecm(i)/g
       enddo
      enddo
c auxiliares
      aux1= xdensi*lam
      aux2= xdensi*(resuaux-resu)/dx
c EKmat  = E*Kmat ; EKmat2 = E*Kmat2 ; EKmat3 = E*Kmat3
c EKmat  = EKmat + aux1*EKmat2 + aux2*EKmat3
      call zzero(ekmat,ndimk**2)
      do i=1,3
       do j=1,3
         do k=1,ndimx
          ekmat(i,j) =ekmat(i,j)+
     .                eemat(i+(k-1)*ndimx)*kmat(k,j)
          ekmat2(i,j)=ekmat2(i,j)+
     .                eemat(i+(k-1)*ndimx)*kmat2(k,j)
          ekmat3(i,j)=ekmat3(i,j)+
     .                eemat(i+(k-1)*ndimx)*kmat3(k,j)
         enddo 
         ekmat(i,j)=ekmat(i,j)+aux1*ekmat2(i,j)+aux2*ekmat3(i,j)
       enddo
      enddo

      if (nescri.eq.1)
     .   write(ues,999)' EKmat_dp ',((ekmat(i,j),j=1,3),i=1,3)
999   format(2x,A12,/,3(3(1x,E12.6),/))
      return
      end