Esope source files

   1 : C IPLNU2    SOURCE    FANDEUR   11/07/19    21:16:12     7039           
   2 : 
   3 :         SUBROUTINE IPLNU2(INUA,ITR1,ITR2,MAX1,IVAL)
   4 : c
   5 : c fonction
   6 : c     cette sous routine interpole une fonction de n variables dont
   7 : c     certains couple (x,f(x)) on ete stockés dans un objet de type
   8 : c     nuage. On utilise une methode des elements diffus. En chaque 
   9 : c     point on cheche la forme lineaire qui minimise la distance ponderée 
  10 : c     de l'hyperplan plan au point du nuage. La valeur obtenue est 
  11 : c     celle de la forme lineaire au point considere.
  12 : c
  13 : c variables
  14 : c     inua pointeur sur l'objet de type nuage
  15 : c     itr1 pointeur sur le segment mtr1 de iplnu2.eso
  16 : c                   qui contient une correspondance entre la position
  17 : c                   des variables dans les tableua du nuage
  18 : c     itr2  pointeur sur le segment mtr2 qui contient les variables x
  19 : c                    dans xi et les valeurs f(x) dans yi
  20 : c     max1  pointeur sur l'objet qui contient le max des composantes du nuage
  21 : c
  22 : c appelé par iplnu2.eso
  23 : c date 05/07/94
  24 : c auteur A de Gayffier 
  25 : c langage fortran 77 + esope
  26 : c                                
  27 :       IMPLICIT INTEGER(I-N)
  28 :       IMPLICIT    REAL*8(A-H,O-Z) 
  29 :       PARAMETER ( GRAND = 1.D50) 
  30 : -INC CCOPTIO
  31 : -INC SMNUAGE
  32 :       SEGMENT MTR1
  33 :          INTEGER IADD(NNMM)
  34 :       ENDSEGMENT
  35 :       SEGMENT MTR2
  36 :          REAL*8  XI(NN),YI(MM)
  37 :       ENDSEGMENT
  38 :       SEGMENT MA 
  39 : *     ce segment contient la matrice du syteme lineaire
  40 : *     a contient les termes et indx contient les indices de permutation
  41 : *     une fois la decomposition effectuée 
  42 :          REAL*8 A(NN+1,NN+1),D
  43 :          INTEGER INDX(NN+1)
  44 :       ENDSEGMENT
  45 :       SEGMENT MB 
  46 :          REAL*8 B(NN+1,MM)
  47 :       ENDSEGMENT
  48 :       SEGMENT MDX
  49 :           REAL*8 DX(NN+1)
  50 :       ENDSEGMENT
  51 :       SEGMENT MDIS
  52 :           REAL*8 DIMINI(NN+1)
  53 :           INTEGER NNBCP(NN+1)
  54 :       ENDSEGMENT
  55 :       SEGMENT MAX1
  56 :           REAL*8 MAXV(NNMM)
  57 :       ENDSEGMENT
  58 : c
  59 : c   on active tous les segments necessaires
  60 : c
  61 :       MNUAGE = INUA   
  62 : c      SEGACT ,MNUAGE
  63 : c
  64 :       NVAR   = NUAPOI(/1)
  65 :       DO 10 I=1,NVAR 
  66 :         NUAVFL = NUAPOI(I)
  67 :         SEGACT ,NUAVFL
  68 :  10   CONTINUE
  69 : *
  70 :       NBCOUP = NUAFLO(/1)
  71 : *
  72 :       IF (NBCOUP .EQ. 0) THEN
  73 : *    nuage vide operation impossible
  73 :         CALL ERREUR(680)
  74 :         SEGDES ,MNUAGE
  75 :         RETURN
  76 :       ENDIF
  77 : *
  78 :       MTR1 = ITR1
  79 :       MTR2 = ITR2
  80 :       SEGACT MTR1,MTR2*MOD
  81 : c      SEGACT ,MAX1
  82 :       NN = XI(/1)
  83 :       MM = YI(/1)
  84 :       NNMM = IADD(/1)
  85 : c
  86 : c      IF ( NVAR .NE. (NN+MM)) THEN
  87 : c        write(ioimp,*) 'ERREUR DANS LA DIMENSION DES MTR'
  88 : c      ENDIF
  89 : c
  90 :       SEGINI ,MA,MB,MDIS
  91 : c    initialisation du tableau pour stocker les distances
  92 : *    on place des valeurs tres grandes
  93 :       DO 1 I=1,NN+1
  94 :         DIMINI(I)= GRAND
  95 :  1    CONTINUE
  96 : c
  97 : c   calcul de la distance cararctéristique du nuage
  98 : c
  99 :       DO 60 I=1,NBCOUP
 100 :         SUM = 0.D0
 101 :         DO 30 J=1,NN  
 102 :            NUAVFL = NUAPOI(IADD(J))
 103 : c        on divise les distances par les maximums de chaque composante
 104 :            DUM = (XI(J)-NUAFLO(I))/MAXV(IADD(J))
 105 :            SUM = SUM + DUM*DUM
 106 :  30     CONTINUE
 107 :         DIST = SQRT(SUM)
 108 : c   stockage des plus petites valeurs
 109 :         DO 50 J=1,NN+1
 110 :            IF ( DIST .LT. DIMINI(J) )THEN
 111 :              DO 40 K=NN+1,J+1,-1
 112 :                DIMINI(K)=DIMINI(K-1)
 113 :                NNBCP(K)=NNBCP(K-1)
 114 :  40          CONTINUE
 115 :              DIMINI(J)=DIST
 116 :              NNBCP(J)=I
 117 :              GOTO 60
 118 :            ENDIF
 119 :  50     CONTINUE
 120 :  60   CONTINUE
 121 : c
 122 :       DCARA =0.D0
 123 :       DO 65 J=1,NN+1
 124 :         DCARA = DCARA + DIMINI(J)
 125 :   65  CONTINUE
 126 :       DCARA = DCARA / (NN+1.D0)
 127 : c
 128 : c   affichage de mise au point avec opti impi 9022 ;
 129 : c
 130 :       IF (IIMPI .EQ. 9022) WRITE(6,1001) DCARA
 131 :       IF (IIMPI .EQ. 9022) THEN
 132 :            WRITE(ioimp,1002) (NN+1)
 133 :            write(ioimp,1005) (DIMINI(I),I=1,NN+1)
 134 : *               boucle sur les composantes
 135 :            DO 67 I=1,NN 
 136 :              NUAVFL = NUAPOI(IADD(I))              
 137 :              WRITE(ioimp,1004) NUANOM(IADD(I))
 138 :              write(ioimp,1003) (NUAFLO(NNBCP(J))/MAXV(IADD(I)),J=1,NN+1)
 139 :              WRITE(ioimp,*)
 140 :    67      CONTINUE
 141 :            DO 68 I=NN+1,NNMM
 142 :              NUAVFL = NUAPOI(IADD(I))              
 143 :              WRITE(ioimp,1004) NUANOM(IADD(I))
 144 :              write(ioimp,1003) (NUAFLO(NNBCP(J))/MAXV(IADD(I)),J=1,NN+1)
 145 :              WRITE(ioimp,*)
 146 :    68      CONTINUE
 147 :       ENDIF  
 148 :  1001 FORMAT('La distance caractéristique vaut ',G10.4)
 149 :  1002 FORMAT('Distance et coordonées des',1X,I2,1X,'points 
 150 :      &        les plus proches sont')
 151 :  1003 FORMAT(30(G10.4,1X))
 152 :  1004 FORMAT('Composante ',A)
 153 :  1005 FORMAT('Distance au point courant:',/,30(G10.4))
 154 : *
 155 :       SEGSUP ,MDIS
 156 : c
 157 : c On sort l'initialisation de ce segment de la boucle
 158 :       SEGINI MDX
 159 : c
 160 : c   boucle sur les points du nuage (nbcoup)
 161 : c
 162 :       DO 110 I=1,NBCOUP
 163 : c  
 164 : c     boucle sur les monomes du polynome d'interpolation (1,x,y,...)
 165 : c  remplissage complet de DX (segment MDX)
 166 :         DX(1)=1.D0
 167 :         SUM = 0.D0
 168 :         DO 70 J=2,(NN+1)  
 169 :            NUAVFL = NUAPOI(IADD(J-1))  
 170 : c        on divise les distances par les maximums de chaque composante
 171 :            DX(J)=(XI(J-1)-NUAFLO(I))/MAXV(IADD(J-1))
 172 :            SUM = SUM + DX(J)*DX(J)
 173 :  70     CONTINUE
 174 :         DIST = SQRT(SUM)
 175 : C
 176 : *  fonction de ponderation de la distance
 177 :         Z = DIST/DCARA
 178 :         IF (IVAL .EQ. 1 ) THEN
 179 :           WDED = 1.D0/(Z+1.D0)
 180 :         ELSE
 181 :           WDED = EXP(-Z*Z)
 182 :         ENDIF
 183 : c       remplissage de a et de b
 184 :         DO 100 J=1,(NN+1)
 185 :           r_z = WDED*DX(J)
 186 :           DO 80 K=1,(NN+1)
 187 : C*            A(J,K) = A(J,K) + WDED*DX(J) * DX(K)
 188 :               A(J,K) = A(J,K) + r_z * DX(K)
 189 :  80       CONTINUE
 190 :           DO 90 K=1,MM
 191 :               NUAVFL = NUAPOI(IADD(NN+K))
 192 : C*            B(J,K)=B(J,K) + WDED*DX(J)*NUAFLO(I)
 193 :               B(J,K)=B(J,K) + r_z * NUAFLO(I)
 194 :  90      CONTINUE
 195 : c
 196 :  100    CONTINUE
 197 : c
 198 :  110  CONTINUE
 199 : c
 200 :       SEGSUP,MDX
 201 : c
 202 : c  phase de resolution
 203 : c
 204 : c  on resoud a.x=b pp fois
 205 : c
 205 :        CALL IPLNU3(MA)
 206 : 
 207 :        DO 120 I=1,MM
 207 :          CALL IPLNU4(MA,MB,I)
 208 :  120   CONTINUE        
 209 : c
 210 : c   on replace le resultat dans mtr2
 211 : c   dans les B(*,I) sont stockees les dérivées partielles de la fonction
 212 : c
 213 :         DO 130 I=1,MM
 214 :           YI(I)=B(1,I)
 215 :  130    CONTINUE
 216 : c
 217 : c  on ferme et on détruit les segments activés
 218 : c
 219 :         DO 140 I=1,NVAR 
 220 :           NUAVFL = NUAPOI(I)
 221 :           SEGDES ,NUAVFL
 222 :  140    CONTINUE
 223 : c
 224 : c        SEGDES MNUAGE
 225 :         SEGSUP MA,MB
 226 : c
 227 :         RETURN
 228 :         END
 229 : 
 230 :