Télécharger iplnu2.eso

Retour à la liste

Numérotation des lignes :

  1. C IPLNU2 SOURCE FD218221 16/08/04 21:15:03 9046
  2.  
  3. SUBROUTINE IPLNU2(INUA,ITR1,XI,YI,MAXI1,IVAL,
  4. & N,INDJ,NN,MM,NNMM,XP)
  5. C
  6. C Fonction
  7. C Cette sous routine interpole une fonction de n variables dont
  8. C certains couples (x,f(x)) on ete stockes dans un objet de type
  9. C nuage.
  10. C 2 techniques sont disponibles :
  11. C - Pour IVAL = 1 ou 2, on utilise une methode des elements diffus.
  12. C En chaque point on cherche la forme lineaire qui minimise la distance
  13. C ponderee de l'hyperplan plan au point du nuage. La valeur obtenue est
  14. C celle de la forme lineaire au point considere.
  15. C - Pour IVAL = 3, on utilise une interpolation exacte par ponderation
  16. C inverse a la distance. Le poids attribue a chaque point du nuage est
  17. C en 1/D**XP ou D est sa distance au point cible et XP le parametre
  18. C puissance fourni par l'utilisateur.
  19. C Variables
  20. C INUA pointeur sur l'objet de type nuage
  21. C NNMM nombe de composantes du nuage
  22. C NN nombre de composantes connues du point cible
  23. C MM nombre de composantes a calculer du point cible
  24. C NBCOUP nombe de points (cad de NNMM-uplets) du nuage
  25. C ITR1.IADD contient la correspondance entre les numeros des composantes
  26. C du CHPOINT/MCHAML et celles du nuage :
  27. C - La composante connue I (entre 1 et NN) du champ est a la position
  28. C IADD(I) dans le nuage
  29. C - La composante a calculer J (entre 1 et MM) est a la position
  30. C IADD(J) dans le nuage
  31. C MAXI1.MAXV contient les max des composantes du nuage (utilise pour normer
  32. C le calcul des distances)
  33. C X tableau contenant les valeurs des NN composantes connues du CHPOINT/MCHAML
  34. C X(i,j) est la valeur de la composante j au i-eme point support du
  35. C champ, j est compris entre 1 et NN, la parallelisation est faite sur i
  36. C Y tableau contenant les valeurs des MM composantes interpolees, c'est
  37. C le resultat de l'interpolation
  38. C Y(i,j) est la valeur de la composante j au i-eme point support du
  39. C champ, j est compris entre 1 et MM, la parallelisation est faite sur i
  40. C IVAL est l'entier indiquant l'option choisie
  41. C IVAL=1 : ponderation par fonction gaussienne EXP(-Z**2) (valeur par defaut)
  42. C IVAL=2 : ponderation par fonction rationnelle 1/(1+Z)
  43. C IVAL=3 : interpolation exacte par ponderation inverse a la distance
  44. C
  45. C Auteur A De Gayffier
  46. C Date 05/07/94
  47. C
  48. C Langage esope+fortran 77
  49. C
  50. IMPLICIT INTEGER(I-N)
  51. IMPLICIT REAL*8(A-H,O-Z)
  52. PARAMETER ( GRAND = 1.D50)
  53. -INC CCOPTIO
  54. -INC SMNUAGE
  55. SEGMENT MTR1
  56. INTEGER IADD(NNMM)
  57. ENDSEGMENT
  58.  
  59. REAL*8 XI,YI
  60. DIMENSION XI(N,*)
  61. DIMENSION YI(N,*)
  62.  
  63. C Ces tableaux contienent la matrice du syteme lineaire
  64. C A contient les termes et INDX contient les indices de permutation
  65. C une fois la decomposition effectuee
  66. REAL*8 A(NN+1,NN+1),D
  67. INTEGER INDX(NN+1)
  68.  
  69. REAL*8 B(NN+1,MM)
  70.  
  71. REAL*8 DX(NN+1)
  72.  
  73. REAL*8 DIMINI(NN+1)
  74. INTEGER NNBCP(NN+1)
  75.  
  76. SEGMENT MAXI1
  77. REAL*8 MAXV(NNMM)
  78. ENDSEGMENT
  79.  
  80. C
  81. MNUAGE = INUA
  82. NUAVFL = NUAPOI(1)
  83. NBCOUP = NUAFLO(/1)
  84.  
  85. C (FDP 2016) Nouvelle option 'PID' : interpolation exacte par ponderation
  86. C inverse a la distance
  87. IF (IVAL.EQ.3) THEN
  88. CALL IPPID(INUA,NBCOUP,ITR1,XI,YI,N,INDJ,NN,MM,NNMM,XP)
  89. RETURN
  90. ENDIF
  91.  
  92.  
  93. C INITIALISATIONS A ZERO
  94. DO INDICE=1,NN+1
  95. INDX(INDICE) =0
  96. DIMINI(INDICE)=0.D0
  97. NNBCP(INDICE) =0
  98. DX(INDICE) =0
  99.  
  100. DO JNDICE=1,NN+1
  101. A(INDICE,JNDICE)=0.D0
  102. ENDDO
  103.  
  104. DO JNDICE=1,MM
  105. B(INDICE,JNDICE)=0.D0
  106. ENDDO
  107. ENDDO
  108.  
  109. C
  110. IF (NBCOUP .EQ. 0) THEN
  111. C nuage vide operation impossible
  112. CALL ERREUR(680)
  113. RETURN
  114. ENDIF
  115. C
  116. MTR1 = ITR1
  117. C initialisation du tableau pour stocker les distances
  118. C on place des valeurs tres grandes
  119. DO 1 I=1,NN+1
  120. DIMINI(I)= GRAND
  121. 1 CONTINUE
  122. C
  123. C calcul de la distance cararcteristique du nuage
  124. C
  125. DO 60 I=1,NBCOUP
  126. SUM = 0.D0
  127. DO 30 J=1,NN
  128. NUAVFL = NUAPOI(IADD(J))
  129. C on divise les distances par les maximums de chaque composante
  130. DUM = (XI(INDJ,J)-NUAFLO(I))/MAXV(IADD(J))
  131. SUM = SUM + DUM*DUM
  132. 30 CONTINUE
  133. DIST = SQRT(SUM)
  134. C stockage des plus petites valeurs
  135. DO 50 J=1,NN+1
  136. IF ( DIST .LT. DIMINI(J) )THEN
  137. DO 40 K=NN+1,J+1,-1
  138. DIMINI(K)=DIMINI(K-1)
  139. NNBCP(K)=NNBCP(K-1)
  140. 40 CONTINUE
  141. DIMINI(J)=DIST
  142. NNBCP(J)=I
  143. GOTO 60
  144. ENDIF
  145. 50 CONTINUE
  146. 60 CONTINUE
  147. C
  148. DCARA =0.D0
  149. DO 65 J=1,NN+1
  150. DCARA = DCARA + DIMINI(J)
  151. 65 CONTINUE
  152. DCARA = DCARA / (NN+1.D0)
  153. C
  154. C affichage de mise au point avec opti impi 9022 ;
  155. C
  156. IF (IIMPI .EQ. 9022) WRITE(6,1001) DCARA
  157. IF (IIMPI .EQ. 9022) THEN
  158. WRITE(ioimp,1002) (NN+1)
  159. write(ioimp,1005) (DIMINI(I),I=1,NN+1)
  160. C boucle sur les composantes
  161. DO 67 I=1,NN
  162. NUAVFL = NUAPOI(IADD(I))
  163. WRITE(ioimp,1004) NUANOM(IADD(I))
  164. write(ioimp,1003) (NUAFLO(NNBCP(J))/MAXV(IADD(I)),J=1,NN+1)
  165. WRITE(ioimp,*)
  166. 67 CONTINUE
  167. DO 68 I=NN+1,NNMM
  168. NUAVFL = NUAPOI(IADD(I))
  169. WRITE(ioimp,1004) NUANOM(IADD(I))
  170. write(ioimp,1003) (NUAFLO(NNBCP(J))/MAXV(IADD(I)),J=1,NN+1)
  171. WRITE(ioimp,*)
  172. 68 CONTINUE
  173. ENDIF
  174. 1001 FORMAT('La distance caracteristique vaut ',G10.4)
  175. 1002 FORMAT('Distance et coordonees des',1X,I2,1X,'points
  176. & les plus proches sont')
  177. 1003 FORMAT(30(G10.4,1X))
  178. 1004 FORMAT('Composante ',A)
  179. 1005 FORMAT('Distance au point courant:',/,30(G10.4))
  180. C
  181. C boucle sur les points du nuage (nbcoup)
  182. C
  183. DO 110 I=1,NBCOUP
  184. C
  185. C boucle sur les monomes du polynome d'interpolation (1,x,y,...)
  186. C remplissage complet de DX (segment MDX)
  187. DX(1)=1.D0
  188. SUM = 0.D0
  189. DO 70 J=2,(NN+1)
  190. NUAVFL = NUAPOI(IADD(J-1))
  191. C on divise les distances par les maximums de chaque composante
  192. DX(J)=(XI(INDJ,J-1)-NUAFLO(I))/MAXV(IADD(J-1))
  193. SUM = SUM + DX(J)*DX(J)
  194. 70 CONTINUE
  195. DIST = SQRT(SUM)
  196. C
  197. C fonction de ponderation de la distance
  198. Z = DIST/DCARA
  199. IF (IVAL .EQ. 1 ) THEN
  200. WDED = EXP(-Z*Z)
  201. ELSE
  202. WDED = 1.D0/(1.D0+Z)
  203. ENDIF
  204. C remplissage de a et de b
  205. DO 100 J=1,(NN+1)
  206. r_z = WDED*DX(J)
  207. DO 80 K=1,(NN+1)
  208. C* A(J,K) = A(J,K) + WDED*DX(J) * DX(K)
  209. A(J,K) = A(J,K) + r_z * DX(K)
  210. 80 CONTINUE
  211. DO 90 K=1,MM
  212. NUAVFL = NUAPOI(IADD(NN+K))
  213. C* B(J,K)=B(J,K) + WDED*DX(J)*NUAFLO(I)
  214. B(J,K)=B(J,K) + r_z * NUAFLO(I)
  215. 90 CONTINUE
  216. C
  217. 100 CONTINUE
  218. C
  219. 110 CONTINUE
  220. C
  221. C SEGSUP,MDX
  222. C
  223. C phase de resolution
  224. C
  225. C on resoud a.x=b pp fois
  226. C
  227. CALL IPLNU3(A,D,INDX,NN+1)
  228.  
  229. DO 120 I=1,MM
  230. CALL IPLNU4(A,D,INDX,NN+1,B,MM,I)
  231. 120 CONTINUE
  232. C
  233. C on replace le resultat dans mtr2
  234. C dans les B(*,I) sont stockees les derivees partielles de la fonction
  235. C
  236. DO 130 I=1,MM
  237. YI(INDJ,I)=B(1,I)
  238. 130 CONTINUE
  239. C
  240. C
  241. RETURN
  242. END
  243.  
  244.  

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