2dimaging.mpi.c

主从模式粗粒级并行算法C程序：这是我以前研究生期间编写的叠前地震成像C源码

          fclose(fp);
        }
        fclose(iwfp);
        MPI_Barrier( MPI_COMM_WORLD);

        fp=fopen(tmpfile, "rb");

        ncig=nmx/ncdpstep_cig;
        if(nmx%ncdpstep_cig==0&&ncdpstep_cig!=1) ncig=nmx/ncdpstep_cig+1;

	if(myid==0) {
	   printf("ncig=%d \n",ncig);
	   printf("ncdpstep_cig=%d \n",ncdpstep_cig);
	}

        vz=alloc1float(nz);
        dss=alloc3float(nhx, nmx, nz);
        dsg=alloc3float(nhx, nmx, nz);

        if(myid==0) memory_used(nxv, nmx, nhx, nz, ncig, operator);

        zero1float(vz, nz);
        zero3float(dss, nhx, nmx, nz);
        zero3float(dsg, nhx, nmx, nz);

        velocity(vv, vz, dss, dsg, nxv, nzv, nmx, nhx, nz, xmmin, offsetmin,
                 xv0, dxv, dzv, dz, dmx, dhx);

        free2float(vv);

        wave=alloc2complex(nhx, nmx);
        image=alloc2float(nz, nmx);
        cimage=alloc3float(nz, nhx, ncig);
        taperx=alloc1float(nmx);
        taperh=alloc1float(nhx);

        zero2complex(wave, nhx, nmx);
        zero2float(image, nz, nmx);
        zero3float(cimage, nz, nhx, ncig);
        zero1float(taperx, nmx);
        zero1float(taperh, nhx);

	if(myid==0) {
          printf("***************************************************\n");
          printf("         Wavefield Extrapolating & Imaging         \n");
          printf("***************************************************\n");
	}

        tapfunc(taperx, nmx, 15, 15, 2);
        tapfunc(taperh, nhx, 6, 6, +1);

	if(np==1){
           //for(iw=myid;iw<nw;iw=iw+np){
           for(iw=nw/3+myid;iw<nw/3+1;iw=iw+np){
              if(myid==0) printf("Myid= %d _______ iw/nw= %d/%d\n", myid,iw,nw);
             
              fseek(fp, (long)(sizeof(complex)*nmx*nhx*iw), 0);
              for(imx=0;imx<nmx;imx++) 
                 fread(wave[imx],sizeof(complex),nhx,fp);

              w=(nf1+iw)*dw;
              presdm(wave, cimage, image, dss, dsg, vz, taperx, taperh,
                  nmx, nhx, nz, dmx, dhx, dz, w, operator, xmmin,
		  offsetmin, ncig, ncdpstep_cig, myid); 
           }
	 }else{ /* np>1 */
	   if(myid==0) {
	      num=1;

	      nnp=np-1<nw?np-1:nw;

	      for(i=1;i<=nnp;i++){
	        iw=num;
	        MPI_Send(&iw,1,MPI_INT,i,i,MPI_COMM_WORLD);
                num++;
	      }
	      for(i=1;i<=nw;i++){
                MPI_Recv(&iww,1,MPI_INT,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&status);
		mpinode = status.MPI_SOURCE;
                printf("Myid= %d ___ iw/nw= %d/%d finished !\n", mpinode,iww,nw);
		if(num<=nw){
		   iw=num;
		   MPI_Send(&iw, 1, MPI_INT, mpinode, iw, MPI_COMM_WORLD); 
                   num++;
		}else
		   MPI_Send(&unit, 0, MPI_INT, mpinode, 0, MPI_COMM_WORLD);
	      }/* i loop */

	   } else{
loop:
             MPI_Recv(&iw, 1, MPI_INT, 0, MPI_ANY_TAG, MPI_COMM_WORLD, &status);
             ip=status.MPI_TAG;
	     if(ip!=0) {
             
              fseek(fp, (long)(sizeof(complex)*nmx*nhx*iw), 0);
              for(imx=0;imx<nmx;imx++) 
                 fread(wave[imx],sizeof(complex),nhx,fp);

              w=(nf1+iw)*dw;
              presdm(wave, cimage, image, dss, dsg, vz, taperx, taperh,
                  nmx, nhx, nz, dmx, dhx, dz, w, operator, xmmin,
		  offsetmin, ncig, ncdpstep_cig, myid); 

	      MPI_Send(&iw,1,MPI_INT,0,ip,MPI_COMM_WORLD);
	      goto loop;

	     }/* ip!=0 */
	   }/* myid==0 */

	 }/* np>1 */

        MPI_Barrier(MPI_COMM_WORLD);

        free2complex(wave);
        free1float(taperx);
        free1float(taperh);
        free1float(vz);
        free3float(dss);
        free3float(dsg);

	if(myid==0)printf("---------- Stack_ImageCIGs_allnodes ----------\n");
        Stack_ImageCIGs_allnodes(owfile1, owfile2, cimage, image, seismic, 
             nmx, nhx, nz, ncdp0, offsetmin, dhx, xmmin, dmx,
	     dz, ncig, ncdpstep_cig, myid);

        free2float(image);
        free3float(cimage);
        fclose(fp);

        MPI_Barrier(MPI_COMM_WORLD);
      
        if(myid==0){
	  fgama=gamamin;
	  dgama=(gamamax-gamamin)/(ngama-1);

	  printf("---------- odcig2adcig_gama  ----------\n");
          odcig2adcig_gama(owfile2, owfile3, seismic, nmx, nhx, ngama, nz, ncdp0,
                      ncig, ncdpstep_cig, xmmin, dmx, offsetmin, dhx, fgama, dgama, dz);
	}

        MPI_Finalize();
}

void presdm(complex **wave, float ***cimage, float **image,  
     float ***dss, float ***dsg, float *vz, float *taperx, float *taperh,
     int nmx, int nhx, int nz, float dmx, float dhx, float dz, float w, 
     int operator, float xmmin, float offsetmin, 
     int ncig, int ncdpstep_cig, int myid)
{
     int    iz;

     for(iz=0;iz<nz;iz++){
       if(operator==1){   /* Split-step Fourier */
          ssf(wave, cimage, image, dss, dsg, vz, taperx, taperh,
          nmx, nhx, iz, dmx, dhx, dz, w, xmmin, offsetmin, 
	  ncig, ncdpstep_cig, myid); 
       }
       if(operator==2){   /* Quasi-linear Fourier/Born       */
          LQELBF(wave, cimage, image, dss, dsg, vz, taperx, taperh, 
          nmx, nhx, iz, dmx, dhx, dz, w, xmmin, offsetmin, 
	  ncig, ncdpstep_cig, myid);
       }
    }
}
void filter(complex *cq,int nw,int nf1,int nf2,int nf3,int nf4)
{
    int iw;
    float tmpp=0.0;
    for (iw=0; iw<nw; iw++)
    {
       if(iw>=nf1&&iw<=nf2)
       {
          tmpp=0.54+0.46*cos(PI*(iw-nf1)/(nf2-nf1)-PI);
          cq[iw]=crmul(cq[iw],tmpp);
       }
       else if(iw>=nf3&&iw<=nf4)
       {
          tmpp=0.54+0.46*cos(PI*(nf3-iw)/(nf4-nf3));
          cq[iw]=crmul(cq[iw],tmpp);
       }
       else if(iw<nf1||iw>nf4)
          cq[iw]=cmplx(0.0, 0.0);
    }
}

void myid_node_delay(int myid)
{
     int   i, j;
      
     for(i=0;i<myid*500000000;i++)  j=100*100*100;
}

/********************************************************************/
/*     calculating taper coeficient                                 */
/********************************************************************/
void tapfunc(float *taper, int nx, int nxleft, int nxright, int iflag)
{
     /*****************************************/ 
     /*  iflag =-1 (left-sided taper)          */
     /*  iflag =+1 (right-sided taper)          */
     /*  iflag = 2 (two-sided taper)          */
     /*****************************************/ 

     int   ntpleft, ntpright;
     int   kk, kk1;
     float aa;

     ntpleft = nxleft - 1;
     ntpright = nxright - 1;

     for(kk=0;kk<nx;kk++)
        taper[kk]=1.0;

     for(kk1=0;kk1<nx;kk1++){
       if(iflag==2||iflag==-1){
          if(kk1<ntpleft-1&&ntpleft>3){
            aa=0.5-0.5*cos(PI*kk1/ntpleft);
            taper[kk1]=aa;
          }
       }
       if(iflag==2||iflag==+1){
          if(kk1>nx-nxright&&ntpright>3){
            aa=0.5-0.5*cos(PI*(nx-kk1)/ntpright);
            taper[kk1]=aa;
          }
       }
     }
}
void memory_used(int nxv, int nmx, int nhx, int nz, int ncig, int operator)
{
     int mem=0, mem0;

         mem+=nxv*nz*FSIZE;
         mem+=nz*FSIZE;
         mem+=nmx*nhx*nz*2*FSIZE;
         mem+=nmx*nhx*sizeof(complex); 
         mem+=nmx*nz*FSIZE;
         mem+=ncig*nhx*nz*FSIZE;

         if(operator==1)
           mem+=nmx*nhx*sizeof(complex); 
         else
           mem+=nmx*nhx*sizeof(complex)*3; 

         mem0=mem/1000000+5;

         printf("\n");
         printf("@@@@: Running this program, need memory >%d M-byte\n",mem0);
         printf("\n");
}
/*****************************************************************************
XINDEX - determine index of x with respect to an array of x values

xindex		determine index of x with respect to an array of x values

******************************************************************************
Input:
nx		number of x values in array ax
ax		array[nx] of monotonically increasing or decreasing x values
x		the value for which index is to be determined
index		index determined previously (used to begin search)

Output:
index		for monotonically increasing ax values, the largest index
		for which ax[index]<=x, except index=0 if ax[0]>x;
		for monotonically decreasing ax values, the largest index
		for which ax[index]>=x, except index=0 if ax[0]<x

******************************************************************************
Notes:
This function is designed to be particularly efficient when called
repeatedly for slightly changing x values; in such cases, the index 
returned from one call should be used in the next.

******************************************************************************
Author:  Dave Hale, Colorado School of Mines, 12/25/89
*****************************************************************************/
/**************** end self doc ********************************/

void xindex (int nx, float ax[], float x, int *index)
/*****************************************************************************
determine index of x with respect to an array of x values
******************************************************************************
Input:
nx		number of x values in array ax
ax		array[nx] of monotonically increasing or decreasing x values
x		the value for which index is to be determined
index		index determined previously (used to begin search)

Output:
index		for monotonically increasing ax values, the largest index
		for which ax[index]<=x, except index=0 if ax[0]>x;
		for monotonically decreasing ax values, the largest index
		for which ax[index]>=x, except index=0 if ax[0]<x
******************************************************************************
Notes:
This function is designed to be particularly efficient when called
repeatedly for slightly changing x values; in such cases, the index 
returned from one call should be used in the next.
******************************************************************************
Author:  Dave Hale, Colorado School of Mines, 12/25/89
*****************************************************************************/
{
	int lower,upper,middle,step;

	/* initialize lower and upper indices and step */
	lower = *index;
	if (lower<0) lower = 0;
	if (lower>=nx) lower = nx-1;
	upper = lower+1;