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

CMakeLists.txt

@@ -11,16 +11,32 @@ option(BUILD_SHARED_LIBS "Build shared libraries." OFF)
 option(BUILD_STATIC_LIBS "Build static libraries." ON)
 
 find_path(NETCDF_INCLUDE_PATH NAMES netcdf.h)
+find_path(NETCDFCPP_INCLUDE_PATH NAMES netcdfcpp.h netcdf)
 find_path(GSL_INCLUDE_PATH NAMES gsl/gsl_blas.h)
 
+IF(EXISTS ${NETCDFCPP_INCLUDE_PATH}/netcdf)
+  SET(FOUND_NETCDF4 1)
+  FILE(WRITE ${CMAKE_BINARY_DIR}/src/ctool_netcdf_ver.hpp "#define NETCDFCPP4 1")
+ELSE(EXISTS ${NETCDFCPP_INCLUDE_PATH}/netcdf)
+  SET(FOUND_NETCDF3 1)
+  FILE(WRITE ${CMAKE_BINARY_DIR}/src/ctool_netcdf_ver.hpp "#undef NETCDFCPP4")
+ENDIF(EXISTS ${NETCDFCPP_INCLUDE_PATH}/netcdf)
+
 find_library(NETCDF_LIBRARY netcdf)
-find_library(NETCDFCPP_LIBRARY netcdf_c++)
+find_library(NETCDFCPP_LIBRARY NAMES netcdf_c++ netcdf_c++4)
 find_library(GSL_LIBRARY gsl)
 find_library(GSLCBLAS_LIBRARY gslcblas)
 
-set(HDF5_FIND_COMPONENTS CXX)
+set(HDF5_FIND_COMPONENTS HL CXX)
+if(HDF5_ROOTDIR)
+  SET(ENV{HDF5_ROOT} ${HDF5_ROOTDIR})
+endif(HDF5_ROOTDIR)
 include(FindHDF5)
 
+include(FindPkgConfig)
+
+pkg_check_modules(FFTW3 fftw3>=3.3)
+pkg_check_modules(EIGEN3 eigen3)
 
 include(FindPackageHandleStandardArgs)
 set(NETCDF_FIND_REQUIRED TRUE)

sample/CMakeLists.txt

@@ -1,5 +1,5 @@
 SET(tolink ${GSL_LIBRARIES} CosmoTool ${CosmoTool_LIBS})
-include_directories(${CMAKE_SOURCE_DIR}/src ${NETCDF_INCLUDE_PATH} ${GSL_INCLUDE_PATH})
+include_directories(${CMAKE_SOURCE_DIR}/src ${FFTW3_INCLUDE_DIRS} ${EIGEN3_INCLUDE_DIRS} ${NETCDF_INCLUDE_PATH} ${GSL_INCLUDE_PATH})
 
 add_executable(testBQueue testBQueue.cpp)
 target_link_libraries(testBQueue ${tolink})
@@ -36,3 +36,11 @@ target_link_libraries(testEskow ${tolink})
 
 add_executable(testAlgo testAlgo.cpp)
 target_link_libraries(testAlgo ${tolink})
+
+add_executable(testBSP testBSP.cpp)
+target_link_libraries(testBSP ${tolink})
+
+if (FFTW3_FOUND AND EIGEN3_FOUND)
+  add_executable(test_fft_calls test_fft_calls)
+  target_link_libraries(test_fft_calls ${tolink} ${FFTW3_LIBRARIES})
+endif (FFTW3_FOUND AND EIGEN3_FOUND)

sample/testBSP.cpp

@@ -1,3 +1,4 @@
+#if 0
 #include "bsp_simple.hpp"
 
 int main(int argc, char **argv)
@@ -10,3 +11,5 @@ int main(int argc, char **argv)
 
   return 0;
 }
+#endif
+int main() {}

sample/testEskow.cpp

@@ -50,6 +50,8 @@ int main(int argc, char **argv)
 
   chol.cholesky(M, M.N, norm_E);
 
+  cout << "norm_E = " << norm_E << endl;
+
   for (int i = 0; i < M.N; i++)
     {
       for (int j = 0; j < M.N; j++)

sample/testReadFlash.cpp

@@ -9,7 +9,7 @@ using namespace std;
 
 int main () {
  
-  char* filename = "lss_read_hdf5_chk_0000";
+  const char* filename = "lss_read_hdf5_chk_0000";
  
   SimuData* data = CosmoTool::loadFlashMulti(filename, 0, 0);
 

sample/test_fft_calls.cpp

@@ -0,0 +1,12 @@
+#include "fourier/euclidian.hpp"
+
+using namespace CosmoTool;
+
+int main()
+{
+  EuclidianFourierTransform_2d<double> dft(128,128,1.0,1.0);
+
+  dft.realSpace().eigen().setRandom();
+  dft.analysis();
+  return 0;
+}

src/CMakeLists.txt

@@ -6,11 +6,19 @@ SET(CosmoTool_SRCS
 	loadRamses.cpp
 	octTree.cpp
 	powerSpectrum.cpp
-	yorick.cpp
 	miniargs.cpp
 	growthFactor.cpp
 )
 
+IF(FOUND_NETCDF3)
+	SET(CosmoTool_SRCS ${CosmoTool_SRCS} yorick_nc3.cpp)
+ELSE(FOUND_NETCDF3)
+IF(FOUND_NETCDF4)
+	SET(CosmoTool_SRCS ${CosmoTool_SRCS} yorick_nc4.cpp)
+ENDIF(FOUND_NETCDF4)
+ENDIF(FOUND_NETCDF3)
+
+
 if (HDF5_FOUND)
   set(CosmoTool_SRCS ${CosmoTool_SRCS}
     h5_readFlash.cpp
@@ -46,15 +54,14 @@ SET(CosmoTool_SRCS ${CosmoTool_SRCS}
 	growthFactor.hpp
 )
 
-include_directories(${GSL_INCLUDE_PATH} ${NETCDF_INCLUDE_PATH})
+include_directories(${GSL_INCLUDE_PATH} ${NETCDF_INCLUDE_PATH} ${NETCDFCPP_INCLUDE_PATH} ${CMAKE_BINARY_DIR}/src)
 
-set(CosmoTool_LIBS  ${NETCDF_LIBRARY} ${NETCDFCPP_LIBRARY} ${GSL_LIBRARIES})
+set(CosmoTool_LIBS  ${NETCDFCPP_LIBRARY} ${NETCDF_LIBRARY} ${GSL_LIBRARIES})
 if (HDF5_FOUND)
   set(CosmoTool_LIBS ${CosmoTool_LIBS} ${HDF5_CXX_LIBRARIES} ${HDF5_LIBRARIES})
   include_directories(${HDF5_INCLUDE_DIRS})
 endif (HDF5_FOUND)
 
-message("${CosmoTool_LIBS}")
 set(CosmoTool_LIBS ${CosmoTool_LIBS} PARENT_SCOPE)
 
 if (BUILD_SHARED_LIBS)

src/bsp_simple.hpp

@@ -154,6 +154,8 @@ namespace CosmoTool
 	f.data = data;
 	insert(f);
       }
+
+      bool inside(const typename space_t::coord_t& p) const;
     };
 
   };

src/bsp_simple.tcc

@@ -98,7 +98,19 @@ namespace CosmoTool
 	*(*i) = current;
 
       allocated.push(current);
-   }
+    }
+
+    template<typename T, typename PType, int N>
+    bool BSP<T,PType,N>::inside(const typename space_t::coord_t& p) const
+    {
+      node_t *current = root;
+
+      do
+	{
+	}
+      while();
+      current
+    }
 
   };
 

src/fourier/base_types.hpp

@@ -0,0 +1,102 @@
+#ifndef __BASE_FOURIER_TYPES_HPP
+#define __BASE_FOURIER_TYPES_HPP
+
+#include <string>
+#include <Eigen/Dense>
+#include <complex>
+
+namespace CosmoTool
+{
+  template<typename T>
+  class FourierMap
+  {
+  protected:
+    FourierMap() {}
+    
+  public:
+    typedef Eigen::Matrix<T, 1, Eigen::Dynamic> VecType;
+    typedef Eigen::Map<VecType, Eigen::Aligned> MapType;
+    typedef Eigen::Map<VecType const, Eigen::Aligned> ConstMapType;
+
+    virtual ~FourierMap() {}
+
+    virtual const T* data() const = 0;
+    virtual T* data()  = 0;
+
+    virtual long size() const = 0;
+
+    MapType eigen()
+    {
+      return MapType(data(), size());
+    }
+
+    ConstMapType eigen() const
+    {
+      return ConstMapType(data(), size());
+    }
+ 
+    void scale(const T& factor)
+    {
+      MapType m(data(), size());
+      m *= factor;
+    }
+
+    void scale(const FourierMap<T> *map2)
+    {
+      assert(size() == map2->size());
+      MapType m(data(), size());
+      MapType m2(map2->data(), map2->size());
+      m *= m2;
+    }
+
+    void add(const T&  factor)
+    {
+      eigen() += factor;
+    }
+
+    void add(const FourierMap<T> *map2)
+    {
+      assert(size() == map2->size());
+      MapType m(data(), size());
+      MapType m2(map2->data(), map2->size());
+
+      eigen() += map2->eigen();
+    }
+
+    virtual FourierMap<T> *copy() const
+    {
+      FourierMap<T> *m = this->mimick();
+      
+      m->eigen() = this->eigen();
+      return m;
+    }
+
+    virtual FourierMap<T> *mimick() const = 0;
+  };
+
+  template<typename T>
+  class FourierTransform
+  {
+  protected:
+    FourierTransform() {}
+  public:
+    virtual ~FourierTransform() { }
+    
+    virtual const FourierMap<std::complex<T> >& fourierSpace() const = 0;
+    virtual FourierMap<std::complex<T> >& fourierSpace() = 0;
+
+    virtual const FourierMap<T>& realSpace() const = 0;
+    virtual FourierMap<T>& realSpace() = 0;
+
+    virtual FourierTransform<T> *mimick() const = 0;
+
+    virtual void analysis() = 0;
+    virtual void synthesis() = 0;
+    virtual void analysis_conjugate() = 0;
+    virtual void synthesis_conjugate() = 0;
+  };
+
+
+};
+
+#endif

src/fourier/euclidian.hpp

@@ -0,0 +1,200 @@
+#ifndef __COSMOTOOL_FOURIER_EUCLIDIAN_HPP
+#define __COSMOTOOL_FOURIER_EUCLIDIAN_HPP
+
+#include <vector>
+#include <boost/shared_ptr.hpp>
+#include "base_types.hpp"
+#include "fft/fftw_calls.hpp"
+
+namespace CosmoTool
+{
+  
+  template<typename T>
+  class EuclidianFourierMap: public FourierMap<T>
+  {
+  private:
+    boost::shared_ptr<T> m_data;
+    std::vector<int> m_dims;
+    long m_size;
+  public:
+    EuclidianFourierMap(boost::shared_ptr<T> indata, std::vector<int> indims)
+    {
+      m_data = indata;
+      m_dims = indims;
+      m_size = 1;
+      for (int i = 0; i < m_dims.size(); i++)
+	m_size *= m_dims[i];
+    }
+
+    virtual ~EuclidianFourierMap()
+    {
+    }
+
+    virtual const T *data() const { return m_data.get(); }
+    virtual T *data() { return m_data.get(); }
+    virtual long size() const  { return m_size; } 
+    
+    virtual FourierMap<T> *copy() const
+    {
+      FourierMap<T> *m = this->mimick();
+      m->eigen() = this->eigen();
+      return m;
+    }
+
+    virtual FourierMap<T> *mimick() const
+    {
+      return new EuclidianFourierMap<T>(
+		      boost::shared_ptr<T>((T *)fftw_malloc(sizeof(T)*m_size), 
+					   std::ptr_fun(fftw_free)),
+		       m_dims);
+    }
+  };
+
+  
+  template<typename T>
+  class EuclidianFourierTransform: public FourierTransform<T>
+  {
+  private:
+    typedef FFTW_Calls<T> calls;
+    EuclidianFourierMap<T> *realMap;
+    EuclidianFourierMap<std::complex<T> > *fourierMap;
+    typename calls::plan_type m_analysis, m_synthesis;
+    double volume;
+    long N;
+    std::vector<int> m_dims;
+    std::vector<double> m_L;
+  public:
+    EuclidianFourierTransform(const std::vector<int>& dims, const std::vector<double>& L)
+    {
+      assert(L.size() == dims.size());
+
+      m_dims = dims;
+      m_L = L;
+
+      N = 1;
+      volume = 1;
+      for (int i = 0; i < dims.size(); i++)
+	{
+	  N *= dims[i];
+	  volume *= L[i];
+	}
+
+      realMap = new EuclidianFourierMap<T>(
+		   boost::shared_ptr<T>(calls::alloc_real(N),
+					std::ptr_fun(calls::free)),
+		   dims);
+      fourierMap = new EuclidianFourierMap<std::complex<T> >(
+		   boost::shared_ptr<std::complex<T> >((std::complex<T>*)calls::alloc_complex(N),
+						       std::ptr_fun(calls::free)), 
+		   dims);
+      m_analysis = calls::plan_dft_r2c(dims.size(), &dims[0],
+				       realMap->data(), (typename calls::complex_type *)fourierMap->data(),
+				       FFTW_MEASURE);
+      m_synthesis = calls::plan_dft_c2r(dims.size(), &dims[0],
+					(typename calls::complex_type *)fourierMap->data(), realMap->data(),
+					FFTW_MEASURE);
+    }
+    
+    virtual ~EuclidianFourierTransform()
+    {
+      delete realMap;
+      delete fourierMap;
+      calls::destroy_plan(m_synthesis);
+      calls::destroy_plan(m_analysis);
+    }
+
+    void synthesis()
+    {
+      calls::execute(m_synthesis);
+      realMap->scale(1/volume);
+    }
+    
+    void analysis()
+    {
+      calls::execute(m_analysis);
+      fourierMap->scale(volume/N);
+    }
+    
+    void synthesis_conjugate()
+    {
+      calls::execute(m_analysis);
+      fourierMap->scale(1/volume);
+    }
+    
+    void analysis_conjugate()
+    {
+      calls::execute(m_synthesis);
+      realMap->scale(volume/N);
+    }
+
+    const FourierMap<std::complex<T> >& fourierSpace() const
+    {
+      return *fourierMap;
+    }
+
+    FourierMap<std::complex<T> >& fourierSpace()
+    {
+      return *fourierMap;
+    }
+
+    const FourierMap<T>& realSpace() const
+    {
+      return *realMap;
+    }
+
+    FourierMap<T>& realSpace()
+    {
+      return *realMap;
+    }
+
+    FourierTransform<T> *mimick() const
+    {
+      return new EuclidianFourierTransform(m_dims, m_L);
+    }
+  };
+
+  template<typename T>
+  class EuclidianFourierTransform_2d: public EuclidianFourierTransform<T>
+  {
+  private:
+    template<typename T2>
+    static std::vector<T2> make_2d_vector(T2 a, T2 b)
+    {
+      T2 arr[2] = { a, b};
+      return std::vector<T2>(&arr[0], &arr[2]);
+    }
+  public:
+    EuclidianFourierTransform_2d(int Nx, int Ny, double Lx, double Ly)
+      : EuclidianFourierTransform<T>(make_2d_vector<int>(Nx, Ny), make_2d_vector<double>(Lx, Ly))
+    {
+    }
+
+    virtual ~EuclidianFourierTransform_2d() {}
+
+  };
+
+  template<typename T>
+  class EuclidianFourierTransform_3d: public EuclidianFourierTransform<T>
+  {
+  private:
+    template<typename T2>
+    static std::vector<T2> make_3d_vector(T2 a, T2 b, T2 c)
+    {
+      T2 arr[2] = { a, b, c};
+      return std::vector<T2>(&arr[0], &arr[3]);
+    }
+
+  public:
+    EuclidianFourierTransform_3d(int Nx, int Ny, int Nz, double Lx, double Ly, double Lz)
+      : EuclidianFourierTransform<T>(make_3d_vector<int>(Nx, Ny, Nz), make_3d_vector<double>(Lx, Ly, Lz))
+    {
+    }
+
+    virtual ~EuclidianFourierTransform_3d() {}
+  };
+
+
+
+};
+
+#endif

src/fourier/fft/fftw_calls.hpp

@@ -0,0 +1,75 @@
+#ifndef __FFTW_UNIFIED_CALLS_HPP
+#define __FFTW_UNIFIED_CALLS_HPP
+
+#include <fftw3.h>
+
+namespace CosmoTool
+{
+
+static inline void init_fftw_wisdom()
+{
+  fftw_import_system_wisdom();
+  fftw_import_wisdom_from_filename("fft_wisdom");
+}
+
+static inline void save_fftw_wisdom()
+{
+  fftw_export_wisdom_to_filename("fft_wisdom");
+}
+
+template<typename T> class FFTW_Calls {};
+
+
+#define FFTW_CALLS_BASE(rtype, prefix) \
+  template<>				\
+class FFTW_Calls<rtype> {		\
+public: \
+  typedef rtype real_type; \
+  typedef prefix ## _complex complex_type; \
+  typedef prefix ## _plan plan_type; \
+  \
+  static complex_type *alloc_complex(int N) { return prefix ## _alloc_complex(N); } \
+  static real_type *alloc_real(int N) { return prefix ## _alloc_real(N); } \
+  static void free(void *p) { fftw_free(p); } \
+\
+  static void execute(plan_type p) { prefix ## _execute(p); } \
+  static plan_type plan_dft_r2c_2d(int Nx, int Ny,  \
+		       real_type *in, complex_type *out, \
+		       unsigned flags) \
+  { \
+    return prefix ## _plan_dft_r2c_2d(Nx, Ny, in, out, \
+				flags); \
+  } \
+  static plan_type plan_dft_c2r_2d(int Nx, int Ny,  \
+                       complex_type *in, real_type *out, \
+                       unsigned flags) \
+  { \
+    return prefix ## _plan_dft_c2r_2d(Nx, Ny, in, out, \
+                                flags); \
+  } \
+  static plan_type plan_dft_r2c_3d(int Nx, int Ny, int Nz, \
+                                   real_type *in, complex_type *out, \
+                                   unsigned flags) \
+  { \
+    return prefix ## _plan_dft_r2c_3d(Nx, Ny, Nz, in, out, flags); \
+  } \
+  static plan_type plan_dft_r2c(int rank, const int *n, real_type *in, \
+                                complex_type *out, unsigned flags) \
+  { \
+    return prefix ## _plan_dft_r2c(rank, n, in, out, flags); \
+  } \
+  static plan_type plan_dft_c2r(int rank, const int *n, complex_type *in, \
+                                real_type *out, unsigned flags) \
+  { \
+    return prefix ## _plan_dft_c2r(rank, n, in, out, flags); \
+  } \
+  static void destroy_plan(plan_type plan) { prefix ## _destroy_plan(plan); } \
+}
+
+
+FFTW_CALLS_BASE(double, fftw);
+FFTW_CALLS_BASE(float, fftwf);
+
+};
+
+#endif

src/fourier/healpix.hpp

@@ -0,0 +1,135 @@
+#ifndef __COSMOTOOL_FOURIER_HEALPIX_HPP
+#define __COSMOTOOL_FOURIER_HEALPIX_HPP
+
+#include <boost/bind.hpp>
+#include <boost/shared_ptr.hpp>
+#include <alm.h>
+#include <healpix_base.h>
+#include <psht_cxx.h>
+
+namespace CosmoTool
+{
+    
+  template<typename T>
+  class HealpixFourierMap: public FourierMap<T>, public Healpix_Base
+  {
+  private:
+    T *m_data;
+    Eigen::aligned_allocator<T> alloc;
+  public:
+    HealpixFourierMap(long nSide)
+      : Healpix_Base(RING, nSide, SET_NSIDE)
+    {
+      m_data = alloc.allocate(Npix);
+    }
+
+    virtual ~HealpixFourierMap()
+    {
+      alloc.deallocate(m_data);
+    }
+
+    virtual const T* data() const { return m_data; }
+    virtual T *data() { return m_data; }
+    virtual long size() const { return Npix(); }
+    
+    virtual FourierMap<T> *mimick() const
+    {
+      return new HealpixFourierMap<T>(Nside());
+    }
+  };
+
+
+  template<typename T>
+  class HealpixFourierALM: public FourierMap<std::complex<T> >, public Alm_Base
+  {
+  private:
+    std::complex<T> *alms;
+    long size;
+    Eigen::aligned_allocator<std::complex<T> > alloc;
+  public:
+    HealpixFourierALM(long Lmax, long Mmax)
+      : Alm_Base(Lmax, Mmax)
+    {
+      size = Num_Alms(Lmax, Mmax);
+      alms = alloc.allocate(size);
+    }
+
+    virtual ~HealpixFourierALM()
+    {
+      alloc.deallocate(alms);
+    }
+
+    virtual const T* data() const { return alms; }
+    virtual T * data() { return alms;} 
+    virtual long size() const { return size; }
+
+    virtual FourierMap<T> *mimick() const
+    {
+      return new HealpixFourierALM<T>(Lmax(), Mmax());
+    }
+  };
+
+  template<typename T>
+  class HealpixFourierTransform: public FourierTransform<T>
+  {
+  private:
+    HealpixFourierMap<T> realMap;
+    HealpixFourierALM<T> fourierMap;
+    psht_joblist<T> jobs;
+  public:
+    HealpixFourierTransform(long nSide, long Lmax, long Mmax)
+      : realMap(nSide), fourierMap(Lmax, Mmax)
+    {
+      jobs.set_Healpix_geometry(nSide);
+      jobs.set_triangular_alm_info(Lmax, Mmax);
+    }
+
+    virtual ~HealpixFourierTransform() {}
+
+    virtual const FourierMap<std::complex<T> >& fourierSpace() const { return fourierMap; }
+
+    virtual FourierMap<std::complex<T> >& fourierSpace() { return fourierMap; }
+
+    virtual const FourierMap<T>& realSpace() const { return realMap; }
+
+    virtual FourierMap<T>& realSpace() { return realMap; }
+    
+    virtual FourierTransform<T> *mimick() const 
+    {
+      return new HealpixFourierTransform<T>(realMap.Nside(), fourierMap.Lmax(), fourierMap.Mmax());
+    }
+
+    virtual void analysis()
+    {
+      jobs.add_map2alm(realMap.data(), 
+		       reinterpret_cast<xcomplex<T> *>(fourierMap.data()), 
+		       false);
+      jobs.execute();
+      jobs.clear_jobs();
+    }
+
+    virtual void synthesis()
+    {
+      jobs.add_alm2map(reinterpret_cast<xcomplex<T> *>(fourierMap.data()), 
+		       realMap.data(),
+		       false);
+      jobs.execute();
+      jobs.clear_jobs();
+    }
+
+    virtual void analysis_conjugate()
+    {
+      synthesis();
+      realMap.scale(4*M_PI/realMap.Npix());
+    }
+
+    virtual void synthesis_conjugate()
+    {
+      analysis();
+      fourierMap.scale(realMap.Npix()/(4*M_PI));
+    }
+    
+  };
+};
+
+#endif

src/yorick_nc3.cpp

@@ -1,5 +1,11 @@
+#include "ctool_netcdf_ver.hpp"
 #include "config.hpp"
+#ifdef NETCDFCPP4
+#include <netcdf>
+using namespace netCDF
+#else
 #include <netcdfcpp.h>
+#endif
 #include <fstream>
 #include "yorick.hpp"
 #include <assert.h>

src/yorick_nc4.cpp

@@ -0,0 +1,237 @@
+#include "ctool_netcdf_ver.hpp"
+#include "config.hpp"
+#include <netcdf>
+using namespace netCDF;
+#include <fstream>
+#include "yorick.hpp"
+#include <assert.h>
+
+using namespace CosmoTool;
+using namespace std;
+
+class NetCDF_handle
+{
+public:
+  NcFile *outFile;
+  NcVar curVar;
+  vector<size_t> curPos;
+  vector<size_t> counts;
+  vector<NcDim> dimList;
+  uint32_t rank;
+  
+  NetCDF_handle(NcFile *f, NcVar v, vector<NcDim>& dimList, uint32_t rank);
+  virtual ~NetCDF_handle();
+};
+
+NetCDF_handle::NetCDF_handle(NcFile *f, NcVar v, vector<NcDim>& dimList, uint32_t rank)
+{
+  this->outFile = f;
+  this->curVar = v;
+  this->dimList = dimList;
+  this->rank = rank;
+  this->counts.resize(rank);
+  this->curPos.resize(rank);
+
+  for (long i = 0; i < rank; i++)
+    this->curPos[i] = 0;
+
+  for (long i = 0; i < rank; i++)
+    this->counts[i] = 1;
+}
+
+NetCDF_handle::~NetCDF_handle()
+{
+  delete outFile;
+}
+
+template<typename T>
+class InputGenCDF: public NetCDF_handle, public ProgressiveInputImpl<T>
+{
+public:
+  InputGenCDF(NcFile *f, NcVar v, vector<NcDim>& dimList, uint32_t rank)
+    : NetCDF_handle(f,v,dimList,rank)
+  {}
+  virtual ~InputGenCDF() {} 
+
+  virtual T read()
+  {
+    T a;
+
+    curVar.getVar(curPos, counts, &a);
+    
+    curPos[rank-1]++;
+    for (long i = rank-1; i >= 1; i--)
+      {
+	if (curPos[i] == dimList[i].getSize())
+	  {
+	    curPos[i-1]++;
+	  curPos[i] = 0;
+	}
+      }    
+    return a;
+  }
+
+  virtual void seek(uint32_t *pos)
+  {
+    for (long i = rank-1; i >= 0; i--)
+      curPos[i] = pos[rank-1-i];
+  }
+};
+
+template<typename T>
+class OutputGenCDF: public NetCDF_handle, public ProgressiveOutputImpl<T>
+{
+public:
+  OutputGenCDF(NcFile *f, NcVar v, vector<NcDim>& dimList, uint32_t rank)
+    : NetCDF_handle(f,v,dimList,rank)
+  {}
+  virtual ~OutputGenCDF() {} 
+
+  virtual void put(T a)
+  {
+    curVar.putVar(curPos, counts, &a);
+    
+    curPos[rank-1]++;
+    for (long i = rank-1; i >= 1; i--)
+      {
+	if (curPos[i] == dimList[i].getSize())
+	  {
+	    curPos[i-1]++;
+	  curPos[i] = 0;
+	}
+      }    
+  }
+};
+
+template<typename T> NcType& get_NetCDF_type();
+
+#define IMPL_TYPE(T,ncT) \
+template<> \
+NcType& get_NetCDF_type<T>() \
+{ \
+  return ncT; \
+}
+
+IMPL_TYPE(int,ncInt);
+IMPL_TYPE(unsigned int,ncInt);
+IMPL_TYPE(double,ncDouble);
+IMPL_TYPE(float,ncFloat);
+
+namespace CosmoTool {
+  template<typename T>
+  ProgressiveOutput<T>
+  ProgressiveOutput<T>::saveArrayProgressive(const std::string& fname, uint32_t *dimList,
+					     uint32_t rank)
+  {
+    NcFile *f = new NcFile(fname, NcFile::replace);
+    
+    vector<NcDim> dimArray;
+    for (uint32_t i = 0; i < rank; i++)
+      {
+	char dimName[255];
+	
+	sprintf(dimName, "dim%d", i);
+	dimArray.push_back(f->addDim(dimName, dimList[rank-1-i]));
+      }
+    
+    NcVar v = f->addVar("array", get_NetCDF_type<T>(), dimArray);  
+
+    vector<NcDim> ldimList;
+
+    for (uint32_t i = 0; i < rank; i++)
+      ldimList.push_back(dimArray[rank-1-i]);
+
+    OutputGenCDF<T> *impl = new OutputGenCDF<T>(f, v, ldimList, rank);
+    return ProgressiveOutput<T>(impl);   
+  }
+
+  template<typename T>
+  ProgressiveInput<T>
+  ProgressiveInput<T>::loadArrayProgressive(const std::string& fname, uint32_t *&dimList,
+					      uint32_t& rank)
+  {
+    NcFile *f = new NcFile(fname, NcFile::read);
+
+    NcVar v = f->getVar("array");
+
+    rank = v.getDimCount();
+    vector<NcDim> vdimlist = v.getDims();
+    dimList = new uint32_t[rank];
+    for (uint32_t i = 0; i < rank; i++)
+      {
+	dimList[rank-i-1] = vdimlist[i].getSize();
+      }
+    InputGenCDF<T> *impl = new InputGenCDF<T>(f, v, vdimlist, rank);
+
+    return ProgressiveInput<T>(impl);
+  }
+
+  template<typename T>
+  void saveArray(const std::string& fname,
+		 T *array, uint32_t *dimList, uint32_t rank)
+  {
+    NcFile f(fname.c_str(), NcFile::replace);
+    
+    vector<NcDim> dimArray;
+    for (uint32_t i = 0; i < rank; i++)
+      {
+	char dimName[255];
+	
+	sprintf(dimName, "dim%d", i);
+	dimArray.push_back(f.addDim(dimName, dimList[i]));
+      }
+
+    NcVar v = f.addVar("array", get_NetCDF_type<T>(), dimArray);
+    
+    v.putVar(array);
+  }
+
+  template<typename T>
+  void loadArray(const std::string& fname,
+		 T*&array, uint32_t *&dimList, uint32_t& rank)
+	throw (NoSuchFileException)
+  {
+    NcFile f(fname.c_str(), NcFile::read);
+    
+    //if (!f.is_valid())
+    //  throw NoSuchFileException(fname);
+    
+    NcVar v = f.getVar("array"); 
+    vector<NcDim> dims = v.getDims();
+    rank = v.getDimCount();
+    uint32_t fullSize = 1;
+    dimList = new uint32_t[rank];
+    for (int i = 0; i < rank; i++)
+      {
+	dimList[i] = dims[i].getSize();
+	fullSize *= dimList[i];
+      }
+    if (fullSize != 0) {
+      array = new T[fullSize];
+      v.getVar(array);
+    }
+  }
+
+  template class ProgressiveInput<int>;
+  template class ProgressiveInput<float>;
+  template class ProgressiveInput<double>;
+
+  template class ProgressiveOutput<int>;
+  template class ProgressiveOutput<float>;
+  template class ProgressiveOutput<double>;
+
+  template void loadArray<int>(const std::string& fname,
+			       int*& array, uint32_t *&dimList, uint32_t& rank);
+  template void loadArray<float>(const std::string& fname,
+				 float*& array, uint32_t *&dimList, uint32_t& rank);
+  template void loadArray<double>(const std::string& fname,
+				  double*& array, uint32_t *&dimList, uint32_t& rank);
+
+  template void saveArray<int>(const std::string& fname,
+			       int *array, uint32_t *dimList, uint32_t rank);
+  template void saveArray<float>(const std::string& fname,
+				 float *array, uint32_t *dimList, uint32_t rank);
+  template void saveArray<double>(const std::string& fname,
+				  double *array, uint32_t *dimList, uint32_t rank);
+  
+}