Télécharger iplnu1.eso

Retour à la liste

Numérotation des lignes :

iplnu1
  1. C IPLNU1 SOURCE BP208322 21/01/28 21:15:10 10868
  2. SUBROUTINE IPLNU1(INUA)
  3. C
  4. C Fonction
  5. C Cette sous routine interpole une fonction de n variables dont
  6. C certains couples (x,f(x)) ont ete stockes dans un objet de type
  7. C nuage
  8. C
  9. C Variables
  10. C INUA pointeur sur l'objet de type nuage
  11. C NN nombre de composantes du CHPOINT/MCHAML en entree
  12. C MM nombre de composantes du CHPOINT/MCHAML a calculer
  13. C NNMM nombe de composantes du nuage (= NN + MM)
  14. C IADD contient la correspondance entre les numeros des composantes
  15. C du CHPOINT/MCHAML et celles du nuage :
  16. C - La composante connue I (entre 1 et NN) du champ est a la position
  17. C IADD(I) dans le nuage
  18. C - La composante a calculer J (entre 1 et MM) est a la position
  19. C IADD(J) dans le nuage
  20. C MAXV contient les max des composantes du nuage (utilise pour normer
  21. C le calcul des distances)
  22. C
  23. C Appele par interp qui est le point d'entree de l'operateur IPOL
  24. C
  25. C Auteur A De Gayffier
  26. C Date 05/07/94
  27. C
  28. C Evolution 2015 : nouvelle option GRILL, interpolation dans une grille
  29. C Evolution 2016 : nouvelle option PID, interpolation par ponderation
  30. C inverse a la distance
  31. C
  32. C Langage esope+fortran 77
  33. C
  34.  
  35. IMPLICIT INTEGER(I-N)
  36. IMPLICIT REAL*8(A-H,O-Z)
  37.  
  38. LOGICAL FLAG
  39.  
  40. -INC SMNUAGE
  41. -INC SMCHAML
  42. -INC CCREEL
  43.  
  44. -INC PPARAM
  45. -INC CCOPTIO
  46. -INC SMCHPOI
  47. -INC CCASSIS
  48.  
  49. C Introduction d'un COMMON pour la parallelisation
  50. C Il est a repercuter aussi dans les sources suivantes :
  51. C - ipln2i.eso
  52. COMMON/iplnuc/XP1,NBTH1,INUA1,ITR1,IMAX1,IVAL1,N1,NN1,MM1,NNMM1,
  53. & IMPOV1,IMPOV2,
  54. & IMCHA1,IMCHA2,N3EL,N3PTEL,
  55. & ITR2,XELIM1
  56.  
  57. EXTERNAL ipln2i
  58.  
  59.  
  60. SEGMENT MTR1
  61. INTEGER IADD(NNMM)
  62. ENDSEGMENT
  63.  
  64. SEGMENT MTR2
  65. REAL*8 XI(N,NN)
  66. REAL*8 YI(N,MM)
  67. ENDSEGMENT
  68.  
  69.  
  70. SEGMENT MAXI1
  71. REAL*8 MAXV(NNMM)
  72. ENDSEGMENT
  73.  
  74. CHARACTER*5 CLE(4)
  75. DATA CLE /'GAUSS','RATIO','PID','GRILL'/
  76. CHARACTER*4 CLE2(1)
  77. DATA CLE2 /'ELIM'/
  78.  
  79. PARAMETER (TINY=1.D-20)
  80.  
  81.  
  82. INUAG = INUA
  83. IMPOV1 = 0
  84. IMPOV2 = 0
  85. IMCHA1 = 0
  86. IMCHA2 = 0
  87. N3EL = 0
  88. N3PTEL = 0
  89. ITR2 = 0
  90.  
  91. ITHRD = 0
  92. C
  93. C Option pour la fonction de ponderation
  94. CALL LIRMOT(CLE,4,IVAL,0)
  95. IF (IVAL .EQ. 0) IVAL = 1
  96. C
  97. C (FDP 2016) Acquisition d'un flottant/entier pour le mot clef PID
  98. XP=1.D0
  99. XELIM=XPETIT
  100. IF (IVAL.EQ.3) THEN
  101. CALL LIRREE(XP,0,IRETOU)
  102. IF (XP.LT.0.D0) THEN
  103. REAERR(1)=REAL(0.D0)
  104. REAERR(2)=REAL(XP)
  105. CALL ERREUR(191)
  106. RETURN
  107. ENDIF
  108. c + un flottant supplementaire pour la tolerance a l'ELIMination
  109. c (i.e. distance en deca de laquelle les noeuds sont confondus)
  110. CALL LIRMOT(CLE2,1,IVAL2,0)
  111. IF(IVAL2.EQ.1) THEN
  112. CALL LIRREE(XELIM,0,IRETOU)
  113. IF(IRETOU.EQ.0) THEN
  114. c XP est en fait XELIM ! -> on intervertit
  115. XELIM=XP
  116. XP=1.D0
  117. ELSE
  118. c on a lu XP et XELIM (positionnel) -> on teste XELIM
  119. IF (XELIM.LT.0.D0) THEN
  120. REAERR(1)=REAL(0.D0)
  121. REAERR(2)=REAL(XELIM)
  122. CALL ERREUR(191)
  123. RETURN
  124. ENDIF
  125. ENDIF
  126. ENDIF
  127. ENDIF
  128. c WRITE(IOIMP,*) IVAL,'>>> XP,XELIM=',XP,XELIM
  129. C
  130. C (FDP 2015) Interpolation dans un nuage representant une grille de valeurs
  131. C La parallelisation est faite dans IPGRIL !
  132. IF (IVAL.EQ.4) THEN
  133. CALL IPGRIL(INUA)
  134. RETURN
  135. ENDIF
  136.  
  137. C
  138. MNUAGE = INUA
  139.  
  140. NNMM = NUANOM(/2)
  141. SEGINI ,MAXI1
  142.  
  143. C recherche du maximum de chaque parametre
  144.  
  145. DO 20 K=1,NNMM
  146. MAXV(K) = 0
  147. NUAVFL = NUAPOI(K)
  148. NBCOUP = NUAFLO(/1)
  149. MAXV(K)=TINY
  150. DO 10 J=1,NBCOUP
  151. IF(MAXV(K).LT.ABS(NUAFLO(J))) THEN
  152. MAXV(K) = ABS(NUAFLO(J))
  153. ENDIF
  154. 10 CONTINUE
  155. 20 CONTINUE
  156. C
  157. CALL LIROBJ('CHPOINT',IPOI,0,IRETOU)
  158. IF ( IRETOU .EQ. 1 ) THEN
  159. CALL ACTOBJ('CHPOINT',IPOI,1)
  160. GOTO 210
  161. ENDIF
  162.  
  163. CALL LIROBJ('MCHAML',ICHML,0,IRETOU)
  164. IF ( IRETOU .EQ. 1 ) THEN
  165. CALL ACTOBJ('MCHAML',ICHML,1)
  166. ELSE
  167. C pas d'argument correct trouve: erreur
  168. CALL QUETYP(MOTERR(1:8),0,IRETOU)
  169. IF(IRETOU.NE.0) THEN
  170. CALL ERREUR(39)
  171. ELSE
  172. CALL ERREUR(533)
  173. ENDIF
  174. RETURN
  175. ENDIF
  176.  
  177. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  178. C Cas d'un MCHAML C
  179. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  180.  
  181. MCHELM = ICHML
  182. C creation du chapeau du chamelem resultat
  183. SEGINI ,MCHEL1 = MCHELM
  184. C
  185. N1 = ICHAML(/1)
  186. C
  187. C boucle sur les chamelems elementaires
  188. C
  189. DO 200 I=1,N1
  190. MCHAML = ICHAML(I)
  191. C calcul des dimensions nn et pp
  192. NNMM = 0
  193. SEGINI MTR1
  194. NN = NOMCHE(/2)
  195. MM = NUANOM(/2) - NN
  196. C initialisation du nouveau sous chamelem
  197. N2 = MM
  198. SEGINI ,MCHAM1
  199. MCHEL1.ICHAML(I) = MCHAM1
  200. C
  201. C on verifie que les noms des composantes du chamelem sont
  202. C dans le nuage et on rempli iadd
  203. C
  204. DO 110 J=1,NOMCHE(/2)
  205. IF ( TYPCHE(J) .NE. 'REAL*8' ) THEN
  206. MOTERR(1:8) = NOMCHE(J)
  207. MOTERR(9:13) = 'CHAMP'
  208. CALL ERREUR(681)
  209. C
  210. C 'la composante',nomche(j),'du chamelem''n est pas scalaire'*
  211. SEGSUP MTR1,MCHAM1,MCHEL1
  212. RETURN
  213. ENDIF
  214. C
  215. DO 100 K=1,NUANOM(/2)
  216. IF (NUANOM(K) .EQ. NOMCHE(J)) THEN
  217. IADD(**)=K
  218. NNMM = IADD(/1)
  219. ENDIF
  220. 100 CONTINUE
  221. IF ( NNMM .NE. J ) THEN
  222. C une des composantes du champ n'est pas un composante du nuage
  223. CALL ERREUR(682)
  224. SEGSUP MTR1,MCHAM1,MCHEL1
  225. C return
  226. ENDIF
  227. 110 CONTINUE
  228. C
  229. C on recupere les composantes du nuage qui ne sont pas celles
  230. C du chamelem
  231. C
  232. DO 130 J=1,NUANOM(/2)
  233. IF ( NUATYP(J) .NE. 'FLOTTANT' ) THEN
  234. MOTERR(1:8) = NUANOM(J)
  235. MOTERR(9:13) = 'NUAGE'
  236. CALL ERREUR(681)
  237. C 'la composante', nuanom(j), 'du nuage', 'n est pas scalaire'
  238. SEGSUP MTR1,MCHAM1,MCHEL1
  239. RETURN
  240. ENDIF
  241. FLAG = .TRUE.
  242. DO 120 K=1,NOMCHE(/2)
  243. IF (NUANOM(J) .EQ. NOMCHE(K)) THEN
  244. FLAG = .FALSE.
  245. ENDIF
  246. 120 CONTINUE
  247. IF (FLAG) THEN
  248. C la composante du nuage n'est pas dans le chamelem
  249. IADD(**) = J
  250. NNMM = IADD(/1)
  251. ENDIF
  252. 130 CONTINUE
  253. IF ( NNMM .NE. NUANOM(/2)) THEN
  254. CALL ERREUR(682)
  255. C 'incoherence entre les composantes du nuage et du champ'
  256. SEGSUP MTR1,MCHAM1,MCHEL1
  257. RETURN
  258. ENDIF
  259. C
  260. C Remplissage des types et du nom du nouveau MCHAML
  261. DO 135 J=1,MM
  262. MCHAM1.TYPCHE(J)='REAL*8'
  263. MCHAM1.NOMCHE(J)=NUANOM(IADD(NN+J))
  264. 135 CONTINUE
  265. C
  266. C on recupere le nombre d'element et le nombre de points par element
  267. MELVAL = MCHAML.IELVAL(1)
  268. C
  269. N1PTEL = VELCHE(/1)
  270. N1EL = VELCHE(/2)
  271.  
  272. C Creation des nouveaux MELVAL pour ecrire le resultat
  273. N2PTEL = 0
  274. N2EL = 0
  275. DO 138 J=1,MM
  276. SEGINI,MELVA1
  277. MCHAM1.IELVAL(J) = MELVA1
  278. 138 CONTINUE
  279.  
  280. N=N1EL*N1PTEL
  281. NT1= N / (100 * nbthrs)
  282.  
  283. C
  284. C Pour la paralelisation de l'interpolation
  285. C
  286. ITH = 0
  287. IF (NBESC .NE. 0) ith=oothrd
  288. IF ((NT1 .LE. 1) .OR. (NBTHRS .EQ. 1) .OR. (ITH .GT. 0)) THEN
  289. NBTHR = 1
  290. ITHRD = 0
  291. ELSE
  292. ithrd=1
  293. NBTHR = MIN(NT1,NBTHRS)
  294. call threadii
  295. ENDIF
  296.  
  297. SEGINI,MTR2
  298.  
  299. C Remplissage du 'COMMON/iplnuc' apres THREADII : pthread_mutex_lock
  300. C sinon soucis de cohabitation entre les ASSISTANTS qui ecrivent tous dans le meme COMMON...
  301. XP1 = XP
  302. NBTH1 = NBTHR
  303. INUA1 = INUAG
  304. ITR1 = MTR1
  305. IMAX1 = MAXI1
  306. IVAL1 = IVAL
  307. N1 = N
  308. NN1 = NN
  309. MM1 = MM
  310. NNMM1 = NNMM
  311. IMCHA1= MCHAML
  312. IMCHA2= MCHAM1
  313. N3EL = N1EL
  314. N3PTEL= N1PTEL
  315. ITR2 = MTR2
  316. XELIM1= XELIM
  317.  
  318. C Lancement des Threads
  319. if ((nbthr .gt. 1) .AND. (ithrd .eq. 1)) then
  320. do ith=2,nbthr
  321. call threadid(ith,ipln2i)
  322. enddo
  323. call ipln2i(1)
  324. do ith=2,nbthr
  325. call threadif(ith)
  326. enddo
  327. else
  328. call ipln2i(1)
  329. endif
  330.  
  331. if (ithrd .eq. 1) call threadis
  332.  
  333.  
  334. C SEGINI ,MTR2
  335. C ITR2=MTR2
  336.  
  337. C
  338. C boucle sur chaque point de gauss
  339. C
  340.  
  341. C DO 170 J=1,N1EL
  342. C DO 160 K=1,N1PTEL
  343. C
  344. C boucle sur les composantes pour remplir le tableau xi
  345. C qui contient les valeurs des variables maitres
  346. C DO 140 L=1,NN
  347. C MELVAL = IELVAL(L)
  348. C XI((J-1)*N1EL+N1PTEL,L) = VELCHE(K,J)
  349. C 140 CONTINUE
  350. C
  351. C on interpole le nuage pour les valeurs xi
  352. C
  353. C CALL IPLNU2(INUA,ITR1,ITR2,MAXI1,IVAL)
  354. C IF ( IERR .NE. 0 ) THEN
  355. C SEGDES MCHAML,MCHELM
  356. C DO 145 L=1,MCHAM1.IELVAL(/1)
  357. C MELVAL = MCHAM1.IELVAL(L)
  358. C SEGSUP MELVAL
  359. C 145 CONTINUE
  360. C SEGSUP MCHAM1,MCHEL1,MTR1
  361. C SEGSUP,MTR2
  362. C if (ithrd.eq.1) call threadis
  363. C RETURN
  364. C ENDIF
  365. C
  366. C remplissage des resultats stockes dans yi
  367. C
  368. C DO 150 L=1,MM
  369. C MELVAL = MCHAM1.IELVAL(L)
  370. C VELCHE(K,J) = YI((J-1)*N1EL+N1PTEL,L)
  371. C 150 CONTINUE
  372. C SEGSUP ,MTR2
  373. C
  374. C 160 CONTINUE
  375. C 170 CONTINUE
  376. C
  377.  
  378. SEGSUP ,MTR1,MTR2
  379. 200 CONTINUE
  380. C Fin de la boucle sur les MCHAML
  381.  
  382.  
  383. ICH1 = MCHEL1
  384. SEGSUP MAXI1
  385.  
  386. CALL ACTOBJ('MCHAML ',ICH1,1)
  387. CALL ECROBJ('MCHAML ',ICH1)
  388.  
  389. RETURN
  390.  
  391.  
  392.  
  393. 210 CONTINUE
  394. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  395. C Cas d'un CHPOINT C
  396. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  397.  
  398. MCHPOI = IPOI
  399.  
  400. SEGINI ,MCHPO1 = MCHPOI
  401. NSOUPO = IPCHP(/1)
  402. C
  403. C boucle sur les CHPOINTS elementaire
  404. DO 280 I=1,NSOUPO
  405. MSOUPO = IPCHP(I)
  406. C
  407. C calcul de nn et mm
  408. C
  409. NNMM=0
  410. SEGINI,MTR1
  411. ITR1 = MTR1
  412. NN = NOCOMP(/2)
  413. MM = NUANOM(/2) - NN
  414. C initialisation du nouveau soupo
  415. NC = MM
  416. SEGINI MSOUP1
  417. MCHPO1.IPCHP(I)=MSOUP1
  418. MSOUP1.IGEOC = IGEOC
  419. C
  420. C on verifie que les noms des composantes du CHPOINT sont
  421. C dans le nuage et on rempli iadd
  422. C
  423. DO 230 J=1,NOCOMP(/2)
  424. DO 220 K=1,NUANOM(/2)
  425. IF (NUANOM(K) .EQ. NOCOMP(J)) THEN
  426. IADD(**)=K
  427. NNMM = IADD(/1)
  428. ENDIF
  429. 220 CONTINUE
  430. IF ( NNMM .NE. J ) THEN
  431. C une des composantes du champ n'est pas un composante du nuage
  432. CALL ERREUR(682)
  433. C 'incoherence entre les composantes du nuage et du champ'
  434. SEGSUP MTR1,MCHPO1,MSOUP1
  435. RETURN
  436. ENDIF
  437. 230 CONTINUE
  438. C
  439. C on recupere les composantes du nuage qui ne sont pas dans le CHPOINT
  440. C
  441. DO 250 J=1,NUANOM(/2)
  442. IF ( NUATYP(J) .NE. 'FLOTTANT' ) THEN
  443. C la composante du nuanom(j) du nuage n'estpas scalaire
  444. MOTERR(1:8) = NUANOM(J)
  445. MOTERR(9:13) = 'NUAGE'
  446. CALL ERREUR(681)
  447. SEGSUP MTR1,MCHPO1,MSOUP1
  448. RETURN
  449. ENDIF
  450. FLAG = .TRUE.
  451. DO 240 K=1,NOCOMP(/2)
  452. IF (NUANOM(J) .EQ. NOCOMP(K)) THEN
  453. FLAG = .FALSE.
  454. ENDIF
  455. 240 CONTINUE
  456. IF (FLAG) THEN
  457. C la composante du nuage n'est pas dans le CHPOINT
  458. IADD(**) = J
  459. NNMM = NNMM +1
  460. ENDIF
  461. 250 CONTINUE
  462. IF ( NNMM .NE. NUANOM(/2)) THEN
  463. MOTERR='NUAGE'
  464. CALL ERREUR(980)
  465. RETURN
  466. ENDIF
  467. C
  468. C on rempli le noms des composantes du nouveau champ
  469. C
  470. DO 255 J=1,MM
  471. MSOUP1.NOCOMP(J)=NUANOM(IADD(J+NN))
  472. 255 CONTINUE
  473. MPOVAL = IPOVAL
  474.  
  475. C Creation du nouveau segment MPOVAL pour ecrire les resultats
  476. N = VPOCHA(/1)
  477. NC = MM
  478. SEGINI,MPOVA1
  479. MSOUP1.IPOVAL = MPOVA1
  480.  
  481.  
  482. C Preparation pour le calcul en parallele
  483. NT1= N / (100 * nbthrs)
  484. C
  485. C Pour la paralelisation de l'interpolation
  486. C
  487. ITH = 0
  488. IF (NBESC .NE. 0) ith=oothrd
  489. IF ((NT1 .LE. 1) .OR. (NBTHRS .EQ. 1) .OR. (ITH .GT. 0)) THEN
  490. NBTHR = 1
  491. ITHRD = 0
  492. ELSE
  493. ithrd=1
  494. NBTHR = MIN(NT1,NBTHRS)
  495. call threadii
  496. ENDIF
  497.  
  498. C Remplissage du 'COMMON/iplnuc' apres THREADII : pthread_mutex_lock
  499. C sinon soucis de cohabitation entre les ASSISTANTS qui ecrivent tous dans le meme COMMON...
  500. XP1 = XP
  501. NBTH1 = NBTHR
  502. INUA1 = INUAG
  503. ITR1 = MTR1
  504. IMAX1 = MAXI1
  505. IVAL1 = IVAL
  506. N1 = N
  507. NN1 = NN
  508. MM1 = MM
  509. NNMM1 = NNMM
  510. IMPOV1= MPOVAL
  511. IMPOV2= MPOVA1
  512. XELIM1= XELIM
  513.  
  514. C Lancement des Threads
  515. if ((nbthr.gt.1) .AND. (ithrd.eq.1)) then
  516. do ith=2,nbthr
  517. call threadid(ith,ipln2i)
  518. enddo
  519. call ipln2i(1)
  520. do ith=2,nbthr
  521. call threadif(ith)
  522. enddo
  523. else
  524. call ipln2i(1)
  525. endif
  526.  
  527. if (ithrd .eq. 1) call threadis
  528.  
  529.  
  530. C DO 270 J=1,N
  531. C Boucle sur les POINTS dont les valeurs sont a evaluer
  532.  
  533. C SEGINI,MTR2
  534. C ITR2 = MTR2
  535. C DO 260 K=1,NN
  536. C Boucle sur les composantes CONNUES du CHPOINT donne en argument
  537. C XI(K)=VPOCHA(J,K)
  538. C 260 CONTINUE
  539. C
  540. C Interpolation
  541. C
  542. C CALL IPLNU2(INUA,ITR1,ITR2,MAXI1,IVAL)
  543. CC en cas d'erreur on sort proprement
  544. C IF ( IERR .NE. 0 ) THEN
  545. C SEGDES MCHPOI,MSOUPO
  546. C SEGSUP MCHPO1,MSOUP1,MPOVAL,MTR1,MTR2,MAXI1
  547. C if (ithrd.eq.1) call threadis
  548. C RETURN
  549. C ENDIF
  550. C
  551. C DO 265 K=1,MM
  552. C Boucle sur les composantes MANQUANTES du CHPOINT resultat
  553. C MPOVA1.VPOCHA(J,K)=YI(K)
  554. C 265 CONTINUE
  555. C SEGSUP,MTR2
  556. C
  557. C 270 CONTINUE
  558.  
  559. C Desactivation des SEGMENTS
  560. SEGSUP,MTR1
  561. C
  562. 280 CONTINUE
  563. C
  564. IP1 = MCHPO1
  565. SEGSUP ,MAXI1
  566.  
  567. CALL ACTOBJ('CHPOINT ',IP1,1)
  568. CALL ECROBJ('CHPOINT ',IP1)
  569.  
  570. END
  571.  
  572.  
  573.  
  574.  
  575.  
  576.  
  577.  

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