Télécharger la plaquette Cast3M


fr en

Télécharger agreg1.eso

Retour à la liste

Numérotation des lignes : 

agreg1

  1. C AGREG1    SOURCE    FD218221  25/12/24    21:15:01     12435          
  2.       SUBROUTINE AGREG1(MLREE1,ICAS,XP,IDERI,IROBU,AGR,MLREE2)
  3.  
  4. C     Application d'une fonction d'agregation sur un
  5. C     objet LISTREEL
  6.  
  7. C     Entrees :
  8. C     ---------
  9. C       MLREE1 : Pointeur vers le LISTREEL (suppose actif)
  10. C       ICAS   : Entier, code de la fonction
  11. C                = 1  'SOMM'  Somme
  12. C                = 2  'PROD'  Produit
  13. C                = 3  'MOYE'  Moyenne arithmetique (moment d'ordre 1)
  14. C                = 4  'MOHA'  Moyenne harmonique
  15. C                = 5  'MOGE'  Moyenne geometrique
  16. C                = 6  'VARI'  Variance (moment centre d'ordre 2)
  17. C                = 7  'ECTY'  Ecart type
  18. C                = 8  'ASYM'  Coefficient d'asymetrie (moment centre reduit d'ordre 3)
  19. C                = 9  'KURT'  Kurtosis (moment centre reduit d'ordre 4)
  20. C                = 10 'MEDI'  Mediane
  21. C                = 11 'PMOM'  Moment d'ordre P
  22. C                = 12 'PMOY'  Moyenne generalise d'ordre P
  23. C                = 13 'PNOR'  Norme generalisee d'ordre P
  24. C                = 14 'LEHM'  Fonction de Lehmer d'ordre P
  25. C                = 15 'KSL'   Fonction de Kreisselmeir Steinhauser inferieure d'ordre P (MellowMax)
  26. C                = 16 'KSU'   Fonction de Kreisselmeir Steinhauser superieure d'ordre P (LogSumExp)
  27. C                = 17 'BOLT'  Fonction de Boltzmann d'ordre P
  28. C       XP    : Flottant, parametre pour les fonctions 'PMOY' 'PMOM' 'PNOR' 'LEHM'
  29. C                                                      'KSL'  'KSL'  'BOLT'
  30. C       IDERI : Entier, pour calcul des derivees partielles
  31. C                = 0 pas de calcul des derivees
  32. C               != 0 calcul des derivees partielles
  33. C       IROBU : Entier, pour calcul robuste au overflow
  34. C                = 0 pour un calcul "naif"
  35. C               != 0 pour un calcul "robuste" en normalisant les valeurs avec
  36. C                    la norme infinie ou bien le maximum, selon la fonction choisie
  37.  
  38. C     Sorties :
  39. C     ---------
  40. C       AGR    : Flottant, valeur de la fonction d'agregation
  41. C       MLREE2 : Pointeur vers le LISTREEL, de meme dimentsion que MLREE1, contenant les valeurs des
  42. C                derivees partielles de la fonction d'agregation par rapport a chaque terme de MLREE1
  43.  
  44.  
  45. C     Typages implicites habituels
  46.       IMPLICIT INTEGER(I-N)
  47.       IMPLICIT REAL*8(A-H,O-Z)
  48.  
  49. C     Les includes necessaires
  50. -INC CCREEL
  51. -INC SMLREEL
  52.  
  53. C     Quelques objets
  54.       PARAMETER(UN=1.D0)
  55.       LOGICAL ROBU,DERI
  56.  
  57. C     Taille du LISTREEL
  58.       NX=MLREE1.PROG(/1)
  59.  
  60. C     Initialisation du resultat
  61.       AGR=XZERO
  62.       MLREE2=0
  63.       IF (IDERI.EQ.1) THEN
  64.         JG=NX
  65.         SEGINI MLREE2
  66.       ENDIF
  67.  
  68. C     Cas trivial non traite
  69.       IF (NX.LT.1) THEN
  70.         CALL ERREUR(21)
  71.         RETURN
  72.       ENDIF
  73.  
  74. C     Calcul robuste ?
  75.       ROBU=.FALSE.
  76.       IF (IROBU.NE.0) THEN
  77.         ROBU=.TRUE.
  78.         IF (ICAS.GE.12) THEN
  79.           CALL MAXIN3(MLREE1,IMAX,VMAX,1,0)
  80.           VINF=ABS(VMAX)
  81.           VPMAX=XP*VMAX
  82.         ENDIF
  83.       ENDIF
  84.  
  85. C     Calcul des derivees partielles ?
  86.       DERI=.FALSE.
  87.       IF (IDERI.NE.0) THEN
  88.         DERI=.TRUE.
  89.       ENDIF
  90.  
  91. C     Cas de la P moyenne avec P=0 --> on utilise MOGE (moyenne geometrique)
  92.       IF ((ICAS.EQ.12).AND.(ABS(XP).LT.XPETIT)) THEN
  93.         ICAS=5
  94.       ENDIF
  95.  
  96. C     1. Somme (aussi utilise pour la moyenne, la variance, l'ecart type,
  97. C               l'asymetrie, le kurtosis)
  98.       IF ((ICAS.EQ.1).OR.(ICAS.EQ.3).OR.(ICAS.EQ.6).OR.(ICAS.EQ.7).OR.
  99.      &    (ICAS.EQ.8).OR.(ICAS.EQ.9)) THEN
  100.         SUM=XZERO
  101.         DO I=1,NX
  102.           XI=MLREE1.PROG(I)
  103.           SUM=SUM+XI
  104.         ENDDO
  105.         IF (ICAS.EQ.1) THEN
  106.           AGR=SUM
  107.           IF (DERI) THEN
  108.             MLREE2.PROG(I)=UN
  109.           ENDIF
  110.         ENDIF
  111.       ENDIF
  112.  
  113. C     2. Produit (aussi utilise pour la moyenne geometrique)
  114.       IF ((ICAS.EQ.2).OR.(ICAS.EQ.5)) THEN
  115.         XPRO=UN
  116.         DO I=1,NX
  117.           XI=MLREE1.PROG(I)
  118.           XPRO=XPRO*XI
  119.         ENDDO
  120.         IF (DERI) THEN
  121.           DO I=1,NX
  122.             XI=MLREE1.PROG(I)
  123.             IF (ABS(XI).LT.XPETIT) THEN
  124.               MLREE2.PROG(I)=XZERO
  125.             ELSE
  126.               MLREE2.PROG(I)=XPRO/XI
  127.             ENDIF
  128.           ENDDO
  129.         ENDIF
  130.         IF (ICAS.EQ.2) AGR=XPRO
  131.       ENDIF
  132.  
  133. C     3. Moyenne (aussi utilise pour la variance, l'ecart type, l'asymetrie, le kurtosis)
  134.       IF ((ICAS.EQ.3).OR.(ICAS.EQ.6).OR.(ICAS.EQ.7).OR.(ICAS.EQ.8).OR.
  135.      &    (ICAS.EQ.9)) THEN
  136.         XMOY=SUM/NX
  137.         IF (ICAS.EQ.3) THEN
  138.           AGR=XMOY
  139.           IF (DERI) THEN
  140.             DO I=1,NX
  141.               IF (NX.EQ.0) THEN
  142.                 MLREE2.PROG(I)=XZERO
  143.               ELSE
  144.                 MLREE2.PROG(I)=UN/NX
  145.               ENDIF
  146.             ENDDO
  147.           ENDIF
  148.         ENDIF
  149.       ENDIF
  150.  
  151. C     4. Moyenne harmonique
  152.       IF (ICAS.EQ.4) THEN
  153.         DO I=1,NX
  154.           XI=MLREE1.PROG(I)
  155.           AGR=AGR+(UN/XI)
  156.         ENDDO
  157.         XTMP=AGR
  158.         AGR=NX/XTMP
  159.         IF (DERI) THEN
  160.           DO I=1,NX
  161.             XI=MLREE1.PROG(I)
  162.             MLREE2.PROG(I)=AGR/(XI*XI*XTMP)
  163.           ENDDO
  164.         ENDIF
  165.       ENDIF
  166.  
  167. C     5. Moyenne geometrique
  168.       IF (ICAS.EQ.5) THEN
  169.         AGR=XPRO**(UN/NX)
  170.         IF (DERI) THEN
  171.           DO I=1,NX
  172.             MLREE2.PROG(I)=MLREE2.PROG(I)*(XPRO**((UN/NX)-1))/NX
  173.           ENDDO
  174.         ENDIF
  175.       ENDIF
  176.  
  177. C     6. Variance (aussi utilise pour l'ecart type, l'asymetrie, le kurtosis)
  178.       IF ((ICAS.EQ.6).OR.(ICAS.EQ.7).OR.(ICAS.EQ.8).OR.
  179.      &    (ICAS.EQ.9)) THEN
  180.         VAR=XZERO
  181.         DO I=1,NX
  182.           XI=MLREE1.PROG(I)
  183.           XM=XI-XMOY
  184.           VAR=VAR+(XM*XM)
  185.           IF (DERI) THEN
  186.             MLREE2.PROG(I)=2.D0*XM/NX
  187.           ENDIF
  188.         ENDDO
  189.         VAR=VAR/NX
  190.         IF (ICAS.EQ.6) AGR=VAR
  191.       ENDIF
  192.  
  193. C     7. Ecart type (aussi utilise pour l'asymetrie, le kurtosis)
  194.       IF ((ICAS.EQ.7).OR.(ICAS.EQ.8).OR.(ICAS.EQ.9)) THEN
  195.         SIG=SQRT(VAR)
  196.         IF (ICAS.EQ.7) THEN
  197.           AGR=SIG
  198.           IF (DERI) THEN
  199.             DO i=1,NX
  200.               MLREE2.PROG(I)=MLREE2.PROG(I)/(2.D0*SIG)
  201.             ENDDO
  202.           ENDIF
  203.         ENDIF
  204.       ENDIF
  205.  
  206. C     8. Coefficient d'asymetrie (utile aussi pour la derivee du kurtosis)
  207.       IF ((ICAS.EQ.8).OR.((ICAS.EQ.9).AND.(DERI))) THEN
  208.         ASYM=XZERO
  209.         DO I=1,NX
  210.           XI=MLREE1.PROG(I)
  211.           ASYM=ASYM+(((XI-XMOY)/SIG)**3)
  212.         ENDDO
  213.         ASYM=ASYM/NX
  214.         AGR=ASYM
  215.         IF (DERI) THEN
  216.           DO I=1,NX
  217.             XI=MLREE1.PROG(I)
  218.             ZI=(XI-XMOY)/SIG
  219.             MLREE2.PROG(I)=3.D0*(ZI**2-(ASYM*ZI)-UN)/(NX*SIG)
  220.           ENDDO
  221.         ENDIF
  222.       ENDIF
  223.  
  224. C     Kurtosis
  225.       IF (ICAS.EQ.9) THEN
  226.         XKUR=XZERO
  227.         DO I=1,NX
  228.           XI=MLREE1.PROG(I)
  229.           XKUR=XKUR+(((XI-XMOY)/SIG)**4)
  230.         ENDDO
  231.         XKUR=XKUR/NX
  232.         AGR=XKUR
  233.         IF (DERI) THEN
  234.           DO I=1,NX
  235.             XI=MLREE1.PROG(I)
  236.             ZI=(XI-XMOY)/SIG
  237.             MLREE2.PROG(I)=4.D0*(ZI**3-(XKUR*ZI)-ASYM)/(NX*SIG)
  238.           ENDDO
  239.         ENDIF
  240.       ENDIF
  241.  
  242. C     Mediane (par de derivee pour cette fonction)
  243.       IF (ICAS.EQ.10) THEN
  244. C       Tri des valeurs en ordre croissant (par insertion)
  245.         SEGINI,MLREE3=MLREE1
  246.         CALL ORDON1(MLREE3,.TRUE.,.FALSE.,0)
  247. C       Obtention de la valeur mediane
  248.         IF (MOD(NX,2).EQ.0) THEN
  249.           XMED=0.5D0*(MLREE3.PROG(NX/2)+MLREE3.PROG(NX/2+1))
  250.         ELSE
  251.           XMED=MLREE3.PROG(NX/2+1)
  252.         ENDIF
  253.         AGR=XMED
  254.         SEGSUP MLREE3
  255.       ENDIF
  256.  
  257. C     Moment d'ordre P
  258.       IF (ICAS.EQ.11) THEN
  259.         XMOMP=XZERO
  260.         DO I=1,NX
  261.           XI=MLREE1.PROG(I)
  262.           XMOMP=XMOMP+(XI**XP)
  263.           IF (DERI) THEN
  264.             MLREE2.PROG(I)=XP*(XI**(XP-UN))
  265.           ENDIF
  266.         ENDDO
  267.         AGR=XMOMP
  268.       ENDIF
  269.  
  270. C     P moyenne (aussi utlise pour LEHM)
  271.       IF ((ICAS.EQ.12).OR.(ICAS.EQ.14)) THEN
  272.         SUMP=XZERO
  273.         DO I=1,NX
  274.           XI=MLREE1.PROG(I)
  275.           IF (ROBU) XI=XI/VINF
  276.           SUMP=SUMP+(XI**XP)
  277.         ENDDO
  278.         IF (ICAS.EQ.12) THEN
  279.           XMOYP=(SUMP/NX)**(UN/XP)
  280.           IF (ROBU) XMOYP=VINF*XMOYP
  281.           AGR=XMOYP
  282.           IF (DERI) THEN
  283.             DO I=1,NX
  284.               XI=MLREE1.PROG(I)
  285.               MLREE2.PROG(I)=((XI/XMOYP)**(XP-UN))/NX
  286.             ENDDO
  287.           ENDIF
  288.         ENDIF
  289.       ENDIF
  290.  
  291. C     P norme
  292.       IF (ICAS.EQ.13) THEN
  293.         XNORP=XZERO
  294.         DO I=1,NX
  295.           XI=MLREE1.PROG(I)
  296.           IF (ROBU) XI=XI/VINF
  297.           XNORP=XNORP+((ABS(XI))**XP)
  298.         ENDDO
  299.         XNORP=XNORP**(UN/XP)
  300.         IF (ROBU) XNORP=VINF*XNORP
  301.         AGR=XNORP
  302.         IF (DERI) THEN
  303.           DO I=1,NX
  304.             XI=MLREE1.PROG(I)
  305.             MLREE2.PROG(I)=((ABS(XI)/XNORP)**(XP-UN))*SIGN(UN,XI)
  306.           ENDDO
  307.         ENDIF
  308.       ENDIF
  309.  
  310. C     Fonction de Lehmer
  311.       IF (ICAS.EQ.14) THEN
  312.         XPM1=XP-UN
  313.         SUMPM1=XZERO
  314.         DO I=1,NX
  315.           XI=MLREE1.PROG(I)
  316.           IF (ROBU) XI=XI/VINF
  317.           SUMPM1=SUMPM1+(XI**XPM1)
  318.         ENDDO
  319.         XLEHM=SUMP/SUMPM1
  320.         IF (ROBU) XLEHM=VINF*XLEHM
  321.         AGR=XLEHM
  322.         IF (DERI) THEN
  323.           IF (ROBU) SUMP=(VINF**XP)*SUMP
  324.           IF (ROBU) SUMPM1=(VINF**XPM1)*SUMPM1
  325.           XPM2=XPM1-UN
  326.           DO I=1,NX
  327.             XI=MLREE1.PROG(I)
  328.             MLREE2.PROG(I)=(XI**XPM2)*((XP*XI*SUMPM1)-(XPM1*SUMP))/
  329.      &                     (SUMPM1**2)
  330.           ENDDO
  331.         ENDIF
  332.       ENDIF
  333.  
  334. C     Fonctions de Kreisselmeir Steinhauser (aussi utilise pour Boltzmann)
  335.       IF ((ICAS.EQ.15).OR.(ICAS.EQ.16).OR.(ICAS.EQ.17)) THEN
  336.         SUMEP=XZERO
  337.         DO I=1,NX
  338.           XI=MLREE1.PROG(I)
  339.           IF (ROBU) XI=XI-VMAX
  340.           SUMEP=SUMEP+(EXP(XP*XI))
  341.         ENDDO
  342.       ENDIF
  343.       IF ((ICAS.EQ.15).OR.(ICAS.EQ.16)) THEN
  344.         IF (ICAS.EQ.15) XKS=(LOG(SUMEP/NX))/XP
  345.         IF (ICAS.EQ.16) XKS=(LOG(SUMEP))/XP
  346.         IF (ROBU) XKS=VMAX+XKS
  347.         AGR=XKS
  348.         IF (DERI) THEN
  349.           IF (ROBU) SUMEP=SUMEP*EXP(XP*VMAX)
  350.           DO I=1,NX
  351.             XI=MLREE1.PROG(I)
  352.             MLREE2.PROG(I)=EXP(XP*XI)/SUMEP
  353.           ENDDO
  354.         ENDIF
  355.       ENDIF
  356.  
  357. C     Moyenne de Boltzmann
  358.       IF (ICAS.EQ.17) THEN
  359.         SUMXEP=XZERO
  360.         DO I=1,NX
  361.           XI=MLREE1.PROG(I)
  362.           IF (ROBU) XI=XI-VMAX
  363.           SUMXEP=SUMXEP+(XI*EXP(XP*XI))
  364.         ENDDO
  365.         BOL=SUMXEP/SUMEP
  366.         IF (ROBU) BOL=VMAX+BOL
  367.         AGR=BOL
  368.         IF (DERI) THEN
  369.           IF (ROBU) SUMEP=SUMEP*EXP(XP*VMAX)
  370.           DO I=1,NX
  371.             XI=MLREE1.PROG(I)
  372.             MLREE2.PROG(I)=(UN+XP*(XI-BOL))*EXP(XP*XI)/SUMEP
  373.           ENDDO
  374.         ENDIF
  375.       ENDIF
  376.  
  377.       RETURN
  378.       END
  379.  
  380.  
  381.  

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