Télécharger arpitq.eso

Retour à la liste

Numérotation des lignes :

  1. C ARPITQ SOURCE PV 20/05/14 21:15:05 10615
  2. SUBROUTINE ARPITQ (IPRTRA,IPMAUP,SIGMA,INVER,EPSI,OUT)
  3.  
  4.  
  5. ***********************************************************************
  6. *
  7. * A R P I T Q
  8. *
  9. * FONCTION:
  10. * ---------
  11. *
  12. * STEP DE LA FACTORISATION D'ARNOLDI POUR UN PROBLEME QUADRATIQUE.
  13. *
  14. * REMARQUE:
  15. * ---------
  16. *
  17. * LE PROBLEME EST DE DIMENSION DOUBLE PAR RAPPORT AU NOMBRE
  18. * D'INCONNUES: CEPENDANT LES FONCTIONS DISPONIBLES POUR LA
  19. * RESOLUTION DE PROBLEMES D'ORDRE 1 SERONT UTILISES CAR LES
  20. * OPERATIONS MATRICE-VECTEUR SE FONT PAR BLOCS
  21. *
  22. * RIGI(1) : RIGIDITE
  23. * RIGI(3) : AMORTISSEMENT
  24. * RIGI(2) : MASSE
  25. *
  26. *
  27. * PARAMETRES: (E)=ENTREE (S)=SORTIE
  28. * -----------
  29. *
  30. *
  31. * IPRTRA ENTIER (E) OPERATEURS DE TRAVAIL
  32. *
  33. * IPMAUP ENTIER (E/S) POINTEUR VARIABLES ARPACK
  34. *
  35. * SIGMA COMPLEX DP (E) VALEUR DU SHIFT (PEUT ETRE NUL)
  36. *
  37. * INVER LOGIQUE (E) .TRUE. -> PRODUIT SCALAIRE X'KX
  38. * .FALSE. -> PRODUIT SCALAIRE X'MX
  39. *
  40. * EPSI REEL DP (E) TOLERANCE
  41. *
  42. * OUT LOGIQUE (S) FLAG DE CONVERGENCE
  43. *
  44. *
  45. * SOUS-PROGRAMMES APPELES:
  46. * ------------------------
  47. *
  48. * DNAUPD,ADCHPO,ARPCH2,MUCPRI,RESOU1,LDMT,DTCHPO
  49. *
  50. * AUTEUR,DATE DE CREATION:
  51. * -------------------------
  52. *
  53. * PASCAL BOUDA 17 JUILLET 2015
  54. *
  55. * LANGAGE:
  56. * --------
  57. *
  58. * FORTRAN 77 & 90
  59. *
  60. ***********************************************************************
  61.  
  62. -INC PPARAM
  63. -INC CCOPTIO
  64. -INC SMLCHPO
  65. -INC TARWORK
  66. c -INC TARTRAK
  67.  
  68. SEGMENT idemen(0)
  69.  
  70.  
  71. INTEGER IPRTRA
  72. INTEGER IPMAUP
  73. LOGICAL INVER
  74. COMPLEX*16 SIGMA
  75. REAL*8 EPSI
  76. LOGICAL OUT
  77.  
  78.  
  79. INTEGER IPRIGI,IPMASS,IPAMOR,IPSCAL
  80. INTEGER TEST
  81. INTEGER IPLCHP
  82. INTEGER IPCTRA(9)
  83. INTEGER NOID,NOEN
  84. INTEGER ndim,ncv,lworkl
  85.  
  86.  
  87. SEGINI IDEMEN
  88. IDEMEN(**)=0
  89. NOID=0
  90. NOEN=1
  91.  
  92. OUT=.FALSE.
  93.  
  94.  
  95. MRITRA=IPRTRA
  96. SEGACT MRITRA
  97. IPRIGI=RIGI(1)
  98. IPMASS=RIGI(2)
  99. IPAMOR=RIGI(3)
  100. IPCSIM=RIGI(5)
  101. IPKCSM=RIGI(6)
  102.  
  103.  
  104. MAUP=IPMAUP
  105. SEGACT MAUP*MOD
  106.  
  107. *Récupération de la dimension des tableaux
  108. ndim=resid(/1)
  109. ncv=v(/2)
  110. lworkl=workl(/1)
  111.  
  112. *Le probleme est non symetrique. Appel de la routine ARPACK specifique
  113. CALL DNAUPD (ido,bmat,ndim,which,nev,EPSI,resid,
  114. & ncv,v,ldv,iparam,ipntr,workd,workl,lworkl,ITRAK,info)
  115.  
  116.  
  117. *Reverse communication: On récupère les paramètres de sortie et on
  118. *effectue des actions en fonction de leurs valeurs
  119. TEST=ido
  120.  
  121. IPMAUP=MAUP
  122. c SEGDES MAUP
  123.  
  124.  
  125. *On verifie si on doit stopper le processus
  126. CALL ARPERR (IPMAUP,.FALSE.,OUT)
  127. IF (OUT) RETURN
  128.  
  129.  
  130. *Initialisation des chpoints de travail
  131. DO i=1,9
  132. IPCTRA(i)=0
  133. ENDDO
  134.  
  135.  
  136.  
  137. IF (TEST .EQ. -1 .OR. TEST .EQ. 1) THEN
  138.  
  139. * &--------------------------------------------------&
  140. * | Calcul du produit matrice vecteur |
  141. * | Y <---- RE( inv(A-SIGMA*B)*B*X ) |
  142. * | |
  143. * | X : workd(ipntr(1)) |
  144. * | Y : workd(ipntr(2)) |
  145. * &--------------------------------------------------&
  146.  
  147. *Note: Initialement, ARPACK est prevu pour distinguer test=-1 ou on
  148. *fournit inv(A-SIGMA*B)*B*X et test=1 ou on fournit inv(A-SIGMA*B)*X
  149. *(economie du produit par la matrice masse deja effectue).
  150. *Cependant,la "logique" cast3m nous retire la garantie de fournir le
  151. *produit pour tout type de matrice.
  152. *En outre, meme si la reussite de l'operation etait garantie, elle
  153. *serait plus lourde a effectuer, cette derniere impliquant l'inversion
  154. *d'une matrice de masse
  155.  
  156. IF (AIMAG(SIGMA) .EQ. 0.D0) THEN
  157.  
  158. *Shift reel : RE( inv(A-SIGMA*B)*B*X ) = inv(A-REAL(SIGMA)*B)*B*X
  159.  
  160. CALL ARPCH2 (IPRIGI,IPRIGI,IPMAUP,IPLCHP,1,3)
  161. MLCHPO=IPLCHP
  162. SEGACT MLCHPO
  163.  
  164.  
  165. IPCTRA(3)=ICHPOI(1)
  166.  
  167. CALL MUCPRI (IPCTRA(3),IPMASS,IPCTRA(6))
  168.  
  169. IPCTRA(1)=ICHPOI(2)
  170.  
  171. CALL MUCPRI (IPCTRA(1),IPCSIM,IPCTRA(7))
  172.  
  173. CALL ADCHPO (IPCTRA(6),IPCTRA(7),IPCTRA(2),-1.D0,-1.D0)
  174.  
  175. SEGDES MLCHPO
  176.  
  177. IDEMEN(1)=IPCTRA(2)
  178.  
  179. IF (SYME(6) .EQ. 0) THEN
  180. CALL RESOU1 (IPKCSM,IDEMEN,NOID,NOEN,1.D-18,0,0)
  181. ELSEIF (SYME(6) .EQ. 1) THEN
  182. CALL LDMT (IPKCSM,IDEMEN,NOID,NOEN,1.D-18,0)
  183. ENDIF
  184. IF (IERR.NE.0) RETURN
  185. C RESOU desactive IDEMEN on le reactive
  186. SEGACT,IDEMEN
  187.  
  188.  
  189. IPCTRA(4)=IDEMEN(1)
  190.  
  191. CALL ADCHPO(IPCTRA(4),IPCTRA(1),IPCTRA(5),REAL(SIGMA),1.D0)
  192.  
  193. MLCHPO=IPLCHP
  194. SEGACT MLCHPO*MOD
  195. ICHPOI(1)=IPCTRA(5)
  196. ICHPOI(2)=IPCTRA(4)
  197. IPLCHP=MLCHPO
  198. SEGDES MLCHPO
  199.  
  200. CALL ARPCH2 (IPRIGI,IPRIGI,IPMAUP,IPLCHP,2,1)
  201.  
  202. ELSE
  203.  
  204. *Shift complexe : RE( inv(A-SIGMA*B)*B*X ) a fournir (non implemente)
  205.  
  206. ENDIF
  207.  
  208.  
  209. ELSEIF (TEST .EQ. 2) THEN
  210.  
  211. * &------------------------------------------&
  212. * | Calcul du produit matrice vecteur |
  213. * | |
  214. * | Si INVER |
  215. * | Y <---- DIAG(IPRIGI,IPRIGI)*X |
  216. * | |
  217. * | Sinon |
  218. * | Y <---- DIAG(IPMASS,IPMASS)*X |
  219. * | |
  220. * | X : workd(ipntr(1)) |
  221. * | Y : workd(ipntr(2)) |
  222. * &------------------------------------------&
  223.  
  224. IF (INVER) THEN
  225. IPSCAL=IPRIGI
  226. ELSE
  227. IPSCAL=IPMASS
  228. ENDIF
  229.  
  230. CALL ARPCH2 (IPRIGI,IPRIGI,IPMAUP,IPLCHP,1,3)
  231.  
  232. MLCHPO=IPLCHP
  233. SEGACT MLCHPO
  234.  
  235. IPCTRA(1)=ICHPOI(1)
  236.  
  237. CALL MUCPRI (IPCTRA(1),IPSCAL,IPCTRA(3))
  238.  
  239. IPCTRA(2)=ICHPOI(2)
  240.  
  241. CALL MUCPRI (IPCTRA(2),IPSCAL,IPCTRA(4))
  242.  
  243. SEGDES MLCHPO
  244.  
  245. MLCHPO=IPLCHP
  246. SEGACT MLCHPO*MOD
  247. ICHPOI(1)=IPCTRA(3)
  248. ICHPOI(2)=IPCTRA(4)
  249.  
  250. IPLCHP=MLCHPO
  251. SEGDES MLCHPO
  252.  
  253. CALL ARPCH2 (IPRIGI,IPRIGI,IPMAUP,IPLCHP,2,2)
  254.  
  255.  
  256. ENDIF
  257.  
  258.  
  259. *Destruction des chpoints de travail
  260. DO i=1,9
  261. IF (IPCTRA(i) .NE. 0) THEN
  262. CALL DTCHPO (IPCTRA(i))
  263. ENDIF
  264. ENDDO
  265.  
  266. SEGDES MRITRA
  267.  
  268.  
  269. END
  270.  
  271.  
  272.  
  273.  
  274.  
  275.  

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