add precompute for m2l

3dpp/include/kernel.h

@@ -3,6 +3,8 @@
 #include "body.h"
 #include "tree.h"
 #include <map>
+#include "traverser.h"
+#include <fftw3.h>
 
 namespace rtfmm
 {
@@ -32,6 +34,12 @@ class LaplaceKernel
 
     void m2l_fft(int P, Cell3& cell_src, Cell3& cell_tar, vec3r offset = vec3r(0,0,0));
 
+    void m2l_fft_precompute_naive(int P, Cells3& cs, PeriodicInteractionMap& m2l_map, PeriodicInteractionPairs& m2l_pairs);
+
+    void m2l_fft_precompute_advanced(int P, Cells3& cs, PeriodicInteractionMap& m2l_map, PeriodicInteractionPairs& m2l_pairs);
+
+    void m2l_fft_precompute_advanced2(int P, Cells3& cs, PeriodicInteractionMap& m2l_map, PeriodicInteractionPairs& m2l_pairs);
+
     void l2l(int P, Cell3& cell_parent, Cells3& cs);
 
     void l2l_precompute(int P, Cell3& cell_parent, Cells3& cs);
@@ -48,16 +56,18 @@ class LaplaceKernel
 
     void precompute(int P, real r0);
 
+    void precompute_m2l(int P, real r0, Cells3 cs, PeriodicInteractionMap m2l_map, int images);
+
 private:
     Matrix get_p2p_matrix(
         std::vector<vec3r>& x_src, 
         std::vector<vec3r>& x_tar
     );
 
     Matriv get_force_naive(
-        std::vector<vec3r> x_src, 
-        std::vector<vec3r> x_tar, 
-        Matrix q_src
+        std::vector<vec3r>& x_src, 
+        std::vector<vec3r>& x_tar, 
+        Matrix& q_src
     );
 
     std::vector<rtfmm::real> get_G_matrix(std::vector<rtfmm::vec3r>& grid, int N);
@@ -68,6 +78,7 @@ class LaplaceKernel
     Matrix UT_p2m_precompute;
     Matrix V_p2m_precompute;
     Matrix Sinv_p2m_precompute;
+    Matrix VSinv_p2m_precompute;
 
     Matrix UT_l2p_precompute;
     Matrix V_l2p_precompute;
@@ -77,5 +88,9 @@ class LaplaceKernel
     Matrix matrix_m2m_img[27];
     Matrix matrix_l2l[8];
     Matrix matrix_l2l_img;
+
+    std::vector<int> m2l_tars;
+    std::vector<int> m2l_srcs;
+    std::map<int, std::vector<complexr>> m2l_Gks;
 };
 };

3dpp/include/mathfunc.h

@@ -15,13 +15,14 @@ namespace rtfmm
 void mat_scale(Matrix& A, real scale);
 
 /**
- * @brief c = A x b
+ * @brief c = k x A x b
  * 
  * @param A m by n matrix
  * @param b n by 1 vector
+ * @param k real number
  * @return m by 1 vector
  */
-Matrix mat_vec_mul(Matrix A, Matrix b); 
+Matrix mat_vec_mul(const Matrix& A, const Matrix& b, real k = 1.0); 
 
 /**
  * @brief C = A + B
@@ -39,7 +40,7 @@ Matrix mat_mat_add(Matrix A, Matrix B);
  * @param B n by k matrix
  * @return m by k matrix
  */
-Matrix mat_mat_mul(Matrix A, Matrix B); 
+Matrix mat_mat_mul(const Matrix& A, const Matrix& B); 
 
 /**
  * @brief solve A = U x S x VT

3dpp/include/type.h

@@ -25,6 +25,12 @@ using Indices = std::vector<int>;
 
 using real = double;
 
+struct complexr
+{
+    real r;
+    real i;
+};
+
 /**
  * @brief alert error message and exit if conditon is false
 */

3dpp/src/fmm.cpp

@@ -30,6 +30,7 @@ rtfmm::Bodies3 rtfmm::LaplaceFMM::solve()
     {
         TIME_BEGIN(precompute);
         kernel.precompute(args.P, args.r);
+        kernel.precompute_m2l(args.P, args.r, cs, traverser.get_map(OperatorType::M2L), args.images);
         if(args.timing) {TIME_END(precompute);}
     }
 
@@ -134,18 +135,33 @@ void rtfmm::LaplaceFMM::M2L()
 {
     TIME_BEGIN(M2L);
     PeriodicInteractionMap m2l_map = traverser.get_map(OperatorType::M2L);
-    //#pragma omp parallel for
-    for(int i = 0; i < m2l_map.size(); i++)
+    PeriodicInteractionPairs m2l_pairs = traverser.get_pairs(OperatorType::M2L);
+    if(args.use_precompute)
+    {
+        TIME_BEGIN(M2L_kernel);
+        kernel.m2l_fft_precompute_advanced2(args.P, cs, m2l_map, m2l_pairs);
+        TIME_END(M2L_kernel);
+    }
+    else
     {
-        auto m2l_list = m2l_map.begin();
-        std::advance(m2l_list, i);
-        Cell3& ctar = cs[m2l_list->first];
-        for(int j = 0; j < m2l_list->second.size(); j++)
+        for(int i = 0; i < m2l_map.size(); i++)
         {
-            Cell3& csrc = cs[m2l_list->second[j].first];
-            vec3r offset_src = m2l_list->second[j].second;
-            if(args.use_fft) kernel.m2l_fft(args.P, csrc, ctar, offset_src);
-            else kernel.m2l(args.P, csrc, ctar, offset_src);
+            auto m2l_list = m2l_map.begin();
+            std::advance(m2l_list, i);
+            Cell3& ctar = cs[m2l_list->first];
+            for(int j = 0; j < m2l_list->second.size(); j++)
+            {
+                Cell3& csrc = cs[m2l_list->second[j].first];
+                vec3r offset_src = m2l_list->second[j].second;
+                if(args.use_fft) 
+                {
+                    kernel.m2l_fft(args.P, csrc, ctar, offset_src);
+                }
+                else
+                {
+                    kernel.m2l(args.P, csrc, ctar, offset_src);
+                }
+            }
         }
     }
     if(args.timing)