sf_grad_regulation_eb_3

/* A demo of how to calculate the gradient of FWI
Note: surface-free and effective boundary for reconstructing the back wave field.
*/
/*
Important references:
[1] Clayton, Robert, and Bj?rn Engquist. "Absorbing boundary
conditions for acoustic and elastic wave equations." Bulletin
of the Seismological Society of America 67.6 (1977): 1529-1540.
[2] Tarantola, Albert. "Inversion of seismic reflection data in the
acoustic approximation." Geophysics 49.8 (1984): 1259-1266.
[3] Pica, A., J. P. Diet, and A. Tarantola. "Nonlinear inversion
of seismic reflection data in a laterally invariant medium."
Geophysics 55.3 (1990): 284-292.
[4] Dussaud, E., Symes, W. W., Williamson, P., Lemaistre, L.,
Singer, P., Denel, B., & Cherrett, A. (2008). Computational
strategies for reverse-time migration. In SEG Technical Program
Expanded Abstracts 2008 (pp. 2267-2271).
[5] Hager, William W., and Hongchao Zhang. "A survey of nonlinear
conjugate gradient methods." Pacific journal of Optimization
2.1 (2006): 35-58.
*/

#include <rsf.h>
#include <time.h>

#ifdef _OPENMP
#include <omp.h>
#include <time.h>
#endif

static int nb, nz, nx, nt, ng, ns, nzpad, nxpad;
static float dz, dx, dt, fm,amp, c0, c11, c12, c21, c22;
static float *bndr, *wlt;
static float **vv, **sp0, **sp1, **gp0, **gp1, **ptr = NULL, *trans, *dobs;
static int *sxz, *gxz;

void expand2d(float** b, float** a)
/*< expand domain of 'a' to 'b': source(a)-->destination(b) >*/
{
	int iz, ix;
	for (ix = 0; ix<nx; ix++) {
		for (iz = 0; iz<nz; iz++) {
			b[nb + ix][iz] = a[ix][iz];
		}
	}

	for (ix = 0; ix<nxpad; ix++) {
		for (iz = 0; iz<nb; iz++) {
			//b[ix][iz] = b[ix][nb];
			b[ix][nzpad - iz - 1] = b[ix][nzpad - nb - 1];
		}
	}

	for (ix = 0; ix<nb; ix++) {
		for (iz = 0; iz<nzpad; iz++) {
			b[ix][iz] = b[nb][iz];
			b[nxpad - ix - 1][iz] = b[nxpad - nb - 1][iz];
		}
	}
	}


void window2d(float **a, float **b)
/*< window 'b' to 'a': source(b)-->destination(a) >*/
{
	int iz, ix;


	for (ix = 0; ix<nx; ix++) {
		for (iz = 0; iz<nz; iz++) {
			a[ix][iz] = b[nb + ix][iz];
		}
	}
}

void apply_sponge(float**p0, float **p1)
/*< apply absorbing boundary condition >*/
{
	int ix, iz;

#ifdef _OPENMP
#pragma omp parallel for			\
    private(ix,iz)				\
    shared(bndr,p0,p1)
#endif
	for (ix = 0; ix<nxpad; ix++)
	{
		//for (iz = 0; iz<nb; iz++) {	// top ABC			
		//	p0[ix][iz] = bndr[iz] * p0[ix][iz];
		//	p1[ix][iz] = bndr[iz] * p1[ix][iz];
		//}
		for (iz = nz; iz<nzpad; iz++) {// bottom ABC			
			p0[ix][iz] = bndr[nzpad - iz - 1] * p0[ix][iz];
			p1[ix][iz] = bndr[nzpad - iz - 1] * p1[ix][iz];
		}
	}

#ifdef _OPENMP
#pragma omp parallel for			\
    private(ix,iz)				\
    shared(bndr,p0,p1)
#endif
	for (iz = 0; iz<nzpad; iz++)
	{
		for (ix = 0; ix<nb; ix++) {	// left ABC			
			p0[ix][iz] = bndr[ix] * p0[ix][iz];
			p1[ix][iz] = bndr[ix] * p1[ix][iz];
		}
		for (ix = nx + nb; ix<nxpad; ix++) {// right ABC			
			p0[ix][iz] = bndr[nxpad - ix - 1] * p0[ix][iz];
			p1[ix][iz] = bndr[nxpad - ix - 1] * p1[ix][iz];
		}
	}
}


void step_forward(float **p0, float **p1)
/*< forward modeling step >*/
{
	int ix, iz;
	float tmp;

#ifdef _OPENMP
#pragma omp parallel for				\
    private(ix,iz,tmp)					\
    shared(p0,p1,vv,c0,c11,c12,c21,c22,nzpad,nxpad)
#endif	
	for (ix = 2; ix < nxpad - 2; ix++)
		for (iz = 2; iz < nzpad - 2; iz++)
		{
			tmp = c0*p1[ix][iz] +
				c11*(p1[ix][iz - 1] + p1[ix][iz + 1]) +
				c12*(p1[ix][iz - 2] + p1[ix][iz + 2]) +
				c21*(p1[ix - 1][iz] + p1[ix + 1][iz]) +
				c22*(p1[ix - 2][iz] + p1[ix + 2][iz]);
			p0[ix][iz] = 2 * p1[ix][iz] - p0[ix][iz] + vv[ix][iz] * tmp;
		}
}

void matrix_transpose(float *matrix, float *trans, int n1, int n2)
/*< matrix transpose: matrix tansposed to be trans >*/
{
	int i1, i2;

	for (i2 = 0; i2<n2; i2++)
		for (i1 = 0; i1<n1; i1++)
			trans[i2 + n2*i1] = matrix[i1 + n1*i2];
}

void sg_init(int *sxz, int szbeg, int sxbeg, int jsz, int jsx, int ns, int nz)
/*< shot/geophone position initialize >*/
{
	int is, sz, sx;
	for (is = 0; is<ns; is++)
	{
		sz = szbeg + is*jsz;
		sx = sxbeg + is*jsx;
		sxz[is] = sz + nz*sx;
	}
}

void add_source(float **p, float *source, int *sxz, int ns, int nz, int nb, bool add)
/*< add/subtract seismic sources >*/
{
	int is, sx, sz;
	if (add) {
		for (is = 0; is<ns; is++) {
			sx = sxz[is] / nz;
			sz = sxz[is] % nz;
			p[sx + nb][sz ] += source[is];
		}
	}
	else {
		for (is = 0; is<ns; is++) {
			sx = sxz[is] / nz;
			sz = sxz[is] % nz;
			p[sx + nb][sz ] -= source[is];
		}
	}
}

void record_seis(float *seis_it, int *gxz, float **p, int ng, int nz, int nb)
/*< record seismogram at time it into a vector length of ng >*/
{
	int ig, gx, gz;
	for (ig = 0; ig<ng; ig++)
	{
		gx = gxz[ig] / nz;
		gz = gxz[ig] % nz;
		seis_it[ig] = p[gx + nb][gz ];
	}
}




void step_backward(float **illum, float **lap, float **p0, float **p1)
/*< step backward >*/
{
	int ix, iz;
	float tmp;
	float  diff1, diff2;
#ifdef _OPENMP
#pragma omp parallel for				\
    private(ix,iz,tmp)					\
    shared(p0,p1,vv,c0,c11,c12,c21,c22,nzpad,nxpad)
#endif	
	for (ix = 2; ix < nxpad - 2; ix++)
		for (iz = 2; iz < nzpad - 2; iz++)
		{
			tmp = c0*p1[ix][iz] +
				c11*(p1[ix][iz - 1] + p1[ix][iz + 1]) +
				c12*(p1[ix][iz - 2] + p1[ix][iz + 2]) +
				c21*(p1[ix - 1][iz] + p1[ix + 1][iz]) +
				c22*(p1[ix - 2][iz] + p1[ix + 2][iz]);
			diff1 = diff2 = -2.0*p1[ix][iz];
			diff1 += (iz - 1 >= 0) ? p1[ix][iz - 1] : 0.0;
			diff1 += (iz + 1<nzpad) ? p1[ix][iz + 1] : 0.0;
			diff2 += (ix - 1 >= 0) ? p1[ix - 1][iz] : 0.0;
			diff2 += (ix + 1<nxpad) ? p1[ix + 1][iz] : 0.0;
			lap[ix][iz] = diff1 + diff2;
			illum[ix][iz] += p1[ix][iz] * p1[ix][iz];
			p0[ix][iz] = 2 * p1[ix][iz] - p0[ix][iz] + vv[ix][iz] * tmp;
		}
}



void rw_bndr(float *bndr, float **p, int nz, int nx,int nb,int order, bool write)
/*< if write==true, write/save boundaries out of variables;
else  read boundaries into variables (for 2nd order FD) >*/
{
	int i,ix,iz;
        int mg=order/2;
	if (write){
		// top and bottom
		for (ix = 0; ix<nx; ix++)
		{
			for (iz = 0; iz<order; iz++) {	
				//bndr[iz+order*2*ix] = p[ix+nb][iz+nb];  // top 		
				bndr[iz + order  * ix] = p[ix + nb][nz - order + iz];  // bottom 	
			}
		}
		// left and right
		for (ix = 0; ix<order; ix++)
		{
			for (iz = 0; iz<nz; iz++) {
				bndr[nx * order + iz + ix*nz] = p[ix + nb][iz];  // left
				bndr[nx  * order + iz + (ix+order)*nz] = p[nb+nx-order+ix][iz];  // right	
			}
			
		}
	}
	else{
		for (ix = 0; ix<nx; ix++)
		{
			for (iz = 0; iz<order; iz++) {
				//p[ix + nb][iz + nb]=bndr[iz + order * 2 * ix];  // top 		
				p[ix + nb][nz - order + iz] = bndr[iz + order  * ix];  // bottom 	
			}
		}
		// left and right
		for (ix = 0; ix<order; ix++)
		{
			for (iz = 0; iz<nz; iz++) {
				p[ix + nb][iz] = bndr[nx * order + iz + ix*nz];   // left
				p[nb + nx - order + ix][iz] = bndr[nx  * order + iz + (ix + order)*nz];  // right	
			}

		}
		
	}
}

void cal_residuals(float *dcal, float *dobs, float *dres, int ng)
/*< calculate residual >*/
{
	int ig;
	for (ig = 0; ig<ng; ig++){
		dres[ig] = dcal[ig] - dobs[ig];
	}
}

void cal_gradient(float **grad, float **lap, float **gp)
/*< calculate gradient >*/
{
	int ix, iz;
	for (ix = 0; ix<nx; ix++){
		for (iz = 0; iz<nz; iz++){
			grad[ix][iz] += lap[ix+nb][iz] * gp[ix+nb][iz];
		}
	}
}

void scale_gradient(float **grad, float **vv, float **illum, int nz, int nx, bool precon)
/*< scale gradient >*/
{
	int ix, iz;
	float a;
	for (ix = 1; ix<nx - 1; ix++){
		for (iz = 1; iz<nz - 1; iz++){
			a = vv[ix+nb][iz];
			if (precon) a *= sqrtf(illum[ix+nb][iz] + SF_EPS);/*precondition with residual wavefield illumination*/
			grad[ix][iz] *= 2.0 / a;
		}
	}

	for (ix = 0; ix<nx; ix++){
		grad[ix][0] = grad[ix][1];
		grad[ix][nz - 1] = grad[ix][nz - 2];
	}

	for (iz = 0; iz<nz; iz++){
		grad[0][iz] = grad[1][iz];
		grad[nx - 1][iz] = grad[nx - 2][iz];
	}
}

float cal_objective(float *dres, int ng)
/*< calculate the value of objective function >*/
{
	int i;
	float a, obj = 0;

	for (i = 0; i<ng; i++){
		a = dres[i];
		obj += a*a;
	}
	return obj;
}

float cal_beta(float **g0, float **g1, float **cg, int nz, int nx)
/*< calculate beta >*/
{
	int ix, iz;
	float a, b, c;

	a = b = c = 0;
	for (ix = 0; ix<nx; ix++){
		for (iz = 0; iz<nz; iz++){
			a += g1[ix][iz] * (g1[ix][iz] - g0[ix][iz]);// numerator of HS
			b += cg[ix][iz] * (g1[ix][iz] - g0[ix][iz]);// denominator of HS,DY
			c += g1[ix][iz] * g1[ix][iz];		// numerator of DY
		}
	}

	float	beta_HS = (fabsf(b)>0) ? (a / b) : 0.0;
	float beta_DY = (fabsf(b)>0) ? (c / b) : 0.0;
	return SF_MAX(0.0, SF_MIN(beta_HS, beta_DY));
}


void cal_conjgrad(float **g1, float **cg, float beta, int nz, int nx)
/*< calculate conjugate gradient >*/
{
	int ix, iz;

	for (ix = 0; ix<nx; ix++){
		for (iz = 0; iz<nz; iz++){
			cg[ix][iz] = -g1[ix][iz] + beta*cg[ix][iz];
		}
	}
}

float cal_epsilon(float **vv, float **cg, int nz, int nx)
/*< calculate epsilcon >*/
{
	int ix, iz;
	float vvmax, cgmax;
	vvmax = cgmax = 0.0;

	for (ix = 0; ix<nx; ix++){
		for (iz = 0; iz<nz; iz++){
			vvmax = SF_MAX(vvmax, fabsf(vv[ix][iz]));
			cgmax = SF_MAX(cgmax, fabsf(cg[ix][iz]));
		}
	}

	return 0.01*vvmax / (cgmax + SF_EPS);
}

void cal_vtmp(float **vtmp, float **vv, float **cg, float epsil, int nz, int nx)
/*< calculate temporary velcity >*/
{
	int ix, iz;

	for (ix = 0; ix<nx; ix++){
		for (iz = 0; iz<nz; iz++){
			vtmp[ix][iz] = vv[ix+nb][iz] + epsil*cg[ix][iz];
		}
	}
}

void sum_alpha12(float *alpha1, float *alpha2, float *dcaltmp, float *dobs, float *derr, int ng)
/*< calculate numerator and denominator of alpha >*/
{
	int ig;
	float a, b, c;
	for (ig = 0; ig<ng; ig++){
		c = derr[ig];
		a = dobs[ig] + c;/* since f(mk)-dobs[id]=derr[id], thus f(mk)=b+c; */
		b = dcaltmp[ig] - a;/* f(mk+epsil*cg)-f(mk) */
		alpha1[ig] -= b*c;
		alpha2[ig] += b*b;
	}
}


float cal_alpha(float *alpha1, float *alpha2, float epsil, int ng)
/*< calculate alpha >*/
{
	int ig;
	float a, b;

	a = b = 0;
	for (ig = 0; ig<ng; ig++){
		a += alpha1[ig];
		b += alpha2[ig];
	}

	return (a*epsil / (b + SF_EPS));
}

void update_vel(float **vv, float **cg, float alpha, int nz, int nx)
/*< update velcity >*/
{
	int ix, iz;

	for (ix = 0; ix<nx; ix++){
		for (iz = 0; iz<nz; iz++){
			vv[ix][iz] += alpha*cg[ix][iz];
		}
	}
}

void bell_smoothz(float **g, float **smg, int rbell, int nz, int nx)
/*< gaussian bell smoothing for z-axis >*/
{
	int ix, iz, i;
	float s;

	for (ix = 0; ix<nx; ix++)
	for (iz = 0; iz<nz; iz++)
	{
		s = 0.0;
		for (i = -rbell; i <= rbell; i++) if (iz + i >= 0 && iz + i<nz) s += expf(-(2.0*i*i) / rbell)*g[ix][iz + i];
		smg[ix][iz] = s;
	}
}

void bell_smoothx(float **g, float **smg, int rbell, int nz, int nx)
/*< gaussian bell smoothing for x-axis >*/
{
	int ix, iz, i;
	float s;

	for (ix = 0; ix<nx; ix++)
	for (iz = 0; iz<nz; iz++)
	{
		s = 0.0;
		for (i = -rbell; i <= rbell; i++) if (ix + i >= 0 && ix + i<nx) s += expf(-(2.0*i*i) / rbell)*g[ix + i][iz];
		smg[ix][iz] = s;
	}
}

void set_zero(float **g, int nn)
/*< gaussian bell smoothing for x-axis >*/
{
	int ix, iz, i;
	float s;

	for (ix = 0; ix<nx; ix++)
	for (iz = 0; iz<nn; iz++)
	{
		g[ix][iz]=0;
	}
}

int main(int argc, char *argv[])
{
	/* variables on host */
	bool verb, precon, csdgather,ifsmooth;
	int is, it, iter, niter, distx, distz, csd, rbell;
	int sxbeg, szbeg, gxbeg, gzbeg, jsx, jsz, jgx, jgz;/*  parameters of acquisition geometery */
	float tmp,  obj1, obj, beta, epsil, alpha;
	float  *dcal, *derr,  *brw,  *objval;
	float  **illum, **lap, **vtmp,  **g0, **g1, **cg, *alpha1, *alpha2 ;
	float **v0;
	int order,ib,ix,iz,nn;
	clock_t start, stop;/* timer */
	sf_file vinit, shots, grads, objs, illums;/* I/O files */

	/* initialize Madagascar */
	sf_init(argc, argv);

	/* set up I/O files */
	vinit = sf_input("in");   /* initial velocity model, unit=m/s */
	shots = sf_input("shots"); /* recorded shots from exact velocity model */
	grads = sf_output("out");  /* gradient in iterations */   
	illums = sf_output("illums");/* source illumination in iterations */
	objs = sf_output("objs");/* values of objective function in iterations */
	
	/* get parameters from velocity model and recorded shots */
	if (!sf_getbool("verb", &verb)) verb = true;/* vebosity */
	if (!sf_histint(vinit, "n1", &nz)) sf_error("no n1");/* nz */
	if (!sf_histint(vinit, "n2", &nx)) sf_error("no n2");/* nx */
	if (!sf_histfloat(vinit, "d1", &dz)) sf_error("no d1");/* dz */
	if (!sf_histfloat(vinit, "d2", &dx)) sf_error("no d2");/* dx */
	if (!sf_getbool("precon", &precon)) precon = false;/* precondition or not */
	if (!sf_getint("niter", &niter))   niter = 1;	/* number of iterations */
	if (!sf_getint("order", &order))    sf_error("no order");/* dx */	/* number of iterations */
	if (!sf_getint("rbell", &rbell))   rbell = 1;	/* number of iterations */
        if (!sf_getint("nn", &nn))   nn = 10;	/* number of iterations */
	if (!sf_getbool("ifsmooth", &ifsmooth)) ifsmooth = true;/* precondition or not */
	if (!sf_histint(shots, "n1", &nt)) sf_error("no nt");
	/* total modeling time steps */
	if (!sf_histint(shots, "n2", &ng)) sf_error("no ng");
	if (!sf_histint(shots, "nb", &nb)) sf_error("no nb");
	/* total receivers in each shot */
	if (!sf_histint(shots, "n3", &ns)) sf_error("no ns");
	/* number of shots */
	if (!sf_histfloat(shots, "d1", &dt)) sf_error("no dt");
	/* time sampling interval */
	if (!sf_histfloat(shots, "amp", &amp)) sf_error("no amp");
	/* maximum amplitude of ricker */
	if (!sf_histfloat(shots, "fm", &fm)) sf_error("no fm");
	/* dominant freq of ricker */
	if (!sf_histint(shots, "sxbeg", &sxbeg)) sf_error("no sxbeg");
	/* x-begining index of sources, starting from 0 */
	if (!sf_histint(shots, "szbeg", &szbeg)) sf_error("no szbeg");
	/* x-begining index of sources, starting from 0 */
	if (!sf_histint(shots, "gxbeg", &gxbeg)) sf_error("no gxbeg");
	/* x-begining index of receivers, starting from 0 */
	if (!sf_histint(shots, "gzbeg", &gzbeg)) sf_error("no gzbeg");
	/* x-begining index of receivers, starting from 0 */
	if (!sf_histint(shots, "jsx", &jsx)) sf_error("no jsx");
	/* source x-axis  jump interval  */
	if (!sf_histint(shots, "jsz", &jsz)) sf_error("no jsz");
	/* source z-axis jump interval  */
	if (!sf_histint(shots, "jgx", &jgx)) sf_error("no jgx");
	/* receiver x-axis jump interval  */
	if (!sf_histint(shots, "jgz", &jgz)) sf_error("no jgz");
	/* receiver z-axis jump interval  */
	if (!sf_histint(shots, "csdgather", &csd)) sf_error("csdgather or not required");
	/* default, common shot-gather; if n, record at every point*/
	nzpad = nz + nb ;
	nxpad = nx + 2 * nb;

	sf_putint(grads, "n1", nz);
	sf_putint(grads, "n2", nx);
	sf_putint(grads, "n3", niter);
	sf_putfloat(grads, "d1", dz);
	sf_putfloat(grads, "d2", dx);
	sf_putint(grads, "d3", 1);
	sf_putint(grads, "o3", 1);
	sf_putstring(grads, "label1", "Depth");
	sf_putstring(grads, "label2", "Distance");
	sf_putstring(grads, "label3", "Iteration");
	sf_putint(illums, "n1", nz);
	sf_putint(illums, "n2", nx);
	sf_putfloat(illums, "d1", dz);
	sf_putfloat(illums, "d2", dx);
	sf_putint(illums, "n3", niter);
	sf_putint(illums, "d3", 1);
	sf_putint(illums, "o3", 1);
	sf_putint(objs, "n1", niter);
	sf_putint(objs, "n2", 1);
	sf_putfloat(objs, "d1", 1);
	sf_putfloat(objs, "o1", 1);
	



	csdgather = (csd>0) ? true : false;
	tmp = 1.0 / (dz*dz);
	c11 = 4.0*tmp / 3.0;
	c12 = -tmp / 12.0;
	tmp = 1.0 / (dx*dx);
	c21 = 4.0*tmp / 3.0;
	c22 = -tmp / 12.0;
	c0 = -2.0*(c11 + c12 + c21 + c22);

	vv = sf_floatalloc2(nzpad, nxpad);/* updated velocity */
	v0 = sf_floatalloc2(nz, nx);/* updated velocity */
	vtmp = sf_floatalloc2(nz, nx);/* temporary velocity computed with epsil */
	sp0 = sf_floatalloc2(nzpad, nxpad);
	sp1 = sf_floatalloc2(nzpad, nxpad);
	gp0 = sf_floatalloc2(nzpad, nxpad);
	gp1 = sf_floatalloc2(nzpad, nxpad);
	bndr = sf_floatalloc(nb);
	g0 = sf_floatalloc2(nz, nx);/* gradient at previous step */
	g1 = sf_floatalloc2(nz, nx);/* gradient at curret step */
	cg = sf_floatalloc2(nz, nx);/* conjugate gradient */
	lap = sf_floatalloc2(nzpad, nxpad);/* laplace of the source wavefield */
	illum = sf_floatalloc2(nzpad, nxpad);/* illumination of the source wavefield */
	objval = (float*)malloc(niter*sizeof(float));/* objective/misfit function */
	wlt = (float*)malloc(nt*sizeof(float));/* ricker wavelet */
	sxz = (int*)malloc(ns*sizeof(int)); /* source positions */
	gxz = (int*)malloc(ng*sizeof(int)); /* geophone positions */
	brw = (float*)malloc(nt*(2*order* nz + order*nx)*sizeof(float));/* boundaries for wavefield reconstruction */
	trans = (float*)malloc(ng*nt*sizeof(float));/* transposed one shot */
	dobs = (float*)malloc(ng*nt*sizeof(float));/* observed seismic data */
	dcal = (float*)malloc(ng*sizeof(float));	/* calculated/synthetic seismic data */
	derr = (float*)malloc(ns*ng*nt*sizeof(float));/* residual/error between synthetic and observation */
	alpha1 = (float*)malloc(ng*sizeof(float));/* numerator of alpha, length=ng */
	alpha2 = (float*)malloc(ng*sizeof(float));/* denominator of alpha, length=ng */
	
	for (it = 0; it<nt; it++) {
		tmp = SF_PI*fm*(it*dt - 1.0 / fm); tmp *= tmp;
		wlt[it] = (1.0 - 2.0*tmp)*expf(-tmp);
	}	
	for (ib = 0; ib<nb; ib++) {
		tmp = 0.015*(nb - ib);
		bndr[ib] = expf(-tmp*tmp);
	}
	/* initialize varibles */
 FILE *fp_snap,*fp_v,*fp_snap2;
        fp_snap=fopen("nsap.bin","wb"); fp_snap2=fopen("nsap2.bin","wb");
        fp_v=fopen("vv.bin","wb");
	sf_floatread(v0[0], nz*nx, vinit);
	expand2d(vv, v0);
	for (ix = 0; ix<nxpad; ix++) {
		for (iz = 0; iz<nzpad; iz++) {
			tmp = vv[ix][iz] * dt;
			vv[ix][iz] = tmp*tmp;
		}
	}

	memset(sp0[0], 0, nzpad*nxpad *sizeof(float));
	memset(sp1[0], 0, nzpad*nxpad *sizeof(float));
	memset(gp0[0], 0, nzpad*nxpad *sizeof(float));
	memset(gp1[0], 0, nzpad*nxpad *sizeof(float));
	memset(g0[0], 0, nz*nx*sizeof(float));
	memset(g1[0], 0, nz*nx*sizeof(float));
	memset(cg[0], 0, nz*nx*sizeof(float));
	memset(lap[0], 0, nzpad*nxpad *sizeof(float));
	memset(vtmp[0], 0, nz*nx*sizeof(float));
	memset(illum[0], 0, nzpad*nxpad *sizeof(float));

	if (!(sxbeg >= 0 && szbeg >= 0 && sxbeg + (ns - 1)*jsx<nx && szbeg + (ns - 1)*jsz<nz))
	{
		sf_warning("sources exceeds the computing zone!\n"); exit(1);
	}
	sg_init(sxz, szbeg, sxbeg, jsz, jsx, ns, nz);
	distx = sxbeg - gxbeg;
	distz = szbeg - gzbeg;
	if (csdgather)	{
		if (!(gxbeg >= 0 && gzbeg >= 0 && gxbeg + (ng - 1)*jgx<nx && gzbeg + (ng - 1)*jgz<nz &&
			(sxbeg + (ns - 1)*jsx) + (ng - 1)*jgx - distx <nx && (szbeg + (ns - 1)*jsz) + (ng - 1)*jgz - distz <nz))
		{
			sf_warning("geophones exceeds the computing zone!\n"); exit(1);
		}
	}
	else{
		if (!(gxbeg >= 0 && gzbeg >= 0 && gxbeg + (ng - 1)*jgx<nx && gzbeg + (ng - 1)*jgz<nz))
		{
			sf_warning("geophones exceeds the computing zone!\n"); exit(1);
		}
	}
	sg_init(gxz, gzbeg, gxbeg, jgz, jgx, ng, nz);
	memset(brw, 0, nt*(2 * order* nz + order*nx) *sizeof(float));
	memset(dobs, 0, ng*nt*sizeof(float));
	memset(dcal, 0, ng*sizeof(float));
	memset(derr, 0, ns*ng*nt*sizeof(float));
	memset(alpha1, 0, ng*sizeof(float));
	memset(alpha2, 0, ng*sizeof(float));
	memset(dobs, 0, ng*nt*sizeof(float));
	memset(objval, 0, niter*sizeof(float));

	FILE *fp_err,*fp_g;
	fp_err = fopen("errors.bin", "wb");
     //   fp_g = fopen("g.bin", "wb");
	for (iter = 0; iter<niter; iter++)
	{
		if (verb) {
			start = clock();// record starting time
			sf_warning("iter=%d", iter);
		}
		sf_seek(shots, 0L, SEEK_SET);
		memcpy(g0[0], g1[0], nz*nx*sizeof(float));
		memset(g1[0], 0, nz*nx*sizeof(float));
		memset(illum[0], 0, nzpad*nxpad *sizeof(float));
		memset(lap[0], 0, nzpad*nxpad * sizeof(float));
		for (is = 0; is<ns; is++)
		{
			sf_floatread(trans, ng*nt, shots);
			matrix_transpose(trans, dobs, nt, ng);
			if (csdgather)	{
				gxbeg = sxbeg + is*jsx - distx;
				sg_init(gxz, gzbeg, gxbeg, jgz, jgx, ng, nz);
			}
			memset(sp0[0], 0, nzpad*nxpad*sizeof(float));
			memset(sp1[0], 0, nzpad*nxpad*sizeof(float));
			for (it = 0; it<nt; it++)
			{
				add_source(sp1, &wlt[it], &sxz[is], 1, nz,nb, true);
				step_forward(sp0, sp1);
				apply_sponge(sp0, sp1);
				ptr = sp0; sp0 = sp1; sp1 = ptr;
				record_seis(dcal, gxz, sp0, ng, nz,nb);
				cal_residuals(dcal, &dobs[it*ng], &derr[is*ng*nt + it*ng], ng);
				//fwrite(dcal, sizeof(float), ng, fp_err);	
				rw_bndr(&brw[it*(2*order*nz + order*nx)], sp0, nz, nx,nb,order, true);
				//rw_bndr(float *bndr, float **p, int nz, int nx, int nb, int order, bool write)
    if(it==50000)
                               {  window2d(v0, sp0);
                              fwrite(v0[0],sizeof(float),nx*nz,fp_snap2); }
			}
			
                     //   matrix_transpose(&derr[is*ng*nt],trans, nt, ng);
           // sf_floatwrite(&derr[is*ng*nt], ng*nt, errors);
		    fwrite(&derr[is*ng*nt], sizeof(float), nt*ng, fp_err);	

			ptr = sp0; sp0 = sp1; sp1 = ptr;
			memset(gp0[0], 0, nzpad*nxpad*sizeof(float));
			memset(gp1[0], 0, nzpad*nxpad*sizeof(float));
			for (it = nt - 1; it>-1; it--)
			{
				rw_bndr(&brw[it*(2*order*nz + order*nx)], sp1, nz, nx,nb,order, false);
				step_backward(illum, lap, sp0, sp1);
				//step_backward(float **illum, float **lap, float **p0, float **p1)
				add_source(sp1, &wlt[it], &sxz[is], 1, nz,nb, false);
				ptr = sp0; sp0 = sp1; sp1 = ptr;
 
				add_source(gp1, &derr[is*ng*nt + it*ng], gxz, ng, nz,nb, true);
				step_forward(gp0, gp1);
				apply_sponge(gp0, gp1);
				cal_gradient(g1, lap, gp1);
				  if(it==500000)
                               {  window2d(v0, sp0);
                              fwrite(v0[0],sizeof(float),nx*nz,fp_snap); }
				ptr = gp0; gp0 = gp1; gp1 = ptr;
			}

		}

		obj = cal_objective(derr, ng*nt*ns);
//if(it== 1000)

		scale_gradient(g1, vv, illum, nz, nx, precon);
                window2d(v0, illum);
		sf_floatwrite(v0[0], nz*nx, illums);
                if(ifsmooth){
		bell_smoothz(g1, v0, rbell, nz, nx);
		bell_smoothx(v0, g1, rbell, nz, nx);
		
//fwrite(g1[0], sizeof(float), nx*nz, fp_g);
 }     
        set_zero(g1,nn);
        sf_floatwrite(g1[0], nz*nx, grads);
	}
	sf_floatwrite(&obj, niter, objs);

	free(*vv); free(vv);
	free(*vtmp); free(vtmp);
	free(*sp0); free(sp0);
	free(*sp1); free(sp1);
	
	free(*gp0); free(gp0);
	free(*gp1); free(gp1);
	
	free(*g0); free(g0);
	free(*g1); free(g1);
	free(*cg); free(cg);
	free(*lap); free(lap);
	free(*illum); free(illum);
	free(objval);
	free(wlt);
	free(sxz);
	free(gxz);
	free(bndr);
	free(trans);
	free(dobs);
	free(dcal);
	free(derr);
	free(alpha1);
	free(alpha2);

	exit(0);
}