added MDD application

MDD/Makefile

+# Makefile
+
+include ../Make_include
+
+########################################################################
+# define general include and system library
+ALLINC  = -I.
+LIBS += -mkl -L$L -lgenfft $(LIBSM)
+CFLAGS += -I$(MKLROOT)/include  
+
+#LIBS += -lblas -llapack -L$L -lgenfft $(LIBSM) -lc -lm
+
+all: mdd 
+
+PRG = mdd
+
+SRCC =  $(PRG).c \
+	atopkge.c \
+	docpkge.c \
+	getpars.c \
+	readShotData.c \
+	writeEigen.c \
+	deconvolve.c \
+	computeMatrixInverse.c \
+	getFileInfo.c \
+	verbosepkg.c \
+	name_ext.c \
+	wallclock_time.c
+
+OBJC	= $(SRCC:%.c=%.o)
+
+$(PRG):	$(OBJC) 
+	$(CC) $(LDFLAGS) $(CFLAGS) $(OPTC) -o $(PRG) $(OBJC) $(LIBS)
+
+install: $(PRG) 
+	cp $(PRG) $B
+
+clean:
+		rm -f core $(OBJC) $(OBJM) $(PRG) 
+
+realclean:
+		rm -f core $(OBJC) $(OBJM) $(PRG) $B/$(PRG) 
+
+
+print:	Makefile $(SRC)
+	$(PRINT) $?
+	@touch print
+
+count:
+	@wc $(SRC)
+
+tar:
+	@tar cf $(PRG).tar Makefile $(SRC) && compress $(PRG).tar
+
+
+

MDD/atopkge.c

+../utils/atopkge.c
\ No newline at end of file

MDD/computeMatrixInverse.c

+#include<math.h>
+#include<stdlib.h>
+#include<stdio.h>
+#include <assert.h>
+
+#define MAX(x,y) ((x) > (y) ? (x) : (y))
+
+/* Cholesky based inverse */
+void cpotrf_(char *uplo, int *N, float *A, int *lda, int *info);
+void cpotri_(char *uplo, int *N, float *A, int *lda, int *info);
+				
+/* LU based inverse */
+void cgetrf_(int *M, int *N, float *A, int *lda, int *ipvt, int *info);
+void cgetri_(int *N, float *A, int *lda, int *ipvt, float *work, int *lwork, int *info);
+void zgetrf_(int *M, int *N, double *A, int *lda, int *ipvt, int *info);
+void zgetri_(int *N, double *A, int *lda, int *ipvt, double *work, int *lwork, int *info);
+int ilaenv_(int *ispec, char *name, char *opts, int *n1, int *n2, int *n3, int *n4);
+
+/* SVD based inverse */
+void cgesvd_(char *jobu, char *jobvt, int *M, int *N, float *A, int *lda, float *S, float *U, int *ldu, float *vt, int *ldvt, float *work, int *lwork, float *rwork, int *info);
+void zgesvd_(char *jobu, char *jobvt, int *M, int *N, double *A, int *lda, double *S, double *U, int *ldu, double *vt, int *ldvt, double *work, int *lwork, double *rwork, int *info);
+void cgesdd_(char *jobz, int *M, int *N, float *A, int *lda, float *S, float *U, int *ldu, float *vt, int *ldvt, float *work, int *lwork, float *rwork, int *iwork, int *info);
+
+/* Eigenvalues */
+void zgeev_(char *jobvl, char *jobvr, int *N, double *A, int *lda, double *S, double *vl, int *ldvl, double *vr, int *ldvr, 
+			double *work, int *lwork, double *rwork, int *info);
+
+typedef struct { /* complex number */
+	float r,i;
+} complex;
+
+void computeMatrixInverse(complex *matrix, int nxm, int rthm, float eps_a, float eps_r, float numacc, int eigenvalues, float *eigen, int iw, int verbose)
+{
+	int i,j,k,N,lda,info,lwork,*ipvt;
+	float energy;
+	complex tmp, one, *work;
+	char *uplo;
+
+	uplo = "U";
+	lda = N = nxm;
+	one.r = 1.0;
+	one.i = 0.0;
+
+	if (rthm==0) {
+		energy=0.0;
+		if (eps_r != 0.0) {
+			for (i=0; i<nxm; i++) {
+				for (j=0; j<nxm; j++) {
+					tmp = matrix[i*nxm+j];
+					energy += sqrt(tmp.r*tmp.r+tmp.i*tmp.i);
+				}
+//		fprintf(stderr,"i=%d energy=%e\n", i, energy);
+			}
+		}
+		if (verbose>1) fprintf(stderr,"energy=%e eps_r=%e eps_a=%e\n", energy, eps_r*energy, eps_a);
+		/* add small value at diagonal */
+#pragma ivdep
+		for (i=0; i<nxm; i++) {
+			tmp.r = eps_r*energy+eps_a;
+			matrix[i*nxm+i].r+=tmp.r;
+		}
+		/* Cholesky based matrix inversion */
+		cpotrf_(uplo, &N, &matrix[0].r, &lda, &info);
+		assert (info == 0);
+		cpotri_(uplo, &N, &matrix[0].r, &lda, &info);
+		assert (info == 0);
+		/* fill lower part of inverse matrix */
+		for (i=0; i<nxm; i++) {
+#pragma ivdep
+			for (j=i+1; j<nxm; j++) {
+					matrix[i*nxm+j].r=matrix[j*nxm+i].r;
+					matrix[i*nxm+j].i=-1.0*matrix[j*nxm+i].i;
+			}
+		}
+
+	}
+	else if (rthm==1) {
+		int ispec, n1, nb;
+		char *name , *opts;
+
+		ispec = 1;
+		name = "CGETRI";
+		n1 = nxm;
+		nb = ilaenv_(&ispec, name, opts, &n1, &n1, &n1, &n1);
+		nb = MAX(1,nb);
+		lwork = nb*nxm;
+		ipvt = (int *)malloc(nxm*sizeof(int));
+		work = (complex *)malloc(lwork*sizeof(complex));
+
+		energy=0.0;
+		if (eps_r != 0.0) {
+			for (i=0; i<nxm; i++) {
+				for (j=0; j<nxm; j++) {
+					tmp = matrix[i*nxm+j];
+					energy += sqrt(tmp.r*tmp.r+tmp.i*tmp.i);
+				}
+			}
+		}
+		if (verbose>1) fprintf(stderr,"eps_r=%e eps_a=%e\n", eps_r*energy, eps_a);
+		/* add small value at diagonal */
+		for (i=0; i<nxm; i++) {
+			tmp.r = eps_r*energy+eps_a;
+			matrix[i*nxm+i].r+=tmp.r;
+		}
+		/* LU based matrix inversion */
+		cgetrf_(&nxm, &nxm, &matrix[0].r, &nxm, ipvt, &info);
+		assert (info == 0);
+		cgetri_(&nxm, &matrix[0].r, &nxm, ipvt, &work[0].r, &lwork, &info);
+		assert (info == 0);
+
+		free(ipvt);
+		free(work);
+	}
+	else if (rthm==2) { /* SVD general algorithm most accurate */
+		float *rwork, *S;
+		double S0,Si;
+		complex *U, *VT, a, b;
+		char *jobu, *jobvt;
+		int neig;
+
+		energy=0.0;
+		if (eps_r != 0.0) {
+			for (i=0; i<nxm; i++) {
+				for (j=0; j<nxm; j++) {
+					tmp = matrix[i*nxm+j];
+					energy += sqrt(tmp.r*tmp.r+tmp.i*tmp.i);
+				}
+			}
+			fprintf(stderr,"energy = %e\n", energy);
+		}
+		if (verbose>1) fprintf(stderr,"eps_r=%e eps_a=%e\n", eps_r*energy, eps_a);
+		/* add small value at diagonal */
+		for (i=0; i<nxm; i++) {
+			tmp.r = eps_r*energy+eps_a;
+			matrix[i*nxm+i].r+=tmp.r;
+		}
+
+		jobu = "A";
+		jobvt = "A";
+		lda = N = nxm;
+		lwork = N*8;
+		S = (float *)malloc(N*sizeof(float));
+		U = (complex *)malloc(N*N*sizeof(complex));
+		VT = (complex *)malloc(N*N*sizeof(complex));
+		work = (complex *)malloc(lwork*sizeof(complex));
+		rwork = (float *)malloc(5*N*sizeof(float));
+
+		/* Compute SVD */
+		cgesvd_(jobu, jobvt, &N, &N, &matrix[0].r, &lda, S, &U[0].r, &lda, &VT[0].r, 
+			&lda, &work[0].r, &lwork, rwork, &info);
+		assert (info == 0);
+
+		if (eigenvalues) {
+			for (i=0; i<N; i++) {
+					eigen[i] = S[i];
+			}
+		}
+
+		/* Compute inverse */
+		S0 = S[0];
+		neig = 0;
+		for (i=0; i<N; i++) {
+/*			fprintf(stderr,"S[%d] = %e ",i,S[i]);*/
+			Si = S[i];
+			if ((Si/S0) > numacc) { S[i]=1.0/S[i]; neig++; }
+			else S[i] = 0.0;
+			/*S[i]=1.0/(S[i]+eps_r*S[0]);*/
+/*			fprintf(stderr,"S^-1[%d] = %e\n",i,S[i]);*/
+		}
+		if(verbose) fprintf(stderr,"fraction of eigenvalues used = %.3f\n",(float)(neig/((float)N)));
+
+		for (j=0; j<N; j++) {
+			for (i=0; i<N; i++) {
+				U[j*N+i].r=S[j]*U[j*N+i].r;
+				U[j*N+i].i=-1.0*S[j]*U[j*N+i].i;
+			}
+		}
+		for (j=0; j<N; j++) {
+			for (i=0; i<N; i++) {
+				tmp.r = tmp.i = 0.0;
+				for (k=0; k<N; k++) {
+					a = U[k*N+j];
+					b.r = VT[i*N+k].r;
+					b.i = -1.0*VT[i*N+k].i;
+					tmp.r += (a.r*b.r-a.i*b.i);
+					tmp.i += (a.r*b.i+a.i*b.r);
+				}
+				matrix[j*nxm+i] = tmp;
+			}
+		}
+
+		free(U);
+		free(VT);
+		free(S);
+		free(work);
+		free(rwork);
+	}
+	else if (rthm==3) { /* SVD algorithm Divide and Conquerer less accurate */
+		/* CGESDD*/
+		int *iwork;
+		int neig;
+		float *rwork, *S;
+		double S0,Si;
+		complex *U, *VT, a, b;
+		char *jobz;
+
+		energy=0.0;
+		if (eps_r != 0.0) {
+			for (i=0; i<nxm; i++) {
+				for (j=0; j<nxm; j++) {
+					tmp = matrix[i*nxm+j];
+					energy += sqrt(tmp.r*tmp.r+tmp.i*tmp.i);
+				}
+			}
+		}
+		if (verbose>1) fprintf(stderr,"eps_r=%e eps_a=%e\n", eps_r*energy, eps_a);
+		/* add small value at diagonal */
+		for (i=0; i<nxm; i++) {
+			tmp.r = eps_r*energy+eps_a;
+			matrix[i*nxm+i].r+=tmp.r;
+		}
+
+		jobz = "A";
+		lda = N = nxm;
+		lwork = N*N+4*N;
+		S = (float *)malloc(N*sizeof(float));
+		U = (complex *)malloc(N*N*sizeof(complex));
+		VT = (complex *)malloc(N*N*sizeof(complex));
+		work = (complex *)malloc(lwork*sizeof(complex));
+		rwork = (float *)malloc(5*(N*N+N)*sizeof(float));
+		iwork = (int *)malloc(8*N*sizeof(int));
+
+		/* Compute SVD */
+		cgesdd_(jobz, &N, &N, &matrix[0].r, &lda, S, &U[0].r, &lda, &VT[0].r, 
+			&lda, &work[0].r, &lwork, rwork, iwork, &info);
+		assert (info == 0);
+
+		if (eigenvalues) {
+			for (i=0; i<N; i++) {
+					eigen[i] = S[i];
+			}
+		}
+
+		/* Compute inverse */
+		S0 = S[0];
+		neig = 0;
+		for (i=0; i<N; i++) {
+/*			fprintf(stderr,"S[%d] = %e S0 = %e\n ",i,S[i], S0);*/
+			Si = S[i];
+			if ((Si/S0) > numacc) { S[i]=1.0/S[i]; neig++; }
+			else S[i] = 0.0;
+/*			fprintf(stderr,"S^-1[%d] = %e\n",i,S[i]);*/
+		}
+		if(verbose) fprintf(stderr,"fraction of eigenvalues used = %.3f\n",(float)(neig/((float)N)));
+
+		for (j=0; j<N; j++) {
+			for (i=0; i<N; i++) {
+				U[j*N+i].r=S[j]*U[j*N+i].r;
+				U[j*N+i].i=-1.0*S[j]*U[j*N+i].i;
+			}
+		}
+		for (j=0; j<N; j++) {
+			for (i=0; i<N; i++) {
+				tmp.r = tmp.i = 0.0;
+				for (k=0; k<N; k++) {
+					a = U[k*N+j];
+					b.r = VT[i*N+k].r;
+					b.i = -1.0*VT[i*N+k].i;
+					tmp.r += (a.r*b.r-a.i*b.i);
+					tmp.i += (a.r*b.i+a.i*b.r);
+				}
+				matrix[j*nxm+i] = tmp;
+			}
+		}
+
+		free(U);
+		free(VT);
+		free(S);
+		free(work);
+		free(rwork);
+		free(iwork);
+	}
+	else if (rthm==4) { /* SVD general algorithm double precission most accurate */
+		double *rwork, *S, *U, *VT, ar, ai, br, bi, tmpr, tmpi;
+		double S0,Si,*Mat,*dwork;
+		int neig;
+		char *jobu, *jobvt;
+
+		energy=0.0;
+		if (eps_r != 0.0) {
+			for (i=0; i<nxm; i++) {
+				for (j=0; j<nxm; j++) {
+					tmp = matrix[i*nxm+j];
+					energy += sqrt(tmp.r*tmp.r+tmp.i*tmp.i);
+				}
+			}
+		}
+		if (verbose>1) fprintf(stderr,"eps_r=%e eps_a=%e\n", eps_r*energy, eps_a);
+		/* add small value at diagonal */
+		for (i=0; i<nxm; i++) {
+			tmp.r = eps_r*energy+eps_a;
+			matrix[i*nxm+i].r+=tmp.r;
+		}
+
+		Mat = (double *)malloc(2*N*N*sizeof(double));
+		/* convert to doubles */
+		for (i=0; i<nxm; i++) {
+			for (j=0; j<nxm; j++) {
+				Mat[i*2*nxm+j*2] = (double)matrix[i*nxm+j].r;
+				Mat[i*2*nxm+j*2+1] = (double)matrix[i*nxm+j].i;
+			}
+		}
+		jobu = "A";
+		jobvt = "A";
+		lda = N = nxm;
+		lwork = N*8;
+		S = (double *)malloc(N*sizeof(double));
+		U = (double *)malloc(2*N*N*sizeof(double));
+		VT = (double *)malloc(2*N*N*sizeof(double));
+		dwork = (double *)malloc(2*lwork*sizeof(double));
+		rwork = (double *)malloc(5*N*sizeof(double));
+
+		/* Compute SVD */
+		zgesvd_(jobu, jobvt, &N, &N, &Mat[0], &lda, S, &U[0], &lda, &VT[0], 
+			&lda, &dwork[0], &lwork, rwork, &info);
+		assert (info == 0);
+
+		if (eigenvalues) {
+			for (i=0; i<N; i++) {
+					eigen[i] = (float)S[i];
+			}
+		}
+
+		/* Compute inverse */
+		S0 = S[0];
+		neig = 0;
+		for (i=0; i<N; i++) {
+			if (verbose=4) fprintf(stderr,"S[%d] = %e ",i,S[i]);
+			Si = S[i];
+			if ((Si/S0) > numacc) { S[i]=1.0/S[i]; neig++; }
+			else S[i] = 0.0;
+			/*S[i]=1.0/(S[i]+eps_r*S[0]);*/
+/*			fprintf(stderr,"S^-1[%d] = %e\n",i,S[i]);*/
+		}
+		if(verbose) fprintf(stderr,"fraction of eigenvalues used = %.3f\n",(float)(neig/((float)N)));
+
+		for (j=0; j<N; j++) {
+			for (i=0; i<N; i++) {
+				U[j*2*N+2*i]=S[j]*U[j*2*N+2*i];
+				U[j*2*N+2*i+1]=-1.0*S[j]*U[j*2*N+2*i+1];
+			}
+		}
+		for (j=0; j<N; j++) {
+			for (i=0; i<N; i++) {
+				tmpr = tmpi = 0.0;
+				for (k=0; k<N; k++) {
+					ar = U[k*2*N+2*j];
+					ai = U[k*2*N+2*j+1];
+					br = VT[i*2*N+2*k];
+					bi = -1.0*VT[i*2*N+2*k+1];
+					tmpr += (ar*br-ai*bi);
+					tmpi += (ar*bi+ai*br);
+				}
+				matrix[j*nxm+i].r = (float)tmpr;
+				matrix[j*nxm+i].i = (float)tmpi;
+			}
+		}
+
+		free(U);
+		free(VT);
+		free(S);
+		free(dwork);
+		free(rwork);
+		free(Mat);
+	}
+	else if (rthm==5) { /* double precission LU decomposition */
+		int ispec, n1, nb;
+		char *name , *opts;
+		double *Mat, *dwork;
+
+		ispec = 1;
+		name = "ZGETRI";
+		n1 = nxm;
+		nb = ilaenv_(&ispec, name, opts, &n1, &n1, &n1, &n1);
+		nb = MAX(1,nb);
+		lwork = nb*nxm;
+		ipvt = (int *)malloc(nxm*sizeof(int));
+		dwork = (double *)malloc(2*lwork*sizeof(double));
+		Mat = (double *)malloc(2*N*N*sizeof(double));
+
+		energy=0.0;
+		if (eps_r != 0.0) {
+			for (i=0; i<nxm; i++) {
+				for (j=0; j<nxm; j++) {
+					tmp = matrix[i*nxm+j];
+					energy += sqrt(tmp.r*tmp.r+tmp.i*tmp.i);
+				}
+			}
+		}
+		if (verbose>1) fprintf(stderr,"eps_r=%e eps_a=%e\n", eps_r*energy, eps_a);
+		/* convert to doubles */
+		for (i=0; i<nxm; i++) {
+			for (j=0; j<nxm; j++) {
+				Mat[i*2*nxm+j*2] = (double)matrix[i*nxm+j].r;
+				Mat[i*2*nxm+j*2+1] = (double)matrix[i*nxm+j].i;
+			}
+		}
+
+		/* add small value at diagonal */
+		for (i=0; i<nxm; i++) {
+			Mat[i*2*nxm+i*2]  +=eps_r*energy+eps_a;
+//			Mat[i*2*nxm+i*2+1]+=eps_r*energy+eps_a;
+		}
+
+		/* LU based matrix inversion */
+		zgetrf_(&nxm, &nxm, &Mat[0], &nxm, ipvt, &info);
+		if (info != 0) fprintf(stderr,"error in zgetrf %d at frequency %d\n", info, iw);
+		assert (info == 0);
+		zgetri_(&nxm, &Mat[0], &nxm, ipvt, &dwork[0], &lwork, &info);
+		if (info != 0) fprintf(stderr,"error in zgetri %d at frequency %d\n", info, iw);
+		assert (info == 0);
+
+		/* convert back to floats */
+		for (i=0; i<nxm; i++) {
+			for (j=0; j<nxm; j++) {
+				matrix[i*nxm+j].r = (float)Mat[i*2*nxm+j*2];
+				matrix[i*nxm+j].i = (float)Mat[i*2*nxm+j*2+1];
+			}
+		}
+
+		free(ipvt);
+		free(dwork);
+		free(Mat);
+	}
+	else if (rthm==6) { /* eigenvalue decomposition */
+		int *iwork;
+		int neig;
+		double *work, *vr, *vl;
+		double *rwork, *S, *U, *VT, ar, ai, br, bi, tmpr, tmpi;
+		double S0,Si,nxi,*Mat;
+		char *jobvl, *jobvr;
+
+		jobvl = "V";
+		jobvr = "V";
+		lwork = N*N+2*N;
+		work  = (double *)malloc(2*lwork*sizeof(double));
+		rwork = (double *)malloc(N*2*sizeof(double));
+		vr    = (double *)malloc(2*N*N*sizeof(double));
+		vl    = (double *)malloc(2*N*N*sizeof(double));
+		S = (double *)malloc(2*N*sizeof(double));
+		U = (double *)malloc(2*N*N*sizeof(double));
+
+		Mat = (double *)malloc(2*N*N*sizeof(double));
+		/* convert to doubles */
+		for (i=0; i<nxm; i++) {
+			for (j=0; j<nxm; j++) {
+				Mat[i*2*nxm+j*2] = (double)matrix[i*nxm+j].r;
+				Mat[i*2*nxm+j*2+1] = (double)matrix[i*nxm+j].i;
+			}
+		}
+
+		zgeev_(jobvl, jobvr, &N, Mat, &N, S, vl, &N, vr, &N, 
+			work, &lwork, rwork, &info);
+		assert (info == 0);
+
+		nxi = 1.0/N;
+		for (i=0; i<N; i++) {
+			S[2*i] = (float)S[2*i]*nxi;
+			S[2*i+1] = (float)S[2*i+1]*nxi;
+		}
+		
+		for (i=0; i<N; i++) {
+			for (j=0; j<N; j++) {
+				U[i*2*N+2*j]  = (float)vr[(j)*2*N+2*i];
+				U[i*2*N+2*j+1]  = (float)vr[(i)*2*N+2*j+1];
+			}
+		}
+
+		/* Compute inverse */
+		S0 = S[0];
+		neig = 0;
+		for (i=0; i<N; i++) {
+/*			fprintf(stderr,"S[%d] = %e ",i,S[i]);*/
+			Si = S[i];
+			if ((Si/S0) > numacc) { S[i]=1.0/S[i]; neig++; }
+			else S[i] = 0.0;
+/*			fprintf(stderr,"S^-1[%d] = %e\n",i,S[i]);*/
+		}
+		if(verbose) fprintf(stderr,"fraction of eigenvalues used = %.3f\n",(float)(neig/((float)N)));
+
+		for (j=0; j<N; j++) {
+			for (i=0; i<N; i++) {
+				U[j*2*N+2*i]=S[j]*U[j*2*N+2*i];
+				U[j*2*N+2*i+1]=-1.0*S[j]*U[j*2*N+2*i+1];
+			}
+		}
+		for (j=0; j<N; j++) {
+			for (i=0; i<N; i++) {
+				tmpr = tmpi = 0.0;
+				for (k=0; k<N; k++) {
+					ar = U[k*2*N+2*j];
+					ai = U[k*2*N+2*j+1];
+					br = U[i*2*N+2*k];
+					bi = U[i*2*N+2*k+1];
+					tmpr += (ar*br-ai*bi);
+					tmpi += (ar*bi+ai*br);
+				}
+				matrix[j*nxm+i].r = (float)tmpr;
+				matrix[j*nxm+i].i = (float)tmpi;
+			}
+		}
+
+
+		free(work);
+		free(rwork);
+		free(vr);
+		free(Mat);
+		free(S);
+		free(U);
+	}
+
+	return;
+}
+

MDD/deconvolve.c

+#include <stdlib.h>
+#include <stdio.h>
+#include <assert.h>
+#include <math.h>
+#include <string.h>
+#include<mkl_cblas.h>
+
+typedef struct { /* complex number */
+	float r,i;
+} complex;
+
+/*
+cblas interface
+void cgemm(const char *transa, const char *transb, const MKL_INT *m, const MKL_INT *n, const MKL_INT *k,
+           const MKL_Complex8 *alpha, const MKL_Complex8 *a, const MKL_INT *lda,
+           const MKL_Complex8 *b, const MKL_INT *ldb, const MKL_Complex8 *beta,
+           MKL_Complex8 *c, const MKL_INT *ldc);
+*/
+
+void cgemm_(char *transA, char *transb, int *M, int *N, int *K, float *alpha, float *A, int *lda, float *B, int *ldb, float *beta, float *C, int *ldc);
+/*
+CGEMM - perform one of the matrix-matrix operations C := alpha*op( A )*op( B ) + beta*C,
+
+Synopsis
+
+SUBROUTINE CGEMM ( TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC )
+
+CHARACTER*1 TRANSA, TRANSB
+
+INTEGER M, N, K, LDA, LDB, LDC
+
+COMPLEX ALPHA, BETA
+
+COMPLEX A( LDA, * ), B( LDB, * ), C( LDC, * )
+
+TRANSA - CHARACTER*1. On entry, TRANSA specifies the form of op( A ) to be used in the matrix multiplication as follows:
+
+TRANSA = 'N' or 'n', op( A ) = A.
+
+TRANSA = 'T' or 't', op( A ) = A'.
+
+TRANSA = 'C' or 'c', op( A ) = conjg( A' ).
+
+Unchanged on exit.
+
+TRANSB - CHARACTER*1. On entry, TRANSB specifies the form of op( B ) to be used in the matrix multiplication as follows:
+
+TRANSB = 'N' or 'n', op( B ) = B.
+
+TRANSB = 'T' or 't', op( B ) = B'.
+
+TRANSB = 'C' or 'c', op( B ) = conjg( B' ).
+
+Unchanged on exit.
+
+M - INTEGER.
+On entry, M specifies the number of rows of the matrix op( A ) and of the matrix C. M must be at least zero. Unchanged on exit.
+
+N - INTEGER.
+On entry, N specifies the number of columns of the matrix op( B ) and the number of columns of the matrix C. N must be at least zero. Unchanged on exit.
+
+K - INTEGER.
+On entry, K specifies the number of columns of the matrix op( A ) and the number of rows of the matrix op( B ). K must be at least zero. Unchanged on exit.
+
+ALPHA - COMPLEX .
+On entry, ALPHA specifies the scalar alpha. Unchanged on exit.
+
+A - COMPLEX array of DIMENSION ( LDA, ka ), where ka is k when TRANSA = 'N' or 'n', and is m otherwise. Before entry with TRANSA = 'N' or 'n', the leading m by k part of the array A must contain the matrix A, otherwise the leading k by m part of the array A must contain the matrix A. Unchanged on exit.
+
+LDA - INTEGER.
+On entry, LDA specifies the first dimension of A as declared in the calling (sub) program. When TRANSA = 'N' or 'n' then LDA must be at least max( 1, m ), otherwise LDA must be at least max( 1, k ). Unchanged on exit.
+
+B - COMPLEX array of DIMENSION ( LDB, kb ), where kb is n when TRANSB = 'N' or 'n', and is k otherwise. Before entry with TRANSB = 'N' or 'n', the leading k by n part of the array B must contain the matrix B, otherwise the leading n by k part of the array B must contain the matrix B. Unchanged on exit.
+
+LDB - INTEGER.
+On entry, LDB specifies the first dimension of B as declared in the calling (sub) program. When TRANSB = 'N' or 'n' then LDB must be at least max( 1, k ), otherwise LDB must be at least max( 1, n ). Unchanged on exit.
+
+BETA - COMPLEX .
+On entry, BETA specifies the scalar beta. When BETA is supplied as zero then C need not be set on input. Unchanged on exit.
+
+C - COMPLEX array of DIMENSION ( LDC, n ).
+Before entry, the leading m by n part of the array C must contain the matrix C, except when beta is zero, in which case C need not be set on entry. On exit, the array C is overwritten by the m by n matrix ( alpha*op( A )*op( B ) + beta*C ).
+
+LDC - INTEGER.
+On entry, LDC specifies the first dimension of C as declared in the calling (sub) program. LDC must be at least max( 1, m ).  Unchanged on exit.
+
+*/
+
+void computeMatrixInverse(complex *matrix, int nxm, int rthm, float eps_a, float eps_r, float numacc, int eigenvalues, float *eigen, int iw, int verbose);
+
+int deconvolve(complex *cA, complex *cB, complex *cC, complex *oBB, int nfreq, int nblock, size_t nstationA, size_t nstationB, float eps_a, float eps_r, float numacc, int eigenvalues, float *eigen, int rthm, int mdd, int conjgA, int conjgB, int verbose)
+{
+	int istation, jstation, i, j, k, icc, ibb, NA, NB, NC, nshots;
+	size_t  iwnA, iw, iwnB, iwAB, iwBB;
+	complex *AB, *BB;
+	char *transa, *transb,*transN;
+	complex beta, alpha, tmp, a, b;
+
+	AB = (complex *)calloc(nstationA*nstationB,sizeof(complex));
+	BB = (complex *)calloc(nstationB*nstationB,sizeof(complex));
+
+	if (conjgA == 1) transa = "C";
+	else if (conjgA == 0) transa = "N";
+	else transa = "T";
+	if (conjgB == 1) transb = "C";
+	else if(conjgB ==0) transb = "N";
+	else transb = "T";
+	transN = "N";
+	alpha.r = 1.0; alpha.i = 0.0;
+	beta.r = 0.0; beta.i = 0.0;
+	nshots = nblock;
+	NA = nstationA;
+	NB = nstationB;
+	if (conjgA) NC = nshots;
+	else NC = nstationB;
+
+//	if (verbose) fprintf(stderr,"transa=%s transb=%s %d %d %d\n", transa, transb, NA, NB, nshots);
+
+#pragma omp for schedule(static) \
+private(iw, iwnA, iwnB, iwAB, iwBB) 
+	for (iw=0; iw< nfreq; iw++) {
+
+		iwnA = iw*nstationA*nshots;
+		iwnB = iw*nstationB*nshots;
+		iwAB = iw*NC*NC;
+		if (mdd==0) { /* Correlation */
+				/* cblas_cgemm(CblasRowMajor,CblasNoTrans, CblasConjTrans, NA, NB, nshots, &alpha.r, 
+				&cA[iwnA].r, NA, 
+				&cB[iwnB].r, NB, &beta.r,
+				&cC[iwAB].r, NC); */
+			cgemm_(transa, transb, &NA, &NB, &nshots, &alpha.r, 
+				&cA[iwnA].r, &NA, 
+				&cB[iwnB].r, &NB, &beta.r, 
+				&cC[iwAB].r, &NC); 	
+//				memcpy(&cC[iwAB].r, &cB[iwnA].r, sizeof(float)*2*nstationA*nshots);
+		}
+		else if (mdd==1) { /* Multi Dimensional deconvolution */
+            /* compute AB^h and BB^h */
+			iwBB = iw*nstationB*nstationB;
+			cgemm_(transa, transb, &NA, &NB, &nshots, &alpha.r, 
+				&cA[iwnA].r, &NA, 
+				&cB[iwnB].r, &NB, &beta.r, 
+				&AB[0].r, &NA);
+	
+			cgemm_(transa, transb, &NB, &NB, &nshots, &alpha.r, 
+				&cB[iwnB].r, &NB, 
+				&cB[iwnB].r, &NB, &beta.r, 
+				&BB[0].r, &NB);
+	
+			if (oBB!=NULL) memcpy(&oBB[iwBB].r, &BB[0].r, nstationB*nstationB*sizeof(complex));
+
+			/* compute inverse of BB^h as [BB^h+eps]^-1 */
+			computeMatrixInverse(BB, NB, rthm, eps_a, eps_r, numacc, eigenvalues, &eigen[iw*NB], iw, verbose);
+	
+			/* multiply with AB to get Least Squares inversion */
+			/* C = A/B => AB^h/(BB^h+eps) */
+			cgemm_(transa, transa, &NA, &NB, &NB, &alpha.r, 
+				&AB[0].r, &NA, 
+				&BB[0].r, &NB, &beta.r, 
+				&cC[iwAB].r, &NA);
+		}
+		else if (mdd==2) { /* Multi Dimensional deconvolution, but AB^H en BB^H already computed */
+
+			memcpy(&BB[0].r, &cB[iwnB].r, nstationB*nshots*sizeof(complex));
+
+			computeMatrixInverse(BB, NB, rthm, eps_a, eps_r, numacc, eigenvalues, &eigen[iw*NB], iw, verbose);
+	
+			transN = "N";
+			transN = "N";
+			cgemm_(transN, transN, &NA, &NB, &NB, &alpha.r, 
+				&cA[iwnA].r, &NA, 
+				&BB[0].r, &NB, &beta.r, 
+				&cC[iwAB].r, &NA);
+		}
+		else if (mdd==3) { /* Copy matrix A or B to memory for testing purposes */
+			memcpy(&cC[iwAB].r, &cA[iwnA].r, sizeof(complex)*nstationA*nshots);
+		}
+		else if (mdd==4) {
+			memcpy(&cC[iwAB].r, &cB[iwnB].r, sizeof(complex)*nstationB*nshots);
+		}
+		else if (mdd==5) {
+			cblas_cdotu_sub(nshots, &cA[iwnA].r, NA, &cB[iwnB].r, NB, &cC[iwnA].r);
+		}
+
+	}
+
+	free(AB);
+	free(BB);
+
+	return 0;
+}
+

MDD/docpkge.c

+../utils/docpkge.c
\ No newline at end of file

MDD/getFileInfo.c

+../utils/getFileInfo.c
\ No newline at end of file

MDD/getpars.c

+../utils/getpars.c
\ No newline at end of file

MDD/mdd.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))
+
+#ifdef _OPENMP
+int omp_get_thread_num(void);
+#endif
+double wallclock_time(void);
+void name_ext(char *filename, char *extension);
+
+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 readShotData(char *filename, float xmin, float dx, float *xrcv, float *xsrc, int *xnx, complex *cdata, int nw, int nw_low, int ngath, int nx, int nxm, int ntfft, float alpha, float scl, float conjg, int transpose, int verbose);
+
+int deconvolve(complex *cA, complex *cB, complex *cC, complex *oBB, int nfreq, int nblock, size_t nstationA, size_t nstationB, float eps_a, float eps_r, float numacc, int eigenvalues, float *eigen, int rthm, int mdd, int conjgA, int conjgB, int verbose);
+
+void writeEigen(char *file_out, float df, int nw_low, int nw_high, int nw, float *eigen, int nx, float dx, float xmin);
+void writeDatamatrix(char *file_out, complex *P, int ntfft, int ntc, int Nrec, int Nshot, int nfreq, int nw_low, float dt, int verbose);
+
+void gausstaper(float *taper, float dx, int n, float enddecay);
+
+/**************
+* ntc output samples of deconvolution result
+* note that nt (the number of samples read by the IO routine)
+* should be 2*ntc and a number efficient for FFT's
+*/
+
+/*********************** self documentation **********************/
+char *sdoc[] = {
+" ",
+" mdd - multi-dimensional deconvolution (OpenMP)",
+"  ",
+" mdd file_A= file_B= file_out= [optional parameters]",
+"  ",
+" Required parameters: ",
+" ",
+"   file_A= .................. name of file(s) which store the data in location A",
+"   file_B= .................. name of file(s) which store the data in location B",
+" ",
+" Optional parameters: ",
+" ",
+"   ntc=nt ................... number of output time samples",
+"   ntfft=nt ................. number of samples used in fft",
+"   fmin=0 ................... minimum frequency",
+"   fmax=70 .................. maximum frequency to use in deconvolution",
+" INPUT DEFINITION ",
+"   cjA=1 .................... -1 => apply complex conjugate to A",
+"   sclA=1 ................... apply scaling factor to A",
+"   tranposeA=0 .............. apply transpose to A",
+"   cjB=1 .................... -1 => apply complex conjugate to B",
+"   sclB=1 ................... apply scaling factor to B",
+"   tranposeB=0 .............. apply transpose to B",
+" MATRIX INVERSION CALCULATION ",
+"   conjgA=0 ................. apply complex conjugate-transpose to A",
+"   conjgB=1 ................. apply complex conjugate-transpose to B",
+"   rthm=0 ................... see below for options",
+"   eps_a=1e-5 ............... absolute stabilization factor for LS",
+"   eps_r=1e-4 ............... relative stabilization factor for LS",
+"   numacc=1e-6 .............. numerical accurary for SVD",
+"   ntap=0 ................... number of taper points matrix",
+"   ftap=0 ................... percentage for tapering",
+"   tap=0 .................... type of taper: 0=cos 1=exp",
+"   eigenvalues= ............. write SVD eigenvalues to file ",
+"   mdd=1 .................... mdd=0 => computes correlation ",
+" OUTPUT DEFINITION ",
+"   file_out= ................ output base name ",
+"   causal=1 ................. output causal(1), non-causal(2), both(3), or summed(4)",
+"   one_file=1 ............... write all shots into one file ",
+"   file_dmat= ............... if defined writes matrix in frequency domain",
+"   verbose=0 ................ silent option; >0 displays info",
+" ",
+" Notes: ",
+"    ntc output samples of deconvolution result",
+"    nt (the number of samples read by the IO routine)",
+" ",
+" Options for mdd= ",
+"     2 = A/(B + eps) ",
+"     1 = A*B^H/(B*B^H + eps) ",
+"     0 = A*B^H ",
+" ",
+" Option for rthm= ",
+"     0 = Least Squares QR based inversion",
+"     1 = Least Squares LU based inversion",
+"     2 = SVD inversion single precision",
+"     3 = SVD divide-and-conquer method",
+"     4 = SVD inversion double precision",
+"     5 = Least Squares LU based inversion double precision",
+"     6 = Eigenvalue based (not yet working)",
+" ",
+" author  : Jan Thorbecke : 2008 (j.w.thorbecke@tudelft.nl)",
+" ",
+NULL};
+/**************** end self doc ***********************************/
+
+int main (int argc, char **argv)
+{
+	FILE    *fpin, *fpout;
+	int		i, j, k, ret, nshots, ntraces;
+	int		size, n1, n2, ntfft, nf, causal;
+	int     verbose, fullcorr, ncorstat, err;
+	int     nt, nc, ncc, ntc, nshotA, nshotB;
+	size_t  nstationA, nstationB, nfreq, istation, jstation, iw;
+	int     pgsz, istep,jstep;
+	int     mdd;
+	int		conjgA, conjgB;
+	int 	ntap, nxm, ngath, nw, nw_low, nw_high, eigenvalues, rthm, combine, distance;
+	size_t  nwrite, cdatainSize, datainSize, cdataoutSize, stationSize, is;
+	float	dx, dt, fmin, fmax, df, eps_r, eps_a, ftap, numacc;
+	float	*rC, scl, *rl, *eigen;
+	float   f1, f2, d1, d2, sclsxgx, xmin, xmax, alpha, wshot, wpi, wrec;
+ 	float   *xrcvA, *xsrcA, *xrcvB, *xsrcB;
+	float	*taper;
+    int     *xnx;
+	float sclA,sclB, cjA, cjB;
+	int  transposeA, transposeB;
+
+
+	complex *cdataout;
+	double  t0, t1, t2, t3, tinit, twrite, tread, tdec, tfft;
+	char	*file_A, *file_B, *file_out, *file_dmat, filename[1024], number[128], *rthmName;
+	int     pe=0, root_pe=0, npes=1, ipe, size_s, one_file;
+	complex *cA, *cB, *oBB;
+	segy *hdr;
+
+	t0 = wallclock_time();
+	initargs(argc, argv);
+	requestdoc(1);
+
+	if (!getparint("verbose", &verbose)) verbose = 0;
+	if (!getparstring("file_A", &file_A)) file_A=NULL;
+	assert(file_A != NULL);
+	if (!getparstring("file_B", &file_B)) file_B=NULL;
+	assert(file_B != NULL);
+	if (!getparstring("file_out", &file_out)) file_out=NULL;
+	if (!getparstring("file_dmat", &file_dmat)) file_dmat=NULL;
+	if (!getparint("one_file", &one_file)) one_file = 1;
+
+	if (!getparfloat("fmin", &fmin)) fmin = 0.0;
+	if (!getparint("rthm", &rthm)) rthm = 0;
+	if (!getparint("combine", &combine)) combine = 0;
+	if (!getparint("causal", &causal)) causal = 1;
+	if (!getparint("ntap", &ntap)) ntap = 0;
+	if (!getparfloat("ftap", &ftap)) ftap = 0.;
+	if (!getparfloat("eps_r", &eps_r)) eps_r = 1e-4;
+	if (!getparfloat("eps_a", &eps_a)) eps_a = 1e-5;
+	if (!getparfloat("numacc", &numacc)) numacc = 1e-6;
+	if (!getparint("eigenvalues", &eigenvalues)) eigenvalues = 0;
+	if (!getparint("mdd", &mdd)) mdd = 1;
+
+	if (!getparint("transposeA", &transposeA)) transposeA = 0;
+	if (!getparfloat("sclA", &sclA)) sclA = 1.;
+	if (!getparfloat("cjA", &cjA)) cjA = 1.;
+	if (!getparint("transposeB", &transposeB)) transposeB = 0;
+	if (!getparfloat("sclB", &sclB)) sclB = 1.;
+	if (!getparfloat("cjB", &cjB)) cjB = 1.;
+
+#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_A, &n1, &n2, &nshotA, &d1, &d2, &f1, &f2, &xmin, &xmax, &sclsxgx, &nxm);
+	if (!getparint("nt", &nt)) nt=n1;
+	if (!getparint("ntc", &ntc)) ntc = n1;
+	if (!getparint("conjgA", &conjgA)) conjgA = 0;
+	if (!getparint("conjgB", &conjgB)) conjgB = 1;
+	if (!getparfloat("dt", &dt)) dt = d1;
+	if (!getparfloat("dx", &dx)) dx = d2;
+	if (!getparfloat("fmax", &fmax)) fmax = 1.0/(2.0*dt);
+
+	nstationA = n2;
+
+	nshotB = 0;
+	getFileInfo(file_B, &n1, &n2, &nshotB, &d1, &d2, &f1, &f2, &xmin, &xmax, &sclsxgx, &nxm);
+	assert( n1 == nt);
+	nstationB = n2;
+	assert( nshotA == nshotB);
+
+/*================ initializations ================*/
+
+	tinit = 0.0;
+	tfft  = 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);
+
+/* scaling of the results by Johno van IJsseldijk */ 
+    if (mdd == 0) scl = dx*dt/((float)ntfft); //correlation
+    else if (mdd==1) scl = 1/((float)ntfft)/dx/dt; // MDD
+    else if (mdd==2) scl = 1/((float)ntfft)/dx/dt; // MDD with A and B already computed (NOT TESTED)
+    else scl = 1.0/((float)ntfft); // Passing A or B through
+
+/* allocate in shared memory the in- and output data */
+
+	jstep        = nfreq*nshotA;
+	cdatainSize  = nfreq*nshotA*sizeof(complex);
+	cdataoutSize = nstationA*nstationB*nfreq*sizeof(complex);
+	cdataout     = (complex *)malloc(cdataoutSize);
+	cA           = (complex *)malloc(nstationA*cdatainSize);
+	cB           = (complex *)malloc(nstationB*cdatainSize);
+	taper        = (float *)malloc(2*nstationB*sizeof(float));
+	if (file_dmat!=NULL) oBB = (complex *)malloc(nstationB*nstationB*nfreq*sizeof(complex));
+	else oBB = NULL;
+	assert(cdataout != NULL);
+	assert(cA != NULL);
+	assert(cB != NULL);
+
+
+/* for first touch binding of allocated memory */
+#pragma omp parallel for schedule(static) private(jstation,is) default(shared)
+	for (jstation=0; jstation<nstationB; jstation++) {
+		stationSize=nstationA*nfreq*sizeof(complex);
+		is = jstation*nstationA*nfreq;
+		memset(&cdataout[is],0,stationSize);
+		memset(&cB[jstation*jstep],0,jstep*sizeof(complex));
+	}
+
+#pragma omp parallel for schedule(static) private(jstation) default(shared)
+	for (jstation=0; jstation<nstationA; jstation++) {
+		memset(&cA[jstation*jstep],0,jstep*sizeof(complex));
+	}
+
+    if (verbose) {
+        if (rthm==0) rthmName="Cholesky";
+        else if (rthm==1) rthmName="LU";
+        else if (rthm==2) rthmName="SVD single precision";
+        else if (rthm==3) rthmName="SVD divide-and-conquer";
+        else if (rthm==4) rthmName="SVD double precision";
+        else if (rthm==5) rthmName="LU double precision";
+        else if (rthm==6) rthmName="Eigenvalue double precision";
+        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,"  number t-fft ..... : %d\n", ntfft);
+		fprintf(stderr,"  Input  size ...... : %ld MB\n", (nstationA+nstationB)*cdatainSize/(1024*1024));
+		fprintf(stderr,"  Output size ...... : %ld MB\n", (cdataoutSize/((size_t)1024*1024)));
+        fprintf(stderr,"  taper points ..... : %d (%.2f %%)\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);
+        if (mdd) fprintf(stderr,"  Matrix inversion . : %s\n", rthmName);
+        else  fprintf(stderr,"  Correlation ...... : \n");
+        fprintf(stderr,"  eps_r ............ : %e\n", eps_r);
+        fprintf(stderr,"  eps_a ............ : %e\n", eps_a);
+        fprintf(stderr,"  mdd .............. : %d\n", mdd);
+    }
+
+	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;
+    readShotData(file_A, xmin, dx, xrcvA, xsrcA, xnx, cA, nw, nw_low, nshotA, nstationA, nstationA, ntfft, alpha, sclA, cjA, transposeA, verbose);
+
+    xsrcB     = (float *)calloc(nshotB,sizeof(float));
+    xrcvB     = (float *)calloc(nshotB*nstationB,sizeof(float));
+	alpha = 0.0;
+    readShotData(file_B, xmin, dx, xrcvB, xsrcB, xnx, cB, nw, nw_low, nshotB, nstationB, nstationB, ntfft, alpha, sclB, cjB, transposeB, verbose);
+
+	//cB = cA;
+
+	eigen = (float *)malloc(nfreq*nstationB*sizeof(float));
+
+	t2 = wallclock_time();
+	tread += t2-t1;
+
+#pragma omp parallel default(none) \
+	private(t1,t2,pe) \
+	shared(cA,cB,eigen,eigenvalues,numacc,eps_r,eps_a) \
+	shared(nstationA,nstationB,verbose,cdatainSize) \
+	shared(rthm,mdd,nfreq,nshotA,conjgA,conjgB) \
+	shared(cdataout,oBB)
+{ /* start of OpenMP parallel part */
+
+
+#ifdef _OPENMP
+	pe = omp_get_thread_num();
+#endif
+
+	/* compute deconvolution */
+	deconvolve(cA, cB, cdataout, oBB, nfreq, nshotA, nstationA, nstationB, 
+		eps_a, eps_r, numacc, eigenvalues, eigen, rthm, mdd, conjgA, conjgB, verbose);
+
+} /*end of parallel OpenMP part */
+
+	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 deconvolution     = %.3f\n", tdec);
+	}
+
+/* for writing out combined shots cA */
+	free(cA);
+	free(cB);
+
+/* Inverse FFT of deconvolution results */
+/* This is done for every deconvolution component seperately */
+
+	rC = (float *)malloc(nstationA*ntc*sizeof(float));
+	assert(rC != NULL);
+
+/*
+#pragma omp parallel default(none) \
+	private(istation,jstation,pe,j,i,t1,t2,t3,hdr,rl) \
+	private(filename, k, fpout, nwrite, cA, iw,number) \
+	shared(tfft) \
+	shared(rC,dt,ntc,file_out) \
+	shared(nt,nstationA,nstationB,verbose,err,ntfft,t0,twrite) \
+	shared(nfreq,stderr,stdout, nshotA, nshotB, nw_low, causal) \
+	shared(cdataout,istep,jstep,one_file)
+*/
+//{ /* start of OpenMP parallel part */
+//#ifdef _OPENMP
+//	pe = omp_get_thread_num();
+//#else 
+    pe = 0;
+//#endif
+
+	rl  = (float *)calloc(ntfft,sizeof(float));
+	cA  = (complex *)calloc(ntfft,sizeof(complex));
+	hdr = (segy *)calloc(1,sizeof(segy));
+
+/* for writing out combined shots cA */
+
+	tfft   = 0.0;
+	twrite = 0.0;
+	if (one_file && pe==0) {
+		strcpy(filename, file_out);
+		if (verbose>2) fprintf(stderr,"writing all output shot into file %s\n", filename);
+		fpout = fopen( filename, "w+" );
+	}
+//#pragma omp for
+	for (jstation=0; jstation<nstationB; jstation++) {
+		/* FFT */
+		t1 = wallclock_time();
+		for (istation=0; istation<nstationA; istation++) {
+			memset(cA,0,ntfft*sizeof(complex));
+			for (iw=0;iw<nfreq;iw++) {
+				cA[iw+nw_low].r = cdataout[(iw*nstationB+jstation)*nstationA+istation].r*scl;
+				cA[iw+nw_low].i = cdataout[(iw*nstationB+jstation)*nstationA+istation].i*scl;
+			}
+			cr1fft(cA, rl, ntfft, 1);
+			memcpy(&rC[istation*ntc],rl,ntc*sizeof(float));
+
+			if (causal==1) {
+				memcpy(&rC[istation*ntc],rl,ntc*sizeof(float));
+			}
+			else if (causal==2) {
+				rC[istation*ntc] = rl[0];
+				for (j=1;j<ntc; j++) {
+					rC[istation*ntc+j] = rl[ntfft-j];
+				}
+			}
+			else if (causal==3) {
+				for (j=1;j<=(ntc/2); j++) {
+					rC[istation*ntc+ntc/2-j] = rl[ntfft-j];
+				}
+				for (j=ntc/2;j<ntc; j++) {
+					rC[istation*ntc+j] = rl[j-ntc/2];
+				}
+			}
+			else if (causal==4) {
+				rC[istation*ntc] = rl[0];
+				for (j=1;j<ntc; j++) {
+					rC[istation*ntc+j] = rl[ntfft-j] + rl[j];
+				}
+			}
+		}
+		t2 = wallclock_time();
+		tfft += t2-t1;
+
+		if (pe == 0) {
+			/* write data to file */
+			hdr[0].d1  = dt;
+			if (causal == 3) hdr[0].f1=-0.5*ntc*dt;
+			else hdr[0].f1=0.0;
+			hdr[0].dt  = (int)(dt*1000000);
+			hdr[0].ns  = ntc;
+			hdr[0].fldr  = jstation+1;
+			hdr[0].scalco = -1000;
+			hdr[0].scalel = -1000;
+			hdr[0].trid = 1;
+			hdr[0].f2 = f2;
+			hdr[0].d2 = dx;
+//			hdr[0].trwf = nstationA;
+			hdr[0].sx = NINT((f2+dx*jstation)*1000);
+			hdr[0].ntr = nstationA*nstationB;
+			if (!one_file) {
+				strcpy(filename, file_out);
+				sprintf(number,"Station%03d\0",jstation+1);
+				name_ext(filename, number);
+				if (verbose>3) fprintf(stderr,"writing to file %s\n", filename);
+				fpout = fopen( filename, "w+" );
+			}
+			for (istation=0; istation<nstationA; istation++) {
+				hdr[0].tracl = istation+1;
+				hdr[0].gx = NINT((f2+dx*istation)*1000);
+				hdr[0].offset = NINT((f2+dx*istation));
+				nwrite = fwrite( hdr, 1, TRCBYTES, fpout );
+				assert (nwrite == TRCBYTES);
+				nwrite = fwrite( &rC[istation*ntc], sizeof(float), ntc, fpout );
+				assert (nwrite == ntc);
+			}
+			if (!one_file) {
+				fflush(fpout);
+				fclose(fpout);
+			}
+			t3 = wallclock_time();
+			twrite += t3-t2;
+//			fprintf(stderr,"write %f and fft %f for %d\n",twrite, tfft, jstation);
+		}
+	}
+	if (one_file && pe==0) {
+		fflush(fpout);
+		fclose(fpout);
+	}
+	free(cA);
+	free(rl);
+//}
+
+	free(rC);
+	free(cdataout);
+
+	if (eigenvalues) {
+		writeEigen(file_out, df, nw_low, nw_high, nfreq, eigen, nstationB, dx, f2);
+	}
+	free(eigen);
+
+	/* if file_dmat write frequency slices of matrix */
+	if (file_dmat!=NULL) {
+		t2 = wallclock_time();
+		strcpy(filename, file_dmat);
+		fpout = fopen( filename, "w+" );
+		hdr[0].d1  = df;
+		hdr[0].dt  = (int)(df*1000000);
+		hdr[0].ns  = nfreq;
+		hdr[0].trid  = 111;
+/*
+		for (iw=0;iw<nfreq;iw++) {
+			hdr[0].fldr  = iw+1;
+//			sprintf(number,"Station%03d\0",jstation+1);
+//			name_ext(filename, number);
+//			if (verbose>3) fprintf(stderr,"writing to file %s\n", filename);
+//			fpout = fopen( filename, "w+" );
+			twrite = 0.0;
+			for (istation=0; istation<nstationB; istation++) {
+				hdr[0].tracl = istation+1;
+				nwrite = fwrite( hdr, 1, TRCBYTES, fpout );
+				assert (nwrite == TRCBYTES);
+//				nwrite = fwrite( &oBB[iw*nstationB*nstationB+istation].r, sizeof(complex), nfreq, fpout );
+//				assert (nwrite == nfreq);
+			}
+		}
+*/
+		fflush(fpout);
+		fclose(fpout);
+		t3 = wallclock_time();
+		twrite += t3-t2;
+		free(oBB);
+	}
+	free(hdr);
+
+/*================ end ================*/
+
+	if (verbose) {
+		t3 = wallclock_time();
+		fprintf(stderr,"CPU-time inverse FFT's     = %.3f\n", tfft);
+		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);
+	}
+
+	return 0;
+}
+
+void gausstaper(float *taper, float dx, int n, float enddecay)
+{
+	int 	ix, hn;
+	float 	dist, sigma2;
+
+	if (enddecay > 0.999) {
+		for (ix = 0; ix < n; ix++) taper[ix] = 1.0;
+		return;
+	}
+
+	hn = (n-1)/2;
+	sigma2 = (hn*dx*hn*dx)/(log(enddecay));
+
+	for (ix = 0; ix <= hn; ix++) {
+		dist = ix*dx;
+		taper[hn+ix] = exp(dist*dist/sigma2);
+	}
+
+	for (ix = 0; ix < hn; ix++) 
+		taper[ix] = taper[n-1-ix];
+
+	return;
+}
+
+void writeDatamatrix(char *file_out, complex *P, int ntfft, int ntc, int Nrec, int Nshot, int nfreq, int nw_low, float dt, int verbose)
+{
+	FILE *fpout;
+	char filename[1024];
+	size_t  nwrite;
+	int jstation, istation, iw;
+	float *rl, *rC;
+	complex *cA;
+	segy *hdr;
+
+	rC = (float *)malloc(Nrec*ntc*sizeof(float));
+	rl  = (float *)calloc(ntfft,sizeof(float));
+	cA  = (complex *)calloc(ntfft,sizeof(complex));
+	hdr = (segy *)calloc(1,sizeof(segy));
+
+/* for writing out combined shots cA */
+
+	strcpy(filename, file_out);
+	if (verbose>2) fprintf(stderr,"writing all output shot into file %s\n", filename);
+	fpout = fopen( file_out, "w+" );
+	for (jstation=0; jstation<Nshot; jstation++) {
+
+		/* FFT */
+		for (istation=0; istation<Nrec; istation++) {
+			memset(cA,0,ntfft*sizeof(complex));
+			for (iw=0;iw<nfreq;iw++) {
+				cA[iw+nw_low] = P[(iw*Nshot+jstation)*Nrec+istation];
+			}
+			cr1fft(cA, rl, ntfft, 1);
+			memcpy(&rC[istation*ntc],rl,ntc*sizeof(float));
+		}
+
+		/* write data to file */
+		hdr[0].d1  = dt;
+		hdr[0].dt  = (int)(dt*1000000);
+		hdr[0].ns  = ntc;
+		hdr[0].fldr  = jstation+1;
+		for (istation=0; istation<Nrec; istation++) {
+			hdr[0].tracl = istation+1;
+			nwrite = fwrite( hdr, 1, TRCBYTES, fpout );
+			assert (nwrite == TRCBYTES);
+			nwrite = fwrite( &rC[istation*ntc], sizeof(float), ntc, fpout );
+			assert (nwrite == ntc);
+		}
+	}
+
+	free(cA);
+	free(rl);
+	free(rC);
+	return;
+}
+

MDD/name_ext.c

+../utils/name_ext.c
\ No newline at end of file

MDD/par.h

+../utils/par.h
\ No newline at end of file

MDD/readShotData.c

+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include "segy.h"
+#include <assert.h>
+
+extern FILE *fopen64 (__const char *__restrict __filename,
+                      __const char *__restrict __modes);
+
+typedef struct { /* complex number */
+        float r,i;
+} complex;
+
+#define NINT(x) ((int)((x)>0.0?(x)+0.5:(x)-0.5))
+
+int optncr(int n);
+void cc1fft(complex *data, int n, int sign);
+void rc1fft(float *rdata, complex *cdata, int n, int sign);
+
+int compare(const void *a, const void *b) 
+{ return (*(float *)b-*(float *)a); }
+
+int readShotData(char *filename, float xmin, float dx, float *xrcv, float *xsrc, int *xnx, complex *cdata, int nw, int nw_low, int ngath, int nx, int nxm, int ntfft, float alpha, float scale, float conjg, int transpose, int verbose)
+{
+	FILE *fp;
+	segy hdr;
+	size_t nread;
+	int fldr_shot, sx_shot, itrace, one_shot, igath, iw, i, j, k;
+	int end_of_file, nt, ir, is;
+	float scl, dt, *trace;
+	complex *ctrace;
+
+	/* Reading first header  */
+
+	if (filename == NULL) fp = stdin;
+	else fp = fopen64( filename, "r" );
+	if ( fp == NULL ) {
+		fprintf(stderr,"input file %s has an error\n", filename);
+		perror("error in opening file: ");
+		fflush(stderr);
+		return -1;
+	}
+
+	fseek(fp, 0, SEEK_SET);
+	nread = fread( &hdr, 1, TRCBYTES, fp );
+	assert(nread == TRCBYTES);
+	if (hdr.scalco < 0) scl = 1.0/fabs(hdr.scalco);
+	else if (hdr.scalco == 0) scl = 1.0;
+	else scl = hdr.scalco;
+	fseek(fp, 0, SEEK_SET);
+
+	nt        = hdr.ns;
+
+	trace  = (float *)calloc(nx*ntfft,sizeof(float));
+	ctrace = (complex *)malloc(ntfft*sizeof(complex));
+
+	end_of_file = 0;
+	one_shot    = 1;
+	igath       = 0;
+
+	/* Read shots in file */
+
+	while (!end_of_file) {
+
+		/* start reading data (shot records) */
+		itrace = 0;
+		nread = fread( &hdr, 1, TRCBYTES, fp );
+		if (nread != TRCBYTES) { /* no more data in file */
+			break;
+		}
+
+		sx_shot  = hdr.sx;
+		fldr_shot  = hdr.fldr;
+		xsrc[igath] = sx_shot*scl;
+		xnx[igath]=0;
+		/* read in all traces within a shot */
+		while (one_shot) {
+			xrcv[igath*nxm+itrace] = hdr.gx*scl;
+			nread = fread( &trace[itrace*ntfft], sizeof(float), nt, fp );
+			assert (nread == hdr.ns);
+			itrace++;
+			xnx[igath]+=1;
+
+			/* read next hdr of next trace */
+			nread = fread( &hdr, 1, TRCBYTES, fp );
+			if (nread != TRCBYTES) { 
+				one_shot = 0;
+				end_of_file = 1;
+				break;
+			}
+			if ((sx_shot != hdr.sx) || (fldr_shot != hdr.fldr) ) break;
+		}
+
+		for (i=0; i<itrace; i++) {
+			/* apply alpha factor */
+			if (alpha != 0.0) {
+        		for (j=0; j<nt; j++) {
+					trace[i*ntfft+j] *= exp(alpha*j*dt);
+				}
+			}
+        	for (j=nt; j<ntfft; j++) {
+				trace[i*ntfft+j] = 0.0;
+			}
+
+			/* transform to frequency domain */
+        	rc1fft(&trace[i*ntfft],ctrace,ntfft,-1);
+
+			if (transpose == 0) {
+        		for (iw=0; iw<nw; iw++) {
+        			cdata[iw*ngath*nx+igath*nx+i].r = scale*ctrace[nw_low+iw].r;
+        			cdata[iw*ngath*nx+igath*nx+i].i = conjg*scale*ctrace[nw_low+iw].i;
+        		}
+			}
+			else {
+        		for (iw=0; iw<nw; iw++) {
+        			cdata[iw*ngath*nx+i*ngath+igath].r = scale*ctrace[nw_low+iw].r;
+        			cdata[iw*ngath*nx+i*ngath+igath].i = conjg*scale*ctrace[nw_low+iw].i;
+        		}
+			}
+		}
+
+		if (verbose>2) {
+			fprintf(stderr,"finished reading shot %d (%d) with %d traces\n",sx_shot,igath,itrace);
+		}
+
+		if (itrace != 0) { /* end of shot record */
+			fseek( fp, -TRCBYTES, SEEK_CUR );
+			igath++;
+		}
+		else {
+			end_of_file = 1;
+		}
+	}
+
+	free(ctrace);
+	free(trace);
+
+	return 0;
+}
+
+

MDD/segy.h

+/* Copyright (c) Colorado School of Mines, 2011.*/
+/* All rights reserved.                       */
+
+/* segy.h - include file for SEGY traces
+ *
+ * declarations for:
+ *	typedef struct {} segy - the trace identification header
+ *	typedef struct {} bhed - binary header
+ *
+ * Note:
+ *	If header words are added, run the makefile in this directory
+ *	to recreate hdr.h.
+ *
+ * Reference:
+ *	K. M. Barry, D. A. Cavers and C. W. Kneale, "Special Report:
+ *		Recommended Standards for Digital Tape Formats",
+ *		Geophysics, vol. 40, no. 2 (April 1975), P. 344-352.
+ *	
+ * $Author: john $
+ * $Source: /usr/local/cwp/src/su/include/RCS/segy.h,v $
+ * $Revision: 1.33 $ ; $Date: 2011/11/11 23:56:14 $
+ */ 
+
+#include <limits.h>
+#include "par.h"
+
+#ifndef SEGY_H
+#define SEGY_H
+#define TRCBYTES		240
+
+#define SU_NFLTS	32767	/* Arbitrary limit on data array size	*/
+
+
+/* TYPEDEFS */
+typedef struct {	/* segy - trace identification header */
+
+	int tracl;	/* Trace sequence number within line
+			   --numbers continue to increase if the
+			   same line continues across multiple
+			   SEG Y files.
+			   byte# 1-4
+			 */
+
+	int tracr;	/* Trace sequence number within SEG Y file
+			   ---each file starts with trace sequence
+			   one
+			   byte# 5-8
+			 */
+
+	int fldr;	/* Original field record number
+			   byte# 9-12 
+			*/
+
+	int tracf;	/* Trace number within original field record 
+			   byte# 13-16
+			*/
+
+	int ep;		/* energy source point number
+			   ---Used when more than one record occurs
+			   at the same effective surface location.
+			   byte# 17-20
+			 */
+
+	int cdp;	/* Ensemble number (i.e. CDP, CMP, CRP,...) 
+			   byte# 21-24
+			*/
+
+	int cdpt;	/* trace number within the ensemble
+			   ---each ensemble starts with trace number one.
+			   byte# 25-28
+			 */
+
+	short trid;	/* trace identification code:
+			-1 = Other
+		         0 = Unknown
+			 1 = Seismic data
+			 2 = Dead
+			 3 = Dummy
+			 4 = Time break
+			 5 = Uphole
+			 6 = Sweep
+			 7 = Timing
+			 8 = Water break
+			 9 = Near-field gun signature
+			10 = Far-field gun signature
+			11 = Seismic pressure sensor
+			12 = Multicomponent seismic sensor
+				- Vertical component
+			13 = Multicomponent seismic sensor
+				- Cross-line component 
+			14 = Multicomponent seismic sensor
+				- in-line component 
+			15 = Rotated multicomponent seismic sensor
+				- Vertical component
+			16 = Rotated multicomponent seismic sensor
+				- Transverse component
+			17 = Rotated multicomponent seismic sensor
+				- Radial component
+			18 = Vibrator reaction mass
+			19 = Vibrator baseplate
+			20 = Vibrator estimated ground force
+			21 = Vibrator reference
+			22 = Time-velocity pairs
+			23 ... N = optional use 
+				(maximum N = 32,767)
+
+			Following are CWP id flags:
+
+			109 = autocorrelation
+			110 = Fourier transformed - no packing
+			     xr[0],xi[0], ..., xr[N-1],xi[N-1]
+			111 = Fourier transformed - unpacked Nyquist
+			     xr[0],xi[0],...,xr[N/2],xi[N/2]
+			112 = Fourier transformed - packed Nyquist
+	 		     even N:
+			     xr[0],xr[N/2],xr[1],xi[1], ...,
+				xr[N/2 -1],xi[N/2 -1]
+				(note the exceptional second entry)
+			     odd N:
+			     xr[0],xr[(N-1)/2],xr[1],xi[1], ...,
+				xr[(N-1)/2 -1],xi[(N-1)/2 -1],xi[(N-1)/2]
+				(note the exceptional second & last entries)
+			113 = Complex signal in the time domain
+			     xr[0],xi[0], ..., xr[N-1],xi[N-1]
+			114 = Fourier transformed - amplitude/phase
+			     a[0],p[0], ..., a[N-1],p[N-1]
+			115 = Complex time signal - amplitude/phase
+			     a[0],p[0], ..., a[N-1],p[N-1]
+			116 = Real part of complex trace from 0 to Nyquist
+			117 = Imag part of complex trace from 0 to Nyquist
+			118 = Amplitude of complex trace from 0 to Nyquist
+			119 = Phase of complex trace from 0 to Nyquist
+			121 = Wavenumber time domain (k-t)
+			122 = Wavenumber frequency (k-omega)
+			123 = Envelope of the complex time trace
+			124 = Phase of the complex time trace
+			125 = Frequency of the complex time trace
+			130 = Depth-Range (z-x) traces
+			201 = Seismic data packed to bytes (by supack1)
+			202 = Seismic data packed to 2 bytes (by supack2)
+			   byte# 29-30
+			*/
+
+	short nvs;	/* Number of vertically summed traces yielding
+			   this trace. (1 is one trace, 
+			   2 is two summed traces, etc.)
+			   byte# 31-32
+			 */
+
+	short nhs;	/* Number of horizontally summed traces yielding
+			   this trace. (1 is one trace
+			   2 is two summed traces, etc.)
+			   byte# 33-34
+			 */
+
+	short duse;	/* Data use:
+				1 = Production
+				2 = Test
+			   byte# 35-36
+			 */
+
+	int offset;	/* Distance from the center of the source point 
+			   to the center of the receiver group 
+			   (negative if opposite to direction in which 
+			   the line was shot).
+			   byte# 37-40
+			 */
+
+	int gelev;	/* Receiver group elevation from sea level
+			   (all elevations above the Vertical datum are 
+			   positive and below are negative).
+			   byte# 41-44
+			 */
+
+	int selev;	/* Surface elevation at source.
+			   byte# 45-48
+			 */
+
+	int sdepth;	/* Source depth below surface (a positive number).
+			   byte# 49-52
+			 */
+
+	int gdel;	/* Datum elevation at receiver group.
+			   byte# 53-56
+			*/
+
+	int sdel;	/* Datum elevation at source.
+			   byte# 57-60
+			*/
+
+	int swdep;	/* Water depth at source.
+			   byte# 61-64
+			*/
+
+	int gwdep;	/* Water depth at receiver group.
+			   byte# 65-68
+			*/
+
+	short scalel;	/* Scalar to be applied to the previous 7 entries
+			   to give the real value. 
+			   Scalar = 1, +10, +100, +1000, +10000.
+			   If positive, scalar is used as a multiplier,
+			   if negative, scalar is used as a divisor.
+			   byte# 69-70
+			 */
+
+	short scalco;	/* Scalar to be applied to the next 4 entries
+			   to give the real value. 
+			   Scalar = 1, +10, +100, +1000, +10000.
+			   If positive, scalar is used as a multiplier,
+			   if negative, scalar is used as a divisor.
+			   byte# 71-72
+			 */
+
+	int  sx;	/* Source coordinate - X 
+			   byte# 73-76
+			*/
+
+	int  sy;	/* Source coordinate - Y 
+			   byte# 77-80
+			*/
+
+	int  gx;	/* Group coordinate - X 
+			   byte# 81-84
+			*/
+
+	int  gy;	/* Group coordinate - Y 
+			   byte# 85-88
+			*/
+
+	short counit;	/* Coordinate units: (for previous 4 entries and
+				for the 7 entries before scalel)
+			   1 = Length (meters or feet)
+			   2 = Seconds of arc
+			   3 = Decimal degrees
+			   4 = Degrees, minutes, seconds (DMS)
+
+			In case 2, the X values are longitude and 
+			the Y values are latitude, a positive value designates
+			the number of seconds east of Greenwich
+				or north of the equator
+
+			In case 4, to encode +-DDDMMSS
+			counit = +-DDD*10^4 + MM*10^2 + SS,
+			with scalco = 1. To encode +-DDDMMSS.ss
+			counit = +-DDD*10^6 + MM*10^4 + SS*10^2 
+			with scalco = -100.
+			   byte# 89-90
+			*/
+
+	short wevel;	/* Weathering velocity. 
+			   byte# 91-92
+			*/
+
+	short swevel;	/* Subweathering velocity. 
+			   byte# 93-94
+			*/
+
+	short sut;	/* Uphole time at source in milliseconds. 
+			   byte# 95-96
+			*/
+
+	short gut;	/* Uphole time at receiver group in milliseconds. 
+			   byte# 97-98
+			*/
+
+	short sstat;	/* Source static correction in milliseconds. 
+			   byte# 99-100
+			*/
+
+	short gstat;	/* Group static correction  in milliseconds.
+			   byte# 101-102
+			*/
+
+	short tstat;	/* Total static applied  in milliseconds.
+			   (Zero if no static has been applied.)
+			   byte# 103-104
+			*/
+
+	short laga;	/* Lag time A, time in ms between end of 240-
+			   byte trace identification header and time
+			   break, positive if time break occurs after
+			   end of header, time break is defined as
+			   the initiation pulse which maybe recorded
+			   on an auxiliary trace or as otherwise
+			   specified by the recording system 
+			   byte# 105-106
+			*/
+
+	short lagb;	/* lag time B, time in ms between the time break
+			   and the initiation time of the energy source,
+			   may be positive or negative 
+			   byte# 107-108
+			*/
+
+	short delrt;	/* delay recording time, time in ms between
+			   initiation time of energy source and time
+			   when recording of data samples begins
+			   (for deep water work if recording does not
+			   start at zero time) 
+			   byte# 109-110
+			*/
+
+	short muts;	/* mute time--start 
+			   byte# 111-112
+			*/
+
+	short mute;	/* mute time--end 
+			   byte# 113-114
+			*/
+
+	unsigned short ns;	/* number of samples in this trace 
+			   byte# 115-116
+			*/
+
+	unsigned short dt;	/* sample interval; in micro-seconds
+			   byte# 117-118
+			*/
+
+	short gain;	/* gain type of field instruments code:
+				1 = fixed
+				2 = binary
+				3 = floating point
+				4 ---- N = optional use 
+			   byte# 119-120
+			*/
+
+	short igc;	/* instrument gain constant 
+			   byte# 121-122
+			*/
+
+	short igi;	/* instrument early or initial gain 
+			   byte# 123-124
+			*/
+
+	short corr;	/* correlated:
+				1 = no
+				2 = yes 
+			   byte# 125-126
+			*/
+
+	short sfs;	/* sweep frequency at start 
+			   byte# 127-128
+			*/
+
+	short sfe;	/* sweep frequency at end
+			   byte# 129-130
+			*/
+
+	short slen;	/* sweep length in ms 
+			   byte# 131-132
+			*/
+
+	short styp;	/* sweep type code:
+				1 = linear
+				2 = cos-squared
+				3 = other
+			   byte# 133-134
+			*/
+
+	short stas;	/* sweep trace length at start in ms
+			   byte# 135-136
+			*/
+
+	short stae;	/* sweep trace length at end in ms 
+			   byte# 137-138
+			*/
+
+	short tatyp;	/* taper type: 1=linear, 2=cos^2, 3=other 
+			   byte# 139-140
+			*/
+
+	short afilf;	/* alias filter frequency if used 
+			   byte# 141-142
+			*/
+
+	short afils;	/* alias filter slope
+			   byte# 143-144
+			*/
+
+	short nofilf;	/* notch filter frequency if used
+			   byte# 145-146
+			*/
+
+	short nofils;	/* notch filter slope
+			   byte# 147-148
+			*/
+
+	short lcf;	/* low cut frequency if used
+			   byte# 149-150
+			*/
+
+	short hcf;	/* high cut frequncy if used
+			   byte# 151-152
+			*/
+
+	short lcs;	/* low cut slope
+			   byte# 153-154
+			*/
+
+	short hcs;	/* high cut slope
+			   byte# 155-156
+			*/
+
+	short year;	/* year data recorded
+			   byte# 157-158
+			*/
+
+	short day;	/* day of year
+			   byte# 159-160
+			*/
+
+	short hour;	/* hour of day (24 hour clock) 
+			   byte# 161-162
+			*/
+
+	short minute;	/* minute of hour
+			   byte# 163-164
+			*/
+
+	short sec;	/* second of minute
+			   byte# 165-166
+			*/
+
+	short timbas;	/* time basis code:
+				1 = local
+				2 = GMT
+				3 = other
+			   byte# 167-168
+			*/
+
+	short trwf;	/* trace weighting factor, defined as 1/2^N
+			   volts for the least sigificant bit
+			   byte# 169-170
+			*/
+
+	short grnors;	/* geophone group number of roll switch
+			   position one
+			   byte# 171-172
+			*/
+
+	short grnofr;	/* geophone group number of trace one within
+			   original field record
+			   byte# 173-174
+			*/
+
+	short grnlof;	/* geophone group number of last trace within
+			   original field record
+			   byte# 175-176
+			*/
+
+	short gaps;	/* gap size (total number of groups dropped)
+			   byte# 177-178
+			*/
+
+	short otrav;	/* overtravel taper code:
+				1 = down (or behind)
+				2 = up (or ahead)
+			   byte# 179-180
+			*/
+
+#ifdef SLTSU_SEGY_H  /* begin Unocal SU segy.h differences */
+
+
+	/* cwp local assignments */
+	float d1;	/* sample spacing for non-seismic data
+			   byte# 181-184
+			*/
+
+	float f1;	/* first sample location for non-seismic data
+			   byte# 185-188
+			*/
+
+	float d2;	/* sample spacing between traces
+			   byte# 189-192
+			*/
+
+	float f2;	/* first trace location
+			   byte# 193-196
+			*/
+
+	float ungpow;	/* negative of power used for dynamic
+			   range compression
+			   byte# 197-200
+			*/
+
+	float unscale;	/* reciprocal of scaling factor to normalize
+			   range
+			   byte# 201-204
+			*/
+
+	short mark;	/* mark selected traces
+			   byte# 205-206
+			*/
+
+	/* SLTSU local assignments */ 
+	short mutb;	/* mute time at bottom (start time)
+			   bottom mute ends at last sample
+			   byte# 207-208
+			*/
+	float dz;	/* depth sampling interval in (m or ft)
+			if =0.0, input are time samples
+			   byte# 209-212
+			*/
+
+	float fz;	/* depth of first sample in (m or ft)
+			   byte# 213-116
+			*/
+
+	short n2;	/* number of traces per cdp or per shot
+			   byte# 217-218
+			*/
+
+        short shortpad; /* alignment padding
+			   byte# 219-220
+			*/
+
+	int ntr; 	/* number of traces
+			   byte# 221-224
+			*/
+
+	/* SLTSU local assignments end */ 
+
+	short unass[8];	/* unassigned
+			   byte# 225-240
+			*/
+
+#else
+
+	/* cwp local assignments */
+	float d1;	/* sample spacing for non-seismic data
+			   byte# 181-184
+			*/
+
+	float f1;	/* first sample location for non-seismic data
+			   byte# 185-188
+			*/
+
+	float d2;	/* sample spacing between traces
+			   byte# 189-192
+			*/
+
+	float f2;	/* first trace location
+			   byte# 193-196
+			*/
+
+	float ungpow;	/* negative of power used for dynamic
+			   range compression
+			   byte# 197-200
+			*/
+
+	float unscale;	/* reciprocal of scaling factor to normalize
+			   range
+			   byte# 201-204
+			*/
+
+	int ntr; 	/* number of traces
+			   byte# 205-208
+			*/
+
+	short mark;	/* mark selected traces
+			   byte# 209-210
+			*/
+
+        short shortpad; /* alignment padding
+			   byte# 211-212
+			*/
+
+
+	short unass[14];	/* unassigned--NOTE: last entry causes 
+			   a break in the word alignment, if we REALLY
+			   want to maintain 240 bytes, the following
+			   entry should be an odd number of short/UINT2
+			   OR do the insertion above the "mark" keyword
+			   entry
+			   byte# 213-240
+			*/
+#endif
+
+} segy;
+
+
+typedef struct {	/* bhed - binary header */
+
+	int jobid;	/* job identification number */
+
+	int lino;	/* line number (only one line per reel) */
+
+	int reno;	/* reel number */
+
+	short ntrpr;	/* number of data traces per record */
+
+        short nart;	/* number of auxiliary traces per record */
+
+	unsigned short hdt; /* sample interval in micro secs for this reel */
+
+	unsigned short dto; /* same for original field recording */
+
+	unsigned short hns; /* number of samples per trace for this reel */
+
+	unsigned short nso; /* same for original field recording */
+
+	short format;	/* data sample format code:
+				1 = floating point, 4 byte (32 bits)
+				2 = fixed point, 4 byte (32 bits)
+				3 = fixed point, 2 byte (16 bits)
+				4 = fixed point w/gain code, 4 byte (32 bits)
+				5 = IEEE floating point, 4 byte (32 bits)
+				8 = two's complement integer, 1 byte (8 bits)
+			*/
+
+	short fold;	/* CDP fold expected per CDP ensemble */
+
+	short tsort;	/* trace sorting code: 
+				1 = as recorded (no sorting)
+				2 = CDP ensemble
+				3 = single fold continuous profile
+				4 = horizontally stacked */
+
+	short vscode;	/* vertical sum code:
+				1 = no sum
+				2 = two sum ...
+				N = N sum (N = 32,767) */
+
+	short hsfs;	/* sweep frequency at start */
+
+	short hsfe;	/* sweep frequency at end */
+
+	short hslen;	/* sweep length (ms) */
+
+	short hstyp;	/* sweep type code:
+				1 = linear
+				2 = parabolic
+				3 = exponential
+				4 = other */
+
+	short schn;	/* trace number of sweep channel */
+
+	short hstas;	/* sweep trace taper length at start if
+			   tapered (the taper starts at zero time
+			   and is effective for this length) */
+
+	short hstae;	/* sweep trace taper length at end (the ending
+			   taper starts at sweep length minus the taper
+			   length at end) */
+
+	short htatyp;	/* sweep trace taper type code:
+				1 = linear
+				2 = cos-squared
+				3 = other */
+
+	short hcorr;	/* correlated data traces code:
+				1 = no
+				2 = yes */
+
+	short bgrcv;	/* binary gain recovered code:
+				1 = yes
+				2 = no */
+
+	short rcvm;	/* amplitude recovery method code:
+				1 = none
+				2 = spherical divergence
+				3 = AGC
+				4 = other */
+
+	short mfeet;	/* measurement system code:
+				1 = meters
+				2 = feet */
+
+	short polyt;	/* impulse signal polarity code:
+				1 = increase in pressure or upward
+				    geophone case movement gives
+				    negative number on tape
+				2 = increase in pressure or upward
+				    geophone case movement gives
+				    positive number on tape */
+
+	short vpol;	/* vibratory polarity code:
+				code	seismic signal lags pilot by
+				1	337.5 to  22.5 degrees
+				2	 22.5 to  67.5 degrees
+				3	 67.5 to 112.5 degrees
+				4	112.5 to 157.5 degrees
+				5	157.5 to 202.5 degrees
+				6	202.5 to 247.5 degrees
+				7	247.5 to 292.5 degrees
+				8	293.5 to 337.5 degrees */
+
+	short hunass[170];	/* unassigned */
+
+} bhed;
+
+/* DEFINES */
+#define gettr(x)	fgettr(stdin, (x))
+#define vgettr(x)	fvgettr(stdin, (x))
+#define puttr(x)	fputtr(stdout, (x))
+#define vputtr(x)	fvputtr(stdout, (x))
+#define gettra(x, y)    fgettra(stdin, (x), (y))
+
+
+/* TOTHER represents "other"					*/
+#define		TOTHER		-1	
+/* TUNK represents time traces of an unknown type		*/
+#define		TUNK		0
+/* TREAL represents real time traces 				*/
+#define		TREAL		1
+/* TDEAD represents dead time traces 				*/
+#define		TDEAD		2
+/* TDUMMY represents dummy time traces 				*/
+#define		TDUMMY		3
+/* TBREAK represents time break traces 				*/
+#define		TBREAK		4
+/* UPHOLE represents uphole traces 				*/
+#define		UPHOLE		5
+/* SWEEP represents sweep traces 				*/
+#define		SWEEP		6
+/* TIMING represents timing traces 				*/
+#define		TIMING		7
+/* WBREAK represents timing traces 				*/
+#define		WBREAK		8
+/* NFGUNSIG represents near field gun signature 		*/
+#define		NFGUNSIG	9	
+/* FFGUNSIG represents far field gun signature	 		*/
+#define		FFGUNSIG	10
+/* SPSENSOR represents seismic pressure sensor	 		*/
+#define		SPSENSOR	11
+/* TVERT represents multicomponent seismic sensor 
+	- vertical component */
+#define		TVERT		12
+/* TXLIN represents multicomponent seismic sensor 
+	- cross-line component */
+#define		TXLIN		13
+/* TINLIN represents multicomponent seismic sensor 
+	- in-line component */
+#define		TINLIN	14
+/* ROTVERT represents rotated multicomponent seismic sensor 
+	- vertical component */
+#define		ROTVERT		15
+/* TTRANS represents rotated multicomponent seismic sensor 
+	- transverse component */
+#define		TTRANS		16
+/* TRADIAL represents rotated multicomponent seismic sensor 
+	- radial component */
+#define		TRADIAL		17
+/* VRMASS represents vibrator reaction mass */
+#define		VRMASS		18
+/* VBASS represents vibrator baseplate */
+#define		VBASS		19
+/* VEGF represents vibrator estimated ground force */
+#define		VEGF		20
+/* VREF represents vibrator reference */
+#define		VREF		21
+
+/*** CWP trid assignments ***/
+/* ACOR represents autocorrelation  */
+#define		ACOR		109
+/* FCMPLX represents fourier transformed - no packing 
+   xr[0],xi[0], ..., xr[N-1],xi[N-1] */
+#define		FCMPLX		110
+/* FUNPACKNYQ represents fourier transformed - unpacked Nyquist
+   xr[0],xi[0],...,xr[N/2],xi[N/2] */
+#define		FUNPACKNYQ	111
+/* FTPACK represents fourier transformed - packed Nyquist
+   even N: xr[0],xr[N/2],xr[1],xi[1], ...,
+	xr[N/2 -1],xi[N/2 -1]
+   (note the exceptional second entry)
+    odd N:
+     xr[0],xr[(N-1)/2],xr[1],xi[1], ...,
+     xr[(N-1)/2 -1],xi[(N-1)/2 -1],xi[(N-1)/2]
+   (note the exceptional second & last entries)
+*/
+#define		FTPACK		112
+/* TCMPLX represents complex time traces 			*/
+#define		TCMPLX		113
+/* FAMPH represents freq domain data in amplitude/phase form	*/
+#define		FAMPH		114
+/* TAMPH represents time domain data in amplitude/phase form	*/
+#define		TAMPH		115
+/* REALPART represents the real part of a trace to Nyquist	*/
+#define		REALPART	116
+/* IMAGPART represents the real part of a trace to Nyquist	*/
+#define		IMAGPART	117
+/* AMPLITUDE represents the amplitude of a trace to Nyquist	*/
+#define		AMPLITUDE	118
+/* PHASE represents the phase of a trace to Nyquist		*/
+#define		PHASE		119
+/* KT represents wavenumber-time domain data 			*/
+#define		KT		121
+/* KOMEGA represents wavenumber-frequency domain data		*/
+#define		KOMEGA		122
+/* ENVELOPE represents the envelope of the complex time trace	*/
+#define		ENVELOPE	123
+/* INSTPHASE represents the phase of the complex time trace	*/
+#define		INSTPHASE	124
+/* INSTFREQ represents the frequency of the complex time trace	*/
+#define		INSTFREQ	125
+/* DEPTH represents traces in depth-range (z-x)			*/
+#define		TRID_DEPTH	130
+/* 3C data...  v,h1,h2=(11,12,13)+32 so a bitmask will convert  */
+/* between conventions */
+/* CHARPACK represents byte packed seismic data from supack1	*/
+#define		CHARPACK	201
+/* SHORTPACK represents 2 byte packed seismic data from supack2	*/
+#define		SHORTPACK	202
+
+
+#define ISSEISMIC(id) (( (id)==TUNK || (id)==TREAL || (id)==TDEAD || (id)==TDUMMY || (id)==TBREAK || (id)==UPHOLE || (id)==SWEEP || (id)==TIMING || (id)==WBREAK || (id)==NFGUNSIG || (id)==FFGUNSIG || (id)==SPSENSOR || (id)==TVERT || (id)==TXLIN || (id)==TINLIN || (id)==ROTVERT || (id)==TTRANS || (id)==TRADIAL || (id)==ACOR ) ? cwp_true : cwp_false ) 
+
+/* FUNCTION PROTOTYPES */
+#ifdef __cplusplus /* if C++, specify external linkage to C functions */
+extern "C" {
+#endif
+
+/* get trace and put trace */
+int fgettr(FILE *fp, segy *tp);
+int fvgettr(FILE *fp, segy *tp);
+void fputtr(FILE *fp, segy *tp);
+void fvputtr(FILE *fp, segy *tp);
+int fgettra(FILE *fp, segy *tp, int itr);
+
+/* get gather and put gather */
+segy **fget_gather(FILE *fp, cwp_String *key,cwp_String *type,Value *n_val,
+                        int *nt,int *ntr, float *dt,int *first);
+segy **get_gather(cwp_String *key, cwp_String *type, Value *n_val,
+			int *nt, int *ntr, float *dt, int *first);
+segy **fput_gather(FILE *fp, segy **rec,int *nt, int *ntr);
+segy **put_gather(segy **rec,int *nt, int *ntr);
+
+/* hdrpkge */
+void gethval(const segy *tp, int index, Value *valp);
+void puthval(segy *tp, int index, Value *valp);
+void getbhval(const bhed *bhp, int index, Value *valp);
+void putbhval(bhed *bhp, int index, Value *valp);
+void gethdval(const segy *tp, char *key, Value *valp);
+void puthdval(segy *tp, char *key, Value *valp);
+char *hdtype(const char *key);
+char *getkey(const int index);
+int getindex(const char *key);
+void swaphval(segy *tp, int index);
+void swapbhval(bhed *bhp, int index);
+void printheader(const segy *tp);
+
+void tabplot(segy *tp, int itmin, int itmax);
+
+#ifdef __cplusplus /* if C++, end external linkage specification */
+}
+#endif
+
+#endif

MDD/verbosepkg.c

+../utils/verbosepkg.c
\ No newline at end of file

MDD/wallclock_time.c

+../utils/wallclock_time.c
\ No newline at end of file

MDD/writeEigen.c

+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "segy.h"
+#include <assert.h>
+
+
+void writeEigen(char *file_out, float df, int nw_low, int nw_high, int nw, float *eigen, int nx, float dx, float xmin)
+{
+    static FILE *out_file;
+	float *trace, scl, re, im;
+	int sign, ntfft, i, j, ie, iw, count;
+	segy *hdrs_out;
+	size_t nwrite;
+    char filename[256], ext[32];
+
+    trace  = (float *)malloc(nx*sizeof(float));
+	hdrs_out  = (segy *)calloc(TRCBYTES,1);
+    
+	hdrs_out[0].dt=df*1000000;
+	hdrs_out[0].trid = 1;
+	hdrs_out[0].ns = nx;
+	hdrs_out[0].d1 = 1;
+	hdrs_out[0].f1 = 1;
+	hdrs_out[0].f2 = nw_low*df;
+	hdrs_out[0].d2 = df;
+	hdrs_out[0].trwf = nw;
+	hdrs_out[0].fldr = 1;
+
+	strcpy(filename, file_out);
+	sprintf(ext,"%s.su", "_eigen");
+	strcpy(strstr(filename, ".su"), ext);
+	out_file = fopen(filename, "w+"); assert( out_file );
+	fprintf(stderr,"writing eigenvalues of matrix to %s\n", filename);
+	count=1;
+
+	for (iw=0; iw<nw; iw++) {
+		hdrs_out[0].tracl = iw+1;
+		for (i = 0; i < nx; i++) {
+           	trace[i] = eigen[iw*nx+i];
+		}
+		nwrite = fwrite(&hdrs_out[0], 1, TRCBYTES, out_file);
+		assert( nwrite == TRCBYTES );
+		nwrite = fwrite(trace, sizeof(float), nx, out_file);
+		assert( nwrite == nx );
+    }
+    fflush(out_file);
+   	fclose(out_file);
+
+    free(hdrs_out);
+    free(trace);
+
+	return;
+}
+

Makefile

@@ -7,6 +7,7 @@ all: mkdirs
 	cd marchenko	; $(MAKE) install
 	cd corrvir		; $(MAKE) install
 	cd raytime		; $(MAKE) install
+	cd MDD			; $(MAKE) install
 
 mkdirs:
 	-mkdir -p lib
@@ -20,6 +21,7 @@ clean:
 	cd marchenko	; $(MAKE) $@
 	cd corrvir		; $(MAKE) $@
 	cd raytime		; $(MAKE) $@
+	cd MDD			; $(MAKE) $@
 
 realclean:
 	cd FFTlib 		; $(MAKE) $@
@@ -28,6 +30,7 @@ realclean:
 	cd marchenko	; $(MAKE) $@
 	cd corrvir		; $(MAKE) $@
 	cd raytime		; $(MAKE) $@
+	cd MDD			; $(MAKE) $@
 	rm -f lib/*
 	rm -f include/*
 	rm -f bin/*

utils/getFileInfo.c

-#define _FILE_OFFSET_BITS 64
-#define _LARGEFILE_SOURCE
-
 #include <assert.h>
 #include <stdio.h>
 #include <stdlib.h>