Télécharger arpitl.eso

Retour à la liste

Numérotation des lignes :

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

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