Télécharger la plaquette Cast3M


fr en

Télécharger didpra.eso

Retour à la liste

Numérotation des lignes : 

didpra

  1. C DIDPRA    SOURCE    PV        22/04/19    16:18:02     11344          
  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(1)
  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.  
  16.       call zzero(amat ,16)
  17.       call zzero(aamat,16)
  18.       call zzero(bbmat,16)
  19.       call zzero(eemat,16)
  20.  
  21.       void(1)=0.D0
  22.       ndimv=1
  23.       ndims=3
  24. c n
  25. c A
  26.       if (ndimx.eq.3) then
  27.        call vyisig (sig,ndims,void,ndimv,vecnaux,nmodel)
  28.        call HessFlsig (sig,ndims,void,ndimv,amat,ndimx,nmodel)
  29.       else if (ndimx.eq.4) then
  30.        call vyisig (sig,ndims,x(ndimx),ndimv,vecnaux,nmodel)
  31.        call vyivar (sig,ndims,x(ndimx),ndimv,vecnaux(ndimx),nmodel)
  32.        call HessFlsig (sig,ndims,x(ndimx),ndimv,amat,ndimx,nmodel)
  33.       endif
  34. c E=d2_ener (ampliada de 3 a ndimx con 1 en la diagonal)
  35.       call der2_enerelas_dpral(x,eemat,ndimx,nmodel)
  36.  
  37. c n^T=n^T*E
  38. c AA=A*E
  39.       do i=1,ndimx
  40.        vecn(i)=0.D0
  41.        do j=1,ndimx
  42.          vecn(i)=vecn(i)+vecnaux(j)*eemat(j+(i-1)*ndimx)
  43.          aamat(i+(j-1)*ndimx)=0.D0
  44.          i2=i+(j-1)*ndimx
  45.          do k=1,ndimx
  46.            aamat(i2)=aamat(i2)+
  47.      .                         amat(i+(k-1)*ndimx)*eemat(k+(j-1)*ndimx)
  48.          enddo
  49.        enddo
  50.       enddo
  51.       res=0.d0
  52.       if ((augla.eq.1)) then
  53. c residuo y v flujo
  54.       if (ndimx.eq.3) then
  55.        call yieldd(sig,ndims,void,ndimv,res,nmodel)
  56. c       call vyisig(sig,ndims,void,ndimv,vecn,nmodel)
  57.       else if (ndimx.eq.4) then
  58.        call yieldd(sig,ndims,x(ndimx),ndimv,res,nmodel)
  59. c       call vyisig(sig,ndims,x(ndimx),ndimv,vecn,nmodel)
  60. c       call vyivar(sig,ndims,x(ndimx),ndimv,vecn(ndimx),nmodel)
  61.       endif
  62. c BB=n*n^T*E
  63.       do i=1,ndimx
  64.        do j=1,ndimx
  65.          bbmat(i+(j-1)*ndimx)=0.D0
  66.          do k=1,ndimx
  67.            bbmat(i+(j-1)*ndimx)=bbmat(i+(j-1)*ndimx)+
  68.      .                         vecm(i)*vecm(k)*eemat(k+(j-1)*ndimx)
  69.          enddo
  70.        enddo
  71.       enddo
  72.       endif
  73.       call zzero(amat,16 )
  74. c A=I+l*AA
  75.       do i=1,ndimx
  76.        amat(i+(i-1)*ndimx)=1.D0
  77.        do j=1,ndimx
  78.          amat(i+(j-1)*ndimx)=amat(i+(j-1)*ndimx)+
  79.      .                       ABS(x(n)+c*res)*aamat(i+(j-1)*ndimx)+
  80.      .                       c*bbmat(i+(j-1)*ndimx)
  81.        enddo
  82.       enddo
  83.       call DescLU(Amat,ndimx)
  84.       dx(n)=(r(n)-vtLUiw(vecn,Amat,r,ndimx))
  85.      .      /vtLUiw(vecn,Amat,vecm,ndimx)
  86.       do i=1,ndimx
  87.          r(i)=-r(i)-dx(n)*vecm(i)
  88.       enddo
  89.       call LUiw(Amat,r,dx,ndimx)
  90.       return
  91.       end
  92.  
  93.  
  94.  
  95.  
  96.  
  97.  
  98.  
  99.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales