!*************************************************************** !* Calculate the median value of an array without sorting * !* ----------------------------------------------------------- * !* REFERENCE: * !* "NUMERICAL RECIPES By W.H. Press, B.P. Flannery, * !* S.A. Teukolsky and W.T. Vetterling, Cambridge * !* University Press, 1986" [BIBLI 08]. * !* ----------------------------------------------------------- * !* SAMPLE RUN: * !* * !* Given table: * !* * !* 407.8 192.8 851.1 604.4 932.3 799.4 914.5 965.8 453.7 295.1 * !* 154.5 977.4 410.2 916.2 934.7 504.8 823.2 225.2 456.6 49.0 * !* 933.5 663.0 335.3 346.6 568.7 956.1 654.7 300.7 379.6 591.9 * !* 992.9 689.6 644.7 305.4 148.2 257.2 664.6 612.1 713.0 99.7 * !* 46.5 167.6 984.6 847.2 55.4 82.7 999.0 10.7 877.7 929.4 * !* 398.1 972.8 874.1 755.1 472.1 122.8 671.4 35.5 128.8 76.8 * !* 454.2 959.2 510.1 791.3 122.8 176.6 237.9 995.8 548.3 309.8 * !* 162.6 996.5 750.0 250.6 577.7 761.1 101.9 797.1 539.0 723.5 * !* * !* Median value: 558.5251 * !* * !* F90 Release By J-P Moreau, Paris.* !* (www.jpmoreau.fr) * !*************************************************************** PROGRAM TMDIAN1 real A(100) !Table to be sorted real MAX_VALUE !Maximum value of table integer Count(1) !Current value of system clock !initialize random number generator call System_Clock(Count(1)) call Random_Seed(Put=Count) N=80 !initialize size of table MAX_VALUE = 1000 !generate random table of numbers (from 0 to 1000) do i=1, N call Random_Number(x) !returns random x >= 0 and <1 A(i)=MAX_VALUE*x end do print *,'Given table:' call TWRIT(N,A) !call MDIAN1 subroutine call MDIAN1(A,N,AMED) print *,' ' print *,'Median value: ', AMED print *,' ' stop END !******************************************************* !* Given an array X of N numbers, returns their median * !* value XMED. The array X is not modified, and is * !* accessed sequentially in each consecutive pass. * !******************************************************* SUBROUTINE MDIAN1(X,N,XMED) real X(N) real, parameter :: BIG = 1.e30, AFAC=1.5, AMP=1.5 ! Here, AMP is an overconvergence factor: on each iteration, ! we move the guess by this factor. AFAC is a factor used to ! optimize the size of the "smoothing constant" EPS at each ! iteration. A=0.5*(X(1)+X(N)) EPS=ABS(X(N)-X(1)) AP=BIG AM=-BIG 1 SUM=0. SUMX=0. NP=0 NM=0 XP=BIG XM=-BIG do J=1, N XX=X(J) if(XX.ne.A)then if(XX.gt.A)then NP=NP+1 if(XX.lt.XP) XP=XX else if(XX.lt.A)then NM=NM+1 if(XX.gt.XM) XM=XX endif DUM=1./(EPS+ABS(XX-A)) SUM=SUM+DUM SUMX=SUMX+XX*DUM endif end do print *,' NP=',NP,' NM=',NM if(NP-NM.ge.2)then !guess is too low, make another pass AM=A AA=XP+MAX(0.,SUMX/SUM-A)*AMP if(AA.gt.AP) AA=0.5*(A+AP) EPS=AFAC*ABS(AA-A) A=AA goto 1 else if(NM-NP.ge.2)then !guess is too high AM=A AA=XM+MIN(0.,SUMX/SUM-A)*AMP if(AA.lt.AM) AA=0.5*(A+AM) EPS=AFAC*ABS(AA-A) A=AA goto 1 else !guess is ok if(MOD(N,2).eq.0)then !for even N median is an average if(NP.eq.NM)then XMED=0.5*(XP+XM) else if(NP.gt.NM)then XMED=0.5*(A+XP) else XMED=0.5*(XM+A) endif else if(NP.eq.NM)then XMED=A else if(NP.gt.NM)then XMED=XP else XMED=XM endif endif endif return END !write table of size N to standard output SUBROUTINE TWRIT(N,ARR) real ARR(N) print *,' ' WRITE(*,10) (ARR(I),I=1,N) return 10 FORMAT(10F6.1) END !end of file tmdian1.f90