Fix bug in periodic BC

The bug affected some scale-invariant kernels.

include/fmm_pts.txx

@@ -454,8 +454,8 @@ Matrix<typename FMMNode::Real_t>& FMM_Pts<FMMNode>::Precomp(int level, Mat_Type
       if(M0.Dim(0)==0 || M0.Dim(1)==0) return M_;
 
       for(size_t i=0;i<Perm_Count;i++) this->PrecompPerm(type, (Perm_Type) i);
-      Permutation<Real_t>& Pr = this->interac_list.Perm_R(level, type, mat_indx);
-      Permutation<Real_t>& Pc = this->interac_list.Perm_C(level, type, mat_indx);
+      Permutation<Real_t>& Pr = this->interac_list.Perm_R(abs(level), type, mat_indx);
+      Permutation<Real_t>& Pc = this->interac_list.Perm_C(abs(level), type, mat_indx);
       if(Pr.Dim()>0 && Pc.Dim()>0 && M0.Dim(0)>0 && M0.Dim(1)>0) return M_;
     }
   }
@@ -843,18 +843,24 @@ Matrix<typename FMMNode::Real_t>& FMM_Pts<FMMNode>::Precomp(int level, Mat_Type
         for(int level=0; level>=-BC_LEVELS; level--){
           { // Compute M_l2l
             this->Precomp(level, D2D_Type, 0);
-            Permutation<Real_t>& Pr = this->interac_list.Perm_R(level, D2D_Type, 0);
-            Permutation<Real_t>& Pc = this->interac_list.Perm_C(level, D2D_Type, 0);
+            Permutation<Real_t> Pr = this->interac_list.Perm_R(abs(level), D2D_Type, 0);
+            Permutation<Real_t> Pc = this->interac_list.Perm_C(abs(level), D2D_Type, 0);
+            { // Invert scaling because level<0
+              for(long i=0;i<Pr.Dim();i++) Pr.scal[i]=1.0/Pr.scal[i];
+              for(long i=0;i<Pc.Dim();i++) Pc.scal[i]=1.0/Pc.scal[i];
+            }
             M_l2l[-level] = M_check_zero_avg * Pr * this->Precomp(level, D2D_Type, this->interac_list.InteracClass(D2D_Type, 0)) * Pc * M_check_zero_avg;
             assert(M_l2l[-level].Dim(0)>0 && M_l2l[-level].Dim(1)>0);
           }
 
           // Compute M_m2m
           for(size_t mat_indx=0; mat_indx<mat_cnt_m2m; mat_indx++){
-            this->Precomp(level, U2U_Type, mat_indx);
-            Permutation<Real_t>& Pr = this->interac_list.Perm_R(level, U2U_Type, mat_indx);
-            Permutation<Real_t>& Pc = this->interac_list.Perm_C(level, U2U_Type, mat_indx);
-            Matrix<Real_t> M = Pr * this->Precomp(level, U2U_Type, this->interac_list.InteracClass(U2U_Type, mat_indx)) * Pc;
+            this->Precomp(level-1, U2U_Type, mat_indx);
+            Permutation<Real_t> Pr = this->interac_list.Perm_R(abs(level-1), U2U_Type, mat_indx);
+            Permutation<Real_t> Pc = this->interac_list.Perm_C(abs(level-1), U2U_Type, mat_indx);
+            for(long i=0;i<Pr.Dim();i++) Pr.scal[i]=1.0/Pr.scal[i];
+            for(long i=0;i<Pc.Dim();i++) Pc.scal[i]=1.0/Pc.scal[i];
+            Matrix<Real_t> M = Pr * this->Precomp(level-1, U2U_Type, this->interac_list.InteracClass(U2U_Type, mat_indx)) * Pc;
             assert(M.Dim(0)>0 && M.Dim(1)>0);
 
             if(mat_indx==0) M_m2m[-level] = M_equiv_zero_avg*M*M_equiv_zero_avg;
@@ -1113,8 +1119,8 @@ void FMM_Pts<FMMNode>::PrecompAll(Mat_Type type, int level){
     //#pragma omp parallel for //lets use fine grained parallelism
     for(size_t mat_indx=0;mat_indx<mat_cnt;mat_indx++){
       Matrix<Real_t>& M0=interac_list.ClassMat(level,(Mat_Type)type,mat_indx);
-      Permutation<Real_t>& pr=interac_list.Perm_R(level, (Mat_Type)type, mat_indx);
-      Permutation<Real_t>& pc=interac_list.Perm_C(level, (Mat_Type)type, mat_indx);
+      Permutation<Real_t>& pr=interac_list.Perm_R(abs(level), (Mat_Type)type, mat_indx);
+      Permutation<Real_t>& pc=interac_list.Perm_C(abs(level), (Mat_Type)type, mat_indx);
       if(pr.Dim()!=M0.Dim(0) || pc.Dim()!=M0.Dim(1)) Precomp(level, (Mat_Type)type, mat_indx);
     }
   }

include/interac_list.txx

@@ -285,6 +285,7 @@ Matrix<typename Node_t::Real_t>& InteracList<Node_t>::ClassMat(int l, Mat_Type t
 
 template <class Node_t>
 Permutation<typename Node_t::Real_t>& InteracList<Node_t>::Perm_R(int l, Mat_Type type, size_t indx){
+  assert(l>=0);
   size_t indx0=InteracClass(type, indx);
   Matrix     <Real_t>& M0      =mat->Mat   (l, type, indx0);
   Permutation<Real_t>& row_perm=mat->Perm_R(l, type, indx );
@@ -313,6 +314,7 @@ Permutation<typename Node_t::Real_t>& InteracList<Node_t>::Perm_R(int l, Mat_Typ
 
 template <class Node_t>
 Permutation<typename Node_t::Real_t>& InteracList<Node_t>::Perm_C(int l, Mat_Type type, size_t indx){
+  assert(l>=0);
   size_t indx0=InteracClass(type, indx);
   Matrix     <Real_t>& M0      =mat->Mat   (l, type, indx0);
   Permutation<Real_t>& col_perm=mat->Perm_C(l, type, indx );

include/kernel.txx

@@ -922,6 +922,7 @@ void Kernel<T>::BuildMatrix(T* r_src, int src_cnt,
                  T* r_trg, int trg_cnt, T* k_out) const{
   int dim=3; //Only supporting 3D
   memset(k_out, 0, src_cnt*ker_dim[0]*trg_cnt*ker_dim[1]*sizeof(T));
+  #pragma omp parallel for
   for(int i=0;i<src_cnt;i++) //TODO Optimize this.
     for(int j=0;j<ker_dim[0];j++){
       std::vector<T> v_src(ker_dim[0],0);

include/legendre_rule.hpp

@@ -1,5 +1,5 @@
 
-# include <cstring>
+# include <string>
 
 #ifndef _LEGENDRE_RULE_HPP_
 #define _LEGENDRE_RULE_HPP_

include/pvfmm_common.hpp

@@ -28,7 +28,7 @@
 
 #define MAX_DEPTH 15
 
-#define BC_LEVELS 60
+#define BC_LEVELS 30
 
 #define RAD0 1.05 //Radius of upward equivalent (downward check) surface.
 #define RAD1 2.95 //Radius of downward equivalent (upward check) surface.