Télécharger itgl2l.eso

Retour à la liste

Numérotation des lignes :

  1. C ITGL2L SOURCE CHAT 05/01/13 00:44:39 5004
  2. CCC
  3. C **********************************************************************
  4. CCC
  5. SUBROUTINE INTEGRA_LS_2LEVELS(XTRI,XINI,NDIMX,LAM,NMODEL,TOLREL,
  6. . jmax,nescri,ues,nnumer,deltax,kerre,jtot)
  7.  
  8. IMPLICIT INTEGER(I-N)
  9. integer nmodel,nitmax,nescri,nnumer,ues,j,kerre,
  10. . ndimv,jmax,i,ii,k,n,ndimx,uesb,jl,jtot
  11. real*8 xtri(ndimx),xini(ndimx),lam,tolrel,deltax,
  12. . x(5),dx(5),r(5),err,normainf,norm_l2_2,xnew(5),rnew(5),
  13. . cotaere,lam0,lam1,lam2,f0,fp0,f1,f2,f3,s1,s2,aa,bb,cc,
  14. . sig(4),vecm(4),err0,dla,dlb,ddll,beta
  15. logical escr,last
  16. integer level,kmax
  17. common /linesearch/ level,kmax
  18. integer augla
  19. real*8 c
  20. common /auglagrang1/ augla
  21. common /auglagrang2/ c
  22.  
  23. c muestra tensiones
  24. c if (nescri.eq.1) then
  25. c call der_enerelas_dpral(xini,x,nmodel)
  26. c write(ues,*)' Tensiones iniciales'
  27. c write(ues,*)(x(i),i=1,3)
  28. c if (ndimx.eq.4) write(ues,*)' Variable interna inicial',x(4)
  29. c endif
  30.  
  31. c muestra iteraciones de line-search
  32. uesb=299
  33. escr=.false.
  34. if (nescri.eq.1) escr=.true.
  35.  
  36. c parametros generales
  37. kerre=0
  38. n=ndimx+1
  39. ndimv=1
  40. call equv(x,xini,ndimx)
  41. x(n)=lam
  42. beta=1.D-4
  43. j=0
  44.  
  45. jtot=-1
  46. last=.false.
  47. jl=0
  48.  
  49. C CALCULOS PARA LAMBDA; SUPONEMOS SIG_INI=SIG_TRIAL
  50. 2 continue
  51. if ((augla.eq.1).and.(jl.eq.0)) then
  52. jl=1
  53. goto 10
  54. endif
  55. jl=jl+1
  56. C A; B
  57. call actualiza_lambda(x,ndimx,nmodel,lam,dla,dlb)
  58. ddll=-dla/dlb
  59. lam=lam+ddll
  60. if (lam.lt.0) then
  61. write(ues,*)' GAMMA negative !',lam,' => ',0.D0
  62. write(ues,*)' dll:',ddll
  63. lam=0
  64. endif
  65. C CHECK DL => GOTO END
  66. jtot=jtot+1
  67. c if (nescri.eq.1) write(ues,999)jtot,jl,lam,ddll,dla
  68. if (last) goto 100
  69. if (abs(ddll).lt.tolrel) last=.true.
  70. if (jl.gt.50+1) then
  71. write(ues,*)' ERROR: loop gamma no ha convergido en', 50
  72. kerre=1
  73. goto 100
  74. endif
  75.  
  76. C COMPUTE SIGMA_J FROM LAMBDA_J WITH LINE-SEARCH
  77. 10 continue
  78. j=j+1
  79. call residuo_2levels(x,r,ndimx,ndimx,nmodel,xtri,sig,vecm,lam)
  80. f0=norm_l2_2(r,ndimx)
  81. fp0=-2.D0*f0
  82. call direction_2levels
  83. . (x,r,dx,ndimx,ndimx,nmodel,sig,vecm,nnumer,deltax,lam)
  84. err0=normainf(dx,ndimx)
  85. lam0=1.D0
  86. call addvec(xnew,x,lam0,dx,ndimx)
  87. call residuo_2levels
  88. . (xnew,rnew,ndimx,ndimx,nmodel,xtri,sig,vecm,lam)
  89. f1=norm_l2_2(rnew,ndimx)
  90. cotaere=f0+beta*lam0*fp0
  91. lam1=lam0
  92. if (escr) write(uesb,*)' more ',0,lam0
  93. if (f1.lt.cotaere) goto 80
  94. if ((f0.lt.tolrel).and.(f1.lt.tolrel)) goto 80
  95. if (kmax.le.0) goto 70
  96.  
  97. k=1
  98. lam1=(-fp0)/2.D0/(f1-f0-fp0)
  99. if (lam1.lt.0.1D0) lam1=0.1D0
  100. call addvec(xnew,x,lam1,dx,ndimx)
  101. call residuo_2levels(xnew,rnew,n,ndimx,nmodel,xtri,sig,vecm,lam)
  102. f2=norm_l2_2(rnew,ndimx)
  103. cotaere=f0+beta*lam1*fp0
  104. if (escr) write(uesb,*)' more ',k,lam1
  105. if (f2.lt.cotaere) goto 80
  106. if (kmax.le.1) goto 70
  107.  
  108. 20 continue
  109. k=k+1
  110. lam2=(-fp0*lam1**2)/2.D0/(f2-f0-fp0*lam1)
  111. if (lam2.lt.(0.1D0*lam1)) lam2=0.1D0*lam1
  112. if (lam2.gt.(0.5D0*lam1)) lam2=0.5D0*lam1
  113. call addvec(xnew,x,lam2,dx,ndimx)
  114. call residuo_2levels(xnew,rnew,n,ndimx,nmodel,xtri,sig,vecm,lam)
  115. f3=norm_l2_2(rnew,ndimx)
  116. cotaere=f0+beta*lam2*fp0
  117. f1=f2
  118. f2=f3
  119. lam0=lam1
  120. lam1=lam2
  121. if (escr) write(uesb,*)' more ',k,lam2
  122. if (f3.lt.cotaere) goto 80
  123. if (k.lt.kmax) goto 20
  124. 70 continue
  125. if (escr) write(uesb,*)' more line search',kmax
  126. 80 continue
  127. err=err0*lam1/normainf(xnew,ndimx)
  128. call equv(x,xnew,ndimx)
  129. c MUESTRA ITERACIONES EN SIGMA_J
  130. c if (nescri.eq.1) then
  131. c call Invari_t(x,3,aa)
  132. c call Invari_q(x,3,bb)
  133. c write(ues,999)jtot,j,(x(i),i=1,3),lam,err
  134. c write(ues,999)j,aa,bb,(x(i),i=4,n),err
  135. c 999 format(2I4,' ',6(E10.4,1x))
  136. c call der_enerelas_dpral(x,xnew,nmodel)
  137. c call Invari_t(xnew,3,cc)
  138. c call Invari_q(xnew,3,bb)
  139. c call Invari_p(xnew,3,aa)
  140. c write(ues,1999)j,aa,bb,cc*45.D0/atan(1.D0),err,rnew(n)
  141. c1999 format(I3,' ',5(E14.8,1x))
  142. c endif
  143. if (err.lt.tolrel) goto 99
  144. if (j.lt.jmax) goto 10
  145. kerre=1
  146. write(ues,*)' ERROR: Loop interno sigmas no ha convergido'
  147. 99 continue
  148. c FIN ITERACIONES EN SIGMA_J
  149. jtot=jtot+j
  150. j=0
  151. goto 2
  152.  
  153. c RESULTADOS FINALES
  154. 100 continue
  155. call equv(xini,x,ndimx)
  156.  
  157. c MUESTRA RESULTADOS
  158. c if (nescri.eq.1) then
  159. c write(ues,*)' Deform elas finales'
  160. c write(ues,*)(x(i),i=1,3)
  161. c write(ues,*)' Deform plas finales'
  162. c write(ues,*)(xtri(i)-x(i),i=1,3)
  163. c if (ndimx.eq.4) write(ues,*)' Variable interna final',x(4)
  164. c call der_enerelas_dpral(x,xini,nmodel)
  165. c write(ues,*)' Tensiones finales'
  166. c write(ues,*)(xini(i),i=1,3)
  167. c if (ndimx.eq.4) write(ues,*)' Variable interna inicial',x(4)
  168. c call equv(xini,x,ndimx)
  169. c endif
  170. end
  171.  
  172.  
  173.  
  174.  
  175.  

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