source: pcp/src/riemann.c@ 1f9f1b

Last change on this file since 1f9f1b was 79290f, checked in by Frederik Heber <heber@…>, 17 years ago

config.h is included in each and every file. After trying to compile on JUMP (with xlc).
  • Property mode set to 100644
File size: 21.4 KB
Line 
1/** \file riemann.c
2 * RiemannTensor-transformations.
3 *
4 * This is all still not really working, thus so far untouched.
5 *
6 Project: ParallelCarParrinello
7 \author Jan Hamaekers
8 \date 2000
9
10 File: riemann.c
11 $Id: riemann.c,v 1.8 2006/03/30 22:19:52 foo Exp $
12*/
13
14#ifdef HAVE_CONFIG_H
15#include <config.h>
16#endif
17
18#include <stdlib.h>
19#include <stdio.h>
20#include <stddef.h>
21#include <math.h>
22#include <string.h>
23
24#include"data.h"
25#include"errors.h"
26#include"helpers.h"
27#include"mymath.h"
28#include"riemann.h"
29#include"myfft.h"
30
31static void RiemannRTransformPosNFRto0(const struct RiemannTensor *RT, fftw_real *source, fftw_real *dest, const int NF)
32{
33 struct LatticeLevel *Lev = RT->LevR;
34 int es = Lev->NUp0[2]*NF;
35 unsigned int cpyes = sizeof(fftw_real)*es;
36 int nx=Lev->Plan0.plan->local_nx,ny=Lev->N[1],nz=Lev->N[2],Nx=Lev->NUp0[0],Ny=Lev->NUp0[1];
37 int lx,ly,z,Lx,Ly;
38 for(lx=0; lx < nx; lx++)
39 for(Lx=0; Lx < Nx; Lx++)
40 for(ly=0; ly < ny; ly++)
41 for(Ly=0; Ly < Ny; Ly++)
42 for(z=0; z < nz; z++) {
43 memcpy( &dest[es*(z+nz*(Ly+Ny*(ly+ny*(Lx+Nx*lx))))],
44 &source[es*(Ly+Ny*(Lx+Nx*(z+nz*(ly+ny*lx))))],
45 cpyes);
46 }
47}
48
49static void RiemannRTransformPosNFRtoS(const struct RiemannTensor *RT, fftw_real *source, fftw_real *dest, const int NF)
50{
51 struct LatticeLevel *Lev = RT->LevR;
52 int es = Lev->NUp1[2]*NF;
53 unsigned int cpyes = sizeof(fftw_real)*es;
54 int nx=Lev->Plan0.plan->local_nx,ny=Lev->N[1],nz=Lev->N[2],Nx=Lev->NUp1[0],Ny=Lev->NUp1[1];
55 int lx,ly,z,Lx,Ly;
56 for(lx=0; lx < nx; lx++)
57 for(Lx=0; Lx < Nx; Lx++)
58 for(ly=0; ly < ny; ly++)
59 for(Ly=0; Ly < Ny; Ly++)
60 for(z=0; z < nz; z++) {
61 memcpy( &dest[es*(z+nz*(Ly+Ny*(ly+ny*(Lx+Nx*lx))))],
62 &source[es*(Ly+Ny*(Lx+Nx*(z+nz*(ly+ny*lx))))],
63 cpyes);
64 }
65}
66
67static void RiemannRTransformPos(const struct LatticeLevel *Lev, fftw_real *source, fftw_real *dest, const int NF)
68{
69 int es = Lev->NUp[2]*NF;
70 unsigned int cpyes = sizeof(fftw_real)*es;
71 int nx=Lev->Plan0.plan->local_nx,ny=Lev->N[1],nz=Lev->N[2],Nx=Lev->NUp[0],Ny=Lev->NUp[1];
72 int lx,ly,z,Lx,Ly;
73 for(lx=0; lx < nx; lx++)
74 for(Lx=0; Lx < Nx; Lx++)
75 for(ly=0; ly < ny; ly++)
76 for(Ly=0; Ly < Ny; Ly++)
77 for(z=0; z < nz; z++) {
78 memcpy( &dest[es*(z+nz*(Ly+Ny*(ly+ny*(Lx+Nx*lx))))],
79 &source[es*(Ly+Ny*(Lx+Nx*(z+nz*(ly+ny*lx))))],
80 cpyes);
81 }
82}
83
84static void CalculateLapiGcgS(const struct Lattice *Lat, const struct RiemannTensor *RT, fftw_complex *source)
85{
86 int i,Index,j;
87 struct fft_plan_3d *plan = Lat->plan;
88 struct LatticeLevel *LevS = RT->LevS;
89 fftw_complex *LapiGcgC = RT->DensityC[RTAiGcg];
90 fftw_complex *LapcgC = RT->DensityC[RTAcg];
91 fftw_complex *work = RT->TempC;
92 fftw_complex *destCS, *destCD;
93 int TotalSize = RT->TotalSize[RTAiGcg]/LevS->MaxNUp;
94 double *G;
95 SetArrayToDouble0((double *)LapiGcgC, TotalSize*2);
96 SetArrayToDouble0((double *)LapcgC, (TotalSize/NDIM)*2);
97 for (i=0;i<LevS->MaxG;i++) {
98 Index = LevS->GArray[i].Index;
99 LapcgC[Index].re = source[i].re;
100 LapcgC[Index].im = source[i].im;
101 G = LevS->GArray[i].G;
102 for (j=0;j<NDIM;j++) {
103 LapiGcgC[Index*NDIM+j].re = -source[i].im*G[j];
104 LapiGcgC[Index*NDIM+j].im = source[i].re*G[j];
105 }
106 }
107 for (i=0; i<LevS->MaxDoubleG; i++) {
108 destCS = &LapiGcgC[LevS->DoubleG[2*i]*NDIM];
109 destCD = &LapiGcgC[LevS->DoubleG[2*i+1]*NDIM];
110 for (j=0; j < NDIM; j++) {
111 destCD[j].re = destCS[j].re;
112 destCD[j].im = -destCS[j].im;
113 }
114 }
115 fft_3d_complex_to_real(plan, LevS->LevelNo, FFTNFSVec, LapiGcgC, work);
116 for (i=0; i<LevS->MaxDoubleG; i++) {
117 destCS = &LapcgC[LevS->DoubleG[2*i]];
118 destCD = &LapcgC[LevS->DoubleG[2*i+1]];
119 destCD[0].re = destCS[0].re;
120 destCD[0].im = -destCS[0].im;
121 }
122 fft_3d_complex_to_real(plan, LevS->LevelNo, FFTNF1, LapcgC, work);
123}
124
125static void CalculateLapiGcg0(const struct Lattice *Lat, const struct RiemannTensor *RT, fftw_complex *source/*, const double FactorC_R, const double FactorR_C*/)
126{
127 int i,Index,j;
128 struct fft_plan_3d *plan = Lat->plan;
129 struct LatticeLevel *Lev0 = RT->Lev0;
130 fftw_complex *LapiGcgC = RT->DensityC[RTAiGcg];
131 fftw_complex *LapcgC = RT->DensityC[RTAcg];
132 fftw_complex *work = RT->TempC;
133 fftw_complex *destCS, *destCD;
134 int TotalSize = RT->TotalSize[RTAiGcg];
135 double *G;
136 SetArrayToDouble0((double *)LapiGcgC, TotalSize*2);
137 SetArrayToDouble0((double *)LapcgC, (TotalSize/NDIM)*2);
138 for (i=0;i<Lev0->MaxG;i++) {
139 Index = Lev0->GArray[i].Index;
140 LapcgC[Index].re = source[i].re;
141 LapcgC[Index].im = source[i].im;
142 G = Lev0->GArray[i].G;
143 for (j=0;j<NDIM;j++) {
144 LapiGcgC[Index*NDIM+j].re = -source[i].im*G[j];
145 LapiGcgC[Index*NDIM+j].im = source[i].re*G[j];
146 }
147 }
148 for (i=0; i<Lev0->MaxDoubleG; i++) {
149 destCS = &LapiGcgC[Lev0->DoubleG[2*i]*NDIM];
150 destCD = &LapiGcgC[Lev0->DoubleG[2*i+1]*NDIM];
151 for (j=0; j < NDIM; j++) {
152 destCD[j].re = destCS[j].re;
153 destCD[j].im = -destCS[j].im;
154 }
155 }
156 fft_3d_complex_to_real(plan, Lev0->LevelNo, FFTNF0Vec, LapiGcgC, work);
157 for (i=0; i<Lev0->MaxDoubleG; i++) {
158 destCS = &LapcgC[Lev0->DoubleG[2*i]];
159 destCD = &LapcgC[Lev0->DoubleG[2*i+1]];
160 destCD[0].re = destCS[0].re;
161 destCD[0].im = -destCS[0].im;
162 }
163 fft_3d_complex_to_real(plan, Lev0->LevelNo, FFTNF1, LapcgC, work);
164}
165
166void CalculateRiemannLaplace0(const struct Lattice *Lat, struct RiemannTensor *RT, fftw_complex *src, fftw_complex *Lap/*, const double FactorC_R, const double FactorR_C*/)
167{
168 int i,j,Index;
169 struct fft_plan_3d *plan = Lat->plan;
170 struct LatticeLevel *Lev0 = RT->Lev0;
171 fftw_real **RTAR = RT->DensityR;
172 const int *LocalSizeR = RT->LocalSizeR;
173 const int *LocalSizeC = RT->LocalSizeC;
174 fftw_complex *work = RT->TempC;
175 fftw_real *LapiGcgR = RTAR[RTAiGcg], *LapcgR = RTAR[RTAcg];
176 fftw_real *LapValR = RTAR[RTAARTA];
177 fftw_real *LapVecR = RTAR[RTARTA];
178 fftw_real *LapMatR = RTAR[RTAIRT];
179 fftw_real *MulValR = RTAR[RTAcg];
180 fftw_complex *MulValC = (fftw_complex *)MulValR;
181 fftw_real *MulVecR = RTAR[RTAiGcg];
182 fftw_complex *MulVecC = (fftw_complex *)MulVecR;
183 fftw_real ValR, VecR[NDIM];
184 double *G;
185 double factorR_C_0 = /*FactorR_C*/1.0/Lev0->MaxN;
186 CalculateLapiGcg0(Lat,RT,src);
187 for (i=0; i < LocalSizeR[RTAcg]; i++) {
188 ValR = LapValR[i]*LapcgR[i] + RSP3(&LapVecR[i*NDIM],&LapiGcgR[i*NDIM]);
189 RMat33Vec3(VecR,&LapMatR[i*NDIM_NDIM],&LapiGcgR[i*NDIM]);
190 for (j=0;j<NDIM;j++)
191 MulVecR[i*NDIM+j] = VecR[j]+LapVecR[i*NDIM+j]*LapcgR[i];
192 MulValR[i] = ValR;
193 }
194 fft_3d_real_to_complex(plan, Lev0->LevelNo, FFTNF1, MulValC, work);
195 for (i=0; i < LocalSizeC[RTAcg]; i++) {
196 MulValC[i].re *= factorR_C_0;
197 MulValC[i].im *= factorR_C_0;
198 }
199 fft_3d_real_to_complex(plan, Lev0->LevelNo, FFTNF0Vec, MulVecC, work);
200 for (i=0; i < LocalSizeC[RTAiGcg]; i++) {
201 MulVecC[i].re *= factorR_C_0;
202 MulVecC[i].im *= factorR_C_0;
203 }
204 for (i=0;i<Lev0->MaxG;i++) {
205 Index = Lev0->GArray[i].Index;
206 G = Lev0->GArray[i].G;
207 Lap[i].re = MulValC[Index].re;
208 Lap[i].im = MulValC[Index].im;
209 for (j=0;j<NDIM;j++) {
210 Lap[i].im -= G[j]*(MulVecC[Index*NDIM+j].re);
211 Lap[i].re += G[j]*(MulVecC[Index*NDIM+j].im);
212 }
213 }
214}
215
216void CalculateRiemannLaplaceS(const struct Lattice *Lat, struct RiemannTensor *RT, fftw_complex *src, fftw_complex *Lap/*, const double FactorC_R, const double FactorR_C*/)
217{
218 int i,i0,iS,j,Index;
219 struct fft_plan_3d *plan = Lat->plan;
220 struct LatticeLevel *LevS = RT->LevS;
221 fftw_real **RTAR = RT->DensityR;
222 const int *TotalSize = RT->TotalSize;
223 fftw_complex *work = RT->TempC;
224 fftw_real *LapiGcgR = RTAR[RTAiGcg], *LapcgR = RTAR[RTAcg];
225 fftw_real *LapValR = RTAR[RTAARTA];
226 fftw_real *LapVecR = RTAR[RTARTA];
227 fftw_real *LapMatR = RTAR[RTAIRT];
228 fftw_real *MulValR = RTAR[RTAcg];
229 fftw_complex *MulValC = (fftw_complex *)MulValR;
230 fftw_real *MulVecR = RTAR[RTAiGcg];
231 fftw_complex *MulVecC = (fftw_complex *)MulVecR;
232 fftw_real ValR, VecR[NDIM];
233 int nx,ny,nz;
234 const int Nx = LevS->Plan0.plan->local_nx;
235 const int Ny = LevS->Plan0.plan->N[1];
236 const int Nz = LevS->Plan0.plan->N[2];
237 const int NUpx = LevS->NUp[0];
238 const int NUpy = LevS->NUp[1];
239 const int NUpz = LevS->NUp[2];
240 double *G;
241 double factorR_C_S = /*FactorR_C*/1.0/LevS->MaxN;
242 CalculateLapiGcgS(Lat,RT,src);
243 for (nx=0;nx<Nx;nx++)
244 for (ny=0;ny<Ny;ny++)
245 for (nz=0;nz<Nz;nz++) {
246 i0 = nz*NUpz+Nz*NUpz*(ny*NUpy+Ny*NUpy*nx*NUpx);
247 iS = nz+Nz*(ny+Ny*nx);
248 ValR = LapValR[i0]*LapcgR[iS] + RSP3(&LapVecR[i0*NDIM],&LapiGcgR[iS*NDIM]);
249 RMat33Vec3(VecR,&LapMatR[i0*NDIM_NDIM],&LapiGcgR[iS*NDIM]);
250 for (j=0;j<NDIM;j++) {
251 MulVecR[iS*NDIM+j] = VecR[j]+LapVecR[i0*NDIM+j]*LapcgR[iS];
252 }
253 MulValR[iS] = ValR;
254 }
255 fft_3d_real_to_complex(plan, LevS->LevelNo, FFTNF1, MulValC, work);
256 for (i=0; i < TotalSize[RTAcg]/LevS->MaxNUp; i++) {
257 MulValC[i].re *= factorR_C_S;
258 MulValC[i].im *= factorR_C_S;
259 }
260 fft_3d_real_to_complex(plan, LevS->LevelNo, FFTNFSVec, MulVecC, work);
261 for (i=0; i < TotalSize[RTAiGcg]/LevS->MaxNUp; i++) {
262 MulVecC[i].re *= factorR_C_S;
263 MulVecC[i].im *= factorR_C_S;
264 }
265 for (i=0;i<LevS->MaxG;i++) {
266 Index = LevS->GArray[i].Index;
267 G = LevS->GArray[i].G;
268 Lap[i].re = MulValC[Index].re;
269 Lap[i].im = MulValC[Index].im;
270 for (j=0;j<NDIM;j++) {
271 Lap[i].im -= G[j]*(MulVecC[Index*NDIM+j].re);
272 Lap[i].re += G[j]*(MulVecC[Index*NDIM+j].im);
273 }
274 }
275}
276
277static void CalculateRTData(const struct Lattice *Lat, struct RiemannTensor *RT/*, const double FactorC_R, const double FactorR_C*/)
278{
279 struct fft_plan_3d *plan = Lat->plan;
280 struct LatticeLevel *LevR = RT->LevR;
281 struct LatticeLevel *LevS = RT->LevS;
282 const double factorR_C_S = /*FactorR_C*/1./LevS->MaxN;
283 fftw_complex **RTAC = RT->DensityC;
284 fftw_real **RTAR = RT->DensityR;
285 fftw_complex **PosFactor = RT->PosFactor;
286 const int *TotalSize = RT->TotalSize;
287 const int *LocalSizeR = RT->LocalSizeR;
288 const int *LocalSizeC = RT->LocalSizeC;
289 const int *NFields = RT->NFields;
290 const int *MaxNUp = RT->MaxNUp;
291 const int MaxNUpS0 = LevS->MaxNUp;
292 fftw_complex *PosFactorS0 = LevS->PosFactorUp;
293 fftw_complex *coeff, *posfac, *destpos, *destCS, *destCD, *srcC, source;
294 fftw_complex *work = RT->TempC;
295 fftw_real *workR = (fftw_real*)work;
296 fftw_real *srcR, *srcR2, *destR3;
297 fftw_real *srcRT = RTAR[RTAIRT], *srcA = workR, *destDet = RTAR[RTADetPreRT], *destRT = RTAR[RTAIRT];
298 fftw_real *destIRT = RTAR[RTAIRT], *destA = RTAR[RTAA], *destRTA = RTAR[RTARTA], *destARTA = RTAR[RTAARTA];
299 int i,i0,iS,j,Index,pos,n,ng,n0,n1,nx,ny,nz;
300 const int Nx = LevS->Plan0.plan->local_nx;
301 const int Ny = LevS->Plan0.plan->N[1];
302 const int Nz = LevS->Plan0.plan->N[2];
303 const int NUpx = LevS->NUp[0];
304 const int NUpy = LevS->NUp[1];
305 const int NUpz = LevS->NUp[2];
306 double *G;
307 double MatR[NDIM_NDIM], MatRT[NDIM_NDIM], MatIG[NDIM_NDIM], VecR[NDIM];
308 SetArrayToDouble0((double *)RTAC[RTAPreA], TotalSize[RTAPreA]*2);
309 SetArrayToDouble0((double *)RTAC[RTAIRT], TotalSize[RTAIRT]*2);
310 for (i=0;i < LevR->MaxG;i++) {
311 Index = LevR->GArray[i].Index;
312 coeff = &RT->Coeff[i];
313 G = LevR->GArray[i].G;
314 /* RTAIRT */
315 posfac = &(PosFactor[RTPFRto0][MaxNUp[RTPFRto0]*i]);
316 destpos = &(RTAC[RTAIRT][MaxNUp[RTPFRto0]*Index*NFields[RTAIRT]]);
317 for (pos=0; pos < MaxNUp[RTPFRto0]; pos++) {
318 for (n1=0; n1 < NDIM; n1++) {
319 for (n0=0; n0 < NDIM; n0++) {
320 source.re = -coeff[n0].im*G[n1];
321 source.im = coeff[n0].re*G[n1];
322 n = n0+NDIM*n1;
323 destpos[n+NDIM_NDIM*pos].re = source.re*posfac[pos].re-source.im*posfac[pos].im;
324 destpos[n+NDIM_NDIM*pos].im = source.re*posfac[pos].im+source.im*posfac[pos].re;
325 }
326 }
327 }
328 /* RTAPreA */
329 posfac = &(PosFactor[RTPFRtoS][MaxNUp[RTPFRtoS]*i]);
330 destpos = &(RTAC[RTAPreA][MaxNUp[RTPFRtoS]*Index*NFields[RTAPreA]]);
331 for (pos=0; pos < MaxNUp[RTPFRtoS]; pos++) {
332 for (n1=0; n1 < NDIM; n1++) {
333 for (n0=0; n0 < NDIM; n0++) {
334 for (ng=0; ng < NDIM; ng++) {
335 source.re = -coeff[n0].re*G[n1]*G[ng];
336 source.im = -coeff[n0].im*G[n1]*G[ng];
337 n = ng+NDIM*(n0+NDIM*n1);
338 destpos[n+NDIM_NDIM*NDIM*pos].re = source.re*posfac[pos].re-source.im*posfac[pos].im;
339 destpos[n+NDIM_NDIM*NDIM*pos].im = source.re*posfac[pos].im+source.im*posfac[pos].re;
340 }
341 }
342 }
343 }
344 }
345 /* RTAIRT */
346 for (i=0; i < LevR->MaxDoubleG; i++) {
347 destCS = &(RTAC[RTAIRT][LevR->DoubleG[2*i]*MaxNUp[RTPFRto0]*NFields[RTAIRT]]);
348 destCD = &(RTAC[RTAIRT][LevR->DoubleG[2*i+1]*MaxNUp[RTPFRto0]*NFields[RTAIRT]]);
349 for (pos=0; pos < MaxNUp[RTPFRto0]; pos++) {
350 for (n=0; n < NFields[RTAIRT]; n++) {
351 destCD[n+NFields[RTAIRT]*pos].re = destCS[n+NFields[RTAIRT]*pos].re;
352 destCD[n+NFields[RTAIRT]*pos].im = -destCS[n+NFields[RTAIRT]*pos].im;
353 }
354 }
355 }
356 fft_3d_complex_to_real(plan, LevR->LevelNo, (int)FFTNFRMatUp0, RTAC[RTAIRT], work);
357 RiemannRTransformPosNFRto0(RT, RTAR[RTAIRT], workR, NFields[RTAIRT]);
358 memcpy(RTAR[RTAIRT],workR,sizeof(fftw_real)*LocalSizeR[RTAIRT]*NFields[RTAIRT]);
359 /* RTAPreA */
360 for (i=0; i < LevR->MaxDoubleG; i++) {
361 destCS = &(RTAC[RTAPreA][LevR->DoubleG[2*i]*MaxNUp[RTPFRtoS]*NFields[RTAPreA]]);
362 destCD = &(RTAC[RTAPreA][LevR->DoubleG[2*i+1]*MaxNUp[RTPFRtoS]*NFields[RTAPreA]]);
363 for (pos=0; pos < MaxNUp[RTPFRtoS]; pos++) {
364 for (n=0; n < NFields[RTAPreA]; n++) {
365 destCD[n+NFields[RTAPreA]*pos].re = destCS[n+NFields[RTAPreA]*pos].re;
366 destCD[n+NFields[RTAPreA]*pos].im = -destCS[n+NFields[RTAPreA]*pos].im;
367 }
368 }
369 }
370 fft_3d_complex_to_real(plan, LevR->LevelNo, (int)FFTNFRMatVecUpS, RTAC[RTAPreA], work);
371 RiemannRTransformPosNFRtoS(RT, RTAR[RTAPreA], workR, NFields[RTAPreA]);
372 memcpy(RTAR[RTAPreA],workR,sizeof(fftw_real)*LocalSizeR[RTAPreA]*NFields[RTAPreA]);
373 /* RTAA(S) - RTAIRT */
374 for (i=0; i< LocalSizeR[RTAIRT];i++) {
375 srcR = &(RTAR[RTAIRT][i*NFields[RTAIRT]]);
376 srcR[0+NDIM*0] += 1.; /* Diag plus eins */
377 srcR[1+NDIM*1] += 1.;
378 srcR[2+NDIM*2] += 1.;
379 }
380 for (nx=0;nx<Nx;nx++)
381 for (ny=0;ny<Ny;ny++)
382 for (nz=0;nz<Nz;nz++) {
383 i0 = nz*NUpz+Nz*NUpz*(ny*NUpy+Ny*NUpy*nx*NUpx);
384 iS = nz+Nz*(ny+Ny*nx);
385 srcR = &(RTAR[RTAIRT][i0*NFields[RTAIRT]]);
386 srcR2 = &(RTAR[RTAPreA][iS*NFields[RTAPreA]]);
387 destR3 = &(RTAR[RTAA][iS*NFields[RTAA]]);
388 for (ng=0; ng<NDIM;ng++) {
389 destR3[ng]
390 = srcR2[ng+NDIM*(0+NDIM*0)] * srcR[1+NDIM*1] * srcR[2+NDIM*2]
391 + srcR[0+NDIM*0] * srcR2[ng+NDIM*(1+NDIM*1)] * srcR[2+NDIM*2]
392 + srcR[0+NDIM*0] * srcR[1+NDIM*1] * srcR2[ng+NDIM*(2+NDIM*2)];
393 destR3[ng]
394 += srcR2[ng+NDIM*(0+NDIM*1)] * srcR[1+NDIM*2] * srcR[2+NDIM*0]
395 + srcR[0+NDIM*1] * srcR2[ng+NDIM*(1+NDIM*2)] * srcR[2+NDIM*0]
396 + srcR[0+NDIM*1] * srcR[1+NDIM*2] * srcR2[ng+NDIM*(2+NDIM*0)];
397 destR3[ng]
398 += srcR2[ng+NDIM*(0+NDIM*2)] * srcR[1+NDIM*0] * srcR[2+NDIM*1]
399 + srcR[0+NDIM*2] * srcR2[ng+NDIM*(1+NDIM*0)] * srcR[2+NDIM*1]
400 + srcR[0+NDIM*2] * srcR[1+NDIM*0] * srcR2[ng+NDIM*(2+NDIM*1)];
401 destR3[ng]
402 +=-srcR2[ng+NDIM*(2+NDIM*0)] * srcR[1+NDIM*1] * srcR[0+NDIM*2]
403 - srcR[2+NDIM*0] * srcR2[ng+NDIM*(1+NDIM*1)] * srcR[0+NDIM*2]
404 - srcR[2+NDIM*0] * srcR[1+NDIM*1] * srcR2[ng+NDIM*(0+NDIM*2)];
405 destR3[ng]
406 +=-srcR2[ng+NDIM*(2+NDIM*1)] * srcR[1+NDIM*2] * srcR[0+NDIM*0]
407 - srcR[2+NDIM*1] * srcR2[ng+NDIM*(1+NDIM*2)] * srcR[0+NDIM*0]
408 - srcR[2+NDIM*1] * srcR[1+NDIM*2] * srcR2[ng+NDIM*(0+NDIM*0)];
409 destR3[ng]
410 +=-srcR2[ng+NDIM*(2+NDIM*2)] * srcR[1+NDIM*0] * srcR[0+NDIM*1]
411 - srcR[2+NDIM*2] * srcR2[ng+NDIM*(1+NDIM*0)] * srcR[0+NDIM*1]
412 - srcR[2+NDIM*2] * srcR[1+NDIM*0] * srcR2[ng+NDIM*(0+NDIM*1)];
413 }
414 }
415 fft_3d_real_to_complex(plan, LevS->LevelNo, FFTNFSVec, RTAC[RTAA], work);
416 for (i=0; i< LocalSizeC[RTAPreA]*NFields[RTAA]; i++) {
417 work[i].re = RTAC[RTAA][i].re * factorR_C_S;
418 work[i].im = RTAC[RTAA][i].im * factorR_C_S;
419 }
420 /* RTAA (0) */
421 SetArrayToDouble0((double *)RTAC[RTAA], TotalSize[RTAA]*2);
422 for (i=0;i < LevS->MaxG;i++) {
423 Index = LevS->GArray[i].Index;
424 posfac = &(PosFactorS0[MaxNUpS0*i]);
425 destpos = &(RTAC[RTAA][MaxNUpS0*Index*NFields[RTAA]]);
426 srcC = &(work[Index*NFields[RTAA]]);
427 for (pos=0; pos < MaxNUpS0; pos++) {
428 for (n=0; n < NDIM; n++) {
429 destpos[n+NDIM*pos].re = srcC[n].re*posfac[pos].re-srcC[n].im*posfac[pos].im;
430 destpos[n+NDIM*pos].im = srcC[n].re*posfac[pos].im+srcC[n].im*posfac[pos].re;
431 }
432 }
433 }
434 for (i=0; i < LevS->MaxDoubleG; i++) {
435 destCS = &(RTAC[RTAA][LevS->DoubleG[2*i]*MaxNUpS0*NFields[RTAA]]);
436 destCD = &(RTAC[RTAA][LevS->DoubleG[2*i+1]*MaxNUpS0*NFields[RTAA]]);
437 for (pos=0; pos < MaxNUpS0; pos++) {
438 for (n=0; n < NFields[RTAA]; n++) {
439 destCD[n+NFields[RTAA]*pos].re = destCS[n+NFields[RTAA]*pos].re;
440 destCD[n+NFields[RTAA]*pos].im = -destCS[n+NFields[RTAA]*pos].im;
441 }
442 }
443 }
444 fft_3d_complex_to_real(plan, LevS->LevelNo, (int)FFTNFSVecUp, RTAC[RTAA], work);
445 RiemannRTransformPos(LevS, RTAR[RTAA], workR, NFields[RTAA]);
446 /* All RTAï¿œs in R*/
447 for (i=0; i < LocalSizeR[RTAIRT]; i++) {
448 destDet[i] = 1./RDET3(&srcRT[NDIM_NDIM*i]);
449 for(j=0;j<NDIM_NDIM;j++) {
450 MatR[j] = srcRT[NDIM_NDIM*i+j];
451 MatRT[j] = MatR[j];
452 }
453 RTranspose3(MatRT);
454 RMatMat33(&destRT[NDIM_NDIM*i],MatRT,MatR);
455 RMatReci3(MatIG,&destRT[NDIM_NDIM*i]);
456 for (j=0;j<NDIM_NDIM;j++)
457 destIRT[NDIM_NDIM*i+j] = MatIG[j];
458 for(j=0;j<NDIM;j++) {
459 VecR[j] = srcA[j]*destDet[i]*0.5;
460 destA[NDIM*i+j] = VecR[j];
461 }
462 RMat33Vec3(&destRTA[i*NDIM],MatIG,VecR);
463 destARTA[i]=RSP3(VecR,&destRTA[i*NDIM]);
464 }
465}
466
467
468void CalculateRiemannTensorData(struct Problem *const P)
469{
470 if (P->Lat.RT.ActualUse != active) return;
471 CalculateRTData(&P->Lat, &P->Lat.RT/*, P->Lat.FactorDensityR, P->Lat.FactorDensityC*/);
472}
473
474void InitRiemannTensor(struct Problem *const P)
475{
476 int i, MaxSize=0;
477 struct Lattice *Lat = &P->Lat;
478 struct RunStruct *R = &P->R;
479 struct LatticeLevel *LevR;
480 struct LatticeLevel *LevS = NULL;
481 struct LatticeLevel *Lev0 = NULL;
482 struct RiemannTensor *RT = &Lat->RT;
483 if (P->Lat.RT.Use == UseNotRT) return;
484 if (P->Call.out[NormalOut]) fprintf(stderr, "(%i)InitRiemannTensor\n", P->Par.me);
485 RT->LevR = LevR = &Lat->Lev[R->LevRNo];
486 RT->LevS = LevS = &Lat->Lev[R->LevRSNo];
487 RT->Lev0 = Lev0 = &Lat->Lev[R->LevR0No];
488 RT->Coeff = (fftw_complex*)
489 Malloc(sizeof(fftw_complex)*LevR->MaxG*3,"InitRiemannTensor: RT->Coeff");
490 SetArrayToDouble0((double *)RT->Coeff, LevR->MaxG*3*2);
491 RT->RTLaplaceS = (fftw_complex*)
492 Malloc(sizeof(fftw_complex)*LevS->MaxG,"InitRiemannTensor: RT->Laplace");
493 RT->RTLaplace0 = (fftw_complex*)
494 Malloc(sizeof(fftw_complex)*Lev0->MaxG,"InitRiemannTensor: RT->Laplace");
495 RT->MaxNUp[RTPFRtoS] = LevR->NUp1[0]*LevR->NUp1[1]*LevR->NUp1[2];
496 RT->MaxNUp[RTPFRto0] = LevR->NUp0[0]*LevR->NUp0[1]*LevR->NUp0[2];
497 /* RTADetPreRT : */
498 RT->TotalSize[RTADetPreRT] = Lev0->Plan0.plan->LocalSizeC;
499 RT->LocalSizeC[RTADetPreRT] = Lev0->Plan0.plan->LocalSizeC;
500 RT->LocalSizeR[RTADetPreRT] = Lev0->Plan0.plan->LocalSizeR;
501 RT->NFields[RTADetPreRT] = 1;
502 /* RTAPreA : */
503 RT->TotalSize[RTAPreA] = LevR->Plan0.plan->LocalSizeC*RT->MaxNUp[RTPFRtoS]*NDIM_NDIM*NDIM;
504 RT->LocalSizeC[RTAPreA] = LevS->Plan0.plan->LocalSizeC;
505 RT->LocalSizeR[RTAPreA] = LevS->Plan0.plan->LocalSizeR;
506 RT->NFields[RTAPreA] = NDIM_NDIM*NDIM;
507 /* RTAA : */
508 RT->TotalSize[RTAA] = LevS->Plan0.plan->LocalSizeC*LevS->MaxNUp*NDIM;
509 RT->LocalSizeC[RTAA] = Lev0->Plan0.plan->LocalSizeC;
510 RT->LocalSizeR[RTAA] = Lev0->Plan0.plan->LocalSizeR;
511 RT->NFields[RTAA] = NDIM;
512 /* RTAIRT : */
513 RT->TotalSize[RTAIRT] = LevR->Plan0.plan->LocalSizeC*RT->MaxNUp[RTPFRto0]*NDIM_NDIM;
514 RT->LocalSizeC[RTAIRT] = Lev0->Plan0.plan->LocalSizeC;
515 RT->LocalSizeR[RTAIRT] = Lev0->Plan0.plan->LocalSizeR;
516 RT->NFields[RTAIRT] = NDIM_NDIM;
517 /* RTARTA : */
518 RT->TotalSize[RTARTA] = Lev0->Plan0.plan->LocalSizeC*NDIM;
519 RT->LocalSizeC[RTARTA] = Lev0->Plan0.plan->LocalSizeC;
520 RT->LocalSizeR[RTARTA] = Lev0->Plan0.plan->LocalSizeR;
521 RT->NFields[RTARTA] = NDIM;
522 /* RTAARTA : */
523 RT->TotalSize[RTAARTA] = Lev0->Plan0.plan->LocalSizeC;
524 RT->LocalSizeC[RTAARTA] = Lev0->Plan0.plan->LocalSizeC;
525 RT->LocalSizeR[RTAARTA] = Lev0->Plan0.plan->LocalSizeR;
526 RT->NFields[RTAARTA] = 1;
527 /* RTAiGcg : fuer Vektoren*/
528 RT->TotalSize[RTAiGcg] = LevS->Plan0.plan->LocalSizeC*LevS->MaxNUp*NDIM;
529 RT->LocalSizeC[RTAiGcg] = Lev0->Plan0.plan->LocalSizeC;
530 RT->LocalSizeR[RTAiGcg] = Lev0->Plan0.plan->LocalSizeR;
531 RT->NFields[RTAiGcg] = NDIM;
532 /* RTAcg : fuer values*/
533 RT->TotalSize[RTAcg] = LevS->Plan0.plan->LocalSizeC*LevS->MaxNUp;
534 RT->LocalSizeC[RTAcg] = Lev0->Plan0.plan->LocalSizeC;
535 RT->LocalSizeR[RTAcg] = Lev0->Plan0.plan->LocalSizeR;
536 RT->NFields[RTAcg] = 1;
537 for (i=0; i< MAXRTARRAYS; i++) {
538 switch (i) {
539 case RTAiGcg:
540 RT->DensityC[i] = (fftw_complex *)
541 Malloc(sizeof(fftw_complex)*MAX(RT->TotalSize[i],RT->TotalSize[RTAIRT]/NDIM),"InitRiemannTensor: RT->Density");
542 RT->DensityR[i] = (fftw_real *)RT->DensityC[i];
543 if ( MAX(RT->TotalSize[i],RT->TotalSize[RTAIRT]/NDIM) > MaxSize) MaxSize = MAX(RT->TotalSize[i],RT->TotalSize[RTAIRT]/NDIM);
544 break;
545 default:
546 RT->DensityC[i] = (fftw_complex *)
547 Malloc(sizeof(fftw_complex)*RT->TotalSize[i],"InitRiemannTensor: RT->Density");
548 RT->DensityR[i] = (fftw_real *)RT->DensityC[i];
549 if (RT->TotalSize[i] > MaxSize) MaxSize = RT->TotalSize[i];
550 break;
551 }
552 }
553 RT->TempC = (fftw_complex *)
554 Malloc(sizeof(fftw_complex)*MaxSize,"InitRiemannTensor: RT->TempC");
555 RT->TempTotalSize = sizeof(fftw_complex)*MaxSize;
556}
557
Note: See TracBrowser for help on using the repository browser.