-
Notifications
You must be signed in to change notification settings - Fork 6
/
mo_orderpack.f90
16157 lines (15986 loc) · 501 KB
/
mo_orderpack.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
!> \file mo_orderpack.f90
!> \brief Sort and ranking routines
!> \details This module is the Orderpack 2.0 from Michel Olagnon.
!> It provides order and unconditional, unique, and partial
!> ranking, sorting, and permutation.
!> \authors Michel Olagnon
!> \date 2000-2012
MODULE mo_orderpack
! Written, Matthias Cuntz, Apr 2014 - adapted to JAMS library
! - one module, cleaned all warnings
! Modified, Juliane Mai, Nov 2014 - replaced floating comparison by ne(), eq(), etc. from mo_utils
! Modified, Matthias Cuntz, Jul 2015 - median -> omedian
! License
! -------
! This file is part of the JAMS Fortran package, distributed under the MIT License.
!
! Copyright (c) 2014-2015 Michel Olagnon, Matthias Cuntz - mc (at) macu (dot) de
!
! Permission is hereby granted, free of charge, to any person obtaining a copy
! of this software and associated documentation files (the "Software"), to deal
! in the Software without restriction, including without limitation the rights
! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is
! furnished to do so, subject to the following conditions:
!
! The above copyright notice and this permission notice shall be included in all
! copies or substantial portions of the Software.
!
! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
! SOFTWARE.
USE mo_kind, ONLY: i4, sp, dp
USE mo_utils, ONLY: ne, eq, le
IMPLICIT NONE
! Public routines
public :: sort ! alias for refsor
public :: sort_index ! wrapper for mrgrnk
public :: ctrper
public :: fndnth
public :: indmed
public :: indnth
public :: inspar
public :: inssor
public :: omedian
public :: mrgref
public :: mrgrnk
public :: mulcnt
public :: rapknr
public :: refpar
#ifndef __PYTHON__
public :: refsor
#endif
public :: rinpar
public :: rnkpar
public :: uniinv
public :: unipar
public :: unirnk
public :: unista
public :: valmed
public :: valnth
!> Unconditional ranking
!>
!> Subroutine MRGRNK (XVALT, IMULT)
!> Ranks array XVALT into index array IRNGT, using merge-sort
!> For performance reasons, the first 2 passes are taken out of the
!> standard loop, and use dedicated coding.
!>
!> Subroutine MRGREF (XVALT, IRNGT) Ranks array XVALT into index array
!> IRNGT, using merge-sort This version is not optimized for performance,
!> and is thus not as difficult to read as the previous one.
!>
!> Partial ranking
!>
!> Subroutine RNKPAR (XVALT, IRNGT, NORD) Ranks partially XVALT by IRNGT,
!> up to order NORD (refined for speed) This routine uses a pivoting
!> strategy such as the one of finding the median based on the quicksort
!> algorithm, but we skew the pivot choice to try to bring it to NORD as
!> fast as possible. It uses 2 temporary arrays, one where it stores the
!> indices of the values smaller than the pivot, and the other for the
!> indices of values larger than the pivot that we might still need later
!> on. It iterates until it can bring the number of values in ILOWT to
!> exactly NORD, and then uses an insertion sort to rank this set, since
!> it is supposedly small.
!>
!> Subroutine RAPKNR (XVALT, IRNGT, NORD) Same as RNKPAR, but in
!> decreasing order (RAPKNR = RNKPAR spelt backwards).
!>
!> Subroutine REFPAR (XVALT, IRNGT, NORD) Ranks partially XVALT by IRNGT,
!> up to order NORD This version is not optimized for performance, and is
!> thus not as difficult to read as some other ones. It uses a pivoting
!> strategy such as the one of finding the median based on the quicksort
!> algorithm. It uses a temporary array, where it stores the partially
!> ranked indices of the values. It iterates until it can bring the
!> number of values lower than the pivot to exactly NORD, and then uses
!> an insertion sort to rank this set, since it is supposedly small.
!>
!> Subroutine RINPAR (XVALT, IRNGT, NORD) Ranks partially XVALT by IRNGT,
!> up to order NORD This version is not optimized for performance, and is
!> thus not as difficult to read as some other ones. It uses insertion
!> sort, limiting insertion to the first NORD values. It does not use any
!> work array and is faster when NORD is very small (2-5), but worst case
!> behavior (intially inverse sorted) can easily happen. In many cases,
!> the refined quicksort method is faster.
!>
!> Integer Function INDNTH (XVALT, NORD) Returns the index of the NORDth
!> value of XVALT (in increasing order) This routine uses a pivoting
!> strategy such as the one of finding the median based on the quicksort
!> algorithm, but we skew the pivot choice to try to bring it to NORD as
!> fast as possible. It uses 2 temporary arrays, one where it stores the
!> indices of the values smaller than the pivot, and the other for the
!> indices of values larger than the pivot that we might still need later
!> on. It iterates until it can bring the number of values in ILOWT to
!> exactly NORD, and then takes out the original index of the maximum
!> value in this set.
!>
!> Subroutine INDMED (XVALT, INDM) Returns the index of the median
!> (((Size(XVALT)+1))/2th value) of XVALT This routine uses the recursive
!> procedure described in Knuth, The Art of Computer Programming, vol. 3,
!> 5.3.3 - This procedure is linear in time, and does not require to be
!> able to interpolate in the set as the one used in INDNTH. It also has
!> better worst case behavior than INDNTH, but is about 10% slower in
!> average for random uniformly distributed values.
!>
!> Note that in Orderpack 1.0, this routine was a Function procedure, and
!> is now changed to a Subroutine.
!>
!> Unique ranking
!>
!> Subroutine UNIRNK (XVALT, IRNGT, NUNI) Ranks an array, removing
!> duplicate entries (uses merge sort). The routine is similar to pure
!> merge-sort ranking, but on the last pass, it discards indices that
!> correspond to duplicate entries. For performance reasons, the first 2
!> passes are taken out of the standard loop, and use dedicated coding.
!>
!> Subroutine UNIPAR (XVALT, IRNGT, NORD) Ranks partially XVALT by IRNGT,
!> up to order NORD at most, removing duplicate entries This routine uses
!> a pivoting strategy such as the one of finding the median based on the
!> quicksort algorithm, but we skew the pivot choice to try to bring it
!> to NORD as quickly as possible. It uses 2 temporary arrays, one where
!> it stores the indices of the values smaller than the pivot, and the
!> other for the indices of values larger than the pivot that we might
!> still need later on. It iterates until it can bring the number of
!> values in ILOWT to exactly NORD, and then uses an insertion sort to
!> rank this set, since it is supposedly small. At all times, the NORD
!> first values in ILOWT correspond to distinct values of the input
!> array.
!>
!> Subroutine UNIINV (XVALT, IGOEST) Inverse ranking of an array, with
!> removal of duplicate entries The routine is similar to pure merge-sort
!> ranking, but on the last pass, it sets indices in IGOEST to the rank
!> of the original value in an ordered set with duplicates removed. For
!> performance reasons, the first 2 passes are taken out of the standard
!> loop, and use dedicated coding.
!>
!> Subroutine MULCNT (XVALT, IMULT) Gives, for each array value, its
!> multiplicity The number of times that a value appears in the array is
!> computed by using inverse ranking, counting for each rank the number
!> of values that ``collide'' to this rank, and returning this sum to the
!> locations in the original set. Uses subroutine UNIINV.
!>
!> Random permutation: an interesting use of ranking
!>
!> A variation of the following problem was raised on the internet
!> sci.math.num-analysis news group: Given an array, I would like to find
!> a random permutation of this array that I could control with a
!> ``nearbyness'' parameter so that elements stay close to their initial
!> locations. The ``nearbyness'' parameter ranges from 0 to 1, with 0
!> such that no element moves from its initial location, and 1 such that
!> the permutation is fully random.
!>
!> Subroutine CTRPER (XVALT, PCLS) Permute array XVALT randomly, but
!> leaving elements close to their initial locations The routine takes
!> the 1...size(XVALT) index array as real values, takes a combination of
!> these values and of random values as a perturbation of the index
!> array, and sorts the initial set according to the ranks of these
!> perturbated indices. The relative proportion of initial order and
!> random order is 1-PCLS / PCLS, thus when PCLS = 0, there is no change
!> in the order whereas the new order is fully random when PCLS = 1. Uses
!> subroutine MRGRNK.
!>
!> The above solution found another application when I was asked the
!> following question: I am given two arrays, representing parents'
!> incomes and their children's incomes, but I do not know which parents
!> correspond to which children. I know from an independent source the
!> value of the correlation coefficient between the incomes of the
!> parents and of their children. I would like to pair the elements of
!> these arrays so that the given correlation coefficient is attained,
!> i.e. to reconstruct a realistic dataset, though very likely not to be
!> the true one.
!>
!> Program GIVCOR Given two arrays of equal length of unordered values,
!> find a "matching value" in the second array for each value in the
!> first so that the global correlation coefficient reaches exactly a
!> given target The routine first sorts the two arrays, so as to get the
!> match of maximum possible correlation. It then iterates, applying the
!> random permutation algorithm of controlled disorder ctrper to the
!> second array. When the resulting correlation goes beyond (lower than)
!> the target correlation, one steps back and reduces the disorder
!> parameter of the permutation. When the resulting correlation lies
!> between the current one and the target, one replaces the array with
!> the newly permuted one. When the resulting correlation increases from
!> the current value, one increases the disorder parameter. That way, the
!> target correlation is approached from above, by a controlled increase
!> in randomness. Since full randomness leads to zero correlation, the
!> iterations meet the desired coefficient at some point. It may be noted
!> that there could be some cases when one would get stuck in a sort of
!> local minimum, where local perturbations cannot further reduce the
!> correlation and where global ones lead to overpass the target. It
!> seems easier to restart the program with a different seed when this
!> occurs than to design an avoidance scheme. Also, should a negative
!> correlation be desired, the program should be modified to start with
!> one array in reverse order with respect to the other, i.e. coorelation
!> as close to -1 as possible.
!>
!>
!> Sorting
!>
!> Full sorting
!>
!> Subroutine INSSOR (XVALT) Sorts XVALT into increasing order (Insertion
!> sort) This subroutine uses insertion sort. It does not use any work
!> array and is faster when XVALT is of very small size (< 20), or
!> already almost sorted, but worst case behavior (intially inverse
!> sorted) can easily happen. In most cases, the quicksort or merge sort
!> method is faster.
!>
!> Subroutine REFSOR (XVALT) Sorts XVALT into increasing order (Quick
!> sort) This version is not optimized for performance, and is thus not
!> as difficult to read as some other ones. This subroutine uses
!> quicksort in a recursive implementation, and insertion sort for the
!> last steps with small subsets. It does not use any work array
!>
!> Partial sorting
!>
!> Subroutine INSPAR (XVALT, NORD) Sorts partially XVALT, bringing the
!> NORD lowest values at the begining of the array. This subroutine uses
!> insertion sort, limiting insertion to the first NORD values. It does
!> not use any work array and is faster when NORD is very small (2-5),
!> but worst case behavior can happen fairly probably (initially inverse
!> sorted). In many cases, the refined quicksort method is faster.
!>
!> Function FNDNTH (XVALT, NORD) Finds out and returns the NORDth value
!> in XVALT (ascending order) This subroutine uses insertion sort,
!> limiting insertion to the first NORD values, and even less when one
!> can know that the value that is considered will not be the NORDth. It
!> uses only a work array of size NORD and is faster when NORD is very
!> small (2-5), but worst case behavior can happen fairly probably
!> (initially inverse sorted). In many cases, the refined quicksort
!> method implemented by VALNTH / INDNTH is faster, though much more
!> difficult to read and understand.
!>
!> Function VALNTH (XVALT, NORD) Finds out and returns the NORDth value
!> in XVALT (ascending order) This subroutine simply calls INDNTH.
!>
!> Function VALMED (XVALT) Finds out and returns the median
!> (((Size(XVALT)+1))/2th value) of XVALT This routine uses the recursive
!> procedure described in Knuth, The Art of Computer Programming, vol. 3,
!> 5.3.3 - This procedure is linear in time, and does not require to be
!> able to interpolate in the set as the one used in VALNTH/INDNTH. It
!> also has better worst case behavior than VALNTH/INDNTH, and is about
!> 20% faster in average for random uniformly distributed values.
!>
!> Function OMEDIAN (XVALT) It is a modified version of VALMED that
!> provides the average between the two middle values in the case
!> Size(XVALT) is even.
!>
!> Unique sorting
!>
!> Subroutine UNISTA (XVALT, NUNI) Removes duplicates from an array This
!> subroutine uses merge sort unique inverse ranking. It leaves in the
!> initial set only those entries that are unique, packing the array, and
!> leaving the order of the retained values unchanged.
! aliases/wrapper
interface sort
module procedure d_refsor, r_refsor, i_refsor
end interface sort
interface sort_index
module procedure sort_index_dp, sort_index_sp, sort_index_i4
end interface sort_index
interface ctrper
module procedure d_ctrper, r_ctrper, i_ctrper
end interface ctrper
interface fndnth
module procedure d_fndnth, r_fndnth, i_fndnth
end interface fndnth
interface indmed
module procedure d_indmed, r_indmed, i_indmed
end interface indmed
interface indnth
module procedure d_indnth, r_indnth, i_indnth
end interface indnth
interface inspar
module procedure d_inspar, r_inspar, i_inspar
end interface inspar
interface inssor
module procedure d_inssor, r_inssor, i_inssor
end interface inssor
interface omedian
module procedure d_median, r_median, i_median
end interface omedian
interface mrgref
module procedure d_mrgref, r_mrgref, i_mrgref
end interface mrgref
interface mrgrnk
module procedure D_mrgrnk, R_mrgrnk, I_mrgrnk
end interface mrgrnk
interface mulcnt
module procedure d_mulcnt, r_mulcnt, i_mulcnt
end interface mulcnt
interface rapknr
module procedure d_rapknr, r_rapknr, i_rapknr
end interface rapknr
interface refpar
module procedure d_refpar, r_refpar, i_refpar
end interface refpar
#ifndef __PYTHON__
interface refsor
module procedure d_refsor, r_refsor, i_refsor
end interface refsor
#endif
interface rinpar
module procedure d_rinpar, r_rinpar, i_rinpar
end interface rinpar
interface rnkpar
module procedure d_rnkpar, r_rnkpar, i_rnkpar
end interface rnkpar
interface uniinv
module procedure d_uniinv, r_uniinv, i_uniinv
end interface uniinv
interface nearless
module procedure D_nearless, R_nearless, I_nearless
end interface nearless
interface unipar
module procedure d_unipar, r_unipar, i_unipar
end interface unipar
interface unirnk
module procedure D_unirnk, R_unirnk, I_unirnk
end interface unirnk
interface unista
module procedure d_unista, r_unista, i_unista
end interface unista
interface valmed
module procedure d_valmed, r_valmed, i_valmed
end interface valmed
interface valnth
module procedure d_valnth, r_valnth, i_valnth
end interface valnth
private :: R_ctrper, I_ctrper, D_ctrper
private :: R_fndnth, I_fndnth, D_fndnth
private :: R_indmed, I_indmed, D_indmed
private :: R_indnth, I_indnth, D_indnth
private :: R_inspar, I_inspar, D_inspar
private :: R_inssor, I_inssor, D_inssor
private :: R_median, I_median, D_median
private :: R_mrgref, I_mrgref, D_mrgref
private :: R_mrgrnk, I_mrgrnk, D_mrgrnk
private :: R_mulcnt, I_mulcnt, D_mulcnt
private :: R_nearless, I_nearless, D_nearless, nearless
private :: R_rapknr, I_rapknr, D_rapknr
private :: R_refpar, I_refpar, D_refpar
private :: R_refsor, I_refsor, D_refsor
private :: R_rinpar, I_rinpar, D_rinpar
private :: R_rnkpar, I_rnkpar, D_rnkpar
private :: R_subsor, I_subsor, D_subsor
private :: R_uniinv, I_uniinv, D_uniinv
private :: R_unipar, I_unipar, D_unipar
private :: R_unirnk, I_unirnk, D_unirnk
private :: R_unista, I_unista, D_unista
private :: R_valmed, I_valmed, D_valmed
private :: R_valnth, I_valnth, D_valnth
private :: r_med, i_med, d_med
PRIVATE
Integer(kind=i4), Allocatable, Dimension(:), Save :: IDONT
CONTAINS
! ------------------------------------------------------------------
FUNCTION sort_index_dp(arr)
IMPLICIT NONE
REAL(dp), DIMENSION(:), INTENT(IN) :: arr
INTEGER(i4), DIMENSION(size(arr)) :: sort_index_dp
call mrgrnk(arr, sort_index_dp)
END FUNCTION sort_index_dp
FUNCTION sort_index_sp(arr)
IMPLICIT NONE
REAL(sp), DIMENSION(:), INTENT(IN) :: arr
INTEGER(i4), DIMENSION(size(arr)) :: sort_index_sp
call mrgrnk(arr, sort_index_sp)
END FUNCTION sort_index_sp
FUNCTION sort_index_i4(arr)
IMPLICIT NONE
integer(i4), DIMENSION(:), INTENT(IN) :: arr
INTEGER(i4), DIMENSION(size(arr)) :: sort_index_i4
call mrgrnk(arr, sort_index_i4)
END FUNCTION sort_index_i4
! ------------------------------------------------------------------
Subroutine D_ctrper (XDONT, PCLS)
! Permute array XVALT randomly, but leaving elements close
! to their initial locations (nearbyness is controled by PCLS).
! _________________________________________________________________
! The routine takes the 1...size(XVALT) index array as real
! values, takes a combination of these values and of random
! values as a perturbation of the index array, and sorts the
! initial set according to the ranks of these perturbated indices.
! The relative proportion of initial order and random order
! is 1-PCLS / PCLS, thus when PCLS = 0, there is no change in
! the order whereas the new order is fully random when PCLS = 1.
! Michel Olagnon - May 2000.
! _________________________________________________________________
! __________________________________________________________
real(kind=dp), Dimension (:), Intent (InOut) :: XDONT
Real(kind=dp), Intent (In) :: PCLS
! __________________________________________________________
!
Real(kind=dp), Dimension (Size(XDONT)) :: XINDT
Integer(kind=i4), Dimension (Size(XDONT)) :: JWRKT
Real(kind=dp) :: PWRK
Integer(kind=i4) :: I
Real(kind=dp), Dimension (Size(XDONT)) :: II
!
Call Random_number (XINDT(:))
PWRK = Min (Max (0.0_dp, PCLS), 1.0_dp)
XINDT = Real(Size(XDONT),dp) * XINDT
forall(I=1:size(XDONT)) II(I) = real(I,dp) ! for gnu compiler to be initialised
XINDT = PWRK*XINDT + (1.0_dp-PWRK)*II
Call MRGRNK (XINDT, JWRKT)
XDONT = XDONT (JWRKT)
!
End Subroutine D_ctrper
Subroutine R_ctrper (XDONT, PCLS)
! Permute array XVALT randomly, but leaving elements close
! to their initial locations (nearbyness is controled by PCLS).
! _________________________________________________________________
! The routine takes the 1...size(XVALT) index array as real
! values, takes a combination of these values and of random
! values as a perturbation of the index array, and sorts the
! initial set according to the ranks of these perturbated indices.
! The relative proportion of initial order and random order
! is 1-PCLS / PCLS, thus when PCLS = 0, there is no change in
! the order whereas the new order is fully random when PCLS = 1.
! Michel Olagnon - May 2000.
! _________________________________________________________________
! _________________________________________________________
Real(kind=sp), Dimension (:), Intent (InOut) :: XDONT
Real(kind=sp), Intent (In) :: PCLS
! __________________________________________________________
!
Real(kind=sp), Dimension (Size(XDONT)) :: XINDT
Integer(kind=i4), Dimension (Size(XDONT)) :: JWRKT
Real(kind=sp) :: PWRK
Integer(kind=i4) :: I
Real(kind=sp), Dimension (Size(XDONT)) :: II
!
Call Random_number (XINDT(:))
PWRK = Min (Max (0.0, PCLS), 1.0)
XINDT = Real(Size(XDONT),sp) * XINDT
forall(I=1:size(XDONT)) II(I) = real(I,sp) ! for gnu compiler to be initialised
XINDT = PWRK*XINDT + (1.0-PWRK)*II
Call MRGRNK (XINDT, JWRKT)
XDONT = XDONT (JWRKT)
!
End Subroutine R_ctrper
Subroutine I_ctrper (XDONT, PCLS)
! Permute array XVALT randomly, but leaving elements close
! to their initial locations (nearbyness is controled by PCLS).
! _________________________________________________________________
! The routine takes the 1...size(XVALT) index array as real
! values, takes a combination of these values and of random
! values as a perturbation of the index array, and sorts the
! initial set according to the ranks of these perturbated indices.
! The relative proportion of initial order and random order
! is 1-PCLS / PCLS, thus when PCLS = 0, there is no change in
! the order whereas the new order is fully random when PCLS = 1.
! Michel Olagnon - May 2000.
! _________________________________________________________________
! __________________________________________________________
Integer(kind=i4), Dimension (:), Intent (InOut) :: XDONT
Real(kind=sp), Intent (In) :: PCLS
! __________________________________________________________
!
Real(kind=sp), Dimension (Size(XDONT)) :: XINDT
Integer(kind=i4), Dimension (Size(XDONT)) :: JWRKT
Real(kind=sp) :: PWRK
Integer(kind=i4) :: I
Real(kind=sp), Dimension (Size(XDONT)) :: II
!
Call Random_number (XINDT(:))
PWRK = Min (Max (0.0, PCLS), 1.0)
XINDT = Real(Size(XDONT),sp) * XINDT
forall(I=1:size(XDONT)) II(I) = real(I,sp) ! for gnu compiler to be initialised
XINDT = PWRK*XINDT + (1.0-PWRK)*II
Call MRGRNK(XINDT, JWRKT)
XDONT = XDONT(JWRKT)
!
End Subroutine I_ctrper
Function D_fndnth (XDONT, NORD) Result (FNDNTH)
! Return NORDth value of XDONT, i.e fractile of order NORD/SIZE(XDONT).
! ______________________________________________________________________
! This subroutine uses insertion sort, limiting insertion
! to the first NORD values. It is faster when NORD is very small (2-5),
! and it requires only a workarray of size NORD and type of XDONT,
! but worst case behavior can happen fairly probably (initially inverse
! sorted). In many cases, the refined quicksort method is faster.
! Michel Olagnon - Aug. 2000
! __________________________________________________________
! __________________________________________________________
real(Kind=dp), Dimension (:), Intent (In) :: XDONT
real(Kind=dp) :: FNDNTH
Integer(kind=i4), Intent (In) :: NORD
! __________________________________________________________
real(Kind=dp), Dimension (NORD) :: XWRKT
real(Kind=dp) :: XWRK, XWRK1
!
!
Integer(kind=i4) :: ICRS, IDCR, ILOW, NDON
!
XWRKT (1) = XDONT (1)
Do ICRS = 2, NORD
XWRK = XDONT (ICRS)
Do IDCR = ICRS - 1, 1, - 1
If (XWRK >= XWRKT(IDCR)) Exit
XWRKT (IDCR+1) = XWRKT (IDCR)
End Do
XWRKT (IDCR+1) = XWRK
End Do
!
NDON = SIZE (XDONT)
XWRK1 = XWRKT (NORD)
ILOW = 2*NORD - NDON
Do ICRS = NORD + 1, NDON
If (XDONT(ICRS) < XWRK1) Then
XWRK = XDONT (ICRS)
Do IDCR = NORD - 1, MAX (1, ILOW) , - 1
If (XWRK >= XWRKT(IDCR)) Exit
XWRKT (IDCR+1) = XWRKT (IDCR)
End Do
XWRKT (IDCR+1) = XWRK
XWRK1 = XWRKT(NORD)
End If
ILOW = ILOW + 1
End Do
FNDNTH = XWRK1
!
End Function D_fndnth
Function R_fndnth (XDONT, NORD) Result (FNDNTH)
! Return NORDth value of XDONT, i.e fractile of order NORD/SIZE(XDONT).
! ______________________________________________________________________
! This subroutine uses insertion sort, limiting insertion
! to the first NORD values. It is faster when NORD is very small (2-5),
! and it requires only a workarray of size NORD and type of XDONT,
! but worst case behavior can happen fairly probably (initially inverse
! sorted). In many cases, the refined quicksort method is faster.
! Michel Olagnon - Aug. 2000
! __________________________________________________________
! _________________________________________________________
Real(kind=sp), Dimension (:), Intent (In) :: XDONT
Real(kind=sp) :: FNDNTH
Integer(kind=i4), Intent (In) :: NORD
! __________________________________________________________
Real(kind=sp), Dimension (NORD) :: XWRKT
Real(kind=sp) :: XWRK, XWRK1
!
!
Integer(kind=i4) :: ICRS, IDCR, ILOW, NDON
!
XWRKT (1) = XDONT (1)
Do ICRS = 2, NORD
XWRK = XDONT (ICRS)
Do IDCR = ICRS - 1, 1, - 1
If (XWRK >= XWRKT(IDCR)) Exit
XWRKT (IDCR+1) = XWRKT (IDCR)
End Do
XWRKT (IDCR+1) = XWRK
End Do
!
NDON = SIZE (XDONT)
XWRK1 = XWRKT (NORD)
ILOW = 2*NORD - NDON
Do ICRS = NORD + 1, NDON
If (XDONT(ICRS) < XWRK1) Then
XWRK = XDONT (ICRS)
Do IDCR = NORD - 1, MAX (1, ILOW) , - 1
If (XWRK >= XWRKT(IDCR)) Exit
XWRKT (IDCR+1) = XWRKT (IDCR)
End Do
XWRKT (IDCR+1) = XWRK
XWRK1 = XWRKT(NORD)
End If
ILOW = ILOW + 1
End Do
FNDNTH = XWRK1
!
End Function R_fndnth
Function I_fndnth (XDONT, NORD) Result (FNDNTH)
! Return NORDth value of XDONT, i.e fractile of order NORD/SIZE(XDONT).
! ______________________________________________________________________
! This subroutine uses insertion sort, limiting insertion
! to the first NORD values. It is faster when NORD is very small (2-5),
! and it requires only a workarray of size NORD and type of XDONT,
! but worst case behavior can happen fairly probably (initially inverse
! sorted). In many cases, the refined quicksort method is faster.
! Michel Olagnon - Aug. 2000
! __________________________________________________________
! __________________________________________________________
Integer(kind=i4), Dimension (:), Intent (In) :: XDONT
Integer(kind=i4) :: fndnth
Integer(kind=i4), Intent (In) :: NORD
! __________________________________________________________
Integer(kind=i4), Dimension (NORD) :: XWRKT
Integer(kind=i4) :: XWRK, XWRK1
!
!
Integer(kind=i4) :: ICRS, IDCR, ILOW, NDON
!
XWRKT (1) = XDONT (1)
Do ICRS = 2, NORD
XWRK = XDONT (ICRS)
Do IDCR = ICRS - 1, 1, - 1
If (XWRK >= XWRKT(IDCR)) Exit
XWRKT (IDCR+1) = XWRKT (IDCR)
End Do
XWRKT (IDCR+1) = XWRK
End Do
!
NDON = SIZE (XDONT)
XWRK1 = XWRKT (NORD)
ILOW = 2*NORD - NDON
Do ICRS = NORD + 1, NDON
If (XDONT(ICRS) < XWRK1) Then
XWRK = XDONT (ICRS)
Do IDCR = NORD - 1, MAX (1, ILOW) , - 1
If (XWRK >= XWRKT(IDCR)) Exit
XWRKT (IDCR+1) = XWRKT (IDCR)
End Do
XWRKT (IDCR+1) = XWRK
XWRK1 = XWRKT(NORD)
End If
ILOW = ILOW + 1
End Do
FNDNTH = XWRK1
!
End Function I_fndnth
Subroutine D_indmed (XDONT, INDM)
! Returns index of median value of XDONT.
! __________________________________________________________
real(kind=dp), Dimension (:), Intent (In) :: XDONT
Integer(kind=i4), Intent (Out) :: INDM
! __________________________________________________________
Integer(kind=i4) :: IDON
!
Allocate (IDONT (SIZE(XDONT)))
Do IDON = 1, SIZE(XDONT)
IDONT (IDON) = IDON
End Do
!
Call d_med (XDONT, IDONT, INDM)
!
Deallocate (IDONT)
End Subroutine D_indmed
Recursive Subroutine d_med (XDATT, IDATT, ires_med)
! Finds the index of the median of XDONT using the recursive procedure
! described in Knuth, The Art of Computer Programming,
! vol. 3, 5.3.3 - This procedure is linear in time, and
! does not require to be able to interpolate in the
! set as the one used in INDNTH. It also has better worst
! case behavior than INDNTH, but is about 30% slower in
! average for random uniformly distributed values.
! __________________________________________________________
real(kind=dp), Dimension (:), Intent (In) :: XDATT
Integer(kind=i4), Dimension (:), Intent (In) :: IDATT
Integer(kind=i4), Intent (Out):: ires_med
! __________________________________________________________
!
real(kind=dp), Parameter :: XHUGE = HUGE (XDATT)
real(kind=dp) :: XWRK, XWRK1, XMED7, XMAX, XMIN
!
Integer(kind=i4), Dimension (7*(((Size (IDATT)+6)/7+6)/7)) :: ISTRT, IENDT, IMEDT
Integer(kind=i4), Dimension (7*((Size(IDATT)+6)/7)) :: IWRKT
Integer(kind=i4) :: NTRI, NMED, NORD, NEQU, NLEQ, IMED, IDON, IDON1
Integer(kind=i4) :: IDEB, ITMP, IDCR, ICRS, ICRS1, ICRS2, IMAX, IMIN
Integer(kind=i4) :: IWRK, IWRK1, IMED1, IMED7, NDAT
!
NDAT = Size (IDATT)
NMED = (NDAT+1) / 2
IWRKT = IDATT
!
! If the number of values is small, then use insertion sort
!
If (NDAT < 35) Then
!
! Bring minimum to first location to save test in decreasing loop
!
IDCR = NDAT
If (XDATT (IWRKT (1)) < XDATT (IWRKT (IDCR))) Then
IWRK = IWRKT (1)
Else
IWRK = IWRKT (IDCR)
IWRKT (IDCR) = IWRKT (1)
Endif
XWRK = XDATT (IWRK)
Do ITMP = 1, NDAT - 2
IDCR = IDCR - 1
IWRK1 = IWRKT (IDCR)
XWRK1 = XDATT (IWRK1)
If (XWRK1 < XWRK) Then
IWRKT (IDCR) = IWRK
XWRK = XWRK1
IWRK = IWRK1
Endif
End Do
IWRKT (1) = IWRK
!
! Sort the first half, until we have NMED sorted values
!
Do ICRS = 3, NMED
XWRK = XDATT (IWRKT (ICRS))
IWRK = IWRKT (ICRS)
IDCR = ICRS - 1
Do
If (XWRK >= XDATT (IWRKT(IDCR))) Exit
IWRKT (IDCR+1) = IWRKT (IDCR)
IDCR = IDCR - 1
End Do
IWRKT (IDCR+1) = IWRK
End Do
!
! Insert any value less than the current median in the first half
!
XWRK1 = XDATT (IWRKT (NMED))
Do ICRS = NMED+1, NDAT
XWRK = XDATT (IWRKT (ICRS))
IWRK = IWRKT (ICRS)
If (XWRK < XWRK1) Then
IDCR = NMED - 1
Do
If (XWRK >= XDATT (IWRKT(IDCR))) Exit
IWRKT (IDCR+1) = IWRKT (IDCR)
IDCR = IDCR - 1
End Do
IWRKT (IDCR+1) = IWRK
XWRK1 = XDATT (IWRKT (NMED))
End If
End Do
ires_med = IWRKT (NMED)
Return
End If
!
! Make sorted subsets of 7 elements
! This is done by a variant of insertion sort where a first
! pass is used to bring the smallest element to the first position
! decreasing disorder at the same time, so that we may remove
! remove the loop test in the insertion loop.
!
IMAX = 1
IMIN = 1
XMAX = XDATT (IWRKT(IMAX))
XMIN = XDATT (IWRKT(IMIN))
DO IDEB = 1, NDAT-6, 7
IDCR = IDEB + 6
If (XDATT (IWRKT(IDEB)) < XDATT (IWRKT(IDCR))) Then
IWRK = IWRKT(IDEB)
Else
IWRK = IWRKT (IDCR)
IWRKT (IDCR) = IWRKT(IDEB)
Endif
XWRK = XDATT (IWRK)
Do ITMP = 1, 5
IDCR = IDCR - 1
IWRK1 = IWRKT (IDCR)
XWRK1 = XDATT (IWRK1)
If (XWRK1 < XWRK) Then
IWRKT (IDCR) = IWRK
IWRK = IWRK1
XWRK = XWRK1
Endif
End Do
IWRKT (IDEB) = IWRK
If (XWRK < XMIN) Then
IMIN = IWRK
XMIN = XWRK
End If
Do ICRS = IDEB+1, IDEB+5
IWRK = IWRKT (ICRS+1)
XWRK = XDATT (IWRK)
IDON = IWRKT(ICRS)
If (XWRK < XDATT(IDON)) Then
IWRKT (ICRS+1) = IDON
IDCR = ICRS
IWRK1 = IWRKT (IDCR-1)
XWRK1 = XDATT (IWRK1)
Do
If (XWRK >= XWRK1) Exit
IWRKT (IDCR) = IWRK1
IDCR = IDCR - 1
IWRK1 = IWRKT (IDCR-1)
XWRK1 = XDATT (IWRK1)
End Do
IWRKT (IDCR) = IWRK
EndIf
End Do
If (XWRK > XMAX) Then
IMAX = IWRK
XMAX = XWRK
End If
End Do
!
! Add-up alternatively MAX and MIN values to make the number of data
! an exact multiple of 7.
!
IDEB = 7 * (NDAT/7)
NTRI = NDAT
If (IDEB < NDAT) Then
!
Do ICRS = IDEB+1, NDAT
XWRK1 = XDATT (IWRKT (ICRS))
IF (XWRK1 > XMAX) Then
IMAX = IWRKT (ICRS)
XMAX = XWRK1
End If
IF (XWRK1 < XMIN) Then
IMIN = IWRKT (ICRS)
XMIN = XWRK1
End If
End Do
IWRK1 = IMAX
Do ICRS = NDAT+1, IDEB+7
IWRKT (ICRS) = IWRK1
If (IWRK1 == IMAX) Then
IWRK1 = IMIN
Else
NMED = NMED + 1
IWRK1 = IMAX
End If
End Do
!
Do ICRS = IDEB+2, IDEB+7
IWRK = IWRKT (ICRS)
XWRK = XDATT (IWRK)
Do IDCR = ICRS - 1, IDEB+1, - 1
If (XWRK >= XDATT (IWRKT(IDCR))) Exit
IWRKT (IDCR+1) = IWRKT (IDCR)
End Do
IWRKT (IDCR+1) = IWRK
End Do
!
NTRI = IDEB+7
End If
!
! Make the set of the indices of median values of each sorted subset
!
IDON1 = 0
Do IDON = 1, NTRI, 7
IDON1 = IDON1 + 1
IMEDT (IDON1) = IWRKT (IDON + 3)
End Do
!
! Find XMED7, the median of the medians
!
Call d_med (XDATT, IMEDT(1:IDON1), IMED7)
XMED7 = XDATT (IMED7)
!
! Count how many values are not higher than (and how many equal to) XMED7
! This number is at least 4 * 1/2 * (N/7) : 4 values in each of the
! subsets where the median is lower than the median of medians. For similar
! reasons, we also have at least 2N/7 values not lower than XMED7. At the
! same time, we find in each subset the index of the last value < XMED7,
! and that of the first > XMED7. These indices will be used to restrict the
! search for the median as the Kth element in the subset (> or <) where
! we know it to be.
!
IDON1 = 1
NLEQ = 0
NEQU = 0
Do IDON = 1, NTRI, 7
IMED = IDON+3
If (XDATT (IWRKT (IMED)) > XMED7) Then
IMED = IMED - 2
If (XDATT (IWRKT (IMED)) > XMED7) Then
IMED = IMED - 1
Else If (XDATT (IWRKT (IMED)) < XMED7) Then
IMED = IMED + 1
Endif
Else If (XDATT (IWRKT (IMED)) < XMED7) Then
IMED = IMED + 2
If (XDATT (IWRKT (IMED)) > XMED7) Then
IMED = IMED - 1
Else If (XDATT (IWRKT (IMED)) < XMED7) Then
IMED = IMED + 1
Endif
Endif
If (XDATT (IWRKT (IMED)) > XMED7) Then
NLEQ = NLEQ + IMED - IDON
IENDT (IDON1) = IMED - 1
ISTRT (IDON1) = IMED
Else If (XDATT (IWRKT (IMED)) < XMED7) Then
NLEQ = NLEQ + IMED - IDON + 1
IENDT (IDON1) = IMED
ISTRT (IDON1) = IMED + 1
Else ! If (XDATT (IWRKT (IMED)) == XMED7)
NLEQ = NLEQ + IMED - IDON + 1
NEQU = NEQU + 1
IENDT (IDON1) = IMED - 1
Do IMED1 = IMED - 1, IDON, -1
If (eq(XDATT (IWRKT (IMED1)) , XMED7)) Then
NEQU = NEQU + 1
IENDT (IDON1) = IMED1 - 1
Else
Exit
End If
End Do
ISTRT (IDON1) = IMED + 1
Do IMED1 = IMED + 1, IDON + 6
If (eq(XDATT (IWRKT (IMED1)) , XMED7)) Then
NEQU = NEQU + 1
NLEQ = NLEQ + 1
ISTRT (IDON1) = IMED1 + 1
Else
Exit
End If
End Do
Endif
IDON1 = IDON1 + 1
End Do
!
! Carry out a partial insertion sort to find the Kth smallest of the
! large values, or the Kth largest of the small values, according to
! what is needed.
!
!
If (NLEQ - NEQU + 1 <= NMED) Then