C++ sf_cmul函数代码示例

void cbanded_solve (sf_complex *b)
/*< multiply by inverse (in place) >*/
{
    int k, m;
    sf_complex t;

    for (k = 1; k < band; k++) {
	t = b[k];
	for (m = 1; m <= k; m++) {
#ifdef SF_HAS_COMPLEX_H
	    t -= o[m-1][k-m] * b[k-m];
#else
	    t = sf_csub(t,sf_cmul(o[m-1][k-m],b[k-m]));
#endif
	}
	b[k] = t;
    }
    for (k = band; k < n; k++) {
	t = b[k];
	for (m = 1; m <= band; m++) {
#ifdef SF_HAS_COMPLEX_H
	    t -= o[m-1][k-m] * b[k-m];
#else
	    t = sf_csub(t,sf_cmul(o[m-1][k-m],b[k-m]));
#endif
	}
	b[k] = t;
    }
    for (k = n-1; k >= n - band; k--) {
#ifdef SF_HAS_COMPLEX_H
	t = b[k]/d[k];
#else
	t = sf_crmul(b[k],1./d[k]);
#endif
	for (m = 0; m < n - k - 1; m++) {
#ifdef SF_HAS_COMPLEX_H
	    t -= conjf(o[m][k]) * b[k+m+1];
#else
	    t = sf_csub(t,sf_cmul(conjf(o[m][k]),b[k+m+1]));
#endif
	}
	b[k] = t;
    }
    for (k = n - band - 1; k >= 0; k--) {
#ifdef SF_HAS_COMPLEX_H
	t = b[k]/d[k];
#else
	t = sf_crmul(b[k],1./d[k]);
#endif
	for (m = 0; m < band; m++) {
#ifdef SF_HAS_COMPLEX_H
	    t -= conjf(o[m][k]) * b[k+m+1];
#else
	    t = sf_csub(t,sf_cmul(conjf(o[m][k]),b[k+m+1]));
#endif
	}
	b[k] = t;
    }
}

void ctoeplitz_solve (const sf_complex *r /* top row of the matrix */, 
		      sf_complex *f       /* inverted in place */)
/*< apply the solver >*/
{    
    int i,j;
    sf_complex e,c,w, bot;
    float v;
    
    v=crealf(r[0]);
#ifdef SF_HAS_COMPLEX_H
    f[0] /= v;
#else
    f[0] = sf_crmul(f[0],1./v);
#endif
    
    for (j=1; j < n; j++) {
	e = cdprod(j,a,r);
#ifdef SF_HAS_COMPLEX_H
	c = -e/v;
#else
	c = sf_crmul(e,-1./v);
#endif

	v += crealf(c)*crealf(e) + cimagf(c)*cimagf(e);
       
	for (i=1; i <= j/2; i++) {
#ifdef SF_HAS_COMPLEX_H
	    bot  = a[j-i] + c*conjf(a[i]);
	    a[i] += c*conjf(a[j-i]);
#else
	    bot  = sf_cadd(a[j-i],sf_cmul(c,conjf(a[i])));
	    a[i] = sf_cadd(a[i],sf_cmul(c,conjf(a[j-i])));
#endif
	    a[j-i] = bot;
	}
	a[j] = c;
       
	w = cdprod(j,f,r);
#ifdef SF_HAS_COMPLEX_H
	c = (f[j]-w)/v;
#else
	c = sf_crmul(sf_csub(f[j],w),1./v);
#endif
       
	for (i=0; i < j; i++) {
#ifdef SF_HAS_COMPLEX_H
	    f[i] += c*conjf(a[j-i]);
#else
	    f[i] = sf_cadd(f[i],sf_cmul(c,conjf(a[j-i])));
#endif
	}
	f[j] = c;
    }
}

void radon_lop (bool adj, bool add, int nm, int nd, 
		sf_complex *mm, sf_complex *dd)
/*< linear operator >*/
{
    int ix, ip;
    sf_complex c, d;

    if (nm != np || nd != nx) 
	sf_error("%s: mismatched data sizes",__FILE__);

    sf_cadjnull(adj, add, nm, nd, mm, dd);
    
    for (ix=0; ix < nx; ix++) {
	if (adj == true) {
	    c = dc[ix];
#ifdef SF_HAS_COMPLEX_H
	    d = dd[ix]*c0[ix];
#else
	    d = sf_cmul(dd[ix],c0[ix]);
#endif
	    for (ip=0; ip < np-1; ip++) {
#ifdef SF_HAS_COMPLEX_H
		mm[ip] += d;
		d *= c; 
#else
		mm[ip] = sf_cadd(mm[ip],d);
		d = sf_cmul(d,c); 
#endif
	    }
#ifdef SF_HAS_COMPLEX_H
	    mm[np-1] += d;
#else
	    mm[np-1] = sf_cadd(mm[np-1],d);
#endif
	} else {
	    c = conjf(dc[ix]);
	    d = mm[np-1];
	    for (ip=np-2; ip >= 0; ip--) {
#ifdef SF_HAS_COMPLEX_H
		d = d*c + mm[ip];
#else
		d = sf_cadd(sf_cmul(d,c),mm[ip]);
#endif
	    }
#ifdef SF_HAS_COMPLEX_H
	    dd[ix] += d*conjf(c0[ix]);
#else
	    dd[ix] = sf_cadd(dd[ix],sf_cmul(d,conjf(c0[ix])));
#endif
	}
    }
}

void BornScatteredField(float xx,float xy,float xz,float *output)
/*< Born scattering field >*/
{
    int iw;
    double omega,scale;

    sf_complex  U,dU;
    sf_complex val,fkern;

    ZeroArray(u,  nt); 
    ZeroArray(du, nt); 
    *output=0.;

    for (iw=0; iw<nw; iw++) {
	omega=ow+dw*iw;

	scale =cos((SF_PI/2)*(((double) iw+1)/((double) nw+1)));
	scale*=scale*omega;

#ifdef SF_HAS_COMPLEX_H
	/* background field */
	U=Green(xx,xy,xz,sx,sy,sz,omega)*scale;
	    
	/* scattered field */
	val=Green(px,py,pz,sx,sy,sz,omega)*scale;
	val *= Green(xx,xy,xz,px,py,pz,omega);
	dU = val*(-omega*omega*dv);

	U += dU;
	fkern=cexpf(sf_cmplx(0.,-omega*0.6));
	val=fkern*U;
#else
	/* background field */
	U=sf_crmul(Green(xx,xy,xz,sx,sy,sz,omega),scale);
	    
	/* scattered field */
	val=sf_crmul(Green(px,py,pz,sx,sy,sz,omega),scale);
	val=sf_cmul(val,Green(xx,xy,xz,px,py,pz,omega));
	dU=sf_crmul(val,-omega*omega*dv);

	U=sf_cadd(U,dU);
	fkern=cexpf(sf_cmplx(0.,-omega*0.6));
	val=sf_cmul(fkern,U);
#endif

	*output+=crealf(val);
    }
    return;
}    

void radon_toep (sf_complex *qq /* Toeplitz row */, 
		 float eps         /* regularization */)
/*< fill Toeplitz matrix for inversion >*/
{
    int ix, ip;
    sf_complex c, q;
    
    qq[0] = sf_cmplx(eps*eps,0.);
    for (ip=1; ip < np; ip++) {
	qq[ip] = czero;
    }

    for (ix=0; ix < nx; ix++) {
	c = conjf(dc[ix]);
	q = sf_cmplx(1.,0.);
	for (ip=0; ip < np-1; ip++) {
#ifdef SF_HAS_COMPLEX_H
	    qq[ip] += q;
	    q *= c;
#else
	    qq[ip] = sf_cadd(qq[ip],q);
	    q = sf_cmul(q,c);
#endif
	}
#ifdef SF_HAS_COMPLEX_H
	qq[np-1] += q;
#else
	qq[np-1] = sf_cadd(qq[np-1],q);
#endif
    }
}

/*------------------------------------------------------------*/
void cam3_ssf(
    sf_complex       w /* frequency */,
    sf_complex   ***wx /* wavefield */,
    cub3d          cub,
    cam3d          cam,
    tap3d          tap,
    slo3d          slo,
    int            imz,
    int         ompith
    )
/*< Wavefield extrapolation by SSF >*/
{
    sf_complex co,cs,cr,w2,khy,kss,krr;
    float s, kmy, d,dsc,drc;
    int imy,imx,ihx,jmy, js,jr;

    co = sf_cmplx(0,0);
#ifdef SF_HAS_COMPLEX_H
    w2 = w*w;
#else
    w2 = sf_cmul(w,w);
#endif
    
#ifdef SF_HAS_COMPLEX_H
    LOOP( s = slo->so[ompith][ cam->jy[imy] ][ cam->is[ihx][imx] ] 
	  +   slo->so[ompith][ cam->jy[imy] ][ cam->ir[ihx][imx] ];
	  wx[ihx][imy][imx] *= cexpf(-w*s* cub->amz.d/2); );

/*------------------------------------------------------------*/
void sft2(sf_complex **pp)
/*< apply shift >*/
{
    int i1,i2;

    for    (i2=0; i2<n2; i2++) {
	for(i1=0; i1<n1; i1++) {
#ifdef SF_HAS_COMPLEX_H
	    pp[i2][i1] *= shf1[i1]*shf2[i2];
#else
	    pp[i2][i1] = sf_cmul(pp[i2][i1],
				 sf_cmul(shf1[i1],shf2[i2])
		);
#endif
	}
    }    
}

/*------------------------------------------------------------*/
void ompsft2(sf_complex **pp,
	     fft2d fft)
/*< apply shift >*/
{
    int i1,i2;

    for     (i2=0; i2<fft->n2; i2++) {
	for (i1=0; i1<fft->n1; i1++) {
#ifdef SF_HAS_COMPLEX_H
	    pp[i2][i1] *= fft->shf1[i1]*fft->shf2[i2];
#else
	    pp[i2][i1] = sf_cmul(pp[i2][i1],
				 sf_cmul(fft->shf1[i1],fft->shf2[i2]));
#endif
	}
    }    
}

void fft1_2D_fwd (float *input, float *output, int n2)
/*< Fast Fourier Transform along the first axis forward>*/
{
    /*bool inv, sym, opt;*/
    int n1, i1, i2;
    float  *p,  wt, shift;
    kiss_fft_cpx *pp, ce;
    char *label;
    sf_file out=NULL;
    kiss_fftr_cfg cfg;
	bool verb;

	verb=1;

    
    
  wt = sym? 1./sqrtf((float) ntopt): 1.0/ntopt;

  p = sf_floatalloc(ntopt);
  pp = (kiss_fft_cpx*) sf_complexalloc(nw);

  cfg = kiss_fftr_alloc(ntopt,0,NULL,NULL);
    

	for (i2=0; i2 < n2; i2++) {

    /*sf_floatread (p,n1,in);*/
    memcpy(p,&input[i2*nt],nt*sizeof(float));
    if (sym) {
      for (i1=0; i1 < nt; i1++) {
    	  p[i1] *= wt;
      }
    }
    for (i1=nt; i1 < ntopt; i1++) {
      p[i1]=0.0; // padding
    }
    
      kiss_fftr (cfg,p,pp);
    
    if (0. != ot) {
      for (i1=0; i1 < nw; i1++) {
          shift = -2.0*SF_PI*i1*dw*ot;
          ce.r = cosf(shift);
          ce.i = sinf(shift);
          pp[i1]=sf_cmul(pp[i1],ce);
      }
    }
    
    memcpy(&output[i2*2*nw],pp,2*nw*sizeof(float));
  }
    free(p);
    free(pp);

	/*exit (0);*/
}

void rweone_phs(
    sf_complex *v,
    float w,
    float a0,
    float b0
    )
/*< Fourier-domain phase shift >*/
{
    int ig,ikg;
    float kg;

    float a2,b2,k2;
    sf_complex iw,ikt,w2;

    a2 = a0*a0;
    b2 = b0*b0;

    iw = sf_cmplx(2e-3,-w);
#ifdef SF_HAS_COMPLEX_H
    w2 = iw*iw;
#else
    w2 = sf_cmul(iw,iw);
#endif

    rweone_fft(false,(kiss_fft_cpx*) v);

    for(ig=0;ig<ag.n;ig++) {
	ikg = KMAP(ig,ag.n);
	kg = okg + ikg * dkg;
	k2 = kg*kg;
	
#ifdef SF_HAS_COMPLEX_H
	ikt = csqrtf( w2*a2 + k2*b2 );
	v[ig] *= cexpf(-ikt*at.d);
#else
	ikt = csqrtf(sf_cadd(sf_crmul(w2,a2),sf_cmplx(k2*b2,0.)));
	v[ig] = sf_cmul(v[ig],cexpf(sf_crmul(ikt,-at.d)));
#endif
    }
    
    rweone_fft( true,(kiss_fft_cpx*) v);
}

void cmatmul(kiss_fft_cpx *a,kiss_fft_cpx *b,int m,int n,int k,kiss_fft_cpx *c) 
/*< complex matrix multiplication: a->m*n, b->n*k >*/
{
    int i,j,l,u;
    for (i=0; i<=m-1; i++)
	for (j=0; j<=k-1; j++)	
	{ 
	    u=i*k+j; c[u].r=0.0; c[u].i=0.0;
	    for (l=0; l<=n-1; l++)
		c[u]=sf_cadd(c[u],sf_cmul(a[i*n+l],b[l*k+j]));
	}
}

void cicai1_lop (bool adj, bool add, int nx, int ny, sf_complex* xx, sf_complex* yy) 
/*< linear operator >*/
{
    int b, x, y;

    if(ny != nx) sf_error("%s: size problem: %d != %d",__FILE__,ny,nx);
    sf_cadjnull (adj, add, nx, ny, xx, yy);
    
    for( b=0; b < nb; b++) {
	for( y = nb - lg; y <= ny - lg; y++) {
	    x = y - b + lg - 1;
#ifdef SF_HAS_COMPLEX_H
	    if( adj) xx[x] += yy[y] * conjf(bb[b]);
	    else     yy[y] += xx[x] * bb[b];
#else
	    if( adj) xx[x] = sf_cadd(xx[x],sf_cmul(yy[y],conjf(bb[b])));
	    else     yy[y] = sf_cadd(yy[y],sf_cmul(xx[x],bb[b]));
#endif
	}
    }
}

void rweone_thr(
    sf_complex *a,
    sf_complex *b,
    sf_complex *c,
    sf_complex *v,
    int n)
/*< tridiagonal solver >*/
{
    int i;

#ifdef SF_HAS_COMPLEX_H
    b[0]/=a[0];
    v[0]/=a[0];
#else
    b[0] = sf_cdiv(b[0],a[0]);
    v[0] = sf_cdiv(v[0],a[0]);
#endif
    for(i=1;i<n;i++) {
#ifdef SF_HAS_COMPLEX_H
	a[i] -= c[i-1]*b[i-1];
	if(i<n-1) b[i]/=a[i];
	v[i]=( v[i]-c[i-1]*v[i-1] ) / a[i];
#else
	a[i] = sf_csub(a[i],sf_cmul(c[i-1],b[i-1]));
	if(i<n-1) b[i] = sf_cdiv(b[i],a[i]);
	v[i]=sf_cdiv(sf_csub(v[i],sf_cmul(c[i-1],v[i-1])),a[i]);
#endif
    }

    for(i=n-2;i>=0;i--) {
#ifdef SF_HAS_COMPLEX_H
	v[i] -= b[i]*v[i+1];
#else
	v[i] = sf_csub(v[i],sf_cmul(b[i],v[i+1]));
#endif
    }
}

static sf_complex cdprod (int j, 
			  const sf_complex* a, const sf_complex* b) 
/* complex dot product */
{
    int i;
    sf_complex c;
    c = sf_cmplx(0.,0.);
    for (i=1; i <= j; i++) {
#ifdef SF_HAS_COMPLEX_H
	c += a[j-i]*conjf(b[i]);
#else
	c = sf_cadd(c,sf_cmul(a[j-i],conjf(b[i])));
#endif
    }
    return c;
}

sf_complex plyr_val(int n, float *r, sf_complex x)
/*< polynomial value >*/
{
    int i;
    sf_complex val;

    val = sf_cmplx(r[n],0.0);
    for(i=0; i<n; i++)
    {
#ifdef SF_HAS_COMPLEX_H
        val = val*(x-r[i]);
#else
        val = sf_cmul(val,sf_cadd(x,sf_cmplx(-r[i],0.0)));
#endif
    }
    return val;
}