C++ sf_alloc函数代码示例

void psrtm_lop(bool adj, bool add, int nm, int nd, sf_complex *mod, sf_complex *dat)
    /*< lowrank onestep psrtm linear operator >*/
{
    sf_complex ***rec, **img;
    int nz, nx, nt, gpl, shtnum;
    int ix,is;

    nz = geop->nz;
    nx = geop->nx;
    nt = geop->nt;
    gpl = geop->gpl;
    shtnum = geop->shtnum;

    /*check the dimension*/
    if (nm!=nz*nx || nd!=nt*gpl*shtnum) sf_error("%s: wrong dimensions",__FILE__);
    sf_cadjnull(adj, add, nm, nd, mod, dat);

    /*allocate memory and transform the dimension*/
    img = (sf_complex**) sf_alloc (nx,sizeof(sf_complex*)); 
    img[0] = mod;
    for (ix=1; ix<nx; ix++) img[ix] = img[0]+ix*nz;

    rec = (sf_complex***) sf_alloc (shtnum,sizeof(sf_complex**));
    rec[0] = (sf_complex**) sf_alloc (gpl*shtnum,sizeof(sf_complex*));
    rec[0][0] = dat;
    for (ix=1; ix<gpl*shtnum; ix++) rec[0][ix] = rec[0][0]+ix*nt; 
    for (is=1; is<shtnum; is++) rec[is] = rec[0]+is*gpl;

    /* prestack migration */
    geop->adj=adj;
    psrtm(rec, img, geop);

    free(img);
    free(*rec); free(rec);
}

int main (int argc, char **argv)
{
	int np;
	char buf[BUFSZ*PARSZ], **para, *p, *sep;

	sf_init(argc, argv);
	if ((sep=sf_getstring("sep"))==NULL) sep="";
	/* cmd separater */

	para = (char **)sf_alloc(PARSZ, sizeof(char*));
	*para = (char *)sf_alloc(PARSZ*BUFSZ, sizeof(char));
	for(np=1;np<PARSZ; np++) para[np] = para[np-1] + BUFSZ;

    vp_init();

	p=buf;
	while(1)
	{
		*p=fgetc(stdin);
		if(issep(*p, sep))
		{
			p = cmdline(buf, p, para);
		}else if(*p==EOF) 
		{
			p = cmdline(buf, p, para);
			break;
		}else	p++;
	}
	free(*para);
	free(para);
    return 0;
}

void layer_init(int n_in        /* number of layers */,
		char *types     /* [n] velocity model types */, 
		float *v0_in    /* [n] reference velocity  */,
		float **g_in    /* [n] velocity gradient */,
		float **x0_in   /* [n] gradient reference point */,
		int order       /* interpolation order */,
		int nx          /* horizontal sampling */,
		float **lays, 
		float **dips, 
		float **curs    /* [n+1][nx] interfaces */) 
/*< initialize >*/
{
    int k;

    n = n_in;

    y = sf_floatalloc(n-1);
    type = types;
    lay = (sf_eno*) sf_alloc(n,sizeof(sf_eno));
    dip = (sf_eno*) sf_alloc(n,sizeof(sf_eno));
    cur = (sf_eno*) sf_alloc(n,sizeof(sf_eno));
    table = (ktable) sf_alloc(1,sizeof(*table));

    for (k=0; k < n+1; k++) {
	lay[k] = sf_eno_init (order,nx); sf_eno_set (lay[k],lays[k]);
	dip[k] = sf_eno_init (order,nx); sf_eno_set (dip[k],dips[k]);
	cur[k] = sf_eno_init (order,nx); sf_eno_set (cur[k],curs[k]);
    }    

    v0 = v0_in;
    g = g_in;
    x0 = x0_in;
}

agrid2 agrid2_init (int order      /* interpolation order */, 
		    int n1, int n2 /* data dimensions */, 
		    int nd, int max)
/*< Initialize interpolation object >*/
{
    agrid2 pnt;
    int i2;

    pnt = (agrid2) sf_alloc(1,sizeof(*pnt));
    pnt->order = order; 
    pnt->n1 = n1; 
    pnt->n2 = n2;
    pnt->nd = nd;
    pnt->ng = 2*order-2;
    if (pnt->ng > pnt->n2) sf_error("ng is too big in agrid2");
    pnt->jnt = sf_eno_init (order, pnt->ng);
    pnt->f  = sf_floatalloc2(pnt->nd,pnt->ng);
    pnt->f1 = sf_floatalloc(pnt->ng);
    pnt->ent = (agrid*) sf_alloc(n2,sizeof(agrid));
    for (i2 = 0; i2 < n2; i2++) {
	pnt->ent[i2] = agrid_init (n1, nd, max);
    }

    return pnt;
}

/*------------------------------------------------------------*/
ompfft3d sf_ompfft3a3_init(int n1_, 
			   int n2_, 
			   int n3_,
			   int ompnth_)
/*< initialize FFT on axis 3 >*/
{    
    int ompith;

    ompfft3d fft;
    fft = (ompfft3d) sf_alloc(1,sizeof(*fft));

    fft->n1 = n1_; 
    fft->n2 = n2_;
    fft->n3 = n3_;
    fft->ompnth=ompnth_;

    fft->forw = (kiss_fft_cfg*) sf_alloc(fft->ompnth,sizeof(kiss_fft_cfg));
    fft->invs = (kiss_fft_cfg*) sf_alloc(fft->ompnth,sizeof(kiss_fft_cfg));

    for(ompith=0; ompith<fft->ompnth; ompith++) {
	fft->forw[ompith] = kiss_fft_alloc(fft->n3,0,NULL,NULL);
	fft->invs[ompith] = kiss_fft_alloc(fft->n3,1,NULL,NULL);
    
	if (NULL == fft->forw[ompith] || NULL == fft->invs[ompith]) 
	    sf_error("%s: KISS FFT allocation error",__FILE__);
    
	fft->trace = (kiss_fft_cpx**) sf_complexalloc2(fft->n3,fft->ompnth);
    }

    fft->scale = 1./sqrtf(fft->n3);

    return fft;
}

nfilter nallocate(int np   /* number of patches */, 
		  int nd   /* data size */, 
		  int *nh  /* filter size [np] */, 
		  int *pch /* patching [nd] */) 
/*< allocate >*/
{
    nfilter aa;
    int ip, id;
    
    aa = (nfilter) sf_alloc(1,sizeof(*aa));

    aa->np = np;
    aa->hlx = (sf_filter*) sf_alloc(np,sizeof(sf_filter));

    for (ip=0; ip < np; ip++) {
	aa->hlx[ip] = sf_allocatehelix(nh[ip]);
    }
    
    aa->pch = sf_intalloc(nd);
    for (id=0; id < nd; id++) {
	aa->pch[id] = pch[id];
    }

    aa->mis = NULL;
    return aa;
}

eno3 eno3_init (int order              /* interpolation order */, 
		int n1, int n2, int n3 /* data dimensions */)
/*< Initialize interpolation object >*/
{
    eno3 pnt;
    int i2, i3;
    
    pnt = (eno3) sf_alloc(1,sizeof(*pnt));
    pnt->order = order; 
    pnt->n1 = n1; 
    pnt->n2 = n2;
    pnt->n3 = n3;
    pnt->ng = 2*order-2;
    if (pnt->ng > n2 || pnt->ng > n3) 
	sf_error("%s: ng=%d is too big",__FILE__,pnt->ng);
    pnt->jnt = sf_eno2_init (order, pnt->ng, pnt->ng);
    pnt->f  = sf_floatalloc2(pnt->ng,pnt->ng);
    pnt->f1 = sf_floatalloc2(pnt->ng,pnt->ng);
    pnt->ent = (sf_eno**) sf_alloc(n3,sizeof(sf_eno*));
    for (i3 = 0; i3 < n3; i3++) {
	pnt->ent[i3] = (sf_eno*) sf_alloc(n2,sizeof(sf_eno));
	for (i2 = 0; i2 < n2; i2++) {
	    pnt->ent[i3][i2] = sf_eno_init (order, n1);
	}
    }

    return pnt;
}

int fft2_init(bool cmplx1        /* if complex transform */,
              int pad1           /* padding on the first axis */,
              int nx,   int ny   /* input data size */,
              int *nx2, int *ny2 /* padded data size */)
/*< initialize >*/
{
#ifndef SF_HAS_FFTW
    int i2;
#endif

    cmplx = cmplx1;

    if (cmplx) {
        nk = n1 = kiss_fft_next_fast_size(nx*pad1);

#ifndef SF_HAS_FFTW
        cfg1  = kiss_fft_alloc(n1,0,NULL,NULL);
        icfg1 = kiss_fft_alloc(n1,1,NULL,NULL);
#endif
    } else {
        nk = kiss_fft_next_fast_size(pad1*(nx+1)/2)+1;
        n1 = 2*(nk-1);

#ifndef SF_HAS_FFTW
        cfg  = kiss_fftr_alloc(n1,0,NULL,NULL);
        icfg = kiss_fftr_alloc(n1,1,NULL,NULL);
#endif
    }

    n2 = kiss_fft_next_fast_size(ny);

    if (cmplx) {
        cc = sf_complexalloc2(n1,n2);
    } else {
        ff = sf_floatalloc2(n1,n2);
    }

#ifndef SF_HAS_FFTW
    cfg2  = kiss_fft_alloc(n2,0,NULL,NULL);
    icfg2 = kiss_fft_alloc(n2,1,NULL,NULL);

    tmp =    (kiss_fft_cpx **) sf_alloc(n2,sizeof(*tmp));
    tmp[0] = (kiss_fft_cpx *)  sf_alloc(nk*n2,sizeof(kiss_fft_cpx));
    for (i2=0; i2 < n2; i2++) {
        tmp[i2] = tmp[0]+i2*nk;
    }

    trace2 = sf_complexalloc(n2);
    ctrace2 = (kiss_fft_cpx *) trace2;
#endif

    *nx2 = n1;
    *ny2 = n2;

    wt =  1.0/(n1*n2);

    return (nk*n2);
}

sf_stack sf_stack_init (size_t size)
/*< create a stack >*/
{
    sf_stack s;

    s = (sf_stack) sf_alloc (1, sizeof(*s));
    s->size = size;
    s->entry = (struct entry*) sf_alloc(size,sizeof(struct entry));
    s->top = s->entry-1;

    return s;
}

int cfft2_init(int pad1           /* padding on the first axis */,
               int nx,   int ny   /* input data size */,
               int *nx2, int *ny2 /* padded data size */)
/*< initialize >*/
{

#ifdef SF_HAS_FFTW
#ifdef _OPENMP
    fftw_init_threads();
    sf_warning("Using threaded FFTW3! \n");
    fftw_plan_with_nthreads(omp_get_max_threads());
#endif
#endif

#ifndef SF_HAS_FFTW
    int i2;
#endif

    nk = n1 = kiss_fft_next_fast_size(nx*pad1);

#ifndef SF_HAS_FFTW
    cfg1  = kiss_fft_alloc(n1,0,NULL,NULL);
    icfg1 = kiss_fft_alloc(n1,1,NULL,NULL);
#endif

    n2 = kiss_fft_next_fast_size(ny);

    cc = sf_complexalloc2(n1,n2);
    dd = sf_complexalloc2(nk,n2);

#ifndef SF_HAS_FFTW
    cfg2  = kiss_fft_alloc(n2,0,NULL,NULL);
    icfg2 = kiss_fft_alloc(n2,1,NULL,NULL);

    tmp =    (kiss_fft_cpx **) sf_alloc(n2,sizeof(*tmp));
    tmp[0] = (kiss_fft_cpx *)  sf_alloc(nk*n2,sizeof(kiss_fft_cpx));
    for (i2=0; i2 < n2; i2++) {
        tmp[i2] = tmp[0]+i2*nk;
    }

    trace2 = sf_complexalloc(n2);
    ctrace2 = (kiss_fft_cpx *) trace2;
#endif

    *nx2 = n1;
    *ny2 = n2;

    wt =  1.0/(n1*n2);

    return (nk*n2);
}

/*------------------------------------------------------------*/
wexcub3d wex_cube(bool    verb_,
		  sf_axis amx_,
		  sf_axis amy_,
		  sf_axis az_,
		  sf_axis alx_,
		  sf_axis aly_,
		  sf_axis aw_,
		  sf_axis ae_,
		  float   eps_,
		  int     ompnth_
    )
/*< initialize SR migration space >*/
{
    wexcub3d cub;
    cub = (wexcub3d) sf_alloc(1,sizeof(*cub));

    cub->verb=verb_;
    cub->amx = sf_nod(amx_);
    cub->amy = sf_nod(amy_);
    cub->az  = sf_nod(az_);

    cub->alx = sf_nod(alx_);
    cub->aly = sf_nod(aly_);
    cub->ae  = sf_nod(ae_);

    cub->aw    = sf_nod(aw_);
    cub->aw.d *= 2.*SF_PI; /* from hertz to radians */
    cub->aw.o *= 2.*SF_PI;

    cub->eps      = eps_;

    cub->ompnth   = ompnth_;

    return cub;
}

sf_map sf_stretch_init (int n1, float o1, float d1 /* regular axis */, 
			int nd                     /* data length */, 
			float eps                  /* regularization */, 
			bool narrow                /* if zero boundary */)
/*< initialize >*/
{
    sf_map str;
    
    str = (sf_map) sf_alloc (1, sizeof(*str));

    str->nt = n1; 
    str->t0 = o1; 
    str->dt = d1; 
    str->nd = nd; 
    str->eps = eps;
    str->narrow = narrow;
    
    str->x = sf_intalloc (nd);
    str->m = sf_boolalloc (nd);
    str->w = sf_floatalloc (nd);
    str->diag = sf_floatalloc (n1);
    str->offd = sf_floatalloc (n1-1);
  
    str->slv = sf_tridiagonal_init (n1);

    return str;
}

/*------------------------------------------------------------*/
fdm2d fdutil_init(bool verb_, 
		  bool free_,
		  sf_axis a1_, 
		  sf_axis a2_, 
		  int     nb_,
		  int ompchunk_) 
/*< init fdm utilities >*/
{ 
    fdm2d fdm;
    fdm = (fdm2d) sf_alloc(1,sizeof(*fdm));

    fdm->free=free_;
    fdm->verb=verb_;

    fdm->nb=nb_;

    fdm->n1=sf_n(a1_);
    fdm->n2=sf_n(a2_);

    fdm->d1=sf_d(a1_);
    fdm->d2=sf_d(a2_);

    fdm->o1=sf_o(a1_);
    fdm->o2=sf_o(a2_);

    fdm->n1pad=sf_n(a1_)+2*fdm->nb;
    fdm->n2pad=sf_n(a2_)+2*fdm->nb;
	
    fdm->o1pad=sf_o(a1_)-fdm->nb*fdm->d1;
    fdm->o2pad=sf_o(a2_)-fdm->nb*fdm->d2;

    fdm->ompchunk=ompchunk_;

    return fdm;
}

/*------------------------------------------------------------*/
abcone2d abcone2d_make(int     nop,
		       float    dt,
		       float**  vp,
		       bool   free, 
		       fdm2d   fdm)
/*< init 2D ABC >*/
{
    abcone2d abc;
    int i1,i2;
    float d;

    abc = (abcone2d) sf_alloc(1,sizeof(*abc));

    abc->free = free;

    abc->b1l = sf_floatalloc(fdm->n2pad);
    abc->b1h = sf_floatalloc(fdm->n2pad);
    abc->b2l = sf_floatalloc(fdm->n1pad);
    abc->b2h = sf_floatalloc(fdm->n1pad);

    for (i2=0;i2<fdm->n2pad;i2++) {
	d = vp[i2][           nop  ] *dt/fdm->d1; abc->b1l[i2] = (1-d)/(1+d);
	d = vp[i2][fdm->n1pad-nop-1] *dt/fdm->d1; abc->b1h[i2] = (1-d)/(1+d);
    }
    for (i1=0;i1<fdm->n1pad;i1++) {
	d = vp[           nop  ][i1] *dt/fdm->d2; abc->b2l[i1] = (1-d)/(1+d);
	d = vp[fdm->n2pad-nop-1][i1] *dt/fdm->d2; abc->b2h[i1] = (1-d)/(1+d);
    }

    return abc;
}

/*------------------------------------------------------------*/
ompsft3d sf_ompsft3_init(int n,
			 float o,
			 float d,
			 int ompnth)
/*< origin shift (assumes no centering) >*/
{
    int ompith;

    int k,i;
    float w,s;
    
    ompsft3d sft;
    sft = (ompsft3d) sf_alloc(1,sizeof(*sft));
    
    w=2.0*SF_PI/(n*d) * o;
    k=n/2;
    
    sft->www = sf_complexalloc2(n,ompnth);
    for(ompith=0; ompith<ompnth; ompith++) {
	for(i=0; i<n; i++) { sft->www[ompith][i]=sf_cmplx(1.0,0.0); }
    }

    for(ompith=0; ompith<ompnth; ompith++) {
	for(i=0; i<k; i++) {
	    s = w * i;
	    sft->www[ompith][i]   = sf_cmplx(cosf(s),sinf(s));
	    
	    s = w * (-k+i);
	    sft->www[ompith][k+i] = sf_cmplx(cosf(s),sinf(s));
	}
    }
    
    return sft;
}