C DI2LEV    SOURCE    PV        22/04/19    16:18:01     11344          

      SUBROUTINE DIRECTION_2LEVELS(X,R,DX,N,NDIMX,NMODEL,SIG,VECM,
     .                             nnumer,deltax,lam)
      IMPLICIT INTEGER(I-N)
      integer n,i,j,k,ndims,ndimv,nmodel,nnumer,ndimx
      real*8  r(n),x(n),dx(n),deltax,sig(3),vecm(4)
      real*8  vecnaux(4),vecn(4),amat(16),vtLUiw,void(1)
      real*8  aamat(16),eemat(16),lam
      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 HessFlsig(sig,ndims,void,ndimv,amat,ndimx,nmodel)
      else if (ndimx.eq.4) then
       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 AA=A*E
      do i=1,ndimx
       do j=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
**
      if ((augla.eq.1)) then
c residuo y v flujo
      if (ndimx.eq.3) then
       call yieldd(sig,ndims,void,ndimv,res,nmodel)
       call vyisig(sig,ndims,void,ndimv,vecn,nmodel)
      else if (ndimx.eq.4) then
       call yieldd(sig,ndims,x(ndimx),ndimv,res,nmodel)
       call vyisig(sig,ndims,x(ndimx),ndimv,vecn,nmodel)
       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)+
     .                         vecn(i)*vecn(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
c         amat(i+(j-1)*ndimx)=amat(i+(j-1)*ndimx)+
c     .                       lam*aamat(i+(j-1)*ndimx)
**
         amat(i+(j-1)*ndimx)=amat(i+(j-1)*ndimx)+
     .                       ABS(lam+c*res)*aamat(i+(j-1)*ndimx)+
     .                       c*bbmat(i+(j-1)*ndimx)
**
       enddo
      enddo
      call DescLU(Amat,ndimx)
      call LUiw(Amat,r,dx,ndimx)
      do i=1,ndimx
        dx(i)=-1.D0*dx(i)
      enddo
      return
      end