#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);
}