Merge branch 'master' of file:///home/lavaux/Dropbox/gitRoot/CosmoToolbox

sample/CMakeLists.txt

@@ -26,3 +26,7 @@ if (HDF5_FOUND)
 	  add_executable(testReadFlash testReadFlash.cpp)
 	  target_link_libraries(testReadFlash ${tolink})
 endif (HDF5_FOUND)
+
+
+add_executable(testEskow testEskow.cpp)
+target_link_libraries(testEskow ${tolink})

sample/Hartmann_Matrix.txt

@@ -0,0 +1,7 @@
+6
+ 14.8253   -6.4243   7.8746  -1.2498  10.2733  10.2733 
+ -6.4243   15.1024  -1.1155  -0.2761  -8.2117  -8.2117 
+  7.8746   -1.1155  51.8519 -23.3482  12.5902  12.5902 
+ -1.2498   -0.2761 -23.3482  22.7962  -9.8958  -9.8958 
+ 10.2733   -8.2117  12.5902  -9.8958  21.0656  21.0656 
+ 10.2733   -8.2117  12.5902  -9.8958  21.0656  21.0656 

sample/testEskow.cpp

@@ -0,0 +1,65 @@
+#include <fstream>
+#include <cstring>
+#include <iostream>
+#include <iomanip>
+#include <vector>
+#include "eskow.hpp"
+
+using namespace std;
+
+double Hartmann_Matrix[6][6] = {
+  { 14.8253,   -6.4243,   7.8746,  -1.2498,  10.2733,  10.2733 },
+  { -6.4243,   15.1024,  -1.1155,  -0.2761,  -8.2117,  -8.2117 },
+  {  7.8746,   -1.1155,  51.8519, -23.3482,  12.5902,  12.5902 },
+  { -1.2498,   -0.2761, -23.3482,  22.7962,  -9.8958,  -9.8958 },
+  { 10.2733,   -8.2117,  12.5902,  -9.8958,  21.0656,  21.0656 },
+  { 10.2733,   -8.2117,  12.5902,  -9.8958,  21.0656,  21.0656 }
+};
+
+struct MatrixOp
+{
+  vector<double> M;
+  int N;
+
+  double& operator()(int i, int j)
+  {
+    return M[i*N + j];
+  }
+};
+
+int main(int argc, char **argv)
+{
+  MatrixOp M;
+  double norm_E;
+  ifstream fi(argv[1]);
+  ofstream f("eskowed.txt");
+  CholeskyEskow<double,MatrixOp> chol;
+
+  fi >> M.N;
+  M.M.resize(M.N*M.N);
+
+  for (int i = 0; i < M.N; i++)
+    {
+      for (int j = 0; j < M.N; j++)
+	{
+	  fi >> M(i,j);
+	  if (j > i)
+	    M(i,j) =0;
+	}
+    }
+
+  chol.cholesky(M, M.N, norm_E);
+
+  for (int i = 0; i < M.N; i++)
+    {
+      for (int j = 0; j < M.N; j++)
+	{
+	  if (j > i) 
+            f << "0 ";
+          else
+ 	    f << setprecision(25) << M(i,j) << " ";
+	}
+      f << endl;
+    }
+  return 0;
+}

src/cic.cpp

@@ -0,0 +1,205 @@
+#include <assert.h>
+#include <math.h>
+#include <inttypes.h>
+#include "cic.hpp"
+
+CICFilter::CICFilter(uint32_t N, double len)
+{
+  spatialLen = len;
+  szGrid = N;
+  totalSize = N*N*N;
+  densityGrid = new CICType[totalSize];
+  resetMesh();
+}
+
+CICFilter::~CICFilter()
+{
+  delete[] densityGrid;
+}
+  
+void CICFilter::resetMesh()
+{
+  for (uint32_t i = 0; i < totalSize; i++)
+    densityGrid[i] = 0;
+}
+
+void CICFilter::putParticles(CICParticles *particles, uint32_t N)
+{
+#if 0
+  uint32_t numCorners = 1 << NUMDIMS;
+
+  for (uint32_t i = 0; i < N; i++)
+    {
+      Coordinates xyz;
+      int32_t ixyz[NUMDIMS];
+      int32_t rxyz[NUMDIMS];
+      CICType alpha[NUMDIMS];
+      CICType beta[NUMDIMS];
+      for (int j = 0; j < NUMDIMS; j++)
+	{
+	  xyz[j] = (particles[i].coords[j] / spatialLen * szGrid);
+	  ixyz[j] = (int32_t)floor(xyz[j] - 0.5);
+	  beta[j] = xyz[j] - ixyz[j] - 0.5;
+	  alpha[j] = 1 - beta[j];
+	  if (ixyz[j] < 0)
+	    ixyz[j] = szGrid-1;
+	}
+
+      CICType tot_mass = 0;
+      for (int j = 0; j < numCorners; j++)
+	{
+	  CICType rel_mass = 1;
+	  uint32_t idx = 0;
+	  uint32_t mul = 1;
+	  uint32_t mul2 = 1;
+
+	  for (int k = 0; k < NUMDIMS; k++)
+	    {
+	      uint32_t ipos = ((j & mul2) != 0);
+
+	      if (ipos == 1)
+		{
+		  rel_mass *= beta[k];
+		}
+	      else
+		{
+		  rel_mass *= alpha[k];
+		}
+
+	      rxyz[k] = ixyz[k] + ipos;
+	      
+	      if (rxyz[k] >= szGrid)
+		idx += (rxyz[k] - szGrid) * mul;
+	      else
+		idx += rxyz[k] * mul;
+
+	      mul2 *= 2;
+	      mul *= szGrid;
+	    }
+
+	  assert(rel_mass > 0);
+	  assert(rel_mass < 1);
+	  assert(idx < totalSize);
+	  densityGrid[idx] += rel_mass * particles[i].mass;
+	  tot_mass += rel_mass;
+	}
+      assert(tot_mass < 1.1);
+      assert(tot_mass > 0.9);
+    }
+#endif
+#if 0  
+  for (uint32_t i = 0; i < N; i++)
+    {
+      Coordinates xyz;
+      int32_t ixyz[NUMDIMS];
+      for (int j = 0; j < NUMDIMS; j++)
+	{
+	  xyz[j] = (particles[i].coords[j] / spatialLen * szGrid);
+	  ixyz[j] = (int32_t)round(xyz[j] - 0.5);
+	  if (ixyz[j] < 0)
+	    ixyz[j] = szGrid-1;
+	  else if (ixyz[j] >= szGrid)
+	    ixyz[j] = 0;
+	}
+
+      uint32_t idx = ixyz[0] + ixyz[1] * szGrid + ixyz[2] * szGrid * szGrid;
+      densityGrid[idx] += particles[i].mass;      
+    }
+
+#endif
+
+  for (uint32_t i = 0; i < N; i++)
+    {
+      CICType x, y, z;
+      int32_t ix, iy, iz;
+      int32_t ix2, iy2, iz2;
+
+      x = particles[i].coords[0] / spatialLen * szGrid + 0.5;
+      y = particles[i].coords[1] / spatialLen * szGrid + 0.5;
+      z = particles[i].coords[2] / spatialLen * szGrid + 0.5;
+
+      if (x < 0)
+	x += szGrid;
+      if (y < 0)
+	y += szGrid;
+      if (z < 0)
+	z += szGrid;
+
+      ix = ((int32_t)floor(x));
+      iy = ((int32_t)floor(y));
+      iz = ((int32_t)floor(z));
+      
+      ix2 = (ix + 1) % szGrid;
+      iy2 = (iy + 1) % szGrid;
+      iz2 = (iz + 1) % szGrid;
+
+      CICType alpha_x = x - ix;
+      CICType alpha_y = y - iy;
+      CICType alpha_z = z - iz;
+
+      ix %= szGrid;
+      iy %= szGrid;
+      iz %= szGrid;
+      
+      assert(alpha_x >= 0);
+      assert(alpha_y >= 0);
+      assert(alpha_z >= 0);
+
+      CICType beta_x = 1 - alpha_x;
+      CICType beta_y = 1 - alpha_y;
+      CICType beta_z = 1 - alpha_z;
+
+      assert(beta_x >= 0);
+      assert(beta_y >= 0);
+      assert(beta_z >= 0);
+
+      CICType mass = particles[i].mass;
+      uint32_t idx;
+
+      // 000
+      idx = ix + (iy + iz * szGrid) * szGrid;
+      densityGrid[idx] +=
+	mass * beta_x * beta_y * beta_z;
+      
+      // 100
+      idx = ix2 + (iy + iz * szGrid) * szGrid;
+      densityGrid[idx] +=
+	mass * alpha_x * beta_y * beta_z;
+      
+      // 010
+      idx = ix + (iy2 + iz * szGrid) * szGrid;
+      densityGrid[idx] +=
+	mass * beta_x * alpha_y * beta_z;
+      
+      // 110
+      idx = ix2 + (iy2 + iz * szGrid) * szGrid;
+      densityGrid[idx] +=
+	mass * alpha_x * alpha_y * beta_z;
+
+      // 001
+      idx = ix + (iy + iz2 * szGrid) * szGrid;
+      densityGrid[idx] +=
+	mass * beta_x * beta_y * alpha_z;
+
+      // 101
+      idx = ix2 + (iy + iz2 * szGrid) * szGrid;
+      densityGrid[idx] +=
+	mass * alpha_x * beta_y * alpha_z;
+
+      // 011
+      idx = ix + (iy2 + iz2 * szGrid) * szGrid;
+      densityGrid[idx] +=
+	mass * beta_x * alpha_y * alpha_z;
+
+      // 111
+      idx = ix2 + (iy2 + iz2 * szGrid) * szGrid;
+      densityGrid[idx] +=
+	mass * alpha_x * alpha_y * alpha_z;
+    }
+}
+  
+void CICFilter::getDensityField(CICType*& field, uint32_t& res)
+{
+  field = densityGrid;
+  res = totalSize;
+}

src/cic.hpp

@@ -0,0 +1,35 @@
+#ifndef __CICFILTER_HPP
+#define __CICFILTER_HPP
+
+#include "CosmoTool/config.hpp"
+#include <inttypes.h>
+
+using namespace CosmoTool;
+
+typedef float CICType;
+
+  typedef struct
+  {
+    float mass;
+    Coordinates coords;
+  } CICParticles;
+
+  class CICFilter
+  {
+  public:
+    CICFilter(uint32_t resolution, double spatialLen);
+    ~CICFilter();
+
+    void resetMesh();
+    void putParticles(CICParticles *particles, uint32_t N);    
+
+    void getDensityField(CICType*& field, uint32_t& res);
+
+  protected:
+    CICType *densityGrid;
+    double spatialLen;
+    uint32_t totalSize;
+    uint32_t szGrid;
+  };
+
+#endif

src/eskow.hpp

@@ -0,0 +1,271 @@
+#ifndef __ESKOW_CHOLESKY_HPP
+#define __ESKOW_CHOLESKY_HPP
+
+#include <cmath>
+#include <vector>
+#include "mach.hpp"
+
+/* Implementation of Schnabel & Eskow, 1999, Vol. 9, No. 4, pp. 1135-148, SIAM J. OPTIM. */
+
+template<typename T, typename A>
+class CholeskyEskow
+{
+private:
+  static const bool verbose_eskow = true;
+  T tau, tau_bar, mu;
+
+  void print_matrix(A& m, int N)
+  {
+    using std::cout;
+    using std::endl;
+    using std::setprecision;
+
+    if (verbose_eskow)
+      {
+
+	for (int i = 0; i < N; i++)
+	  {
+	    for (int j = 0; j < N; j++)
+	      {
+		cout.width(6);
+		cout << setprecision(5) << m(i,j) << " ";
+	      }
+	    cout << endl;
+	  }
+	cout << endl;
+      }
+  }
+
+  T max_diag(A& m, int j, int N)
+  {
+    T maxval = std::abs(m(j,j));
+
+    for (int k = j+1; k < N; k++)
+      {
+	maxval = std::max(maxval, std::abs(m(k,k)));
+      }
+    return maxval;
+  }
+
+  void minmax_diag(A& m, int j, int N, T& minval, T& maxval, int& i_min, int& i_max)
+  {
+    i_min = i_max = j;
+    minval = maxval = m(j,j);
+
+    for (int k = j+1; k < N; k++)
+      {
+	maxval = std::max(maxval, m(k,k));
+	minval = std::min(minval, m(k,k));
+      }
+  
+    for (int k = j; k < N; k++)
+      {
+	if (m(k,k) == minval && i_min < 0)
+	  i_min = k;
+	if (m(k,k) == maxval && i_max < 0)
+	  i_max = k;
+      }
+  }
+
+  void swap_rows(A& m, int N, int i0, int i1)
+  {
+    for (int r = 0; r < N; r++)
+      std::swap(m(r,i0), m(r,i1));
+  }
+
+  void swap_cols(A& m, int N, int i0, int i1)
+  {
+    for (int c = 0; c < N; c++)
+      std::swap(m(i0,c), m(i1,c));
+  }
+
+  T square(T x)
+  {
+    return x*x;
+  }
+
+  T min_row(A& m, int j, int N)
+  {
+    T a = 1/m(j,j);
+    T v = m(j+1,j+1) - square(m(j+1,j))*a;
+    
+    for (int i = j+2; i < N; i++)
+      {
+	v = std::min(v, m(i, i) - square(m(i,j))*a);
+      }
+
+    return v;
+  }
+
+  int g_max(const std::vector<T>& g, int j, int N)
+  {
+    T a = g[j];
+    int k = j;
+
+    for (int i = j+1; i < N; i++)
+      {
+	if (a < g[i])
+	  {
+	    a = g[i];
+	    k = i;
+	  }
+      }
+    return k;
+  }
+
+public:
+  CholeskyEskow()
+  {
+    tau = std::pow(mach_epsilon<T>(), 1./3);
+    tau_bar = std::pow(mach_epsilon<T>(), 2./3);
+    mu=0.1;
+  }
+
+  void cholesky(A& m, int N, T& norm_E)
+  {
+    bool phaseone = true;
+    T gamma = max_diag(m, 0, N);
+    int j;
+
+    norm_E = 0;
+    
+    for (j = 0; j < N && phaseone; j++)
+      {
+	T minval, maxval;
+	int i_min, i_max;
+	
+	print_matrix(m, N);
+     
+	minmax_diag(m, j, N, minval, maxval, i_min, i_max);
+	if (maxval < tau_bar*gamma || minval < -mu*maxval)
+	  {
+	    phaseone = false;
+	    break;
+	  }
+      
+	if (i_max != j)	  
+	  {
+	    std::cout << "Have to swap i=" << i_max << " and j=" << j << std::endl;
+	    swap_cols(m, N, i_max, j);
+	    swap_rows(m, N, i_max, j);
+
+	  }
+      
+	if (min_row(m, j, N) < -mu*gamma)
+	  {
+	    phaseone = false;
+	    break;
+	  }
+      
+	T L_jj = std::sqrt(m(j,j));
+
+	m(j,j) = L_jj;
+	for (int i = j+1; i < N; i++)
+	  {
+	    m(i,j) /= L_jj;
+	    for (int k = j+1; k <= i; k++)
+	      m(i,k) -= m(i,j)*m(k,j);
+	  }
+      }
+
+
+    if (!phaseone && j == N-1)
+      {
+	T A_nn = m(N-1,N-1);
+	T delta = -A_nn + std::max(tau*(-A_nn)/(1-tau), tau_bar*gamma);
+      
+	m(N-1,N-1) = std::sqrt(m(N-1,N-1) + delta);	  	  
+      }
+
+    
+
+    if (!phaseone && j < (N-1))
+      {
+	std::cout << "Phase two ! (j=" << j << ")" << std::endl;
+
+	int k = j-1;
+	std::vector<T> g(N);
+
+	for (int i = k+1; i < N; i++)
+	  {
+	    g[i] = m(i,i);
+	    for (int j = k+1; j < i; j++)
+	      g[i] -= std::abs(m(i,j));
+	    for (int j = i+1; j < N; j++)
+	      g[i] -= std::abs(m(j,i));
+	  }
+
+	T delta, delta_prev = 0;
+	
+	for (int j = k+1; j < N-2; j++)
+	  {
+	    int i = g_max(g, j, N);
+	    T norm_j;
+
+	    print_matrix(m, N);
+
+	    if (i != j)
+	      {
+		swap_cols(m, N, i, j);
+		swap_rows(m, N, i, j);		
+	      }
+
+
+	    for (int i = j+1; j < N; j++)
+	      {
+		norm_j += std::abs(m(i,j));
+	      }
+	    
+	    delta = std::max(delta_prev, std::max((T)0, -m(j,j) + std::max(norm_j,tau_bar*gamma)));
+	    if (delta > 0)
+	      {
+		m(j,j) += delta;
+		delta_prev = delta;
+	      }
+	    
+	    if (m(j,j) != norm_j)
+	      {
+		T temp = 1 - norm_j/m(j,j);
+		
+		for (int i = j+1; j < N; j++)
+		  {
+		    g[i] += std::abs(m(i,j))*temp;
+		  }
+	      }
+
+	    // Now we do the classic cholesky iteration
+	    T L_jj = std::sqrt(m(j,j));
+	    
+	    m(j,j) = L_jj;
+	    for (int i = j+1; i < N; i++)
+	      {
+		m(i,j) /= L_jj;
+		for (int k = j+1; k <= i; k++)
+		  m(i,k) -= m(i,j)*m(k,j);
+	      }
+	  }
+
+	// The final 2x2 submatrix is special
+	T A00 = m(N-2, N-2), A01 = m(N-2, N-1), A11 = m(N-1,N-1);
+	T sq_DELTA = std::sqrt(square(A00-A11) + square(A01));
+	T lambda_hi = 0.5*((A00+A11) + sq_DELTA);
+	T lambda_lo = 0.5*((A00+A11) - sq_DELTA);
+       
+	delta = std::max(std::max((T)0, -lambda_lo + std::max(tau*sq_DELTA/(1-tau), tau_bar*gamma)),delta_prev);
+	if (delta > 0)
+	  {
+	    m(N-1,N-1) += delta;
+	    m(N,N) += delta;
+	    delta_prev = delta;
+	  }
+	m(N-2,N-2) = A00 = std::sqrt(A00);
+	m(N-1,N-2) = (A01 /= A00);
+	m(N-1,N-1) = std::sqrt(A11-A01*A01);
+	norm_E = delta_prev;
+      }
+  }
+
+};
+
+
+#endif

src/loadGadget.cpp

@@ -273,3 +273,73 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id, int loadflags, i
 }
 
 
+
+void CosmoTool::writeGadget(const char *fname, SimuData *data, int GadgetFormat)
+{
+  UnformattedWrite *f;
+  int npart[6];
+  float mass[6];
+
+  if (data->Pos[0] == 0 || data->Vel[0] == 0 || data->Id == 0)
+    return;
+
+  f = new UnformattedWrite(fname);
+  if (f == 0)
+    return;
+
+  for (int i = 0; i < 6; i++)
+    {
+      npart[i] = 0;
+      mass[i] = 0;
+    }
+
+  npart[1] = data->NumPart;
+  
+  f->beginCheckpoint();
+  for (int i = 0; i < 6; i++)
+    f->writeInt32(npart[i]);
+  for (int i = 0; i < 6; i++)
+    f->writeReal64(mass[i]);
+
+  f->writeReal64(data->time);
+  f->writeReal64(1/data->time-1);
+  f->writeInt32(0);
+  f->writeInt32(0);
+
+  for (int i = 0; i < 6; i++)
+    f->writeInt32(npart[i]);
+  f->writeInt32(0);
+  f->writeInt32(1);
+  f->writeReal64(data->BoxSize);
+  f->writeReal64(data->Omega_M);
+  f->writeReal64(data->Omega_Lambda);
+  f->writeReal64(data->Hubble);
+  f->endCheckpoint();
+ 
+  f->beginCheckpoint();
+  for(int n = 0; n < data->NumPart; n++) {
+    for (int k = 0; k < 3; k++)
+      f->writeReal32(data->Pos[k][n]);
+  }
+  f->endCheckpoint();
+
+  float velmul = 1.0;
+  if (GadgetFormat == 1)
+    velmul = sqrt(data->time);
+
+  f->beginCheckpoint();
+  for(int n = 0; n < data->NumPart; n++) {
+    for (int k = 0; k < 3; k++)
+      f->writeReal32(data->Vel[k][n]/velmul);
+  }
+  f->endCheckpoint();
+
+  f->beginCheckpoint();
+  for(int n = 0; n < data->NumPart; n++)
+    {
+      f->writeReal32(data->Id[n]);
+    }
+  f->endCheckpoint();
+}
+
+  

src/loadGadget.hpp

@@ -9,7 +9,10 @@ namespace CosmoTool {
   PurePositionData *loadGadgetPosition(const char *fname);
 
   SimuData *loadGadgetMulti(const char *fname, int id, int flags, int GadgetFormat = 1);
-  
+ 
+  // Only single snapshot supported
+  void writeGadget(const char *fname, SimuData *data, int GadgetFormat = 1);
+ 
 };
 
 #endif

src/mach.hpp

@@ -0,0 +1,20 @@
+#ifndef __COSMO_MACHINE_TEST_HPP
+#define __COSMO_MACHINE_TEST_HPP
+
+#include <iostream>
+
+template<typename T>
+T mach_epsilon()
+{
+  T eps = (T)1;
+
+  do
+    {
+      eps /= 2;
+    }
+  while ((T)(1 + (eps/2)) != (T)1);
+
+  return eps;
+}
+
+#endif