1
2! CONVENTIONS AND CHANGES:
3! 1) All the periods have been changed by underscores
4! 2) Instead using a derive type in fortran, fields of variable "vgm" has been split
5!    into two new arrays: "vgm_ibin" and "vgm_ipair"
6! 3) The integers n,p,n_out has been introduced to define the array sizes
7!	 B= size of gamma_hat
8!        n= size of vgm_pair
9!        p= size of rho_grid
10! 4) The output of the original R-function is included as the last argument in the subrutine
11!------------------- hg function ---------------------------------------
12
13subroutine hg(rho,m,value)
14
15 implicit none
16 integer,parameter:: dp = KIND(1.0d0)
17 integer,intent(in)::m
18 real(kind=dp),intent(in)::rho(m)
19 real(kind=dp)::value(m),a,b,cc,fn(m),fnold(m),facn
20 integer::n
21 real(kind=dp),external::fgamma
22
23 ! Note: the Gamma function was included as intrinsic in Gfortran 4.3
24
25   a     = 3.0_dp/4.0_dp
26   b     = 3.0_dp/4.0_dp
27   cc    = 1.0_dp/2.0_dp
28   fn    = fgamma(a)*fgamma(b)/fgamma(cc)
29   facn  = 1.0_dp
30   n     = 1
31   fnold = 0.1_dp
32
33   do while (maxval(abs(fn-fnold)/fnold) > 0.0001_dp)
34      facn  = facn*n
35      fnold = fn
36      fn    = fn + fgamma(a + n)*fgamma(b + n)*rho**n/(fgamma(cc+n)*facn)
37      n     = n + 1
38    end do
39
40   value=fn*fgamma(cc)/(fgamma(a)*fgamma(b))
41
42end subroutine hg
43
44!------------------- cor_sqrtabs function ------------------------------
45subroutine cor_sqrtabs(rho,m,value)
46
47 implicit none
48
49 interface
50   subroutine hg(rho,m,value)
51      integer,parameter:: dp = KIND(1.0d0)
52      integer,intent(in)::m
53      real(kind=dp),intent(in)::rho(m)
54      real(kind=dp)::value(m)
55   end subroutine hg
56 end interface
57
58 integer,parameter:: dp = KIND(1.0d0)
59 integer,intent(in)::m
60 real(kind=dp),intent(in)::rho(m)
61 real(kind=dp)::pi,value(m),val(m)
62 real(kind=dp),external::fgamma
63 pi=acos(-1.0_dp)
64
65
66 ! Note: the Gamma function was included as intrinsi! in Gfortran 4.3
67 call hg(rho**2,m,val)
68 value=fgamma(0.75_dp)**2*((1.0_dp-rho**2)*val-1.0_dp)/(sqrt(pi)-fgamma(0.75_dp)**2)
69
70
71end subroutine cor_sqrtabs
72
73!------------------- approx_linear function ------------------------------
74subroutine approx_linear(x,y,n,v,m,yleft,yright,value)
75
76 implicit none
77 integer,parameter:: dp = KIND(1.0d0)
78 integer,intent(in)::n,m
79 integer::j,l
80 real(kind=dp)::value(m)
81 real(kind=dp),intent(in)::x(n),y(n),v(m),yleft,yright
82
83 ! Loop in the components of v
84 do j=1,m
85
86    if(v(j)<=x(1)) then
87       value(j)=yleft
88    elseif(v(j)>=x(n)) then
89       value(j)=yright
90    else
91       l=count(v(j)>x)
92       value(j)=y(l)+(y(l+1)-y(l))*((v(j)-x(l))/(x(l+1)-x(l)))
93    endif
94
95 enddo
96end subroutine approx_linear
97
98
99subroutine cov_bin_fun(B,n,p,i,j,vgm_ibin,vgm_ipair,gamma_hat,mean_cv)
100
101 implicit none
102
103 interface
104   subroutine approx_linear(x,y,n,v,m,yleft,yright,value)
105      integer,parameter:: dp = KIND(1.0d0)
106      integer,intent(in)::n,m
107      real(kind=dp)::value(m)
108      real(kind=dp),intent(in)::x(n),y(n),v(m),yleft,yright
109   end subroutine approx_linear
110 end interface
111
112 interface
113   subroutine cor_sqrtabs(rho,m,value)
114     integer,parameter:: dp = KIND(1.0d0)
115     integer,intent(in)::m
116     real(kind=dp),intent(in)::rho(m)
117     real(kind=dp)::value(m)
118   end subroutine cor_sqrtabs
119 end interface
120
121! Simple precision
122!  integer,parameter:: dp = KIND(1.0)
123! Double precision
124  integer,parameter:: dp = KIND(1.0d0)
125
126 integer,intent(in)::B,n,p,i,j
127 integer,intent(in):: vgm_ipair(n,2),vgm_ibin(n)
128 real(kind=dp),intent(in):: gamma_hat(B)
129 real(kind=dp),intent(inout):: mean_cv
130
131 integer::iter,k,l,m,m1,m2,m3,m4,n_bin_i,n_bin_j,nmax,ib,jb,n_out,unique_vgm_ibin(B)
132 integer,dimension(:),allocatable::vrow1,vrow2,aux, &
133                      vind1,vind2,vind3,vind4,ind,ind1,ind2,ind3,ind4
134 integer,dimension(:,:),allocatable::aux1,aux2,mp,mp1,mp2, &
135                      mpair1(:,:),mpair2(:,:),mpair3(:,:),mpair4(:,:)
136
137 real(kind=dp):: rho_grid(p+1),f_grid(p+1)
138 real(kind=dp),allocatable::cr(:),cv(:),tmp(:),mpair1_hat(:),mpair2_hat(:), &
139                                    mpair3_hat(:),mpair4_hat(:),auxcv(:)
140
141 logical,allocatable::ind_bin_i(:),ind_bin_j(:), &
142                      mask(:),mask1(:),mask2(:),mask3(:),mask4(:)
143
144 ! Memory allocation for ind_bin_i,ind_bin_j
145 allocate(ind_bin_i(n),ind_bin_j(n))
146
147! ------------ rho evaluation ------------------------------------------
148 k=0; rho_grid(1:p) = (/((0.96_dp/(p-1))*(k-1),k=1,p)/)
149 call cor_sqrtabs(rho_grid(1:p),p,f_grid(1:p))
150 rho_grid(p+1) = 1.0_dp
151 f_grid(p+1) = 1.0_dp
152
153! ------------- Sort computations --------------------------------------
154 ind_bin_i = (vgm_ibin == i)
155 ind_bin_j = (vgm_ibin == j)
156 n_bin_i = count(ind_bin_i)
157 n_bin_j = count(ind_bin_j)
158
159 ! Memory allocation for ind_bin_i,ind_bin_j
160 n_out=n_bin_i*n_bin_j
161 allocate(vrow1(n_out),vrow2(n_out))
162
163 ! Computing the values for vrow1,vrow2
164 do k=1,n_bin_j
165    vrow1(1+(k-1)*n_bin_i:k*n_bin_i) = (/(l,l=1,n_bin_i)/)
166 enddo
167 do k=1,n_bin_j
168    vrow2(1+(k-1)*n_bin_i:k*n_bin_i) = k
169 enddo
170
171 ! Memory allocation for mp1,mp2 and mp
172 allocate(aux1(n_bin_i,2),aux2(n_bin_j,2), &
173          mp1(n_out,2),mp2(n_out,2))
174
175 ! Construction mp1
176 aux1=vgm_ipair(pack((/(l,l=1,n)/),ind_bin_i),:)
177 mp1=aux1(vrow1,:)
178
179 ! Construction mp2
180 aux2=vgm_ipair(pack((/(l,l=1,n)/),ind_bin_j),:)
181 mp2=aux2(vrow2,:)
182
183 ! Construction mp
184 allocate(mp(n_out,4))
185 mp(:,1:2)=mp1
186 mp(:,3:4)=mp2
187
188 ! Memory deallocation
189 deallocate(aux1,aux2,vrow1,vrow2,mp1,mp2,ind_bin_i,ind_bin_j)
190
191!---------------- mpair* matrices: --------------------------------------
192! 1: pair(1)
193! 2: pair(2)
194! 3: pair(1)-pair(2)
195! 4: index
196! 5: gamma_hat, if different pairs --> go to mpair*_hat(:)
197
198 ! Memory allocation
199 allocate(mpair1(n_out,4),mpair2(n_out,4), &
200          mpair3(n_out,4),mpair4(n_out,4))
201
202 ! Initialize to zero
203 mpair1 = 0; mpair2 = 0; mpair3 = 0; mpair4 = 0
204
205 ! Compute the values for mpair*(:,1:3)
206 mpair1(:,1)= merge(mp(:,2),mp(:,3),mp(:,2)<=mp(:,3))
207 mpair1(:,2)= merge(mp(:,2),mp(:,3),mp(:,2)>=mp(:,3))
208 mpair2(:,1)= merge(mp(:,1),mp(:,4),mp(:,1)<=mp(:,4))
209 mpair2(:,2)= merge(mp(:,1),mp(:,4),mp(:,1)>=mp(:,4))
210 mpair3(:,1)= merge(mp(:,1),mp(:,3),mp(:,1)<=mp(:,3))
211 mpair3(:,2)= merge(mp(:,1),mp(:,3),mp(:,1)>=mp(:,3))
212 mpair4(:,1)= merge(mp(:,2),mp(:,4),mp(:,2)<=mp(:,4))
213 mpair4(:,2)= merge(mp(:,2),mp(:,4),mp(:,2)>=mp(:,4))
214
215 mpair1(:,3)= mpair1(:,2)-mpair1(:,1)
216 mpair2(:,3)= mpair2(:,2)-mpair2(:,1)
217 mpair3(:,3)= mpair3(:,2)-mpair3(:,1)
218 mpair4(:,3)= mpair4(:,2)-mpair4(:,1)
219
220 ! Memory allocation for temporal arrays
221 allocate(mask1(n_out),mask2(n_out),aux(n_out),mask3(n_out),mask4(n_out))
222
223 ! Compute the mask for the mpair* (the pairs outside the diagonal)
224 mask1=(mpair1(:,3) /= 0)
225 mask2=(mpair2(:,3) /= 0)
226 mask3=(mpair3(:,3) /= 0)
227 mask4=(mpair4(:,3) /= 0)
228
229 m1=count(mask1)
230 m2=count(mask2)
231 m3=count(mask3)
232 m4=count(mask4)
233
234 ! Memory allocation for vind*
235 allocate(vind1(m1),vind2(m2),vind3(m3),vind4(m4))
236
237 aux=(/(k,k=1,n_out)/)
238 vind1= pack(aux,mask1)
239 vind2= pack(aux,mask2)
240 vind3= pack(aux,mask3)
241 vind4= pack(aux,mask4)
242
243 ! Memory deallocation
244 deallocate(mask1,mask2,mask3,mask4)
245
246 nmax= maxval(vgm_ipair)
247 mpair1(vind1,4) = (mpair1(vind1,1)-1)*nmax- &
248                   ((mpair1(vind1,1)-1)*mpair1(vind1,1))/2 + &
249                    mpair1(vind1,2)-mpair1(vind1,1)
250
251 mpair2(vind2,4) = (mpair2(vind2,1)-1)*nmax- &
252                   ((mpair2(vind2,1)-1)*mpair2(vind2,1))/2 + &
253                    mpair2(vind2,2)-mpair2(vind2,1)
254
255 mpair3(vind3,4) = (mpair3(vind3,1)-1)*nmax- &
256                   ((mpair3(vind3,1)-1)*mpair3(vind3,1))/2 + &
257                    mpair3(vind3,2)-mpair3(vind3,1)
258
259 mpair4(vind4,4) = (mpair4(vind4,1)-1)*nmax- &
260                   ((mpair4(vind4,1)-1)*mpair4(vind4,1))/2 + &
261                    mpair4(vind4,2)-mpair4(vind4,1)
262
263 ! For gamma_hat at each bin
264 allocate(ind1(m1),ind2(m2),ind3(m3),ind4(m4))
265
266 ! Construction of the uniquenness for vgm_ibin
267 unique_vgm_ibin=maxval(vgm_ibin)+1
268 iter=1
269 do k=1,int(maxval(vgm_ibin))
270     if(any(vgm_ibin==k)) then
271        unique_vgm_ibin(iter)=k
272        iter=iter+1
273     endif
274 enddo
275
276 do k=1,m1
277   ind1(k)=count(unique_vgm_ibin<vgm_ibin(mpair1(vind1(k),4)))+1
278 enddo
279 do k=1,m2
280   ind2(k)=count(unique_vgm_ibin<vgm_ibin(mpair2(vind2(k),4)))+1
281 enddo
282 do k=1,m3
283   ind3(k)=count(unique_vgm_ibin<vgm_ibin(mpair3(vind3(k),4)))+1
284 enddo
285 do k=1,m4
286   ind4(k)=count(unique_vgm_ibin<vgm_ibin(mpair4(vind4(k),4)))+1
287 enddo
288
289 ib = count(unique_vgm_ibin<i)+1
290 jb = count(unique_vgm_ibin<j)+1
291
292 ! Memory deallocation
293 deallocate(mpair1,mpair2,mpair3,mpair4)
294
295 ! Memory allocation
296 allocate(cr(n_out),cv(n_out),tmp(n_out),mpair1_hat(n_out),mpair2_hat(n_out), &
297                      mpair3_hat(n_out),mpair4_hat(n_out))
298
299 ! Initialize to zero
300 mpair1_hat = 0.0_dp; mpair2_hat = 0.0_dp
301 mpair3_hat = 0.0_dp; mpair4_hat = 0.0_dp
302
303 mpair1_hat(vind1) = gamma_hat(ind1)
304 mpair2_hat(vind2) = gamma_hat(ind2)
305 mpair3_hat(vind3) = gamma_hat(ind3)
306 mpair4_hat(vind4) = gamma_hat(ind4)
307
308 ! Memory deallocation
309 deallocate(vind1,vind2,vind3,vind4,ind1,ind2,ind3,ind4)
310
311 ! Compute cr
312 cr = (-mpair1_hat-mpair2_hat+mpair3_hat+mpair4_hat)/ &
313         (2.0_dp*sqrt(abs(gamma_hat(ib)*gamma_hat(jb))))!<--- (2.0_dp*sqrt((gamma_hat(ib)*gamma_hat(jb))))
314
315 tmp=merge(abs(cr),(/(1.0_dp,k=1,n_out)/),abs(cr)<1.0_dp)
316 cr=merge(tmp,-tmp,cr>=0.0_dp)
317
318 ! Memory deallocation
319 deallocate(mpair1_hat,mpair2_hat,mpair3_hat,mpair4_hat)
320
321 ! Compute mult (using the function "approx" to compute a linear interpolation)
322 allocate(mask(n_out))
323
324 ! Memory allocation and compute the index array "ind"
325 mask=(abs(cr)<=0.96_dp)
326 m=count(mask)
327 allocate(ind(m),auxcv(m))
328 ind=pack(aux,mask)
329 cv=1.0_dp
330 call approx_linear(rho_grid,f_grid,p+1,abs(cr(ind)),m,f_grid(1),f_grid(p+1),auxcv)
331 cv(ind)=auxcv
332 !cv=cv*0.172402_dp*sqrt(sqrt((gamma_hat(ib)* gamma_hat(jb)))) !<-------------
333 cv=cv*0.172402_dp*sqrt(sqrt(abs(gamma_hat(ib)* gamma_hat(jb))))
334 mean_cv=sum(cv)/n_out
335
336 ! Memory deallocation
337 deallocate(mask,ind,cr,cv,aux,auxcv)
338
339end subroutine cov_bin_fun
340
341
342subroutine diag_cov_bin_fun(B,n,p,vgm_ibin,vgm_ipair,gamma_hat,mean_cv)
343
344 implicit none
345
346 interface
347   subroutine cov_bin_fun(B,n,p,i,j,vgm_ibin,vgm_ipair,gamma_hat,mean_cv)
348     integer,parameter:: dp = KIND(1.0d0)
349     integer,intent(in)::B,n,p,i,j
350     integer,intent(in):: vgm_ipair(n,2),vgm_ibin(n)
351     real(kind=dp),intent(in):: gamma_hat(B)
352     real(kind=dp),intent(inout):: mean_cv
353   end subroutine cov_bin_fun
354 end interface
355
356 integer,parameter:: dp = KIND(1.0d0)
357 integer,intent(in)::B,n,p
358 integer,intent(in):: vgm_ipair(n,2),vgm_ibin(n)
359 real(kind=dp),intent(in):: gamma_hat(B)
360 real(kind=dp),intent(inout):: mean_cv(B)
361 integer:: i
362
363 do i=1,B
364   call cov_bin_fun(B,N,p,i,i,vgm_ibin,vgm_ipair,gamma_hat,mean_cv(i))
365 enddo
366
367end subroutine diag_cov_bin_fun
368
369
370subroutine full_cov_bin_fun(B,n,p,vgm_ibin,vgm_ipair,gamma_hat,mean_cv)
371
372 implicit none
373
374 interface
375   subroutine cov_bin_fun(B,n,p,i,j,vgm_ibin,vgm_ipair,gamma_hat,mean_cv)
376     integer,parameter:: dp = KIND(1.0d0)
377     integer,intent(in)::B,n,p,i,j
378     integer,intent(in):: vgm_ipair(n,2),vgm_ibin(n)
379     real(kind=dp),intent(in):: gamma_hat(B)
380     real(kind=dp),intent(inout):: mean_cv
381   end subroutine cov_bin_fun
382 end interface
383
384 integer,parameter:: dp = KIND(1.0d0)
385 integer,intent(in)::B,n,p
386 integer,intent(in):: vgm_ipair(n,2),vgm_ibin(n)
387 real(kind=dp),intent(in):: gamma_hat(B)
388 real(kind=dp),intent(inout):: mean_cv(B,B)
389 integer:: i,j
390
391 do i=1,B
392   do j=i,B
393     call cov_bin_fun(B,N,p,i,j,vgm_ibin,vgm_ipair,gamma_hat,mean_cv(i,j))
394     mean_cv(j,i)=mean_cv(i,j)
395   enddo
396 enddo
397
398end subroutine full_cov_bin_fun
399
400
401
402