/** \file myfft.c
 * FastFourierTransformations for parallel computing.
 * Implements FFT on multiple machines using fftw package. The most important
 * routines are the hither and back transformations: fft_3d_real_to_complex() and
 * fft_3d_complex_to_real(). For fftw plans are needed, which are created by
 * fft_3d_create_plan() and destroyed by fft_3d_destroy_plan(). Plan weights
 * can be copied CopyElementOnepw() and there are sort critera GetKeyOnepw(),
 * GetKeyOnepw2PZI().\n
 * Wave functions can - for the purpose of testing - be filled with random values
 * by InitDataTestR().
 * 
  Project: ParallelCarParrinello
 \author Jan Hamaekers
 \date 2000

  File: myfft.c
  $Id: myfft.c,v 1.23 2007/02/05 15:28:39 foo Exp $
*/

#include<stdlib.h>
#include<stdio.h>
#include<unistd.h>
#include<math.h>
#include<string.h>

// use double precision fft when we have it
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#ifdef HAVE_DFFTW_H
#include "dfftw.h"
#else
#include "fftw.h"
#endif

#ifdef HAVE_DRFFTW_H
#include "drfftw.h"
#else
#include "rfftw.h"
#endif

#include"data.h"
#include"helpers.h"
#include"output.h"
#include"errors.h"
#include"mergesort2.h"
#include"mymath.h"
#include"myfft.h"


/** Returns negative value of the weight of a plan_weight as sort critera.
 * \param *a plan_weight array
 * \param i index of element
 * \param *Args unused (only for contingency)
 */
static double GetKeyOnepw(void *a, int i, void *Args) {
  return(-((struct plan_weight *)a)[i].weight);
}

/** Returns "index + twice the process number + \a Args[2]" of a plan_weight as sort critera.
 * \a Args[2] is only added if: Index % \a Args[1] != Args[1] - 1 
 * \param *a plan_weight array
 * \param i index of element
 * \param *Args 
 */
static double GetKeyOnepw2PZI(void *a, int i, void *Args) {
  return(((struct plan_weight *)a)[i].Index+((int *)Args)[2]*( ( ((struct plan_weight *)a)[i].Index%((int *)Args)[1]!=((int *)Args)[1]-1)+2*(((struct plan_weight *)a)[i].PE)));
}

/** Copies one element from a plan weight array \a *b to another \a *a.
 * Copies from \a i-th element of \a *a index, weight and PE to \a *b
 * \param *a source plan weight array
 * \param i index of source element
 * \param *b destination plan weight array
 * \param j index of destination element
 */
static void CopyElementOnepw(void *a, int i, void *b, int j) {
  ((struct plan_weight *)a)[i].Index = ((struct plan_weight *)b)[j].Index;
  ((struct plan_weight *)a)[i].weight = ((struct plan_weight *)b)[j].weight;
  ((struct plan_weight *)a)[i].PE = ((struct plan_weight *)b)[j].PE;
}

/** Creates a 3d plan for the FFT to use in transforms.
 * \sa fft_plan_3d
 * \param *P Problem at hand
 * \param comm MPI Communicator
 * \param E_cut Energy cutoff
 * \param MaxLevel number of lattice levels
 * \param *LevelSizes size increas factor from level to level for MaxLevel's
 * \param GBasisSQ[] scalar product of reciprocal basis vectors
 * \param N[] mesh points per dimension NDIM 
 * \param *MaxNoOfnFields maximum number of NFields per level
 * \param **NFields NField per level per field
 */
struct fft_plan_3d *fft_3d_create_plan(struct Problem *P, MPI_Comm comm, const double E_cut, const int MaxLevel, const int *LevelSizes, const double GBasisSQ[NDIM], const int N[NDIM], const int *MaxNoOfnFields, int **NFields) {
  struct fft_plan_3d *plan;
  int *MaxNFields;
  int i;
  int j;
  int ny;
  int Index;
  int LevelSize=1;      // factorial level increase from 0th to maxlevel
  int n[NDIM];
  int Args[3];
  double weight;
  // Initialization of plan
  plan = (struct fft_plan_3d *)
    Malloc(sizeof(struct fft_plan_3d),"create_plan");
  plan->P = P;
  plan->local_data = NULL;
  plan->work = NULL;
  plan->MaxLevel = MaxLevel;
  plan->LevelSizes = (int *)
    Malloc(sizeof(int)*MaxLevel,"create_plan: MaxLevels");
  for (i=0; i < MaxLevel; i++) {
    LevelSize *= LevelSizes[i];
    plan->LevelSizes[i] = LevelSizes[i];
  }
  plan->E_cut = E_cut;
  plan->LevelSize = LevelSize;
  if (N[0] % LevelSize != 0 || N[1] % LevelSize != 0 || N[2] % LevelSize != 0)
    Error(SomeError,"create_plan: N[0] % LevelSize != 0 || N[1] % LevelSize != 0 || N[2] % LevelSize != 0");
  // MPI part of plan
  MPI_Comm_dup(comm, &plan->comm);	       // Duplicates an existing communicator with all its cached information
  MPI_Comm_size(plan->comm, &plan->MaxPE); // Determines the size of the group associated with a communictor
  MPI_Comm_rank(plan->comm, &plan->myPE);  // Determines the rank of the calling process in the communicator
  plan->PENo = (int *)
    Malloc(2*sizeof(int)*plan->MaxPE,"create_plan: PENo");	// Malloc 2 ints per process in communicator group
  for (i=0; i < 2*plan->MaxPE; i++)	// all set to zero
    plan->PENo[i] = 0;

  for (i=0; i < NDIM; i++) {
    plan->N[i] = N[i];
    plan->GBasisSQ[i] = GBasisSQ[i];
    plan->NLSR[i] = N[i] / LevelSize;
  }
  plan->Maxpw = (N[1] / LevelSize) * ((N[2] / LevelSize)/2+1);
  if (P->Call.out[LeaderOut]) fprintf(stderr,"Create_Plan: E_cut %g, N[] %i %i %i, LS %i, Maxpw %i\n",E_cut, N[0],N[1],N[2], LevelSize, plan->Maxpw);
  if (plan->NLSR[0] % plan->MaxPE != 0 || plan->Maxpw % plan->MaxPE != 0) // check if dividable among processes
    Error(SomeError,"create_plan: plan->NLSR[0] % plan->MaxPE != 0 || plan->Maxpw % plan->MaxPE != 0");
  plan->pw = (struct plan_weight *)
    Malloc(plan->Maxpw*sizeof(struct plan_weight),"create_plan_weight"); // Allocating plan weight
  plan->pwFS = (struct plan_weight *)
    Malloc((plan->Maxpw+1)*sizeof(struct plan_weight),"create_plan_weight_FS");
  for (n[1]=0; n[1] < plan->NLSR[1]; n[1]++) {
    for (n[2]=0; n[2] < (plan->NLSR[2]/2+1); n[2]++) {
      ny = (n[1] >= (plan->NLSR[1]/2+1) ? n[1]-plan->NLSR[1] : n[1]);
      Index = CalcRowMajor2D(n[1],n[2],plan->NLSR[1],plan->NLSR[2]/2+1);
      plan->pw[Index].Index = Index;
      //fprintf(stderr,"(%i) plan->pw[%i] = %i\n",P->Par.me,Index, plan->pw[Index].Index);
      if (GBasisSQ[0] < MYEPSILON) fprintf(stderr,"fft_3d_create_plan: GBasisSQ[0] = %lg\n",GBasisSQ[0]);
      weight = (2.*E_cut/(LevelSize*LevelSize)-ny*ny*GBasisSQ[1]-n[2]*n[2]*GBasisSQ[2])/GBasisSQ[0];
      if (weight > 0.0)
	      plan->pw[Index].weight = sqrt(weight);
      else
	      plan->pw[Index].weight = 0.0;
      if (n[2] == plan->NLSR[2]/2) plan->pw[Index].weight += 4.*sqrt(2.*E_cut/(LevelSize*LevelSize)/GBasisSQ[0]);
      if (n[2] == 0) plan->pw[Index].weight += 2.*sqrt(2.*E_cut/(LevelSize*LevelSize)/GBasisSQ[0]);
    }
  }
  naturalmergesort(plan->pw,plan->pwFS,0,plan->Maxpw-1,&GetKeyOnepw,NULL,&CopyElementOnepw);
  if (plan->Maxpw < plan->MaxPE) Error(SomeError,"create_plan: plan->Maxpw < plan->MaxPE");
  if (plan->Maxpw / plan->MaxPE <  plan->NLSR[1] + plan->NLSR[1]/plan->MaxPE)
    Error(SomeError,"create_plan: plan->Maxpw / plan->MaxPE <  plan->NLSR[1] + plan->NLSR[1]/plan->MaxPE");
  for (i=0; i < plan->Maxpw; i++) {
    if (i < plan->NLSR[1]) {
      plan->pw[i].PE = i % plan->MaxPE;
    } else {
      if (i-plan->NLSR[1] < plan->NLSR[1]) {
      	plan->pw[i].PE = 0;
      } else {
      	if (i-2*plan->NLSR[1] < (plan->MaxPE-1)*plan->NLSR[1]) {
      	  plan->pw[i].PE = ((i%(plan->MaxPE-1))+1);
      	} else {
      	  plan->pw[i].PE = i%plan->MaxPE;
      	}
      }
    }
    plan->PENo[2*(i % plan->MaxPE)]++;
    if (plan->pw[i].Index%(plan->NLSR[2]/2+1) == plan->NLSR[2]/2)
      plan->pw[i].weight -= 4.*sqrt(2.*E_cut/(LevelSize*LevelSize)/GBasisSQ[0]);
    if (plan->pw[i].Index%(plan->NLSR[2]/2+1) == 0) {
      plan->pw[i].weight -= 2.*sqrt(2.*E_cut/(LevelSize*LevelSize)/GBasisSQ[0]);
      if (plan->pw[i].Index%(plan->NLSR[2]/2+1) >= plan->NLSR[2]/2+1) {
      	plan->pw[i].weight = 0.0;
      } else {
      	if (plan->pw[i].Index % (plan->NLSR[2]/2+1) == 0) {
      	  plan->PENo[2*(plan->pw[i].PE)+1] += floor(plan->pw[i].weight)+1;
      	} else {
      	  plan->PENo[2*(plan->pw[i].PE)+1] += 2*floor(plan->pw[i].weight)+1;
      	}
      }
    } else {
      plan->PENo[2*(plan->pw[i].PE)+1] += 2*floor(plan->pw[i].weight)+1;
    }
  }
  plan->LocalMaxpw = plan->PENo[2*plan->myPE];
  Args[0]=plan->MaxPE; Args[1]=plan->NLSR[2]/2+1; Args[2]=plan->Maxpw;
  naturalmergesort(plan->pw,plan->pwFS,0,plan->Maxpw-1,&GetKeyOnepw2PZI,Args,&CopyElementOnepw);
  for (j=0; j < plan->MaxPE; j++) {
    for (i= plan->NLSR[1]/plan->MaxPE-1; i >= 0; i--) {
      /*if (i != 0 && i+plan->LocalMaxpw*j >= i*(plan->NLSR[1]/2+1)+plan->LocalMaxpw*j)
      	Error(SomeError,"Createplan: i != 0 && i+plan->LocalMaxpw*j ? i*(plan->NLSR[2]/2+1)+plan->LocalMaxpw*j");*/
      CopyElementOnepw(plan->pwFS,plan->LocalMaxpw,plan->pw,i+plan->LocalMaxpw*j);
      CopyElementOnepw(plan->pw,i+plan->LocalMaxpw*j,plan->pw,(i+1)*(plan->NLSR[2]/2+1)-1+plan->LocalMaxpw*j);
      CopyElementOnepw(plan->pw,(i+1)*(plan->NLSR[2]/2+1)-1+plan->LocalMaxpw*j,plan->pwFS,plan->LocalMaxpw);
    }
  }
  plan->Localpw = &plan->pw[plan->myPE*plan->LocalMaxpw];
  if (P->Call.out[PsiOut]) {
    fprintf(stderr,"pw\n");
    for (i=plan->myPE*plan->LocalMaxpw; i < (plan->myPE+1)*plan->LocalMaxpw; i++) // print my local weights
      fprintf(stderr,"(%i)W(%g)I(%i)P(%i) ",i,plan->pw[i].weight,plan->pw[i].Index,plan->pw[i].PE);
    fprintf(stderr,"\n");
  }
  if (P->Call.out[LeaderOut]) {
     fprintf(stderr,"pwPENo\n");
    for (i=0; i < plan->MaxPE; i++)
      fprintf(stderr,"(%i)W(%i)X(%i) ",i,plan->PENo[2*i+1],plan->PENo[2*i]);
    fprintf(stderr,"\n");
  }
  Free(plan->pwFS);
  plan->pwFS = NULL;

  plan->Lev_CR_RC = (int *)
    Malloc(2*sizeof(int)*MaxLevel,"create_plan: Lev_CR_RC");
  for (i=0; i < 2*MaxLevel; i++) plan->Lev_CR_RC[i] = 1;
  plan->Levplan = (struct LevelPlan *)
    Malloc(sizeof(struct LevelPlan)*MaxLevel,"create_plan:Levplan");
  plan->MaxNoOfnFields = (int *)
    Malloc(sizeof(int)*MaxLevel,"create_plan:MaxNoOfnFields");
  for (i=0; i<MaxLevel;i++)
    plan->MaxNoOfnFields[i] = MaxNoOfnFields[i];
  plan->NFields = (int **)
    Malloc(sizeof(int *)*MaxLevel,"create_plan:NFields");
  MaxNFields = (int *)
    Malloc(sizeof(int)*MaxLevel,"create_plan:MaxNFields");
  for (i=0; i<MaxLevel;i++) {
    plan->NFields[i] = (int *)
      Malloc(sizeof(int )*MaxNoOfnFields[i],"create_plan:NFields");
    MaxNFields[i] = 0;
    for (j=0; j <  plan->MaxNoOfnFields[i]; j++) {
      plan->NFields[i][j] = NFields[i][j];
      if (NFields[i][j] < 1) Error(SomeError,"Create_plan: NFields[j] < 1");
      if (NFields[i][j] > MaxNFields[i]) MaxNFields[i] = NFields[i][j];
    }
  }
  plan->LevelBlockSize = plan->LevelSize*plan->LevelSize*plan->LevelSize;
  for (i=MaxLevel-1; i >= 0; i--) {
    if (i == MaxLevel-1)
      plan->Levplan[i].LevelFactor = LevelSizes[MaxLevel-1];
    else
      plan->Levplan[i].LevelFactor = plan->Levplan[i+1].LevelFactor*LevelSizes[i];
    for (j=0; j < NDIM; j++) 
      plan->Levplan[i].N[j] = plan->NLSR[j]*plan->Levplan[i].LevelFactor;
    plan->Levplan[i].LevelNo = i;
    plan->LocalDataSize = (plan->Levplan[i].N[0]*plan->Levplan[i].N[1]*(plan->Levplan[i].N[2]/2+1)*MaxNFields[i])/plan->MaxPE;
    plan->LocalWorkSize = plan->LocalDataSize;
    plan->local_data = (fftw_complex *)
      Realloc(plan->local_data, sizeof(fftw_complex)*plan->LocalDataSize,"create_plan: localdata");
    plan->work = (fftw_complex *)
      Realloc(plan->work, sizeof(fftw_complex)*plan->LocalWorkSize,"create_plan: work");
    plan->Levplan[i].ny = plan->Levplan[i].N[1];
    plan->Levplan[i].local_ny = plan->Levplan[i].N[1]/plan->MaxPE;
    plan->Levplan[i].nz = plan->Levplan[i].N[2]/2+1;
    plan->Levplan[i].local_nx = plan->Levplan[i].N[0]/plan->MaxPE;
    plan->Levplan[i].nx = plan->Levplan[i].N[0];
    plan->Levplan[i].start_nx =  plan->Levplan[i].local_nx * plan->myPE;
    plan->Levplan[i].local_nyz = plan->Levplan[i].local_ny*plan->Levplan[i].nz;
    plan->Levplan[i].nyz = plan->Levplan[i].ny*plan->Levplan[i].nz;
    fprintf(stderr,"(%i) nlocal (x%i, y%i, yz%i), n (%i, %i, %i), N (%i, %i, %i)\n", P->Par.me, plan->Levplan[i].local_nx, plan->Levplan[i].local_ny, plan->Levplan[i].local_nyz, plan->Levplan[i].nx, plan->Levplan[i].ny, plan->Levplan[i].nz, plan->Levplan[i].N[0], plan->Levplan[i].N[1], plan->Levplan[i].N[2]);  
    plan->Levplan[i].LocalSizeC = plan->Levplan[i].nx*plan->Levplan[i].local_nyz;
    plan->Levplan[i].LocalSizeR = plan->Levplan[i].local_nx*plan->Levplan[i].ny*plan->Levplan[i].N[2];
    plan->Levplan[i].LevelBlockSize = plan->Levplan[i].LevelFactor* plan->Levplan[i].LevelFactor*plan->Levplan[i].LevelFactor;
    plan->Levplan[i].SendRecvBlockSize = plan->Levplan[i].local_nyz*plan->Levplan[i].local_nx;
    plan->Levplan[i].PermBlockSize = plan->Levplan[i].LevelFactor;
    if (P->Call.out[LeaderOut])
      fprintf(stderr,"Create_Plan: Level(%i) LevelFactor(%i) LevelBlockSize(%i)\n",i,plan->Levplan[i].LevelFactor, plan->Levplan[i].LevelBlockSize);
    if (plan->Lev_CR_RC[2*i]) {
      plan->Levplan[i].xplanCR = (fftw_plan *)
	      Malloc(sizeof(fftw_plan)*plan->MaxNoOfnFields[i],"Create_Plan: xplanCR");
      plan->Levplan[i].yzplanCR = (rfftwnd_plan *)
	      Malloc(sizeof(rfftwnd_plan)*plan->MaxNoOfnFields[i],"Create_Plan: yzplanCR");
      for (j=0; j<plan->MaxNoOfnFields[i]; j++) {
      	plan->Levplan[i].xplanCR[j] = fftw_create_plan_specific(plan->Levplan[i].N[0], FFTW_BACKWARD, fftw_flag | FFTW_OUT_OF_PLACE,plan->local_data,plan->NFields[i][j],plan->work,plan->NFields[i][j]*plan->Levplan[i].local_nyz);
      	plan->Levplan[i].yzplanCR[j] = rfftw2d_create_plan_specific(plan->Levplan[i].N[1], plan->Levplan[i].N[2], FFTW_COMPLEX_TO_REAL, fftw_flag | FFTW_OUT_OF_PLACE,(fftw_real *)plan->local_data,plan->NFields[i][j],(fftw_real *)plan->work,plan->NFields[i][j]);
      }
    } else {
      plan->Levplan[i].xplanCR = NULL;
      plan->Levplan[i].yzplanCR = NULL;
    }
    if (plan->Lev_CR_RC[2*i+1]) {
      plan->Levplan[i].xplanRC = (fftw_plan *)
      	Malloc(sizeof(fftw_plan)*plan->MaxNoOfnFields[i],"Create_Plan: xplanRC");
      plan->Levplan[i].yzplanRC = (rfftwnd_plan *)
      	Malloc(sizeof(rfftwnd_plan)*plan->MaxNoOfnFields[i],"Create_Plan: yzplanRC");
      for (j=0; j<plan->MaxNoOfnFields[i]; j++) {
      	plan->Levplan[i].xplanRC[j] = fftw_create_plan_specific(plan->Levplan[i].N[0], FFTW_FORWARD, fftw_flag | FFTW_OUT_OF_PLACE,plan->local_data,plan->NFields[i][j]*plan->Levplan[i].local_nyz,plan->work,plan->NFields[i][j]);
      	plan->Levplan[i].yzplanRC[j] = rfftw2d_create_plan_specific(plan->Levplan[i].N[1], plan->Levplan[i].N[2], FFTW_REAL_TO_COMPLEX, fftw_flag | FFTW_OUT_OF_PLACE,(fftw_real *)plan->local_data,plan->NFields[i][j],(fftw_real *)plan->work,plan->NFields[i][j]);
      }
    } else {
      plan->Levplan[i].xplanRC = NULL;
      plan->Levplan[i].yzplanRC = NULL;
    }
  }
  if (plan->local_data != NULL) {
    Free(plan->local_data);
    plan->local_data = NULL;
  }
  if (plan->work != NULL) {
    Free(plan->work);
    plan->work = NULL;
  }
  Free(MaxNFields);
  return plan;
}

void fft_3d_destroy_plan(const struct Problem *P, struct fft_plan_3d *plan) {
  int i,j;
  //if (P->Call.out[LeaderOut]) fprintf(stderr,"Destroy_Plan\n");
  for (i=0; i< plan->MaxLevel; i++) {
    for (j=0; j< plan->MaxNoOfnFields[i]; j++) {
      fftw_destroy_plan(plan->Levplan[i].xplanCR[j]);
      fftw_destroy_plan(plan->Levplan[i].xplanRC[j]);
      rfftwnd_destroy_plan(plan->Levplan[i].yzplanCR[j]);
      rfftwnd_destroy_plan(plan->Levplan[i].yzplanRC[j]);
    }
    Free(plan->Levplan[i].xplanCR);
    Free(plan->Levplan[i].xplanRC);
    Free(plan->Levplan[i].yzplanCR);
    Free(plan->Levplan[i].yzplanRC);
    Free(plan->NFields[i]);
  }
  Free(plan->MaxNoOfnFields);
  Free(plan->NFields);
  Free(plan->Levplan);
  Free(plan->PENo);
  plan->PENo = NULL;
  Free(plan->LevelSizes);
  plan->LevelSizes = NULL;
  Free(plan->Lev_CR_RC);
  plan->Lev_CR_RC = NULL;
  Free(plan->pw);
  plan->pw = NULL;
  plan->Levplan = NULL;
  MPI_Comm_free(&plan->comm);
  Free(plan);
}

static void FirstDimTransCR(const struct fft_plan_3d *plan,  const struct LevelPlan *Lplan, const int nFieldsNo, fftw_complex *local_data, fftw_complex *work) { /* local_data -> work */
  int local_nyz = Lplan->local_nyz;
  int nx = Lplan->N[0];
  int fft_iter;
  int nFields = plan->NFields[Lplan->LevelNo][nFieldsNo];
  fftw_plan fft_x = Lplan->xplanCR[nFieldsNo];
  if (nFields > 1) {
    for (fft_iter = 0; fft_iter < local_nyz; ++fft_iter)
      fftw(fft_x, nFields,
	   local_data + (nx * nFields) * fft_iter, nFields, 1,
	   work+nFields*fft_iter, nFields*local_nyz, 1);
  } else
    fftw(fft_x, local_nyz,
	 local_data, 1, nx,
	 work, local_nyz, 1);
}

static void OtherDimTransCR(const struct fft_plan_3d *plan, const struct LevelPlan *Lplan, const int nFieldsNo, fftw_complex *local_data, fftw_complex *work) { /* work -> local_data */
  int local_nx = Lplan->local_nx;
  int nyz = Lplan->nyz;
  int fft_iter;
  int rnyz = Lplan->N[1]*Lplan->N[2];
  int nFields = plan->NFields[Lplan->LevelNo][nFieldsNo];
  rfftwnd_plan fft = Lplan->yzplanCR[nFieldsNo];
  if (nFields > 1) {
    for (fft_iter = 0; fft_iter < local_nx; ++fft_iter)
      rfftwnd_complex_to_real(fft, nFields,
			      ((fftw_complex *) work)+ (nyz * nFields) * fft_iter,
			      nFields, 1,
			      ((fftw_real *)local_data)+ (rnyz * nFields) * fft_iter
			      , nFields, 1);
  } else {
    rfftwnd_complex_to_real(fft, local_nx,
			    (fftw_complex *) work,
			    1, nyz,
			    (fftw_real *)local_data
			    , 1, rnyz);
  }
}

static void TransposeMPICR(const struct fft_plan_3d *plan, const struct LevelPlan *Lplan, const int el_size, double *local_data,  double *work) { /* work -> local_data */
  MPI_Alltoall(work, Lplan->SendRecvBlockSize * el_size, MPI_DOUBLE,
	       local_data, Lplan->SendRecvBlockSize * el_size, MPI_DOUBLE,
	       plan->comm);
}

static void LocalPermutationCR(const struct fft_plan_3d *plan, const struct LevelPlan *Lplan, const int el_size, double *local_data,  double *work) { /* data -> work */
  int x,PEy,ly,z,Z,Y,srcIndex,destIndex,Index,cpyes,Ly;
  int lnx=Lplan->local_nx,lny=plan->NLSR[1]/plan->MaxPE,nz=plan->NLSR[2]/2+1,MaxPE=plan->MaxPE,ny=plan->NLSR[1];
  int es = el_size*Lplan->PermBlockSize;
  int LevelFactor = Lplan->LevelFactor;
  int snz = (nz-1)*es+el_size;
  struct plan_weight *pw = plan->pw;
  for (x=0; x < lnx; x++)
    for (PEy=0; PEy < MaxPE; PEy++)
      for (ly=0; ly < lny; ly++)
	for (Ly=0; Ly < LevelFactor;Ly++)
	  for (z=0; z < nz; z++) {
	    Z = z*es;
	    srcIndex = Z+snz*(Ly+LevelFactor*(ly+lny*(x+lnx*(PEy))));
	    Index = pw[z+nz*(ly+lny*PEy)].Index;
	    Z = Index%nz;
	    Z = Z*es;
	    Y = ((Index/nz)*LevelFactor+Ly);
	    destIndex = Z+snz*(Y+LevelFactor*ny*(x));
	    cpyes = (z == nz-1 ? el_size : es);
	    memcpy( &work[destIndex],
		    &local_data[srcIndex],
		    cpyes * sizeof(double));
	  }
}

// does complex to real 3d fftransformation (reciprocal base -> real base)
/** Inverse Fast Fourier Transformation of plane wave coefficients.
 * If given state is \f$\langle \chi_{i,G} | \psi_i (r) \rangle \f$, result is \f$| \psi_i (r) \rangle\f$
 * Calls subsequently (reverse order of calling as in fft_3d_real_to_complex())
 * -# FirstDimTransCR()\n
 * -# TransposeMPICR()\n
 * -# LocalPermutationCR()\n
 * -# OtherDimTransCR()\n
 * \param *plan the transformation plan
 * \param Level current level number
 * \param nFieldsNo number of parallel ffts (in case of LevelUp)
 * \param *local_data real base wave function coefficients, \f$\psi_i (G)\f$, on return \f$\psi_i (r)\f$
 * \param *work temporary array for coefficients
 * \warning coefficients in \a *local_date and *work are destroyed!
 * \note Due to the symmetry of the reciprocal coefficients, not all of them are stored in memory. Thus
 *       the whole set must be rebuild before this inverse transformation can be called. From 0 to
 *       LatticeLevel::MaxG copied to negative counterpart and for 0 to LatticeLevel::MaxDoubleG the
 *       complex conjugate must be formed (see \ref density.c).
 */
void fft_3d_complex_to_real(const struct fft_plan_3d *plan, const int Level, const int nFieldsNo, fftw_complex *local_data, fftw_complex *work) { /* local data -> local_data, work destroyed */
  int el_size = (sizeof(fftw_complex) / sizeof(double));
  int nFields = plan->NFields[Level][nFieldsNo];
  struct LevelPlan *Lplan = &plan->Levplan[Level];
  if (!plan->Lev_CR_RC[2*Level]) Error(SomeError,"fft_3d_complex_to_real: !plan->Lev_CR_RC[2*i]");
  if (nFieldsNo < 0 || nFieldsNo >= plan->MaxNoOfnFields[Level]) Error(SomeError,"fft_3d_complex_to_real: nFieldsNo < 0 || nFieldsNo >= plan->MaxNoOfnFields[Level]");
  el_size *= nFields;
  //if (isnan(work[0].re)) 
  //fprintf(stderr,"fft_3d_complex_to_real,FirstDimTransCR: before work %lg\n",work[0].re);
  FirstDimTransCR(plan, Lplan, nFieldsNo, local_data, work);
  //if (isnan(work[0].re)) 
  //fprintf(stderr,"fft_3d_complex_to_real,FirstDimTransCR: after work %lg\n",work[0].re);
  TransposeMPICR(plan, Lplan, el_size, (double *)local_data, (double *)work);
  LocalPermutationCR(plan, Lplan, el_size, (double *)local_data, (double *)work);
  OtherDimTransCR(plan, Lplan, nFieldsNo, local_data, work);
}

static void FirstDimTransRC(const struct fft_plan_3d *plan, const struct LevelPlan *Lplan, const int nFieldsNo, fftw_complex *local_data, fftw_complex *work) { /* work -> local_data */
  int local_nyz = Lplan->local_nyz;
  int nx = Lplan->N[0];
  int fft_iter;
  int nFields = plan->NFields[Lplan->LevelNo][nFieldsNo];
  fftw_plan fft_x = Lplan->xplanRC[nFieldsNo];
  if (nFields > 1) {
    for (fft_iter = 0; fft_iter < local_nyz; ++fft_iter)
      fftw(fft_x, nFields,
	   work+nFields*fft_iter, nFields*local_nyz, 1,
	   local_data + (nx * nFields) * fft_iter, nFields, 1);
  } else
    fftw(fft_x, local_nyz,
	 work, local_nyz, 1,
	 local_data, 1, nx);
}

static void OtherDimTransRC(const struct fft_plan_3d *plan, const struct LevelPlan *Lplan, const int nFieldsNo, fftw_complex *local_data, fftw_complex *work) { /* local_data -> work */
  int local_nx = Lplan->local_nx;
  int nyz = Lplan->nyz;
  int rnyz = Lplan->N[1]*Lplan->N[2];
  int fft_iter;
  int nFields = plan->NFields[Lplan->LevelNo][nFieldsNo];
  rfftwnd_plan fft = Lplan->yzplanRC[nFieldsNo];
  if (nFields > 1) {
    for (fft_iter = 0; fft_iter < local_nx; ++fft_iter)
      rfftwnd_real_to_complex(fft, nFields,
			      ((fftw_real *) local_data) + (rnyz * nFields) * fft_iter,
			      nFields, 1,
			      ((fftw_complex *)work)+ (nyz * nFields) * fft_iter
			      , nFields, 1);
  } else {
    rfftwnd_real_to_complex(fft, local_nx,
			    (fftw_real *) local_data,
			    1, rnyz,
			    (fftw_complex *) work
			    , 1, nyz);
  }
}

static void TransposeMPIRC(const struct fft_plan_3d *plan, const struct LevelPlan *Lplan, const int el_size, double *local_data,  double *work) { /* local_data -> work */
  MPI_Alltoall(local_data, Lplan->SendRecvBlockSize * el_size, MPI_DOUBLE,
	       work, Lplan->SendRecvBlockSize * el_size, MPI_DOUBLE,
	       plan->comm);
}

static void LocalPermutationRC(const struct fft_plan_3d *plan, const struct LevelPlan *Lplan, const int el_size, double *local_data,  double *work) { /* work -> data*/
  int x,PEy,ly,z,Z,Y,srcIndex,destIndex,Index,cpyes,Ly;
  int lnx=Lplan->local_nx,lny=plan->NLSR[1]/plan->MaxPE,nz=plan->NLSR[2]/2+1,MaxPE=plan->MaxPE,ny=plan->NLSR[1];
  int es = el_size*Lplan->PermBlockSize;
  int LevelFactor = Lplan->LevelFactor;
  int snz = (nz-1)*es+el_size;
  struct plan_weight *pw = plan->pw;
  for (x=0; x < lnx; x++)
    for (PEy=0; PEy < MaxPE; PEy++)
      for (ly=0; ly < lny; ly++)
	for (Ly=0; Ly < LevelFactor;Ly++)
	  for (z=0; z < nz; z++) {
	    Z = z*es;
	    destIndex = Z+snz*(Ly+LevelFactor*(ly+lny*(x+lnx*(PEy))));
	    Index = pw[z+nz*(ly+lny*PEy)].Index;
	    Z = Index%nz;
	    Z = Z*es;
	    Y = ((Index/nz)*LevelFactor+Ly);
	    srcIndex = Z+snz*(Y+LevelFactor*ny*(x));
	    cpyes = (z == nz-1 ? el_size : es);
	    memcpy( &local_data[destIndex],
		    &work[srcIndex],
		    cpyes * sizeof(double));
	  }
}

/** Normal Fast Fourier Transformation of plane wave coefficients.
 * If given state is \f$| \psi_i (r) \rangle \f$, result is \f$\langle \chi_{i,G} | \psi_i (r) \rangle\f$
 * Calls subsequently
 * -# OtherDimTransRC()\n
 * -# LocalPermutationRC()\n
 * -# TransposeMPIRC()\n
 * -# FirstDimTransRC()\n
 * \param *plan the transformation plan
 * \param Level current level number
 * \param nFieldsNo number of parallel ffts (in case of LevelUp)
 * \param *local_data real base wave function coefficients, \f$\psi_i (r)\f$, on return \f$\psi_i (G)\f$
 * \param *work temporary array for coefficients
 * \warning coefficients in \a *local_date and *work are destroyed!
 */
void fft_3d_real_to_complex(const struct fft_plan_3d *plan, const int Level, const int nFieldsNo, fftw_complex *local_data, fftw_complex *work) { /* local data -> local_data, work destroyed */
  int el_size = (sizeof(fftw_complex) / sizeof(double));
  int nFields = plan->NFields[Level][nFieldsNo];
  struct LevelPlan *Lplan = &plan->Levplan[Level];
  if (!plan->Lev_CR_RC[2*Level+1]) Error(SomeError,"fft_3d_real_to_complex: !plan->Lev_CR_RC[2*i+1]");
  if (nFieldsNo < 0 || nFieldsNo >= plan->MaxNoOfnFields[Level]) Error(SomeError,"fft_3d_real_to_complex: nFieldsNo < 0 || nFieldsNo >= plan->MaxNoOfnFields[Level]");
  el_size *= nFields;
  OtherDimTransRC(plan, Lplan, nFieldsNo, local_data, work);
  LocalPermutationRC(plan, Lplan, el_size, (double *)local_data, (double *)work);
  TransposeMPIRC(plan, Lplan, el_size, (double *)local_data, (double *)work);
  FirstDimTransRC(plan, Lplan, nFieldsNo, local_data, work);
}

/** Fills \a *local_data with random values.
 * \param *plan fft plan
 * \param Level current level number
 * \param nFields number of nFields (of parallel ffts)
 * \param *local_data array of coefficients to be set at random values
 */
void InitDataTestR(const struct fft_plan_3d *plan, const int Level, const int nFields,fftw_real *local_data) {
  int x,y,z,n,Index,i;
  int Nx = plan->Levplan[Level].N[0]/plan->MaxPE;
  int Ny = plan->Levplan[Level].N[1];
  int Nz = plan->Levplan[Level].N[2];
  srand((unsigned int)Level);
  for (i=0;i<plan->myPE*Nx*Ny*Nz*nFields;i++) rand();
  for (x=0; x < Nx; x++) {
    for (y = 0; y < Ny; y++) {
      for (z = 0; z < Nz; z++) {
      	for (n=0; n < nFields; n++) {
      	  Index = n+nFields*(z+Nz*(y+Ny*(x)));
      	  local_data[Index]=1.-2.*(rand()/(double)RAND_MAX);
      	}
      }
    }
  }
}

/*
static void OutputDataTestWR(const struct fft_plan_3d *plan, const int Level, const int nFields,fftw_real *local_data,FILE *target) {
  int x,y,z,n,Index;
  int Nx = plan->Levplan[Level].N[0]/plan->MaxPE;
  int Ny = plan->Levplan[Level].N[1];
  int Nz = plan->Levplan[Level].N[2];
  fprintf(target,"Level=%i\n",Level);
  for (x=0; x < Nx; x++) {
    fprintf(target,"x=%i\n",x+Nx*plan->myPE);
    for (y = 0; y < Ny; y++) {
      fprintf(target,"y=%5i     ",y);
      for (z = 0; z < Nz; z++) {
	for (n=0; n < nFields; n++) {
	  Index = n+nFields*(z+2*(Nz/2+1)*(y+Ny*(x)));
	  fprintf(target,"%10.5f ",local_data[Index]);
	}
      }
      fprintf(target,"\n");
    }
  }
}

static void OutputDataTestR(const struct fft_plan_3d *plan, const int Level, const int nFields,fftw_real *local_data, double factor,FILE *target) {
  int x,y,z,n,Index;
  int Nx = plan->Levplan[Level].N[0]/plan->MaxPE;
  int Ny = plan->Levplan[Level].N[1];
  int Nz = plan->Levplan[Level].N[2];
  fprintf(target,"Level=%i\n",Level);
  for (x=0; x < Nx; x++) {
    fprintf(target,"x=%i\n",x+Nx*plan->myPE);
    for (y = 0; y < Ny; y++) {
      fprintf(target,"y=%5i     ",y);
      for (z = 0; z < Nz; z++) {
	for (n=0; n < nFields; n++) {
	  Index = n+nFields*(z+Nz*(y+Ny*(x)));
	  fprintf(target,"%10.5f ",local_data[Index]*factor);
	}
      }
      fprintf(target,"\n");
    }
  }
}

static void InitDataTestWR(const struct fft_plan_3d *plan, const int Level, const int nFields,fftw_real *local_data) {
  int i,x,y,z,n,Index,Value=0;
  int Nx = plan->Levplan[Level].N[0]/plan->MaxPE;
  int Ny = plan->Levplan[Level].N[1];
  int Nz = plan->Levplan[Level].N[2];
  srand(Level);
  for (i=0;i<plan->myPE*Nx*Ny*Nz*nFields;i++) rand();
  for (x=0; x < Nx; x++) {
    for (y = 0; y < Ny; y++) {
      for (z = 0; z < Nz; z++) {
	for (n=0; n < nFields; n++) {
	  Index = n+nFields*(z+2*(Nz/2+1)*(y+Ny*(x)));
	  local_data[Index] = 1.-2.*(rand()/(double)RAND_MAX);
	  Value++;
	}
      }
    }
  }
}

static void OutputDataTestWCA(const struct fft_plan_3d *plan, const int Level, const int nFields,fftw_complex *local_data,double factor,FILE *target) {
  int x,y,z,yz,n,Index,Y,Z,YZ;
  int Nx = plan->Levplan[Level].N[0];
  int Ny = plan->Levplan[Level].N[1]/plan->MaxPE;
  int Nz = plan->Levplan[Level].N[2]/2+1;
  int NZ = plan->NLSR[2]/2+1;
  fprintf(target,"Level=%i\n",Level);
  for (y = 0; y < Ny; y++) {
    for (z = 0; z < Nz; z++) {
      Y = (y+plan->myPE*Ny)/plan->Levplan[Level].LevelFactor;
      Z = z/plan->Levplan[Level].LevelFactor;
      YZ = Z+NZ*Y;
      yz = z+Nz*y;
      fprintf(target,"yz=%i(%i)\n",yz+plan->myPE*Ny*Nz,YZ);
      for (x=0; x < Nx; x++) {
	for (n=0; n < nFields; n++) {
	  Index = n+nFields*(z+Nz*(x+Nx*(y)));
	  fprintf(target,"%10.5f %10.5f ",local_data[Index].re*factor,local_data[Index].im*factor);
	}
      }
      fprintf(target,"\n");
    }

  }
}

static void OutputDataTestCA(const struct fft_plan_3d *plan, const int Level, const int nFields,fftw_complex *local_data,double factor,FILE *target) {
  int x,yz,n,Index,Y,Z,YZ,y,z;
  int Nx = plan->Levplan[Level].nx;
  int Ny = plan->Levplan[Level].local_ny ;
  int Nz = plan->Levplan[Level].nz;
  int Nyz = Ny*Nz;
  int NZ = plan->NLSR[2]/2+1;
  fprintf(target,"Level=%i\n",Level);
  for (y = 0; y < Ny; y++) {
    for (z = 0; z < Nz; z++) {
      Y = (y+plan->myPE*Ny)/plan->Levplan[Level].LevelFactor;
      Z = z/plan->Levplan[Level].LevelFactor;
      YZ = Z+NZ*Y;
      YZ = plan->pw[YZ].Index;
      yz = z+Nz*y;
      fprintf(target,"yz=%i(%i)\n",yz+plan->myPE*Nyz,YZ);
      for (x=0; x < Nx; x++) {
	for (n=0; n < nFields; n++) {
	  Index = n+nFields*(x+Nx*(yz));
	  fprintf(target,"%10.5f %10.5f ",((double *)local_data)[2*Index]*factor,((double *)local_data)[2*Index+1]*factor);
	}
      }
      fprintf(target,"\n");
    }
  }
}

static void InitDataTestCA(const struct fft_plan_3d *plan, const int Level, const int nFields,fftw_complex *local_data) {
  int x,yz,n,Index,Y,Z,YZ,y,z,Value=0;
  int Nx = plan->Levplan[Level].nx;
  int Ny = plan->Levplan[Level].local_ny ;
  int Nz = plan->Levplan[Level].nz;
  int Nyz = plan->Levplan[Level].ny*Nz;
  int NZ = plan->NLSR[2]/2+1;
  int YY;
  for (y = 0; y < Ny; y++) {
    for (z = 0; z < Nz; z++) {
      Y = (y+plan->myPE*Ny)/plan->Levplan[Level].LevelFactor;
      YY = (y+plan->myPE*Ny)%plan->Levplan[Level].LevelFactor;
      Z = z/plan->Levplan[Level].Level;
      YZ = Z+NZ*Y;
      YZ = plan->pw[YZ].Index;
      Z = (YZ%NZ)*plan->Levplan[Level].LevelFactor+(z)%plan->Levplan[Level].LevelFactor;
      yz = Z+Nz*(((YZ/NZ)*plan->Levplan[Level].LevelFactor+YY));
      for (x=0; x < Nx; x++) {
	for (n=0; n < nFields; n++) {
	  Index = n+nFields*(x+Nx*(z+Nz*y));
	  Value = n+nFields*(yz+Nyz*x);
	  local_data[Index].re = Value*2;
	  local_data[Index].im = Value*2+1;
	}
      }
    }
  }
}

static void OutputDataTestCB(const struct fft_plan_3d *plan, const int Level, const int nFields,fftw_complex *local_data,double factor,FILE *target) {
  int x,yz,n,Index,Y,Z,YZ,y,z;
  int Nx = plan->Levplan[Level].local_nx;
  int Ny = plan->Levplan[Level].ny ;
  int Nz = plan->Levplan[Level].nz;
  int NZ = plan->NLSR[2]/2+1;
  fprintf(target,"Level=%i\n",Level);
  for (x=0; x < Nx; x++) {
    fprintf(target,"x=%i\n",x+plan->myPE*Nx);
    for (y = 0; y < Ny; y++) {
      for (z = 0; z < Nz; z++) {
	Y = (y)/plan->Levplan[Level].LevelFactor;
	Z = z/plan->Levplan[Level].LevelFactor;
	YZ = Z+NZ*Y;
	YZ = plan->pw[YZ].Index;
	yz = z+Nz*y;
	for (n=0; n < nFields; n++) {
	  Index = n+nFields*(yz+Ny*Nz*(x));
	  fprintf(target,"%10.5f %10.5f ",local_data[Index].re*factor,local_data[Index].im*factor);
	}
      }
      fprintf(target,"\n");
    }
  }
}

void fft_3d_c2r_r2c_Test(struct Problem *P, const struct fft_plan_3d *plan, const int nFields) {
  int i,j,Max=10;
  double time1, time2, timeCR, timeRC, time1CR, time1RC, factor =1.0;
  int local_nx, local_x_start,	local_ny_after_transpose,local_y_start_after_transpose,total_local_size;
  FILE *data0,*data1,*data0C,*data1C;
  OpenFileNo2(P,&data0,".data0",plan->myPE,"w",P->Call.out[ReadOut]);
  OpenFileNo2(P,&data1,".data1",plan->myPE,"w",P->Call.out[ReadOut]);
  OpenFileNo2(P,&data0C,".data0C",plan->myPE,"w",P->Call.out[ReadOut]);
  OpenFileNo2(P,&data1C,".data1C",plan->myPE,"w",P->Call.out[ReadOut]);
  fprintf(stderr, "(%i)C2RTest: \n", P->Par.me);
  for (i=plan->MaxLevel-1; i >= 0; i--) {
    timeCR = 0;
    timeRC = 0;
    time1RC = 0;
    time1CR = 0;
    fprintf(stderr,"PE(%i) L(%i) Nx(%i) Nyz(%i) nF(%i)\n",plan->myPE,i,plan->Levplan[i].N[0],plan->Levplan[i].local_nyz,nFields);
    factor = 1./(plan->Levplan[i].N[0]*plan->Levplan[i].N[1]*plan->Levplan[i].N[2]);
    for (j=0; j < Max; j++) {
      InitDataTestR(plan,i,nFields,(fftw_real *)plan->local_data);
      time1 = GetTime();
      fft_3d_real_to_complex(plan,i,nFields,plan->local_data, plan->work);
      time2 = GetTime();
      timeRC += (time2 - time1);
      time1 = GetTime();
      fft_3d_complex_to_real(plan,i,nFields,plan->local_data, plan->work);
      time2 = GetTime();
      timeCR += (time2 - time1);

    }
    fprintf(stderr,"PE(%i): L(%i) CR:sec(%g) RC:sec(%g)\n",plan->myPE, i, timeCR, timeRC);
  }
  fclose(data0);
  fclose(data1);
  fclose(data0C);
  fclose(data1C);
}
*/