octave_packages/openmpi_ext-1.0.2/montecarlo.m

   1 # Copyright (C) 2006, 2009  Michael Creel <michael.creel@uab.es>
   2 # under the terms of the GNU General Public License.
   3 #
   4 # This program is free software; you can redistribute it and/or modify
   5 # it under the terms of the GNU General Public License as published by
   6 # the Free Software Foundation; either version 2 of the License, or
   7 # (at your option) any later version.
   8 #
   9 # This program is distributed in the hope that it will be useful,
  10 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  11 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12 # GNU General Public License for more details.
  13 #
  14 # You should have received a copy of the GNU General Public License
  15 # along with this program; if not, write to the Free Software
  16 # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  17
  18 # montecarlo.m: generates a specified number of replications of a function's
  19 # output and writes them to a user-specified output file.
  20 #
  21 # USAGE: montecarlo(f,f_args,reps,outfile,n_pooled,n_returns,usempi, verbose)
  22 #
  23 # IMPORTANT: f should return a row vector of output from feval(f,f_args)
  24 #
  25 # For normal evaluation on one core, only the first 4 arguments are required.
  26 # * Arg 1: (required) the function that generates a row vector of output
  27 # * Arg 2: (required) the arguments of the function, in a cell
  28 # * Arg 3: (required) the number of replications to generate
  29 # * Arg 4: (required) the output file name
  30 # * Arg 5 (optional) number of replications to be pooled together between writes
  31 # * Arg 6 (optional) verbose: 1 for on, 0 for off
  32 #
  33 # If using MPI, you should run using ranks equal to number of cores plus 1,
  34 # and should make sure that the core running the frontend is also the one that
  35 # has the second rank. That way the core the frontend is on will also do work.
  36
  37 function n_received = montecarlo(f,f_args,reps,outfile,n_pooled,verbose)
  38
  39         t0 = clock(); # initialize timing
  40
  41         # defaults for optional arguments
  42         if (nargin < 6) verbose = false; endif
  43         if (nargin < 5) n_pooled = 1; endif;
  44
  45         if MPI_Initialized      # check if doing this parallel or serial
  46                 use_mpi = true;
  47                 CW = MPI_Comm_Load("NEWORLD");
  48                 is_node = MPI_Comm_rank(CW);
  49                 nodes = MPI_Comm_size(CW);
  50                 mytag = 48;
  51         else
  52                 use_mpi = false;
  53                 is_node = 0;
  54         endif
  55
  56         if is_node # compute nodes
  57                 more_please = 1;
  58                 while more_please
  59                         for i = 1:n_pooled
  60                                 contrib = feval(f, f_args);
  61                                 contribs(i,:) = contrib;
  62                         endfor
  63                         MPI_Send(contribs, 0, mytag, CW);
  64                         pause(0.05); # give time for the fronted to send a stop message, if done
  65                         # check if we're done
  66                         if (MPI_Iprobe(0, is_node, CW)) # check for ping from rank 0
  67                                 junk = MPI_Recv(0, is_node, CW);
  68                                 break;
  69                         endif
  70                 endwhile
  71         else # frontend
  72                 received = 0;
  73                 done = false;
  74                 while received < reps
  75                         if use_mpi
  76                                 # retrieve results from compute nodes
  77                                 for i = 1:nodes-1
  78                                         # compute nodes have results yet?
  79                                         ready = false;
  80                                         ready = MPI_Iprobe(i, mytag, CW); # check if message pending
  81                                         if ready
  82                                                 # get it if it's there
  83                                                 contribs = MPI_Recv(i, mytag, CW);
  84                                                 need = reps - received;
  85                                                 received = received + n_pooled;
  86                                                 # truncate?
  87                                                 if n_pooled  >= need
  88                                                                 contribs = contribs(1:need,:);
  89                                                                 done = true;
  90                                                 endif
  91                                                 # write to output file
  92                                                 FN = fopen (outfile, "a");
  93                                                 if (FN < 0) error ("montecarlo: couldn't open output file %s", outfile); endif
  94                                                 t = etime(clock(), t0);
  95                                                 for j = 1:rows(contribs)
  96                                                         fprintf(FN, "%f ", i, t, contribs(j,:));
  97                                                         fprintf(FN, "\n");
  98                                                 endfor
  99                                                 fclose(FN);
 100                                                 if verbose printf("\nContribution received from node%d.  Received so far: %d\n", i, received); endif
 101                                                 if done
 102                                                         # tell compute nodes to stop loop
 103                                                         for j = 1:5
 104                                                                 for i = 1:(nodes-1)
 105                                                                         if (j==1) MPI_Send(" ",i,i,CW); endif # send out message to stop
 106                                                                         ready = MPI_Iprobe(i, mytag, CW); # get last messages
 107                                                                         if ready contribs = MPI_Recv(i, mytag, CW); endif
 108                                                                 endfor
 109                                                         endfor
 110                                                         break;
 111                                                 endif
 112                                         endif
 113                                 endfor
 114                         else
 115                                 for i = 1:n_pooled
 116                                         contrib = feval(f, f_args);
 117                                         contribs(i,:) = contrib;
 118                                 endfor
 119                                 need = reps - received;
 120                                 received = received + n_pooled;
 121                                 # truncate?
 122                                 if n_pooled  >= need
 123                                         contribs = contribs(1:need,:);
 124                                 endif
 125                                 # write to output file
 126                                 FN = fopen (outfile, "a");
 127                                 if (FN < 0) error ("montecarlo: couldn't open output file %s", outfile); endif
 128                                 t = etime(clock(), t0);
 129                                 for j = 1:rows(contribs)
 130                                         fprintf(FN, "%f ", 0, t, contribs(j,:));
 131                                         fprintf(FN, "\n");
 132                                 endfor
 133                                 fclose(FN);
 134                                 if verbose printf("\nContribution received from node 0.  Received so far: %d\n", received); endif
 135                         endif
 136                 endwhile
 137         endif
 138 endfunction