Esope source files

   1 : C DI2LEV    SOURCE    CHAT      05/01/12    22:49:31     5004           
   2 : 
   3 :       SUBROUTINE DIRECTION_2LEVELS(X,R,DX,N,NDIMX,NMODEL,SIG,VECM,
   4 :      .                             nnumer,deltax,lam)
   5 :       IMPLICIT INTEGER(I-N)
   6 :       integer n,i,j,k,ndims,ndimv,nmodel,nnumer,ndimx
   7 :       real*8  r(n),x(n),dx(n),deltax,sig(3),vecm(4)
   8 :       real*8  vecnaux(4),vecn(4),amat(16),vtLUiw,void
   9 :       real*8  aamat(16),eemat(16),lam
  10 :       integer             augla
  11 :       real*8                    c
  12 :       common /auglagrang1/ augla
  13 :       common /auglagrang2/ c
  14 :       real*8  bbmat(16),res
  14 :       call zzero(amat ,16)
  14 :       call zzero(aamat,16)
  14 :       call zzero(bbmat,16)
  14 :       call zzero(eemat,16)
  15 :       void=0.D0
  16 :       ndimv=1
  17 :       ndims=3
  18 : c n
  19 : c A
  20 :       if (ndimx.eq.3) then
  20 :        call HessFlsig(sig,ndims,void,ndimv,amat,ndimx,nmodel)
  21 :       else if (ndimx.eq.4) then 
  21 :        call HessFlsig(sig,ndims,x(ndimx),ndimv,amat,ndimx,nmodel)
  22 :       endif
  23 : c E=d2_ener (ampliada de 3 a ndimx con 1 en la diagonal)
  23 :       call der2_enerelas_dpral(x,eemat,ndimx,nmodel)
  24 : c AA=A*E
  25 :       do i=1,ndimx
  26 :        do j=1,ndimx
  27 :          aamat(i+(j-1)*ndimx)=0.D0
  28 :          do k=1,ndimx
  29 :           aamat(i+(j-1)*ndimx)=aamat(i+(j-1)*ndimx)+
  30 :      .                         amat(i+(k-1)*ndimx)*eemat(k+(j-1)*ndimx)
  31 :          enddo 
  32 :        enddo
  33 :       enddo
  34 : **
  35 :       if ((augla.eq.1)) then 
  36 : c residuo y v flujo
  37 :       if (ndimx.eq.3) then
  37 :        call yieldd(sig,ndims,void,ndimv,res,nmodel)
  37 :        call vyisig(sig,ndims,void,ndimv,vecn,nmodel)
  38 :       else if (ndimx.eq.4) then 
  38 :        call yieldd(sig,ndims,x(ndimx),ndimv,res,nmodel)
  38 :        call vyisig(sig,ndims,x(ndimx),ndimv,vecn,nmodel)
  38 :        call vyivar(sig,ndims,x(ndimx),ndimv,vecn(ndimx),nmodel)
  39 :       endif
  40 : c BB=n*n^T*E
  41 :       do i=1,ndimx
  42 :        do j=1,ndimx
  43 :          bbmat(i+(j-1)*ndimx)=0.D0
  44 :            do k=1,ndimx
  45 :             bbmat(i+(j-1)*ndimx)=bbmat(i+(j-1)*ndimx)+
  46 :      .                         vecn(i)*vecn(k)*eemat(k+(j-1)*ndimx)
  47 :            enddo 
  48 :        enddo
  49 :       enddo
  50 :       endif
  51 : **
  51 :       call zzero(amat,16 )
  52 : c A=I+l*AA
  53 :       do i=1,ndimx
  54 :        amat(i+(i-1)*ndimx)=1.D0
  55 :        do j=1,ndimx
  56 : c         amat(i+(j-1)*ndimx)=amat(i+(j-1)*ndimx)+
  57 : c     .                       lam*aamat(i+(j-1)*ndimx)
  58 : **
  59 :          amat(i+(j-1)*ndimx)=amat(i+(j-1)*ndimx)+
  60 :      .                       dabs(lam+c*res)*aamat(i+(j-1)*ndimx)+
  61 :      .                       c*bbmat(i+(j-1)*ndimx)
  62 : **
  63 :        enddo
  64 :       enddo
  64 :       call DescLU(Amat,ndimx)
  64 :       call LUiw(Amat,r,dx,ndimx)
  65 :       do i=1,ndimx
  66 :         dx(i)=-1.D0*dx(i)
  67 :       enddo
  68 :       return
  69 :       end
  70 :  
  71 :  
  72 :