Skip to content

Commit 9a839b8

Browse files
author
Dag Sverre Seljebotn
committed
Fortran wrapper for libsharp
1 parent c3b6277 commit 9a839b8

3 files changed

Lines changed: 317 additions & 0 deletions

File tree

‎fortran/sharp.f90‎

Lines changed: 245 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,245 @@
1+
module sharp
2+
use iso_c_binding
3+
implicit none
4+
! alm_info flags
5+
integer, parameter :: SHARP_PACKED = 1
6+
7+
! sharp job types
8+
enum, bind(c)
9+
enumerator :: SHARP_YtW = 0
10+
enumerator :: SHARP_Y = 1
11+
enumerator :: SHARP_Yt = 2
12+
enumerator :: SHARP_WY = 3
13+
enumerator :: SHARP_ALM2MAP_DERIV1 = 4
14+
end enum
15+
16+
! sharp job flags
17+
integer, parameter :: SHARP_DP = ISHFT(1, 4)
18+
integer, parameter :: SHARP_ADD = ISHFT(1, 5)
19+
integer, parameter :: SHARP_REAL_HARMONICS = ISHFT(1, 6)
20+
integer, parameter :: SHARP_NO_FFT = ISHFT(1, 7)
21+
22+
type sharp_geom_info
23+
type(c_ptr) :: handle
24+
integer(c_ptrdiff_t) :: n_local
25+
end type sharp_geom_info
26+
27+
type sharp_alm_info
28+
type(c_ptr) :: handle
29+
integer(c_ptrdiff_t) :: n_local
30+
end type sharp_alm_info
31+
32+
interface
33+
34+
! alm_info
35+
subroutine sharp_make_general_alm_info( &
36+
lmax, nm, stride, mval, mvstart, flags, alm_info) bind(c)
37+
use iso_c_binding
38+
integer(c_int), value, intent(in) :: lmax, nm, stride, flags
39+
integer(c_int), intent(in) :: mval(nm)
40+
integer(c_ptrdiff_t), intent(in) :: mvstart(nm)
41+
type(c_ptr), intent(out) :: alm_info
42+
end subroutine sharp_make_general_alm_info
43+
44+
subroutine c_sharp_make_mmajor_real_packed_alm_info( &
45+
lmax, stride, nm, ms, alm_info) bind(c, name='sharp_make_mmajor_real_packed_alm_info')
46+
use iso_c_binding
47+
integer(c_int), value, intent(in) :: lmax, nm, stride
48+
integer(c_int), intent(in), optional :: ms(nm)
49+
type(c_ptr), intent(out) :: alm_info
50+
end subroutine c_sharp_make_mmajor_real_packed_alm_info
51+
52+
function c_sharp_alm_count(alm_info) bind(c, name='sharp_alm_count')
53+
use iso_c_binding
54+
integer(c_ptrdiff_t) :: c_sharp_alm_count
55+
type(c_ptr), value, intent(in) :: alm_info
56+
end function c_sharp_alm_count
57+
58+
subroutine c_sharp_destroy_alm_info(alm_info) bind(c, name='sharp_destroy_alm_info')
59+
use iso_c_binding
60+
type(c_ptr), value :: alm_info
61+
end subroutine c_sharp_destroy_alm_info
62+
63+
! geom_info
64+
subroutine sharp_make_subset_healpix_geom_info ( &
65+
nside, stride, nrings, rings, weight, geom_info) bind(c)
66+
use iso_c_binding
67+
integer(c_int), value, intent(in) :: nside, stride, nrings
68+
integer(c_int), intent(in), optional :: rings(nrings)
69+
real(c_double), intent(in), optional :: weight(2 * nside)
70+
type(c_ptr), intent(out) :: geom_info
71+
end subroutine sharp_make_subset_healpix_geom_info
72+
73+
subroutine c_sharp_destroy_geom_info(geom_info) bind(c, name='sharp_destroy_geom_info')
74+
use iso_c_binding
75+
type(c_ptr), value :: geom_info
76+
end subroutine c_sharp_destroy_geom_info
77+
78+
function c_sharp_map_size(info) bind(c, name='sharp_map_size')
79+
use iso_c_binding
80+
integer(c_ptrdiff_t) :: c_sharp_map_size
81+
type(c_ptr), value :: info
82+
end function c_sharp_map_size
83+
84+
85+
! execute
86+
subroutine c_sharp_execute(type, spin, alm, map, geom_info, alm_info, ntrans, &
87+
flags, time, opcnt) bind(c, name='sharp_execute')
88+
use iso_c_binding
89+
integer(c_int), value :: type, spin, ntrans, flags
90+
type(c_ptr), value :: alm_info, geom_info
91+
real(c_double), intent(out), optional :: time
92+
integer(c_long_long), intent(out), optional :: opcnt
93+
type(c_ptr), intent(in) :: alm(*), map(*)
94+
end subroutine c_sharp_execute
95+
96+
subroutine c_sharp_execute_mpi(comm, type, spin, alm, map, geom_info, alm_info, ntrans, &
97+
flags, time, opcnt) bind(c, name='sharp_execute_mpi_fortran')
98+
use iso_c_binding
99+
integer(c_int), value :: comm, type, spin, ntrans, flags
100+
type(c_ptr), value :: alm_info, geom_info
101+
real(c_double), intent(out), optional :: time
102+
integer(c_long_long), intent(out), optional :: opcnt
103+
type(c_ptr), intent(in) :: alm(*), map(*)
104+
end subroutine c_sharp_execute_mpi
105+
106+
end interface
107+
108+
interface sharp_execute
109+
module procedure sharp_execute_d
110+
end interface
111+
112+
contains
113+
! alm info
114+
115+
! if ms is not passed, we default to using m=0..lmax.
116+
subroutine sharp_make_mmajor_real_packed_alm_info(lmax, ms, alm_info)
117+
use iso_c_binding
118+
integer(c_int), value, intent(in) :: lmax
119+
integer(c_int), intent(in), optional :: ms(:)
120+
type(sharp_alm_info), intent(out) :: alm_info
121+
!--
122+
integer(c_int), allocatable :: ms_copy(:)
123+
integer(c_int) :: nm
124+
125+
if (present(ms)) then
126+
nm = size(ms)
127+
allocate(ms_copy(nm))
128+
ms_copy = ms
129+
call c_sharp_make_mmajor_real_packed_alm_info(lmax, 1, nm, ms_copy, alm_info=alm_info%handle)
130+
deallocate(ms_copy)
131+
else
132+
call c_sharp_make_mmajor_real_packed_alm_info(lmax, 1, lmax + 1, alm_info=alm_info%handle)
133+
end if
134+
alm_info%n_local = c_sharp_alm_count(alm_info%handle)
135+
end subroutine sharp_make_mmajor_real_packed_alm_info
136+
137+
subroutine sharp_destroy_alm_info(alm_info)
138+
use iso_c_binding
139+
type(sharp_alm_info), intent(inout) :: alm_info
140+
call c_sharp_destroy_alm_info(alm_info%handle)
141+
alm_info%handle = c_null_ptr
142+
end subroutine sharp_destroy_alm_info
143+
144+
145+
! geom info
146+
subroutine sharp_make_healpix_geom_info(nside, rings, weight, geom_info)
147+
integer(c_int), value :: nside
148+
integer(c_int), optional :: rings(:)
149+
real(c_double), intent(in), optional :: weight(2 * nside)
150+
type(sharp_geom_info), intent(out) :: geom_info
151+
!--
152+
integer(c_int) :: nrings
153+
integer(c_int), allocatable :: rings_copy(:)
154+
155+
if (present(rings)) then
156+
nrings = size(rings)
157+
allocate(rings_copy(nrings))
158+
rings_copy = rings
159+
call sharp_make_subset_healpix_geom_info(nside, 1, nrings, rings_copy, &
160+
weight, geom_info%handle)
161+
deallocate(rings_copy)
162+
else
163+
call sharp_make_subset_healpix_geom_info(nside, 1, nrings=4 * nside - 1, &
164+
weight=weight, geom_info=geom_info%handle)
165+
end if
166+
geom_info%n_local = c_sharp_map_size(geom_info%handle)
167+
end subroutine sharp_make_healpix_geom_info
168+
169+
subroutine sharp_destroy_geom_info(geom_info)
170+
use iso_c_binding
171+
type(sharp_geom_info), intent(inout) :: geom_info
172+
call c_sharp_destroy_geom_info(geom_info%handle)
173+
geom_info%handle = c_null_ptr
174+
end subroutine sharp_destroy_geom_info
175+
176+
177+
! Currently the only mode supported is stacked (not interleaved) maps.
178+
!
179+
! Note that passing the exact dimension of alm/map is necesarry, it
180+
! prevents the caller from doing too crazy slicing prior to pass array
181+
! in...
182+
!
183+
! Usage:
184+
!
185+
! The alm array must have shape exactly alm(alm_info%n_local, nmaps)
186+
! The maps array must have shape exactly map(map_info%n_local, nmaps).
187+
subroutine sharp_execute_d(type, spin, nmaps, alm, alm_info, map, geom_info, &
188+
add, time, opcnt, comm)
189+
use iso_c_binding
190+
use mpi
191+
implicit none
192+
integer(c_int), value :: type, spin, nmaps
193+
integer(c_int), optional :: comm
194+
logical, value, optional :: add ! should add instead of replace out
195+
196+
type(sharp_alm_info) :: alm_info
197+
type(sharp_geom_info) :: geom_info
198+
real(c_double), intent(out), optional :: time
199+
integer(c_long_long), intent(out), optional :: opcnt
200+
real(c_double), target, intent(inout) :: alm(0:alm_info%n_local - 1, 1:nmaps)
201+
real(c_double), target, intent(inout) :: map(0:geom_info%n_local - 1, 1:nmaps)
202+
!--
203+
integer(c_int) :: mod_flags, ntrans, k
204+
type(c_ptr), target :: alm_ptr(nmaps)
205+
type(c_ptr), target :: map_ptr(nmaps)
206+
207+
mod_flags = SHARP_DP
208+
if (present(add) .and. add) then
209+
mod_flags = or(mod_flags, SHARP_ADD)
210+
end if
211+
212+
if (spin == 0) then
213+
ntrans = nmaps
214+
else
215+
ntrans = nmaps / 2
216+
end if
217+
218+
! Set up pointer table to access maps
219+
do k = 1, nmaps
220+
alm_ptr(k) = c_loc(alm(0, k))
221+
map_ptr(k) = c_loc(map(0, k))
222+
end do
223+
224+
if (present(comm)) then
225+
call c_sharp_execute_mpi(comm, type, spin, alm_ptr, map_ptr, &
226+
geom_info=geom_info%handle, &
227+
alm_info=alm_info%handle, &
228+
ntrans=ntrans, &
229+
flags=mod_flags, &
230+
time=time, &
231+
opcnt=opcnt)
232+
else
233+
call c_sharp_execute(type, spin, alm_ptr, map_ptr, &
234+
geom_info=geom_info%handle, &
235+
alm_info=alm_info%handle, &
236+
ntrans=ntrans, &
237+
flags=mod_flags, &
238+
time=time, &
239+
opcnt=opcnt)
240+
end if
241+
end subroutine sharp_execute_d
242+
243+
244+
245+
end module

‎fortran/test_sharp.f90‎

Lines changed: 57 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,57 @@
1+
program test_sharp
2+
use mpi
3+
use sharp
4+
use iso_c_binding, only : c_ptr, c_double
5+
implicit none
6+
7+
integer, parameter :: lmax = 2, nside = 2
8+
type(sharp_alm_info) :: alm_info
9+
type(sharp_geom_info) :: geom_info
10+
11+
real(c_double), dimension(0:(lmax + 1)**2 - 1, 1:1) :: alm
12+
real(c_double), dimension(0:12*nside**2 - 1, 1:1) :: map
13+
14+
integer(c_int), dimension(1:lmax + 1) :: ms
15+
integer(c_int), dimension(1:4 * nside - 1) :: rings
16+
integer(c_int) :: nm, m, nrings, iring
17+
integer :: nodecount, rank, ierr
18+
19+
call MPI_Init(ierr)
20+
call MPI_Comm_size(MPI_COMM_WORLD, nodecount, ierr)
21+
call MPI_Comm_rank(MPI_COMM_WORLD, rank, ierr)
22+
23+
nm = 0
24+
do m = rank, lmax, nodecount
25+
nm = nm + 1
26+
ms(nm) = m
27+
end do
28+
29+
nrings = 0
30+
do iring = rank + 1, 4 * nside - 1, nodecount
31+
nrings = nrings + 1
32+
rings(nrings) = iring
33+
end do
34+
35+
alm = 0
36+
map = 0
37+
if (rank == 0) then
38+
alm(0, 1) = 1
39+
end if
40+
41+
print *, ms(1:nm)
42+
call sharp_make_mmajor_real_packed_alm_info(lmax, ms=ms(1:nm), alm_info=alm_info)
43+
print *, 'alm_info%n_local', alm_info%n_local
44+
call sharp_make_healpix_geom_info(nside, rings=rings(1:nrings), geom_info=geom_info)
45+
print *, 'geom_info%n_local', geom_info%n_local
46+
print *, 'execute'
47+
call sharp_execute(SHARP_Y, 0, 1, alm, alm_info, map, geom_info, comm=MPI_COMM_WORLD)
48+
49+
print *, alm
50+
print *, map
51+
52+
call sharp_destroy_alm_info(alm_info)
53+
call sharp_destroy_geom_info(geom_info)
54+
print *, 'DONE'
55+
call MPI_Finalize(ierr)
56+
57+
end program test_sharp

‎libsharp/sharp_mpi.c‎

Lines changed: 15 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -326,4 +326,19 @@ void sharp_execute_mpi (MPI_Comm comm, sharp_jobtype type, int spin,
326326
if (opcnt!=NULL) *opcnt = job.opcnt;
327327
}
328328

329+
/* We declare this only in C file to make symbol available for Fortran wrappers;
330+
without declaring it in C header as it should not be available to C code */
331+
void sharp_execute_mpi_fortran(MPI_Fint comm, sharp_jobtype type, int spin,
332+
void *alm, void *map, const sharp_geom_info *geom_info,
333+
const sharp_alm_info *alm_info, int ntrans, int flags, double *time,
334+
unsigned long long *opcnt);
335+
void sharp_execute_mpi_fortran(MPI_Fint comm, sharp_jobtype type, int spin,
336+
void *alm, void *map, const sharp_geom_info *geom_info,
337+
const sharp_alm_info *alm_info, int ntrans, int flags, double *time,
338+
unsigned long long *opcnt)
339+
{
340+
sharp_execute_mpi(MPI_Comm_f2c(comm), type, spin, alm, map, geom_info,
341+
alm_info, ntrans, flags, time, opcnt);
342+
}
343+
329344
#endif

0 commit comments

Comments
 (0)