Adding fast multipole method for Legendre-Chebyshev transforms

.github/workflows/ci.yml

@@ -28,7 +28,7 @@ jobs:
       CC: ${{ matrix.compiler }}
       LD_LIBRARY_PATH: $LD_LIBRARY_PATH:.
     steps:
-      - uses: actions/checkout@v3
+      - uses: actions/checkout@v4
       - name: Install dependencies on ubuntu gcc
         if: startsWith(matrix.os, 'ubuntu') && startsWith(matrix.compiler, 'gcc')
         run: |
@@ -65,13 +65,13 @@ jobs:
     name: Documentation
     runs-on: macOS-latest
     steps:
-      - uses: actions/checkout@v3
+      - uses: actions/checkout@v4
       - name: Install Doxygen
         run: brew install doxygen
       - name: Build site
         run: doxygen Doxyfile
       - if: ${{ github.ref == 'refs/heads/master' }}
-        uses: peaceiris/actions-gh-pages@v3
+        uses: peaceiris/actions-gh-pages@v4
         with:
           github_token: ${{ secrets.GITHUB_TOKEN }}
           publish_dir: docs/html

Make.inc

@@ -84,7 +84,7 @@ ASM = src/recurrence/recurrence_default.s \
       src/rotations/rotations_AVX512F.s \
       src/rotations/rotations_NEON.s
 
-SRC = src/recurrence.c src/transforms.c src/rotations.c src/permute.c src/tdc.c src/drivers.c src/fftw.c src/isometries.c src/helmholtzhodge.c
+SRC = src/fmm.c src/recurrence.c src/transforms.c src/rotations.c src/permute.c src/tdc.c src/drivers.c src/fftw.c src/isometries.c src/helmholtzhodge.c
 
 machine := $(shell $(CC) -dumpmachine | cut -d'-' -f1)
 
@@ -120,23 +120,28 @@ ifdef FT_PREFIX
     endif
 else
     ifeq ($(UNAME), Darwin)
-        CFLAGS += -I/Library/Developer/CommandLineTools/SDKs/MacOSX.sdk/usr/include
-        LDFLAGS += -L/Library/Developer/CommandLineTools/SDKs/MacOSX.sdk/usr/lib
-        ifeq ($(FT_USE_APPLEBLAS), 1)
-            CFLAGS += -I/Library/Developer/CommandLineTools/SDKs/MacOSX.sdk/System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/Headers
-            LDFLAGS += -L/Library/Developer/CommandLineTools/SDKs/MacOSX.sdk/System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A
+        ifeq (,$(shell which conda))
+            CFLAGS += -I/Library/Developer/CommandLineTools/SDKs/MacOSX.sdk/usr/include
+            LDFLAGS += -L/Library/Developer/CommandLineTools/SDKs/MacOSX.sdk/usr/lib
+            ifeq ($(FT_USE_APPLEBLAS), 1)
+                CFLAGS += -I/Library/Developer/CommandLineTools/SDKs/MacOSX.sdk/System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/Headers
+                LDFLAGS += -L/Library/Developer/CommandLineTools/SDKs/MacOSX.sdk/System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A
+            else
+                CFLAGS += -I/opt/homebrew/opt/openblas/include
+                LDFLAGS += -L/opt/homebrew/opt/openblas/lib
+            endif
+            CFLAGS += -I/opt/homebrew/opt/libomp/include
+            LDFLAGS += -L/opt/homebrew/opt/libomp/lib
+            CFLAGS += -I/opt/homebrew/opt/fftw/include
+            LDFLAGS += -L/opt/homebrew/opt/fftw/lib
+            CFLAGS += -I/opt/homebrew/opt/mpfr/include
+            LDFLAGS += -L/opt/homebrew/opt/mpfr/lib
+            CFLAGS += -I/opt/homebrew/opt/gmp/include
+            LDFLAGS += -L/opt/homebrew/opt/gmp/lib
         else
-            CFLAGS += -I/usr/local/opt/openblas/include
-            LDFLAGS += -L/usr/local/opt/openblas/lib
+            CFLAGS += -I$(CONDA_PREFIX)/include
+            LDFLAGS += -L$(CONDA_PREFIX)/lib
         endif
-        CFLAGS += -I/usr/local/opt/libomp/include
-        LDFLAGS += -L/usr/local/opt/libomp/lib
-        CFLAGS += -I/usr/local/opt/fftw/include
-        LDFLAGS += -L/usr/local/opt/fftw/lib
-        CFLAGS += -I/usr/local/opt/mpfr/include
-        LDFLAGS += -L/usr/local/opt/mpfr/lib
-        CFLAGS += -I/usr/local/opt/gmp/include
-        LDFLAGS += -L/usr/local/opt/gmp/lib
     else ifeq ($(UNAME), Windows)
         CFLAGS += -Ic:\tools\vcpkg\installed\x64-windows\include\openblas
         CFLAGS += -Ic:\tools\vcpkg\installed\x64-windows\include
@@ -164,7 +169,7 @@ ifndef FT_OPENMP
     endif
 endif
 
-LDLIBS += $(FT_OPENMP) -lm -lmpfr -l$(FT_BLAS) -lfftw3
+LDLIBS += $(FT_OPENMP) -lm -lmpfr -l$(FT_BLAS) -lfftw3 -lfftw3f
 
 ifdef FT_QUADMATH
     AFLAGS += -DFT_QUADMATH
@@ -174,5 +179,5 @@ ifdef FT_QUADMATH
 endif
 
 ifneq ($(FT_FFTW_WITH_COMBINED_THREADS), 1)
-    LDLIBS += -lfftw3_threads
+    LDLIBS += -lfftw3_threads -lfftw3f_threads
 endif

Makefile

@@ -33,6 +33,7 @@ lib:
 	$(CC) $(CFLAGS) $(LIBFLAGS) $(ASM) $(SRC) $(LDFLAGS) $(LDLIBS) -o lib$(LIB).$(SLIB)
 
 tests:
+	$(CC) src/ftutilities.c test/test_fmm.c $(CFLAGS) -L$(LIBDIR) -l$(LIB) $(LDFLAGS) $(LDLIBS) -o test_fmm
 	$(CC) src/ftutilities.c test/test_recurrence.c $(CFLAGS) -L$(LIBDIR) -l$(LIB) $(LDFLAGS) $(LDLIBS) -o test_recurrence
 	$(CC) src/ftutilities.c test/test_transforms.c $(CFLAGS) -L$(LIBDIR) -l$(LIB) $(LDFLAGS) $(LDLIBS) -o test_transforms
 	$(CC) src/ftutilities.c test/test_permute.c $(CFLAGS) -L$(LIBDIR) -l$(LIB) $(LDFLAGS) $(LDLIBS) -o test_permute
@@ -58,6 +59,7 @@ examples:
 	$(CC) src/ftutilities.c examples/subspaceangles.c $(CFLAGS) -L$(LIBDIR) -l$(LIB) $(LDFLAGS) $(LDLIBS) -o subspaceangles
 
 runtests:
+	./test_fmm
 	./test_recurrence
 	./test_transforms
 	./test_permute

src/fmm.c

@@ -0,0 +1,67 @@
+#include "fasttransforms.h"
+#include "ftinternal.h"
+#include "fmm.h"
+
+size_t get_number_of_blocks(const size_t level) { return pow(2, level + 1) - 1; }
+
+size_t get_h(const size_t level, const size_t L) { return pow(2, L - level - 1); }
+
+size_t get_number_of_submatrices(const size_t level) { return 3 * get_number_of_blocks(level); }
+
+size_t get_total_number_of_blocks(const size_t L) { return pow(2, L + 1) - (L + 2); }
+
+size_t get_total_number_of_submatrices(const size_t L) { return 3 * get_total_number_of_blocks(L); }
+
+void get_ij(size_t *ij, const size_t level, const size_t block, const size_t s, const size_t L)
+{
+  size_t h = get_h(level, L);
+  ij[0] = 2 * block * s * h;
+  ij[1] = ij[0] + 2 * s * h;
+}
+
+/*#if defined(FT_QUADMATH)
+    #define FLT quadruple
+    #define X(name) FT_CONCAT(ft_, name, q)
+    #define Y(name) FT_CONCAT(, name, q)
+    #include "fmm_source.c"
+    #undef FLT
+    #undef X
+    #undef Y
+
+    #define FLT long double
+    #define X(name) FT_CONCAT(ft_, name, l)
+    #define X2(name) FT_CONCAT(ft_, name, q)
+    #define Y(name) FT_CONCAT(, name, l)
+    #include "fmm_source.c"
+    #undef FLT
+    #undef X
+    #undef X2
+    #undef Y
+#else
+    #define FLT long double
+    #define X(name) FT_CONCAT(ft_, name, l)
+    #define X2(name) FT_CONCAT(ft_, name, l)
+    #define Y(name) FT_CONCAT(, name, l)
+    #include "fmm_source.c"
+    #undef FLT
+    #undef X
+    #undef X2
+    #undef Y
+#endif*/
+
+
+#define FLT double
+#define X(name) FT_CONCAT(ft_, name, )
+#define FT_USE_DOUBLE
+#include "fmm_source.c"
+#undef FLT
+#undef X
+#undef FT_USE_DOUBLE
+
+#define FLT float
+#define X(name) FT_CONCAT(ft_, name, f)
+#define FT_USE_SINGLE
+#include "fmm_source.c"
+#undef FLT
+#undef X
+#undef FT_USE_SINGLE

src/fmm.h

@@ -0,0 +1,22 @@
+#ifndef FMM_H
+#define FMM_H
+
+#include <assert.h>
+#include <stdbool.h>
+#include <string.h>
+#include <stdint.h>
+#include <time.h>
+
+enum { L2C = 0, C2L = 1, BOTH = 2 };
+
+#ifdef CLOCK_UPTIME_RAW
+#define tic clock_gettime_nsec_np(CLOCK_UPTIME_RAW)
+#define dtics(a, b) (double)(b - a) / 1.0E9
+#define toc(a) (double)(tic - a) / 1.0E9
+#else
+#define tic clock()
+#define dtics(a, b) (double)(b - a) / (double)CLOCKS_PER_SEC
+#define toc(a) (double)(tic - a) / (double)CLOCKS_PER_SEC
+#endif
+
+#endif