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