C AGREG1    SOURCE    FD218221  25/12/24    21:15:01     12435          
      SUBROUTINE AGREG1(MLREE1,ICAS,XP,IDERI,IROBU,AGR,MLREE2)

C     Application d'une fonction d'agregation sur un
C     objet LISTREEL

C     Entrees :
C     ---------
C       MLREE1 : Pointeur vers le LISTREEL (suppose actif)
C       ICAS   : Entier, code de la fonction
C                = 1  'SOMM'  Somme
C                = 2  'PROD'  Produit
C                = 3  'MOYE'  Moyenne arithmetique (moment d'ordre 1)
C                = 4  'MOHA'  Moyenne harmonique
C                = 5  'MOGE'  Moyenne geometrique
C                = 6  'VARI'  Variance (moment centre d'ordre 2)
C                = 7  'ECTY'  Ecart type
C                = 8  'ASYM'  Coefficient d'asymetrie (moment centre reduit d'ordre 3)
C                = 9  'KURT'  Kurtosis (moment centre reduit d'ordre 4)
C                = 10 'MEDI'  Mediane
C                = 11 'PMOM'  Moment d'ordre P
C                = 12 'PMOY'  Moyenne generalise d'ordre P
C                = 13 'PNOR'  Norme generalisee d'ordre P
C                = 14 'LEHM'  Fonction de Lehmer d'ordre P
C                = 15 'KSL'   Fonction de Kreisselmeir Steinhauser inferieure d'ordre P (MellowMax)
C                = 16 'KSU'   Fonction de Kreisselmeir Steinhauser superieure d'ordre P (LogSumExp)
C                = 17 'BOLT'  Fonction de Boltzmann d'ordre P
C       XP    : Flottant, parametre pour les fonctions 'PMOY' 'PMOM' 'PNOR' 'LEHM'
C                                                      'KSL'  'KSL'  'BOLT'
C       IDERI : Entier, pour calcul des derivees partielles
C                = 0 pas de calcul des derivees
C               != 0 calcul des derivees partielles
C       IROBU : Entier, pour calcul robuste au overflow
C                = 0 pour un calcul "naif"
C               != 0 pour un calcul "robuste" en normalisant les valeurs avec
C                    la norme infinie ou bien le maximum, selon la fonction choisie

C     Sorties :
C     ---------
C       AGR    : Flottant, valeur de la fonction d'agregation
C       MLREE2 : Pointeur vers le LISTREEL, de meme dimentsion que MLREE1, contenant les valeurs des
C                derivees partielles de la fonction d'agregation par rapport a chaque terme de MLREE1


C     Typages implicites habituels
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)

C     Les includes necessaires
-INC CCREEL
-INC SMLREEL

C     Quelques objets
      PARAMETER(UN=1.D0)
      LOGICAL ROBU,DERI

C     Taille du LISTREEL
      NX=MLREE1.PROG(/1)

C     Initialisation du resultat
      AGR=XZERO
      MLREE2=0
      IF (IDERI.EQ.1) THEN
        JG=NX
        SEGINI MLREE2
      ENDIF

C     Cas trivial non traite
      IF (NX.LT.1) THEN
        CALL ERREUR(21)
        RETURN
      ENDIF

C     Calcul robuste ?
      ROBU=.FALSE.
      IF (IROBU.NE.0) THEN
        ROBU=.TRUE.
        IF (ICAS.GE.12) THEN
          CALL MAXIN3(MLREE1,IMAX,VMAX,1,0)
          VINF=ABS(VMAX)
          VPMAX=XP*VMAX
        ENDIF
      ENDIF

C     Calcul des derivees partielles ?
      DERI=.FALSE.
      IF (IDERI.NE.0) THEN
        DERI=.TRUE.
      ENDIF

C     Cas de la P moyenne avec P=0 --> on utilise MOGE (moyenne geometrique)
      IF ((ICAS.EQ.12).AND.(ABS(XP).LT.XPETIT)) THEN
        ICAS=5
      ENDIF

C     1. Somme (aussi utilise pour la moyenne, la variance, l'ecart type,
C               l'asymetrie, le kurtosis)
      IF ((ICAS.EQ.1).OR.(ICAS.EQ.3).OR.(ICAS.EQ.6).OR.(ICAS.EQ.7).OR.
     &    (ICAS.EQ.8).OR.(ICAS.EQ.9)) THEN
        SUM=XZERO
        DO I=1,NX
          XI=MLREE1.PROG(I)
          SUM=SUM+XI
        ENDDO
        IF (ICAS.EQ.1) THEN
          AGR=SUM
          IF (DERI) THEN
            MLREE2.PROG(I)=UN
          ENDIF
        ENDIF
      ENDIF

C     2. Produit (aussi utilise pour la moyenne geometrique)
      IF ((ICAS.EQ.2).OR.(ICAS.EQ.5)) THEN
        XPRO=UN
        DO I=1,NX
          XI=MLREE1.PROG(I)
          XPRO=XPRO*XI
        ENDDO
        IF (DERI) THEN
          DO I=1,NX
            XI=MLREE1.PROG(I)
            IF (ABS(XI).LT.XPETIT) THEN
              MLREE2.PROG(I)=XZERO
            ELSE
              MLREE2.PROG(I)=XPRO/XI
            ENDIF
          ENDDO
        ENDIF
        IF (ICAS.EQ.2) AGR=XPRO
      ENDIF

C     3. Moyenne (aussi utilise pour la variance, l'ecart type, l'asymetrie, le kurtosis)
      IF ((ICAS.EQ.3).OR.(ICAS.EQ.6).OR.(ICAS.EQ.7).OR.(ICAS.EQ.8).OR.
     &    (ICAS.EQ.9)) THEN
        XMOY=SUM/NX
        IF (ICAS.EQ.3) THEN
          AGR=XMOY
          IF (DERI) THEN
            DO I=1,NX
              IF (NX.EQ.0) THEN
                MLREE2.PROG(I)=XZERO
              ELSE
                MLREE2.PROG(I)=UN/NX
              ENDIF
            ENDDO
          ENDIF
        ENDIF
      ENDIF

C     4. Moyenne harmonique
      IF (ICAS.EQ.4) THEN
        DO I=1,NX
          XI=MLREE1.PROG(I)
          AGR=AGR+(UN/XI)
        ENDDO
        XTMP=AGR
        AGR=NX/XTMP
        IF (DERI) THEN
          DO I=1,NX
            XI=MLREE1.PROG(I)
            MLREE2.PROG(I)=AGR/(XI*XI*XTMP)
          ENDDO
        ENDIF
      ENDIF

C     5. Moyenne geometrique
      IF (ICAS.EQ.5) THEN
        AGR=XPRO**(UN/NX)
        IF (DERI) THEN
          DO I=1,NX
            MLREE2.PROG(I)=MLREE2.PROG(I)*(XPRO**((UN/NX)-1))/NX
          ENDDO
        ENDIF
      ENDIF

C     6. Variance (aussi utilise pour l'ecart type, l'asymetrie, le kurtosis)
      IF ((ICAS.EQ.6).OR.(ICAS.EQ.7).OR.(ICAS.EQ.8).OR.
     &    (ICAS.EQ.9)) THEN
        VAR=XZERO
        DO I=1,NX
          XI=MLREE1.PROG(I)
          XM=XI-XMOY
          VAR=VAR+(XM*XM)
          IF (DERI) THEN
            MLREE2.PROG(I)=2.D0*XM/NX
          ENDIF
        ENDDO
        VAR=VAR/NX
        IF (ICAS.EQ.6) AGR=VAR
      ENDIF

C     7. Ecart type (aussi utilise pour l'asymetrie, le kurtosis)
      IF ((ICAS.EQ.7).OR.(ICAS.EQ.8).OR.(ICAS.EQ.9)) THEN
        SIG=SQRT(VAR)
        IF (ICAS.EQ.7) THEN
          AGR=SIG
          IF (DERI) THEN
            DO i=1,NX
              MLREE2.PROG(I)=MLREE2.PROG(I)/(2.D0*SIG)
            ENDDO
          ENDIF
        ENDIF
      ENDIF

C     8. Coefficient d'asymetrie (utile aussi pour la derivee du kurtosis)
      IF ((ICAS.EQ.8).OR.((ICAS.EQ.9).AND.(DERI))) THEN
        ASYM=XZERO
        DO I=1,NX
          XI=MLREE1.PROG(I)
          ASYM=ASYM+(((XI-XMOY)/SIG)**3)
        ENDDO
        ASYM=ASYM/NX
        AGR=ASYM
        IF (DERI) THEN
          DO I=1,NX
            XI=MLREE1.PROG(I)
            ZI=(XI-XMOY)/SIG
            MLREE2.PROG(I)=3.D0*(ZI**2-(ASYM*ZI)-UN)/(NX*SIG)
          ENDDO
        ENDIF
      ENDIF

C     Kurtosis
      IF (ICAS.EQ.9) THEN
        XKUR=XZERO
        DO I=1,NX
          XI=MLREE1.PROG(I)
          XKUR=XKUR+(((XI-XMOY)/SIG)**4)
        ENDDO
        XKUR=XKUR/NX
        AGR=XKUR
        IF (DERI) THEN
          DO I=1,NX
            XI=MLREE1.PROG(I)
            ZI=(XI-XMOY)/SIG
            MLREE2.PROG(I)=4.D0*(ZI**3-(XKUR*ZI)-ASYM)/(NX*SIG)
          ENDDO
        ENDIF
      ENDIF

C     Mediane (par de derivee pour cette fonction)
      IF (ICAS.EQ.10) THEN
C       Tri des valeurs en ordre croissant (par insertion)
        SEGINI,MLREE3=MLREE1
        CALL ORDON1(MLREE3,.TRUE.,.FALSE.,0)
C       Obtention de la valeur mediane
        IF (MOD(NX,2).EQ.0) THEN
          XMED=0.5D0*(MLREE3.PROG(NX/2)+MLREE3.PROG(NX/2+1))
        ELSE
          XMED=MLREE3.PROG(NX/2+1)
        ENDIF
        AGR=XMED
        SEGSUP MLREE3
      ENDIF

C     Moment d'ordre P
      IF (ICAS.EQ.11) THEN
        XMOMP=XZERO
        DO I=1,NX
          XI=MLREE1.PROG(I)
          XMOMP=XMOMP+(XI**XP)
          IF (DERI) THEN
            MLREE2.PROG(I)=XP*(XI**(XP-UN))
          ENDIF
        ENDDO
        AGR=XMOMP
      ENDIF

C     P moyenne (aussi utlise pour LEHM)
      IF ((ICAS.EQ.12).OR.(ICAS.EQ.14)) THEN
        SUMP=XZERO
        DO I=1,NX
          XI=MLREE1.PROG(I)
          IF (ROBU) XI=XI/VINF
          SUMP=SUMP+(XI**XP)
        ENDDO
        IF (ICAS.EQ.12) THEN
          XMOYP=(SUMP/NX)**(UN/XP)
          IF (ROBU) XMOYP=VINF*XMOYP
          AGR=XMOYP
          IF (DERI) THEN
            DO I=1,NX
              XI=MLREE1.PROG(I)
              MLREE2.PROG(I)=((XI/XMOYP)**(XP-UN))/NX
            ENDDO
          ENDIF
        ENDIF
      ENDIF

C     P norme
      IF (ICAS.EQ.13) THEN
        XNORP=XZERO
        DO I=1,NX
          XI=MLREE1.PROG(I)
          IF (ROBU) XI=XI/VINF
          XNORP=XNORP+((ABS(XI))**XP)
        ENDDO
        XNORP=XNORP**(UN/XP)
        IF (ROBU) XNORP=VINF*XNORP
        AGR=XNORP
        IF (DERI) THEN
          DO I=1,NX
            XI=MLREE1.PROG(I)
            MLREE2.PROG(I)=((ABS(XI)/XNORP)**(XP-UN))*SIGN(UN,XI)
          ENDDO
        ENDIF
      ENDIF

C     Fonction de Lehmer
      IF (ICAS.EQ.14) THEN
        XPM1=XP-UN
        SUMPM1=XZERO
        DO I=1,NX
          XI=MLREE1.PROG(I)
          IF (ROBU) XI=XI/VINF
          SUMPM1=SUMPM1+(XI**XPM1)
        ENDDO
        XLEHM=SUMP/SUMPM1
        IF (ROBU) XLEHM=VINF*XLEHM
        AGR=XLEHM
        IF (DERI) THEN
          IF (ROBU) SUMP=(VINF**XP)*SUMP
          IF (ROBU) SUMPM1=(VINF**XPM1)*SUMPM1
          XPM2=XPM1-UN
          DO I=1,NX
            XI=MLREE1.PROG(I)
            MLREE2.PROG(I)=(XI**XPM2)*((XP*XI*SUMPM1)-(XPM1*SUMP))/
     &                     (SUMPM1**2)
          ENDDO
        ENDIF
      ENDIF

C     Fonctions de Kreisselmeir Steinhauser (aussi utilise pour Boltzmann)
      IF ((ICAS.EQ.15).OR.(ICAS.EQ.16).OR.(ICAS.EQ.17)) THEN
        SUMEP=XZERO
        DO I=1,NX
          XI=MLREE1.PROG(I)
          IF (ROBU) XI=XI-VMAX
          SUMEP=SUMEP+(EXP(XP*XI))
        ENDDO
      ENDIF
      IF ((ICAS.EQ.15).OR.(ICAS.EQ.16)) THEN
        IF (ICAS.EQ.15) XKS=(LOG(SUMEP/NX))/XP
        IF (ICAS.EQ.16) XKS=(LOG(SUMEP))/XP
        IF (ROBU) XKS=VMAX+XKS
        AGR=XKS
        IF (DERI) THEN
          IF (ROBU) SUMEP=SUMEP*EXP(XP*VMAX)
          DO I=1,NX
            XI=MLREE1.PROG(I)
            MLREE2.PROG(I)=EXP(XP*XI)/SUMEP
          ENDDO
        ENDIF
      ENDIF

C     Moyenne de Boltzmann
      IF (ICAS.EQ.17) THEN
        SUMXEP=XZERO
        DO I=1,NX
          XI=MLREE1.PROG(I)
          IF (ROBU) XI=XI-VMAX
          SUMXEP=SUMXEP+(XI*EXP(XP*XI))
        ENDDO
        BOL=SUMXEP/SUMEP
        IF (ROBU) BOL=VMAX+BOL
        AGR=BOL
        IF (DERI) THEN
          IF (ROBU) SUMEP=SUMEP*EXP(XP*VMAX)
          DO I=1,NX
            XI=MLREE1.PROG(I)
            MLREE2.PROG(I)=(UN+XP*(XI-BOL))*EXP(XP*XI)/SUMEP
          ENDDO
        ENDIF
      ENDIF

      RETURN
      END