From f408f07c05847ece1b1fea58d1ef1987ccbd22d2 Mon Sep 17 00:00:00 2001 From: Jens Keiner Date: Fri, 4 Mar 2016 21:30:26 +0100 Subject: [PATCH 01/11] #2: Ignore some more files. --- .gitignore | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/.gitignore b/.gitignore index ca52170d..e6b73203 100644 --- a/.gitignore +++ b/.gitignore @@ -110,3 +110,7 @@ examples/nsfft/simple_test examples/solver/glacier examples/solver/simple_test + +CUnitAutomated-Results.xml + +tests/checkall From 39721fdd47e4e953e06081d39cb5187f865af76a Mon Sep 17 00:00:00 2001 From: Jens Keiner Date: Fri, 4 Mar 2016 21:44:28 +0100 Subject: [PATCH 02/11] #2: Remove debug code. --- tests/Makefile.am | 4 ---- 1 file changed, 4 deletions(-) diff --git a/tests/Makefile.am b/tests/Makefile.am index fcc391b5..bb18ae3c 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -19,7 +19,6 @@ else endif check_PROGRAMS = $(CHECK) $(CHECK_THREADS) -#debug TESTS = $(check_PROGRAMS) @@ -46,8 +45,5 @@ if HAVE_OPENMP endif endif -#debug_SOURCES = debug.c nfct.c nfst.c nfct.h nfst.h -#debug_LDADD = $(top_builddir)/libnfft3@PREC_SUFFIX@.la -lcunit -lncurses - clean-local: rm -f CUnitAutomated-Results.xml CUnitAutomated_threads-Results.xml checkall.log checkall.trs checkall_threads.log checkall_threads.trs From 5d20b5994fc4d4e62bcc58fd0c2ce2de9a5be1c8 Mon Sep 17 00:00:00 2001 From: Jens Keiner Date: Fri, 4 Mar 2016 21:44:51 +0100 Subject: [PATCH 03/11] #2: Add openMP to Travis CI configuration. --- .travis.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.travis.yml b/.travis.yml index 71e53be2..b22958a8 100644 --- a/.travis.yml +++ b/.travis.yml @@ -28,7 +28,7 @@ env: - WINDOW=sinc PRECISION=--enable-long-double # Compile, link, and run tests. -script: ./bootstrap.sh && ./configure --with-window=$WINDOW $PRECISION --enable-all && make && make check +script: ./bootstrap.sh && ./configure --with-window=$WINDOW $PRECISION --enable-all --enable-openmp && make && make check ## Print config.log for debugging. after_failure: "cat config.log" From f434d2a5e928087e92fcde23ee20e04134e96b63 Mon Sep 17 00:00:00 2001 From: Jens Keiner Date: Fri, 4 Mar 2016 22:36:05 +0100 Subject: [PATCH 04/11] #2: Disable OpenMP when using clang. --- .travis.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.travis.yml b/.travis.yml index b22958a8..fe4d61ee 100644 --- a/.travis.yml +++ b/.travis.yml @@ -28,7 +28,7 @@ env: - WINDOW=sinc PRECISION=--enable-long-double # Compile, link, and run tests. -script: ./bootstrap.sh && ./configure --with-window=$WINDOW $PRECISION --enable-all --enable-openmp && make && make check +script: ./bootstrap.sh && ./configure --with-window=$WINDOW $PRECISION --enable-all $(if test "$CC" = "clang"; then echo ""; else echo "--enable-openmp"; fi) && make && make check ## Print config.log for debugging. after_failure: "cat config.log" From 008265d3e0cee7a6322056e85d1f245db1c7f1fb Mon Sep 17 00:00:00 2001 From: Toni Volkmer Date: Mon, 7 Mar 2016 14:28:52 +0100 Subject: [PATCH 05/11] removed temporary array for bspline function from function header (R *scratch) and from nfft_plan (R * spline_coeffs), the temporary array is now created within the bspline function on the stack --- applications/iterS2/iterS2.c | 5 +---- include/infft.h | 20 +++++++------------- include/nfft3.h | 6 ------ kernel/mri/mri.c | 2 -- kernel/util/bspline.c | 3 ++- tests/bspline.c | 3 +-- 6 files changed, 11 insertions(+), 28 deletions(-) diff --git a/applications/iterS2/iterS2.c b/applications/iterS2/iterS2.c index 0b12f3f2..51aa80b1 100644 --- a/applications/iterS2/iterS2.c +++ b/applications/iterS2/iterS2.c @@ -200,7 +200,6 @@ int main (int argc, char **argv) double re; double im; double a; - double *scratch; double xs; double *ys; double *temp; @@ -284,7 +283,6 @@ int main (int argc, char **argv) X(next_power_of_2_exp)(N, &npt, &npt_exp); fprintf(stderr, "npt = %d, npt_exp = %d\n", npt, npt_exp); fprintf(stderr,"Optimal interpolation!\n"); - scratch = (double*) nfft_malloc(4*sizeof(double)); ys = (double*) nfft_malloc((N+1)*sizeof(double)); temp = (double*) nfft_malloc((2*N+1)*sizeof(double)); temp2 = (double _Complex*) nfft_malloc((N+1)*sizeof(double _Complex)); @@ -293,7 +291,7 @@ int main (int argc, char **argv) for (j = 0; j <= N; j++) { xs = 2.0 + (2.0*j)/(N+1); - ys[j] = (2.0-((j == 0)?(1.0):(0.0)))*4.0*nfft_bspline(4,xs,scratch); + ys[j] = (2.0-((j == 0)?(1.0):(0.0)))*4.0*nfft_bspline(4,xs); //fprintf(stdout,"%3d: g(%le) = %le\n",j,xs,ys[j]); a += ys[j]; } @@ -369,7 +367,6 @@ int main (int argc, char **argv) fftw_destroy_plan(fplan); - nfft_free(scratch); nfft_free(qweights); nfft_free(ys); nfft_free(temp); diff --git a/include/infft.h b/include/infft.h index 1a5c2941..a8739150 100644 --- a/include/infft.h +++ b/include/infft.h @@ -187,12 +187,9 @@ typedef ptrdiff_t INT; POW(SIN((k) * KPI / n) / ((k) * KPI / n), \ K(2.0) * ths->m)/n)) #define PHI(n,x,d) (Y(bspline)(2*ths->m,((x)*n) + \ - (R)ths->m,ths->spline_coeffs) / n) - #define WINDOW_HELP_INIT \ - { \ - ths->spline_coeffs= (R*)Y(malloc)(2*ths->m*sizeof(R)); \ - } - #define WINDOW_HELP_FINALIZE {Y(free)(ths->spline_coeffs);} + (R)ths->m) / n) + #define WINDOW_HELP_INIT + #define WINDOW_HELP_FINALIZE #if defined(NFFT_LDOUBLE) #define WINDOW_HELP_ESTIMATE_m 11 #elif defined(NFFT_SINGLE) @@ -203,17 +200,14 @@ typedef ptrdiff_t INT; #elif defined(SINC_POWER) #define PHI_HUT(n,k,d) (Y(bspline)(2 * ths->m, (K(2.0) * ths->m*(k)) / \ ((K(2.0) * ths->sigma[(d)] - 1) * n / \ - ths->sigma[(d)]) + (R)ths->m, ths->spline_coeffs)) + ths->sigma[(d)]) + (R)ths->m)) #define PHI(n,x,d) ((R)(n / ths->sigma[(d)] * \ (K(2.0) * ths->sigma[(d)] - K(1.0))/ (K(2.0)*ths->m) * \ POW(Y(sinc)(KPI * n / ths->sigma[(d)] * (x) * \ (K(2.0) * ths->sigma[(d)] - K(1.0)) / (K(2.0)*ths->m)) , 2*ths->m) / \ n)) - #define WINDOW_HELP_INIT \ - { \ - ths->spline_coeffs= (R*)Y(malloc)(2 * ths->m * sizeof(R)); \ - } - #define WINDOW_HELP_FINALIZE {Y(free)(ths->spline_coeffs);} + #define WINDOW_HELP_INIT + #define WINDOW_HELP_FINALIZE #if defined(NFFT_LDOUBLE) #define WINDOW_HELP_ESTIMATE_m 13 #elif defined(NFFT_SINGLE) @@ -1439,7 +1433,7 @@ R Y(lambda2)(R mu, R nu); R Y(bessel_i0)(R x); /* bspline.c: */ -R Y(bspline)(const INT, const R x, R*); +R Y(bspline)(const INT, const R x); /* float.c: */ typedef enum {NFFT_EPSILON = 0, NFFT_SAFE__MIN = 1, NFFT_BASE = 2, diff --git a/include/nfft3.h b/include/nfft3.h index eade8696..2c976194 100644 --- a/include/nfft3.h +++ b/include/nfft3.h @@ -154,8 +154,6 @@ typedef struct\ C *g_hat; /**< Zero-padded vector of Fourier coefficients, size is \ref n_total fftw_complex */\ C *g1; /**< Input of fftw */\ C *g2; /**< Output of fftw */\ -\ - R *spline_coeffs; /**< Input for de Boor algorithm if B_SPLINE or SINC_POWER is defined */\ \ NFFT_INT *index_x; /**< Index array for nodes x used when flag \ref NFFT_SORT_NODES is set. */\ } X(plan); \ @@ -261,7 +259,6 @@ typedef struct\ R *g1; /**< input of fftw */\ R *g2; /**< output of fftw */\ \ - R *spline_coeffs; /**< input for de Boor algorithm, if B_SPLINE or SINC_2m is defined */\ } X(plan);\ \ NFFT_EXTERN void X(init_1d)(X(plan) *ths_plan, int N0, int M_total); \ @@ -340,8 +337,6 @@ typedef struct\ R *g_hat;\ R *g1; /**< input of fftw */\ R *g2; /**< output of fftw */\ -\ - R *spline_coeffs; /**< input for de Boor algorithm, if B_SPLINE or SINC_2m is defined */\ \ R X(full_psi_eps);\ } X(plan);\ @@ -412,7 +407,6 @@ typedef struct\ int *psi_index_g; /**< only for thin B */\ int *psi_index_f; /**< only for thin B */\ C *F;\ - R *spline_coeffs; /**< input for de Boor algorithm, if B_SPLINE or SINC_2m is defined */\ } X(plan);\ \ NFFT_EXTERN void X(init)(X(plan) *ths_plan, int d, int N_total, int M_total, int *N); \ diff --git a/kernel/mri/mri.c b/kernel/mri/mri.c index 0740ff3f..15218b82 100644 --- a/kernel/mri/mri.c +++ b/kernel/mri/mri.c @@ -37,8 +37,6 @@ typedef struct window_funct_plan_ { int n[1]; double sigma[1]; double *b; - double *spline_coeffs; /**< input for de Boor algorithm, if - B_SPLINE or SINC_2m is defined */ } window_funct_plan; /** diff --git a/kernel/util/bspline.c b/kernel/util/bspline.c index 5cc88939..f12b730e 100644 --- a/kernel/util/bspline.c +++ b/kernel/util/bspline.c @@ -34,7 +34,7 @@ static inline void bspline_help(const INT k, const R x, R *scratch, const INT j, } /* Evaluate the cardinal B-Spline B_{n-1} supported on [0,n]. */ -R Y(bspline)(const INT k, const R _x, R *scratch) +R Y(bspline)(const INT k, const R _x) { const R kk = (R)k; R result_value; @@ -43,6 +43,7 @@ R Y(bspline)(const INT k, const R _x, R *scratch) INT j, idx, ug, og; /* indices */ R a; /* Alpha of de Boor scheme*/ R x = _x; + R scratch[k]; result_value = K(0.0); diff --git a/tests/bspline.c b/tests/bspline.c index a185351f..52bd45db 100644 --- a/tests/bspline.c +++ b/tests/bspline.c @@ -6333,14 +6333,13 @@ static const R bound = K(16.0) * EPSILON; static int check_bspline(const unsigned n, const unsigned int m, const R *r) { unsigned int j; - R scratch[100]; R err = K(0.0); int ok; for (j = 0; j < m; j++) { const R x = r[2*j], yr = r[2*j+1]; - R y = X(bspline)((INT)(n + 1), x, scratch); + R y = X(bspline)((INT)(n + 1), x); /*printf("x = " __FE__ ", err = " __FE__ "\n", x, ERR(y,yr));*/ err = FMAX(err, ERR(y, yr)); } From 1ea53b0b9f5cfa90fab054fbd996e2b84f3fa906 Mon Sep 17 00:00:00 2001 From: Toni Volkmer Date: Mon, 7 Mar 2016 17:49:15 +0100 Subject: [PATCH 06/11] replaced creal call by macro CREAL and cimag by CIMAG for multi-precision in nfft kernel --- kernel/nfft/nfft.c | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/kernel/nfft/nfft.c b/kernel/nfft/nfft.c index 28fa8108..d5d7fe64 100644 --- a/kernel/nfft/nfft.c +++ b/kernel/nfft/nfft.c @@ -1440,10 +1440,10 @@ static void nfft_adjoint_B_compute_full_psi(C *g, const INT *psi_index_g, R *gref_real = (R*) gref; #pragma omp atomic - gref_real[0] += creal(val); + gref_real[0] += CREAL(val); #pragma omp atomic - gref_real[1] += cimag(val); + gref_real[1] += CIMAG(val); #else g[psi_index_g[j * lprod + l]] += psi[j * lprod + l] * f[j]; #endif @@ -1881,10 +1881,10 @@ static void nfft_adjoint_1d_compute_omp_atomic(const C f, C *g, R *lhs_real = (R*)lhs; C val = psij_const[l] * f; #pragma omp atomic - lhs_real[0] += creal(val); + lhs_real[0] += CREAL(val); #pragma omp atomic - lhs_real[1] += cimag(val); + lhs_real[1] += CIMAG(val); } } #endif @@ -2726,10 +2726,10 @@ static void nfft_adjoint_2d_compute_omp_atomic(const C f, C *g, C val = psij_const0[l0] * psij_const1[l1] * f; #pragma omp atomic - lhs_real[0] += creal(val); + lhs_real[0] += CREAL(val); #pragma omp atomic - lhs_real[1] += cimag(val); + lhs_real[1] += CIMAG(val); } } } @@ -4132,10 +4132,10 @@ static void nfft_adjoint_3d_compute_omp_atomic(const C f, C *g, C val = psij_const0[l0] * psij_const1[l1] * psij_const2[l2] * f; #pragma omp atomic - lhs_real[0] += creal(val); + lhs_real[0] += CREAL(val); #pragma omp atomic - lhs_real[1] += cimag(val); + lhs_real[1] += CIMAG(val); } } } From f751d1566df31fdb59e9468bbeef7b1767bf81b8 Mon Sep 17 00:00:00 2001 From: Toni Volkmer Date: Mon, 7 Mar 2016 23:04:38 +0100 Subject: [PATCH 07/11] added plan field spline_coeffs again for compability --- include/nfft3.h | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/include/nfft3.h b/include/nfft3.h index 2c976194..eade8696 100644 --- a/include/nfft3.h +++ b/include/nfft3.h @@ -154,6 +154,8 @@ typedef struct\ C *g_hat; /**< Zero-padded vector of Fourier coefficients, size is \ref n_total fftw_complex */\ C *g1; /**< Input of fftw */\ C *g2; /**< Output of fftw */\ +\ + R *spline_coeffs; /**< Input for de Boor algorithm if B_SPLINE or SINC_POWER is defined */\ \ NFFT_INT *index_x; /**< Index array for nodes x used when flag \ref NFFT_SORT_NODES is set. */\ } X(plan); \ @@ -259,6 +261,7 @@ typedef struct\ R *g1; /**< input of fftw */\ R *g2; /**< output of fftw */\ \ + R *spline_coeffs; /**< input for de Boor algorithm, if B_SPLINE or SINC_2m is defined */\ } X(plan);\ \ NFFT_EXTERN void X(init_1d)(X(plan) *ths_plan, int N0, int M_total); \ @@ -337,6 +340,8 @@ typedef struct\ R *g_hat;\ R *g1; /**< input of fftw */\ R *g2; /**< output of fftw */\ +\ + R *spline_coeffs; /**< input for de Boor algorithm, if B_SPLINE or SINC_2m is defined */\ \ R X(full_psi_eps);\ } X(plan);\ @@ -407,6 +412,7 @@ typedef struct\ int *psi_index_g; /**< only for thin B */\ int *psi_index_f; /**< only for thin B */\ C *F;\ + R *spline_coeffs; /**< input for de Boor algorithm, if B_SPLINE or SINC_2m is defined */\ } X(plan);\ \ NFFT_EXTERN void X(init)(X(plan) *ths_plan, int d, int N_total, int M_total, int *N); \ From 84ff14b572019819ef592b39d022d61e984bc9b5 Mon Sep 17 00:00:00 2001 From: Toni Volkmer Date: Tue, 8 Mar 2016 11:39:48 +0100 Subject: [PATCH 08/11] renamed function nfft_bspline to nfft_bsplines due to removal of argument R* scratch --- applications/iterS2/iterS2.c | 2 +- include/infft.h | 6 +++--- kernel/util/bspline.c | 2 +- tests/bspline.c | 2 +- 4 files changed, 6 insertions(+), 6 deletions(-) diff --git a/applications/iterS2/iterS2.c b/applications/iterS2/iterS2.c index 51aa80b1..18c686fd 100644 --- a/applications/iterS2/iterS2.c +++ b/applications/iterS2/iterS2.c @@ -291,7 +291,7 @@ int main (int argc, char **argv) for (j = 0; j <= N; j++) { xs = 2.0 + (2.0*j)/(N+1); - ys[j] = (2.0-((j == 0)?(1.0):(0.0)))*4.0*nfft_bspline(4,xs); + ys[j] = (2.0-((j == 0)?(1.0):(0.0)))*4.0*nfft_bsplines(4,xs); //fprintf(stdout,"%3d: g(%le) = %le\n",j,xs,ys[j]); a += ys[j]; } diff --git a/include/infft.h b/include/infft.h index a8739150..26612458 100644 --- a/include/infft.h +++ b/include/infft.h @@ -186,7 +186,7 @@ typedef ptrdiff_t INT; #define PHI_HUT(n,k,d) ((R)(((k) == 0) ? K(1.0) / n : \ POW(SIN((k) * KPI / n) / ((k) * KPI / n), \ K(2.0) * ths->m)/n)) - #define PHI(n,x,d) (Y(bspline)(2*ths->m,((x)*n) + \ + #define PHI(n,x,d) (Y(bsplines)(2*ths->m,((x)*n) + \ (R)ths->m) / n) #define WINDOW_HELP_INIT #define WINDOW_HELP_FINALIZE @@ -198,7 +198,7 @@ typedef ptrdiff_t INT; #define WINDOW_HELP_ESTIMATE_m 11 #endif #elif defined(SINC_POWER) - #define PHI_HUT(n,k,d) (Y(bspline)(2 * ths->m, (K(2.0) * ths->m*(k)) / \ + #define PHI_HUT(n,k,d) (Y(bsplines)(2 * ths->m, (K(2.0) * ths->m*(k)) / \ ((K(2.0) * ths->sigma[(d)] - 1) * n / \ ths->sigma[(d)]) + (R)ths->m)) #define PHI(n,x,d) ((R)(n / ths->sigma[(d)] * \ @@ -1433,7 +1433,7 @@ R Y(lambda2)(R mu, R nu); R Y(bessel_i0)(R x); /* bspline.c: */ -R Y(bspline)(const INT, const R x); +R Y(bsplines)(const INT, const R x); /* float.c: */ typedef enum {NFFT_EPSILON = 0, NFFT_SAFE__MIN = 1, NFFT_BASE = 2, diff --git a/kernel/util/bspline.c b/kernel/util/bspline.c index f12b730e..77c5774d 100644 --- a/kernel/util/bspline.c +++ b/kernel/util/bspline.c @@ -34,7 +34,7 @@ static inline void bspline_help(const INT k, const R x, R *scratch, const INT j, } /* Evaluate the cardinal B-Spline B_{n-1} supported on [0,n]. */ -R Y(bspline)(const INT k, const R _x) +R Y(bsplines)(const INT k, const R _x) { const R kk = (R)k; R result_value; diff --git a/tests/bspline.c b/tests/bspline.c index 52bd45db..ee012ca8 100644 --- a/tests/bspline.c +++ b/tests/bspline.c @@ -6339,7 +6339,7 @@ static int check_bspline(const unsigned n, const unsigned int m, const R *r) for (j = 0; j < m; j++) { const R x = r[2*j], yr = r[2*j+1]; - R y = X(bspline)((INT)(n + 1), x); + R y = X(bsplines)((INT)(n + 1), x); /*printf("x = " __FE__ ", err = " __FE__ "\n", x, ERR(y,yr));*/ err = FMAX(err, ERR(y, yr)); } From 9da9d1387753f8151271a2efc5880968581e52b0 Mon Sep 17 00:00:00 2001 From: Jens Keiner Date: Tue, 8 Mar 2016 20:44:13 +0100 Subject: [PATCH 09/11] Ignore Eclipse .settings folder. --- .gitignore | 1 + 1 file changed, 1 insertion(+) diff --git a/.gitignore b/.gitignore index e6b73203..541e8cd6 100644 --- a/.gitignore +++ b/.gitignore @@ -34,6 +34,7 @@ # Other .cproject .project +.settings aclocal.m4 Makefile.in INSTALL From 3049064e100e212153de74cc2cb81896bb0feba5 Mon Sep 17 00:00:00 2001 From: Jens Keiner Date: Tue, 8 Mar 2016 21:26:16 +0100 Subject: [PATCH 10/11] Fix URL after repo name change. --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index d743c68c..37327cfa 100644 --- a/README.md +++ b/README.md @@ -1,4 +1,4 @@ -[![Build Status](https://travis-ci.org/NFFT/nfft_new.svg?branch=develop)](https://travis-ci.org/NFFT/nfft_new) +[![Build Status](https://travis-ci.org/NFFT/nfft.svg?branch=develop)](https://travis-ci.org/NFFT/nfft) NFFT - Nonequispaced FFT ========================= From 89d53180cde0ff141654cdbb64d3e92bafda139e Mon Sep 17 00:00:00 2001 From: Jens Keiner Date: Thu, 31 Mar 2016 23:06:11 +0200 Subject: [PATCH 11/11] #11: Temporarily add int version of next_power_of_2_exp() to avoid implicit type casts by the compiler that cause crashes with GCC 5. --- applications/iterS2/iterS2.c | 2 +- include/infft.h | 1 + kernel/fpt/fpt.c | 10 ++++----- kernel/nfsft/nfsft.c | 4 ++-- kernel/util/int.c | 39 ++++++++++++++++++++++++++++++++++++ 5 files changed, 48 insertions(+), 8 deletions(-) diff --git a/applications/iterS2/iterS2.c b/applications/iterS2/iterS2.c index 18c686fd..bc4c5934 100644 --- a/applications/iterS2/iterS2.c +++ b/applications/iterS2/iterS2.c @@ -280,7 +280,7 @@ int main (int argc, char **argv) /* */ if ((N+1)*(N+1) > M) { - X(next_power_of_2_exp)(N, &npt, &npt_exp); + X(next_power_of_2_exp_int)(N, &npt, &npt_exp); fprintf(stderr, "npt = %d, npt_exp = %d\n", npt, npt_exp); fprintf(stderr,"Optimal interpolation!\n"); ys = (double*) nfft_malloc((N+1)*sizeof(double)); diff --git a/include/infft.h b/include/infft.h index 26612458..e02c281f 100644 --- a/include/infft.h +++ b/include/infft.h @@ -1446,6 +1446,7 @@ R Y(prod_real)(R *vec, INT d); /* int.c: */ INT Y(log2i)(const INT m); void Y(next_power_of_2_exp)(const INT N, INT *N2, INT *t); +void Y(next_power_of_2_exp_int)(const int N, int *N2, int *t); /* error.c: */ /* not used */ R Y(error_l_infty_double)(const R *x, const R *y, const INT n); diff --git a/kernel/fpt/fpt.c b/kernel/fpt/fpt.c index 798b7323..4ffb47b9 100644 --- a/kernel/fpt/fpt.c +++ b/kernel/fpt/fpt.c @@ -1045,7 +1045,7 @@ void fpt_precompute(fpt_set set, const int m, double *alpha, double *beta, { /* Stabilize. */ degree_stab = degree*(2*l+1); - X(next_power_of_2_exp)((l+1)*(1<<(tau+1)),&N_stab,&t_stab); + X(next_power_of_2_exp_int)((l+1)*(1<<(tau+1)),&N_stab,&t_stab); /* Old arrays are to small. */ nfft_free(a11); @@ -1177,7 +1177,7 @@ void fpt_trafo_direct(fpt_set set, const int m, const double _Complex *x, double //fprintf(stderr, "Executing dpt.\n"); - X(next_power_of_2_exp)(k_end+1,&Nk,&tk); + X(next_power_of_2_exp_int)(k_end+1,&Nk,&tk); norm = 2.0/(Nk<<1); //fprintf(stderr, "Norm = %e.\n", norm); @@ -1275,7 +1275,7 @@ void fpt_trafo(fpt_set set, const int m, const double _Complex *x, double _Compl return; } - X(next_power_of_2_exp)(k_end,&Nk,&tk); + X(next_power_of_2_exp_int)(k_end,&Nk,&tk); k_start_tilde = K_START_TILDE(data->k_start,Nk); k_end_tilde = K_END_TILDE(k_end,Nk); @@ -1520,7 +1520,7 @@ void fpt_transposed_direct(fpt_set set, const int m, double _Complex *x, int tk; double norm; - X(next_power_of_2_exp)(k_end+1,&Nk,&tk); + X(next_power_of_2_exp_int)(k_end+1,&Nk,&tk); norm = 2.0/(Nk<<1); if (set->flags & FPT_NO_DIRECT_ALGORITHM) @@ -1607,7 +1607,7 @@ void fpt_transposed(fpt_set set, const int m, double _Complex *x, return; } - X(next_power_of_2_exp)(k_end,&Nk,&tk); + X(next_power_of_2_exp_int)(k_end,&Nk,&tk); k_start_tilde = K_START_TILDE(data->k_start,Nk); k_end_tilde = K_END_TILDE(k_end,Nk); diff --git a/kernel/nfsft/nfsft.c b/kernel/nfsft/nfsft.c index 79fe9d59..a66de0fc 100644 --- a/kernel/nfsft/nfsft.c +++ b/kernel/nfsft/nfsft.c @@ -280,7 +280,7 @@ void nfsft_init_guru(nfsft_plan *plan, int N, int M, unsigned int flags, /* Calculate the next greater power of two with respect to the bandwidth N * and the corresponding exponent. */ - //next_power_of_2_exp(plan->N,&plan->NPT,&plan->t); + //next_power_of_2_exp_int(plan->N,&plan->NPT,&plan->t); /* Save length of array of Fourier coefficients. Owing to the data layout the * length is (2N+2)(2N+2) */ @@ -382,7 +382,7 @@ void nfsft_precompute(int N, double kappa, unsigned int nfsft_flags, /* Compute and save N_max = 2^{\ceil{log_2 N}} as next greater * power of two with respect to N. */ - X(next_power_of_2_exp)(N,&wisdom.N_MAX,&wisdom.T_MAX); + X(next_power_of_2_exp_int)(N,&wisdom.N_MAX,&wisdom.T_MAX); /* Check, if precomputation for direct algorithms needs to be performed. */ if (wisdom.flags & NFSFT_NO_DIRECT_ALGORITHM) diff --git a/kernel/util/int.c b/kernel/util/int.c index 5210174e..c4cf51af 100644 --- a/kernel/util/int.c +++ b/kernel/util/int.c @@ -109,3 +109,42 @@ void Y(next_power_of_2_exp)(const INT N, INT *N2, INT *t) *t = logn+1; } } + +void Y(next_power_of_2_exp_int)(const int N, int *N2, int *t) +{ + int n,i,logn; + int N_is_not_power_of_2=0; + + if (N == 0) + { + *N2 = 1; + *t = 0; + } + else + { + n = N; + logn = 0; + while (n != 1) + { + if (n%2 == 1) + { + N_is_not_power_of_2=1; + } + n = n/2; + logn++; + } + + if (!N_is_not_power_of_2) + { + logn--; + } + + for (i = 0; i <= logn; i++) + { + n = n*2; + } + + *N2 = n; + *t = logn+1; + } +}