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