📄 fst_semi_memo.cpp
字号:
}
}
void FST_semi_memo::Conv2Sphere_semi_memo(double *rdata, double *idata, double *rfilter, double *ifilter, double *rres, double *ires, int size, double *workspace)
{
int bw, spharmonic_bound;
double *frres, *fires, *filtrres, *filtires, *trres, *tires;
double **spharmonic_pml_table, **transpose_spharmonic_pml_table;
double *spharmonic_result_space, *transpose_spharmonic_result_space;
double *scratchpad;
bw = size/2;
/* assign space */
spharmonic_bound =cosp.Spharmonic_TableSize(bw);
spharmonic_result_space = workspace;
transpose_spharmonic_result_space = spharmonic_result_space +
spharmonic_bound;
frres = transpose_spharmonic_result_space +
spharmonic_bound; /* needs (bw*bw) */
fires = frres + (bw*bw); /* needs (bw*bw) */
trres = fires + (bw*bw); /* needs (bw*bw) */
tires = trres + (bw*bw); /* needs (bw*bw) */
filtrres = tires + (bw*bw); /* needs bw */
filtires = filtrres + bw; /* needs bw */
scratchpad = filtires + bw; /* needs (8 * bw^2) + (32 * bw) */
spharmonic_pml_table =
cosp.Spharmonic_Pml_Table(bw,
spharmonic_result_space,
scratchpad);
transpose_spharmonic_pml_table =
cosp.Transpose_Spharmonic_Pml_Table(spharmonic_pml_table,
bw,
transpose_spharmonic_result_space,
scratchpad);
FST_semi_memo1(rdata, idata,
frres, fires,
size,
spharmonic_pml_table,
scratchpad,
1,
bw);
FZT_semi_memo(rfilter, ifilter,
filtrres, filtires,
size,
spharmonic_pml_table[0],
scratchpad,
1);
TransMult(frres, fires, filtrres, filtires, trres, tires, size);
InvFST_semi_memo(trres, tires,
rres, ires,
size,
transpose_spharmonic_pml_table,
scratchpad,
1,
bw);
/***
have to free the memory that was allocated in
Spharmonic_Pml_Table() and
Transpose_Spharmonic_Pml_Table()
***/
free(spharmonic_pml_table);
free(transpose_spharmonic_pml_table);
}
void FST_semi_memo::FST_semi_memo2(double *rdata, double *rcoeffs, double *icoeffs, int size, double **seminaive_naive_table, double *workspace, double *wSave, double *CoswSave, int cutoff)
{
int bw, m, i, j;
double *rres, *ires;
double *rdataptr, *idataptr;
double *scratchpad;
int tmpindex[2];
double pow_one;
double *cos_even;
double *tmprdata_ptr;
int l, dummy ;
double tmpA, tmpB ;
bw = size/2;
/* assign space */
rres = workspace; /* needs (size * size) = (4 * bw^2) */
cos_even = rres + (size * size); /* needs bw */
scratchpad = cos_even + (2*bw); /* needs (24 * bw) ?? */
/* total workspace is (4 * bw^2) + (29 * bw) */
/* first save input data: don't want to write over it! */
memcpy(rres, rdata, sizeof(double) * size * size);
tmprdata_ptr = rres; /* to be used for advancing */
/* do the FFTs along phi */
fftg1.grid_fourier_FFTPACK(rres, size, wSave);
/* point to start of output data buffers */
rdataptr = rcoeffs;
idataptr = icoeffs;
for (m=0; m<bw; m++) {
/*
fprintf(stderr,"m = %d\n",m);
*/
/*** test to see if before cutoff or after ***/
if (m < cutoff){
/* do the real part */
semi1.SemiNaiveReduced(tmprdata_ptr,
bw,
m,
rdataptr,
seminaive_naive_table[m],
scratchpad,
cos_even,
CoswSave );
rdataptr += bw-m;
tmprdata_ptr += size;
/* do ONLY IF m != 0; otherwise set imaginary coeffs = 0 */
if (m != 0)
{
/* do imaginary part */
semi1.SemiNaiveReduced(tmprdata_ptr,
bw,
m,
idataptr,
seminaive_naive_table[m],
scratchpad,
cos_even,
CoswSave );
idataptr += bw-m;
tmprdata_ptr += size;
}
else
{
memset(idataptr, 0, sizeof(double) * bw);
idataptr += bw;
}
}
else{
/* do real part */
naive1.Naive_AnalysisX(tmprdata_ptr,
bw,
m,
rdataptr,
seminaive_naive_table[m],
workspace);
rdataptr += bw-m;
tmprdata_ptr += size;
/* do ONLY IF m != 0; otherwise set imaginary coeffs = 0 */
if (m != 0)
{
/* do imaginary part */
naive1.Naive_AnalysisX(tmprdata_ptr,
bw,
m,
idataptr,
seminaive_naive_table[m],
workspace);
idataptr += bw-m;
tmprdata_ptr += size;
}
else
{
memset(idataptr, 0, sizeof(double) * bw);
idataptr += bw;
}
}
}
/*** now do upper coefficients ****/
pow_one = 1.0;
for(i = 1; i < bw; i++)
{
pow_one *= -1.0;
for( j = i; j < bw; j++)
{
prif.seanindex2(i, j, bw, tmpindex);
rcoeffs[tmpindex[1]] = pow_one * rcoeffs[tmpindex[0]];
icoeffs[tmpindex[1]] = -1.0 * pow_one * icoeffs[tmpindex[0]];
}
}
/*****************
New in 2.6: Need to massage the coefficients one more time,
to get that right factor in there (how embarrassing)
******************/
tmpA = 2. * sqrt( PI );
tmpB = sqrt( 2. * PI );
for(m=0;m<bw;m++)
for(l=m;l<bw;l++)
{
dummy = prif.seanindex(m,l,bw);
if ( m == 0 )
{
rcoeffs[dummy] *= tmpA ;
icoeffs[dummy] *= tmpA ;
}
else
{
rcoeffs[dummy] *= tmpB ;
icoeffs[dummy] *= tmpB ;
}
/* now for the negative-order coefficients */
if ( m != 0 )
{
dummy = prif.seanindex(-m,l,bw);
rcoeffs[dummy] *= tmpB ;
icoeffs[dummy] *= tmpB ;
}
}
}
void FST_semi_memo::InvFST_semi_memo2(double *rcoeffs, double *icoeffs, double *rdata, int size, double **transpose_seminaive_naive_table, double *workspace, double *wSave, int cutoff)
{
int bw, m, i, n;
double *rdataptr, *idataptr;
double *scratchpad;
double *sin_values, *eval_pts;
double *tmp_ptr;
int l, dummy ;
double tmpA, tmpB ;
FILE *fp;
bw = size/2;
/* allocate space */
sin_values = workspace; /* needs (2 * bw) */
eval_pts = sin_values + (2 * bw); /* needs (2 * bw) */
scratchpad = eval_pts + (2 * bw); /* needs (16 * bw) */
/* load up the sin_values array */
n = 2*bw;
prim1.ArcCosEvalPts(n, eval_pts);
for (i=0; i<n; i++)
sin_values[i] = sin(eval_pts[i]);
/*****************
New in 2.6: Need to massage the coefficients one more time,
to get that right factor in there (how embarrassing)
******************/
tmpA = 1./(2. * sqrt(PI) );
tmpB = 1./sqrt(2. * PI ) ;
for(m=0;m<bw;m++)
for(l=m;l<bw;l++)
{
dummy = prif.seanindex(m,l,bw);
if ( m == 0 )
{
rcoeffs[dummy] *= tmpA ;
icoeffs[dummy] *= tmpA ;
}
else
{
rcoeffs[dummy] *= tmpB ;
icoeffs[dummy] *= tmpB ;
}
/* now for the negative-order coefficients */
if ( m != 0 )
{
dummy = prif.seanindex(-m,l,bw);
rcoeffs[dummy] *= tmpB ;
icoeffs[dummy] *= tmpB ;
}
}
/* Now do all of the inverse Legendre transforms */
tmp_ptr = rdata;
rdataptr = rcoeffs;
idataptr = icoeffs;
for (m=0; m<bw; m++) {
/*
fprintf(stderr,"m = %d\n",m);
*/
if(m < cutoff){
/* do real part first */
semi1.InvSemiNaiveReduced(rdataptr,
bw,
m,
tmp_ptr,
transpose_seminaive_naive_table[m],
sin_values,
scratchpad);
tmp_ptr += size;
/* move to next set of coeffs */
rdataptr += bw - m;
if ( m != 0)
{
/* now do imaginary part */
semi1.InvSemiNaiveReduced(idataptr,
bw,
m,
tmp_ptr,
transpose_seminaive_naive_table[m],
sin_values,
scratchpad);
tmp_ptr += size;
/* move to next set of coefs */
idataptr += bw - m;
}
else
{
/* move to next set of coeffs */
idataptr += bw - m;
}
}
else
{
/* first do the real part */
naive1.Naive_SynthesizeX(rdataptr,
bw,
m,
tmp_ptr,
transpose_seminaive_naive_table[m]);
tmp_ptr += size;
/* move to next set of coeffs */
rdataptr += bw - m;
if ( m!= 0 )
{
/* now do the imaginary */
naive1.Naive_SynthesizeX(idataptr,
bw,
m,
tmp_ptr,
transpose_seminaive_naive_table[m]);
tmp_ptr += size;
/* move to next set of coefs */
idataptr += bw - m;
}
else
{
/* move to next set of coeffs */
idataptr += bw - m;
}
}
} /* closes m loop */
/* now fill in zero values where m = bw (from problem definition) */
memset(tmp_ptr, 0, sizeof(double) * size);
/** now transpose **/
prif.transpose(rdata, size);
/* now do inverse fourier grid computation */
fftg1.grid_invfourier_FFTPACK(rdata, size, wSave);
/* amscray */
}
void FST_semi_memo::FZT_semi_memo2(double *rdata, double *rres, double *ires, int size, double *cos_pml_table, double *workspace, int dataformat, double *wSave)
{
int i, j, bw;
double *r0, *i0, dsize;
double tmpreal, tmpimag;
double *scratchpad;
double *cos_even;
int l, dummy ;
double tmpA ;
bw = size/2;
/* assign memory */
r0 = workspace; /* needs (2 * bw) */
i0 = r0 + (2 * bw); /* needs (2 * bw) */
cos_even = i0 + (2 * bw); /* needs (1 * bw) */
scratchpad = cos_even + (1 * bw); /* needs (8 * bw) */
/* total workspace = 13*bw */
dsize = (double) size;
/* compute the m = 0 components */
for (i=0; i<size; i++) {
tmpreal = 0.0;
for(j=0; j<size; j++) {
tmpreal += rdata[(i*size)+j];
}
/* normalize */
r0[i] = tmpreal/dsize;
}
/* do the real part */
semi1.SemiNaiveReduced(r0,
bw,
0,
rres,
cos_pml_table,
scratchpad,
cos_even,
wSave);
/* set imag coefficients = 0 */
memset(ires, 0, sizeof(double) * size);
/*****************
New in 2.6: Need to massage the coefficients one more time,
to get that right factor in there (how embarrassing)
******************/
tmpA = 2. * sqrt( PI );
for(l=0;l<bw;l++)
{
dummy = prif.seanindex(0,l,bw);
rres[dummy] *= tmpA ;
}
}
void FST_semi_memo::TransMult2(double *rdatacoeffs, double *idatacoeffs, double *rfiltercoeffs, double *ifiltercoeffs, double *rres, double *ires, int size)
{
int m, l, bw;
double *rdptr, *idptr, *rrptr, *irptr;
bw = size/2;
rdptr = rdatacoeffs;
idptr = idatacoeffs;
rrptr = rres;
irptr = ires;
for (m=0; m<bw; m++) {
for (l=m; l<bw; l++) {
compmult(rfiltercoeffs[l], ifiltercoeffs[l],
rdptr[l-m], idptr[l-m],
rrptr[l-m], irptr[l-m]);
rrptr[l-m] *= sqrt(4*M_PI/(2*l+1));
irptr[l-m] *= sqrt(4*M_PI/(2*l+1));
}
rdptr += bw-m; idptr += bw-m;
rrptr += bw-m; irptr += bw-m;
}
for (m=bw+1; m<size; m++) {
for (l=size-m; l<bw; l++){
compmult(rfiltercoeffs[l], ifiltercoeffs[l],
rdptr[l-size+m], idptr[l-size+m],
rrptr[l-size+m], irptr[l-size+m]);
rrptr[l-size+m] *= sqrt(4*M_PI/(2*l+1));
irptr[l-size+m] *= sqrt(4*M_PI/(2*l+1));
}
rdptr += m-bw; idptr += m-bw;
rrptr += m-bw; irptr += m-bw;
}
}
void FST_semi_memo::Conv2Sphere_semi_memo2(double *rdata, double *rfilter, double *rres, int size, double *workspace, double *wSave, double *CoswSave)
{
int bw, spharmonic_bound;
double *frres, *fires, *filtrres, *filtires, *trres, *tires;
double **spharmonic_pml_table, **transpose_spharmonic_pml_table;
double *spharmonic_result_space, *transpose_spharmonic_result_space;
double *scratchpad;
bw = size/2;
/* assign space */
spharmonic_bound = cosp.Spharmonic_TableSize(bw);
spharmonic_result_space = workspace;
transpose_spharmonic_result_space = spharmonic_result_space +
spharmonic_bound;
frres = transpose_spharmonic_result_space +
spharmonic_bound; /* needs (bw*bw) */
fires = frres + (bw*bw); /* needs (bw*bw) */
trres = fires + (bw*bw); /* needs (bw*bw) */
tires = trres + (bw*bw); /* needs (bw*bw) */
filtrres = tires + (bw*bw); /* needs bw */
filtires = filtrres + bw; /* needs bw */
scratchpad = filtires + bw; /* needs (8 * bw^2) + (32 * bw) */
spharmonic_pml_table =
cosp.Spharmonic_Pml_Table(bw,
spharmonic_result_space,
scratchpad);
transpose_spharmonic_pml_table =
cosp.Transpose_Spharmonic_Pml_Table(spharmonic_pml_table,
bw,
transpose_spharmonic_result_space,
scratchpad);
FST_semi_memo2(rdata,
frres, fires,
size,
spharmonic_pml_table,
scratchpad,
wSave,
CoswSave,
bw);
/*
FZT_semi_memo2(rfilter,
filtrres, filtires,
size,
spharmonic_pml_table[0],
scratchpad,
1,
CoswSave );
TransMult2(frres, fires, filtrres, filtires, trres, tires, size);
InvFST_semi_memo2(trres, tires,
rres,
size,
transpose_spharmonic_pml_table,
scratchpad,
wSave,
bw);*/
/***
have to free the memory that was allocated in
Spharmonic_Pml_Table() and
Transpose_Spharmonic_Pml_Table()
***/
free(spharmonic_pml_table);
free(transpose_spharmonic_pml_table);
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -