Merge pull request #142 from NFFT/feature/nfft_trafo_direct_d_1

Improve efficiency of NFFT direct transformation in one dimension.

Author: michaelquellmalz

kernel/nfft/nfft.c

@@ -146,8 +146,6 @@ void X(trafo_direct)(const X(plan) *ths)
 {
   C *f_hat = (C*)ths->f_hat, *f = (C*)ths->f;
 
-  memset(f, 0, (size_t)(ths->M_total) * sizeof(C));
-
   if (ths->d == 1)
   {
     /* specialize for univariate case, rationale: faster */
@@ -157,16 +155,22 @@ void X(trafo_direct)(const X(plan) *ths)
 #endif
     for (j = 0; j < ths->M_total; j++)
     {
+      C v = K(0.0);
       INT k_L;
       for (k_L = 0; k_L < ths->N_total; k_L++)
       {
         R omega = K2PI * ((R)(k_L - ths->N_total/2)) * ths->x[j];
-        f[j] += f_hat[k_L] * BASE(-II * omega);
+        v += f_hat[k_L] * (COS(omega) - II * SIN(omega));
       }
+
+      f[j] = v;
     }
   }
   else
   {
+
+    memset(f, 0, (size_t)(ths->M_total) * sizeof(C));
+
     /* multivariate case */
     INT j;
 #ifdef _OPENMP

tests/nfft.c

@@ -783,13 +783,13 @@ static R compare_adjoint(check_delegate_t *ego, X(plan) *p, const int NN, const
 static check_delegate_t check_trafo = {prepare_trafo, compare_trafo};
 static check_delegate_t check_adjoint = {prepare_adjoint, compare_adjoint};
 
-static trafo_delegate_t trafo_direct = {"trafo_direct", X(trafo_direct), 0, trafo_direct_cost, err_trafo_direct};
+static trafo_delegate_t trafo_direct = {"trafo_direct", (trafo_t)X(trafo_direct), 0, trafo_direct_cost, err_trafo_direct};
 static trafo_delegate_t trafo = {"trafo", X(trafo), X(check), 0, err_trafo};
 static trafo_delegate_t trafo_1d = {"trafo_1d", X(trafo_1d), X(check), 0, err_trafo};
 static trafo_delegate_t trafo_2d = {"trafo_2d", X(trafo_2d), X(check), 0, err_trafo};
 static trafo_delegate_t trafo_3d = {"trafo_3d", X(trafo_3d), X(check), 0, err_trafo};
 
-static trafo_delegate_t adjoint_direct = {"adjoint_direct", X(adjoint_direct), 0, trafo_direct_cost, err_trafo_direct};
+static trafo_delegate_t adjoint_direct = {"adjoint_direct", (trafo_t)X(adjoint_direct), 0, trafo_direct_cost, err_trafo_direct};
 static trafo_delegate_t adjoint = {"adjoint", X(adjoint), X(check), 0, err_trafo};
 static trafo_delegate_t adjoint_1d = {"adjoint_1d", X(adjoint_1d), X(check), 0, err_trafo};
 static trafo_delegate_t adjoint_2d = {"adjoint_2d", X(adjoint_2d), X(check), 0, err_trafo};