1 # Copyright (C) 2006, 2009 Michael Creel <michael.creel@uab.es>
2 # under the terms of the GNU General Public License.
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.
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.
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
18 # montecarlo.m: generates a specified number of replications of a function's
19 # output and writes them to a user-specified output file.
21 # USAGE: montecarlo(f,f_args,reps,outfile,n_pooled,n_returns,usempi, verbose)
23 # IMPORTANT: f should return a row vector of output from feval(f,f_args)
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
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.
37 function n_received = montecarlo(f,f_args,reps,outfile,n_pooled,verbose)
39 t0 = clock(); # initialize timing
41 # defaults for optional arguments
42 if (nargin < 6) verbose = false; endif
43 if (nargin < 5) n_pooled = 1; endif;
45 if MPI_Initialized # check if doing this parallel or serial
47 CW = MPI_Comm_Load("NEWORLD");
48 is_node = MPI_Comm_rank(CW);
49 nodes = MPI_Comm_size(CW);
56 if is_node # compute nodes
60 contrib = feval(f, f_args);
61 contribs(i,:) = contrib;
63 MPI_Send(contribs, 0, mytag, CW);
64 pause(0.05); # give time for the fronted to send a stop message, if done
66 if (MPI_Iprobe(0, is_node, CW)) # check for ping from rank 0
67 junk = MPI_Recv(0, is_node, CW);
76 # retrieve results from compute nodes
78 # compute nodes have results yet?
80 ready = MPI_Iprobe(i, mytag, CW); # check if message pending
82 # get it if it's there
83 contribs = MPI_Recv(i, mytag, CW);
84 need = reps - received;
85 received = received + n_pooled;
88 contribs = contribs(1:need,:);
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,:));
100 if verbose printf("\nContribution received from node%d. Received so far: %d\n", i, received); endif
102 # tell compute nodes to stop loop
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
116 contrib = feval(f, f_args);
117 contribs(i,:) = contrib;
119 need = reps - received;
120 received = received + n_pooled;
123 contribs = contribs(1:need,:);
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,:));
134 if verbose printf("\nContribution received from node 0. Received so far: %d\n", received); endif