Fix issues with Vector, Matrix, Kernel, MPI_Tree

- For Vector and Matrix classes, it was not clear what Resize(...)
  function should do when own_data=false. The container will now e
  re-initialize with own_data=true when resizing to a larger size.

- Fix issues with detecting scale invariance and symmettry for certain
  kernels (like modified Laplace).

- Fix issue with file names of VTK files.

include/kernel.txx

@@ -66,38 +66,64 @@ void Kernel<T>::Initialize(bool verbose) const{
   T eps=1.0;
   while(eps+(T)1.0>1.0) eps*=0.5;
 
-  if(ker_dim[0]*ker_dim[1]>0){ // Determine scal
+  T scal=1.0;
+  if(ker_dim[0]*ker_dim[1]>0){ // Determine scaling
     Matrix<T> M_scal(ker_dim[0],ker_dim[1]);
     size_t N=1024;
     T eps_=N*eps;
 
     T src_coord[3]={0,0,0};
     std::vector<T> trg_coord1(N*COORD_DIM);
-    std::vector<T> trg_coord2(N*COORD_DIM);
-    for(size_t i=0;i<N;i++){
-      T x,y,z,r;
-      do{
-        x=(drand48()-0.5);
-        y=(drand48()-0.5);
-        z=(drand48()-0.5);
-        r=sqrt(x*x+y*y+z*z);
-      }while(r<0.25);
-      trg_coord1[i*COORD_DIM+0]=x;
-      trg_coord1[i*COORD_DIM+1]=y;
-      trg_coord1[i*COORD_DIM+2]=z;
+    Matrix<T> M1(N,ker_dim[0]*ker_dim[1]);
+    while(true){
+      T abs_sum=0;
+      for(size_t i=0;i<N/2;i++){
+        T x,y,z,r;
+        do{
+          x=(drand48()-0.5);
+          y=(drand48()-0.5);
+          z=(drand48()-0.5);
+          r=sqrt(x*x+y*y+z*z);
+        }while(r<0.25);
+        trg_coord1[i*COORD_DIM+0]=x*scal;
+        trg_coord1[i*COORD_DIM+1]=y*scal;
+        trg_coord1[i*COORD_DIM+2]=z*scal;
+      }
+      for(size_t i=N/2;i<N;i++){
+        T x,y,z,r;
+        do{
+          x=(drand48()-0.5);
+          y=(drand48()-0.5);
+          z=(drand48()-0.5);
+          r=sqrt(x*x+y*y+z*z);
+        }while(r<0.25);
+        trg_coord1[i*COORD_DIM+0]=x*1.0/scal;
+        trg_coord1[i*COORD_DIM+1]=y*1.0/scal;
+        trg_coord1[i*COORD_DIM+2]=z*1.0/scal;
+      }
+      for(size_t i=0;i<N;i++){
+        BuildMatrix(&src_coord [          0], 1,
+                    &trg_coord1[i*COORD_DIM], 1, &(M1[i][0]));
+        for(size_t j=0;j<ker_dim[0]*ker_dim[1];j++){
+          abs_sum+=fabs(M1[i][j]);
+        }
+      }
+      if(abs_sum>sqrt(eps) || scal<eps) break;
+      scal=scal*0.5;
     }
+
+    std::vector<T> trg_coord2(N*COORD_DIM);
+    Matrix<T> M2(N,ker_dim[0]*ker_dim[1]);
     for(size_t i=0;i<N*COORD_DIM;i++){
       trg_coord2[i]=trg_coord1[i]*0.5;
     }
-
-    T max_val=0;
-    Matrix<T> M1(N,ker_dim[0]*ker_dim[1]);
-    Matrix<T> M2(N,ker_dim[0]*ker_dim[1]);
     for(size_t i=0;i<N;i++){
-      BuildMatrix(&src_coord [          0], 1,
-                  &trg_coord1[i*COORD_DIM], 1, &(M1[i][0]));
       BuildMatrix(&src_coord [          0], 1,
                   &trg_coord2[i*COORD_DIM], 1, &(M2[i][0]));
+    }
+
+    T max_val=0;
+    for(size_t i=0;i<N;i++){
       for(size_t j=0;j<ker_dim[0]*ker_dim[1];j++){
         max_val=std::max<T>(max_val,M1[i][j]);
         max_val=std::max<T>(max_val,M2[i][j]);
@@ -170,7 +196,19 @@ void Kernel<T>::Initialize(bool verbose) const{
     T src_coord[3]={0,0,0};
     std::vector<T> trg_coord1(N*COORD_DIM);
     std::vector<T> trg_coord2(N*COORD_DIM);
-    for(size_t i=0;i<N;i++){
+    for(size_t i=0;i<N/2;i++){
+      T x,y,z,r;
+      do{
+        x=(drand48()-0.5);
+        y=(drand48()-0.5);
+        z=(drand48()-0.5);
+        r=sqrt(x*x+y*y+z*z);
+      }while(r<0.25);
+      trg_coord1[i*COORD_DIM+0]=x*scal;
+      trg_coord1[i*COORD_DIM+1]=y*scal;
+      trg_coord1[i*COORD_DIM+2]=z*scal;
+    }
+    for(size_t i=N/2;i<N;i++){
       T x,y,z,r;
       do{
         x=(drand48()-0.5);
@@ -178,9 +216,9 @@ void Kernel<T>::Initialize(bool verbose) const{
         z=(drand48()-0.5);
         r=sqrt(x*x+y*y+z*z);
       }while(r<0.25);
-      trg_coord1[i*COORD_DIM+0]=x;
-      trg_coord1[i*COORD_DIM+1]=y;
-      trg_coord1[i*COORD_DIM+2]=z;
+      trg_coord1[i*COORD_DIM+0]=x*1.0/scal;
+      trg_coord1[i*COORD_DIM+1]=y*1.0/scal;
+      trg_coord1[i*COORD_DIM+2]=z*1.0/scal;
     }
 
     for(size_t p_type=0;p_type<C_Perm;p_type++){ // For each symmetry transform

include/matrix.txx

@@ -119,8 +119,14 @@ void Matrix<T>::Swap(Matrix<T>& M){
 
 template <class T>
 void Matrix<T>::ReInit(size_t dim1, size_t dim2, T* data_, bool own_data_){
-  Matrix<T> tmp(dim1,dim2,data_,own_data_);
-  this->Swap(tmp);
+  if(own_data_ && own_data && dim[0]*dim[1]>=dim1*dim2){
+    if(dim[0]*dim[1]!=dim1*dim2) FreeDevice(false);
+    dim[0]=dim1; dim[1]=dim2;
+    if(data_) mem::memcopy(data_ptr,data_,dim[0]*dim[1]*sizeof(T));
+  }else{
+    Matrix<T> tmp(dim1,dim2,data_,own_data_);
+    this->Swap(tmp);
+  }
 }
 
 template <class T>
@@ -202,56 +208,26 @@ size_t Matrix<T>::Dim(size_t i) const{
 
 template <class T>
 void Matrix<T>::Resize(size_t i, size_t j){
-  if(dim[0]==i && dim[1]==j) return;
-  FreeDevice(false);
-  if(own_data){
-    if(data_ptr!=NULL){
-      mem::aligned_delete(data_ptr);
-#if !defined(__MIC__) || !defined(__INTEL_OFFLOAD)
-      Profile::Add_MEM(-dim[0]*dim[1]*sizeof(T));
-#endif
-
-    }
-  }
-  dim[0]=i;
-  dim[1]=j;
-  if(own_data){
-    if(dim[0]*dim[1]>0){
-      data_ptr=mem::aligned_new<T>(dim[0]*dim[1]);
-#if !defined(__MIC__) || !defined(__INTEL_OFFLOAD)
-      Profile::Add_MEM(dim[0]*dim[1]*sizeof(T));
-#endif
-    }else
-      data_ptr=NULL;
-  }
+  if(dim[0]*dim[1]!=i*j) FreeDevice(false);
+  if(dim[0]*dim[1]>=i*j){
+    dim[0]=i; dim[1]=j;
+  }else ReInit(i,j);
 }
 
 template <class T>
 void Matrix<T>::SetZero(){
-  if(dim[0]*dim[1]>0)
+  if(dim[0]*dim[1])
     memset(data_ptr,0,dim[0]*dim[1]*sizeof(T));
 }
 
 template <class T>
 Matrix<T>& Matrix<T>::operator=(const Matrix<T>& M){
   if(this!=&M){
-    FreeDevice(false);
-    if(own_data && dim[0]*dim[1]!=M.dim[0]*M.dim[1]){
-      if(data_ptr!=NULL){
-        mem::aligned_delete(data_ptr); data_ptr=NULL;
-#if !defined(__MIC__) || !defined(__INTEL_OFFLOAD)
-        Profile::Add_MEM(-dim[0]*dim[1]*sizeof(T));
-#endif
-      }
-      if(M.dim[0]*M.dim[1]>0){
-        data_ptr=mem::aligned_new<T>(M.dim[0]*M.dim[1]);
-#if !defined(__MIC__) || !defined(__INTEL_OFFLOAD)
-        Profile::Add_MEM(M.dim[0]*M.dim[1]*sizeof(T));
-#endif
-      }
+    if(dim[0]*dim[1]!=M.dim[0]*M.dim[1]) FreeDevice(false);
+    if(dim[0]*dim[1]<M.dim[0]*M.dim[1]){
+      ReInit(M.dim[0],M.dim[1]);
     }
-    dim[0]=M.dim[0];
-    dim[1]=M.dim[1];
+    dim[0]=M.dim[0]; dim[1]=M.dim[1];
     mem::memcopy(data_ptr,M.data_ptr,dim[0]*dim[1]*sizeof(T));
   }
   return *this;

include/mpi_tree.txx

@@ -2241,8 +2241,11 @@ void MPI_Tree<TreeNode>::Write2File(const char* fname, int lod){
   {
     // Extract filename from path.
     std::stringstream vtupath;
-    vtupath<<'/'<<fname<<'\0';
-    char *fname_ = (char*)strrchr(vtupath.str().c_str(), '/') + 1;
+    vtupath<<'/'<<fname;
+    std::string pathname = vtupath.str();
+    unsigned found = pathname.find_last_of("/\\");
+    std::string fname_ = pathname.substr(found+1);
+    //char *fname_ = (char*)strrchr(vtupath.str().c_str(), '/') + 1;
     //std::string fname_ = boost::filesystem::path(fname).filename().string().
     for(int i=0;i<np;i++) pvtufile<<"      <Piece Source=\""<<fname_<<std::setfill('0')<<std::setw(6)<<i<<".vtu\"/>\n";
   }

include/vector.hpp

@@ -70,7 +70,7 @@ class Vector{
 
   size_t Capacity() const;
 
-  void Resize(size_t dim_, bool fit_size=false);
+  void Resize(size_t dim_);
 
   void SetZero();
 

include/vector.txx

@@ -124,8 +124,13 @@ void Vector<T>::Swap(Vector<T>& v1){
 
 template <class T>
 void Vector<T>::ReInit(size_t dim_, T* data_, bool own_data_){
-  Vector<T> tmp(dim_,data_,own_data_);
-  this->Swap(tmp);
+  if(own_data_ && own_data && dim_<=capacity){
+    if(dim!=dim_) FreeDevice(false); dim=dim_;
+    if(data_) mem::memcopy(data_ptr,data_,dim*sizeof(T));
+  }else{
+    Vector<T> tmp(dim_,data_,own_data_);
+    this->Swap(tmp);
+  }
 }
 
 template <class T>
@@ -178,31 +183,10 @@ inline size_t Vector<T>::Capacity() const{
 }
 
 template <class T>
-void Vector<T>::Resize(size_t dim_,bool fit_size){
-  ASSERT_WITH_MSG(own_data || capacity>=dim_, "Resizing array beyond capacity when own_data=false.");
+void Vector<T>::Resize(size_t dim_){
   if(dim!=dim_) FreeDevice(false);
-
-  if((capacity>dim_ && !fit_size) || capacity==dim_ || !own_data){
-    dim=dim_;
-    return;
-  }
-
-  {
-    if(data_ptr!=NULL){
-      mem::aligned_delete(data_ptr); data_ptr=NULL;
-#if !defined(__MIC__) || !defined(__INTEL_OFFLOAD)
-      Profile::Add_MEM(-capacity*sizeof(T));
-#endif
-    }
-    capacity=dim_;
-    if(capacity>0){
-      data_ptr=mem::aligned_new<T>(capacity);
-#if !defined(__MIC__) || !defined(__INTEL_OFFLOAD)
-      Profile::Add_MEM(capacity*sizeof(T));
-#endif
-    }
-  }
-  dim=dim_;
+  if(capacity>=dim_) dim=dim_;
+  else ReInit(dim_);
 }
 
 template <class T>
@@ -213,26 +197,9 @@ void Vector<T>::SetZero(){
 
 template <class T>
 Vector<T>& Vector<T>::operator=(const Vector<T>& V){
-  ASSERT_WITH_MSG(own_data || capacity>=V.dim, "Resizing array beyond capacity when own_data=false.");
-
   if(this!=&V){
-    FreeDevice(false);
-    if(own_data && capacity<V.dim){
-      if(data_ptr!=NULL){
-        mem::aligned_delete(data_ptr); data_ptr=NULL;
-#if !defined(__MIC__) || !defined(__INTEL_OFFLOAD)
-        Profile::Add_MEM(-capacity*sizeof(T));
-#endif
-      }
-      capacity=V.dim;
-      if(capacity>0){
-        data_ptr=mem::aligned_new<T>(capacity);
-#if !defined(__MIC__) || !defined(__INTEL_OFFLOAD)
-        Profile::Add_MEM(capacity*sizeof(T));
-#endif
-      }
-    }
-    dim=V.dim;
+    if(dim!=V.dim) FreeDevice(false);
+    if(capacity<V.dim) ReInit(V.dim); dim=V.dim;
     mem::memcopy(data_ptr,V.data_ptr,dim*sizeof(T));
   }
   return *this;
@@ -240,26 +207,9 @@ Vector<T>& Vector<T>::operator=(const Vector<T>& V){
 
 template <class T>
 Vector<T>& Vector<T>::operator=(const std::vector<T>& V){
-  ASSERT_WITH_MSG(own_data || capacity>=V.size(), "Resizing array beyond capacity when own_data=false.");
-
   {
-    FreeDevice(false);
-    if(own_data && capacity<V.size()){
-      if(data_ptr!=NULL){
-        mem::aligned_delete(data_ptr); data_ptr=NULL;
-#if !defined(__MIC__) || !defined(__INTEL_OFFLOAD)
-        Profile::Add_MEM(-capacity*sizeof(T));
-#endif
-      }
-      capacity=V.size();
-      if(capacity>0){
-        data_ptr=mem::aligned_new<T>(capacity);
-#if !defined(__MIC__) || !defined(__INTEL_OFFLOAD)
-        Profile::Add_MEM(capacity*sizeof(T));
-#endif
-      }
-    }
-    dim=V.size();
+    if(dim!=V.size()) FreeDevice(false);
+    if(capacity<V.size()) ReInit(V.size()); dim=V.size();
     mem::memcopy(data_ptr,&V[0],dim*sizeof(T));
   }
   return *this;