From 7ad49118728dc5e98124a86f6cfeffa001ccde81 Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Sun, 11 Nov 2012 15:46:03 -0500 Subject: [PATCH] Changed interface of newRandomFourier to avoid generating copies. --- sample/CMakeLists.txt | 2 +- sample/test_fft_calls.cpp | 3 +-- src/fourier/base_types.hpp | 4 ++-- src/fourier/euclidian.hpp | 33 ++++++++++++++------------------- src/fourier/healpix.hpp | 12 ++++-------- 5 files changed, 22 insertions(+), 32 deletions(-) diff --git a/sample/CMakeLists.txt b/sample/CMakeLists.txt index 9d65cc3..17e01c4 100644 --- a/sample/CMakeLists.txt +++ b/sample/CMakeLists.txt @@ -47,7 +47,7 @@ target_link_libraries(testBSP ${tolink}) if (FFTW3_FOUND AND FFTW3F_FOUND AND EIGEN3_FOUND) add_executable(test_fft_calls test_fft_calls.cpp) target_link_libraries(test_fft_calls ${tolink} ${FFTW3_LIBRARIES} ${FFTW3F_LIBRARIES}) -endif (FFTW3_FOUND AND EIGEN3_FOUND) +endif (FFTW3_FOUND AND FFTW3F_FOUND AND EIGEN3_FOUND) if (SHARP_LIBRARY AND SHARP_INCLUDE_PATH AND EIGEN3_FOUND) include_directories(${SHARP_INCLUDE_PATH}) diff --git a/sample/test_fft_calls.cpp b/sample/test_fft_calls.cpp index 4ba1f91..f62e678 100644 --- a/sample/test_fft_calls.cpp +++ b/sample/test_fft_calls.cpp @@ -37,8 +37,7 @@ void test_2d(int Nx, int Ny) dft.realSpace().scale(2.0); dft.fourierSpace().scale(0.2); - typename SpectrumFunction::FourierMapPtr m = spectrum.newRandomFourier(rng, dft.fourierSpace()); - dft.fourierSpace() = *m.get(); + spectrum.newRandomFourier(rng, dft.fourierSpace()); dft.synthesis(); uint32_t dims[2] = { Ny, Nx }; diff --git a/src/fourier/base_types.hpp b/src/fourier/base_types.hpp index 29a4c87..73d3698 100644 --- a/src/fourier/base_types.hpp +++ b/src/fourier/base_types.hpp @@ -26,8 +26,8 @@ namespace CosmoTool typedef FourierMap > FourierMapType; typedef boost::shared_ptr FourierMapPtr; - virtual FourierMapPtr - newRandomFourier(gsl_rng *rng, const FourierMapType& like_map) const = 0; + virtual void + newRandomFourier(gsl_rng *rng, FourierMapType& out_map) const = 0; virtual void mul(FourierMapType& m) const = 0; virtual void mul_sqrt(FourierMapType& m) const = 0; diff --git a/src/fourier/euclidian.hpp b/src/fourier/euclidian.hpp index 600ebbb..fb92e40 100644 --- a/src/fourier/euclidian.hpp +++ b/src/fourier/euclidian.hpp @@ -27,12 +27,12 @@ namespace CosmoTool { } - ptr_map newRandomFourier(gsl_rng *rng, const FourierMapType& like_map) const; + void newRandomFourier(gsl_rng *rng, FourierMapType& out_map) const; - void mul(FourierMap >& m) const; - void mul_sqrt(FourierMap >& m) const; - void mul_inv(FourierMap >& m) const; - void mul_inv_sqrt(FourierMap >& m) const; + void mul(FourierMapType& m) const; + void mul_sqrt(FourierMapType& m) const; + void mul_inv(FourierMapType& m) const; + void mul_inv_sqrt(FourierMapType& m) const; }; template @@ -388,16 +388,13 @@ namespace CosmoTool template - typename EuclidianSpectrum_1D::ptr_map - EuclidianSpectrum_1D::newRandomFourier(gsl_rng *rng, const FourierMapType& like_map) const + void EuclidianSpectrum_1D::newRandomFourier(gsl_rng *rng, FourierMapType& out_map) const { typedef EuclidianFourierMapComplex MapT; typedef typename EuclidianSpectrum_1D::ptr_map ptr_map; typedef typename MapT::DimArray DimArray; - const MapT& m_c = dynamic_cast(like_map); - ptr_map ret_map = ptr_map(m_c.mimick()); - MapT& rand_map = dynamic_cast(*ret_map.get()); + MapT& rand_map = dynamic_cast(out_map); std::complex *d = rand_map.data(); long idx; @@ -405,7 +402,7 @@ namespace CosmoTool long plane_size; bool alleven = rand_map.allDimensionsEven(); - for (long p = 1; p < m_c.size(); p++) + for (long p = 1; p < rand_map.size(); p++) { double A_k = std::sqrt(0.5*f(rand_map.get_K(p))); d[p] = std::complex(gsl_ran_gaussian(rng, A_k), @@ -415,14 +412,14 @@ namespace CosmoTool d[0] = std::complex(gsl_ran_gaussian(rng, std::sqrt(f(0))), 0); if (!rand_map.firstDimensionEven()) - return ret_map; + return; // Correct the Nyquist plane idx = dims[0]-1; // Stick to the last element of the first dimension d[idx] = std::complex(d[idx].real() + d[idx].imag(), 0); // 1D is special case if (dims.size() == 1) - return ret_map; + return; plane_size = 1; for (int q = 1; q < dims.size(); q++) @@ -447,12 +444,10 @@ namespace CosmoTool q += dims[0]-1; d[q] = std::complex(d[q].real()+d[q].imag(),0); } - - return ret_map; } template - void EuclidianSpectrum_1D::mul(FourierMap >& m) const + void EuclidianSpectrum_1D::mul(FourierMapType& m) const { EuclidianFourierMapComplex& m_c = dynamic_cast&>(m); std::complex *d = m.data(); @@ -462,7 +457,7 @@ namespace CosmoTool } template - void EuclidianSpectrum_1D::mul_sqrt(FourierMap >& m) const + void EuclidianSpectrum_1D::mul_sqrt(FourierMapType& m) const { EuclidianFourierMapComplex& m_c = dynamic_cast&>(m); std::complex *d = m.data(); @@ -472,7 +467,7 @@ namespace CosmoTool } template - void EuclidianSpectrum_1D::mul_inv(FourierMap >& m) const + void EuclidianSpectrum_1D::mul_inv(FourierMapType& m) const { EuclidianFourierMapComplex& m_c = dynamic_cast&>(m); std::complex *d = m.data(); @@ -488,7 +483,7 @@ namespace CosmoTool } template - void EuclidianSpectrum_1D::mul_inv_sqrt(FourierMap >& m) const + void EuclidianSpectrum_1D::mul_inv_sqrt(FourierMapType& m) const { EuclidianFourierMapComplex& m_c = dynamic_cast&>(m); std::complex *d = m.data(); diff --git a/src/fourier/healpix.hpp b/src/fourier/healpix.hpp index 06da0fc..bd868e7 100644 --- a/src/fourier/healpix.hpp +++ b/src/fourier/healpix.hpp @@ -32,7 +32,7 @@ namespace CosmoTool const T *data() const { return &cls[0]; } long size() const { return cls.size(); } - ptr_map newRandomFourier(gsl_rng *rng, const FourierMapType& like_map) const; + void newRandomFourier(gsl_rng *rng, FourierMapType& like_map) const; void mul(FourierMapType& m) const; void mul_sqrt(FourierMapType& m) const; @@ -244,14 +244,11 @@ namespace CosmoTool }; template - typename HealpixSpectrum::ptr_map - HealpixSpectrum::newRandomFourier(gsl_rng *rng, const FourierMapType& like_map) const + void HealpixSpectrum::newRandomFourier(gsl_rng *rng, FourierMapType& out_map) const { - const HealpixFourierALM& alms = dynamic_cast&>(like_map); - HealpixFourierALM *new_alms; - ptr_map r(new_alms = new HealpixFourierALM(alms.Lmax(), alms.Mmax())); + HealpixFourierALM& alms = dynamic_cast&>(out_map); long lmaxGen = std::min(cls.size()-1, alms.Lmax()); - std::complex *new_data = new_alms->data(); + std::complex *new_data = alms.data(); for (long l = 0; l <= lmaxGen; l++) { @@ -266,7 +263,6 @@ namespace CosmoTool c.imag() = gsl_ran_gaussian(rng, Al); } } - return r; } template