source: pcp/src/riemann.c@ e936b3

non_adaptive_molecuilder
Last change on this file since e936b3 was a0bcf1, checked in by Frederik Heber <heber@…>, 17 years ago

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