nfft.c

用于实现非均匀采样FFT（NDFT）

					   (2*this->m+2)*sizeof(double));
  
  if(this->nfft_flags & FFTW_INIT)
    {
      this->g1=(fftw_complex*)fftw_malloc(this->n_L*sizeof(fftw_complex));
      
      if(this->nfft_flags & FFT_OUT_OF_PLACE)
	this->g2 = (fftw_complex*) fftw_malloc(this->n_L*
					       sizeof(fftw_complex));
      else
	this->g2 = this->g1;
      
      this->my_fftw_plan1 = 
	fftw_plan_dft(this->d, this->n, this->g1, this->g2,
		      FFTW_FORWARD, this->fftw_flags);
      this->my_fftw_plan2 = 
	fftw_plan_dft(this->d, this->n, this->g2, this->g1,
		      FFTW_BACKWARD, this->fftw_flags);
    }
}

void nfft_init(nfft_plan *this, int d, int *N, int M)
{
  int t;                                /**< index over all dimensions        */

  this->d = d;

  this->N=(int*) fftw_malloc(d*sizeof(int));
  for(t = 0;t < d; t++)
    this->N[t] = N[t];

  this->M = M;

  this->n = (int*) fftw_malloc(d*sizeof(int));
  for(t = 0;t < d; t++)
    this->n[t] = 2*next_power_of_2(this->N[t]);

  WINDOW_HELP_ESTIMATE_m;

  this->nfft_flags= PRE_PHI_HUT| PRE_PSI| MALLOC_X| MALLOC_F_HAT| MALLOC_F|
                    FFTW_INIT| FFT_OUT_OF_PLACE;
  this->fftw_flags= FFTW_ESTIMATE| FFTW_DESTROY_INPUT;

  nfft_init_help(this);    
}

void nfft_init_specific(nfft_plan *this, int d, int *N, int M, int *n,
			int m, unsigned nfft_flags, unsigned fftw_flags)
{
  int t;                                /**< index over all dimensions        */

  this->d =d;
  this->N= (int*) fftw_malloc(this->d*sizeof(int));
  for(t=0; t<d; t++)
    this->N[t]= N[t];
  this->M= M;
  this->n= (int*) fftw_malloc(this->d*sizeof(int));
  for(t=0; t<d; t++)
    this->n[t]= n[t];
  this->m= m;
  this->nfft_flags= nfft_flags;
  this->fftw_flags= fftw_flags;

  nfft_init_help(this);  
}

void nfft_init_1d(nfft_plan *this, int N1, int M)
{
  int N[1];

  N[0]=N1;
  nfft_init(this,1,N,M);
}

void nfft_init_2d(nfft_plan *this, int N1, int N2, int M)
{
  int N[2];

  N[0]=N1;
  N[1]=N2;
  nfft_init(this,2,N,M);
}

void nfft_init_3d(nfft_plan *this, int N1, int N2, int N3, int M)
{
  int N[3];

  N[0]=N1;
  N[1]=N2;
  N[2]=N3;
  nfft_init(this,3,N,M);
}

void nfft_finalize(nfft_plan *this)
{
  int t; /* index over dimensions */

  if(this->nfft_flags & FFTW_INIT)
    {
      fftw_destroy_plan(this->my_fftw_plan2);
      fftw_destroy_plan(this->my_fftw_plan1);
      
      if(this->nfft_flags & FFT_OUT_OF_PLACE)
	fftw_free(this->g2);
      
      fftw_free(this->g1);
    }

  /* NO FFTW_FREE HERE */
  if(this->nfft_flags & PRE_PSI)
    {
      if(this->nfft_flags & PRE_FULL_PSI)
	{
	  free(this->psi_index_g);
	  free(this->psi_index_f);
	}

      free(this->psi);
    }
      
  if(this->nfft_flags & PRE_PHI_HUT)
    {
      for(t=0; t<this->d; t++)
        fftw_free(this->c_phi_inv[t]);
      fftw_free(this->c_phi_inv);
    }

  if(this->nfft_flags & MALLOC_F)
    fftw_free(this->f);

  if(this->nfft_flags & MALLOC_F_HAT)
    fftw_free(this->f_hat);

  if(this->nfft_flags & MALLOC_X)
  fftw_free(this->x);
 
  WINDOW_HELP_FINALIZE;
 
  fftw_free(this->sigma);
  fftw_free(this->n);
  fftw_free(this->N);
}


/** inverse nfft 
 * -----------------------------------------------------------------------------
 * -----------------------------------------------------------------------------
 */
void infft_init_specific(infft_plan *this_iplan, nfft_plan *direct_plan,
			 int infft_flags)
{
  this_iplan->direct_plan = direct_plan;
  this_iplan->infft_flags = infft_flags;
  
  this_iplan->given_f = (fftw_complex*)
    fftw_malloc(this_iplan->direct_plan->M*sizeof(fftw_complex));

  this_iplan->r_iter = (fftw_complex*)
    fftw_malloc(this_iplan->direct_plan->M*sizeof(fftw_complex));

  this_iplan->f_hat_iter = (fftw_complex*)
    fftw_malloc(this_iplan->direct_plan->N_L*sizeof(fftw_complex));

  this_iplan->p_hat_iter = (fftw_complex*)
    fftw_malloc(this_iplan->direct_plan->N_L*sizeof(fftw_complex));

  if(this_iplan->infft_flags & LANDWEBER)
    {
      this_iplan->z_hat_iter = this_iplan->p_hat_iter;
    }

  if(this_iplan->infft_flags & STEEPEST_DESCENT)
    {
      this_iplan->z_hat_iter = this_iplan->p_hat_iter;

      this_iplan->v_iter = (fftw_complex*)
	fftw_malloc(this_iplan->direct_plan->M*sizeof(fftw_complex));
    }

  if(this_iplan->infft_flags & CGNR_E)
    {
      this_iplan->z_hat_iter = (fftw_complex*)
	fftw_malloc(this_iplan->direct_plan->N_L*sizeof(fftw_complex));

      this_iplan->v_iter = (fftw_complex*)
	fftw_malloc(this_iplan->direct_plan->M*sizeof(fftw_complex));
    }

  if(this_iplan->infft_flags & CGNE_R)
    {
      this_iplan->z_hat_iter = this_iplan->p_hat_iter;
    }

  if(this_iplan->infft_flags & ITERATE_2nd)
    this_iplan->f_hat_iter_2nd = (fftw_complex*) 
      fftw_malloc(this_iplan->direct_plan->N_L*sizeof(fftw_complex));

  if(this_iplan->infft_flags & PRECOMPUTE_WEIGHT)
    this_iplan->w = 
      (double*) fftw_malloc(this_iplan->direct_plan->M*sizeof(double));
  
  if(this_iplan->infft_flags & PRECOMPUTE_DAMP)
    this_iplan->w_hat = 
      (double*) fftw_malloc(this_iplan->direct_plan->N_L*sizeof(double));
}

void infft_init(infft_plan *this_iplan, nfft_plan *direct_plan)
{
  infft_init_specific(this_iplan, direct_plan, CGNR_E);
}

void infft_before_loop_help(infft_plan *this_iplan)
{
  /** step 2
   *  overwrites this_iplan->direct_plan->f_hat 
   *  overwrites this_iplan->r_iter
   */
  copyc(this_iplan->direct_plan->f_hat, this_iplan->f_hat_iter,
	this_iplan->direct_plan->N_L);
  
  SWAPC(this_iplan->r_iter, this_iplan->direct_plan->f);
  nfft_trafo(this_iplan->direct_plan);
  SWAPC(this_iplan->r_iter, this_iplan->direct_plan->f);

  updatec_axpy(this_iplan->r_iter, -1.0, this_iplan->given_f,
	       this_iplan->direct_plan->M);

  if((!(this_iplan->infft_flags & LANDWEBER)) ||
     (this_iplan->infft_flags & NORMS_FOR_LANDWEBER))
    {
      if(this_iplan->infft_flags & PRECOMPUTE_WEIGHT)
	this_iplan->dot_r_iter =
	  dotproductc_w(this_iplan->r_iter, this_iplan->w,
			this_iplan->direct_plan->M);
      else
	this_iplan->dot_r_iter =
	  dotproductc(this_iplan->r_iter, this_iplan->direct_plan->M);
    }
  
  /** step 3
   *  overwrites this_iplan->direct_plan->f
   *  overwrites this_iplan->z_hat_iter resp. this_iplan->z_hat_iter
   */ 

  if(this_iplan->infft_flags & PRECOMPUTE_WEIGHT)
    copyc_w(this_iplan->direct_plan->f, this_iplan->w, this_iplan->r_iter,
	    this_iplan->direct_plan->M);
  else
    copyc(this_iplan->direct_plan->f, this_iplan->r_iter,
	  this_iplan->direct_plan->M);
  
  SWAPC(this_iplan->z_hat_iter, this_iplan->direct_plan->f_hat);
  nfft_adjoint(this_iplan->direct_plan);
  SWAPC(this_iplan->z_hat_iter, this_iplan->direct_plan->f_hat);
  
  if((!(this_iplan->infft_flags & LANDWEBER)) ||
     (this_iplan->infft_flags & NORMS_FOR_LANDWEBER))
    {
      if(this_iplan->infft_flags & PRECOMPUTE_DAMP)
	this_iplan->dot_z_hat_iter = 
	  dotproductc_w(this_iplan->z_hat_iter, this_iplan->w_hat,
			this_iplan->direct_plan->N_L);
      else
	this_iplan->dot_z_hat_iter =
	  dotproductc(this_iplan->z_hat_iter, this_iplan->direct_plan->N_L);
    }

  if(this_iplan->infft_flags & CGNE_R)
    this_iplan->dot_p_hat_iter = this_iplan->dot_z_hat_iter;

} /* void infft_before_loop_help */

void infft_before_loop(infft_plan *this_iplan)
{
  infft_before_loop_help(this_iplan);

  if(this_iplan->infft_flags & CGNR_E)
    {
      /** step 4-6
       *  overwrites this_iplan->f_hat_iter_2nd
       */
      if(this_iplan->infft_flags & ITERATE_2nd)
	copyc(this_iplan->f_hat_iter_2nd, this_iplan->f_hat_iter,
	      this_iplan->direct_plan->N_L);
     
      /** step 7
       *  overwrites this_iplan->p_hat_iter
       */
      copyc(this_iplan->p_hat_iter, this_iplan->z_hat_iter,
	    this_iplan->direct_plan->N_L);
    }

  if(this_iplan->infft_flags & CGNE_R)
    {
      /** step 4-7
       *  overwrites this_iplan->f_hat_iter_2nd
       */
      if(this_iplan->infft_flags & ITERATE_2nd)
	{
	  this_iplan->gamma_iter=1.0;
	  copyc(this_iplan->f_hat_iter_2nd, this_iplan->f_hat_iter,
		this_iplan->direct_plan->N_L);
	}
    }
} /* void infft_before_loop */

void infft_loop_one_step_landweber(infft_plan *this_iplan)
{
  /** step 5
   *  updates this_iplan->f_hat_iter
   */
  if(this_iplan->infft_flags & PRECOMPUTE_DAMP)
    updatec_xpawy(this_iplan->f_hat_iter, this_iplan->alpha_iter,
		  this_iplan->w_hat, this_iplan->z_hat_iter,
		  this_iplan->direct_plan->N_L);
  else
    updatec_xpay(this_iplan->f_hat_iter, this_iplan->alpha_iter,
		this_iplan->z_hat_iter, this_iplan->direct_plan->N_L);
  
  /** step 6
   *  original residual, not the updated residual,
   *  overwrites this_iplan->r_iter
   *  overwrites this_iplan->direct_plan->f_hat 
   */
  copyc(this_iplan->direct_plan->f_hat, this_iplan->f_hat_iter,
	this_iplan->direct_plan->N_L);

  SWAPC(this_iplan->r_iter,this_iplan->direct_plan->f);
  nfft_trafo(this_iplan->direct_plan);
  SWAPC(this_iplan->r_iter,this_iplan->direct_plan->f);
  
  updatec_axpy(this_iplan->r_iter, -1.0, this_iplan->given_f, 
	       this_iplan->direct_plan->M);

  if(this_iplan->infft_flags & NORMS_FOR_LANDWEBER)
    {
      if(this_iplan->infft_flags & PRECOMPUTE_WEIGHT)
	this_iplan->dot_r_iter = dotproductc_w(this_iplan->r_iter,this_iplan->w,
					       this_iplan->direct_plan->M);
      else
	this_iplan->dot_r_iter =
	  dotproductc(this_iplan->r_iter, this_iplan->direct_plan->M);
    }