vmar.c

#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <math.h>
#include "par.h"
#include "segy.h"

#define MIN(x,y) ((x) < (y) ? (x) : (y))
#define MAX(x,y) ((x) > (y) ? (x) : (y))
#define NINT(x) ((int)((x)>0.0?(x)+0.5:(x)-0.5))

#ifndef PI
#define PI (3.141592653589793)
#endif

#ifdef _OPENMP
int omp_get_thread_num(void);
#endif
double wallclock_time(void);
void name_ext(char *filename, char *extension);

int compare(const void *a, const void *b) 
{ return (*(float *)b-*(float *)a); }

typedef struct { /* complex number */
        float r,i;
} complex;

void cr1fft(complex *cdata, float *rdata, int n, int sign);
int optncr(int n);

int getFileInfo(char *filename, int *n1, int *n2, int *ngath, float *d1, float *d2, float *f1, float *f2, float *xmin, float *xmax, float *sclsxgx, int *nxm);

int readDataTD(char *filename, float xmin, float dx, float *xrcv, float *xsrc, int *xnx, float *rdata, int nw, int nw_low, int ngath, int nx, int nxm, int ntfft, float alpha, float scale, int transpose, int verbose);

int readReflData(char *filename, float xmin, float dx, float *xrcv, float *xsrc, int *xnx, complex *cdata, complex *cdata2, int nw, int nw_low, int ngath, int nx, int nxm, int ntfft, float alpha, float scl, float conjg, int transpose, int verbose);

int writeData(FILE *fp, float *data, segy *hdrs, int n1, int n2);

int Marchenko_Iterations(float *inif, float *WinA, float *WinB, float *rdatavp, float *rdatavm, float *rdatagm, float *rdatagp, complex *Reflw, complex *cjReflw, float fftscl, int ntfft, int nw, int nw_low, int nblock, size_t nstationA, size_t nstationB, int niter, int squaremat, int verbose);

/*********************** self documentation **********************/
char *sdoc[] = {
" ",
" Marchenko (f and v) based on CBLAS routines (OpenMP)",
"  ",
" invmar file_inif= file_Refl= file_WinA= fileWinB= file_fm= file_fp= file_gm= file_gp= [optional parameters]",
"  ",
" Required parameters: ",
" ",
"   file_inif= ................... name of file(s) which store the initial focusing function data",
"   file_Refl= .................. name of file(s) which store the Refl data",
"   file_WinA= ............... name of file which store the Mute window A data",
" ",
" Optional parameters: ",
" ",
"   file_WinB= ............... name of file which store the Mute window B data, if empty winB=winA",
" ",
" OUTPUT FILES ",
"   file_fp= ................. output downgoing focusing function ",
"   file_fm= ................. output upgoing focusing function ",
"   file_gp= ................. output downgoing Green's function ",
"   file_gm= ................. output upgoing Green's function ",
" ",
" INPUT DEFINITION ",
"   ntc=nt ................... number of output time samples",
"   ntfft=nt ................. number of samples used in fft",
"   fmin=0 ................... minimum frequency",
"   fmax=125 ................. maximum frequency to use in deconvolution",
" ",
"   scaling=1.0 .............. 1 => pressure norm and 0 => flux norm ",
"   cjRefl=1 ................. -1 => apply complex conjugate to specific input",
"   sclinif/Win*=1 ........... apply scaling factor to inif/WinA/WinB",
"   sclRefl=2.0 / 1.0 ........ Scaling factor R (defaults: 2.0 for pressure / 1.0 for flux norm)",
"   tranposeinif/Refl/Win*=0 . 1 => apply transpose to inif/Refl/WinA/WinB",
" ITERATION PARAMETERS ",
"   niter=20 ................. Number of iterations",
"   ntap=0 ................... number of taper points matrix",
"   ftap=0 ................... percentage for tapering",
"   square=1 ................. square=0 => no. shots =/= no. receivers ",
" OUTPUT DEFINITION ",
"   verbose=0 ................ silent option; >0 displays info",
" ",
" Notes: ",
"    ntc output samples of deconvolution result",
"    nt (the number of samples read by the IO routine)",
" ",
" author  : Johno van IJsseldijk : 2020 (j.e.vanijsseldijk@tudelft.nl)",
" based on MDD scripts by Jan Thorbecke",
" ",
NULL};
/************************** end self doc *************************/



int main (int argc, char **argv)
{
	FILE    *fpin, *fpout, *fmout, *gpout, *gmout;
	int		i, j, k, ret, tracf, nshots, ntraces;
	int		size, n1, n2, ntfft, nf;
	int     verbose, fullcorr, ncorstat, err;
	int     nt, nc, ncc, ntc, nshotA, nshotB, nshotC, nshotD;
	size_t  nstationA, nstationB, nstationC, nstationD, nfreq, istation, jstation, iw, it;
	int 	ntap, nxm, ngath, nw, nw_low, nw_high, distance;
	size_t  nwrite, cdatainSize, datainSize, cdataoutSize, stationSize, is;
	float	dx, dt, fmin, fmax, df, ftap;
	float	scl;
	float   f1, f2, d1, d2, sclsxgx, xmin, xmax, alpha, wshot, wpi, wrec;
 	float   *xrcvA, *xsrcA, *xrcvB, *xsrcB, *xsrcC, *xrcvC, *xsrcD, *xrcvD;
    int     *xnx;
	
	float 	sclinif,sclRefl, cjRefl, sclWinA, sclWinB;
	int  	transposeinif, transposeRefl, transposeWinA, transposeWinB;
	float	scaling;
	int  	npad;

	float   *rdatavp, *rdatavm, *rdatagp, *rdatagm;
	double  t0, t1, t2, t3, tinit, twrite, tread, tdec;
	char	*file_inif, *file_Refl, *file_WinA, *file_WinB, *file_fp, *file_fm, *file_gp, *file_gm, filename[1024], number[128];
	int     pe=0, root_pe=0, npes=1, ipe, size_s, one_file;
	segy 	*hdrs_out;
	
	complex *Reflw, *cjReflw;
	float   *inif, *WinA, *WinB;
	float	fftscl;
	int		niter, squaremat;
	
	t0 = wallclock_time();
	initargs(argc, argv);
	requestdoc(1);
	
	if (!getparfloat("scaling", &scaling)) scaling = 1.;
	if (!getparstring("file_inif", &file_inif)) file_inif=NULL;
	assert(file_inif != NULL);
	if (!getparstring("file_Refl", &file_Refl)) file_Refl=NULL;
	assert(file_Refl != NULL);
	if (!getparstring("file_WinA", &file_WinA)) file_WinA=NULL;
	assert(file_WinA != NULL);
	if (!getparstring("file_WinB", &file_WinB)) file_WinB=NULL;
	
	if (!getparstring("file_fp", &file_fp)) file_fp=NULL;
	if (!getparstring("file_fm", &file_fm)) file_fm=NULL;
	if (!getparstring("file_gp", &file_gp)) file_gp=NULL;
	if (!getparstring("file_gm", &file_gm)) file_gm=NULL;
	
	if (!getparint("one_file", &one_file)) one_file = 1;

	if (!getparfloat("fmin", &fmin)) fmin = 0.0;

	if (!getparint("niter", &niter)) niter = 10;
	if (!getparint("square", &squaremat)) squaremat = 1;

	if (!getparfloat("sclinif", &sclinif)) sclinif = 1.;
	if (!getparint("transposeinif", &transposeinif)) transposeinif = 0;
	if (!getparint("transposeRefl", &transposeRefl)) transposeRefl = 0;
	if (!getparfloat("cjRefl", &cjRefl)) cjRefl = 1.;
	if (!getparint("transposeWinA", &transposeWinA)) transposeWinA = 0;
	if (!getparfloat("sclWinA", &sclWinA)) sclWinA = 1.;
	if (!getparint("transposeWinB", &transposeWinB)) transposeWinB = 0;
	if (!getparfloat("sclWinB", &sclWinB)) sclWinB = 1.;
	if (!getparint("npad", &npad)) npad = 0;
	if (!getparint("verbose", &verbose)) verbose = 0;

#ifdef _OPENMP
	npes = atoi(getenv("OMP_NUM_THREADS"));
	assert(npes != 0);
	if (verbose) fprintf(stderr,"Number of OpenMP thread's is %d\n", npes);
#else
   npes=1;
#endif

/* get information from input files */

	nshotA = 0;
	getFileInfo(file_inif, &n1, &n2, &nshotA, &d1, &d2, &f1, &f2, &xmin, &xmax, &sclsxgx, &nxm);
	if (!getparint("nt", &nt)) nt=n1;
	if (!getparint("ntc", &ntc)) ntc = n1;
	if (!getparfloat("dt", &dt)) dt = d1;
	if (!getparfloat("dx", &dx)) dx = d2;
	if (!getparfloat("fmax", &fmax)) fmax = 125;
	nstationA = n2;

	nshotB = 0;
	getFileInfo(file_Refl, &n1, &n2, &nshotB, &d1, &d2, &f1, &f2, &xmin, &xmax, &sclsxgx, &nxm);
	assert( n1 == nt);
	nstationB = n2;
	if (squaremat) assert( nshotA == nshotB);

	nshotC = 0;
	getFileInfo(file_WinA, &n1, &n2, &nshotC, &d1, &d2, &f1, &f2, &xmin, &xmax, &sclsxgx, &nxm);
	assert( n1 == nt);
	nstationC = n2;
	if (squaremat) assert( nshotA == nshotC);
	
	if (file_WinB != NULL) {
		nshotD = 0;
		getFileInfo(file_WinB, &n1, &n2, &nshotD, &d1, &d2, &f1, &f2, &xmin, &xmax, &sclsxgx, &nxm);
		assert( n1 == nt);
		nstationD = n2;
		if (squaremat) assert( nshotA == nshotD);
	} 
	
	
	if (!getparint("ntap", &ntap)) ntap = 0;
	if (!getparfloat("ftap", &ftap)) ftap = 0.;
	if (ntap != 0) ftap = (float)ntap / (float)nstationA;
	else if (ftap != 0) ntap = NINT(ftap*nstationA);
/*================ initializations ================*/

	tinit = 0.0;
	tread = 0.0;
	tdec  = 0.0;

    if (!getparint("ntfft", &ntfft)) ntfft = nt;
	ntfft = optncr(ntfft);
	nf    = ntfft/2+1;
	df    = 1.0/(ntfft*dt);
    nw_high  = MIN( (int)((fmax)/df), nf );
    nw_low   = MAX( (int)(fmin/df), 1 );
    nw       = nw_high - nw_low + 1;
	nfreq = MIN(nf,nw);
	if (scaling==1.) fftscl=dt*dx/((float)ntfft); //Pressure Normalized
	else if (scaling==0.) fftscl=1/((float)ntfft); //Flux Normalized
		
	if (!getparfloat("sclRefl", &sclRefl)) sclRefl = ((scaling==1.) ? 2. : 1.);
/* allocate in shared memory the in- and output data */

	cdatainSize  = nt*sizeof(float);
	cdataoutSize = 0;
	if (file_fp != NULL) cdataoutSize += nstationA*nshotA*ntfft*sizeof(float);
	if (file_fm != NULL) cdataoutSize += nstationA*nshotA*ntfft*sizeof(float);
	if (file_gp != NULL) cdataoutSize += nstationA*nshotA*ntfft*sizeof(float);
	if (file_gm != NULL) cdataoutSize += nstationA*nshotA*ntfft*sizeof(float);
	rdatavp      = (float *)malloc(nstationA*nshotA*ntfft*sizeof(float));	
	rdatavm      = (float *)malloc(nstationA*nshotA*ntfft*sizeof(float));
	rdatagp      = (float *)malloc(nstationA*nshotA*ntfft*sizeof(float));	
	rdatagm      = (float *)malloc(nstationA*nshotA*ntfft*sizeof(float));	
	Reflw        = (complex *)malloc(nstationB*nfreq*nshotB*sizeof(complex));
	cjReflw        = (complex *)malloc(nstationB*nfreq*nshotB*sizeof(complex));
	inif         = (float *)malloc(nstationA*ntfft*nshotA*sizeof(float));
	WinA         = (float *)malloc(nstationC*ntfft*nshotC*sizeof(float));
	WinB 		 = (float *)malloc(nstationC*ntfft*nshotC*sizeof(float));
//	taper        = (float *)malloc(nstationA*2*sizeof(float));
	assert(rdatavp != NULL);
	assert(rdatavm != NULL);
	assert(rdatagp != NULL);
	assert(rdatagm != NULL);
	assert(inif     != NULL);
	assert(Reflw    != NULL);
	assert(cjReflw    != NULL);
	assert(WinA != NULL);
	assert(WinB != NULL);

/* for first touch binding of allocated memory */
#pragma omp parallel for schedule(static) private(jstation) default(shared)
	for (jstation=0; jstation<nstationB; jstation++) {	
		memset(&Reflw[jstation*nfreq*nshotB],0,nfreq*nshotB*sizeof(complex));
		memset(&cjReflw[jstation*nfreq*nshotB],0,nfreq*nshotB*sizeof(complex));
	}

	  
#pragma omp parallel for schedule(static) private(jstation) default(shared)
	for (jstation=0; jstation<nstationA; jstation++) { 
		memset(&inif[jstation*ntfft*nshotA],0,ntfft*nshotA*sizeof(float));
		memset(&rdatavp[jstation*ntfft*nshotA],0,ntfft*nshotA*sizeof(float));
		memset(&rdatavm[jstation*ntfft*nshotA],0,ntfft*nshotA*sizeof(float));
		memset(&rdatagp[jstation*ntfft*nshotA],0,ntfft*nshotA*sizeof(float));
		memset(&rdatagm[jstation*ntfft*nshotA],0,ntfft*nshotA*sizeof(float));
		memset(&WinA[jstation*ntfft*nshotC],0,ntfft*nshotC*sizeof(float));
		memset(&WinB[jstation*ntfft*nshotC],0,ntfft*nshotC*sizeof(float));
	}
		
    if (verbose) {
		fprintf(stderr,"--- Input Information ---\n");
        fprintf(stderr,"  dt nt ............ : %f : %d\n", dt, nt);
        fprintf(stderr,"  dx ............... : %f\n", dx);
        fprintf(stderr,"  nshotA ........... : %d\n", nshotA );
        fprintf(stderr,"  nstationA ........ : %ld\n", nstationA );
        fprintf(stderr,"  nshotB ........... : %d\n", nshotB );
        fprintf(stderr,"  nstationB ........ : %ld\n", nstationB );
		fprintf(stderr,"  nshotC ........... : %d\n", nshotC );
		fprintf(stderr,"  nstationC ........ : %ld\n", nstationC );
		if (file_WinB != 0) {
			fprintf(stderr,"  nshotD ........... : %d\n", nshotD );
			fprintf(stderr,"  nstationD ........ : %ld\n", nstationD );
		}
		fprintf(stderr,"  Scaling .......... : %e\n", fftscl);
		fprintf(stderr,"  Refl Scaling...... : %.1f\n", sclRefl);
        fprintf(stderr,"  number t-fft ..... : %d\n", ntfft);
		fprintf(stderr,"  Input  size ...... : %ld MB\n", ((file_WinB != NULL) ? (nstationA*nshotA+nstationB*nshotB+nstationC*nshotC+nstationD*nshotD)*cdatainSize/((size_t)1024*1024) : (nstationA*nshotA+nstationB*nshotB+nstationC*nshotC)*cdatainSize/((size_t)1024*1024)));
		fprintf(stderr,"  Output size ...... : %ld MB\n", (cdataoutSize/((size_t)1024*1024)));
        if (ntap != 0) fprintf(stderr,"  taper points ..... : %d (%.0f %%)\n", ntap, ftap*100.0);
        fprintf(stderr,"  process number ... : %d\n", pe);
        fprintf(stderr,"  fmin ............. : %.3f (%d)\n", fmin, nw_low);
        fprintf(stderr,"  fmax ............. : %.3f (%d)\n", fmax, nw_high);
        fprintf(stderr,"  nfreq  ........... : %ld\n", nfreq);
        fprintf(stderr,"  niter ............ : %d\n", niter);
        fprintf(stderr,"  square ........... : %d\n", squaremat);
    }

	t1 = wallclock_time();
	tinit += t1-t0;

/* read in first nt samples, and store in data */

    xsrcA     = (float *)calloc(nshotA,sizeof(float));
    xrcvA     = (float *)calloc(nshotA*nstationA,sizeof(float));
    xnx       = (int *)calloc(nshotA,sizeof(int));
	alpha = 0.0;
    readDataTD(file_inif, xmin, dx, xrcvA, xsrcA, xnx, inif, nw, nw_low, nshotA, nstationA, nstationA, ntfft, alpha, sclinif, transposeinif, verbose);
	if (verbose >= 2) fprintf(stderr," inif data read!!! \n");

    xsrcB     = (float *)calloc(nshotB,sizeof(float));
    xrcvB     = (float *)calloc(nshotB*nstationB,sizeof(float));
	xnx 	  = (int *)calloc(nshotB,sizeof(int));
	alpha = 0.0;
    readReflData(file_Refl, xmin, dx, xrcvB, xsrcB, xnx, Reflw, cjReflw, nw, nw_low, nshotB, nstationB, nstationB, ntfft, alpha, sclRefl, cjRefl, transposeRefl, verbose);
	if (verbose >= 2) fprintf(stderr," Refl data read!!! \n");

	xsrcC    = (float *)calloc(nshotC,sizeof(float));
    xrcvC     = (float *)calloc(nshotC*nstationC,sizeof(float));
    xnx       = (int *)calloc(nshotC,sizeof(int));
	alpha = 0.0;
	readDataTD(file_WinA, xmin, dx, xrcvC, xsrcC, xnx, WinA, nw, nw_low, nshotC, nstationC, nstationC, ntfft, alpha, sclWinA, transposeWinA, verbose);
	if (verbose >= 2) fprintf(stderr," WinA data read!!! \n");
	
	if (file_WinB != NULL) {
		xsrcD    = (float *)calloc(nshotD,sizeof(float));
		xrcvD     = (float *)calloc(nshotD*nstationD,sizeof(float));
		xnx       = (int *)calloc(nshotD,sizeof(int));
		alpha = 0.0;
		readDataTD(file_WinB, xmin, dx, xrcvD, xsrcD, xnx, WinB, nw, nw_low, nshotD, nstationD, nstationD, ntfft, alpha, sclWinB, transposeWinB, verbose);
		if (verbose >= 2) fprintf(stderr," WinB data read!!! \n");
	}
	else {
	#pragma omp parallel for schedule(static) private(jstation) default(shared)
		for (jstation=0; jstation<nshotC; jstation++) { 
			memcpy(&WinB[jstation*ntfft*nstationC], &WinA[jstation*ntfft*nstationC], sizeof(float)*nstationC*ntfft);	
		}
	}
	
#pragma omp parallel for schedule(static) private(jstation) default(shared)
	for (jstation=0; jstation<nshotA; jstation++) { 
		memcpy(&rdatavp[jstation*ntfft*nstationA], &inif[jstation*ntfft*nstationA], sizeof(float)*nstationA*ntfft);		
	}
	
	t2 = wallclock_time();
	tread += t2-t1; 
	
	/* Use CGEMM or CGEMV to iterate multidimensional marchenko equation */
	Marchenko_Iterations(inif, WinA, WinB, rdatavp, rdatavm, rdatagm, rdatagp, Reflw, cjReflw, fftscl, ntfft, nfreq, nw_low, nshotB, nstationA, nstationB, niter, squaremat, verbose);
	
	fflush(stderr);
	fflush(stdout);

	t3 = wallclock_time();
	tdec += t3-t2;
	if (verbose>=1) {
		fprintf(stderr,"************* PE %d ************* \n", pe);
		fprintf(stderr,"CPU-time read data         = %.3f\n", tread);
		fprintf(stderr,"CPU-time Marchenko scheme  = %.3f\n", tdec);
	}

	free(inif);
	free(WinA);
	free(WinB);
	free(Reflw);
	free(cjReflw);	

    pe = 0;
	
	hdrs_out = (segy *)calloc(nstationA,sizeof(segy));
	assert(hdrs_out != NULL);
	
	twrite = 0.0;
	if (one_file && pe==0) {
		if (file_fp != NULL)  { 
			strcpy(filename, file_fp);
			if (verbose>=2) fprintf(stderr,"writing downgoing focusing function into file %s\n", filename);
			fpout = fopen( filename, "w+" );
			assert(fpout != NULL);
		}
		if (file_fm != NULL) {
			strcpy(filename, file_fm);
			if (verbose>=2) fprintf(stderr,"writing upgoing focusing function into file %s\n", filename);
			fmout = fopen( filename, "w+" );
			assert(fmout != NULL);
		}
		if (file_gp != NULL) {
			strcpy(filename, file_gp);
			if (verbose>=2) fprintf(stderr,"writing downgoing Green's function into file %s\n", filename);
			gpout = fopen( filename, "w+" );
			assert(gpout != NULL);
		}
		if (file_gm != NULL) {
			strcpy(filename, file_gm);
			if (verbose>=2) fprintf(stderr,"writing upgoing Green's function into file %s\n", filename);
			gmout = fopen( filename, "w+" );
			assert(gmout != NULL);
		}
	}

	for (i = 0; i < nstationA; i++) {
        hdrs_out[i].ns     = ntfft;
        hdrs_out[i].trid   = 1;
        hdrs_out[i].dt     = dt*1000000;
        hdrs_out[i].f1     = -0.5*ntfft*dt;
        hdrs_out[i].f2     = f2;
        hdrs_out[i].d1     = d1;
        hdrs_out[i].d2     = d2;
        hdrs_out[i].trwf   = nstationA*nshotA;
        hdrs_out[i].scalco = -1000;
        hdrs_out[i].gx = NINT(1000*(f2+i*d2));
        hdrs_out[i].scalel = -1000;
        hdrs_out[i].tracl = i+1;
    }

	tracf = 1;
	for (jstation=0; jstation<nshotA; jstation++) {
		t1 = wallclock_time();
		
		for (i = 0; i < nstationA; i++) {
            hdrs_out[i].fldr   = jstation+1;
            hdrs_out[i].sx     = NINT((f2+dx*jstation)*1000);;
            hdrs_out[i].offset = hdrs_out[i].sx - hdrs_out[i].gx;
            hdrs_out[i].tracf  = tracf++;

        }
		
		if (file_fp != NULL)  { 
			ret = writeData(fpout, (float *)&rdatavp[jstation*nstationA*ntfft], hdrs_out, ntfft, nstationA);
			if (ret < 0 ) verr("error on writing output file.");
		}
		if (file_fm != NULL)  { 
			ret = writeData(fmout, (float *)&rdatavm[jstation*nstationA*ntfft], hdrs_out, ntfft, nstationA);
			if (ret < 0 ) verr("error on writing output file.");
		}
		if (file_gp != NULL)  { 
			ret = writeData(gpout, (float *)&rdatagp[jstation*nstationA*ntfft], hdrs_out, ntfft, nstationA);
			if (ret < 0 ) verr("error on writing output file.");
		}
		if (file_gm != NULL)  { 
			ret = writeData(gmout, (float *)&rdatagm[jstation*nstationA*ntfft], hdrs_out, ntfft, nstationA);
			if (ret < 0 ) verr("error on writing output file.");
		}
		
		t2 = wallclock_time();
		twrite += t2-t1;
	}
//****************************************************************************************************************************************************************************************************************************************************************************************************************//

	if (one_file) {
		if (file_fp != NULL) {fclose(fpout);}
		if (file_fm != NULL) {fclose(fmout);}
		if (file_gp != NULL) {fclose(gpout);}
		if (file_gm != NULL) {fclose(gmout);}
	}
	
	free(rdatavp);
	free(rdatavm);
	free(rdatagp);
	free(rdatagm);
	
	free(hdrs_out);
	
/*================ end ================*/

	if (verbose) {
		t3 = wallclock_time();
		fprintf(stderr,"CPU-time write data        = %.3f\n", twrite);
		fprintf(stderr,"CPU-time initialization    = %.3f\n", tinit);
		fprintf(stderr,"Total CPU-time             = %.3f\n", t3-t0);
	}

	exit(0);
}