Per stampare: Clicca qui oppure seleziona File » Stampa nel menù del tuo browser.
-----------------------------------------------------------------------------------------------
Questo intervento è stato stampato da Guide di Dada.Net
raggiungibile a http://guide.dada.net
-----------------------------------------------------------------------------------------------
By Fortran di Giuseppe Ciaburro
URL: http://guide.dada.net/fortran/interventi/2007/05/293774.shtml
Fortran di Giuseppe Ciaburro guida dal 04-06-2002
Calcolo della mediana
PROGRAM tstmed
INTEGER, DIMENSION (1025) :: jindt
REAL, DIMENSION (1025) :: xvalt
INTERFACE
SUBROUTINE TRIIND (xvalt, IRNGT)
! triind = tri par interclassement suivant xvalt croissant
REAL, DIMENSION (:) :: xvalt
INTEGER, DIMENSION (:) :: IRNGT
END SUBROUTINE TRIIND
FUNCTION FRCTIL (xvalt, IRNGT, XORD)
! frctil = fractile d'ordre XORD de l'ensemble XVALT dont
! les rangs sont donnes par IRNGT
REAL, DIMENSION (:), INTENT (IN) :: xvalt
INTEGER, DIMENSION (:), INTENT (IN) :: IRNGT
REAL, INTENT (IN) :: XORD
REAL :: FRCTIL
END FUNCTION FRCTIL
FUNCTION XMEDIA (XDONT)
! Retourne la valeur mediane de XDONT, dont l'ordre
! est bouleverse par la procedure
REAL, DIMENSION (:), INTENT (INOUT) :: XDONT
REAL :: XMEDIA
END FUNCTION XMEDIA
END INTERFACE
CALL random_seed
DO i1 = 1, 10
CALL random_number (xvalt(:))
OPEN (1, FILE="tmp.tmp", FORM="unformatted")
WRITE (1) xvalt (:)
CLOSE (1)
CALL TRIIND (xvalt(:), jindt(:))
x1 = FRCTIL (xvalt(:), jindt(:), 0.5)
IF (Mod(i1, 2) == 0) THEN
xvalt (:) = xvalt (jindt(:))
ELSE
xvalt (:) = xvalt (jindt(1025:1:-1))
END IF
x2 = XMEDIA (xvalt(:))
IF (x1 /= x2) write (*,*) x1, x2, i1
END DO
END PROGRAM tstmed
SUBROUTINE TRIIND (XVALT, IRNGT)
! triind = tri par interclassement suivant xvalt croissant
REAL, DIMENSION (:) :: XVALT
INTEGER, DIMENSION (:) :: IRNGT
! __________________________________________________________
INTEGER, DIMENSION (:), ALLOCATABLE :: JWRKT
!
NVAL = MIN (SIZE (XVALT), SIZE (IRNGT))
IF (NVAL <= 0) THEN
RETURN
ENDIF
!
! .... initialisation du tableau d'indices
! (on cree des monotonies de longueurs 2)
!
DO IIND = 2, NVAL, 2
IF (XVALT (IIND - 1) < XVALT (IIND)) THEN
IRNGT (IIND - 1) = IIND - 1
IRNGT (IIND) = IIND
ELSE
IRNGT (IIND - 1) = IIND
IRNGT (IIND) = IIND - 1
ENDIF
ENDDO
IF (MOD (NVAL, 2) /= 0) THEN
IRNGT (NVAL) = NVAL
ENDIF
!
! .... initialisation des monotonies 'a' et 'c'
!
ALLOCATE (JWRKT (1:NVAL))
LMTNC = 2
LMTNA = 2
!
! .... initialisation nouvelles monotonies
!
DO
IF (LMTNA >= NVAL) EXIT
IWRKF = 0
LMTNC = 2 * LMTNC
IWRK = 0
!
! .... passage aux monotonies suivantes 'a', 'b' et 'c'
! on regarde si on deborde, auquel cas on regarde
! si on a encore 2 monotonies a interclasser
!
DO
IINDA = IWRKF
IWRKD = IWRKF + 1
IWRKF = IINDA + LMTNC
JINDA = IINDA + LMTNA
IF (IWRKF >= NVAL) THEN
IF (JINDA >= NVAL) EXIT
IWRKF = NVAL
ENDIF
IINDB = JINDA
!
! .... on compare la derniere valeur de 'a' et la 1ere de 'b',
! pour sauter directement si a < b
!
IVALA = IRNGT (JINDA)
IVALB = IRNGT (JINDA + 1)
IF (XVALT (IVALA) <= XVALT (IVALB)) THEN
IWRK = IWRKF
CYCLE
ENDIF
!
! .... boucle sur la nouvelle monotonie 'c'
! qu'on cree dans le tableau wrk
!
DO
IF (IWRK >= IWRKF) THEN
!
! .... on recopie le tableau de travail dans le tableau d'indices
!
IRNGT (IWRKD:IWRKF) = JWRKT (IWRKD:IWRKF)
EXIT
ENDIF
!
IWRK = IWRK + 1
!
! .... monotonie 'a' et 'b' ne sont pas epuisees
!
IF (IINDA < JINDA) THEN
IF (IINDB < IWRKF) THEN
IVALA = IRNGT (IINDA + 1)
IVALB = IRNGT (IINDB + 1)
IF (XVALT (IVALA) > XVALT (IVALB)) THEN
IINDB = IINDB + 1
JWRKT (IWRK) = IRNGT (IINDB)
ELSE
IINDA = IINDA + 1
JWRKT (IWRK) = IRNGT (IINDA)
ENDIF
ELSE
!
! .... monotonie 'b' epuisee
!
IINDA = IINDA + 1
JWRKT (IWRK) = IRNGT (IINDA)
ENDIF
ELSE
!
! .... monotonie 'a' epuisee
!
IRNGT (IWRKD:IINDB) = JWRKT (IWRKD:IINDB)
IWRK = IWRKF
EXIT
ENDIF
!
ENDDO
ENDDO
!
! ... on double les monotonies de depart et on recommence
!
LMTNA = 2 * LMTNA
ENDDO
!
! ... fin
!
DEALLOCATE (JWRKT)
RETURN
END SUBROUTINE TRIIND
FUNCTION FRCTIL (XVALT, IRNGT, XORD)
! frctil = fractile d'ordre XORD de l'ensemble XVALT dont
! les rangs sont donnes par IRNGT
REAL, DIMENSION (:), INTENT (IN) :: XVALT
INTEGER, DIMENSION (:), INTENT (IN) :: IRNGT
REAL, INTENT (IN) :: XORD
REAL :: FRCTIL
! __________________________________________________________
!
! ... nombre de valeurs rangees au total
!
NVAL = SIZE (IRNGT)
IF (NVAL <= 0 .OR. NVAL > SIZE (XVALT)) THEN
WRITE (*, *) 'Dimensions incorrectes'
FRCTIL = 0.0
RETURN
ENDIF
IF (XORD <= 0.0 .OR. XORD >= 1.0) THEN
WRITE (*, *) 'Ordre incorrect ( ]0.0, 1.0[ )'
FRCTIL = 0.0
RETURN
ENDIF
!
! ... encadrement de la valeur
!
IINF = 1 + FLOOR (XORD * REAL (NVAL - 1))
ISUP = 1 + CEILING (XORD * REAL (NVAL - 1))
ITST = MIN ((IINF - 1), (NVAL - ISUP))
IF (ITST == 0) THEN
WRITE (*, *) 'Significativite douteuse'
ENDIF
!
XVAL1 = XVALT (IRNGT (IINF))
XVAL2 = XVALT (IRNGT (ISUP))
IF (XVAL1 == XVAL2) THEN
FRCTIL = XVAL1
ELSE
!
! ... interpolation
!
XORD1 = REAL (IINF - 1) / REAL (NVAL - 1)
XORD2 = REAL (ISUP - 1) / REAL (NVAL - 1)
FRCTIL = ((XORD - XORD1) * XVAL2 + &
(XORD2 - XORD) * XVAL1) / &
(XORD2 - XORD1)
ENDIF
RETURN
END FUNCTION FRCTIL
FUNCTION XMEDIA (XDONT)
! Retourne la valeur mediane de XDONT, dont l'ordre
! est bouleverse par la procedure
REAL, DIMENSION (:), INTENT (INOUT) :: XDONT
REAL :: XMEDIA
! __________________________________________________________
NDON = SIZE (XDONT)
NMED = (NDON + 1) / 2
IDEB = 1
IFIN = NDON
DO
IF (IDEB >= IFIN) EXIT
IMIL = (IDEB + IFIN) / 2
!
! On choisit un pivot, mediane du premier, dernier et milieu
! de l'ensemble auquel on a restreint la recherche
!
IF (XDONT (IMIL) < XDONT (IDEB)) THEN
XWRK = XDONT (IDEB)
XDONT (IDEB) = XDONT (IMIL)
XDONT (IMIL) = XWRK
ENDIF
IF (XDONT (IMIL) > XDONT (IFIN)) THEN
XWRK = XDONT (IFIN)
XDONT (IFIN) = XDONT (IMIL)
XDONT (IMIL) = XWRK
IF (XDONT (IMIL) < XDONT (IDEB)) THEN
XWRK = XDONT (IDEB)
XDONT (IDEB) = XDONT (IMIL)
XDONT (IMIL) = XWRK
ENDIF
ENDIF
IF ((IFIN - IDEB) < 3) EXIT
XPIV = XDONT (IMIL)
!
! On met, par echanges, les valeurs > pivot a la fin
! et les valeurs <= au debut
!
ICRS = IDEB
IDCR = IFIN
ECHG: DO
DO
ICRS = ICRS + 1
IF (ICRS >= IDCR) THEN
!
! le premier > pivot est IDCR
! le dernier <= pivot est ICRS - 1
! NB: si on arrive la au premier tour, le pivot est le max
! de l'ensemble, la derniere valeur lui est egale, et
! on peut aussi reduire de 1 l'ensemble a explorer
! comme si XDONT (IFIN) > XPIV
!
EXIT ECHG
ENDIF
IF (XDONT (ICRS) > XPIV) EXIT
ENDDO
DO
IF (XDONT (IDCR) <= XPIV) EXIT
IDCR = IDCR - 1
IF (ICRS >= IDCR) THEN
!
! le dernier <= pivot est toujours ICRS - 1
!
EXIT ECHG
ENDIF
ENDDO
!
XWRK = XDONT (IDCR)
XDONT (IDCR) = XDONT (ICRS)
XDONT (ICRS) = XWRK
ENDDO ECHG
!
! On restreint la recherche a la partie ou se trouve
! la mediane
!
IF (ICRS <= NMED) IDEB = ICRS
IF (ICRS > NMED) IFIN = ICRS - 1
ENDDO
!
XMEDIA = XDONT (NMED)
RETURN
END FUNCTION XMEDIA