Télécharger arpitl.eso

Retour à la liste

Numérotation des lignes :

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

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