main_mpi.c

用于2维的射线追踪

/*
 * Ray2mesh : software for geophysicists.
 * Compute various scores attached to the mesh cells, based on geometric
   information that rays bring when they traverse cells.
 *
 * Copyright (C) 2003, St閜hane Genaud and Marc Grunberg
 *
 * This tool is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Library General Public
 * License as published by the Free Software Foundation; either
 * version 2 of the License, or (at your option) any later version.
 *
 * This library 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
 * Library General Public License for more details.
 *
 * You should have received a copy of the GNU Library General Public
 * License along with this library; if not, write to the Free
 * Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */
#include <stdio.h>
#include <stddef.h>		/* for offsetof */
#include "util.h"		/* get_date_stamp() */


#ifdef USE_MPILB
#include <mpilb/tag.h>
#endif

#include "r2mfile.h"
#include "merge.h"
#include "save_to_disk.h"

#define  ROOT 0			/* the MPI rank of sender */
#define MPI_MSG_DATA       1001
#define MAX_LENGTH_RANK 10

/*----------------------------------------------------------------------*/
/* send_or_recv_raydata : given a raydata_t array, give each MPI process */
/* a part this array. The root process sends a sizeperproc element part */
/* of raydata to each process, and each slave process receives this part*/
/* that it stores in the recv buffer                                    */
/*----------------------------------------------------------------------*/
struct raydata_t *
send_or_recv_raydata (struct raydata_t *raydata, int sizeperproc, int root)
{
  struct raydata_t *rbuff;	/* receive buffer */
  MPI_Datatype mpi_raydata;
  int blen[10];			/* length in bytes of each field */
  MPI_Datatype oldtypes[10];	/* type of each field */
  MPI_Aint indices[10];		/* offset in bytes of each field */


  /* every process receives data (root as well) */
#ifndef USE_MPILB
  rbuff = (struct raydata_t *)
    malloc (sizeperproc * sizeof (struct raydata_t));
#else
  rbuff = malloc (tag_spp[rank] * sizeof *rbuff);
#endif
  if (rbuff == NULL)
    {
      MPI_Abort (MPI_COMM_WORLD, 1);
      exit (1);
    }

  /* now, root and slaves cast data as MPI types for send/reveive */
  blen[0] = 1;
  blen[1] = 1;
  blen[2] = 1;
  blen[3] = 1;
  blen[4] = 1;
  blen[5] = 1;
  blen[6] = 1;
  blen[7] = RAYCODE_MAX_STRING_LENGTH + 1;
  blen[8] = 1;
  blen[9] = 1;

  oldtypes[0] = MPI_LONG;
  oldtypes[1] = MPI_DOUBLE;
  oldtypes[2] = MPI_DOUBLE;
  oldtypes[3] = MPI_DOUBLE;
  oldtypes[4] = MPI_DOUBLE;
  oldtypes[5] = MPI_DOUBLE;
  oldtypes[6] = MPI_DOUBLE;
  oldtypes[7] = MPI_CHAR;
  oldtypes[8] = MPI_DOUBLE;
  oldtypes[9] = MPI_UB;


  indices[0] = offsetof (struct raydata_t, event_id);
  indices[1] =
    offsetof (struct raydata_t, src) + offsetof (struct coord_geo_t, prof);
  indices[2] =
    offsetof (struct raydata_t, src) + offsetof (struct coord_geo_t, lon);
  indices[3] =
    offsetof (struct raydata_t, src) + offsetof (struct coord_geo_t, lat);
  indices[4] =
    offsetof (struct raydata_t, dest) + offsetof (struct coord_geo_t, prof);
  indices[5] =
    offsetof (struct raydata_t, dest) + offsetof (struct coord_geo_t, lon);
  indices[6] =
    offsetof (struct raydata_t, dest) + offsetof (struct coord_geo_t, lat);
  indices[7] = offsetof (struct raydata_t, phase);
  indices[8] = offsetof (struct raydata_t, ray_travel_time);
  indices[9] = sizeof (struct raydata_t);

  MPI_Type_struct (10, blen, indices, oldtypes, &mpi_raydata);
  MPI_Type_commit (&mpi_raydata);

#ifndef USE_MPILB
  MPI_Scatter (raydata, sizeperproc, mpi_raydata,
	       rbuff, sizeperproc, mpi_raydata, root, MPI_COMM_WORLD);
#else
  MPI_Scatterv (raydata,
		tag_spp,
		tag_spp_sum,
		mpi_raydata,
		rbuff, tag_spp[rank], mpi_raydata, root, MPI_COMM_WORLD);
#endif

  MPI_Type_free (&mpi_raydata);
  return (rbuff);
}


/*----------------------------------------------------------------------*/
/* main                                                                 */
/*----------------------------------------------------------------------*/
int
main (int argc, char **argv)
{				/* return the nb of ray computed */
  /*struct coord_geo_t source;
     struct coord_geo_t dest; */

  struct ray_config_t ray_config;
  struct ray_filter_t ray_filter;

  char *meshfile = NULL;
  char *inputdatafile = NULL;
  char *rayfilter_filename = NULL;
  char *tmpdir = NULL;

  float limit = -1;
  int nb_ray_computed = 0;
  int nb_ray_rejected = 0;
  int nb_ray_total = 0;
  int nbread = 0;
  int nberr = 0;		/* nb of line with error */

  /* res & sparse files */
  FILE *sparse_fd = NULL;
  FILE *res_fd = NULL;
  FILE *event_fd = NULL;
  FILE *filter_fd = NULL;

  /*---------------- MPI related --------------------*/
  int nbprocs;			/* number of processes in the group */
  int size;			/* offset of remainder of raydata not treated */
  int xml_buffer_size;		/* number of bytes of the xml file */
  int sizeperproc;		/* number of raydata for each process to compute */
  long int offset;

  char *xml_buffer = NULL;
  struct raydata_t *recv_data;	/* transmitted ray data buffer */
  struct raydata_t *raydata = NULL;
#ifndef USE_MPILB
  struct raydata_t *ray_remainder = NULL;	/* some rays may remain */
#endif
  double start_time;
  double comm_time;
  double *merge_times;		/* 9 counters */
  double score_time;
  double ray_time_start;
  double ray_time_end;

  char *date_stamp;

  /* catch Ctrl-C signal */
  /* do not cast 2nd arg to __sighandler_t : fatal make error wit irix cc */
  signal (SIGINT, emergency_halt);

  /* MPI initialization */
  MPI_Init (&argc, &argv);
  MPI_Comm_rank (MPI_COMM_WORLD, &rank);
  MPI_Comm_size (MPI_COMM_WORLD, &nbprocs);

  /* cmd line */
  memset (&ray_filter, 0, sizeof (struct ray_filter_t));
  parse_command_line (argc, argv,
		      &ray_config, &ray_filter,
		      &meshfile,
		      &celldatafilename,
		      &inputdatafile,
		      &rayfilter_filename, &limit, &output_format, &tmpdir);
#ifdef DEBUG
  if (rank == ROOT)
    {
      fprintf (stderr, "%s:main():%d  command line parsing done.\n",
	       __FILE__, __LINE__);
    }
#endif

  if (limit > 0 && rank == ROOT)
    {
      if (strchr (output_format, 's'))
	{
	  fprintf (stderr,
		   "[ROOT] can not use *sco* output format with -L option\n");
	  MPI_Abort (MPI_COMM_WORLD, 1);
	  exit (1);
	}
      fprintf (stderr, "[ROOT] info, limit memory usage to %.1fM\n", limit);
    }

  /* checks if the basename for files is provided */
  if (!celldatafilename)
    {
      fprintf (stderr, "No output file specified. Exiting.\n");
      MPI_Abort (MPI_COMM_WORLD, 1);
      exit (1);
    }

  start_time = MPI_Wtime ();

    /*************************************************************/
  /* first check that files are accessible (on root side only) */
    /*************************************************************/
  if (rank == ROOT)
    {
	/********/
      /* ROOT */
	/********/
      fprintf (stdout, "*** [root] MPI process started (MPI_Init).\n");
      fprintf (stdout, "*** [root] %s running on %d processors.\n",
	       PACKAGE, nbprocs);
      ray2mesh_info (&ray_config);
      check_file_access (meshfile);
      check_file_access (inputdatafile);
#ifdef RAYTRACING_ONLY
      fprintf (stdout,
	       "%s : only ray-tracing executed (compiled with -DRAYTRACING_ONLY)\n",
	       PACKAGE);
#endif
#ifdef DEBUG
      fprintf (stderr, "%s:main():%d [root] after check access .\n",
	       __FILE__, __LINE__);
#endif
      /* init of the mesh->parameters struct  */
      if ((xml_buffer =
	   load_file_to_memory (meshfile, &xml_buffer_size)) == NULL)
	{
	  fprintf (stderr,
		   "%s : could not load %s to memory (%d bytes). Exiting.\n",
		   PACKAGE, meshfile, xml_buffer_size);
	  MPI_Abort (MPI_COMM_WORLD, 1);
	  exit (1);
	}

      /* open the file containing the rays */
      if ((fdinput = fopen (inputdatafile, "r")) == NULL)
	{
	  fprintf (stderr,
		   "Can not open input datafile '%s' for reading. Exiting.\n",
		   inputdatafile);
	  MPI_Abort (MPI_COMM_WORLD, 1);
	  exit (1);
	}
    }
  else
    {
	/**********/
      /* SLAVES */
	/**********/
      date_stamp = get_date_stamp ();
      fprintf (stdout,
	       "*** [process %d] MPI process started (MPI_Init) as process %d (%s)\n",
	       rank, getpid (), date_stamp);
      free (date_stamp);

      /* open the filter file if needed */
      if (rayfilter_filename)
	{
	  char *filename;
	  char *ext = "sei";

	  filename =
	    (char *) malloc (sizeof (char) *
			     (strlen (rayfilter_filename) + strlen ("-p") +
			      RANK_SIZE + strlen (".") + strlen ("sei") + 1));
	  assert (filename);
	  sprintf (filename, "%s-p%.3d.%s", rayfilter_filename, rank, ext);

	  if ((filter_fd = fopen (filename, "w")) == NULL)
	    {
	      fprintf (stderr, "Cannot open %s for writing ray filtered.\n",
		       filename);
	      MPI_Abort (MPI_COMM_WORLD, 1);
	      exit (1);
	    }
	  free (filename);
	}
    }


    /**************/
  /* XML buffer */
    /**************/