Télécharger arpitq.eso

Retour à la liste

Numérotation des lignes :

  1. C ARPITQ SOURCE PV 15/11/25 21:15:04 8707
  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 CCOPTIO
  63. -INC SMLCHPO
  64. -INC TARWORK
  65. -INC TARTRAK
  66.  
  67. SEGMENT idemen(0)
  68.  
  69.  
  70. INTEGER IPRTRA
  71. INTEGER IPMAUP
  72. LOGICAL INVER
  73. COMPLEX*16 SIGMA
  74. REAL*8 EPSI
  75. LOGICAL OUT
  76.  
  77.  
  78. INTEGER IPRIGI,IPMASS,IPAMOR,IPSCAL
  79. INTEGER TEST
  80. INTEGER IPLCHP
  81. INTEGER IPCTRA(9)
  82. INTEGER NOID,NOEN
  83. INTEGER ndim,ncv,lworkl
  84.  
  85.  
  86. SEGINI IDEMEN
  87. IDEMEN(**)=0
  88. NOID=0
  89. NOEN=1
  90.  
  91. OUT=.FALSE.
  92.  
  93.  
  94. MRITRA=IPRTRA
  95. SEGACT MRITRA
  96. IPRIGI=RIGI(1)
  97. IPMASS=RIGI(2)
  98. IPAMOR=RIGI(3)
  99. IPCSIM=RIGI(5)
  100. IPKCSM=RIGI(6)
  101.  
  102.  
  103. MAUP=IPMAUP
  104. SEGACT MAUP*MOD
  105.  
  106. *Récupération de la dimension des tableaux
  107. ndim=resid(/1)
  108. ncv=v(/2)
  109. lworkl=workl(/1)
  110.  
  111. *Le probleme est non symetrique. Appel de la routine ARPACK specifique
  112. CALL DNAUPD (ido,bmat,ndim,which,nev,EPSI,resid,
  113. & ncv,v,ldv,iparam,ipntr,workd,workl,lworkl,info)
  114.  
  115.  
  116. *Reverse communication: On récupère les paramètres de sortie et on
  117. *effectue des actions en fonction de leurs valeurs
  118. TEST=ido
  119.  
  120. IPMAUP=MAUP
  121. SEGDES MAUP
  122.  
  123.  
  124. *On verifie si on doit stopper le processus
  125. CALL ARPERR (IPMAUP,.FALSE.,OUT)
  126. IF (OUT) RETURN
  127.  
  128.  
  129. *Initialisation des chpoints de travail
  130. DO i=1,9
  131. IPCTRA(i)=0
  132. ENDDO
  133.  
  134.  
  135.  
  136. IF (TEST .EQ. -1 .OR. TEST .EQ. 1) THEN
  137.  
  138. * &--------------------------------------------------&
  139. * | Calcul du produit matrice vecteur |
  140. * | Y <---- RE( inv(A-SIGMA*B)*B*X ) |
  141. * | |
  142. * | X : workd(ipntr(1)) |
  143. * | Y : workd(ipntr(2)) |
  144. * &--------------------------------------------------&
  145.  
  146. *Note: Initialement, ARPACK est prevu pour distinguer test=-1 ou on
  147. *fournit inv(A-SIGMA*B)*B*X et test=1 ou on fournit inv(A-SIGMA*B)*X
  148. *(economie du produit par la matrice masse deja effectue).
  149. *Cependant,la "logique" cast3m nous retire la garantie de fournir le
  150. *produit pour tout type de matrice.
  151. *En outre, meme si la reussite de l'operation etait garantie, elle
  152. *serait plus lourde a effectuer, cette derniere impliquant l'inversion
  153. *d'une matrice de masse
  154.  
  155. IF (AIMAG(SIGMA) .EQ. 0.D0) THEN
  156.  
  157. *Shift reel : RE( inv(A-SIGMA*B)*B*X ) = inv(A-REAL(SIGMA)*B)*B*X
  158.  
  159. CALL ARPCH2 (IPRIGI,IPRIGI,IPMAUP,IPLCHP,1,3)
  160. MLCHPO=IPLCHP
  161. SEGACT MLCHPO
  162.  
  163.  
  164. IPCTRA(3)=ICHPOI(1)
  165.  
  166. CALL MUCPRI (IPCTRA(3),IPMASS,IPCTRA(6))
  167.  
  168. IPCTRA(1)=ICHPOI(2)
  169.  
  170. CALL MUCPRI (IPCTRA(1),IPCSIM,IPCTRA(7))
  171.  
  172. CALL ADCHPO (IPCTRA(6),IPCTRA(7),IPCTRA(2),-1.D0,-1.D0)
  173.  
  174. SEGDES MLCHPO
  175.  
  176. IDEMEN(1)=IPCTRA(2)
  177.  
  178. IF (SYME(6) .EQ. 0) THEN
  179. CALL RESOU1 (IPKCSM,IDEMEN,NOID,NOEN,1.D-18,0)
  180. ELSEIF (SYME(6) .EQ. 1) THEN
  181. CALL LDMT (IPKCSM,IDEMEN,NOID,NOEN,1.D-18)
  182. ENDIF
  183. IF (IERR.NE.0) RETURN
  184.  
  185.  
  186. IPCTRA(4)=IDEMEN(1)
  187.  
  188. CALL ADCHPO(IPCTRA(4),IPCTRA(1),IPCTRA(5),REAL(SIGMA),1.D0)
  189.  
  190. MLCHPO=IPLCHP
  191. SEGACT MLCHPO*MOD
  192. ICHPOI(1)=IPCTRA(5)
  193. ICHPOI(2)=IPCTRA(4)
  194. IPLCHP=MLCHPO
  195. SEGDES MLCHPO
  196.  
  197. CALL ARPCH2 (IPRIGI,IPRIGI,IPMAUP,IPLCHP,2,1)
  198.  
  199. ELSE
  200.  
  201. *Shift complexe : RE( inv(A-SIGMA*B)*B*X ) a fournir (non implemente)
  202.  
  203. ENDIF
  204.  
  205.  
  206. ELSEIF (TEST .EQ. 2) THEN
  207.  
  208. * &------------------------------------------&
  209. * | Calcul du produit matrice vecteur |
  210. * | |
  211. * | Si INVER |
  212. * | Y <---- DIAG(IPRIGI,IPRIGI)*X |
  213. * | |
  214. * | Sinon |
  215. * | Y <---- DIAG(IPMASS,IPMASS)*X |
  216. * | |
  217. * | X : workd(ipntr(1)) |
  218. * | Y : workd(ipntr(2)) |
  219. * &------------------------------------------&
  220.  
  221. IF (INVER) THEN
  222. IPSCAL=IPRIGI
  223. ELSE
  224. IPSCAL=IPMASS
  225. ENDIF
  226.  
  227. CALL ARPCH2 (IPRIGI,IPRIGI,IPMAUP,IPLCHP,1,3)
  228.  
  229. MLCHPO=IPLCHP
  230. SEGACT MLCHPO
  231.  
  232. IPCTRA(1)=ICHPOI(1)
  233.  
  234. CALL MUCPRI (IPCTRA(1),IPSCAL,IPCTRA(3))
  235.  
  236. IPCTRA(2)=ICHPOI(2)
  237.  
  238. CALL MUCPRI (IPCTRA(2),IPSCAL,IPCTRA(4))
  239.  
  240. SEGDES MLCHPO
  241.  
  242. MLCHPO=IPLCHP
  243. SEGACT MLCHPO*MOD
  244. ICHPOI(1)=IPCTRA(3)
  245. ICHPOI(2)=IPCTRA(4)
  246.  
  247. IPLCHP=MLCHPO
  248. SEGDES MLCHPO
  249.  
  250. CALL ARPCH2 (IPRIGI,IPRIGI,IPMAUP,IPLCHP,2,2)
  251.  
  252.  
  253. ENDIF
  254.  
  255.  
  256. *Destruction des chpoints de travail
  257. DO i=1,9
  258. IF (IPCTRA(i) .NE. 0) THEN
  259. CALL DTCHPO (IPCTRA(i))
  260. ENDIF
  261. ENDDO
  262.  
  263. SEGDES MRITRA
  264.  
  265.  
  266. END
  267.  
  268.  
  269.  

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