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
|
subroutine dsearchc(X, m, val, n, indX, occ, info)
*
*
* PURPOSE
* val(0..n) being an array (with strict increasing order and n >=1)
* representing intervals, this routine, by the mean of a
* dichotomic search, computes :
*
* a/ for each X(i) its interval number indX(i) :
* indX(i) = j if X(i) in (val(j-1), val(j)]
* = 1 if X(i) = val(0)
* = 0 if X(i) is not in [val(0),val(n)]
*
* b/ the number of points falling in the interval j :
*
* occ(j) = # { X(i) such that X(i) in (val(j-1), val(j)] } for j>1
* and occ(1) = # { X(i) such that X(i) in [val(0), val(1)] }
*
* PARAMETERS
* inputs :
* m integer
* X(1..m) double float array
* n integer
* val(0..n) double float array (val(0) < val(1) < ....)
* outputs
* indX(1..m) integer array
* occ(1..n) integer array
* info integer (number of X(i) not in [val(0), val(n)])
*
* AUTHOR
* Bruno Pincon
*
implicit none
integer m, n, info
double precision X(m), val(0:n)
integer occ(n), indX(m)
integer i, j1, j2, j
do j = 1, n
occ(j) = 0
enddo
info = 0
do i = 1, m
if ( val(0) .le. X(i) .and. X(i) .le. val(n) ) then
* X(i) is in [val(0),val(n)] :
* find j such that val(j-1) <= X(i) <= val(j) by a dicho search
j1 = 0
j2 = n
do while ( j2 - j1 .gt. 1 )
j = (j1 + j2)/2
if ( X(i) .le. val(j) ) then
j2 = j
else
j1 = j
endif
enddo
* we have val(j1) < X(i) <= val(j2) if j2 > 1 (j1=j2-1)
* or val(j1) <= X(i) <= val(j2) if j2 = 1 (j1=j2-1)
* so that j2 is the good interval number in all cases
occ(j2) = occ(j2) + 1
indX(i) = j2
else ! X(i) is not in [val(0), val(n)]
info = info + 1
indX(i) = 0
endif
enddo
end
*
**************************************************************************
*
subroutine dsearchd(X, m, val, n, indX, occ, info)
*
* PURPOSE
* val(1..n) being a strictly increasing array, this
* routines by the mean of a dichotomic search computes :
*
* a/ the number of occurences (occ(j)) of each value val(j)
* in the array X :
*
* occ(j) = #{ X(i) such that X(i) = val(j) }
*
* b/ the array indX : if X(i) = val(j) then indX(i) = j
* (if X(i) is not in val then indX(i) = 0)
*
* PARAMETERS
* inputs :
* m integer
* X(1..m) double float array
* n integer
* val(1..n) double float array (must be in a strict increasing order)
* outputs :
* occ(1..n) integer array
* indX(1..m) integer array
* info integer (number of X(i) which are not in val(1..n))
*
* AUTHOR
* Bruno Pincon
*
implicit none
integer m, n, info
double precision X(m), val(n)
integer occ(n), indX(m)
integer i, j1, j2, j
do j = 1, n
occ(j) = 0
enddo
info = 0
do i = 1, m
if ( val(1) .le. X(i) .and. X(i) .le. val(n) ) then
* find j such that X(i) = val(j) by a dicho search
j1 = 1
j2 = n
do while ( j2 - j1 .gt. 1 )
j = (j1 + j2)/2
if ( X(i) .lt. val(j) ) then
j2 = j
else
j1 = j
endif
enddo
* here we know that val(j1) <= X(i) <= val(j2) with j2 = j1 + 1
* (in fact we have exactly val(j1) <= X(i) < val(j2) if j2 < n)
if (X(i) .eq. val(j1)) then
occ(j1) = occ(j1) + 1
indX(i) = j1
else if (X(i) .eq. val(j2)) then ! (note: this case may happen only for j2=n)
occ(j2) = occ(j2) + 1
indX(i) = j2
else ! X(i) is not in {val(1), val(2),..., val(n)}
info = info + 1
indX(i) = 0
endif
else ! X(i) is not in {val(1), val(2),..., val(n)}
info = info + 1
indX(i) = 0
endif
enddo
end
|
