Merge pull request #4 from dagss/legtrans

Legendre transforms

Author: Dag Sverre Seljebotn

.gitignore

@@ -1,6 +1,9 @@
 *.o
+*.so
 #*
 *~
+*.pyc
+*.pyo
 
 /auto
 /autom4te.cache
@@ -9,3 +12,5 @@
 /config/config.auto
 /configure
 /sharp_oracle.inc
+
+/python/libsharp/libsharp.c

Makefile

@@ -53,3 +53,18 @@ perftest: compile_all
 	$(BINDIR)/sharp_testsuite test gauss 2047 -1 -1 4096 0 1 && \
 	$(BINDIR)/sharp_testsuite test gauss 4095 -1 -1 8192 0 1 && \
 	$(BINDIR)/sharp_testsuite test gauss 8191 -1 -1 16384 0 1
+
+%.c: %.c.in
+# Only do this if the md5sum changed, in order to avoid Python and Jinja
+# dependency when not modifying the c.in file
+	grep `md5sum $< | cut -d ' ' -f 1` $@ || ./runjinja.py < $< > $@
+
+genclean:
+	rm libsharp/sharp_legendre.c || exit 0
+
+pytest:
+	rm python/libsharp/libsharp.so || exit 0
+	cd python && LIBSHARP_INCLUDE=$(INCDIR) LIBSHARP_LIB=$(LIBDIR) python setup.py build_ext --inplace
+	cd python && nosetests libsharp
+
+

fortran/sharp.f90

@@ -103,12 +103,32 @@ subroutine c_sharp_execute_mpi(comm, type, spin, alm, map, geom_info, alm_info,
        type(c_ptr), intent(in)                      :: alm(*), map(*)
      end subroutine c_sharp_execute_mpi
 
+     ! Legendre transforms
+     subroutine c_sharp_legendre_transform(bl, recfac, lmax, x, out, nx) &
+          bind(c, name='sharp_legendre_transform')
+       use iso_c_binding
+       integer(c_ptrdiff_t), value :: lmax, nx
+       real(c_double) :: bl(lmax + 1), x(nx), out(nx)
+       real(c_double), optional :: recfac(lmax + 1)
+     end subroutine c_sharp_legendre_transform
+
+     subroutine c_sharp_legendre_transform_s(bl, recfac, lmax, x, out, nx) &
+          bind(c, name='sharp_legendre_transform_s')
+       use iso_c_binding
+       integer(c_ptrdiff_t), value :: lmax, nx
+       real(c_float) :: bl(lmax + 1), x(nx), out(nx)
+       real(c_float), optional :: recfac(lmax + 1)
+     end subroutine c_sharp_legendre_transform_s
   end interface
 
   interface sharp_execute
      module procedure sharp_execute_d
   end interface
 
+  interface sharp_legendre_transform
+     module procedure sharp_legendre_transform_d, sharp_legendre_transform_s
+  end interface sharp_legendre_transform
+
 contains
   ! alm info
 
@@ -240,6 +260,25 @@ subroutine sharp_execute_d(type, spin, nmaps, alm, alm_info, map, geom_info, &
    end if
   end subroutine sharp_execute_d
 
+  subroutine sharp_legendre_transform_d(bl, x, out)
+    use iso_c_binding
+    real(c_double) :: bl(:)
+    real(c_double) :: x(:), out(size(x))
+    !--
+    integer(c_ptrdiff_t) :: lmax, nx
+    call c_sharp_legendre_transform(bl, lmax=int(size(bl) - 1, c_ptrdiff_t), &
+                                    x=x, out=out, nx=int(size(x), c_ptrdiff_t))
+  end subroutine sharp_legendre_transform_d
+
+  subroutine sharp_legendre_transform_s(bl, x, out)
+    use iso_c_binding
+    real(c_float) :: bl(:)
+    real(c_float) :: x(:), out(size(x))
+    !--
+    integer(c_ptrdiff_t) :: lmax, nx
+    call c_sharp_legendre_transform_s(bl, lmax=int(size(bl) - 1, c_ptrdiff_t), &
+                                      x=x, out=out, nx=int(size(x), c_ptrdiff_t))
+  end subroutine sharp_legendre_transform_s
 
 
 end module

fortran/test_sharp.f90

@@ -54,4 +54,31 @@ program test_sharp
   print *, 'DONE'
   call MPI_Finalize(ierr)
 
+  print *, 'LEGENDRE TRANSFORMS'
+
+  call test_legendre_transforms()
+
+contains
+  subroutine test_legendre_transforms()
+    integer, parameter :: lmax = 20, nx=10
+    real(c_double) :: bl(0:lmax)
+    real(c_double) :: x(nx), out(nx)
+    real(c_float) :: out_s(nx)
+    !--
+    integer :: l, i
+
+    do l = 0, lmax
+       bl(l) = 1.0 / real(l + 1, c_double)
+    end do
+    do i = 1, nx
+       x(i) = 1 / real(i, c_double)
+    end do
+    out = 0
+    call sharp_legendre_transform(bl, x, out)
+    print *, out
+    call sharp_legendre_transform(real(bl, c_float), real(x, c_float), out_s)
+    print *, out_s
+  end subroutine test_legendre_transforms
+
+
 end program test_sharp

libsharp/planck.make

@@ -8,7 +8,7 @@ FULL_INCLUDE+= -I$(SD)
 HDR_$(PKG):=$(SD)/*.h
 LIB_$(PKG):=$(LIBDIR)/libsharp.a
 BIN:=sharp_testsuite
-LIBOBJ:=sharp_ylmgen_c.o sharp.o sharp_announce.o sharp_geomhelpers.o sharp_almhelpers.o sharp_core.o
+LIBOBJ:=sharp_ylmgen_c.o sharp.o sharp_announce.o sharp_geomhelpers.o sharp_almhelpers.o sharp_core.o sharp_legendre.o
 ALLOBJ:=$(LIBOBJ) sharp_testsuite.o
 LIBOBJ:=$(LIBOBJ:%=$(OD)/%)
 ALLOBJ:=$(ALLOBJ:%=$(OD)/%)

libsharp/sharp.h

@@ -39,5 +39,6 @@
 #include <complex.h>
 
 #include "sharp_lowlevel.h"
+#include "sharp_legendre.h"
 
 #endif