Skip to content
Permalink
Browse files

Add C and Fortran interfaces.

  • Loading branch information...
dmalhotra committed May 8, 2018
1 parent 7ee4066 commit 9b8393667d3f7a4f63f6f15c2f4509c9d6ec6a75
Showing with 910 additions and 1,169 deletions.
  1. +1 −2 Makefile.am
  2. +15 −15 examples/Makefile
  3. +97 −0 examples/src/example-c.c
  4. +101 −0 examples/src/example-f.f90
  5. +0 −52 include/fmm_gll.hpp
  6. +87 −0 include/pvfmm.f90
  7. +50 −0 include/pvfmm.h
  8. +1 −1 include/pvfmm_common.hpp
  9. +0 −1,099 src/fmm_gll.cpp
  10. +558 −0 src/pvfmm-wrapper.cpp
@@ -52,7 +52,6 @@ lib_libfmm_a_HEADERS = \
include/dtypes.h \
include/fft_wrapper.hpp \
include/fmm_cheb.hpp \
include/fmm_gll.hpp \
include/fmm_node.hpp \
include/fmm_pts.hpp \
include/fmm_pts_gpu.hpp \
@@ -106,13 +105,13 @@ lib_libpvfmm_a_SOURCES = \
$(lib_libpvfmm_a_HEADERS) \
src/cheb_utils.cpp \
src/device_wrapper.cpp \
src/fmm_gll.cpp \
src/legendre_rule.cpp \
src/mat_utils.cpp \
src/math_utils.cpp \
src/mem_mgr.cpp \
src/mortonid.cpp \
src/profile.cpp \
src/pvfmm-wrapper.cpp \
src/tree_node.cpp

if NVCC_OK
@@ -8,6 +8,8 @@ ifndef CXXFLAGS_PVFMM
$(error Cannot find file: MakeVariables)
endif

# FC=$(FC_PVFMM) # TODO: for now, FC must be provided by user
# CC=$(CC_PVFMM) # TODO: for now, CC must be provided by user
CXX=$(CXX_PVFMM)
CXXFLAGS=$(CXXFLAGS_PVFMM)
LDLIBS=$(LDLIBS_PVFMM)
@@ -28,29 +30,27 @@ TARGET_BIN = \

all : $(TARGET_BIN)

ifeq ($(INTEL_OFFLOAD_OK),yes)

$(BINDIR)/%: $(OBJDIR)/%.o
$(BINDIR)/%: $(SRCDIR)/%.f90
-@$(MKDIRS) $(dir $@)
$(CXX) $(CXXFLAGS) -no-offload $^ $(LDLIBS) -o $@
$(CXX) $(CXXFLAGS) $^_async $(LDLIBS) -o $@_async
$(CXX) $(CXXFLAGS) -D__DEVICE_SYNC__=1 $^_mic $(LDLIBS) -o $@_mic
$(FC) $(CXXFLAGS) -I$(INCDIR) $^ $(LDLIBS) -o $@

$(OBJDIR)/%.o: $(SRCDIR)/%.cpp
$(BINDIR)/%: $(SRCDIR)/%.c
-@$(MKDIRS) $(dir $@)
$(CXX) $(CXXFLAGS) -no-offload -I$(INCDIR) -c $^ -o $@
$(CXX) $(CXXFLAGS) -I$(INCDIR) -c $^ -o $@_async
$(CXX) $(CXXFLAGS) -D__DEVICE_SYNC__=1 -I$(INCDIR) -c $^ -o $@_mic
$(CC) $(CXXFLAGS) -I$(INCDIR) $^ $(LDLIBS) -o $@

else
ifeq ($(INTEL_OFFLOAD_OK),yes)

$(BINDIR)/%: $(OBJDIR)/%.o
$(BINDIR)/%: $(SRCDIR)/%.cpp
-@$(MKDIRS) $(dir $@)
$(CXX) $(CXXFLAGS) $^ $(LDLIBS) -o $@
$(CXX) $(CXXFLAGS) -no-offload -I$(INCDIR) $^ $(LDLIBS) -o $@
$(CXX) $(CXXFLAGS) -I$(INCDIR) $^ $(LDLIBS) -o $@_async
$(CXX) $(CXXFLAGS) -D__DEVICE_SYNC__=1 -I$(INCDIR) $^ $(LDLIBS) -o $@_mic

else

$(OBJDIR)/%.o: $(SRCDIR)/%.cpp
$(BINDIR)/%: $(SRCDIR)/%.cpp
-@$(MKDIRS) $(dir $@)
$(CXX) $(CXXFLAGS) -I$(INCDIR) -c $^ -o $@
$(CXX) $(CXXFLAGS) -I$(INCDIR) $^ $(LDLIBS) -o $@

endif

@@ -0,0 +1,97 @@
#include <omp.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <pvfmm.h>

void BiotSavart(const double* src_X, const double* src_V, long Ns, const double* trg_X, double* trg_V, long Nt) { // Direct evaluation
double oofp = 1 / (4 * M_PI);
for (long t = 0; t < Nt; t++) {
trg_V[t*3+0] = 0;
trg_V[t*3+1] = 0;
trg_V[t*3+2] = 0;
for (long s = 0; s < Ns; s++) {
double X[3];
X[0] = trg_X[t*3+0] - src_X[s*3+0];
X[1] = trg_X[t*3+1] - src_X[s*3+1];
X[2] = trg_X[t*3+2] - src_X[s*3+2];
double r2 = X[0]*X[0] + X[1]*X[1] + X[2]*X[2];
double rinv = (r2 > 0 ? 1/sqrt(r2) : 0);
double rinv3 = rinv * rinv * rinv;

trg_V[t*3+0] += (src_V[s*3+1]*X[2] - src_V[s*3+2]*X[1]) * rinv3 * oofp;
trg_V[t*3+1] += (src_V[s*3+2]*X[0] - src_V[s*3+0]*X[2]) * rinv3 * oofp;
trg_V[t*3+2] += (src_V[s*3+0]*X[1] - src_V[s*3+1]*X[0]) * rinv3 * oofp;
}
}
}

void test_FMM(void* ctx) { // Compare FMM and direct evaluation results
int kdim[2] = {3,3};
int64_t Ns = 20000;
int64_t Nt = 20000;

// Initialize target coordinates
double* trg_X = (double*)malloc(Nt * 3 * sizeof(double));
double* trg_V = (double*)malloc(Nt * kdim[1] * sizeof(double));
double* trg_V0 = (double*)malloc(Nt * kdim[1] * sizeof(double));
for (long i = 0; i < Nt * 3; i++) trg_X[i] = drand48();

// Initialize source coordinates and density
double* src_X = (double*)malloc(Ns * 3 * sizeof(double));
double* src_V = (double*)malloc(Ns * kdim[0] * sizeof(double));
for (long i = 0; i < Ns * 3; i++) src_X[i] = drand48();
for (long i = 0; i < Ns * kdim[0]; i++) src_V[i] = drand48() - 0.5;

// FMM evaluation
double tt;
int32_t setup = 1;
tt = -omp_get_wtime();
PVFMMEvalD(src_X, src_V, NULL, Ns, trg_X, trg_V, Nt, ctx, setup);
printf("FMM evaluation time (with setup) : %f\n", tt + omp_get_wtime());

// FMM evaluation (without setup)
setup = 0;
tt = -omp_get_wtime();
PVFMMEvalD(src_X, src_V, NULL, Ns, trg_X, trg_V, Nt, ctx, setup);
printf("FMM evaluation time (without setup) : %f\n", tt + omp_get_wtime());

// Direct evaluation
tt = -omp_get_wtime();
BiotSavart(&src_X[0], &src_V[0], Ns, &trg_X[0], &trg_V0[0], Nt);
printf("Direct evaluation time : %f\n", tt + omp_get_wtime());

double max_err = 0;
for (long i = 0; i < Nt * kdim[1]; i++) { // Compute error
double err = fabs(trg_V[i] - trg_V0[i]);
if (err > max_err) max_err = err;
}
printf("Max-error : %f\n", max_err);

free(src_X);
free(src_V);
free(trg_X);
free(trg_V);
free(trg_V0);
}

int main(int argc, char** argv) {
MPI_Init(&argc, &argv);

void* ctx;
{ // Create FMM context
double box_size = -1;
int32_t points_per_box = 1000;
int32_t multipole_order = 10;
int32_t kernel = PVFMMBiotSavartPotential;
ctx = PVFMMCreateContextD(box_size, points_per_box, multipole_order, kernel, MPI_COMM_WORLD);
}

test_FMM(ctx);

PVFMMDestroyContextD(&ctx);

MPI_Finalize();
return 0;
}

@@ -0,0 +1,101 @@
program main
use iso_c_binding
implicit none
include 'mpif.h'
include 'pvfmm.f90'
integer ierror
integer :: comm
real*8 :: box_size
integer*4 :: points_per_leaf, multipole_order, kernel
type (c_ptr) :: fmm_ctx
integer*8 :: Ns, Nt

call MPI_Init(ierror)

! Create FMM context
box_size = -1.0 ! for periodic boundary conditions
points_per_leaf = 1000 ! tuning parameter
multipole_order = 10 ! accuracy
kernel = PVFMMBiotSavartPotential ! kernel function
comm = MPI_COMM_WORLD
call PVFMMCreateContextD(fmm_ctx, box_size, points_per_leaf, multipole_order, kernel, comm)

! Evaluate FMM
Ns = 20000
Nt = 20000
call test(fmm_ctx, Ns, Nt)

! Destroy FMM context
call PVFMMDestroyContextD(fmm_ctx)

call MPI_Finalize(ierror)
end

subroutine BiotSavart(Xs, Vs, Ns, Xt, Vt, Nt)
implicit none
integer*8 Ns, Nt, s, t
real*8 :: Xs(Ns * 3), Vs(Ns * 3), Xt(Nt * 3), Vt(Nt * 3), Vt_ref(Nt * 3)
real*8 :: oofp, X(3), rinv, rinv3

oofp = 1/(16*atan(1.0))
do t = 0, Nt-1
Vt(t*3+1) = 0
Vt(t*3+2) = 0
Vt(t*3+3) = 0
do s = 0, Ns-1
X(1) = Xt(t*3+1) - Xs(s*3+1)
X(2) = Xt(t*3+2) - Xs(s*3+2)
X(3) = Xt(t*3+3) - Xs(s*3+3)

rinv = X(1)*X(1) + X(2)*X(2) + X(3)*X(3)
if (rinv .gt. 0) then
rinv = 1/sqrt(rinv)
endif
rinv3 = rinv*rinv*rinv

Vt(t*3+1) = Vt(t*3+1) + (Vs(s*3+2)*X(3) - Vs(s*3+3)*X(2)) * rinv3 * oofp;
Vt(t*3+2) = Vt(t*3+2) + (Vs(s*3+3)*X(1) - Vs(s*3+1)*X(3)) * rinv3 * oofp;
Vt(t*3+3) = Vt(t*3+3) + (Vs(s*3+1)*X(2) - Vs(s*3+2)*X(1)) * rinv3 * oofp;
enddo
enddo
endsubroutine

subroutine test(fmm_ctx, Ns, Nt)
use iso_c_binding
implicit none
include 'pvfmm.f90'
integer*8 :: i, Ns, Nt
real*8 :: Xs(Ns * 3), Vs(Ns * 3), Xt(Nt * 3), Vt(Nt * 3), Vt_ref(Nt * 3)
double precision :: omp_get_wtime, tt
type (c_ptr) :: fmm_ctx
integer*4 :: setup

call srand(0)
do i=1, Ns*3
Xs(i) = rand()
Vs(i) = rand()
enddo
do i=1, Nt*3
Xt(i) = rand()
enddo

setup = 1
tt = -omp_get_wtime()
call PVFMMEvalD(Xs, Vs, Ns, Xt, Vt, Nt, fmm_ctx, setup)
tt = tt + omp_get_wtime()
print*, "FMM evaluation time (with setup) : ", tt

setup = 0
tt = -omp_get_wtime()
call PVFMMEvalD(Xs, Vs, Ns, Xt, Vt, Nt, fmm_ctx, setup)
tt = tt + omp_get_wtime()
print*, "FMM evaluation time (without setup) : ", tt

tt = -omp_get_wtime()
call BiotSavart(Xs, Vs, Ns, Xt, Vt_ref, Nt)
tt = tt + omp_get_wtime()
print*, "Direct evaluation time : ", tt

print*, "Maximum relative error : ", maxval(abs(Vt_ref - Vt)) / maxval(abs(Vt_ref));
endsubroutine

This file was deleted.

Oops, something went wrong.
Oops, something went wrong.

0 comments on commit 9b83936

Please sign in to comment.
You can’t perform that action at this time.
You signed in with another tab or window. Reload to refresh your session. You signed out in another tab or window. Reload to refresh your session.