nfft_adjoint() crash with small N_total

[krivenko]

Dear authors,

I have found that the following simple program crashes (segfault) when compiled and linked against NFFT 3.3.2.

#include "nfft3.h"

int main(void) {
  int N[] = {8};        /* Works for >= 9 */
  int M = 10;
  nfft_plan p;
  nfft_init(&p, 1, N, M);

  for(int i = 0; i < M; ++i) {
   p.x[i] = -0.5 + (double)i / M;
   p.f[i][0] = 1.0;     /* Memory corruption? */
   p.f[i][1] = 0.0;
  }

  if(p.flags & PRE_ONE_PSI) nfft_precompute_one_psi(&p);
  nfft_check(&p);
  nfft_adjoint(&p);     /* Segfault if M <= 6 */
  nfft_vpr_complex(p.f_hat, p.N_total, "adjoint nfft, vector f_hat");
  nfft_finalize(&p);    /* Segfault if M > 6 */
  return 0;
}

Some remarks are to be made about the program.

• It successfully runs to completion if N is 9 or larger;
• For M <= 6 the segmentation fault occurs in nfft_adjoint();
• For larger M it occurs in nfft_finalize();
• There is no crash if p.f is untouched.

Is this a bug or am I just misusing NFFT?

[tvolkmer]

If you use the "simple" init for NFFT and the Kaiser-Bessel window function (configure default), then the parameter N should be at least 12. Also for N=9, there may be corrupted memory or segfault depending on the other parameters and the x values.

Moreover, if you use nfft_check please also test the corresponding result:

#include "nfft3.h"

int main(void) {
  int N[] = {8};
  int M = 10;
  nfft_plan p;

  nfft_init(&p, 1, N, M);

  for(int i = 0; i < M; ++i) {
   p.x[i] = -0.5 + (double)i / M;
   p.f[i][0] = 1.0;     /* Memory corruption? */
   p.f[i][1] = 0.0;
  }

  if(p.flags & PRE_ONE_PSI) nfft_precompute_one_psi(&p);

  const char *check_error_msg = nfft_check(&p);
  if (check_error_msg != 0) {
    printf("Invalid nfft parameter: %s\n", check_error_msg);
    return 1;
  }

  nfft_adjoint(&p);     /* Segfault if M <= 6 */
  nfft_vpr_complex(p.f_hat, p.N_total, "adjoint nfft, vector f_hat");
  nfft_finalize(&p);    /* Segfault if M > 6 */
  return 0;
}

This example will produce the output
Invalid nfft parameter: Polynomial degree N is smaller than cut-off m
and return 1.

Whenever the nfft_check function returns something != 0 and nfft_trafo / nfft_adjoint are called afterwards, there may be some segfault or corrupted memory.

I think we should improve our documentation...

[krivenko]

Thank you for the prompt answer, Toni!

Moreover, if you use nfft_check please also test the corresponding result:

Sorry, my bad. Indeed, calling nfft_check() without actually checking its return value is pointless.

If you use the "simple" init for NFFT and the Kaiser-Bessel window function (configure default), then the parameter N should be at least 12. Also for N=9, there may be corrupted memory or segfault depending on the other parameters and the x values.

O.k., what options do I have if I need only a few frequencies?
I guess for N < 12 I could set N to 12 and discard the unneeded frequencies, but this is a bit of overhead...

[tvolkmer]

Dear Igor,

thank you for your feedback. This reminds us that we should write some things more precisely in the documentation.

Please note that the parameter N should always be an even positive integer.

For small even N, please use nfft_trafo_direct / nfft_adjoint_direct. These methods perform the "slow"/direct computation, which might also be faster.

Only for "not too small" N, the fast computation is really faster than the direct one. The exact "cross-over" N depends on many factors like dimension d, number of nonequispaced nodes M, window cut-off parameter m (=12 by default, influences the precision), your computer, ...

@jenskeiner We should think about calling nfft_trafo_direct internally in nfft_trafo and nfft_adjoint_direct internally in nfft_adjoint for very small values of N.

[krivenko]

Thank you very much!
I have updated my code to use nfft_adjoint_direct() for the smaller N_total, and everything works fine now.

This reminds us that we should write some things more precisely in the documentation.

Indeed, nfft_trafo_direct / nfft_adjoint_direct are not mentioned in docs on your website...

Feel free to close this issue, unless you are going to keep it as a reminder.

[jenskeiner]

Yes, we should internally call the direct trafo when the asymptotically faster one doesn't make sense. This is a legacy from before the time I got involved with the project. Best would be an FFTW-style solver concept but that would require a huge re-factoring and add a lot of additional infrastructure. We can start with some smaller improvements...

[jenskeiner]

Closing this one. Will create a new ticket to call the direct trafo automatically for smaller sizes.

[jenskeiner]

Can this be merged to develop? If so, please create a pull request from the feature branch.