parallel.c

蒙特卡罗模拟光子成像C语言版,代码简洁专业

/* Support code for parallel execution (and only for parallel execution) */

/*
 * This file is part of tMCimg.
 * 
 * tMCimg is free software; you can redistribute it and/or modify it under
 * the terms of the GNU General Public License as published by the Free
 * Software Foundation; either version 2 of the License, or (at your option)
 * any later version.
 * 
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
 * or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
 * for more details.
 * 
 * You should have received a copy of the GNU General Public License along
 * with this program; if not, write to the Free Software Foundation, Inc.,
 * 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include "config.h"
#include "mccore.h"

#if (HAVE_ASSERT_H)
#include <assert.h>
#endif

/* Defines specific to the parallel code */
#include "parallel.h"

static int   rank = -1;
static int   nnodes;

char *mpibuff = NULL;
int   mpiblen = 0;

void initParallel(int *argc, char ***argv)
{
  MPI_Init(argc, argv);
  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
  MPI_Comm_size(MPI_COMM_WORLD, &nnodes);

  return;
}

int isMaster(void)
{
  if (rank == -1)
    {
      fprintf(stderr, "Error, isMaster() called before MPI initialized\n");
      exit(1);
    }

  return (rank == 0);
}

int myNode(void)
{
  return rank;
}

int getNodeCount(void)
{
  return nnodes;
}

/* This is just savehis() with the file write replaced by MPI_Send() */

void parallel_savehis(const struct Config *system,
                      double  x, double  y, double  z,
                      double kx, double ky, double kz, double T, float *len)
{
  static float *buff = NULL;
  static int    bufl = 0;
  int i, tindex;

  /* Allocate storage */

  if ((mpibuff == NULL) || (mpiblen != 4 + system->Ntissue))
    {
      if (mpibuff != NULL)
        free_1d_array(mpibuff, 100 * mpiblen);

      mpiblen = 4 + system->Ntissue;

      /* Allocate space for 100 exiting photons */
      mpibuff = (char *)alloc_1d_array(100 * mpiblen, sizeof(float));

#if (PARANOID)
      assert(mpibuff != NULL);
#endif

      MPI_Buffer_attach(mpibuff, 100 * mpiblen * sizeof(float));
    }

  if ((buff == NULL) || (bufl != 4 + system->Ntissue))
    {
      if (buff != NULL)
        free_1d_array(buff, 4 + system->Ntissue);

      bufl = 4 + system->Ntissue;
      buff = alloc_1d_array(bufl, sizeof(float));

#if (PARANOID)
      assert(buff != NULL);
#endif
    }

  /* Current time gate */

  tindex = (int)((T - system->Tmin) / system->stepT);

  /* Rescale back to physical units */

  for (i = 1; i <= system->Ntissue; i++)
    len[i] = unscale_length(len[i], system);

  /* IF NO DETECTORS DEFINED THEN SAVE EXIT POSITION 
   * ELSE LOOP THROUGH NUMBER OF DETECTORS */

  if (system->nDets == 0)
    {
      /* Convert back to normal dimensions [mm] before saving */

      buff[0] = unscale_length(x, system);
      buff[1] = unscale_length(y, system);
      buff[2] = unscale_length(z, system);

      if (system->flags.exact_exit_time != 0)
	buff[3] = (float)T;
      else
	buff[3] = (float)tindex;

      /* SEND TO MASTER NODE */

      memcpy(&buff[4], &len[1], system->Ntissue * sizeof(float));

      /* Asynchronous I/O would probably be faster ... */
      MPI_Send(buff, 4+system->Ntissue, MPI_FLOAT, 
               0, TAG_HISFILE, MPI_COMM_WORLD);
    }
  else				  /* system->nDets > 0 */
    {
#if (PARANOID)
      int found_det = 0;
#endif

      for (i = 0; i < system->nDets; i++)
	{
	  /* dR and system->detRad[] are both in voxel units */

	  double dR = SQR(x - system->detLoc[i][0]) +
		      SQR(y - system->detLoc[i][1]) + 
		      SQR(z - system->detLoc[i][2]);

	  if (dR <= SQR(system->detRad[i]))
	    {
	      buff[0] = (float)i;	/* Detector #  */

#if (PARANOID)
              /* Other half of double-counting check */
              assert(found_det == 0);
#endif

	      /* (kx,ky,kz) is implicitly a unit vector */

#if (HAVE_ASSERT)
	      assert((SQR(kx) + SQR(ky) + SQR(kz) > 0.9999) &&
		     (SQR(kx) + SQR(ky) + SQR(kz) < 1.0001));
#endif

	      /* Accept/reject based on detector NA */

	      /* Accept both NA=0 and NA=1 to specify all photons
	       *  should be accepted (i.e. effectively disable the check) */

	      if ((system->detNA[i] > 0) && (system->detNA[i] < 1))
		{
                  /* Not implemented yet */
                  MPI_Abort(MPI_COMM_WORLD, MPI_ERR_OTHER);
		}

	      /* use either tgate # or exact photon exit time... */

	      if (system->flags.exact_exit_time != 0)
		buff[1] = (float)T;
	      else
		buff[1] = (float)tindex;

              /* SEND TO MASTER NODE */

              memcpy(&buff[2], &len[1], system->Ntissue * sizeof(float));

              MPI_Send(buff, 2+system->Ntissue, MPI_FLOAT, 
                        0, TAG_HISFILE, MPI_COMM_WORLD);

#if (PARANOID)
              /* Check for double-counting */
              found_det += 1;
#endif
	    }
	}
    }

  return;
}

void endParallel(void)
{
  /* Remove buffer if it's still attached, deallocate memory */

  if (mpibuff != NULL)
    {
      void *p;
      int s;

      MPI_Buffer_detach(&p, &s);

      free_1d_array(mpibuff, 100 * mpiblen);

      mpibuff = NULL;
      mpiblen = 0;
    }

  /* Close down MPI connection */

  MPI_Finalize();

  rank   = -1;
  nnodes =  0;

  return;
}