fastgauss.c

/*
 * Copyright (c) 2002, 2012 Jens Keiner, Stefan Kunis, Daniel Potts
 *
 * This program is free software; you can redistribute it and/or modify it under
 * the terms of the GNU General Public License as published by the Free Software
 * Foundation; either version 2 of the License, or (at your option) any later
 * version.
 *
 * This program is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
 * FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
 * details.
 *
 * You should have received a copy of the GNU General Public License along with
 * this program; if not, write to the Free Software Foundation, Inc., 51
 * Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

/* $Id: fastgauss.c 3896 2012-10-10 12:19:26Z tovo $ */

#include "config.h"

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#ifdef HAVE_COMPLEX_H
  #include <complex.h>
#endif

#include "nfft3.h"
#include "nfft3util.h"
#include "infft.h"

#define DGT_PRE_CEXP     (1U<< 0)

#define FGT_NDFT         (1U<< 1)

#define FGT_APPROX_B     (1U<< 2)

typedef struct
{
  int N;                                
  int M;                                
  double _Complex *alpha;                
  double _Complex *f;                    
  unsigned flags;                       
  double _Complex sigma;                 
  double *x;                            
  double *y;                            
  double _Complex *pre_cexp;             
  int n;                                
  double p;                             
  double _Complex *b;                    
  nfft_plan *nplan1;                    
  nfft_plan *nplan2;                    
} fgt_plan;

void dgt_trafo(fgt_plan *ths)
{
  int j,k,l;

  for(j=0; j<ths->M; j++)
    ths->f[j] = 0;

  if(ths->flags & DGT_PRE_CEXP)
    for(j=0,l=0; j<ths->M; j++)
      for(k=0; k<ths->N; k++,l++)
        ths->f[j] += ths->alpha[k]*ths->pre_cexp[l];
  else
    for(j=0; j<ths->M; j++)
      for(k=0; k<ths->N; k++)
        ths->f[j] += ths->alpha[k]*cexp(-ths->sigma*(ths->y[j]-ths->x[k])*
          (ths->y[j]-ths->x[k]));
}

void fgt_trafo(fgt_plan *ths)
{
  int l;

  if(ths->flags & FGT_NDFT)
    {
      nfft_adjoint_direct(ths->nplan1);

      for(l=0; l<ths->n; l++)
        ths->nplan1->f_hat[l] *= ths->b[l];

      nfft_trafo_direct(ths->nplan2);
    }
  else
    {
      nfft_adjoint(ths->nplan1);

      for(l=0; l<ths->n; l++)
        ths->nplan1->f_hat[l] *= ths->b[l];

      nfft_trafo(ths->nplan2);
    }
}

void fgt_init_guru(fgt_plan *ths, int N, int M, double _Complex sigma, int n,
       double p, int m, unsigned flags)
{
  int j,n_fftw;
  fftw_plan fplan;

  ths->M = M;
  ths->N = N;
  ths->sigma = sigma;
  ths->flags = flags;

  ths->x = (double*)nfft_malloc(ths->N*sizeof(double));
  ths->y = (double*)nfft_malloc(ths->M*sizeof(double));
  ths->alpha = (double _Complex*)nfft_malloc(ths->N*sizeof(double _Complex));
  ths->f = (double _Complex*)nfft_malloc(ths->M*sizeof(double _Complex));

  ths->n = n;
  ths->p = p;

  ths->b = (double _Complex*)nfft_malloc(ths->n*sizeof(double _Complex));

  ths->nplan1 = (nfft_plan*) nfft_malloc(sizeof(nfft_plan));
  ths->nplan2 = (nfft_plan*) nfft_malloc(sizeof(nfft_plan));

  n_fftw=X(next_power_of_2)(2*ths->n);

  nfft_init_guru(ths->nplan1, 1, &(ths->n), ths->N, &n_fftw, m, PRE_PHI_HUT|
                 PRE_PSI| MALLOC_X| MALLOC_F_HAT| FFTW_INIT, FFTW_MEASURE);
  nfft_init_guru(ths->nplan2, 1, &(ths->n), ths->M, &n_fftw, m, PRE_PHI_HUT|
                 PRE_PSI| MALLOC_X| FFTW_INIT, FFTW_MEASURE);

  ths->nplan1->f = ths->alpha;
  ths->nplan2->f_hat = ths->nplan1->f_hat;
  ths->nplan2->f = ths->f;

  if(ths->flags & FGT_APPROX_B)
    {
      fplan = fftw_plan_dft_1d(ths->n, ths->b, ths->b, FFTW_FORWARD,
                               FFTW_MEASURE);

      for(j=0; j<ths->n; j++)
  ths->b[j] = cexp(-ths->p*ths->p*ths->sigma*(j-ths->n/2)*(j-ths->n/2)/
                          ((double)ths->n*ths->n)) / ths->n;

      nfft_fftshift_complex(ths->b, 1, &ths->n);
      fftw_execute(fplan);
      nfft_fftshift_complex(ths->b, 1, &ths->n);

      fftw_destroy_plan(fplan);
    }
  else
    {
      for(j=0; j<ths->n; j++)
  ths->b[j] = 1.0/ths->p * csqrt(PI/ths->sigma)*
    cexp(-PI*PI*(j-ths->n/2)*(j-ths->n/2)/
         (ths->p*ths->p*ths->sigma));
    }
}

void fgt_init(fgt_plan *ths, int N, int M, double _Complex sigma, double eps)
{
  double p;
  int n;

  p=0.5+sqrt(-log(eps)/creal(sigma));
  if(p<1)
    p=1;

  n=2*((int)ceil(p*cabs(sigma)/PI * sqrt(-log(eps)/creal(sigma))));

  if(N*M<=((int)(1U<<20)))
    fgt_init_guru(ths, N, M, sigma, n, p, 7, DGT_PRE_CEXP);
  else
    fgt_init_guru(ths, N, M, sigma, n, p, 7, 0);
}

void fgt_init_node_dependent(fgt_plan *ths)
{
  int j,k,l;

  if(ths->flags & DGT_PRE_CEXP)
   {
     ths->pre_cexp=(double _Complex*)nfft_malloc(ths->M*ths->N*
            sizeof(double _Complex));

     for(j=0,l=0; j<ths->M; j++)
       for(k=0; k<ths->N; k++,l++)
         ths->pre_cexp[l]=cexp(-ths->sigma*(ths->y[j]-ths->x[k])*
                                (ths->y[j]-ths->x[k]));
   }

  for(j=0; j<ths->nplan1->M_total; j++)
    ths->nplan1->x[j] = ths->x[j]/ths->p;
  for(j=0; j<ths->nplan2->M_total; j++)
    ths->nplan2->x[j] = ths->y[j]/ths->p;

  if(ths->nplan1->nfft_flags & PRE_PSI)
    nfft_precompute_psi(ths->nplan1);
  if(ths->nplan2->nfft_flags & PRE_PSI)
    nfft_precompute_psi(ths->nplan2);
}

void fgt_finalize(fgt_plan *ths)
{
  nfft_finalize(ths->nplan2);
  nfft_finalize(ths->nplan1);

  nfft_free(ths->nplan2);
  nfft_free(ths->nplan1);

  nfft_free(ths->b);

  nfft_free(ths->f);
  nfft_free(ths->y);

  nfft_free(ths->alpha);
  nfft_free(ths->x);
}

void fgt_test_init_rand(fgt_plan *ths)
{
  int j,k;

  for(k=0; k<ths->N; k++)
    ths->x[k] = (double)rand()/(2.0*RAND_MAX)-1.0/4.0;

  for(j=0; j<ths->M; j++)
    ths->y[j] = (double)rand()/(2.0*RAND_MAX)-1.0/4.0;

  for(k=0; k<ths->N; k++)
    ths->alpha[k] =   (double)rand()/(RAND_MAX)-1.0/2.0
            + _Complex_I*(double)rand()/(RAND_MAX)-I/2.0;
}

double fgt_test_measure_time(fgt_plan *ths, unsigned dgt)
{
  int r;
  ticks t0, t1;
  double t_out;
  double tau=0.01;

  t_out=0;
  r=0;
  while(t_out<tau)
    {
      r++;
      t0 = getticks();
      if (dgt)
        dgt_trafo(ths);
      else
        fgt_trafo(ths);
      t1 = getticks();
      t_out += nfft_elapsed_seconds(t1,t0);
    }
  t_out/=r;

  return t_out;
}

void fgt_test_simple(int N, int M, double _Complex sigma, double eps)
{
  fgt_plan my_plan;
  double _Complex *swap_dgt;

  fgt_init(&my_plan, N, M, sigma, eps);
  swap_dgt = (double _Complex*)nfft_malloc(my_plan.M*sizeof(double _Complex));

  fgt_test_init_rand(&my_plan);
  fgt_init_node_dependent(&my_plan);

  NFFT_SWAP_complex(swap_dgt,my_plan.f);
  dgt_trafo(&my_plan);
  nfft_vpr_complex(my_plan.f,my_plan.M,"discrete gauss transform");
  NFFT_SWAP_complex(swap_dgt,my_plan.f);

  fgt_trafo(&my_plan);
  nfft_vpr_complex(my_plan.f,my_plan.M,"fast gauss transform");

  printf("\n relative error: %1.3e\n", X(error_l_infty_1_complex)(swap_dgt,
         my_plan.f, my_plan.M, my_plan.alpha, my_plan.N));

  nfft_free(swap_dgt);
  fgt_finalize(&my_plan);
}

void fgt_test_andersson(void)
{
  fgt_plan my_plan;
  double _Complex *swap_dgt;
  int N;

  double _Complex sigma=4*(138+ _Complex_I*100);
  int n=128;
  int N_dgt_pre_exp=(int)(1U<<11);
  int N_dgt=(int)(1U<<19);

  printf("n=%d, sigma=%1.3e+i%1.3e\n",n,creal(sigma),cimag(sigma));

  for(N=((int)(1U<<6)); N<((int)(1U<<22)); N=N<<1)
    {
      printf("$%d$\t & ",N);

      if(N<N_dgt_pre_exp)
        fgt_init_guru(&my_plan, N, N, sigma, n, 1, 7, DGT_PRE_CEXP);
      else
        fgt_init_guru(&my_plan, N, N, sigma, n, 1, 7, 0);

      swap_dgt = (double _Complex*)nfft_malloc(my_plan.M*
                sizeof(double _Complex));

      fgt_test_init_rand(&my_plan);

      fgt_init_node_dependent(&my_plan);

      if(N<N_dgt)
  {
          NFFT_SWAP_complex(swap_dgt,my_plan.f);
          if(N<N_dgt_pre_exp)
            my_plan.flags^=DGT_PRE_CEXP;

    printf("$%1.1e$\t & ",fgt_test_measure_time(&my_plan, 1));
          if(N<N_dgt_pre_exp)
            my_plan.flags^=DGT_PRE_CEXP;
          NFFT_SWAP_complex(swap_dgt,my_plan.f);
  }
      else
  printf("\t\t & ");

      if(N<N_dgt_pre_exp)
  printf("$%1.1e$\t & ",fgt_test_measure_time(&my_plan, 1));
      else
  printf("\t\t & ");

      my_plan.flags^=FGT_NDFT;
      printf("$%1.1e$\t & ",fgt_test_measure_time(&my_plan, 0));
      my_plan.flags^=FGT_NDFT;

      printf("$%1.1e$\t & ",fgt_test_measure_time(&my_plan, 0));

      printf("$%1.1e$\t \\\\ \n",
       X(error_l_infty_1_complex)(swap_dgt, my_plan.f, my_plan.M,
            my_plan.alpha, my_plan.N));
      fflush(stdout);

      nfft_free(swap_dgt);
      fgt_finalize(&my_plan);
      fftw_cleanup();
    }
}

void fgt_test_error(void)
{
  fgt_plan my_plan;
  double _Complex *swap_dgt;
  int n,mi;

  double _Complex sigma=4*(138+ _Complex_I*100);
  int N=1000;
  int M=1000;
  int m[2]={7,3};

  printf("N=%d;\tM=%d;\nsigma=%1.3e+i*%1.3e;\n",N,M,creal(sigma),cimag(sigma));
  printf("error=[\n");

  swap_dgt = (double _Complex*)nfft_malloc(M*sizeof(double _Complex));

  for(n=8; n<=128; n+=4)
    {
      printf("%d\t",n);
      for(mi=0;mi<2;mi++)
        {
          fgt_init_guru(&my_plan, N, M, sigma, n, 1, m[mi], 0);
          fgt_test_init_rand(&my_plan);
          fgt_init_node_dependent(&my_plan);

          NFFT_SWAP_complex(swap_dgt,my_plan.f);
          dgt_trafo(&my_plan);
          NFFT_SWAP_complex(swap_dgt,my_plan.f);

          fgt_trafo(&my_plan);

          printf("%1.3e\t", X(error_l_infty_1_complex)(swap_dgt, my_plan.f,
                 my_plan.M, my_plan.alpha, my_plan.N));
          fflush(stdout);

          fgt_finalize(&my_plan);
          fftw_cleanup();
        }
      printf("\n");
    }
  printf("];\n");

  nfft_free(swap_dgt);
}

void fgt_test_error_p(void)
{
  fgt_plan my_plan;
  double _Complex *swap_dgt;
  int n,pi;

  double _Complex sigma=20+40*_Complex_I;
  int N=1000;
  int M=1000;
  double p[3]={1,1.5,2};

  printf("N=%d;\tM=%d;\nsigma=%1.3e+i*%1.3e;\n",N,M,creal(sigma),cimag(sigma));
  printf("error=[\n");

  swap_dgt = (double _Complex*)nfft_malloc(M*sizeof(double _Complex));

  for(n=8; n<=128; n+=4)
    {
      printf("%d\t",n);
      for(pi=0;pi<3;pi++)
        {
          fgt_init_guru(&my_plan, N, M, sigma, n, p[pi], 7, 0);
          fgt_test_init_rand(&my_plan);
          fgt_init_node_dependent(&my_plan);

          NFFT_SWAP_complex(swap_dgt,my_plan.f);
          dgt_trafo(&my_plan);
          NFFT_SWAP_complex(swap_dgt,my_plan.f);

          fgt_trafo(&my_plan);

          printf("%1.3e\t", X(error_l_infty_1_complex)(swap_dgt, my_plan.f,
                 my_plan.M, my_plan.alpha, my_plan.N));
          fflush(stdout);

          fgt_finalize(&my_plan);
          fftw_cleanup();
        }
      printf("\n");
    }
  printf("];\n");
}

int main(int argc,char **argv)
{
  if(argc!=2)
    {
      fprintf(stderr,"fastgauss type\n");
      fprintf(stderr," type\n");
      fprintf(stderr,"  0 - Simple test.\n");
      fprintf(stderr,"  1 - Compares accuracy and execution time.\n");
      fprintf(stderr,"      Pipe to output_andersson.tex\n");
      fprintf(stderr,"  2 - Compares accuracy.\n");
      fprintf(stderr,"      Pipe to output_error.m\n");
      fprintf(stderr,"  3 - Compares accuracy.\n");
      fprintf(stderr,"      Pipe to output_error_p.m\n");
      return -1;
    }

  if(atoi(argv[1])==0)
    fgt_test_simple(10, 10, 5+3*_Complex_I, 0.001);

  if(atoi(argv[1])==1)
    fgt_test_andersson();

  if(atoi(argv[1])==2)
    fgt_test_error();

  if(atoi(argv[1])==3)
    fgt_test_error_p();

  return 1;
}
/* \} */