sudrefsvmaster.c

su 的源代码库

   }
   initF = NINT(info->f1 / info->dF);
   lastF = NINT(info->f2 / info->dF);  

   /* attenuation of wrap-around */
   info->tau = log(info->tau) / tMax;
   if (info->tau > TAUMAX)
      info->tau = TAUMAX;
      
   if (info->verbose)
      fprintf(stderr, "Discrete frequency range to model: [%d, %d]\n", 
	      initF, lastF);
   
   if (nFreqProc == 0)
   {
      nFreqProc = NINT((float) info->nF / (float) nProc + .5);
      if (nFreqProc > info->nF) nFreqProc = info->nF;
   }
   else
      while (nFreqProc > info->nF) nFreqProc /= 2;
   nFreqPart = NINT((float) info->nF / (float) nFreqProc + .5);

   /* memory allocation for frequency arrays */
   uRF = alloc2complex(info->nSamples / 2 + 1, info->nR);
   uZF = alloc2complex(info->nSamples / 2 + 1, info->nR);
   freqPart = alloc2complex(nFreqProc, info->nR);
   statusFreq = alloc2int(3, nFreqPart);

   /* defining frequency partitions */
   for (k = initF, i = 0; i < nFreqPart; i++, k += nFreqProc)
   {
      statusFreq[i][0] = k;
      statusFreq[i][1] = MIN(k + nFreqProc - 1, lastF);
      statusFreq[i][2] = 0;       
   }

   if (info->verbose)
      fprintf(stderr, "Starting communication with PVM\n");
   
   /* starting communication with PVM */
   if ((apl_pid = pvm_mytid()) < 0) 
   {
      pvm_perror("Error enrolling master process");
      exit(-1);
   } 
   processControl = CreateSlaves(processes, PROCESS, nProc);
   if (processControl != nProc)
   {
      fprintf(stderr,"Problem starting PVM daemons\n");
      exit(-1);
   }

   info->nFreqProc = nFreqProc;
   /* Broadcasting all processes common information */
   BroadINFO(info, 1, processes, nProc, GENERAL_INFORMATION);
   
   if (info->verbose)
      fprintf(stderr, "Broadcasting model information to all slaves\n");

   /* sending all profiles */
   BroadFloat(thick, info->nL + 1, processes, nProc, THICKNESS);
   BroadFloat(rho, info->nL + 1, processes, nProc, DENSITY);
   BroadFloat(alpha, info->nL + 1, processes, nProc, ALPHA);
   BroadFloat(qP, info->nL + 1, processes, nProc, QALPHA);
   BroadFloat(beta, info->nL + 1, processes, nProc, BETA);
   BroadFloat(qS, info->nL + 1, processes, nProc, QBETA);

   /* freeing memory */
   free1float(thick);
   free1float(rho);
   free1float(alpha);
   free1float(qP);
   free1float(beta);
   free1float(qS);

   /* sending frequency partitions for each process */
   for (iProc = 0; iProc < nProc; iProc++)
   {
      FInfo[0] = statusFreq[iProc][0];
      FInfo[1] = statusFreq[iProc][1];

      if (info->verbose)
	 fprintf(stderr, 
	 "Master sending frequencies [%d, %d] out of %d to slave %d [id:%d]\n"
	  ,FInfo[0], FInfo[1], info->nF, iProc, processes[iProc]);

      procInfo[iProc][0] = FInfo[0]; procInfo[iProc][1] = FInfo[1];
      SendInt(FInfo, 2, processes[iProc], FREQUENCY_LIMITS);
      statusFreq[iProc][2] = 1;
   }

   /* waiting modelled frequencies */
   /* master process will send more frequencies if there's more work to do */
   /* measuring elapsed time */
   wallcpu = walltime();  
   
   /* reseting frequency counter */
   FReceived = 0;
   
   while (FOREVER)
   {
      pid = RecvCplx(freqPart[0], info->nR * nFreqProc, -1, 
		     FREQUENCY_PARTITION_VERTICAL);

      /* finding the frequency limits of this process */
      iProc = 0;
      while (pid != processes[iProc])
	 iProc++;

      /* copying into proper place of the total frequency array */
      for (iR = 0; iR < info->nR; iR++)
      {
	 for (k = 0, i = procInfo[iProc][0]; i <= procInfo[iProc][1]; i++, k++)
	 {
	    uZF[iR][i] = freqPart[iR][k];
	 }
      }

      pid = RecvCplx(freqPart[0], info->nR * nFreqProc, -1, 
		     FREQUENCY_PARTITION_RADIAL);
      
      /* finding the frequency limits of this process */
      iProc = 0;
      while (pid != processes[iProc])
	 iProc++;
   
      /* copying into proper place of the total frequency array */
      for (iR = 0; iR < info->nR; iR++)
      { 
	 for (k = 0, i = procInfo[iProc][0]; i <= procInfo[iProc][1]; i++, k++)
	 {
	    uRF[iR][i] = freqPart[iR][k];
	 }
      }

      /* summing frequencies that are done */
      FReceived += procInfo[iProc][1] - procInfo[iProc][0] + 1;

      if (info->verbose)
	 fprintf(stderr, "Master received %d frequencies, remaining %d\n", 
	      FReceived, info->nF - FReceived);

/*       if (FReceived >= info->nF) break; */

      /* defining new frequency limits */
      i = 0;
      while (i < nFreqPart && statusFreq[i][2])
	 i++;
      
      if (i < nFreqPart)
      {
	 /* there is still more work to be done */
	 /* tell this process to not die */
	 die = 0;
	 SendInt(&die, 1, processes[iProc], DIE);

	 FInfo[0] = statusFreq[i][0];
	 FInfo[1] = statusFreq[i][1];

	 if (info->verbose)
	    fprintf(stderr, 
		    "Master sending frequencies [%d, %d] to slave %d\n", 
		    FInfo[0], FInfo[1], processes[iProc]);
	 
	 procInfo[iProc][0] = FInfo[0]; procInfo[iProc][1] = FInfo[1];
	 SendInt(FInfo, 2, processes[iProc], FREQUENCY_LIMITS);
	 statusFreq[i][2] = 1;
      }
      else
      {
	 /* tell this process to die since there is no more work to do */
	 if (info->verbose)
	    fprintf(stderr, "Master ''killing'' slave %d\n", processes[iProc]);
	 die = 1;
	 SendInt(&die, 1, processes[iProc], DIE);
      }
      
      /* a check to get out the loop */
      if (FReceived >= info->nF) break; 
   }

   if (info->verbose)
      fprintf(stderr, "Master ''killing'' remaining slaves\n");

   /* getting elapsed time */
   wallcpu = walltime() - wallcpu;
   fprintf(stderr, "Wall clock time = %f seconds\n", wallcpu);  
   
   /* going to time domain */
   memset( (void *) &trZ, (int) '\0', sizeof(trZ));     
   memset( (void *) &trR, (int) '\0', sizeof(trR));     
   trZ.dt = dt * 1000000;
   trZ.ns = nSamplesOrig;
   trR.dt = dt * 1000000;
   trR.ns = nSamplesOrig;
   
   /* z component */
   for (iR = 0; iR < info->nR; iR++)
   {
      trZ.tracl = iR + 1;
      /* inverse FFT */
      pfacr(1, info->nSamples, uZF[iR], trZ.data); 
      for (i = 0; i < info->nSamples; i++)
      {
	 /* compensating for the complex frequency */
	 trZ.data[i] *= exp(info->tau * i * dt);
      }
      puttr(&trZ);
   }

   /* r component */
   for (iR = 0; iR < info->nR; iR++)
   {
      trR.tracl = iR + 1;
      /* inverse FFT */
      pfacr(1, info->nSamples, uRF[iR], trR.data); 
      for (i = 0; i < info->nSamples; i++)
      {
	 /* compensating for the complex frequency */
	 trR.data[i] *= exp(info->tau * i * dt);
      }
      puttr(&trR);
   }
}   
#include <time.h>

float walltime()
/*****************************************************************************
return elapsed time (wall clock time) in seconds using ANSI C built-ins

Notes:
	return value will be an integral number of seconds since t1 and t2,
	as returned by the time() intrinsic, are the number of seconds
	since the epoch
******************************************************************************
Author:		Jack K. Cohen, Colorado School of Mines, 07/27/90
*****************************************************************************/
{
	static int firsttime = 1;	/* First entry?                 */
	static time_t lasttime;		/* Save return for next entry   */
	time_t t1, t2;			/* Resp. last and current times */

	if (firsttime) {
		firsttime = 0;
		lasttime = time(&t1);
		return 0.0;
	} else {
		t1 = lasttime;
		lasttime = time(&t2);
		return (float) difftime(t2, t1);
	}
}