Télécharger proch3.eso

Retour à la liste

Numérotation des lignes :

  1. C PROCH3 SOURCE BP208322 19/04/29 21:15:13 10213
  2. C SUBROUTINE PROCH3( FREQ,NBMOD,IPRIGI,IPMASS,INF0,IPMODE,LIMAGE,
  3. C INSYM, MTAB3, I )
  4. ************************************************************************
  5. *
  6. * P R O C H 3
  7. * -----------
  8. *
  9. * FONCTION:
  10. * ---------
  11. *
  12. * RECHERCHE DES NBMOD MODES PROPRES LES PLUS PROCHES DE FREQ
  13. * ET CONSTRUIT L'OBJET 'SOLUTION' CORRESPONDANT.
  14. *
  15. * MODE D'APPEL:
  16. * -------------
  17. *
  18. * CALL PROCH3( FREQ,NBMOD,IPRIGI,IPMASS,INF0,IPMODE,LIMAGE, INSYM )
  19. *
  20. * PARAMETRES: (E)=ENTREE (S)=SORTIE
  21. * -----------
  22. *
  23. * FREQ REEL DP (E) FREQUENCE A APPROCHER PAR UNE FREQUENCE
  24. * PROPRE.
  25. * NBMOD ENTIER (E) ORDRE DE MULTIPLICITE DE LA FREQ PROPRE
  26. * IPRIGI ENTIER (E) POINTEUR DE L'OBJET 'RIGIDITE' CONTENANT
  27. * LA MATRICE DE RIGIDITE.
  28. * IPMASS ENTIER (E) POINTEUR DE L'OBJET 'RIGIDITE' CONTENANT
  29. * LA MATRICE MASSE.
  30. * IPMODE ENTIER (S) POINTEUR DE L'OBJET 'SOLUTION' CONTENANT
  31. * LE MODE PROPRE TROUVE.
  32. * LIMAGE BOOLEEN (E) VRAI SI ON SOUHAITE DE FREQ. NEGATIVES.
  33. *
  34. * LEXIQUE: (ORDRE ALPHABETIQUE)
  35. * --------
  36. *
  37. * IPKW2M ENTIER POINTEUR DE LA 'RIGIDITE' "DECALEE" K - W2.M
  38. * IPLVER LISTCHPO LISTE DE CHPO A ITERER. ILS FORMENT UN SOUS-
  39. * ESPACE QUI CONVERGE VERS LE SOUS-ESPACE PROPRE
  40. * DOMINANT DE DIMENSION NBMOD. SORTIE: PARTIE REELLE
  41. * DES CHPOINTS SOLUTIONS
  42. * IPLVEI LISTCHPO LISTE DE CHPO.SORTIE: PARTIE IMAGINAIRE
  43. * DES CHPOINTS SOLUTIONS
  44. * IPVECP ENTIER POINTEUR DU 'CHPOINT' QUI CONTIENT DES NOMBRES
  45. * ALEATOIRES, PUIS UN VECTEUR PROPRE.
  46. * OMEGA2 REEL DP PULSATION PROPRE TROUVEE AU CARRE.
  47. * ( NE SERT A RIEN ! )
  48. * PROPRE REEL DP TABLEAU CONTENANT DES CARACTERISTIQUES DU
  49. * MODE PROPRE CALCULE:
  50. * PROPRE(1) = FREQ. PROPRE ,
  51. * PROPRE(2) = MASSE GEN.
  52. * PROPRE(3)ET(4) ET(5) DEPL.GEN. SELON X,Y,Z
  53. * W2 REEL DP PULSATION AU CARRE A APPROCHER.
  54. *
  55. * MODE DE FONCTIONNEMENT:
  56. * -----------------------
  57. *
  58. * LE CALCUL D'UN VECTEUR PROPRE SE FAIT PAR LA METHODE DE
  59. * L'ITERATION D'UN SOUS-ESPACE, AVEC DECALAGE INITIAL ("SHIFTING")
  60. * DE W2
  61. *
  62. * AUTEUR, DATE DE CREATION:
  63. * -------------------------
  64. *
  65. * PASCAL MANIGOT 16 OCTOBRE 1984
  66. *
  67. * MODIF: Benoit PRABEL MARS 2009
  68. *
  69. * LANGAGE:
  70. * --------
  71. *
  72. * ESOPE
  73. *
  74. ************************************************************************
  75.  
  76. SUBROUTINE PROCH3 ( FREQ,NBMOD,IPRIGI,IPMASS,INF0,IPMODE,LIMAGE
  77. & ,INSYM, MTAB3, I)
  78.  
  79. IMPLICIT INTEGER(I-N)
  80. IMPLICIT REAL*8 (A-H,O-Z)
  81.  
  82. -INC CCREEL
  83. -INC CCOPTIO
  84. -INC SMLCHPO
  85. -INC SMLREEL
  86. -INC SMCHPOI
  87.  
  88. * -- CONSTANTES --
  89.  
  90. PARAMETER (LPROPR = 10)
  91. PARAMETER (DEUXPI = (2.D0*XPI))
  92.  
  93.  
  94. * -- ARGUMENTS --
  95.  
  96. REAL*8 FREQ
  97. INTEGER NBMOD, IPRIGI, IPMASS, INF0, IPMODE
  98. LOGICAL LIMAGE
  99. CHARACTER*4 MOTCLE
  100.  
  101.  
  102. * -- VARIABLES LOCALES --
  103.  
  104. POINTEUR IPLVER.MLCHPO, IPLVE1.MLCHPO, IPLVEI.MLCHPO
  105. POINTEUR IPLVAR.MLREEL, IPLVA1.MLREEL, IPLVAI.MLREEL
  106. INTEGER NBVEC
  107. REAL*8 OMEGA2, PROPRE(LPROPR)
  108.  
  109. ******************************************************************
  110. * INITIALISATION *
  111. ******************************************************************
  112.  
  113.  
  114. * -- CREATION DE (K-W2M) --
  115.  
  116. W2 = (FREQ * DEUXPI) ** 2
  117. IF( LIMAGE ) THEN
  118. W2 = SIGN(W2,FREQ)
  119. ENDIF
  120. CALL DECALE ( IPRIGI, IPMASS, W2, IPKW2M )
  121. IF ( IERR .NE. 0 ) RETURN
  122.  
  123.  
  124. * -- INITIALISATION DES ITERATIONS: CREATION D'UNE LISTE
  125. * DE NBVEC 'CHPOINT' ALEATOIRE --
  126.  
  127.  
  128. * bp (03.2009) : si on souhaite eventuellement fournir 1 CHPOINT
  129. call LIROBJ('LISTCHPO',MLCHP1,0,IRET1)
  130. if(IRET1 .ne. 0) then
  131. segact,MLCHP1
  132. NCHPO1 = MLCHP1.ICHPOI(/1)
  133. IPVECP = 0
  134. * call COPIE2(MLCHP1.ICHPOI(1),IPVECP)
  135. ICHP1 = MLCHP1.ICHPOI(1)
  136. call COPIE2(ICHP1,IPVECP)
  137. * A Faire : vérifier l'accord entre chpoint et matrices...?
  138. else
  139. CALL ALEAT1 (IPKW2M, IPVECP )
  140. endif
  141. IF ( IERR .NE. 0 ) RETURN
  142. N1 = 1
  143. SEGINI ,IPLVER
  144. IPLVER.ICHPOI(1) = IPVECP
  145. ** on cherce le nombre de ddl sans relations
  146. nddl = 0
  147. ndlx = 0
  148. mchpoi = ipvecp
  149. segact mchpoi
  150. nsous = ipchp(/1)
  151. do 81 isous = 1, nsous
  152. msoupo = ipchp(isous)
  153. segact msoupo
  154. mpoval = ipoval
  155. segact mpoval
  156. nn = vpocha(/1)
  157. nc1 = vpocha(/2)
  158. nddl = nddl + nn*nc1
  159. do 95 inc = 1,nc1
  160. if (nocomp(inc).eq. 'LX ') ndlx = ndlx + nn
  161. 95 continue
  162. segdes msoupo, mpoval
  163. 81 continue
  164. segdes mchpoi
  165.  
  166. nddl = nddl - ndlx - (ndlx/2)
  167. nbvec = max(2*nbmod, 8)
  168. if (nbvec .gt. nddl) nbvec = nddl
  169.  
  170. * on continue le tirage au sort (nous aurions pu passer par Lanzos!)
  171. DO 100 IB100 = 1, NBVEC-1
  172. if ( (IRET1 .ne. 0) .and. (IB100 .lt. NCHPO1) ) then
  173. IPVECP = 0
  174. * call COPIE2(MLCHP1.ICHPOI(1 + IB100),IPVECP)
  175. ICHP1 = MLCHP1.ICHPOI(1 + IB100)
  176. call COPIE2(ICHP1,IPVECP)
  177. else
  178. CALL ALEAT1 (IPKW2M, IPVECP )
  179. endif
  180. IF ( IERR .NE. 0 ) RETURN
  181. IPLVER.ICHPOI(**) = IPVECP
  182. 100 CONTINUE
  183. SEGDES ,IPLVER
  184. if(IRET1 .ne. 0) SEGDES,MLCHP1
  185.  
  186. * cette orthonormalisation n'est pas indispensable, dans la plupart
  187. * des cas,
  188. * mais, on l'a gardée pour tester la singularité de la masse
  189. * et mettre à jour nbvec si necessaire. On n'a besoin de ce test
  190. * qu'au debut, car si l'on a eliminé les eventuels vecteurs du noyau
  191. * de la masse ils ne peuvent pas ressurgir lors du processus iteratif
  192. * car ils correspondent à des frequences infinies
  193. * IF ( INSYM .EQ. 0) THEN
  194. CALL GRAAMO ( IPLVER, IPMASS )
  195. IF ( IERR .NE. 0 ) RETURN
  196. segact IPLVER
  197. nbvec = IPLVER.ICHPOI(/1)
  198. segdes iplver
  199. if (nbmod .gt. nbvec) nbmod = nbvec
  200. * END IF
  201.  
  202.  
  203. ******************************************************************
  204. * RESOLUTION *
  205. ******************************************************************
  206.  
  207. * -- RESOLUTION PAR ITERATION DU SOUS-ESPACE IPLVER
  208. * ET NORMALISATION DES MODES PROPRES --
  209.  
  210. *on differentie NBMOD=dim du sous espace initialement recherché et
  211. *NBMOD2 qui peut etre modifié par sespac si dernier vecteur complexe
  212. NBMOD2 = NBMOD
  213.  
  214. IF (INSYM .EQ. 0) THEN
  215.  
  216. CALL SESPA ( IPLVAR, IPLVER, NBMOD, IPKW2M, IPMASS )
  217. IF ( IERR .NE. 0 ) RETURN
  218. segact iplver
  219. nblgn = IPLVER.ICHPOI(/1)
  220. ICOMP = 0
  221. mlchpo = iplver
  222. segact mlchpo*mod
  223. do 101 ibmod = 1,nbmod
  224. ix = ichpoi(ibmod)
  225. CALL MOTS1( IPLMOT, MOTCLE )
  226. call norma1 (ix,iplmot,motcle,ixx)
  227. call dtchpo(ix)
  228. ichpoi(ibmod) = ixx
  229. 101 continue
  230. segdes mlchpo
  231.  
  232. ELSE
  233.  
  234. CALL SESPAC(IPLVAR,IPLVAI,IPLVER,NBMOD2,IPKW2M,IPMASS)
  235. ICOMP = 0
  236. *
  237. END IF
  238.  
  239.  
  240. ******************************************************************
  241. * CREATION DE L'OBJET SOLUTION *
  242. ******************************************************************
  243.  
  244. * -- TRANSLATION DES FREQUENCES PROPRES --
  245.  
  246. SEGACT ,IPLVAR*MOD
  247. * IF (ICOMP .EQ. 0) THEN
  248. *---- Cas d'un appel a sespa (matrices symetriques => mode Reel)
  249. IF (INSYM .EQ. 0) THEN
  250. DO 200 IB200 =1, NBMOD
  251. FREQPR = IPLVAR.PROG(IB200)
  252. CALL W2FREQ ( FREQPR, W2, OMEGA2, FREQP1, LIMAGE )
  253. IF ( IERR .NE. 0 ) RETURN
  254. IPLVAR.PROG(IB200) = FREQP1
  255. 200 CONTINUE
  256.  
  257. *---- Cas d'un appel a sespaC (matrices nonsymetriques => mode Reel ou Complexe)
  258. ELSE
  259. SEGACT IPLVAI*MOD
  260. * DO 400 IB400 = 1, NBMOD2
  261. DO 400 IB400 = 1, NBMOD2
  262. FREQPR = IPLVAR.PROG(IB400)
  263. FREQPI = IPLVAI.PROG(IB400)
  264. CALL W2FRQC ( FREQPR, FREQPI, W2,XRW2,XIW2, FREQP1, FREQP2)
  265. IF ( IERR .NE. 0 ) RETURN
  266. * IPLVAR.PROG(IB400) = FREQP1
  267. * IPLVAI.PROG(IB400) = FREQP2
  268. * on se débarasse des erreurs d'arrondis de W2FRQC
  269. if (FREQPI .eq. 0.) then
  270. if ( (FREQPR + W2) .lt. 0.) then
  271. IPLVAR.PROG(IB400)= 0.D0
  272. IPLVAI.PROG(IB400)= FREQP2
  273. else
  274. IPLVAR.PROG(IB400)= FREQP1
  275. IPLVAI.PROG(IB400)= 0.D0
  276. endif
  277. else
  278. IPLVAR.PROG(IB400)= FREQP1
  279. IPLVAI.PROG(IB400)= FREQP2
  280. endif
  281. * WRITE(6,*) IB400,' FREQUENCE =',FREQP1,'+i',FREQP2,
  282. * & ' soit',(IPLVAR.PROG(IB400)),'+i',(IPLVAI.PROG(IB400))
  283. 400 CONTINUE
  284. END IF
  285.  
  286. * -- ON GARDE LES NBMOD2 PREMIERS ELEMENTS PROPRES --
  287. * -- ET ON DETRUIT LES AUTRES --
  288.  
  289. * -- LES VALEURS PROPRES --
  290.  
  291. * JG = NBMOD
  292. JG = NBMOD2
  293. SEGADJ IPLVAR
  294. SEGDES IPLVAR
  295.  
  296. IF (INSYM .NE. 0) THEN
  297. * JG = NBMOD
  298. JG = NBMOD2
  299. SEGADJ IPLVAI
  300. SEGDES IPLVAI
  301. ENDIF
  302.  
  303. * -- LES MODES PROPRES --
  304. SEGACT IPLVER
  305. * si le dernier mode est complexe, alors il nous faut le (nbmod + 1)^ieme
  306. * qui correspond a la partie imaginaire indissociable de la partie reelle,
  307. * ce qui justifie l'introduction de nbmod2
  308. DO 550 IB550 = (NBMOD2 + 1), NBVEC
  309. IPCHPO = IPLVER.ICHPOI(IB550)
  310. CALL DTCHPO ( IPCHPO )
  311. 550 CONTINUE
  312. N1 = NBMOD2
  313. SEGADJ IPLVER
  314. SEGDES IPLVER
  315.  
  316.  
  317. * -- CLASSEMENT DES MODES PROPRES --
  318.  
  319. IF (INSYM .eq. 0) THEN
  320. c on n'a retenu que le nombre de mode demandes,
  321. * on trie par lambda croissant de maniere a avoir le bon numero
  322. CALL ORDVEC ( IPLVAR, IPLVER,.false. )
  323. ENDIF
  324. * rem: mode complexes deja classés par module de la val propre croissant
  325.  
  326.  
  327. * -- AFFICHAGE --
  328.  
  329. IF ( IIMPI .EQ. 2 ) THEN
  330. WRITE ( IOIMP, 1000 ) FREQ
  331. 1000 FORMAT( /8X, 'SOLUTION POUR LA FREQUENCE:', E12.6,
  332. 1 /8X, '---------------------------', / )
  333. ENDIF
  334.  
  335.  
  336. * -- CONSTRUCTION DE LA SOLUTION --
  337.  
  338. IF (INSYM .EQ. 0) THEN
  339. CALL CRSOLU (W2,IPLVAR,IPLVER,NBMOD,IPKW2M,IPMASS,IPMODE)
  340. ELSE
  341. CALL CRSOL1(W2,IPLVAR,IPLVAI,IPLVER,NBMOD,NBMOD2,
  342. & IPKW2M,IPMASS,MTAB3,I,INF0)
  343. END IF
  344.  
  345.  
  346. * -- SUPPRESSION DES OBJETS DE TRAVAIL --
  347.  
  348. CALL DTRIGI ( IPKW2M )
  349. IF ( IERR .NE. 0 ) RETURN
  350. CALL DTLCHP ( IPLVER )
  351. IF ( IERR .NE. 0 ) RETURN
  352. CALL DTLREE ( IPLVAR )
  353. IF ( IERR .NE. 0 ) RETURN
  354. IF (INSYM .NE. 0) CALL DTLREE ( IPLVAI )
  355.  
  356. RETURN
  357. END
  358.  
  359.  
  360.  
  361.  
  362.  
  363.  
  364.  
  365.  
  366.  
  367.  
  368.  
  369.  
  370.  
  371.  
  372.  
  373.  
  374.  
  375.  

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