/* * 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                                                                 *//*----------------------------------------------------------------------*/intmain (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 */    /**************/