Télécharger iplnu1.eso

Retour à la liste

Numérotation des lignes :

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

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