diff options
author | Karl Berry <karl@freefriends.org> | 2013-04-02 23:03:22 +0000 |
---|---|---|
committer | Karl Berry <karl@freefriends.org> | 2013-04-02 23:03:22 +0000 |
commit | 324efbd652f973cd70f42fe08c665c518995bcbe (patch) | |
tree | 818ae382dfa133f3dcceec8e12e4cbd42be8bcf4 /Build/source/utils/asymptote/fftw++.h | |
parent | bbf699abe485e707133ff3ba502119f261f27269 (diff) |
asy 2.21 sources (no runtime update)
git-svn-id: svn://tug.org/texlive/trunk@29618 c570f23f-e606-0410-a88d-b1316a301751
Diffstat (limited to 'Build/source/utils/asymptote/fftw++.h')
-rw-r--r-- | Build/source/utils/asymptote/fftw++.h | 136 |
1 files changed, 118 insertions, 18 deletions
diff --git a/Build/source/utils/asymptote/fftw++.h b/Build/source/utils/asymptote/fftw++.h index 0ea30502b87..be68c899290 100644 --- a/Build/source/utils/asymptote/fftw++.h +++ b/Build/source/utils/asymptote/fftw++.h @@ -1,5 +1,5 @@ /* Fast Fourier transform C++ header class for the FFTW3 Library - Copyright (C) 2004-10 John C. Bowman, University of Alberta + Copyright (C) 2004-12 John C. Bowman, University of Alberta This program is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by @@ -18,7 +18,7 @@ #ifndef __fftwpp_h__ #define __fftwpp_h__ 1 -#define __FFTWPP_H_VERSION__ 1.08 +#define __FFTWPP_H_VERSION__ 1.12 #include <cstdlib> #include <fstream> @@ -26,11 +26,35 @@ #include <fftw3.h> #include <cerrno> +#ifndef FFTWPP_SINGLE_THREAD +#include <omp.h> +#endif + +inline int get_max_threads() +{ +#ifdef FFTWPP_SINGLE_THREAD + return 1; +#else + return omp_get_max_threads(); +#endif +} + +inline int get_thread_num() +{ +#ifdef FFTWPP_SINGLE_THREAD + return 0; +#else + return omp_get_thread_num(); +#endif +} + #ifndef __Complex_h__ #include <complex> typedef std::complex<double> Complex; #endif +#include "seconds.h" + #ifndef HAVE_POSIX_MEMALIGN #ifdef __GLIBC_PREREQ @@ -167,7 +191,7 @@ inline int GetWisdom(std::ifstream& s) {return s.get();} // class fftw { protected: - unsigned int doubles; // number of double words in dataset + unsigned int doubles; // number of double precision values in dataset int sign; double norm; @@ -194,8 +218,14 @@ protected: static std::ofstream ofWisdom; static bool Wise; static const double twopi; + unsigned int threads; public: + static unsigned int maxthreads; + static double testseconds; + + unsigned int Threads() {return threads;} + // Shift the Fourier origin to (nx/2,0). static void Shift(Complex *data, unsigned int nx, unsigned int ny, int sign=0) { @@ -271,8 +301,15 @@ public: fftw(unsigned int doubles, int sign, unsigned int n=0) : doubles(doubles), sign(sign), norm(1.0/(n ? n : (doubles+1)/2)), - plan(NULL) - {} + plan(NULL) { + static bool initialize=true; + if(initialize) { +#ifndef FFTWPP_SINGLE_THREAD + fftw_init_threads(); +#endif + initialize=false; + } + } virtual ~fftw() {if(plan) fftw_destroy_plan(plan);} @@ -284,6 +321,15 @@ public: << "-byte aligned: address " << p << std::endl; } + double stdev(unsigned int N, double sum, double sum2) { + return N > 1 ? sqrt((sum2-sum*sum/N)/(N-1)) : 0.0; + } + + void noplan() { + std::cerr << "Unable to construct FFTW plan" << std::endl; + exit(1); + } + void Setup(Complex *in, Complex *out=NULL) { if(!Wise) LoadWisdom(); bool alloc=!in; @@ -296,11 +342,62 @@ public: if(!out) out=in; #endif inplace=(out==in); - plan=Plan(in,out); - if(!plan) { - std::cerr << "Unable to construct FFTW plan" << std::endl; - exit(1); - } + +#ifndef FFTWPP_SINGLE_THREAD + fftw_plan_with_nthreads(1); +#endif + fftw_plan plan1=Plan(in,out); + if(!plan1) noplan(); + plan=plan1; + + if(maxthreads > 1) { + double sum=0.0; + double sum2=0.0; + unsigned int N=1; + const unsigned int microseconds=1000000; + unsigned int limit=(int) (testseconds*microseconds); + for(; N < limit; ++N) { + seconds(); + fft(in,out); + double t=seconds(); + sum += t; + sum2 += t*t; + if(sum > testseconds) + break; + } + double mean1=sum/N; + double stdev1=stdev(N,sum,sum2); + + threads=get_max_threads(); + if(maxthreads > threads) maxthreads=threads; +#ifndef FFTWPP_SINGLE_THREAD + fftw_plan_with_nthreads(maxthreads); +#endif + plan=Plan(in,out); + if(!plan) noplan(); + + if(plan) { + sum=0.0; + seconds(); + for(unsigned int i=1; i <= N; ++i) { + seconds(); + fft(in,out); + double t=seconds(); + sum += t; + } + double mean2=sum/N; + + if(mean2 > mean1-stdev1) { + threads=1; + fftw_destroy_plan(plan); + plan=plan1; + } else { + fftw_destroy_plan(plan1); + threads=maxthreads; + } + } + } else + threads=1; if(alloc) Array::deleteAlign(in,(doubles+1)/2); SaveWisdom(); @@ -416,6 +513,9 @@ public: else { out=Setout(in,out); Execute(in,out); +#ifndef FFTWPP_SINGLE_THREAD +#pragma omp parallel for num_threads(threads) +#endif for(unsigned int k=0; k < M; k++) { for(unsigned int j=0; j < nx; j++) { out[j*stride+k*dist] *= norm; @@ -500,8 +600,8 @@ public: {Setup(in,out);} fftw_plan Plan(Complex *in, Complex *out) { - int n[1]={nx}; - return fftw_plan_many_dft(1,n,M, + int n=(int) nx; + return fftw_plan_many_dft(1,&n,M, (fftw_complex *) in,NULL,stride,dist, (fftw_complex *) out,NULL,stride,dist, sign,effort); @@ -648,8 +748,8 @@ public: stride(stride), dist(Dist(nx,stride,dist)) {Setup(in,out);} fftw_plan Plan(Complex *in, Complex *out) { - const int n[1]={nx}; - return fftw_plan_many_dft_r2c(1,n,M, + int n=(int) nx; + return fftw_plan_many_dft_r2c(1,&n,M, (double *) in,NULL,stride,2*dist, (fftw_complex *) out,NULL,stride,dist, effort); @@ -699,8 +799,8 @@ public: nx(nx), M(M), stride(stride), dist(Dist(nx,stride,dist)) {Setup(in,out);} fftw_plan Plan(Complex *in, Complex *out) { - const int n[1]={nx}; - return fftw_plan_many_dft_c2r(1,n,M, + int n=(int) nx; + return fftw_plan_many_dft_c2r(1,&n,M, (fftw_complex *) in,NULL,stride,dist, (double *) out,NULL,stride,2*dist, effort); @@ -1019,11 +1119,11 @@ public: #ifdef __Array_h__ crfft3d(unsigned int nz, const Array::array3<double>& out) : fftw(out.Size(),1,out.Nx()*out.Ny()*nz), - nx(out.Nx()), ny(out.Ny()), nz(nz) {Setup(out);} + nx(out.Nx()), ny(out.Ny()), nz(nz) {Setup(out);} crfft3d(unsigned int nz, const Array::array3<Complex>& in) : fftw(2*in.Size(),1,in.Nx()*in.Ny()*nz), - nx(in.Nx()), ny(in.Ny()), nz(nz) {Setup(in);} + nx(in.Nx()), ny(in.Ny()), nz(nz) {Setup(in);} crfft3d(unsigned int nz, const Array::array3<Complex>& in, const Array::array3<double>& out) |