Esope source files

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