Télécharger arsluc.eso

Retour à la liste

Numérotation des lignes :

  1. C ARSLUC SOURCE CB215821 16/04/15 21:15:02 8907
  2. SUBROUTINE ARSLUC (IPRTRA,SIGMA,IPMAUP,QUAD,EPSI,INVER,
  3. & IPLVAR,IPLVER,IPLVAI,IPLVEI)
  4.  
  5. ************************************************************************
  6. *
  7. * A R S L U C
  8. *
  9. *
  10. *
  11. * FONCTION:
  12. * ---------
  13. *
  14. * CONSTRUCTION DE LISTREEL'S SOLUTION DES FREQUENCES (REELLES ET/OU
  15. * COMPLEXES) ET D'UN MLCHPO'S SOLUTION DES MODES CALCULES (REELS
  16. * ET/OU COMPLEXES) A PARTIR DES VARIABLES DE SORTIE ARPACK
  17. * (PROBLEME NON SYMETRIQUE)
  18. *
  19. *
  20. * PARAMETRES: (E)=ENTREE (S)=SORTIE
  21. * -----------
  22. *
  23. * IPRTRA ENTIER (E) POINTEUR DES OPERATEURS DE TRAVAIL
  24. *
  25. * SIGMA COMPLEX DP (E) VALEUR PROPRE DE SHIFT
  26. *
  27. * INVER LOGIQUE (E) .TRUE. -> PRODUIT SCALAIRE X'KX
  28. * .FALSE. -> PRODUIT SCALAIRE X'MX
  29. *
  30. * IPMAUP ENTIER (E) POINTEUR DES VARIABLES ARPACK
  31. *
  32. * QUAD LOGIQUE (E) PROBLEME QUADRATIQUE OU NON
  33. *
  34. * IPLVAR ENTIER (S) POINTEUR DE L'OBJET 'LISTREEL' CONTENANT
  35. * LA SUITE DES FREQUENCES PROPRES RELLES
  36. *
  37. * IPLVER ENTIER (S) POINTEUR DE L'OBJET 'LISTCHPO' CONTENANT
  38. * LA SUITE DES MODES PROPRES REELS
  39. *
  40. * IPLVAI ENTIER (S) POINTEUR DE L'OBJET 'LISTREEL' CONTENANT
  41. * LA SUITE DES FREQUENCES PROPRES IMAGINAIRES
  42. *
  43. * IPLVEI ENTIER (S) POINTEUR DE L'OBJET 'LISTCHPO' CONTENANT
  44. * LA SUITE DES MODES PROPRES IMAGINAIRES
  45. *
  46. * SOUS-PROGRAMMES APPELES:
  47. * ------------------------
  48. *
  49. * MOTS1,MAXIM1,VCH1,ARPSHI,NORMA1,ARPVER,MUCHPO
  50. *
  51. *
  52. * AUTEUR, DATE DE CREATION:
  53. * -------------------------
  54. *
  55. * PASCAL BOUDA 11 JUILLET 2015
  56. *
  57. ************************************************************************
  58.  
  59. IMPLICIT INTEGER(I-N)
  60. IMPLICIT REAL*8 (A-H,O-Z)
  61.  
  62. -INC CCOPTIO
  63. -INC SMRIGID
  64. -INC SMVECTD
  65. -INC SMLCHPO
  66. -INC SMLREEL
  67. -INC TARWORK
  68.  
  69.  
  70.  
  71. INTEGER IPRTRA
  72. COMPLEX*16 SIGMA
  73. LOGICAL INVER
  74. INTEGER IPMEUP
  75. LOGICAL QUAD
  76. INTEGER IPLVAR
  77. INTEGER IPLVER
  78. INTEGER IPLVAI
  79. INTEGER IPLVEI
  80.  
  81.  
  82. INTEGER IPRIGI
  83. INTEGER JG, N1,INC
  84. INTEGER IPCHO
  85. INTEGER IPVEC
  86. INTEGER IPLMOT
  87. INTEGER MOTCLE
  88. INTEGER TYPRO
  89. INTEGER IPMONO,IPMON1,IPMON2
  90. INTEGER IPMOD1,IPMOD2,IPMO22
  91. REAL*8 MAXVAL, MAXVA1, MAXVA2
  92. INTEGER N
  93. COMPLEX*16 VSHIFT,FSHIFT
  94. COMPLEX*16 CONJ
  95. COMPLEX*16 ZERO
  96.  
  97.  
  98. ZERO =CMPLX(0.D0,0.D0)
  99. VSHIFT=CMPLX(0.D0,0.D0)
  100. FSHIFT=CMPLX(0.D0,0.D0)
  101. i=1
  102.  
  103. MRITRA=IPRTRA
  104. SEGACT MRITRA
  105.  
  106. IPRIGI=RIGI(1)
  107.  
  108. SEGDES MRITRA
  109.  
  110.  
  111. MRIGID=IPRIGI
  112. SEGACT MRIGID
  113. IPCHO=ICHOLE
  114. SEGDES MRIGID
  115.  
  116. MAUP=IPMAUP
  117. SEGACT MAUP
  118.  
  119. IF (QUAD) THEN
  120. N=v(/1)/2
  121. INCR=N
  122. ELSE
  123. N=v(/1)
  124. INCR=0
  125. ENDIF
  126.  
  127. JG=dr(/1)
  128. SEGINI MLREEL
  129. IPLVAR=MLREEL
  130. SEGDES MLREEL
  131.  
  132. SEGINI MLREEL
  133. IPLVAI=MLREEL
  134. SEGDES MLREEL
  135.  
  136. N1=dr(/1)
  137. SEGINI MLCHPO
  138. IPLVER=MLCHPO
  139. SEGDES MLCHPO
  140.  
  141. SEGINI MLCHPO
  142. IPLVEI=MLCHPO
  143. SEGDES MLCHPO
  144.  
  145. *Dans le cas de frequences complexes, il se peut que la nev-ieme
  146. *frequence soit complexe; dans ce cas, ce se sont pas nev mais nev+1
  147. *frequences qui sont calculees. Ainsi dans le cas ou la nev-ieme
  148. *frequence est reelle, on ajuste la taille des objets solutions
  149. IF ( ABS(di(nev+1)) .EQ. 0.D0 ) THEN
  150. JG=JG-1
  151. N1=N1-1
  152. MLREEL=IPLVAR
  153. SEGADJ MLREEL
  154. MLREEL=IPLVAI
  155. SEGADJ MLREEL
  156. MLCHPO=IPLVER
  157. SEGADJ MLCHPO
  158. MLCHPO=IPLVEI
  159. SEGADJ MLCHPO
  160. ENDIF
  161.  
  162. DO 100 k=1,nev
  163.  
  164. *Lecture des valeurs propres: Si la valeur de shift est reelle ou
  165. *nulle, ARPACK fournit directement les valeurs propres. Dans le cas
  166. *contraire, il faut les recalculer a partir des quotients de Rayeigh
  167. *
  168. *En notant lambda=a+ib la valeur propre, V=X+iY le vecteur propre
  169. *associe, et B la matrice pour le produit scalaire, on a :
  170. * - a=X*B*X+Y*B*Y
  171. * - b=X*B*Y-Y*B*X
  172. * 18/08/2015 : Le shift complexe n'est pas implemente
  173.  
  174.  
  175. IF (SIGMA .EQ. ZERO .OR. AIMAG(SIGMA) .EQ. 0.D0) THEN
  176. *Cas "simple": pas besoin de calculer les quotients
  177. VSHIFT=CMPLX(dr(i),di(i))
  178. ELSE
  179. WRITE(IOIMP,*) 'Shift complexe non implemente !'
  180. CALL ERREUR(5)
  181. ENDIF
  182.  
  183. IF (IIMPI .GT. 2) THEN
  184. WRITE(IOIMP,*) ' '
  185. WRITE(IOIMP,*) '--- Mode n°',k,' : on lit a l indice',i,' ---'
  186. WRITE(IOIMP,*) 'Valeur propre',VSHIFT
  187. ENDIF
  188.  
  189.  
  190.  
  191. *Conversion de la valeur propre de shift en frequence de shift
  192. IF (IIMPI.GE.3) WRITE(IOIMP,*) '*** APPEL A ARPSHI ***'
  193. CALL ARPSHI (FSHIFT,VSHIFT,QUAD,INVER,2)
  194. *Les valeurs propres sont classes dans les tableaux dr
  195. *(partie reelle) et di (partie imaginaire). Le tableau de sortie des
  196. *vecteurs propres d'ARPACK prend la convention suivante: si le mode est
  197. *reel, il occupe une colonne du tableau. S'il est complexe, il occupe
  198. *deux colonnes (une pour la partie reelle et une autre pour la partie
  199. *imaginaire.
  200.  
  201.  
  202.  
  203.  
  204. IF ( ABS(di(i)) .EQ. 0.D0 ) THEN
  205.  
  206. IF (IIMPI.GE.3) WRITE(IOIMP,*) 'ARSLUC : *** Cas reel ***'
  207.  
  208. *********************
  209. *******Cas reel******
  210. *********************
  211.  
  212. *La valeur propre est reelle, la frequence est reelle ou imaginaire pure
  213. MLREEL=IPLVAR
  214. SEGACT MLREEL*MOD
  215. PROG(i)=REAL(FSHIFT)
  216. IPLVAR=MLREEL
  217. SEGDES MLREEL
  218.  
  219. MLREEL=IPLVAI
  220. SEGACT MLREEL*MOD
  221. PROG(i)=AIMAG(FSHIFT)
  222. IPLVAI=MLREEL
  223. SEGDES MLREEL
  224.  
  225. *Recuperation du vecteur propre reel pur et transfomation en
  226. chpoint
  227. INC=N
  228. SEGINI MVECTD
  229. DO j=1,N
  230. VECTBB(j)=v(j+INCR,i)
  231. *Tracking : affichage du vecteur propre ARPACK
  232. IF(IIMPI .GT. 300) WRITE(IOIMP,*) j,':',VECTBB(j)
  233. ENDDO
  234.  
  235. IPVEC=MVECTD
  236. SEGDES MVECTD
  237.  
  238. CALL VCH1 (IPCHO,IPVEC,IPMOD1,IPRIGI)
  239. SEGSUP MVECTD
  240.  
  241. *La normalisation est impossible si le chpoint est nul
  242. CALL MOTS1 (IPLMOT,MOTCLE)
  243. CALL MAXIM1 (IPMOD1,IPLMOT,MOTCLE,0,MAXVAL)
  244. IF (MAXVAL .LE. XPETIT) THEN
  245. IPMONO=IPMOD1
  246. ELSE
  247. CALL NORMA1 (IPMOD1,IPLMOT,MOTCLE,IPMONO)
  248. ENDIF
  249.  
  250. *Stockage dans le mlchpo
  251. MLCHPO=IPLVER
  252. SEGACT MLCHPO*MOD
  253. ICHPOI(i)=IPMONO
  254. IPLVER=MLCHPO
  255. SEGDES MLCHPO
  256.  
  257. *Creeation d'un chpoint nul pour la partie imaginaire
  258. INC=N
  259. SEGINI MVECTD
  260. DO j=1,N
  261. VECTBB(j)=0.D0
  262. ENDDO
  263. IPVEC=MVECTD
  264. SEGDES MVECTD
  265. CALL VCH1 (IPCHO,IPVEC,IPMOD2,IPRIGI)
  266. SEGSUP MVECTD
  267.  
  268. *pas de normalisation pour un chpoint nul
  269. IPMONO=IPMOD2
  270.  
  271. MLCHPO=IPLVEI
  272. SEGACT MLCHPO*MOD
  273. ICHPOI(i)=IPMONO
  274. IPLVEI=MLCHPO
  275. SEGDES MLCHPO
  276.  
  277. *Calcul de la norme et du residu du mode
  278. TYPRO=iparam(7)
  279. CALL ARPVER (IPRTRA,TYPRO,i,QUAD,.FALSE.,EPSI,INVER,
  280. & IPMOD1,IPMOD2,VSHIFT)
  281.  
  282.  
  283. ELSE
  284.  
  285. IF (IIMPI.GE.3) WRITE(IOIMP,*) 'ARSLUC : *** Cas complexe ***'
  286. *************************
  287. *******Cas complexe******
  288. *************************
  289.  
  290. *Les valeurs propres sont complexes conjuguees. On travaille par paires
  291. *Note: On "conjugue" la partie reelle dans le cas d'un probleme
  292. *quadratique aux valeurs propres
  293. IF (.NOT. QUAD) THEN
  294.  
  295. MLREEL=IPLVAR
  296. SEGACT MLREEL*MOD
  297. PROG(i)=REAL(FSHIFT)
  298. PROG(i+1)=REAL(FSHIFT)
  299. IPLVAR=MLREEL
  300. SEGDES MLREEL
  301.  
  302. MLREEL=IPLVAI
  303. SEGACT MLREEL*MOD
  304. PROG(i)=AIMAG(FSHIFT)
  305. PROG(i+1)=-AIMAG(FSHIFT)
  306. IPLVAI=MLREEL
  307. SEGDES MLREEL
  308.  
  309. ELSE
  310.  
  311. MLREEL=IPLVAR
  312. SEGACT MLREEL*MOD
  313. PROG(i)=REAL(FSHIFT)
  314. PROG(i+1)=-REAL(FSHIFT)
  315. IPLVAR=MLREEL
  316. SEGDES MLREEL
  317.  
  318. MLREEL=IPLVAI
  319. SEGACT MLREEL*MOD
  320. PROG(i)=AIMAG(FSHIFT)
  321. PROG(i+1)=AIMAG(FSHIFT)
  322. IPLVAI=MLREEL
  323. SEGDES MLREEL
  324.  
  325. ENDIF
  326.  
  327. *Recuperation de la partie reelle du vecteur propre. On travaille par
  328. *paires (partie reele identique pour les deux)
  329. INC=N
  330. SEGINI MVECTD
  331. DO j=1,N
  332. VECTBB(j)=v(j+INCR,i)
  333.  
  334. *Tracking : affichage de la partie superieure du vecteur propre ARPACK
  335. IF (IIMPI .GT. 300) THEN
  336. WRITE(IOIMP,*) J,': RELAMX', v(j,i), 'RECALC',
  337. & REAL(VSHIFT)*v(j+INCR,i)-AIMAG(VSHIFT)*
  338. & v(j+INCR,i+1),'IMLAMX', v(j,i+1),'IMCALC',
  339. & AIMAG(VSHIFT)*v(j+INCR,i)+REAL(VSHIFT)*v(j+INCR,i+1)
  340. ENDIF
  341.  
  342. ENDDO
  343. IPVEC=MVECTD
  344. SEGDES MVECTD
  345.  
  346. IF (IIMPI.GE.3) WRITE(IOIMP,*) '*** APPEL A VCH1 ***'
  347. CALL VCH1 (IPCHO,IPVEC,IPMOD1,IPRIGI)
  348. SEGSUP MVECTD
  349.  
  350. *Recuperation de la partie imaginaire du vecteur propre. On travaille
  351. *par paires (partie imaginaires opposees)
  352. INC=N
  353. SEGINI MVECTD
  354. DO j=1,N
  355. VECTBB(j)=v(j+INCR,i+1)
  356.  
  357. *Tracking : affichage de la partie inferieure du vecteur propre ARPACK
  358. IF (IIMPI .GT. 500) THEN
  359. WRITE(IOIMP,*) J,': REX',v(j+INCR,i),'IMX',v(j+INCR,i+1)
  360. ENDIF
  361.  
  362. ENDDO
  363. IPVEC=MVECTD
  364. SEGDES MVECTD
  365.  
  366. CALL VCH1 (IPCHO,IPVEC,IPMOD2,IPRIGI)
  367. SEGSUP MVECTD
  368.  
  369. CALL MUCHPO (IPMOD2,-1.D0,IPMO22,1)
  370.  
  371.  
  372. *On recupere le max de la paire des chpoints
  373. CALL MOTS1 (IPLMOT,MOTCLE)
  374. CALL MAXIM1 (IPMOD1,IPLMOT,MOTCLE,0,MAXVA1)
  375.  
  376. CALL MOTS1 (IPLMOT,MOTCLE)
  377. CALL MAXIM1 (IPMOD2,IPLMOT,MOTCLE,0,MAXVA2)
  378.  
  379. MAXVAL=MAX(MAXVA1,MAXVA2)
  380.  
  381. *Division de la partie reelle par le max
  382. IF (MAXVA1 .EQ. 0.D0) THEN
  383. IPMONO=IPMOD1
  384. ELSE
  385. CALL MUCHPO (IPMOD1,MAXVAL,IPMONO,-1)
  386. ENDIF
  387.  
  388. *Stockage dans le mlchpo
  389. MLCHPO=IPLVER
  390. SEGACT MLCHPO*MOD
  391. ICHPOI(i)=IPMONO
  392. ICHPOI(i+1)=IPMONO
  393. IPLVER=MLCHPO
  394. SEGDES MLCHPO
  395.  
  396. *Division de la partie imaginaire par le max
  397. IF (MAXVA2 .EQ. 0.D0) THEN
  398. IPMON1=IPMONO
  399. IPMON2=IPMONO
  400. ELSE
  401. CALL MUCHPO (IPMOD2,MAXVAL,IPMON1,-1)
  402. *Construction de la partie imaginaire du conjugue
  403. CALL MUCHPO (IPMOD2,-MAXVAL,IPMON2,-1)
  404. ENDIF
  405.  
  406. *Stockage dans de mlchpo
  407. MLCHPO=IPLVEI
  408. SEGACT MLCHPO*MOD
  409. ICHPOI(i)=IPMON1
  410. IPLVEI=MLCHPO
  411. SEGDES MLCHPO
  412.  
  413. *Stockage dans de mlchpo
  414. MLCHPO=IPLVEI
  415. SEGACT MLCHPO*MOD
  416. ICHPOI(i+1)=IPMON2
  417. IPLVEI=MLCHPO
  418. SEGDES MLCHPO
  419.  
  420. *Calcul de la norme et du residu des modes
  421. TYPRO=iparam(7)
  422.  
  423. CALL ARPVER (IPRTRA,TYPRO,i,QUAD,.FALSE.,EPSI,INVER,
  424. & IPMOD1,IPMOD2,VSHIFT)
  425.  
  426. CONJ=CMPLX(REAL(VSHIFT),-AIMAG(VSHIFT))
  427. CALL ARPVER (IPRTRA,TYPRO,(i+1),QUAD,.FALSE.,EPSI,INVER,
  428. & IPMOD1,IPMO22,CONJ)
  429.  
  430.  
  431. ENDIF
  432.  
  433. *Critere de sortie de boucle
  434. IF ( (i .EQ. JG) .OR. (i .EQ. (JG-1) .AND. di(i) .NE. 0) ) THEN
  435. GOTO 101
  436. ENDIF
  437.  
  438. *L'incrementation differe suivant le type de solution (reelle ou
  439. *complexe)
  440. IF (di(i) .EQ. 0.D0) THEN
  441. i=i+1
  442. ELSE
  443. i=i+2
  444. ENDIF
  445.  
  446. 100 CONTINUE
  447. 101 CONTINUE
  448.  
  449. SEGDES MAUP
  450. SEGDES MRITRA
  451.  
  452. END
  453.  
  454.  
  455.  

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