/** \file riemann.c
 * RiemannTensor-transformations.
 * 
 * This is all still not really working, thus so far untouched.
 * 
  Project: ParallelCarParrinello
 \author Jan Hamaekers
 \date 2000

  File: riemann.c
  $Id: riemann.c,v 1.8 2006/03/30 22:19:52 foo Exp $
*/

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

#include"data.h"
#include"errors.h"
#include"helpers.h"
#include"mymath.h"
#include"riemann.h"
#include"myfft.h"

static void RiemannRTransformPosNFRto0(const struct RiemannTensor *RT, fftw_real *source, fftw_real *dest, const int NF)
{
  struct LatticeLevel *Lev = RT->LevR;
  int es = Lev->NUp0[2]*NF;
  unsigned int cpyes = sizeof(fftw_real)*es;
  int nx=Lev->Plan0.plan->local_nx,ny=Lev->N[1],nz=Lev->N[2],Nx=Lev->NUp0[0],Ny=Lev->NUp0[1];
  int lx,ly,z,Lx,Ly;
  for(lx=0; lx < nx; lx++)
    for(Lx=0; Lx < Nx; Lx++)
      for(ly=0; ly < ny; ly++)
	for(Ly=0; Ly < Ny; Ly++)
	  for(z=0; z < nz; z++) {
	    memcpy( &dest[es*(z+nz*(Ly+Ny*(ly+ny*(Lx+Nx*lx))))],
		    &source[es*(Ly+Ny*(Lx+Nx*(z+nz*(ly+ny*lx))))],
		    cpyes);
	  }
}

static void RiemannRTransformPosNFRtoS(const struct RiemannTensor *RT, fftw_real *source, fftw_real *dest, const int NF)
{
  struct LatticeLevel *Lev = RT->LevR;
  int es = Lev->NUp1[2]*NF;
  unsigned int cpyes = sizeof(fftw_real)*es;
  int nx=Lev->Plan0.plan->local_nx,ny=Lev->N[1],nz=Lev->N[2],Nx=Lev->NUp1[0],Ny=Lev->NUp1[1];
  int lx,ly,z,Lx,Ly;
  for(lx=0; lx < nx; lx++)
    for(Lx=0; Lx < Nx; Lx++)
      for(ly=0; ly < ny; ly++)
	for(Ly=0; Ly < Ny; Ly++)
	  for(z=0; z < nz; z++) {
	    memcpy( &dest[es*(z+nz*(Ly+Ny*(ly+ny*(Lx+Nx*lx))))],
		    &source[es*(Ly+Ny*(Lx+Nx*(z+nz*(ly+ny*lx))))],
		    cpyes);
	  }
}

static void RiemannRTransformPos(const struct LatticeLevel *Lev, fftw_real *source, fftw_real *dest, const int NF)
{
  int es = Lev->NUp[2]*NF;
  unsigned int cpyes = sizeof(fftw_real)*es;
  int nx=Lev->Plan0.plan->local_nx,ny=Lev->N[1],nz=Lev->N[2],Nx=Lev->NUp[0],Ny=Lev->NUp[1];
  int lx,ly,z,Lx,Ly;
  for(lx=0; lx < nx; lx++)
    for(Lx=0; Lx < Nx; Lx++)
      for(ly=0; ly < ny; ly++)
	for(Ly=0; Ly < Ny; Ly++)
	  for(z=0; z < nz; z++) {
	    memcpy( &dest[es*(z+nz*(Ly+Ny*(ly+ny*(Lx+Nx*lx))))],
		    &source[es*(Ly+Ny*(Lx+Nx*(z+nz*(ly+ny*lx))))],
		    cpyes);
	  }
}

static void CalculateLapiGcgS(const struct Lattice *Lat, const struct RiemannTensor *RT, fftw_complex *source) 
{
  int i,Index,j;
  struct  fft_plan_3d *plan = Lat->plan; 
  struct LatticeLevel *LevS = RT->LevS;
  fftw_complex *LapiGcgC = RT->DensityC[RTAiGcg];
  fftw_complex *LapcgC = RT->DensityC[RTAcg];
  fftw_complex *work = RT->TempC;
  fftw_complex *destCS, *destCD;
  int TotalSize = RT->TotalSize[RTAiGcg]/LevS->MaxNUp;
  double *G;
  SetArrayToDouble0((double *)LapiGcgC, TotalSize*2);
  SetArrayToDouble0((double *)LapcgC, (TotalSize/NDIM)*2);
  for (i=0;i<LevS->MaxG;i++) {
    Index = LevS->GArray[i].Index;
    LapcgC[Index].re = source[i].re;
    LapcgC[Index].im = source[i].im;
    G = LevS->GArray[i].G;
    for (j=0;j<NDIM;j++) { 
      LapiGcgC[Index*NDIM+j].re = -source[i].im*G[j];
      LapiGcgC[Index*NDIM+j].im = source[i].re*G[j];
    }
  }
  for (i=0; i<LevS->MaxDoubleG; i++) {
    destCS = &LapiGcgC[LevS->DoubleG[2*i]*NDIM];
    destCD = &LapiGcgC[LevS->DoubleG[2*i+1]*NDIM];
    for (j=0; j < NDIM; j++) {
      destCD[j].re =  destCS[j].re;
      destCD[j].im = -destCS[j].im;
    }
  }
  fft_3d_complex_to_real(plan, LevS->LevelNo, FFTNFSVec, LapiGcgC, work);
  for (i=0; i<LevS->MaxDoubleG; i++) {
    destCS = &LapcgC[LevS->DoubleG[2*i]];
    destCD = &LapcgC[LevS->DoubleG[2*i+1]];
    destCD[0].re =  destCS[0].re;
    destCD[0].im = -destCS[0].im;
  }
  fft_3d_complex_to_real(plan, LevS->LevelNo, FFTNF1, LapcgC, work);
}

static void CalculateLapiGcg0(const struct Lattice *Lat, const struct RiemannTensor *RT, fftw_complex *source/*, const double FactorC_R, const double FactorR_C*/) 
{
  int i,Index,j;
  struct  fft_plan_3d *plan = Lat->plan; 
  struct LatticeLevel *Lev0 = RT->Lev0;
  fftw_complex *LapiGcgC = RT->DensityC[RTAiGcg];
  fftw_complex *LapcgC = RT->DensityC[RTAcg];
  fftw_complex *work = RT->TempC;
  fftw_complex *destCS, *destCD;
  int TotalSize = RT->TotalSize[RTAiGcg];
  double *G;
  SetArrayToDouble0((double *)LapiGcgC, TotalSize*2);
  SetArrayToDouble0((double *)LapcgC, (TotalSize/NDIM)*2);
  for (i=0;i<Lev0->MaxG;i++) {
    Index = Lev0->GArray[i].Index;
    LapcgC[Index].re = source[i].re;
    LapcgC[Index].im = source[i].im;
    G = Lev0->GArray[i].G;
    for (j=0;j<NDIM;j++) { 
      LapiGcgC[Index*NDIM+j].re = -source[i].im*G[j];
      LapiGcgC[Index*NDIM+j].im = source[i].re*G[j];
    }
  }
  for (i=0; i<Lev0->MaxDoubleG; i++) {
    destCS = &LapiGcgC[Lev0->DoubleG[2*i]*NDIM];
    destCD = &LapiGcgC[Lev0->DoubleG[2*i+1]*NDIM];
    for (j=0; j < NDIM; j++) {
      destCD[j].re =  destCS[j].re;
      destCD[j].im = -destCS[j].im;
    }
  }
  fft_3d_complex_to_real(plan, Lev0->LevelNo, FFTNF0Vec, LapiGcgC, work);
  for (i=0; i<Lev0->MaxDoubleG; i++) {
    destCS = &LapcgC[Lev0->DoubleG[2*i]];
    destCD = &LapcgC[Lev0->DoubleG[2*i+1]];
    destCD[0].re =  destCS[0].re;
    destCD[0].im = -destCS[0].im;
  }
  fft_3d_complex_to_real(plan, Lev0->LevelNo, FFTNF1, LapcgC, work);
}

void CalculateRiemannLaplace0(const struct Lattice *Lat, struct RiemannTensor *RT, fftw_complex *src, fftw_complex *Lap/*, const double FactorC_R, const double FactorR_C*/) 
{
  int i,j,Index;
  struct  fft_plan_3d *plan = Lat->plan; 
  struct LatticeLevel *Lev0 = RT->Lev0;
  fftw_real **RTAR = RT->DensityR;
  const int *LocalSizeR = RT->LocalSizeR;
  const int *LocalSizeC = RT->LocalSizeC;
  fftw_complex *work = RT->TempC;
  fftw_real *LapiGcgR = RTAR[RTAiGcg], *LapcgR = RTAR[RTAcg];
  fftw_real *LapValR = RTAR[RTAARTA];
  fftw_real *LapVecR = RTAR[RTARTA];
  fftw_real *LapMatR = RTAR[RTAIRT];
  fftw_real *MulValR = RTAR[RTAcg];
  fftw_complex *MulValC = (fftw_complex *)MulValR;
  fftw_real *MulVecR = RTAR[RTAiGcg];
  fftw_complex *MulVecC = (fftw_complex *)MulVecR;
  fftw_real ValR, VecR[NDIM];
  double *G;
  double factorR_C_0 = /*FactorR_C*/1.0/Lev0->MaxN;
  CalculateLapiGcg0(Lat,RT,src);
  for (i=0; i < LocalSizeR[RTAcg]; i++) {
    ValR = LapValR[i]*LapcgR[i] + RSP3(&LapVecR[i*NDIM],&LapiGcgR[i*NDIM]);
    RMat33Vec3(VecR,&LapMatR[i*NDIM_NDIM],&LapiGcgR[i*NDIM]);
    for (j=0;j<NDIM;j++)
      MulVecR[i*NDIM+j] = VecR[j]+LapVecR[i*NDIM+j]*LapcgR[i];
    MulValR[i] = ValR;
  }
  fft_3d_real_to_complex(plan, Lev0->LevelNo, FFTNF1, MulValC, work);
  for (i=0; i < LocalSizeC[RTAcg]; i++) {
    MulValC[i].re *= factorR_C_0;
    MulValC[i].im *= factorR_C_0;
  }
  fft_3d_real_to_complex(plan, Lev0->LevelNo, FFTNF0Vec, MulVecC, work);
  for (i=0; i < LocalSizeC[RTAiGcg]; i++) {
    MulVecC[i].re *= factorR_C_0;
    MulVecC[i].im *= factorR_C_0;
  }
  for (i=0;i<Lev0->MaxG;i++) {
    Index = Lev0->GArray[i].Index;
    G = Lev0->GArray[i].G;
    Lap[i].re = MulValC[Index].re;
    Lap[i].im = MulValC[Index].im;
    for (j=0;j<NDIM;j++) {
      Lap[i].im -= G[j]*(MulVecC[Index*NDIM+j].re);
      Lap[i].re += G[j]*(MulVecC[Index*NDIM+j].im);
    }
  }
}

void CalculateRiemannLaplaceS(const struct Lattice *Lat, struct RiemannTensor *RT, fftw_complex *src, fftw_complex *Lap/*, const double FactorC_R, const double FactorR_C*/) 
{
  int i,i0,iS,j,Index;
  struct  fft_plan_3d *plan = Lat->plan; 
  struct LatticeLevel *LevS = RT->LevS;
  fftw_real **RTAR = RT->DensityR;
  const int *TotalSize = RT->TotalSize;
  fftw_complex *work = RT->TempC;
  fftw_real *LapiGcgR = RTAR[RTAiGcg], *LapcgR = RTAR[RTAcg];
  fftw_real *LapValR = RTAR[RTAARTA];
  fftw_real *LapVecR = RTAR[RTARTA];
  fftw_real *LapMatR = RTAR[RTAIRT];
  fftw_real *MulValR = RTAR[RTAcg];
  fftw_complex *MulValC = (fftw_complex *)MulValR;
  fftw_real *MulVecR = RTAR[RTAiGcg];
  fftw_complex *MulVecC = (fftw_complex *)MulVecR;
  fftw_real ValR, VecR[NDIM];
  int nx,ny,nz;
  const int Nx = LevS->Plan0.plan->local_nx;
  const int Ny = LevS->Plan0.plan->N[1];
  const int Nz = LevS->Plan0.plan->N[2];
  const int NUpx = LevS->NUp[0];
  const int NUpy = LevS->NUp[1];
  const int NUpz = LevS->NUp[2];
  double *G;
  double factorR_C_S = /*FactorR_C*/1.0/LevS->MaxN;
  CalculateLapiGcgS(Lat,RT,src);
  for (nx=0;nx<Nx;nx++)  
    for (ny=0;ny<Ny;ny++) 
      for (nz=0;nz<Nz;nz++) {
	i0 = nz*NUpz+Nz*NUpz*(ny*NUpy+Ny*NUpy*nx*NUpx);
	iS = nz+Nz*(ny+Ny*nx);
	ValR = LapValR[i0]*LapcgR[iS] + RSP3(&LapVecR[i0*NDIM],&LapiGcgR[iS*NDIM]);
	RMat33Vec3(VecR,&LapMatR[i0*NDIM_NDIM],&LapiGcgR[iS*NDIM]);
	for (j=0;j<NDIM;j++) {
	  MulVecR[iS*NDIM+j] = VecR[j]+LapVecR[i0*NDIM+j]*LapcgR[iS];
	}
	MulValR[iS] = ValR;
      }
  fft_3d_real_to_complex(plan, LevS->LevelNo, FFTNF1, MulValC, work);
  for (i=0; i < TotalSize[RTAcg]/LevS->MaxNUp; i++) {
    MulValC[i].re *= factorR_C_S;
    MulValC[i].im *= factorR_C_S;
  }
  fft_3d_real_to_complex(plan, LevS->LevelNo, FFTNFSVec, MulVecC, work);
  for (i=0; i < TotalSize[RTAiGcg]/LevS->MaxNUp; i++) {
    MulVecC[i].re *= factorR_C_S;
    MulVecC[i].im *= factorR_C_S;
  }
  for (i=0;i<LevS->MaxG;i++) {
    Index = LevS->GArray[i].Index;
    G = LevS->GArray[i].G;
    Lap[i].re = MulValC[Index].re;
    Lap[i].im = MulValC[Index].im;
    for (j=0;j<NDIM;j++) {
      Lap[i].im -= G[j]*(MulVecC[Index*NDIM+j].re);
      Lap[i].re += G[j]*(MulVecC[Index*NDIM+j].im);
    }
  }
}

static void CalculateRTData(const struct Lattice *Lat, struct RiemannTensor *RT/*, const double FactorC_R, const double FactorR_C*/) 
{
  struct  fft_plan_3d *plan = Lat->plan; 
  struct LatticeLevel *LevR = RT->LevR;
  struct LatticeLevel *LevS = RT->LevS;
  const double factorR_C_S = /*FactorR_C*/1./LevS->MaxN;
  fftw_complex **RTAC = RT->DensityC;
  fftw_real **RTAR = RT->DensityR;
  fftw_complex **PosFactor = RT->PosFactor;
  const int *TotalSize = RT->TotalSize;
  const int *LocalSizeR = RT->LocalSizeR;
  const int *LocalSizeC = RT->LocalSizeC;
  const int *NFields = RT->NFields;
  const int *MaxNUp = RT->MaxNUp;
  const int MaxNUpS0 = LevS->MaxNUp;
  fftw_complex *PosFactorS0 = LevS->PosFactorUp;
  fftw_complex *coeff, *posfac, *destpos, *destCS, *destCD, *srcC, source;
  fftw_complex *work = RT->TempC;
  fftw_real *workR = (fftw_real*)work;
  fftw_real *srcR, *srcR2, *destR3;
  fftw_real *srcRT = RTAR[RTAIRT], *srcA = workR, *destDet = RTAR[RTADetPreRT], *destRT = RTAR[RTAIRT];
  fftw_real *destIRT = RTAR[RTAIRT], *destA = RTAR[RTAA], *destRTA = RTAR[RTARTA], *destARTA = RTAR[RTAARTA];
  int i,i0,iS,j,Index,pos,n,ng,n0,n1,nx,ny,nz;
  const int Nx = LevS->Plan0.plan->local_nx;
  const int Ny = LevS->Plan0.plan->N[1];
  const int Nz = LevS->Plan0.plan->N[2];
  const int NUpx = LevS->NUp[0];
  const int NUpy = LevS->NUp[1];
  const int NUpz = LevS->NUp[2];
  double *G;
  double MatR[NDIM_NDIM], MatRT[NDIM_NDIM], MatIG[NDIM_NDIM], VecR[NDIM];
  SetArrayToDouble0((double *)RTAC[RTAPreA], TotalSize[RTAPreA]*2);
  SetArrayToDouble0((double *)RTAC[RTAIRT], TotalSize[RTAIRT]*2);
  for (i=0;i < LevR->MaxG;i++) {
    Index = LevR->GArray[i].Index;
    coeff = &RT->Coeff[i];
    G = LevR->GArray[i].G;
    /* RTAIRT */
    posfac = &(PosFactor[RTPFRto0][MaxNUp[RTPFRto0]*i]);
    destpos = &(RTAC[RTAIRT][MaxNUp[RTPFRto0]*Index*NFields[RTAIRT]]);
    for (pos=0; pos < MaxNUp[RTPFRto0]; pos++) {
      for (n1=0; n1 < NDIM; n1++) {
	for (n0=0; n0 < NDIM; n0++) {
	  source.re = -coeff[n0].im*G[n1];
	  source.im = coeff[n0].re*G[n1];
	  n = n0+NDIM*n1;
	  destpos[n+NDIM_NDIM*pos].re = source.re*posfac[pos].re-source.im*posfac[pos].im;
	  destpos[n+NDIM_NDIM*pos].im = source.re*posfac[pos].im+source.im*posfac[pos].re;
	}
      }
    }
    /* RTAPreA */
    posfac = &(PosFactor[RTPFRtoS][MaxNUp[RTPFRtoS]*i]);
    destpos = &(RTAC[RTAPreA][MaxNUp[RTPFRtoS]*Index*NFields[RTAPreA]]);
    for (pos=0; pos < MaxNUp[RTPFRtoS]; pos++) {
      for (n1=0; n1 < NDIM; n1++) {
	for (n0=0; n0 < NDIM; n0++) {
	  for (ng=0; ng < NDIM; ng++) {
	    source.re = -coeff[n0].re*G[n1]*G[ng];
	    source.im = -coeff[n0].im*G[n1]*G[ng];
	    n = ng+NDIM*(n0+NDIM*n1);
	    destpos[n+NDIM_NDIM*NDIM*pos].re = source.re*posfac[pos].re-source.im*posfac[pos].im;
	    destpos[n+NDIM_NDIM*NDIM*pos].im = source.re*posfac[pos].im+source.im*posfac[pos].re;
	  }
	}
      }
    }
  }
  /* RTAIRT */
  for (i=0; i < LevR->MaxDoubleG; i++) {
    destCS = &(RTAC[RTAIRT][LevR->DoubleG[2*i]*MaxNUp[RTPFRto0]*NFields[RTAIRT]]);
    destCD = &(RTAC[RTAIRT][LevR->DoubleG[2*i+1]*MaxNUp[RTPFRto0]*NFields[RTAIRT]]);
    for (pos=0; pos < MaxNUp[RTPFRto0]; pos++) {
      for (n=0; n < NFields[RTAIRT]; n++) {
	destCD[n+NFields[RTAIRT]*pos].re =  destCS[n+NFields[RTAIRT]*pos].re;
	destCD[n+NFields[RTAIRT]*pos].im = -destCS[n+NFields[RTAIRT]*pos].im;
      }
    }
  }
  fft_3d_complex_to_real(plan, LevR->LevelNo, (int)FFTNFRMatUp0, RTAC[RTAIRT], work);
  RiemannRTransformPosNFRto0(RT, RTAR[RTAIRT], workR, NFields[RTAIRT]);
  memcpy(RTAR[RTAIRT],workR,sizeof(fftw_real)*LocalSizeR[RTAIRT]*NFields[RTAIRT]);
  /* RTAPreA */
  for (i=0; i < LevR->MaxDoubleG; i++) {
    destCS = &(RTAC[RTAPreA][LevR->DoubleG[2*i]*MaxNUp[RTPFRtoS]*NFields[RTAPreA]]);
    destCD = &(RTAC[RTAPreA][LevR->DoubleG[2*i+1]*MaxNUp[RTPFRtoS]*NFields[RTAPreA]]);
    for (pos=0; pos < MaxNUp[RTPFRtoS]; pos++) {
      for (n=0; n < NFields[RTAPreA]; n++) {
	destCD[n+NFields[RTAPreA]*pos].re =  destCS[n+NFields[RTAPreA]*pos].re;
	destCD[n+NFields[RTAPreA]*pos].im = -destCS[n+NFields[RTAPreA]*pos].im;
      }
    }
  }
  fft_3d_complex_to_real(plan, LevR->LevelNo, (int)FFTNFRMatVecUpS, RTAC[RTAPreA], work);
  RiemannRTransformPosNFRtoS(RT, RTAR[RTAPreA], workR, NFields[RTAPreA]);
  memcpy(RTAR[RTAPreA],workR,sizeof(fftw_real)*LocalSizeR[RTAPreA]*NFields[RTAPreA]);
  /* RTAA(S) - RTAIRT */
  for (i=0; i< LocalSizeR[RTAIRT];i++) {
    srcR = &(RTAR[RTAIRT][i*NFields[RTAIRT]]);
    srcR[0+NDIM*0] += 1.; /* Diag plus eins */
    srcR[1+NDIM*1] += 1.;
    srcR[2+NDIM*2] += 1.;
  }
  for (nx=0;nx<Nx;nx++)  
    for (ny=0;ny<Ny;ny++) 
      for (nz=0;nz<Nz;nz++) {
	i0 = nz*NUpz+Nz*NUpz*(ny*NUpy+Ny*NUpy*nx*NUpx);
	iS = nz+Nz*(ny+Ny*nx);
	srcR = &(RTAR[RTAIRT][i0*NFields[RTAIRT]]);
	srcR2 = &(RTAR[RTAPreA][iS*NFields[RTAPreA]]);
	destR3 = &(RTAR[RTAA][iS*NFields[RTAA]]);
	for (ng=0; ng<NDIM;ng++) {
	  destR3[ng] 
	    =  srcR2[ng+NDIM*(0+NDIM*0)] * srcR[1+NDIM*1]            * srcR[2+NDIM*2] 
	    +  srcR[0+NDIM*0]            * srcR2[ng+NDIM*(1+NDIM*1)] * srcR[2+NDIM*2]
	    +  srcR[0+NDIM*0]            * srcR[1+NDIM*1]            * srcR2[ng+NDIM*(2+NDIM*2)];
	  destR3[ng] 
	    += srcR2[ng+NDIM*(0+NDIM*1)] * srcR[1+NDIM*2]            * srcR[2+NDIM*0]
	    +  srcR[0+NDIM*1]            * srcR2[ng+NDIM*(1+NDIM*2)] * srcR[2+NDIM*0]
	    +  srcR[0+NDIM*1]            * srcR[1+NDIM*2]            * srcR2[ng+NDIM*(2+NDIM*0)];
	  destR3[ng] 
	    += srcR2[ng+NDIM*(0+NDIM*2)] * srcR[1+NDIM*0]            * srcR[2+NDIM*1]
	    +  srcR[0+NDIM*2]            * srcR2[ng+NDIM*(1+NDIM*0)] * srcR[2+NDIM*1]
	    +  srcR[0+NDIM*2]            * srcR[1+NDIM*0]            * srcR2[ng+NDIM*(2+NDIM*1)];
	  destR3[ng] 
	    +=-srcR2[ng+NDIM*(2+NDIM*0)] * srcR[1+NDIM*1]            * srcR[0+NDIM*2]
	    -  srcR[2+NDIM*0]            * srcR2[ng+NDIM*(1+NDIM*1)] * srcR[0+NDIM*2]
	    -  srcR[2+NDIM*0]            * srcR[1+NDIM*1]            * srcR2[ng+NDIM*(0+NDIM*2)];
	  destR3[ng] 
	    +=-srcR2[ng+NDIM*(2+NDIM*1)] * srcR[1+NDIM*2]            * srcR[0+NDIM*0]
	    -  srcR[2+NDIM*1]            * srcR2[ng+NDIM*(1+NDIM*2)] * srcR[0+NDIM*0]
	    -  srcR[2+NDIM*1]            * srcR[1+NDIM*2]            * srcR2[ng+NDIM*(0+NDIM*0)];
	  destR3[ng] 
	    +=-srcR2[ng+NDIM*(2+NDIM*2)] * srcR[1+NDIM*0]            * srcR[0+NDIM*1]
	    -  srcR[2+NDIM*2]            * srcR2[ng+NDIM*(1+NDIM*0)] * srcR[0+NDIM*1]
	    -  srcR[2+NDIM*2]            * srcR[1+NDIM*0]            * srcR2[ng+NDIM*(0+NDIM*1)];
	}
      }
  fft_3d_real_to_complex(plan, LevS->LevelNo, FFTNFSVec, RTAC[RTAA], work);
  for (i=0; i< LocalSizeC[RTAPreA]*NFields[RTAA]; i++) {    
    work[i].re = RTAC[RTAA][i].re * factorR_C_S;
    work[i].im = RTAC[RTAA][i].im * factorR_C_S;
  }
  /* RTAA (0) */
  SetArrayToDouble0((double *)RTAC[RTAA], TotalSize[RTAA]*2);
  for (i=0;i < LevS->MaxG;i++) {
    Index = LevS->GArray[i].Index;
    posfac = &(PosFactorS0[MaxNUpS0*i]);
    destpos = &(RTAC[RTAA][MaxNUpS0*Index*NFields[RTAA]]);
    srcC = &(work[Index*NFields[RTAA]]);
    for (pos=0; pos < MaxNUpS0; pos++) {
      for (n=0; n < NDIM; n++) {
	destpos[n+NDIM*pos].re = srcC[n].re*posfac[pos].re-srcC[n].im*posfac[pos].im;
	destpos[n+NDIM*pos].im = srcC[n].re*posfac[pos].im+srcC[n].im*posfac[pos].re;
      }
    }
  }
  for (i=0; i < LevS->MaxDoubleG; i++) {
    destCS = &(RTAC[RTAA][LevS->DoubleG[2*i]*MaxNUpS0*NFields[RTAA]]);
    destCD = &(RTAC[RTAA][LevS->DoubleG[2*i+1]*MaxNUpS0*NFields[RTAA]]);
    for (pos=0; pos < MaxNUpS0; pos++) {
      for (n=0; n < NFields[RTAA]; n++) {
	destCD[n+NFields[RTAA]*pos].re =  destCS[n+NFields[RTAA]*pos].re;
	destCD[n+NFields[RTAA]*pos].im = -destCS[n+NFields[RTAA]*pos].im;
      }
    }
  }
  fft_3d_complex_to_real(plan, LevS->LevelNo, (int)FFTNFSVecUp, RTAC[RTAA], work);
  RiemannRTransformPos(LevS, RTAR[RTAA], workR, NFields[RTAA]);
  /* All RTA´s in R*/
  for (i=0; i < LocalSizeR[RTAIRT]; i++) {
    destDet[i] = 1./RDET3(&srcRT[NDIM_NDIM*i]);
    for(j=0;j<NDIM_NDIM;j++) {
      MatR[j] = srcRT[NDIM_NDIM*i+j];
      MatRT[j] = MatR[j];
    }
    RTranspose3(MatRT);
    RMatMat33(&destRT[NDIM_NDIM*i],MatRT,MatR); 
    RMatReci3(MatIG,&destRT[NDIM_NDIM*i]);
    for (j=0;j<NDIM_NDIM;j++)
      destIRT[NDIM_NDIM*i+j] = MatIG[j];
    for(j=0;j<NDIM;j++) {
      VecR[j] = srcA[j]*destDet[i]*0.5;
      destA[NDIM*i+j] = VecR[j];
    }
    RMat33Vec3(&destRTA[i*NDIM],MatIG,VecR);
    destARTA[i]=RSP3(VecR,&destRTA[i*NDIM]);
  }
}
      

void CalculateRiemannTensorData(struct Problem *const P) 
{
  if (P->Lat.RT.ActualUse != active) return;
  CalculateRTData(&P->Lat, &P->Lat.RT/*, P->Lat.FactorDensityR, P->Lat.FactorDensityC*/);
}

void InitRiemannTensor(struct Problem *const P) 
{
  int i, MaxSize=0; 
  struct Lattice *Lat = &P->Lat;
  struct RunStruct *R = &P->R;
  struct LatticeLevel *LevR;
  struct LatticeLevel *LevS = NULL;
  struct LatticeLevel *Lev0 = NULL;
  struct RiemannTensor *RT = &Lat->RT;
  if (P->Lat.RT.Use == UseNotRT) return;
  if (P->Call.out[NormalOut]) fprintf(stderr, "(%i)InitRiemannTensor\n", P->Par.me);
  RT->LevR = LevR = &Lat->Lev[R->LevRNo];
  RT->LevS = LevS = &Lat->Lev[R->LevRSNo];
  RT->Lev0 = Lev0 = &Lat->Lev[R->LevR0No];
  RT->Coeff = (fftw_complex*)
    Malloc(sizeof(fftw_complex)*LevR->MaxG*3,"InitRiemannTensor: RT->Coeff");
  SetArrayToDouble0((double *)RT->Coeff, LevR->MaxG*3*2);
  RT->RTLaplaceS = (fftw_complex*)
    Malloc(sizeof(fftw_complex)*LevS->MaxG,"InitRiemannTensor: RT->Laplace");
  RT->RTLaplace0 = (fftw_complex*)
    Malloc(sizeof(fftw_complex)*Lev0->MaxG,"InitRiemannTensor: RT->Laplace");
  RT->MaxNUp[RTPFRtoS] = LevR->NUp1[0]*LevR->NUp1[1]*LevR->NUp1[2];
  RT->MaxNUp[RTPFRto0] = LevR->NUp0[0]*LevR->NUp0[1]*LevR->NUp0[2];
  /* RTADetPreRT :  */
  RT->TotalSize[RTADetPreRT] = Lev0->Plan0.plan->LocalSizeC;
  RT->LocalSizeC[RTADetPreRT] = Lev0->Plan0.plan->LocalSizeC;
  RT->LocalSizeR[RTADetPreRT] = Lev0->Plan0.plan->LocalSizeR;
  RT->NFields[RTADetPreRT] = 1;
  /* RTAPreA : */
  RT->TotalSize[RTAPreA] = LevR->Plan0.plan->LocalSizeC*RT->MaxNUp[RTPFRtoS]*NDIM_NDIM*NDIM;
  RT->LocalSizeC[RTAPreA] = LevS->Plan0.plan->LocalSizeC;
  RT->LocalSizeR[RTAPreA] = LevS->Plan0.plan->LocalSizeR;
  RT->NFields[RTAPreA] = NDIM_NDIM*NDIM;
  /* RTAA : */
  RT->TotalSize[RTAA] = LevS->Plan0.plan->LocalSizeC*LevS->MaxNUp*NDIM;
  RT->LocalSizeC[RTAA] = Lev0->Plan0.plan->LocalSizeC;
  RT->LocalSizeR[RTAA] = Lev0->Plan0.plan->LocalSizeR;
  RT->NFields[RTAA] = NDIM;
  /* RTAIRT : */
  RT->TotalSize[RTAIRT] = LevR->Plan0.plan->LocalSizeC*RT->MaxNUp[RTPFRto0]*NDIM_NDIM;
  RT->LocalSizeC[RTAIRT] = Lev0->Plan0.plan->LocalSizeC;
  RT->LocalSizeR[RTAIRT] = Lev0->Plan0.plan->LocalSizeR;
  RT->NFields[RTAIRT] = NDIM_NDIM;
  /* RTARTA : */
  RT->TotalSize[RTARTA] = Lev0->Plan0.plan->LocalSizeC*NDIM;
  RT->LocalSizeC[RTARTA] = Lev0->Plan0.plan->LocalSizeC;
  RT->LocalSizeR[RTARTA] = Lev0->Plan0.plan->LocalSizeR;
  RT->NFields[RTARTA] = NDIM;
  /* RTAARTA : */
  RT->TotalSize[RTAARTA] = Lev0->Plan0.plan->LocalSizeC;
  RT->LocalSizeC[RTAARTA] = Lev0->Plan0.plan->LocalSizeC;
  RT->LocalSizeR[RTAARTA] = Lev0->Plan0.plan->LocalSizeR;
  RT->NFields[RTAARTA] = 1;
  /* RTAiGcg : fuer Vektoren*/
  RT->TotalSize[RTAiGcg] = LevS->Plan0.plan->LocalSizeC*LevS->MaxNUp*NDIM;
  RT->LocalSizeC[RTAiGcg] = Lev0->Plan0.plan->LocalSizeC;
  RT->LocalSizeR[RTAiGcg] = Lev0->Plan0.plan->LocalSizeR;
  RT->NFields[RTAiGcg] = NDIM;
  /* RTAcg : fuer values*/
  RT->TotalSize[RTAcg] = LevS->Plan0.plan->LocalSizeC*LevS->MaxNUp;
  RT->LocalSizeC[RTAcg] = Lev0->Plan0.plan->LocalSizeC;
  RT->LocalSizeR[RTAcg] = Lev0->Plan0.plan->LocalSizeR;
  RT->NFields[RTAcg] = 1;
  for (i=0; i< MAXRTARRAYS; i++) {
    switch (i) {
    case RTAiGcg:
      RT->DensityC[i] = (fftw_complex *)
	Malloc(sizeof(fftw_complex)*MAX(RT->TotalSize[i],RT->TotalSize[RTAIRT]/NDIM),"InitRiemannTensor: RT->Density");
      RT->DensityR[i] = (fftw_real *)RT->DensityC[i];
      if ( MAX(RT->TotalSize[i],RT->TotalSize[RTAIRT]/NDIM) > MaxSize) MaxSize = MAX(RT->TotalSize[i],RT->TotalSize[RTAIRT]/NDIM);
      break;
    default:
      RT->DensityC[i] = (fftw_complex *)
	Malloc(sizeof(fftw_complex)*RT->TotalSize[i],"InitRiemannTensor: RT->Density");
      RT->DensityR[i] = (fftw_real *)RT->DensityC[i];
      if (RT->TotalSize[i] > MaxSize) MaxSize = RT->TotalSize[i];
      break;
    }
  }
  RT->TempC = (fftw_complex *)
    Malloc(sizeof(fftw_complex)*MaxSize,"InitRiemannTensor: RT->TempC");
  RT->TempTotalSize = sizeof(fftw_complex)*MaxSize;
}