Skip to content

Commit

Permalink
upate to v0.8
Browse files Browse the repository at this point in the history
  • Loading branch information
samotracio committed Feb 28, 2018
1 parent afc145f commit a4b62a9
Show file tree
Hide file tree
Showing 8 changed files with 235 additions and 415 deletions.
108 changes: 99 additions & 9 deletions cflibfor.f90
Original file line number Diff line number Diff line change
Expand Up @@ -437,7 +437,7 @@ subroutine rppi_A(nt,npt,dec,dc,x,y,z,nsepp,sepp,nsepv,sepv,sbound, &
!$omp parallel do reduction(+:aapv) default(shared) &
!$omp& private(iq1,iq2,iq3,jq1,jq2,jq3,rcl,stm2,dltdec,jq1m,jq1min) &
!$omp& private(dltra,jq2m,jq2min,jq2max,jq2t,j,i,xi,yi,zi,dci,rv,rp2,ii,jj) &
!$omp& schedule(guided) if(nthr>1)
!$omp& schedule(dynamic) if(nthr>1)
do iq1=1,nc1
!$omp critical
write(*,fmt="(i4)",advance='no') iq1 ! for screen
Expand Down Expand Up @@ -2080,7 +2080,15 @@ subroutine s_A(nt,npt,dec,dc,x,y,z,nseps,seps,sbound,mxh1,mxh2,mxh3,cntid,logf,s
aas(nseps-1) = aas(nseps-1) + 1.0d0
goto 74
endif
do ii=nseps-2,1,-1
if(r2>seps2(nseps-2)) then
aas(nseps-2) = aas(nseps-2) + 1.0d0
goto 74
endif
if(r2>seps2(nseps-3)) then
aas(nseps-3) = aas(nseps-3) + 1.0d0
goto 74
endif
do ii=nseps-4,1,-1
if(r2>seps2(ii)) then
aas(ii) = aas(ii) + 1.0d0
goto 74
Expand Down Expand Up @@ -2271,7 +2279,15 @@ subroutine s_A_wg(nt,npt,dec,dc,wei,x,y,z,nseps,seps,sbound,mxh1,mxh2,mxh3,&
aas(nseps-1) = aas(nseps-1) + wpp
goto 74
endif
do ii=nseps-2,1,-1
if(r2>seps2(nseps-2)) then
aas(nseps-2) = aas(nseps-2) + wpp
goto 74
endif
if(r2>seps2(nseps-3)) then
aas(nseps-3) = aas(nseps-3) + wpp
goto 74
endif
do ii=nseps-4,1,-1
if(r2>seps2(ii)) then
aas(ii) = aas(ii) + wpp
goto 74
Expand Down Expand Up @@ -2465,7 +2481,17 @@ subroutine s_Ab(nt,npt,dec,dc,x,y,z,nseps,seps,sbound,mxh1,mxh2,mxh3,nbts,bseed,
baas(:,nseps-1) = baas(:,nseps-1) + wbts(:,i)*wbts(:,j)
goto 74
endif
do ii=nseps-2,1,-1
if(r2>seps2(nseps-2)) then
aas(nseps-2) = aas(nseps-2) + 1.0d0
baas(:,nseps-2) = baas(:,nseps-2) + wbts(:,i)*wbts(:,j)
goto 74
endif
if(r2>seps2(nseps-3)) then
aas(nseps-3) = aas(nseps-3) + 1.0d0
baas(:,nseps-3) = baas(:,nseps-3) + wbts(:,i)*wbts(:,j)
goto 74
endif
do ii=nseps-4,1,-1
if(r2>seps2(ii)) then
aas(ii) = aas(ii) + 1.0d0
baas(:,ii) = baas(:,ii) + wbts(:,i)*wbts(:,j)
Expand Down Expand Up @@ -2669,7 +2695,17 @@ subroutine s_Ab_wg(nt,npt,dec,dc,wei,x,y,z,nseps,seps,sbound,mxh1,mxh2,mxh3,&
baas(:,nseps-1) = baas(:,nseps-1) + wpp*wbts(:,i)*wbts(:,j)
goto 74
endif
do ii=nseps-2,1,-1
if(r2>seps2(nseps-2)) then
aas(nseps-2) = aas(nseps-2) + wpp
baas(:,nseps-2) = baas(:,nseps-2) + wpp*wbts(:,i)*wbts(:,j)
goto 74
endif
if(r2>seps2(nseps-3)) then
aas(nseps-3) = aas(nseps-3) + wpp
baas(:,nseps-3) = baas(:,nseps-3) + wpp*wbts(:,i)*wbts(:,j)
goto 74
endif
do ii=nseps-4,1,-1
if(r2>seps2(ii)) then
aas(ii) = aas(ii) + wpp
baas(:,ii) = baas(:,ii) + wbts(:,i)*wbts(:,j)
Expand Down Expand Up @@ -2855,7 +2891,19 @@ subroutine s_C(nt,npt,ra,dec,dc,x,y,z,npt1,dc1,x1,y1,z1, &
cds(nseps) = cds(nseps) + 1.0d0
goto 70
endif
do ii=nseps-1,1,-1
if(r2>seps2(nseps-1)) then
cds(nseps-1) = cds(nseps-1) + 1.0d0
goto 70
endif
if(r2>seps2(nseps-2)) then
cds(nseps-2) = cds(nseps-2) + 1.0d0
goto 70
endif
if(r2>seps2(nseps-3)) then
cds(nseps-3) = cds(nseps-3) + 1.0d0
goto 70
endif
do ii=nseps-4,1,-1
if(r2>seps2(ii)) then
cds(ii) = cds(ii) + 1.0d0
goto 70
Expand Down Expand Up @@ -3045,7 +3093,19 @@ subroutine s_C_wg(nt,npt,ra,dec,dc,wei,x,y,z,npt1,dc1,wei1,x1,y1,z1, &
cds(nseps) = cds(nseps) + wpp
goto 70
endif
do ii=nseps-1,1,-1
if(r2>seps2(nseps-1)) then
cds(nseps-1) = cds(nseps-1) + wpp
goto 70
endif
if(r2>seps2(nseps-2)) then
cds(nseps-2) = cds(nseps-2) + wpp
goto 70
endif
if(r2>seps2(nseps-3)) then
cds(nseps-3) = cds(nseps-3) + wpp
goto 70
endif
do ii=nseps-4,1,-1
if(r2>seps2(ii)) then
cds(ii) = cds(ii) + wpp
goto 70
Expand Down Expand Up @@ -3237,7 +3297,22 @@ subroutine s_Cb(nt,npt,ra,dec,dc,x,y,z,npt1,dc1,x1,y1,z1, &
bcds(:,nseps) = bcds(:,nseps) + wbts(:,i)*wbts1(:,j)
goto 70
endif
do ii=nseps-1,1,-1
if(r2>seps2(nseps-1)) then
cds(nseps-1) = cds(nseps-1) + 1.0d0
bcds(:,nseps-1) = bcds(:,nseps-1) + wbts(:,i)*wbts1(:,j)
goto 70
endif
if(r2>seps2(nseps-2)) then
cds(nseps-2) = cds(nseps-2) + 1.0d0
bcds(:,nseps-2) = bcds(:,nseps-2) + wbts(:,i)*wbts1(:,j)
goto 70
endif
if(r2>seps2(nseps-3)) then
cds(nseps-3) = cds(nseps-3) + 1.0d0
bcds(:,nseps-3) = bcds(:,nseps-3) + wbts(:,i)*wbts1(:,j)
goto 70
endif
do ii=nseps-4,1,-1
if(r2>seps2(ii)) then
cds(ii) = cds(ii) + 1.0d0
bcds(:,ii) = bcds(:,ii) + wbts(:,i)*wbts1(:,j)
Expand Down Expand Up @@ -3440,7 +3515,22 @@ subroutine s_Cb_wg(nt,npt,ra,dec,dc,wei,x,y,z,npt1,dc1,wei1,x1,y1,z1,nseps,seps,
bcds(:,nseps) = bcds(:,nseps) + wpp*wbts(:,i)*wbts1(:,j)
goto 70
endif
do ii=nseps-1,1,-1
if(r2>seps2(nseps-1)) then
cds(nseps-1) = cds(nseps-1) + wpp
bcds(:,nseps-1) = bcds(:,nseps-1) + wpp*wbts(:,i)*wbts1(:,j)
goto 70
endif
if(r2>seps2(nseps-2)) then
cds(nseps-2) = cds(nseps-2) + wpp
bcds(:,nseps-2) = bcds(:,nseps-2) + wpp*wbts(:,i)*wbts1(:,j)
goto 70
endif
if(r2>seps2(nseps-3)) then
cds(nseps-3) = cds(nseps-3) + wpp
bcds(:,nseps-3) = bcds(:,nseps-3) + wpp*wbts(:,i)*wbts1(:,j)
goto 70
endif
do ii=nseps-4,1,-1
if(r2>seps2(ii)) then
cds(ii) = cds(ii) + wpp
bcds(:,ii) = bcds(:,ii) + wpp*wbts(:,i)*wbts1(:,j)
Expand Down
Binary file added docs/benchMP_plot1.png
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file modified docs/bench_rppiA_plot1.png
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file modified docs/bench_rppiA_plot2.png
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
8 changes: 7 additions & 1 deletion docs/performance.rst
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,13 @@ Now we compare with one of the fastest publicity available codes
.. image:: bench_rppiA_plot1.png
:scale: 60%
:alt: DD(rp,pi) performance1


And now we compare the performance with *n* threads running in parallel

.. image:: benchMP_plot1.png
:scale: 60%
:alt: DD(rp,pi) performance3

The equipment employed for testing is based in a four core i7-3770K 3.5GHz CPU
(L1 cache: 32KB data + 32KB instruction, L2 cache: 256KB, L3 cache: 8MB shared)
with 16GB RAM, running with OpenSuse Linux, GNU Fortran 6.1.1 compiler and Python
Expand Down
6 changes: 2 additions & 4 deletions docs/pxorder.rst
Original file line number Diff line number Diff line change
Expand Up @@ -20,11 +20,9 @@ list algorithm. The parameter :code:`pxorder` controls the method to sort the
pixels. Available options are:

- **'natural'** |br|
Pixels are sorted in lexical order by three keys, DEC-RA-Z, in that order.
Pixels are sorted in lexical order by three keys, Z-DEC-RA, in that order.
This is usually the fastest option as it more closely matches the mechanics
of the counting algorithm. The number of pixels are set to :code:`(mxh1,mxh2,10*mxh3)`,
where the factor of 10 creates a slighly more fine grained ordering along the
redshift dimension
of the counting algorithm. The number of pixels are set to :code:`(mxh3,mxh1,mxh2)`

- **'morton_dir'** |br|
2D Morton-curve ordering using the RA-DEC values directly, i.e. no division
Expand Down
8 changes: 4 additions & 4 deletions example_pcf.py
Original file line number Diff line number Diff line change
Expand Up @@ -10,10 +10,10 @@
ranf = './examples_data/RAND.fits' # Random sample
outfn = './examples_data/ex_pcf_01' # Name for output files
par = gun.packpars(kind='pcf', file=galf, file1=ranf, outfn=outfn)
par.autogrid = False # Automatic SK grid size
par.mxh1 = 60 # SK size in dec
par.mxh2 = 240 # SK size in ra
par.mxh3 = 24 # SK size in z
par.autogrid = True # Automatic SK grid size
par.mxh1 = 20 # SK size in dec
par.mxh2 = 100 # SK size in ra
par.mxh3 = 10 # SK size in z
par.nsepp = 28 # Number of bins of projected separation rp
par.seppmin = 0.02 # Minimum rp in Mpc/h
par.dsepp = 0.12 # Bin size of rp (in log space)
Expand Down

0 comments on commit a4b62a9

Please sign in to comment.