#4: Merge from develop. Resolve conflicts.

.gitignore

@@ -37,6 +37,7 @@
 # Other
 .cproject
 .project
+.settings
 aclocal.m4
 Makefile.in
 INSTALL

.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 $(if test "$CC" = "clang"; then echo ""; else echo "--enable-openmp"; fi) && make && make check
 
 ## Print config.log for debugging.
 after_failure: "cat config.log"

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
 =========================

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;
@@ -281,10 +280,9 @@ 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");
-      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_bsplines(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);

include/infft.h

@@ -187,13 +187,10 @@ 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) + \
-    (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);}
+  #define PHI(n,x,d) (Y(bsplines)(2*ths->m,((x)*n) + \
+    (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)
@@ -202,19 +199,16 @@ 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, 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)
@@ -1440,7 +1434,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(bsplines)(const INT, const R x);
 
 /* float.c: */
 typedef enum {NFFT_EPSILON = 0, NFFT_SAFE__MIN = 1, NFFT_BASE = 2,
@@ -1453,6 +1447,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);

kernel/fpt/fpt.c

@@ -1049,7 +1049,7 @@ void fpt_precompute(fpt_set set, const int m, R *alpha, R *beta,
         {
           /* Stabilize. */
           degree_stab = degree*(2*l+1);
-          Y(next_power_of_2_exp)((l+1)*(1<<(tau+1)),&N_stab,&t_stab);
+          Y(next_power_of_2_exp_int)((l+1)*(1<<(tau+1)),&N_stab,&t_stab);
 
           /* Old arrays are to small. */
           nfft_free(a11);
@@ -1181,7 +1181,7 @@ void fpt_trafo_direct(fpt_set set, const int m, const C *x, C *y,
 
     //fprintf(stderr, "Executing dpt.\n");  
 
-  Y(next_power_of_2_exp)(k_end+1,&Nk,&tk);
+  Y(next_power_of_2_exp_int)(k_end+1,&Nk,&tk);
   norm = 2.0/(Nk<<1);
 
     //fprintf(stderr, "Norm = %e.\n", norm);  
@@ -1279,7 +1279,7 @@ void fpt_trafo(fpt_set set, const int m, const C *x, C *y,
     return;
   }
 
-  Y(next_power_of_2_exp)(k_end,&Nk,&tk);
+  Y(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);
 
@@ -1524,7 +1524,7 @@ void fpt_transposed_direct(fpt_set set, const int m, C *x,
   int tk;
   R norm;
 
-  Y(next_power_of_2_exp)(k_end+1,&Nk,&tk);
+  Y(next_power_of_2_exp_int)(k_end+1,&Nk,&tk);
   norm = 2.0/(Nk<<1);
 
   if (set->flags & FPT_NO_DIRECT_ALGORITHM)
@@ -1611,7 +1611,7 @@ void fpt_transposed(fpt_set set, const int m, C *x,
     return;
   }
 
-  Y(next_power_of_2_exp)(k_end,&Nk,&tk);
+  Y(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);
 

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;
 
 /**

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);
       }
     }
   }

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)

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(bsplines)(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);
 

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;
+  }
+}

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@[email protected] -lcunit -lncurses
-
 clean-local:
 	rm -f CUnitAutomated-Results.xml CUnitAutomated_threads-Results.xml checkall.log checkall.trs checkall_threads.log checkall_threads.trs

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(bsplines)((INT)(n + 1), x);
     /*printf("x = " __FE__ ", err = " __FE__ "\n", x, ERR(y,yr));*/
     err = FMAX(err, ERR(y, yr));
   }