Télécharger arpitl.eso

Retour à la liste

Numérotation des lignes :

  1. C ARPITL SOURCE BP208322 17/10/03 21:15:05 9580
  2. SUBROUTINE ARPITL (IPRTRA,IPMAUP,SIGMA,INVER,SYM,EPSI,OUT)
  3.  
  4.  
  5. ***********************************************************************
  6. *
  7. * A R P I T L
  8. *
  9. * FONCTION:
  10. * ---------
  11. *
  12. * STEP DE LA FACTORISATION D'ARNOLDI POUR UN PROBLEME LINEAIRE.
  13. *
  14.  
  15. * REMARQUE:
  16. * ---------
  17. *
  18. * ON NOTE:
  19. *
  20. * A=IPRIGI
  21. * B=IPMASS
  22. *
  23. * IPRIGI : RIGIDITE
  24. * IPMASS : MASSE OU KSIGMA
  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 SU SHIFT (PEUT ETRE NUL)
  36. *
  37. * INVER LOGIQUE (E) .TRUE. -> PRODUIT SCALAIRE X'KX
  38. * .FALSE. -> PRODUIT SCALAIRE X'MX
  39. *
  40. * SYM LOGIQUE (E) PROBLEME SYMETRIQUE OU NON
  41. *
  42. * EPSI REEL DP (E) TOLERANCE EIGENPAIRS
  43. *
  44. * OUT LOGIQUE (S) FLAG DE CONVERGENCE
  45. *
  46. *
  47. * SOUS-PROGRAMMES APPELES:
  48. * ------------------------
  49. *
  50. * DSAUPD,DNAUPD,ARPCH1,MUCPRI,RESOU1,LDMT,DECALE,DTCHPO,ARPERR
  51. *
  52. * AUTEUR,DATE DE CREATION:
  53. * -------------------------
  54. *
  55. * PASCAL BOUDA 29 JUIN 2015
  56. *
  57. * LANGAGE:
  58. * --------
  59. *
  60. * FORTRAN 77 & 90
  61. *
  62. ***********************************************************************
  63.  
  64. IMPLICIT INTEGER(I-N)
  65. IMPLICIT REAL*8 (A-H,O-Z)
  66.  
  67. -INC CCOPTIO
  68. -INC CCHAMP
  69. -INC SMRIGID
  70. -INC TARWORK
  71. -INC TARTRAK
  72.  
  73. SEGMENT idemen(0)
  74.  
  75. INTEGER IPRTRA
  76. INTEGER IPMAUP
  77. COMPLEX*16 SIGMA
  78. LOGICAL INVER
  79. LOGICAL SYM
  80. REAL*8 EPSI
  81. LOGICAL OUT
  82.  
  83.  
  84. INTEGER IPRIGI,IPMASS,IPKSIM
  85. INTEGER TEST
  86. CHARACTER*1 SCAL
  87. INTEGER OPT
  88. INTEGER IPCTRA(4)
  89. INTEGER NOID,NOEN
  90. INTEGER ndim,ncv,lworkl
  91.  
  92.  
  93. SEGINI IDEMEN
  94. IDEMEN(**)=0
  95. NOID=0
  96. NOEN=1
  97.  
  98. OUT=.FALSE.
  99.  
  100. MAUP=IPMAUP
  101. SEGACT MAUP*MOD
  102.  
  103. MRITRA=IPRTRA
  104. SEGACT MRITRA
  105.  
  106. *Recuperation des operateurs de travail
  107. IPRIGI=RIGI(1)
  108. IPMASS=RIGI(2)
  109. IPKSIM=RIGI(4)
  110.  
  111.  
  112. *Récupération de la dimension des tableaux
  113. ndim=resid(/1)
  114. ncv=v(/2)
  115. lworkl=workl(/1)
  116.  
  117. *Si le probleme est symétrique, on appelle la routine spécifique aux
  118. *problemes symetriques, sinon on appelle celle pour les problemes
  119. *quelconques
  120.  
  121. IF (SYM) THEN
  122. CALL DSAUPD (ido,bmat,ndim,which,nev,EPSI,resid,
  123. & ncv,v,ldv,iparam,ipntr,workd,workl,lworkl,info)
  124.  
  125. ELSE
  126.  
  127. CALL DNAUPD (ido,bmat,ndim,which,nev,EPSI,resid,
  128. & ncv,v,ldv,iparam,ipntr,workd,workl,lworkl,info)
  129.  
  130. ENDIF
  131.  
  132. *Reverse communication: On récupère les paramètres de sortie et on
  133. *effectue des actions en fonction de leurs valeurs
  134. TEST=ido
  135. SCAL=bmat
  136. OPT=iparam(7)
  137.  
  138. IPMAUP=MAUP
  139. SEGDES MAUP
  140.  
  141. *On verifie si on doit stopper le processus
  142. CALL ARPERR (IPMAUP,SYM,OUT)
  143. IF (OUT) RETURN
  144.  
  145.  
  146. *Initialisation des chpoints de travail
  147. DO i=1,4
  148. IPCTRA(i)=0
  149. ENDDO
  150.  
  151.  
  152. *SCAL: type de probleme
  153. *'I' si standard
  154. *'G' si generalise
  155.  
  156. IF (SCAL .EQ. 'I') THEN
  157.  
  158. IF (IIMPI.GE.10) WRITE(IOIMP,*) '* PB AUX V.P. STANDARD *'
  159.  
  160. IF (TEST .EQ. -1 .OR. TEST .EQ. 1) THEN
  161.  
  162. * &---------------------------------------------------&
  163. * | Calcul du produit matrice vecteur |
  164. * | Y <---- inv(inv(B)*A-SIGMA*I)*X |
  165. * | |
  166. * | X : workd(ipntr(1)) |
  167. * | Y : workd(ipntr(2)) |
  168. * &---------------------------------------------------&
  169.  
  170. ************************************************************************
  171. * 28/08/2015 : Dans ce cas, le shift est obligatoirement nul
  172. * decalage spectral avec une matrice identite non implemente
  173. ************************************************************************
  174.  
  175. CALL ARPCH1 (IPRIGI,IPRIGI,IPMAUP,IPCTRA(3),1,3)
  176.  
  177. CALL MUCPRI (IPCTRA(3),IPMASS,IPCTRA(2))
  178.  
  179. CALL ARPCH1 (IPRIGI,IPRIGI,IPMAUP,IPCTRA(2),1,2)
  180.  
  181. *Mise a sero des inconnues en FLX
  182. CALL ARCORC (IPCTRA(2),10)
  183. *Mise a zero des inconnues en PI inutile ?
  184. CALL ARCORC (IPCTRA(2),15)
  185.  
  186. IDEMEN(1)=IPCTRA(2)
  187.  
  188. IF (SYME(1) .EQ. 0) THEN
  189. CALL RESOU1 (IPRIGI,IDEMEN,NOID,NOEN,1.D-18,0)
  190. ELSEIF (SYME(1) .EQ. 1) THEN
  191. CALL LDMT (IPRIGI,IDEMEN,NOID,NOEN,1.D-18)
  192. ENDIF
  193. IF(IERR.NE.0) RETURN
  194.  
  195. IPCTRA(1)=IDEMEN(1)
  196.  
  197. CALL ARPCH1 (IPRIGI,IPRIGI,IPMAUP,IPCTRA(1),2,1)
  198.  
  199.  
  200. ENDIF
  201.  
  202. ELSEIF (SCAL .EQ. 'G') THEN
  203.  
  204. IF (IIMPI.GE.10)
  205. & WRITE(IOIMP,*) '* PB AUX V.P. GENERALISE *',TEST,INVER
  206.  
  207.  
  208. IF (TEST .EQ. -1) THEN
  209.  
  210. * &--------------------------------------------------&
  211. * | Calcul du produit matrice vecteur |
  212. * | |
  213. * | Y <---- inv(A-SIGMA*B)*B*X |
  214. * | |
  215. * | X : workd(ipntr(1)) |
  216. * | Y : workd(ipntr(2)) |
  217. * &--------------------------------------------------&
  218.  
  219. c WRITE(*,*) 'X1 :'
  220. CALL ARPCH1 (IPRIGI,IPRIGI,IPMAUP,IPCTRA(1),1,3)
  221.  
  222. CALL MUCPRI (IPCTRA(1),IPMASS,IPCTRA(3))
  223. *Mise a sero des inconnues en FLX
  224. CALL ARCORC (IPCTRA(3),10)
  225. *Mise a zero des inconnues en PI inutile ?
  226. CALL ARCORC (IPCTRA(3),15)
  227.  
  228. c WRITE(*,*) '{B*X1} :'
  229. CALL ARPCH1 (IPRIGI,IPRIGI,IPMAUP,IPCTRA(3),1,2)
  230.  
  231. c WRITE(*,*) 'avant RESOU :',SYME(4)
  232. IDEMEN(1)=IPCTRA(3)
  233. IF (SYME(4) .EQ. 0) THEN
  234. CALL RESOU1 (IPKSIM,IDEMEN,NOID,NOEN,1.D-18,0)
  235. ELSEIF (SYME(4) .EQ. 1) THEN
  236. CALL LDMT (IPKSIM,IDEMEN,NOID,NOEN,1.D-18)
  237. ENDIF
  238. IF (IERR.NE.0) RETURN
  239.  
  240. c WRITE(*,*) 'Y1=[OP^-1]*{B*X1} :'
  241. IPCTRA(2)=IDEMEN(1)
  242. cbp CALL ARPCH1 (IPKSIM,IPRIGI,IPMAUP,IPCTRA(2),2,1)
  243. CALL ARPCH1 (IPRIGI,IPRIGI,IPMAUP,IPCTRA(2),2,1)
  244.  
  245.  
  246. ELSEIF (TEST .EQ. 1) THEN
  247.  
  248. * &--------------------------------------------------&
  249. * | Calcul du produit matrice vecteur |
  250. * | |
  251. * | si INVER : |
  252. * | Y <---- inv(A-SIGMA*B)*B*X |
  253. * | |
  254. * | X : workd(ipntr(1)) |
  255. * | Y : workd(ipntr(2)) |
  256. * | |
  257. * | sinon : |
  258. * | Y <---- inv(A-SIGMA*B)*X |
  259. * | |
  260. * | X : workd(ipntr(3)) |
  261. * | Y : workd(ipntr(2)) |
  262. * &--------------------------------------------------&
  263.  
  264. IF (INVER) THEN
  265.  
  266. CALL ARPCH1(IPRIGI,IPRIGI,IPMAUP,IPCTRA(1),1,3)
  267.  
  268. CALL MUCPRI (IPCTRA(1),IPMASS,IPCTRA(3))
  269. *Mise a sero des inconnues en FLX
  270. CALL ARCORC (IPCTRA(3),10)
  271. *Mise a zero des inconnues en PI
  272. CALL ARCORC (IPCTRA(3),15)
  273.  
  274. CALL ARPCH1 (IPRIGI,IPRIGI,IPMAUP,IPCTRA(3),1,2)
  275.  
  276. ELSE
  277.  
  278. c WRITE(*,*) 'X2 :'
  279. cbp CALL ARPCH1 (IPKSIM,IPRIGI,IPMAUP,IPCTRA(3),3,4)
  280. CALL ARPCH1 (IPRIGI,IPRIGI,IPMAUP,IPCTRA(3),3,4)
  281.  
  282. *Mise a sero des inconnues en FLX
  283. CALL ARCORC (IPCTRA(3),10)
  284. *Mise a zero des inconnues en PI
  285. CALL ARCORC (IPCTRA(3),15)
  286. c WRITE(*,*) '{X2} : uniquement le chpoint '
  287.  
  288. ENDIF
  289.  
  290. IDEMEN(1)=IPCTRA(3)
  291.  
  292. IF (SYME(4) .EQ. 0) THEN
  293. CALL RESOU1 (IPKSIM,IDEMEN,NOID,NOEN,1.D-18,0)
  294. ELSEIF (SYME(4) .EQ. 1) THEN
  295. CALL LDMT (IPKSIM,IDEMEN,NOID,NOEN,1.D-18)
  296. ENDIF
  297. IF (IERR.NE.0) RETURN
  298.  
  299. c WRITE(*,*) 'Y2=[OP^-1]*{X2} :'
  300. IPCTRA(2)=IDEMEN(1)
  301. cbp CALL ARPCH1 (IPKSIM,IPRIGI,IPMAUP,IPCTRA(2),2,1)
  302. CALL ARPCH1 (IPRIGI,IPRIGI,IPMAUP,IPCTRA(2),2,1)
  303.  
  304.  
  305. ELSEIF (TEST .EQ. 2) THEN
  306.  
  307. * &-------------------------------------&
  308. * | Calcul du produit matrice vecteur |
  309. * | |
  310. * | Si INVER |
  311. * | Y <---- A*X |
  312. * | |
  313. * | Sinon |
  314. * | Y <---- B*X |
  315. * | |
  316. * | X : workd(ipntr(1)) |
  317. * | Y : workd(ipntr(2)) |
  318. * &-------------------------------------&
  319.  
  320. CALL ARPCH1 (IPRIGI,IPRIGI,IPMAUP,IPCTRA(1),1,3)
  321.  
  322. IF (.NOT. INVER) THEN
  323.  
  324. CALL MUCPRI (IPCTRA(1),IPMASS,IPCTRA(2))
  325. c WRITE(*,*) 'Y3=B*X3 :'
  326.  
  327. ELSE
  328.  
  329. CALL MUCPRI (IPCTRA(1),IPRIGI,IPCTRA(2))
  330.  
  331. *Mise a sero des inconnues en FLX
  332. CALL ARCORC (IPCTRA(2),10)
  333. *Mise a zero des inconnues en PI
  334. CALL ARCORC (IPCTRA(2),15)
  335. c WRITE(*,*) 'Y3={B*X3} :'
  336. ENDIF
  337.  
  338. CALL ARPCH1 (IPRIGI,IPRIGI,IPMAUP,IPCTRA(2),2,2)
  339.  
  340.  
  341. ENDIF
  342.  
  343. ENDIF
  344.  
  345. *Destruction des chpoints de travail
  346. DO i=1,4
  347. IF (IPCTRA(i) .NE. 0) THEN
  348. CALL DTCHPO (IPCTRA(i))
  349. ENDIF
  350. ENDDO
  351.  
  352. SEGDES MRITRA
  353.  
  354. END
  355.  
  356.  
  357.  
  358.  
  359.  

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