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

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