本文整理汇总了C++中creal函数的典型用法代码示例。如果您正苦于以下问题：C++ creal函数的具体用法？C++ creal怎么用？C++ creal使用的例子？那么恭喜您, 这里精选的函数代码示例或许可以为您提供帮助。
在下文中一共展示了creal函数的15个代码示例，这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞，您的评价将有助于系统推荐出更棒的C++代码示例。
示例1: test_13
JL_DLLEXPORT struct13 test_13(struct13 a, double b) {
//Unpack a nested ComplexPair{Float64} struct
if (verbose) fprintf(stderr,"%g + %g i & %g\n", creal(a.x), cimag(a.x), b);
a.x += b*1 - (b*2.0*I);
return a;
}示例2: TEST
TEST(complex, creal) {
ASSERT_EQ(0.0, creal(0));
}示例3: cimag
VrArrayPtrCF32 BlasComplexSingle::scal_minus(VrArrayPtrCF32 A, float complex scal) {
//scal[0]=-scal[0];
//scal[1]=-scal[1];
scal = -creal(scal) - cimag(scal)*I;
return scal_add(VR_GET_NDIMS_CF32(A), A,scal);
}示例4: testBSplineHAtom
int testBSplineHAtom() {
PrintTimeStamp(PETSC_COMM_SELF, "H atom", NULL);
MPI_Comm comm = PETSC_COMM_SELF;
BPS bps; BPSCreate(comm, &bps); BPSSetExp(bps, 30.0, 61, 3.0);
int order = 5;
BSS bss; BSSCreate(comm, &bss); BSSSetKnots(bss, order, bps); BSSSetUp(bss);
Mat H; BSSCreateR1Mat(bss, &H);
Mat S; BSSCreateR1Mat(bss, &S);
Mat V; BSSCreateR1Mat(bss, &V);
BSSD2R1Mat(bss, H);
MatScale(H, -0.5);
BSSENR1Mat(bss, 0, 0.0, V);
MatAXPY(H, -1.0, V, DIFFERENT_NONZERO_PATTERN);
BSSSR1Mat(bss, S);
// -- initial space --
Pot psi0; PotCreate(comm, &psi0); PotSetSlater(psi0, 2.0, 1, 1.1);
int n_init_space = 1;
Vec *xs; PetscMalloc1(n_init_space, &xs);
MatCreateVecs(H, &xs[0], NULL);
BSSPotR1Vec(bss, psi0, xs[0]);
EEPS eps; EEPSCreate(comm, &eps);
EEPSSetOperators(eps, H, S);
// EPSSetType(eps->eps, EPSJD);
EPSSetInitialSpace(eps->eps, 1, xs);
EEPSSetTarget(eps, -0.6);
// EPSSetInitialSpace(eps->eps, 1, xs);
EEPSSolve(eps);
int nconv;
PetscScalar kr;
EPSGetConverged(eps->eps, &nconv);
ASSERT_TRUE(nconv > 0);
EPSGetEigenpair(eps->eps, 0, &kr, NULL, NULL, NULL);
ASSERT_DOUBLE_NEAR(-0.5, kr, pow(10.0, -6.0));
Vec cs;
MatCreateVecs(H, &cs, NULL);
EEPSGetEigenvector(eps, 0, cs);
PetscReal x=1.1;
PetscScalar y=0.0;
PetscScalar dy=0.0;
BSSPsiOne(bss, cs, x, &y);
BSSDerivPsiOne(bss, cs, x, &dy);
ASSERT_DOUBLE_NEAR(creal(y), 2.0*x*exp(-x), pow(10.0, -6));
ASSERT_DOUBLE_NEAR(creal(dy), 2.0*exp(-x)-2.0*x*exp(-x), pow(10.0, -6));
VecDestroy(&xs[0]);
PetscFree(xs);
PFDestroy(&psi0);
BSSDestroy(&bss);
MatDestroy(&H);
MatDestroy(&V);
MatDestroy(&S);
EEPSDestroy(&eps);
VecDestroy(&cs);
return 0;
}示例5: remixmatrixu_
double remixmatrixu_(int*id, int*i,int*j){return creal(cMixMatrixU(*id,*i,*j));}示例6: c
/**
This code computes the integtated autocorrelation time :
It expects a file of length N double precision measurements
Defining c(t) = ( Y(t) - \bar{Y} )
C(T) = \sum_{t=0}^{N} ( c(t) c(t+T) )
R(T) = C(T) / C(0)
Setting a cutoff point "n"
Tau(n) = 0.5 + Nsep * \sum_{T}^{n} R(T)
We estimate the error on Tau(T) by
S_E(n) = n * \sqrt( ( 0.5 + \sum_{t}^{n} R(T) ) / N )
The computation of C(T) is performed by convolution with FFTs
The autocorrelation time is written to the file specified
*/
int
autocorrelation( const struct resampled RAW ,
const int NSEP ,
const char *output )
{
// openmp'd fftws
parallel_ffts( ) ;
if( RAW.restype != RAWDATA ) {
printf( "Resampled data is not RAW ... Cannot compute autocorrelation\n" ) ;
return FAILURE ;
}
// some constants
const int N = RAW.NSAMPLES ;
const int N2 = 2 * N ;
printf( "RAWDATA has %d samples\n" , N ) ;
printf( "Measurement separation %d\n\n" , NSEP ) ;
// allocate memory
double complex *in = calloc( N2 , sizeof( double complex ) ) ;
double complex *out = calloc( N2 , sizeof( double complex ) ) ;
// subtract the average from each data point
int i ;
#pragma omp parallel for private(i)
for( i = 0 ; i < N ; i++ ) {
in[ i ] = ( RAW.resampled[ i ] - RAW.avg ) ;
}
message( "FFT planning" ) ;
// are we doing this using openmp ffts?
#if ( defined OMP_FFTW ) && ( defined HAVE_OMP_H )
if( parallel_ffts( ) == FAILURE ) {
printf( "Parallel FFT setting failed \n" ) ;
return FAILURE ;
}
#endif
// create the plans
const fftw_plan forward = fftw_plan_dft_1d( N2 , in , out ,
FFTW_FORWARD , FFTW_ESTIMATE ) ;
const fftw_plan backward = fftw_plan_dft_1d( N2 , out , in ,
FFTW_BACKWARD , FFTW_ESTIMATE ) ;
fftw_execute( forward ) ;
// convolve
#pragma omp parallel for private(i)
for( i = 0 ; i < N2 ; i++ ) {
out[i] = creal( out[i] ) * creal( out[i] ) +
cimag( out[i] ) * cimag( out[i] ) ;
}
fftw_execute( backward ) ;
// normalise
const double zeropoint = 1.0 / in[ 0 ] ;
#pragma omp parallel for private(i)
for( i = 0 ; i < N2 ; i++ ) {
in[ i ] *= zeropoint ;
}
// summy the lags
message( "Computing tau(n)" ) ;
FILE *output_file = fopen( output , "w" ) ;
printf( "Writing tau(n) to file %s \n" , output ) ;
int n ;
for( n = 0 ; n < 30 ; n++ ) {
register double sum = 0.5 ;
//.........这里部分代码省略.........
示例7: VectorSphericalWaveFunctions
//.........这里部分代码省略.........
}else{
cph=*x/rho;
sph=*y/rho;
}
// Spherical Bessel Funtions
double JLM[*lmax+2];
// double *JLM=&JLM0[0];
gsl_sf_bessel_jl_steed_array(*lmax+1,*k*r,JLM);
/* CALCULATIONS OK
for(l=0;l<(*lmax+2);l++){
printf("%d\t%f\t%E\n",l,*k*r,JLM[l]);
}
*/
// Qlm - primeiros 4 termos
double Qlm[LMAXE];
Qlm[jlm(0, 0)]=1/sqrt(4*M_PI);
Qlm[jlm(1, 1)]=-gammaQ(1)*sth*Qlm[jlm(0,0)]; // Q11
Qlm[jlm(1, 0)]=sqrt(3.0)*cth*Qlm[jlm(0,0)]; // Q10
Qlm[jlm(1,-1)]=-Qlm[jlm(1,1)]; // Q11*(-1)
// Complex Exponencial for m=-1,0,1
double complex Eim[2*(*lmax)+3];
Eim[*lmax-1]=(cph-I*sph);
Eim[*lmax ]=1+I*0;
Eim[*lmax+1]=(cph+I*sph);
// Ylm - primeiros 4 termos
double complex Ylm[LMAXE];
Ylm[jlm(0, 0)]=Qlm[jlm(0, 0)];
Ylm[jlm(1,-1)]=Qlm[jlm(1,-1)]*Eim[*lmax-1];
Ylm[jlm(1, 0)]=Qlm[jlm(1, 0)];
Ylm[jlm(1, 1)]=Qlm[jlm(1, 1)]*Eim[*lmax+1];
/* OK jl, Qlm, Ylm
for(l=0;l<2;l++){
for(m=-l;m<=l;m++){
printf("%d\t%d\t%d\t%f\t%f\t%f+%fi\n",l,m,jlm(l,m),JLM[jlm(l,m)],Qlm[jlm(l,m)],creal(Ylm[jlm(l,m)]),cimag(Ylm[jlm(l,m)]));
}
}
printf("======================================================================\n");
*/
// VECTOR SPHERICAL HARMONICS
double complex XM; //[LMAX];
double complex XZ; //[LMAX];
double complex XP; //[LMAX];
double complex YM; //[LMAX];
double complex YZ; //[LMAX];
double complex YP; //[LMAX];
double complex VM; //[LMAX];
double complex VZ; //[LMAX];
double complex VP; //[LMAX];
// HANSEN MULTIPOLES
double complex MM; //[LMAX];
double complex MZ; //[LMAX];
double complex MP; //[LMAX];
double complex NM; //[LMAX];
double complex NZ; //[LMAX];
double complex NP; //[LMAX];
// OTHERS
double kl;
// MAIN LOOP
for(l=1;l<=(*lmax);l++){
//------------------------------------------------------------------------------
//Qlm extremos positivos um passo a frente
Qlm[jlm(l+1, l+1)]=-gammaQ(l+1)*sth*Qlm[jlm(l,l)];
Qlm[jlm(l+1, l )]= deltaQ(l+1)*cth*Qlm[jlm(l,l)];
//Qlm extremos negativos um passo a frente
Qlm[jlm(l+1,-l-1)]=pow(-1,l+1)*Qlm[jlm(l+1, l+1)];
Qlm[jlm(l+1,-l )]=pow(-1,l )*Qlm[jlm(l+1, l )];示例8: cerfc
cmplxd cerfc(cmplxd x)
{
static const double pv= 1.27813464856668857e+01;
static const double ph= 6.64067324283344283e+00;
static const double p0= 2.94608570191793668e-01;
static const double p1= 1.81694307871527086e-01;
static const double p2= 6.91087778921425355e-02;
static const double p3= 1.62114197106901582e-02;
static const double p4= 2.34533471539159422e-03;
static const double p5= 2.09259199579049675e-04;
static const double p6= 1.15149016557480535e-05;
static const double p7= 3.90779571296927748e-07;
static const double p8= 8.17898509247247602e-09;
static const double p9= 1.05575446466983499e-10;
static const double p10= 8.40470321453263734e-13;
static const double p11= 4.12646136715431977e-15;
static const double p12= 1.24947948599560084e-17;
static const double q0= 6.04152433382652546e-02;
static const double q1= 5.43737190044387291e-01;
static const double q2= 1.51038108345663136e+00;
static const double q3= 2.96034692357499747e+00;
static const double q4= 4.89363471039948562e+00;
static const double q5= 7.31024444393009580e+00;
static const double q6= 1.02101761241668280e+01;
static const double q7= 1.35934297511096823e+01;
static const double q8= 1.74600053247586586e+01;
static const double q9= 2.18099028451137569e+01;
static const double q10= 2.66431223121749773e+01;
static const double q11= 3.19596637259423197e+01;
static const double q12= 3.77595270864157841e+01;
static const double r0= 1.56478036351085356e-01;
static const double r1= 2.45771407110492625e-01;
static const double r2= 1.19035163906534275e-01;
static const double r3= 3.55561834455977740e-02;
static const double r4= 6.55014550718381002e-03;
static const double r5= 7.44188068433574137e-04;
static const double r6= 5.21447276257559040e-05;
static const double r7= 2.25337799750608244e-06;
static const double r8= 6.00556181041662576e-08;
static const double r9= 9.87118243564461826e-10;
static const double r10= 1.00064645539515792e-11;
static const double r11= 6.25587539334288736e-14;
static const double r12= 2.41207864479170276e-16;
static const double s1= 2.41660973353061018e-01;
static const double s2= 9.66643893412244073e-01;
static const double s3= 2.17494876017754917e+00;
static const double s4= 3.86657557364897629e+00;
static const double s5= 6.04152433382652546e+00;
static const double s6= 8.69979504071019666e+00;
static const double s7= 1.18413876942999899e+01;
static const double s8= 1.54663022945959052e+01;
static const double s9= 1.95745388415979425e+01;
static const double s10= 2.41660973353061018e+01;
static const double s11= 2.92409777757203832e+01;
static const double s12= 3.47991801628407866e+01;
cmplxd y=x*x;
if(cabs(creal(x))+cabs(cimag(x)) < ph) {
const cmplxd z=cexp(pv*x);
if(creal(z) >= 0.) {
y = cexp(-y)*x*(p12/(y+q12)+p11/(y+q11)
+p10/(y+q10)+p9/(y+q9)+p8/(y+q8)+p7/(y+q7)
+p6/(y+q6)+p5/(y+q5)+p4/(y+q4)+p3/(y+q3)
+p2/(y+q2)+p1/(y+q1)+p0/(y+q0))+2./(1.+z);
} else {
y = cexp(-y)*x*(r12/(y+s12)+r11/(y+s11)
+r10/(y+s10)+r9/(y+s9)+r8/(y+s8)+r7/(y+s7)
+r6/(y+s6)+r5/(y+s5)+r4/(y+s4)+r3/(y+s3)
+r2/(y+s2)+r1/(y+s1)+r0/y)+2./(1.-z);
}
} else {
y = cexp(-y)*x*(p12/(y+q12)+p11/(y+q11)
+p10/(y+q10)+p9/(y+q9)+p8/(y+q8)+p7/(y+q7)
+p6/(y+q6)+p5/(y+q5)+p4/(y+q4)+p3/(y+q3)
+p2/(y+q2)+p1/(y+q1)+p0/(y+q0));
if(creal(x) <= 0) y=y+2.;
}
return y;
} 示例9: XLALCompareCOMPLEX8Vectors
/** Compare two COMPLEX8 vectors using various different comparison metrics
*/
int
XLALCompareCOMPLEX8Vectors ( VectorComparison *result, ///< [out] return comparison results
const COMPLEX8Vector *x, ///< [in] first input vector
const COMPLEX8Vector *y, ///< [in] second input vector
const VectorComparison *tol ///< [in] accepted tolerances on comparisons, or NULL for no check
)
{
XLAL_CHECK ( result != NULL, XLAL_EINVAL );
XLAL_CHECK ( x != NULL, XLAL_EINVAL );
XLAL_CHECK ( y != NULL, XLAL_EINVAL );
XLAL_CHECK ( x->data != NULL, XLAL_EINVAL );
XLAL_CHECK ( y->data != NULL, XLAL_EINVAL );
XLAL_CHECK ( x->length > 0, XLAL_EINVAL );
XLAL_CHECK ( x->length == y->length, XLAL_EINVAL );
REAL8 x_L1 = 0, x_L2 = 0;
REAL8 y_L1 = 0, y_L2 = 0;
REAL8 diff_L1 = 0, diff_L2 = 0;
COMPLEX16 scalar = 0;
REAL8 maxAbsx = 0, maxAbsy = 0;
COMPLEX8 x_atMaxAbsx = 0, y_atMaxAbsx = 0;
COMPLEX8 x_atMaxAbsy = 0, y_atMaxAbsy = 0;
UINT4 numSamples = x->length;
for ( UINT4 i = 0; i < numSamples; i ++ )
{
COMPLEX8 x_i = x->data[i];
COMPLEX8 y_i = y->data[i];
REAL8 xAbs_i = cabs ( x_i );
REAL8 yAbs_i = cabs ( y_i );
XLAL_CHECK ( isfinite ( xAbs_i ), XLAL_EFPINVAL, "non-finite element: x(%d) = %g + I %g\n", i, crealf(x_i), cimagf(x_i) );
XLAL_CHECK ( isfinite ( yAbs_i ), XLAL_EFPINVAL, "non-finite element: y(%d) = %g + I %g\n", i, crealf(y_i), cimagf(y_i) );
REAL8 absdiff = cabs ( x_i - y_i );
diff_L1 += absdiff;
diff_L2 += SQ(absdiff);
x_L1 += xAbs_i;
y_L1 += yAbs_i;
x_L2 += SQ(xAbs_i);
y_L2 += SQ(yAbs_i);
scalar += x_i * conj(y_i);
if ( xAbs_i > maxAbsx ) {
maxAbsx = xAbs_i;
x_atMaxAbsx = x_i;
y_atMaxAbsx = y_i;
}
if ( yAbs_i > maxAbsy ) {
maxAbsy = yAbs_i;
x_atMaxAbsy = x_i;
y_atMaxAbsy = y_i;
}
} // for i < numSamples
// complete L2 norms by taking sqrt
x_L2 = sqrt ( x_L2 );
y_L2 = sqrt ( y_L2 );
diff_L2 = sqrt ( diff_L2 );
// compute and return comparison results
result->relErr_L1 = diff_L1 / ( 0.5 * (x_L1 + y_L1 ) );
result->relErr_L2 = diff_L2 / ( 0.5 * (x_L2 + y_L2 ) );
REAL8 cosTheta = fmin ( 1, creal ( scalar ) / (x_L2 * y_L2) );
result->angleV = acos ( cosTheta );
result->relErr_atMaxAbsx = cRELERR ( x_atMaxAbsx, y_atMaxAbsx );
result->relErr_atMaxAbsy = cRELERR ( x_atMaxAbsy, y_atMaxAbsy );;
XLAL_CHECK ( XLALCheckVectorComparisonTolerances ( result, tol ) == XLAL_SUCCESS, XLAL_EFUNC );
return XLAL_SUCCESS;
} // XLALCompareCOMPLEX8Vectors()示例10: simmap_covar_matrix_complex
// stationary covariance for Complex eigenvalues
static void simmap_covar_matrix_complex (int *nchar,
double *bt,
Rcomplex *lambda_val,
Rcomplex *S_val,
Rcomplex *S1_val,
double *sigmasq,
int *nterm,
double *V) {
// complex version
double complex *eltj, *elti, *W, *U, *tmp1, *lambda, *S, *S1;
double sij, ti, tj, tmp2;
int n = *nchar, nt = *nterm;
int i, j, k, l, r, s;
// alloc complex vectors
U = Calloc(n*n,double complex);
W = Calloc(n*n,double complex);
tmp1 = Calloc(n*n,double complex);
eltj = Calloc(n,double complex);
elti = Calloc(n,double complex);
S = Calloc(n*n,double complex);
S1 = Calloc(n*n,double complex);
lambda = Calloc(n,double complex);
//zeroing vectors & transform to C complex structure
for(i = 0; i<n; i++){
lambda[i]=comp(lambda_val[i]);
for(j =0; j<n; j++){
S[i+j*n]=comp(S_val[i+j*n]);
S1[i+j*n]=comp(S1_val[i+j*n]);
U[i+j*n] = 0;
W[i+j*n] = 0;
}
}
// computing the P^-1%*%Sigma%*%P^-T
for (i = 0; i < n; i++) {
for (j = 0; j < n; j++) {
for (k = 0; k < n; k++) {
for (l = 0; l < n; l++) {
U[i+j*n] += S1[i+k*n]*sigmasq[k+l*n]*S1[j+l*n];
}
}
}
}
// fill in the covariance
for (i = 0; i < nt; i++) {
for (j = 0; j <= i; j++) {
ti = bt[i+i*nt];
sij = bt[i+j*nt];
tj = bt[j+j*nt];
// complex exponential with time
for (k = 0; k < n; k++) {
elti[k] = cexp(-lambda[k]*(ti-sij));
eltj[k] = cexp(-lambda[k]*(tj-sij));
}
// Integral parts
for (k = 0; k < n; k++) {
for (l = 0; l < n; l++) {
W[k+l*n] = elti[k]*U[k+l*n]*eltj[l]/(lambda[k]+lambda[l]);
}
}
// computing the P%*%Sigma%*%P^T
for (r = 0; r < n; r++) {
for (s = 0; s < n; s++) {
tmp1[r+s*n] = 0;
for (k = 0; k < n; k++) {
for (l = 0; l < n; l++) {
tmp1[r+s*n] += S[r+k*n]*W[k+l*n]*S[s+l*n];
}
}
}
}
for (k = 0; k < n; k++) {
for (l = 0; l < n; l++) {
// Save the real parts
tmp2 = creal(tmp1[k+l*n]);
V[i+nt*(k+n*(j+nt*l))] = tmp2;
if (j != i)
V[j+nt*(l+n*(i+nt*k))] = tmp2;
}
}
// End
}
}
Free(lambda);
Free(S);
Free(S1);
//.........这里部分代码省略.........
示例11:
JL_DLLEXPORT complex float *cfptest(complex float *a) {
//Unpack a ComplexPair{Float64} struct
if (verbose) fprintf(stderr,"%g + %g i\n", creal(*a), cimag(*a));
*a += 1 - (2.0*I);
return a;
}示例12: cftest
JL_DLLEXPORT complex float cftest(complex float a) {
//Unpack a ComplexPair{Float32} struct
if (verbose) fprintf(stderr,"%g + %g i\n", creal(a), cimag(a));
a += 1 - (2.0*I);
return a;
}示例13: CreateGeometry
Geometry CreateGeometry(int N[3], double h[3], int Npml[3], int Nc, int LowerPML, double *eps, double *epsI, double *fprof, double wa, double y){
int i;
Geometry geo = (Geometry) malloc(sizeof(struct Geometry_s));
geo->Nc = Nc;
geo->LowerPML = LowerPML;
geo->interference = 0.0; // default no interference
for(i=0; i<3; i++){
geo->h[i] = h[i];
geo->Npml[i] = Npml[i];
}
CreateGrid(&geo->gN, N, geo->Nc, 2);
CreateGrid(&geo->gM, N, 1, 1); // 3/3/14: set M = N as per Steven
CreateVec(2*Nxyzc(geo)+2, &geo->vepspml);
int manual_epspml = 0;
PetscOptionsGetInt(PETSC_NULL,PETSC_NULL,"-manual_epspml", &manual_epspml, NULL);
if(manual_epspml == 0){
Vecfun pml;
CreateVecfun(&pml,geo->vepspml);
for(i=pml.ns; i<pml.ne; i++){
Point p;
CreatePoint_i(&p, i, &geo->gN);
project(&p, 3);
dcomp eps_geoal;
eps_geoal = pmlval(xyzc(&p), N, geo->Npml, geo->h, geo->LowerPML, 0);
setr(&pml, i, p.ir? cimag(eps_geoal) : creal(eps_geoal) );
}
DestroyVecfun(&pml);
}
CreateVec(Mxyz(geo), &geo->vMscratch[0]);
for(i=0; i<SCRATCHNUM; i++){
geo->vNhscratch[i] = 0; // allows checking whether vN created or not
if(i>0)VecDuplicate(geo->vMscratch[0], &geo->vMscratch[i]);
}
double *scratch;
int ms, me;
VecGetOwnershipRange(geo->vMscratch[0], &ms, &me);
if( !manual_epspml){
VecGetArray(geo->vMscratch[0], &scratch);
for(i=ms; i<me;i++)
scratch[i-ms] = eps[i-ms];
VecRestoreArray(geo->vMscratch[0], &scratch);
}
CreateVec(2*Nxyzc(geo)+2, &geo->vH);
VecDuplicate(geo->vH, &geo->veps);
VecDuplicate(geo->vH, &geo->vIeps);
for(i=0; i<SCRATCHNUM; i++) VecDuplicate(geo->vH, &geo->vscratch[i]);
VecSet(geo->vH, 1.0);
if( !manual_epspml){
VecShift(geo->vMscratch[0], -1.0); //hack, for background dielectric
InterpolateVec(geo, geo->vMscratch[0], geo->vscratch[1]);
VecShift(geo->vscratch[1], 1.0);
VecPointwiseMult(geo->veps, geo->vscratch[1], geo->vepspml);
if(epsI != NULL){ // imaginary part of passive dielectric
VecGetArray(geo->vMscratch[0], &scratch);
for(i=ms; i<me; i++){
scratch[i-ms] = epsI[i-ms];
}
VecRestoreArray(geo->vMscratch[0], &scratch);
InterpolateVec(geo, geo->vMscratch[0], geo->vscratch[1]);
VecPointwiseMult(geo->vscratch[1], geo->vscratch[1], geo->vepspml);
TimesI(geo, geo->vscratch[1], geo->vscratch[2]);
VecAXPY(geo->veps, 1.0, geo->vscratch[2]);
}
}
if(manual_epspml){
char epsManualfile[PETSC_MAX_PATH_LEN];
PetscOptionsGetString(PETSC_NULL,PETSC_NULL,"-epsManualfile", epsManualfile, PETSC_MAX_PATH_LEN, NULL);
FILE *fp = fopen(epsManualfile, "r");
ReadVectorC(fp, 2*Nxyzc(geo)+2, geo->veps);
// 07/11/15: if manual_epspml, then directly read in the Nxyzcr+2 vector
fclose(fp);
}
TimesI(geo, geo->veps, geo->vIeps);
// vIeps for convenience only, make sure to update it later if eps ever changes!
geo->D = 0.0;
geo->wa = wa;
geo->y = y;
VecDuplicate(geo->veps, &geo->vf);
VecDuplicate(geo->vMscratch[0], &geo->vfM);
//.........这里部分代码省略.........
示例14: main
int
main(void)
{
const int n = 1000;
int i;
double _Complex vresult, result, array[n];
bool lvresult, lresult;
for (i = 0; i < n; i++)
array[i] = i;
result = 0;
vresult = 0;
/* '+' reductions. */
#pragma acc parallel vector_length (vl)
#pragma acc loop reduction (+:result)
for (i = 0; i < n; i++)
result += array[i];
/* Verify the reduction. */
for (i = 0; i < n; i++)
vresult += array[i];
if (result != vresult)
abort ();
result = 0;
vresult = 0;
/* Needs support for complex multiplication. */
// /* '*' reductions. */
// #pragma acc parallel vector_length (vl)
// #pragma acc loop reduction (*:result)
// for (i = 0; i < n; i++)
// result *= array[i];
//
// /* Verify the reduction. */
// for (i = 0; i < n; i++)
// vresult *= array[i];
//
// if (fabs(result - vresult) > .0001)
// abort ();
// result = 0;
// vresult = 0;
// /* 'max' reductions. */
// #pragma acc parallel vector_length (vl)
// #pragma acc loop reduction (+:result)
// for (i = 0; i < n; i++)
// result = result > array[i] ? result : array[i];
//
// /* Verify the reduction. */
// for (i = 0; i < n; i++)
// vresult = vresult > array[i] ? vresult : array[i];
//
// printf("%d != %d\n", result, vresult);
// if (result != vresult)
// abort ();
//
// result = 0;
// vresult = 0;
//
// /* 'min' reductions. */
// #pragma acc parallel vector_length (vl)
// #pragma acc loop reduction (+:result)
// for (i = 0; i < n; i++)
// result = result < array[i] ? result : array[i];
//
// /* Verify the reduction. */
// for (i = 0; i < n; i++)
// vresult = vresult < array[i] ? vresult : array[i];
//
// printf("%d != %d\n", result, vresult);
// if (result != vresult)
// abort ();
result = 5;
vresult = 5;
lresult = false;
lvresult = false;
/* '&&' reductions. */
#pragma acc parallel vector_length (vl)
#pragma acc loop reduction (&&:lresult)
for (i = 0; i < n; i++)
lresult = lresult && (creal(result) > creal(array[i]));
/* Verify the reduction. */
for (i = 0; i < n; i++)
lvresult = lresult && (creal(result) > creal(array[i]));
if (lresult != lvresult)
abort ();
result = 5;
vresult = 5;
//.........这里部分代码省略.........
示例15: fftMeasure
static void
fftMeasure (int nframes, int overlap, float *indata)
{
#ifdef USE_FFTW
int i, stepSize = fftSize/overlap;
double freqPerBin = rate/(double)fftSize,
phaseDifference = 2.*M_PI*(double)stepSize/(double)fftSize;
if (!fftSample) fftSample = fftSampleBuffer + (fftSize-stepSize);
// bzero(fftGraphData.db, GRAPH_MAX_FREQ);
for (i=0; i<nframes; i++) {
*fftSample++ = indata[i];
if (fftSample-fftSampleBuffer >= fftSize) {
int k;
Peak peaks[MAX_PEAKS];
for (k=0; k<MAX_PEAKS; k++) {
peaks[k].db = -200.;
peaks[k].freq = 0.;
}
fftSample = fftSampleBuffer + (fftSize-stepSize);
for (k=0; k<fftSize; k++) {
double window = -.5*cos(2.*M_PI*(double)k/(double)fftSize)+.5;
fftIn[k] = fftSampleBuffer[k] * window;
}
fftwf_execute(fftPlan);
for (k=0; k<=fftSize/2; k++) {
long qpd;
float
real = creal(fftOut[k]),
imag = cimag(fftOut[k]),
magnitude = 20.*log10(2.*sqrt(real*real + imag*imag)/fftSize),
phase = atan2(imag, real),
tmp, freq;
/* compute phase difference */
tmp = phase - fftLastPhase[k];
fftLastPhase[k] = phase;
/* subtract expected phase difference */
tmp -= (double)k*phaseDifference;
/* map delta phase into +/- Pi interval */
qpd = tmp / M_PI;
if (qpd >= 0) qpd += qpd&1;
else qpd -= qpd&1;
tmp -= M_PI*(double)qpd;
/* get deviation from bin frequency from the +/- Pi interval */
tmp = overlap*tmp/(2.*M_PI);
/* compute the k-th partials' true frequency */
freq = (double)k*freqPerBin + tmp*freqPerBin;
#ifdef USE_GRAPH
int fi = (int)round(freq);
if (fi < GRAPH_MAX_FREQ) {
fftGraphData.db[fi] = magnitude;
if (magnitude > fftGraphData.dbMax)
fftGraphData.dbMax = magnitude;
if (magnitude < fftGraphData.dbMin)
fftGraphData.dbMin = magnitude;
}
/*
printf("%+8.3f % .5f %i\n", freq, fftGraphData.scale_freq,
(int)round(freq * fftGraphData.scale_freq));
*/
#endif
if (freq > 0.0 && magnitude > peaks[0].db) {
memmove(peaks+1, peaks, sizeof(Peak)*(MAX_PEAKS-1));
peaks[0].freq = freq;
peaks[0].db = magnitude;
}
}
fftFrameCount++;
if (fftFrameCount > 0 && fftFrameCount % overlap == 0) {
int l, maxharm = 0;
k = 0;
for (l=1; l<MAX_PEAKS && peaks[l].freq > 0.0; l++) {
int harmonic;
for (harmonic=5; harmonic>1; harmonic--) {
if (peaks[0].freq / peaks[l].freq < harmonic+.02 &&
peaks[0].freq / peaks[l].freq > harmonic-.02) {
if (harmonic > maxharm &&
peaks[0].db < peaks[l].db/2) {
maxharm = harmonic;
k = l;
}
}
}
//displayFrequency(&peaks[l], lblFreq[l]);
}
//.........这里部分代码省略.........