C MTCLSD    SOURCE    PV        22/04/19    16:18:08     11344          
CCC
C **********************************************************************
CCC
      SUBROUTINE MTC_LS_DPRAL(EKMAT,NDIMK,X,NDIMX,LAM,
     .                        nmodel,nescri,ues,nnumer,deltax)

      IMPLICIT INTEGER(I-N)
      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(1),sig(3)

      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)
      void(1)=0.D0
      ndimv=1
      ndims=3
      call der_enerelas_dpral(x,sig,nmodel)
c m,n,A
      if (ndimx.eq.3) then
       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 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
c EK= E*K
      call zzero(ekmat,ndimk**2)
      do i=1,3
       do j=1,3
         ekmat(i,j)=0.D0
         do k=1,ndimx
          ekmat(i,j)=ekmat(i,j)+
     .               eemat(i+(k-1)*ndimx)*kmat(k,j)
         enddo
       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