Télécharger iplnu1.eso

Retour à la liste

Numérotation des lignes :

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

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