Télécharger itinvc.eso

Retour à la liste

Numérotation des lignes :

  1. C ITINVC SOURCE BP208322 15/06/22 21:19:35 8543
  2. SUBROUTINE ITINVC (IPA,IPB,IPRX,IPIX,PROPRE,CONVRG,ITERMX,IPMX)
  3. IMPLICIT INTEGER(I-N)
  4. IMPLICIT REAL*8 (A-H,O-Z)
  5. ************************************************************************
  6. *
  7. * I T I N V C
  8. * -----------
  9. *
  10. * FONCTION:
  11. * ---------
  12. *
  13. * RESOUDRE, PAR ITERATIONS INVERSES, UN SYSTEME D'EQUATIONS:
  14. * |A|.(X) = V.|B|.(X)
  15. * |A| ET |B| ETANT 2 'RIGIDITE',
  16. * (X) UN 'CHPOINT' A DETERMINER ET
  17. * V UN 'FLOTTANT' EGALEMENT A DETERMINER.
  18. *
  19. * ("ITINVC" VAUT POUR IT-ERATIONS INV-ERSES C-OMPLEXES)
  20. *
  21. * MODE D'APPEL:
  22. * -------------
  23. *
  24. * CALL ITINVC (IPA,IPB,IPRX,IPIX,PROPRE,CONVRG,ITERMX,IPMX)
  25. *
  26. * PARAMETRES: (E)=ENTREE (S)=SORTIE
  27. * -----------
  28. *
  29. * IPA ENTIER (E) POINTEUR DE L'OBJET 'RIGIDITE' |A|.
  30. * IPB ENTIER (E) POINTEUR DE L'OBJET 'RIGIDITE' |B|.
  31. * IPRX ENTIER (E) POINTEUR DE L'OBJET 'CHPOINT' DE DEPART.
  32. * (S) PARTIE REELLE DEL'OBJET 'CHPOINT' SOLUTION.
  33. * IPIX ENTIER (S) PARTIE IMAGINAIRE DE L'OBJET SOLUTION
  34. * PROPRE REEL DP (S) TABLEAU CONTENANT DES CARACTERISTIQUES DU
  35. * MODE PROPRE CALCULE. ACTUELLEMENT,
  36. * PROPRE(1) = PARTIE REELLE DU MODE PRORPE,
  37. * PROPRE(2) = (X)T.|B|.(X) , (X) 'CHPOINT'
  38. * SOLUTION,
  39. * PROPRE(3,4,5) DEPL.GEN. REELS SELON X,Y,Z
  40. * PROPRE(6)= PARTIE IMAGINAIRE DU MODE PROPRE
  41. * PROPRE(7)=PARTIE IM. DE XT.|B|.X
  42. * PROPRE(8,9,10) PARTIE IM DES DEP. GEN.
  43. * CONVRG LOGIQUE (S) INDIQUE PAR .TRUE. OU .FALSE. SI LA
  44. * CONVERGENCE A EU LIEU OU NON.
  45. * ITERMX ENTIER (E) NOMBRE MAXIMUM D'ITERATIONS PERMIS.
  46. *
  47. * LEXIQUE: (ORDRE ALPHABETIQUE)
  48. * --------
  49. *
  50. * DIFREL REEL SP VOIR LE S.P. "ITINV1".
  51. * IACCEL ENTIER NOMBRE D'ITERATIONS CONSECUTIVES EFFECTUEES
  52. * SANS ACCELERATION DE CONVERGENCE.
  53. * IPX0 ENTIER VOIR LE S.P. "ITINV1".
  54. * IPX1 ENTIER VOIR LE S.P. "ITINV1".
  55. * IPX2 ENTIER VOIR LE S.P. "ITINV1".
  56. * NBITER ENTIER NOMBRE D'ITERATIONS EFFECTUEES.
  57. * NUMXBX ENTIER NUMERO DE LA DERNIERE ITERATION OU L'ON A
  58. * CALCULE "XT.B.X" POUR 2 'CHPOINT' ITERES
  59. * CONSECUTIFS.
  60. * VALPP REEL DP VALEUR PROPRE REELLE ASSOCIEE AU 'CHPOINT' SOLUTION.
  61. *
  62. * MODE DE FONCTIONNEMENT:
  63. * -----------------------
  64. *
  65. * METHODE DES ITERATIONS INVERSES:
  66. *
  67. * LA SUITE "(X)I" TELLE QUE:
  68. * |A| . (X)I+1 = |B| . (X)I
  69. * TEND VERS LA (OU UNE DES) SOLUTION(S) DE:
  70. * |A| . (X) = V . |B| . (X)
  71. * CORRESPONDANT AU PLUS PETIT V SOLUTION (EN VALEUR ABSOLUE) SOUS
  72. * RESERVE QUE LE (X)1 DE DEPART N'EST PAS B-ORTHOGONAL AU (X)
  73. * SOLUTION.
  74. *
  75. * SOUS-PROGRAMMES APPELES:
  76. * ------------------------
  77. *
  78. * DTCHPO, ITINV1, XTMX, YTX1, RAYLEI, (?), VRFMOD,DEPGE2 ,DTCHPM
  79. *
  80. * AUTEUR, DATE DE MODIFICATION:
  81. * -------------------------
  82. *
  83. * C. LE BIDEAU JUILLET 2001
  84. * Benoit Prabel mars 2009
  85. *
  86. * LANGAGE:
  87. * --------
  88. *
  89. * FORTRAN77 + ESOPE
  90. *
  91. ************************************************************************
  92. *
  93.  
  94. -INC PPARAM
  95. -INC CCOPTIO
  96. -INC SMLMOTS
  97. -INC SMCHPOI
  98. -INC CCHAMP
  99. *
  100. REAL*8 PROPRE(*), XRVP, XIVP
  101. REAL*8 TRNW1(100),TREPS(100)
  102. *
  103. COMMON/CITINV/ NBITER,IACCEL,NUMAC,IPX2,IPX0,IPX1,IPBX1,
  104. C IBBX1,IBBX2,ITPRO,DIFREL
  105. *
  106. LOGICAL CONVRG
  107. *
  108. PARAMETER (INFINI = 9999)
  109. c precision sur la norme 2
  110. * (rem: pour matrices sym, la precision est de 1.D-5 sur la norme infinie => on est donc bcp + exigeant)
  111. PARAMETER (XPREC1 = 1.D-8)
  112. PARAMETER (XPREC2 = 1.D-8)
  113. *
  114. IF (IIMPI .EQ. 747) THEN
  115. CALL GIBTEM(XKT)
  116. INTERR(1)=XKT
  117. CALL ERREUR(-259)
  118. ENDIF
  119. *
  120. * PREPARATION DES ITERATIONS:
  121. IBBX1=0
  122. IBBX2 = IPMX
  123. IPX2 = IPRX
  124. NUMAC = 100
  125. * bp: on n'accelere pas la convergence pour l instant
  126. * car avec des matrices non-sym, la direction du vecteur tourne
  127. * => il faudra adapter l'acceleration au cas complexe
  128. * car parfois cvge tres lente...
  129. NBITER = 0
  130. IACCEL = 0
  131. NUMXBX = -10
  132. X1BX1 = 1.D10
  133. IPLMOX=0
  134. IPLMOY=0
  135. RNW1=1.
  136. REPS=1.
  137. DIFREL = 1.E10
  138. CONVRG = .false.
  139. C
  140. C PREPARATION DES TABLEAUX DONNANT LA CORRESPONDANCE DES NOMS
  141. C D INCONNUE DANS X ET MX STOCKE DANS UN LIST MOT
  142. C
  143. CALL CORRSP(ipa,IPRX,IPMX,IPLMOX,IPLMOY)
  144. C
  145. * -- DEBUT DES ITERATIONS INVERSES -------------------
  146. *
  147. DO 5000 I = 1,ITERMX
  148.  
  149. IF (IBBX1.NE.0) CALL DTCHPO(IBBX1)
  150. IBBX1 = IBBX2
  151.  
  152. IF (I.GE.2) CALL COPIE2 (IPX1,IPREC)
  153. *** réalisation d'une iteration inverse (avec une eventuelle acceleration)
  154. CALL ITINV1 (IPA,IPB)
  155. IF (IERR .NE. 0) RETURN
  156.  
  157. *** on ne teste la convergence qu'apres 3 iterations
  158. IF (I .GT. 3) THEN
  159. *
  160. * ON RECUPERE DEUX VECTEURS IPREC ET IPX1
  161. * 1ER VECTEUR DE BASE
  162. CALL XTY1 (IPREC,IPREC,IPLMOX,IPLMOX,RNNEC)
  163. RNEC = SQRT( ABS(RNNEC))
  164. CALL MUCHPO (IPREC, RNEC, IV, -1)
  165. *
  166. * 2EME VECTEUR DE BASE
  167. CALL XTY1 (IPX1, IV, IPLMOX, IPLMOX, RNX1)
  168. CALL ADCHPO (IPX1, IV, IW1, 1.D0, -RNX1)
  169. CALL XTY1 (IW1, IW1, IPLMOX, IPLMOX, RNNW1)
  170. RNW1 = SQRT( ABS(RNNW1))
  171. TRNW1(I)=RNW1
  172. *
  173. * SI 2E VECTEUR DE BASE NUL : 1 VP REELLE simple
  174. * IF ((I .GE. 3) .AND. (RNW1 .LT. XPREC1)) THEN
  175. IF ((I .GE. 3) .AND. (RNW1 .LT. (XPREC1*RNEC))) THEN
  176. CONVRG = .TRUE.
  177. GOTO 3001
  178. END IF
  179. CALL MUCHPO (IW1, RNW1, IW, -1)
  180.  
  181. * CRITERE DE CONVERGENCE vers 1 paire de vp complexes conjuguées
  182. CALL XTY1 (IPX2, IV, IPLMOX, IPLMOX, P1)
  183. CALL XTY1 (IPX2, IW, IPLMOX, IPLMOX, P2)
  184. IF (IERR .NE. 0 ) RETURN
  185. CALL ADCHPO (IV, IW, IPX4, P1, P2)
  186. CALL ADCHPO (IPX2, IPX4, IEPS, 1.D0, -1.D0)
  187. CALL XTY1 (IEPS, IEPS, IPLMOX, IPLMOX, RNEPS)
  188. IF (IERR .NE. 0 ) RETURN
  189. REPS = SQRT(ABS(RNEPS))
  190. TREPS(I)=REPS
  191. *
  192. * SI IPX2 appartient au s.e. engendré par (IV,IW) : RETOUR VALEURS PROPRES
  193. * IF (RNEPS .LT. 1.E-6) THEN
  194. * IF (REPS .LT. (XPREC2)) THEN
  195. IF (REPS .LT. (XPREC1*RNEC)) THEN
  196. CALL RAYLE2 (IPA,IPB,IV,IW,XRVP,XIVP,IPLMOX,IPLMOY,DELTA,
  197. & IPRX, IPIX, ICVG)
  198. * IF ( ABS(DELTA) .GE. 1.D-14) THEN
  199. IF ( ABS(DELTA) .GE. 1.D-10) THEN
  200. IF (DELTA .LT. 0.D0) THEN
  201. * cas d une paire de valeur propre Complexe
  202. CONVRG = .TRUE.
  203. GOTO 4000
  204. ELSE
  205. * cas d une valeur propre Réelle simple
  206. ** chgt de stratégie concernant l'acceleration de cvgce
  207. * if(NUMAC .ne. 5) then
  208. * NUMAC = 5
  209. * IACCEL = NUMAC - 1
  210. * endif
  211. * GOTO 4900
  212. * cas d une valeur propre Réelle simple pas encore assez bien calculée!
  213. if(ICVG .eq. 0) goto 4900
  214. CONVRG = .TRUE.
  215. * cas de 2 valeurs propre Réelle Double
  216. if(ICVG .eq. 2) GOTO 3002
  217. * cas de 2 valeurs propre Réelle simples de meme module!
  218. if(ICVG .eq. 3) GOTO 3003
  219. END IF
  220. ELSE
  221. * cas d une valeur propre Réelle double
  222. CONVRG = .TRUE.
  223. GOTO 3002
  224. END IF
  225. END IF
  226.  
  227.  
  228. 4900 CONTINUE
  229. * si l'acceleration de cvgce n'ameliore pas le taux de cvgce, alors on arrete
  230. * IF((NUMAC .ne. 100) .and. (IACCEL .EQ. 0)) then
  231. * XTAU2 = (LOG(TREPS(I))) - (LOG(TREPS(I-1)))
  232. * XTAU1 = (LOG(TREPS(I-1))) - (LOG(TREPS(I-2)))
  233. * if(XTAU2 .ge. XTAU1) then
  234. * NUMAC = 100
  235. * endif
  236. * ENDIF
  237. * segdes,MPOVAL,MSOUPO,MCHPOI
  238.  
  239.  
  240. * si on a atteint le nbre maxi d iterations on retourne ce que l 'on a
  241. IF (NBITER .GE. ITERMX) THEN
  242. CONVRG = .FALSE.
  243. * GOTO 302
  244. CALL RAYLE2 (IPA,IPB,IV,IW,XRVP,XIVP,IPLMOX,IPLMOY,DELTA,
  245. & IPRX, IPIX, ICVG)
  246. IF ( ABS(DELTA) .GE. 1.D-10) THEN
  247. IF (DELTA .LT. 0.D0) THEN
  248. GOTO 4000
  249. ELSE
  250. if(XIVP .ne. 0.) GOTO 3002
  251. GOTO 3001
  252. END IF
  253. ELSE
  254. GOTO 3002
  255. END IF
  256.  
  257. END IF
  258.  
  259. NBITER = I
  260.  
  261. *** cas ou I < 3 iterations
  262. ELSE
  263. TRNW1(I)=RNW1
  264. TREPS(I)=REPS
  265.  
  266. END IF
  267. *** fin des tests de convergence réalisés apres 3 iterations
  268. *
  269. *
  270. 5000 CONTINUE
  271. * write(*,*)'-- FIN DE LA BOUCLE DES ITERATIONS INVERSES --------'
  272. CALL DTCHPO(IPX1)
  273. GOTO 4000
  274.  
  275.  
  276. 3001 CONTINUE
  277. * cas d une valeur propre Réelle (simple?) '
  278. * ----------------------------
  279. TRNW1(I)=RNW1
  280. TREPS(I)=1.
  281. * -on n'est pas passé par Rayleigh : il faut faire le travail ici
  282. * normalisation du vecteur propre IPX2
  283. CALL MOTS1(IPLMOT,MOTCLE)
  284. call NORMA1(IPX2,IPLMOT,MOTCLE,IPX2)
  285. * vecteur propre
  286. IPRX = IPX2
  287. IPIX = 0
  288. * produits
  289. CALL MUCPRI(IPX2, IPA, IPAX2)
  290. CALL XTY1(IPX2, IPAX2, IPLMOX, IPLMOY, X2AX2)
  291. IF (IERR .NE. 0 ) RETURN
  292. CALL MUCPRI(IPX2, IPB, IPBX2)
  293. CALL XTY1(IPX2, IPBX2, IPLMOX, IPLMOY, X2BX2)
  294. IF (IERR .NE. 0 ) RETURN
  295. * valeur propre (simple)
  296. XRVP = X2AX2 / X2BX2
  297. XIVP = 0.D0
  298. PROPRE(1) = XRVP
  299. PROPRE(6) = XIVP
  300. * MASSES GENERALISEES
  301. PROPRE(2) = X2BX2
  302. PROPRE(7) = 0.D0
  303. C INTRODUCTION DES COEF. PI OU 2PI EVENTUELS + calcul DEPGEN
  304. CALL MASGEN(IPX2,PROPRE)
  305. CALL DEPGEN( IPB, IPX2,PROPRE,IBBX2,IPLMOX,IPLMOY)
  306. * write(6,*) 'itinvc trouve lambda=',XRVP
  307. GOTO 302
  308.  
  309. 3002 CONTINUE
  310. * cas d une valeur propre Réelle double
  311. * (calculée par Rayleigh) ----
  312. * valeur propre (double)
  313. PROPRE(1) = XRVP
  314. PROPRE(6) = XIVP
  315. * masse généralisée (2: 1er vecteur reel , 7: 2eme vecteur reel)
  316. CALL MUCPRI(IPRX, IPB, IPBRX)
  317. CALL XTY1(IPRX, IPBRX, IPLMOX, IPLMOY, XRBRX)
  318. CALL MUCPRI(IPIX, IPB, IPBIX)
  319. CALL XTY1(IPIX, IPBIX, IPLMOX, IPLMOY, XIBIX)
  320. PROPRE(2) = XRBRX
  321. PROPRE(7) = XIBIX
  322. * write(6,*) 'itinvc trouve lambda=',XRVP,' et ',XIVP
  323. GOTO 302
  324. * ...-> reste a prevoir ce cas dans crebas.eso...
  325.  
  326. 3003 CONTINUE
  327. * cas de 2 valeur propre Réelle simple de meme module
  328. * (calculée par Rayleigh) ----
  329. PROPRE(1) = XRVP
  330. PROPRE(6) = 0.
  331. * masse généralisée (2: 1er vecteur reel , 7: 2eme vecteur reel)
  332. CALL MUCPRI(IPRX, IPB, IPBRX)
  333. CALL XTY1(IPRX, IPBRX, IPLMOX, IPLMOY, XRBRX)
  334. PROPRE(2) = XRBRX
  335. PROPRE(7) = 0.
  336. C INTRODUCTION DES COEF. PI OU 2PI EVENTUELS + calcul DEPGEN
  337. CALL MASGEN(IPRX,PROPRE)
  338. * CALL DEPGEN( IPB, IPRX,PROPRE,IBBX2,IPLMOX,IPLMOY)
  339. * write(6,*) 'itinvc trouve lambda=',XRVP,' (et ',(-1.*XRVP),')'
  340. GOTO 302
  341.  
  342. 4000 CONTINUE
  343. * cas d une paire de valeur propre Complexe'
  344. * (calculée par Rayleigh) ----
  345. * valeur propre lamdba
  346. PROPRE(1) = XRVP
  347. PROPRE(6) = XIVP
  348. * masse complexe générlaisée
  349. CALL MUCPRI(IPRX, IPB, IPBRX)
  350. CALL XTY1(IPRX, IPBRX, IPLMOX, IPLMOY, XRBRX)
  351. CALL XTY1(IPIX, IPBRX, IPLMOX, IPLMOY, XIBRX)
  352. CALL MUCPRI(IPIX, IPB, IPBIX)
  353. CALL XTY1(IPRX, IPBIX, IPLMOX, IPLMOY, XRBIX)
  354. CALL XTY1(IPIX, IPBIX, IPLMOX, IPLMOY, XIBIX)
  355. PROPRE(2) = XRBRX - XIBIX
  356. PROPRE(7) = XIBRX + XRBIX
  357. * write(*,*) 'itinvc trouve lambda=',XRVP,' + i',XIVP
  358. GOTO 302
  359.  
  360. *
  361. * IMPRESSIONS
  362. *
  363. 302 CONTINUE
  364. *
  365. IF (IIMPI.EQ.2) WRITE (IOIMP,2000) NBITER
  366. 2000 FORMAT (//,1X,I3,' ITERATIONS INVERSES ONT ETE EFFECTUEES.'///)
  367. IF (IIMPI.EQ.30) WRITE(IOIMP,1000) (PROPRE(I),I=1,10)
  368. 1000 FORMAT(/10X,'SBR ITINVC',/10X,5(E12.5,1X))
  369. *
  370. * write(6,*) 'i Residu1 Residu2'
  371. * do ii=1,I
  372. * write(6,*) ii,(TRNW1(ii)),(TREPS(ii))
  373. * enddo
  374.  
  375. *
  376. * SUPPRESSION DES 'CHPOINT' DE TRAVAIL:
  377. IF ( (IACCEL + 1) .EQ. NUMAC) THEN
  378. CALL DTCHPO (IPX0)
  379. END IF
  380. CALL DTCHPO (IBBX2)
  381. MLMOTS =IPLMOX
  382. MLMOT1 =IPLMOY
  383. SEGSUP MLMOTS,MLMOT1
  384. *
  385. IF (IIMPI .EQ. 30) THEN
  386. CALL GIBTEM(XKT)
  387. INTERR(1)=XKT
  388. CALL ERREUR(-259)
  389. CALL GIBTEM(XKT)
  390. INTERR(1)=XKT
  391. CALL ERREUR(-259)
  392. END IF
  393. *
  394. END
  395.  
  396.  
  397.  
  398.  
  399.  

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