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.  
  54. -INC PPARAM
  55. -INC CCOPTIO
  56. -INC SMNUAGE
  57. SEGMENT MTR1
  58. INTEGER IADD(NNMM)
  59. ENDSEGMENT
  60.  
  61. REAL*8 XI,YI
  62. DIMENSION XI(N,*)
  63. DIMENSION YI(N,*)
  64.  
  65. C Ces tableaux contienent la matrice du syteme lineaire
  66. C A contient les termes et INDX contient les indices de permutation
  67. C une fois la decomposition effectuee
  68. REAL*8 A(NN+1,NN+1),D
  69. INTEGER INDX(NN+1)
  70.  
  71. REAL*8 B(NN+1,MM)
  72.  
  73. REAL*8 DX(NN+1)
  74.  
  75. REAL*8 DIMINI(NN+1)
  76. INTEGER NNBCP(NN+1)
  77.  
  78. SEGMENT MAXI1
  79. REAL*8 MAXV(NNMM)
  80. ENDSEGMENT
  81.  
  82. C
  83. MNUAGE = INUA
  84. NUAVFL = NUAPOI(1)
  85. NBCOUP = NUAFLO(/1)
  86.  
  87. C (FDP 2016) Nouvelle option 'PID' : interpolation exacte par ponderation
  88. C inverse a la distance
  89. IF (IVAL.EQ.3) THEN
  90. CALL IPPID(INUA,NBCOUP,ITR1,XI,YI,N,INDJ,NN,MM,NNMM,XP)
  91. RETURN
  92. ENDIF
  93.  
  94.  
  95. C INITIALISATIONS A ZERO
  96. DO INDICE=1,NN+1
  97. INDX(INDICE) =0
  98. DIMINI(INDICE)=0.D0
  99. NNBCP(INDICE) =0
  100. DX(INDICE) =0
  101.  
  102. DO JNDICE=1,NN+1
  103. A(INDICE,JNDICE)=0.D0
  104. ENDDO
  105.  
  106. DO JNDICE=1,MM
  107. B(INDICE,JNDICE)=0.D0
  108. ENDDO
  109. ENDDO
  110.  
  111. C
  112. IF (NBCOUP .EQ. 0) THEN
  113. C nuage vide operation impossible
  114. CALL ERREUR(680)
  115. RETURN
  116. ENDIF
  117. C
  118. MTR1 = ITR1
  119. C initialisation du tableau pour stocker les distances
  120. C on place des valeurs tres grandes
  121. DO 1 I=1,NN+1
  122. DIMINI(I)= GRAND
  123. 1 CONTINUE
  124. C
  125. C calcul de la distance cararcteristique du nuage
  126. C
  127. DO 60 I=1,NBCOUP
  128. SUM = 0.D0
  129. DO 30 J=1,NN
  130. NUAVFL = NUAPOI(IADD(J))
  131. C on divise les distances par les maximums de chaque composante
  132. DUM = (XI(INDJ,J)-NUAFLO(I))/MAXV(IADD(J))
  133. SUM = SUM + DUM*DUM
  134. 30 CONTINUE
  135. DIST = SQRT(SUM)
  136. C stockage des plus petites valeurs
  137. DO 50 J=1,NN+1
  138. IF ( DIST .LT. DIMINI(J) )THEN
  139. DO 40 K=NN+1,J+1,-1
  140. DIMINI(K)=DIMINI(K-1)
  141. NNBCP(K)=NNBCP(K-1)
  142. 40 CONTINUE
  143. DIMINI(J)=DIST
  144. NNBCP(J)=I
  145. GOTO 60
  146. ENDIF
  147. 50 CONTINUE
  148. 60 CONTINUE
  149. C
  150. DCARA =0.D0
  151. DO 65 J=1,NN+1
  152. DCARA = DCARA + DIMINI(J)
  153. 65 CONTINUE
  154. DCARA = DCARA / (NN+1.D0)
  155. C
  156. C affichage de mise au point avec opti impi 9022 ;
  157. C
  158. IF (IIMPI .EQ. 9022) WRITE(6,1001) DCARA
  159. IF (IIMPI .EQ. 9022) THEN
  160. WRITE(ioimp,1002) (NN+1)
  161. write(ioimp,1005) (DIMINI(I),I=1,NN+1)
  162. C boucle sur les composantes
  163. DO 67 I=1,NN
  164. NUAVFL = NUAPOI(IADD(I))
  165. WRITE(ioimp,1004) NUANOM(IADD(I))
  166. write(ioimp,1003) (NUAFLO(NNBCP(J))/MAXV(IADD(I)),J=1,NN+1)
  167. WRITE(ioimp,*)
  168. 67 CONTINUE
  169. DO 68 I=NN+1,NNMM
  170. NUAVFL = NUAPOI(IADD(I))
  171. WRITE(ioimp,1004) NUANOM(IADD(I))
  172. write(ioimp,1003) (NUAFLO(NNBCP(J))/MAXV(IADD(I)),J=1,NN+1)
  173. WRITE(ioimp,*)
  174. 68 CONTINUE
  175. ENDIF
  176. 1001 FORMAT('La distance caracteristique vaut ',G10.4)
  177. 1002 FORMAT('Distance et coordonees des',1X,I2,1X,'points
  178. & les plus proches sont')
  179. 1003 FORMAT(30(G10.4,1X))
  180. 1004 FORMAT('Composante ',A)
  181. 1005 FORMAT('Distance au point courant:',/,30(G10.4))
  182. C
  183. C boucle sur les points du nuage (nbcoup)
  184. C
  185. DO 110 I=1,NBCOUP
  186. C
  187. C boucle sur les monomes du polynome d'interpolation (1,x,y,...)
  188. C remplissage complet de DX (segment MDX)
  189. DX(1)=1.D0
  190. SUM = 0.D0
  191. DO 70 J=2,(NN+1)
  192. NUAVFL = NUAPOI(IADD(J-1))
  193. C on divise les distances par les maximums de chaque composante
  194. DX(J)=(XI(INDJ,J-1)-NUAFLO(I))/MAXV(IADD(J-1))
  195. SUM = SUM + DX(J)*DX(J)
  196. 70 CONTINUE
  197. DIST = SQRT(SUM)
  198. C
  199. C fonction de ponderation de la distance
  200. Z = DIST/DCARA
  201. IF (IVAL .EQ. 1 ) THEN
  202. WDED = EXP(-Z*Z)
  203. ELSE
  204. WDED = 1.D0/(1.D0+Z)
  205. ENDIF
  206. C remplissage de a et de b
  207. DO 100 J=1,(NN+1)
  208. r_z = WDED*DX(J)
  209. DO 80 K=1,(NN+1)
  210. C* A(J,K) = A(J,K) + WDED*DX(J) * DX(K)
  211. A(J,K) = A(J,K) + r_z * DX(K)
  212. 80 CONTINUE
  213. DO 90 K=1,MM
  214. NUAVFL = NUAPOI(IADD(NN+K))
  215. C* B(J,K)=B(J,K) + WDED*DX(J)*NUAFLO(I)
  216. B(J,K)=B(J,K) + r_z * NUAFLO(I)
  217. 90 CONTINUE
  218. C
  219. 100 CONTINUE
  220. C
  221. 110 CONTINUE
  222. C
  223. C SEGSUP,MDX
  224. C
  225. C phase de resolution
  226. C
  227. C on resoud a.x=b pp fois
  228. C
  229. CALL IPLNU3(A,D,INDX,NN+1)
  230.  
  231. DO 120 I=1,MM
  232. CALL IPLNU4(A,D,INDX,NN+1,B,MM,I)
  233. 120 CONTINUE
  234. C
  235. C on replace le resultat dans mtr2
  236. C dans les B(*,I) sont stockees les derivees partielles de la fonction
  237. C
  238. DO 130 I=1,MM
  239. YI(INDJ,I)=B(1,I)
  240. 130 CONTINUE
  241. C
  242. C
  243. RETURN
  244. END
  245.  
  246.  

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