--- /dev/null
+# Copyright (C) 2006, 2009 Michael Creel <michael.creel@uab.es>
+# under the terms of the GNU General Public License.
+#
+# This program 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
+
+# montecarlo.m: generates a specified number of replications of a function's
+# output and writes them to a user-specified output file.
+#
+# USAGE: montecarlo(f,f_args,reps,outfile,n_pooled,n_returns,usempi, verbose)
+#
+# IMPORTANT: f should return a row vector of output from feval(f,f_args)
+#
+# For normal evaluation on one core, only the first 4 arguments are required.
+# * Arg 1: (required) the function that generates a row vector of output
+# * Arg 2: (required) the arguments of the function, in a cell
+# * Arg 3: (required) the number of replications to generate
+# * Arg 4: (required) the output file name
+# * Arg 5 (optional) number of replications to be pooled together between writes
+# * Arg 6 (optional) verbose: 1 for on, 0 for off
+#
+# If using MPI, you should run using ranks equal to number of cores plus 1,
+# and should make sure that the core running the frontend is also the one that
+# has the second rank. That way the core the frontend is on will also do work.
+
+function n_received = montecarlo(f,f_args,reps,outfile,n_pooled,verbose)
+
+ t0 = clock(); # initialize timing
+
+ # defaults for optional arguments
+ if (nargin < 6) verbose = false; endif
+ if (nargin < 5) n_pooled = 1; endif;
+
+ if MPI_Initialized # check if doing this parallel or serial
+ use_mpi = true;
+ CW = MPI_Comm_Load("NEWORLD");
+ is_node = MPI_Comm_rank(CW);
+ nodes = MPI_Comm_size(CW);
+ mytag = 48;
+ else
+ use_mpi = false;
+ is_node = 0;
+ endif
+
+ if is_node # compute nodes
+ more_please = 1;
+ while more_please
+ for i = 1:n_pooled
+ contrib = feval(f, f_args);
+ contribs(i,:) = contrib;
+ endfor
+ MPI_Send(contribs, 0, mytag, CW);
+ pause(0.05); # give time for the fronted to send a stop message, if done
+ # check if we're done
+ if (MPI_Iprobe(0, is_node, CW)) # check for ping from rank 0
+ junk = MPI_Recv(0, is_node, CW);
+ break;
+ endif
+ endwhile
+ else # frontend
+ received = 0;
+ done = false;
+ while received < reps
+ if use_mpi
+ # retrieve results from compute nodes
+ for i = 1:nodes-1
+ # compute nodes have results yet?
+ ready = false;
+ ready = MPI_Iprobe(i, mytag, CW); # check if message pending
+ if ready
+ # get it if it's there
+ contribs = MPI_Recv(i, mytag, CW);
+ need = reps - received;
+ received = received + n_pooled;
+ # truncate?
+ if n_pooled >= need
+ contribs = contribs(1:need,:);
+ done = true;
+ endif
+ # write to output file
+ FN = fopen (outfile, "a");
+ if (FN < 0) error ("montecarlo: couldn't open output file %s", outfile); endif
+ t = etime(clock(), t0);
+ for j = 1:rows(contribs)
+ fprintf(FN, "%f ", i, t, contribs(j,:));
+ fprintf(FN, "\n");
+ endfor
+ fclose(FN);
+ if verbose printf("\nContribution received from node%d. Received so far: %d\n", i, received); endif
+ if done
+ # tell compute nodes to stop loop
+ for j = 1:5
+ for i = 1:(nodes-1)
+ if (j==1) MPI_Send(" ",i,i,CW); endif # send out message to stop
+ ready = MPI_Iprobe(i, mytag, CW); # get last messages
+ if ready contribs = MPI_Recv(i, mytag, CW); endif
+ endfor
+ endfor
+ break;
+ endif
+ endif
+ endfor
+ else
+ for i = 1:n_pooled
+ contrib = feval(f, f_args);
+ contribs(i,:) = contrib;
+ endfor
+ need = reps - received;
+ received = received + n_pooled;
+ # truncate?
+ if n_pooled >= need
+ contribs = contribs(1:need,:);
+ endif
+ # write to output file
+ FN = fopen (outfile, "a");
+ if (FN < 0) error ("montecarlo: couldn't open output file %s", outfile); endif
+ t = etime(clock(), t0);
+ for j = 1:rows(contribs)
+ fprintf(FN, "%f ", 0, t, contribs(j,:));
+ fprintf(FN, "\n");
+ endfor
+ fclose(FN);
+ if verbose printf("\nContribution received from node 0. Received so far: %d\n", received); endif
+ endif
+ endwhile
+ endif
+endfunction
git://git.creatis.insa-lyon.fr / CreaPhase.git / blobdiff