Télécharger la plaquette Cast3M


fr en

Télécharger itinvc.eso

Retour à la liste

Numérotation des lignes : 

itinvc

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

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