C DIDPRA    SOURCE    PV        22/04/19    16:18:02     11344          
C
      SUBROUTINE DIRECTION_DPRAL(X,R,DX,N,NDIMX,NMODEL,SIG,VECM,
     .                           nnumer,deltax)
      IMPLICIT INTEGER(I-N)
      integer n,i,j,k,ndims,ndimv,nmodel,nnumer,ndimx
      real*8  r(n),x(n),dx(n),deltax,sig(4),vecm(4)
      real*8  vecnaux(4),vecn(4),amat(16),vtLUiw,void(1)
      real*8  aamat(16),eemat(16)
      integer             augla
      real*8                    c
      common /auglagrang1/ augla
      common /auglagrang2/ c
      real*8  bbmat(16),res

      call zzero(amat ,16)
      call zzero(aamat,16)
      call zzero(bbmat,16)
      call zzero(eemat,16)

      void(1)=0.D0
      ndimv=1
      ndims=3
c n
c A
      if (ndimx.eq.3) then
       call vyisig (sig,ndims,void,ndimv,vecnaux,nmodel)
       call HessFlsig (sig,ndims,void,ndimv,amat,ndimx,nmodel)
      else if (ndimx.eq.4) then
       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
         i2=i+(j-1)*ndimx
         do k=1,ndimx
           aamat(i2)=aamat(i2)+
     .                         amat(i+(k-1)*ndimx)*eemat(k+(j-1)*ndimx)
         enddo
       enddo
      enddo
      res=0.d0
      if ((augla.eq.1)) then
c residuo y v flujo
      if (ndimx.eq.3) then
       call yieldd(sig,ndims,void,ndimv,res,nmodel)
c       call vyisig(sig,ndims,void,ndimv,vecn,nmodel)
      else if (ndimx.eq.4) then
       call yieldd(sig,ndims,x(ndimx),ndimv,res,nmodel)
c       call vyisig(sig,ndims,x(ndimx),ndimv,vecn,nmodel)
c       call vyivar(sig,ndims,x(ndimx),ndimv,vecn(ndimx),nmodel)
      endif
c BB=n*n^T*E
      do i=1,ndimx
       do j=1,ndimx
         bbmat(i+(j-1)*ndimx)=0.D0
         do k=1,ndimx
           bbmat(i+(j-1)*ndimx)=bbmat(i+(j-1)*ndimx)+
     .                         vecm(i)*vecm(k)*eemat(k+(j-1)*ndimx)
         enddo
       enddo
      enddo
      endif
      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)+
     .                       ABS(x(n)+c*res)*aamat(i+(j-1)*ndimx)+
     .                       c*bbmat(i+(j-1)*ndimx)
       enddo
      enddo
      call DescLU(Amat,ndimx)
      dx(n)=(r(n)-vtLUiw(vecn,Amat,r,ndimx))
     .      /vtLUiw(vecn,Amat,vecm,ndimx)
      do i=1,ndimx
         r(i)=-r(i)-dx(n)*vecm(i)
      enddo
      call LUiw(Amat,r,dx,ndimx)
      return
      end