Télécharger proch3.eso

Retour à la liste

Numérotation des lignes :

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

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