Télécharger simu21.eso

Retour à la liste

Numérotation des lignes :

  1. C SIMU21 SOURCE FANDEUR 10/12/14 21:19:41 6806
  2.  
  3. SUBROUTINE SIMU21(IPLA,IPLB,iplc,NMOD,CONV,nbonve,numv1,
  4. $ iplval,iphi,IPERM,XPREC21)
  5.  
  6. *********************************************************************
  7. *
  8. * S I M U 2 1
  9. * -----------
  10. *
  11. * FONCTION:
  12. * ---------
  13. * Passage au petit probleme (calcul de T et de D),
  14. * Resolution du petit probleme [T-(1/lambda)D]*Phi=0
  15. * Test de convergence
  16. *
  17. * PARAMETRES : (E)=ENTREE (S)=SORTIE
  18. * -----------
  19. * IPLA LISTREEL (E) alpha
  20. * IPLB LISTREEL (E) beta
  21. * iplc LISTREEL (E) sign(d)
  22. * NMOD ENTIER (E) nombre de modes demandés pour ce cycle
  23. * CONV LOGIQUE (S) = vrai si on a trouvé nbonve > NMOD
  24. * nbonve ENTIER (S) nombre de modes ayant convergés
  25. * numv1 ENTIER (E) dimension du probleme reduit
  26. * iplval ENTIER (S) POINTEUR VERS LES VALEURS PROPRES CONVERGEES
  27. * iphi ENTIER (S) POINTEUR VERS LES VECTEURS PROPRES CONVERGES
  28. * IPERM ENTIER (S) POINTEUR VERS LE TABLEAU DES PERMUTATIONS SI
  29. * LES MODES (iplval et iphi) SONT RANGES DANS
  30. * UN ORDRE DIFFERENT DE 1 2 3 ... nbonve
  31. * XPREC21 FLOTTANT (E) PRECISION SUR LES VALEURS PROPRES
  32. *
  33. *********************************************************************
  34.  
  35. IMPLICIT INTEGER(I-N)
  36. IMPLICIT REAL*8 (A-H,O-Z)
  37.  
  38. -INC CCOPTIO
  39. -INC SMLREEL
  40. -INC SMLENTI
  41. -INC SMLCHPO
  42.  
  43. *=======================================================*
  44. * TYPE MATRIX et DECLARATIONS
  45. *=======================================================*
  46.  
  47. SEGMENT MATRIX
  48. REAL*8 A(N,N)
  49. ENDSEGMENT
  50.  
  51. POINTEUR IPRIG1.MATRIX, IPMAS1.MATRIX, IPHI.MATRIX
  52. POINTEUR IPLVAL.MLREEL ,IPLERR.MLREEL
  53. POINTEUR IPLA.MLREEL,IPLB.MLREEL,iplc.mlreel
  54. POINTEUR IPERM.MLENTI
  55. LOGICAL CONV,CONV2
  56. REAL*8 XPREC21
  57.  
  58. mlree1=ipla
  59. mlree2=iplb
  60. ildim=numv1
  61. JG=ILDIM
  62. SEGINI,IPLERR,IPERM
  63. N=ILDIM
  64. SEGINI,iprig1,ipmas1
  65. do ibr=1,ildim
  66. iprig1.a(ibr,ibr)=iplc.prog(ibr)
  67. enddo
  68. do ibr=1,ildim-1
  69. ipmas1.a(ibr,ibr)=ipla.prog(ibr)
  70. ipmas1.a(ibr,ibr+1)=iplb.prog(ibr)
  71. ipmas1.a(ibr+1,ibr)=iplb.prog(ibr)
  72. enddo
  73. ipmas1.a(ildim,ildim)=ipla.prog(ildim)
  74. *
  75. c write(6,*) 'simu21: appel sespa3 pour resoudre [T-(1/l)D]Phi=0'
  76. c write(6,*) ' D=',(iplc.prog(iuy),iuy=1,ildim)
  77. c write(6,*) ' alpha=',(ipla.prog(iuy),iuy=1,ildim)
  78. c write(6,*) ' beta=',(iplb.prog(iuy),iuy=1,ildim)
  79.  
  80. *=======================================================*
  81. * Resolution du systeme reduit tridiagonal
  82. *=======================================================*
  83. CALL SESPA3 ( IPMAS1, IPRIG1, IPHI, IPLVAL)
  84. IF(IERR .NE. 0) RETURN
  85.  
  86.  
  87. *=======================================================*
  88. * Test d'arret basé sur un estimateur du residu qui est egalement lié
  89. * a l'erreur relative des valeurs propres
  90. *=======================================================*
  91. *
  92. * REMARQUE IMPORTANTE :
  93. * bp 25/10/2010 : on ne parvient pas a demontrer (ni analytiquement et
  94. * ni numeriquement) ce lien avec l algo de Lanczo.eso tel quil est ecrit
  95. * (utilisation de la K-norme) alors qu il est evident avec la M-norme
  96. * cclusion : il faudrait passer a la M-norme un jour...
  97. * pour l instant, on se contente d etre + severe sur la precision demandee
  98. XPREC0 = XPREC21
  99. XPREC21 = XPREC21**1.5
  100.  
  101. XBETA=IPLB.PROG(ILDIM)
  102. DO 100 IB100 = 1 , ILDIM
  103. XLAMB=IPLVAL.PROG(IB100)
  104. XPHIN= IPHI.A( ILDIM, IB100 )
  105. * erreur absolue
  106. * XERR = ABS( ( XBETA * XPHIN ) )
  107. * erreur relative (sur XLAMB=1/lambda)
  108. XERR = ABS( ( XBETA * XPHIN )/ XLAMB )
  109. IPLERR.PROG( IB100 ) = XERR
  110. 100 CONTINUE
  111.  
  112. * lambda = 1 / (1/ lambda)
  113. DO 200 IB200=1,ILDIM
  114. iplval.prog(ib200)=1.D0/IPLVAL.PROG(IB200)
  115. 200 CONTINUE
  116. c if (iimpi.ge.6) then
  117. c write(6,*) 'simu21: iplval=',(iplval.prog(iou),iou=1,ILDIM)
  118. c write(6,*) 'simu21: IPLERR=',(IPLERR.prog(iou),iou=1,ILDIM)
  119. c endif
  120.  
  121.  
  122. *=========================================================*
  123. * On regarde si les NMODs premieres valeurs ont convergé
  124. *=========================================================*
  125.  
  126. jg=ildim
  127. segini,mlree1,mlree2
  128.  
  129. c----- Stratégie 1 : On trie par module croissant ---------
  130. call trifla(iplval.prog(1),mlree1.prog(1),iplerr.prog(1),
  131. $ mlree2.prog(1),jg,IPERM.lect(1))
  132.  
  133.  
  134. c----- Stratégie 2 : On trie par erreur croissante ---------
  135. c Possible car tri des vp fait par simul7 + tard
  136. c On obtient + de modes par cycle au prix d un + grand nombre de trous
  137. c => implique une modification de strate(gie) pour aller moins loin
  138. c et pour controler le spectre qui est tres diffus...
  139. c call trifla(iplerr.prog(1),mlree2.prog(1),
  140. c $ iplval.prog(1),mlree1.prog(1),jg,IPERM.lect(1))
  141. cbp 10/2010: finalement abandonné car difficile a mettre en oeuvre...
  142.  
  143. segsup,mlree1,mlree2
  144. IF (IERR .NE. 0) RETURN
  145. if (iimpi.ge.6) then
  146. write(6,*) 'simu21: iplval=',(iplval.prog(iou),iou=1,ILDIM)
  147. write(6,*) 'simu21: IPLERR=',(IPLERR.prog(iou),iou=1,ILDIM)
  148. endif
  149.  
  150. c convergence sur l'erreur relative
  151. CONV = .TRUE.
  152. nbonve0=nbonve
  153. nbonve=NMOD
  154.  
  155.  
  156. c------boucle sur les valeur propres
  157. DO 300 IB300 = 1,ildim
  158.  
  159. c -erreur relative < tolerance?
  160. EPS1 = IPLERR.PROG( IB300 )
  161. EPS=ABS(EPS1)
  162. IF ( EPS .GT. XPREC21 ) THEN
  163. if (IB300.eq.1) then
  164. nbonve=ib300-1
  165. if(ib300.le.NMOD+1) CONV=.FALSE.
  166. GOTO 310
  167. endif
  168. c pour eviter de redemarrer a cause d un mode multiple,
  169. c puisque cet estimateur d erreur est assez mechant,
  170. c on va tolerer une erreur relative plus grande
  171. EVP1 = (iplval.prog(IB300)) / (iplval.prog(IB300-1))
  172. EVP1 = ABS (EVP1 - 1.D0)
  173. if (EVP1.LT.XPREC21.AND.EPS.LT.(2.D2*XPREC21)) then
  174. if (iimpi.ge.6) then
  175. write(6,108) IB300,EPS
  176. 108 FORMAT(2X,'on sauve le candidat',I5,
  177. $ ' car mode multiple et erreur acceptable=',E10.5)
  178. endif
  179. elseif(EVP1.LT.(2.D2*XPREC21).AND.IB300.gt.2) then
  180. if (iimpi.ge.6) then
  181. write(6,109) IB300,EPS
  182. 109 FORMAT(2X,'on annule le candidat',I5,' car ',
  183. $ 'probable mode multiple et erreur inacceptable=',E10.5)
  184. endif
  185. nbonve=ib300-2
  186. if(ib300.le.NMOD+1) CONV=.FALSE.
  187. GOTO 310
  188. else
  189. * a priori ce n est pas un mode multiple et erreur inacceptable
  190. nbonve=ib300-1
  191. if(ib300.le.NMOD+1) CONV=.FALSE.
  192. GOTO 310
  193. endif
  194. ENDIF
  195. 300 CONTINUE
  196. c------fin de boucle sur les valeur propres
  197. 310 CONTINUE
  198.  
  199. * on remet la valeur initiale
  200. XPREC21 = XPREC0
  201. *=========================================================*
  202. * destruction des objets de travail
  203. *=========================================================*
  204. mlreel=IPLERR
  205. segsup,mlreel
  206. mlree1=iplval
  207. mlree2=IPERM
  208. SEGDES,mlree1,mlree2
  209. SEGSUP,IPRIG1,IPMAS1
  210.  
  211. RETURN
  212. END
  213.  
  214.  
  215.  

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