--- /dev/null
+## Copyright (C) 2004-2007 Javier Fernández Baldomero, Mancia Anguita López
+## This code has been adjusted for octave3.2.3 and octave 3.3.50+ in
+## 2009 by Riccardo Corradini <riccardocorradini@yahoo.it>
+##
+## 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, see <http://www.gnu.org/licenses/>.
+
+
+# mpirun -np 5 octave -q --eval "Pi(2E7,'s')"
+
+function Pi(N,mod)
+# Pi: Classic PI computation by numeric integration of arctan'(x) in [0..1]
+#
+# Pi [ ( N [ ,mod ] ) ]
+#
+# N [1E7] #subdivisions of the [0, 1] interval
+# mod ['s'] communication modality: (s)end (r)educe
+#
+# printed results struct contains
+# pi estimated pi value
+# err error
+# time from argument xmit to pi computed
+#
+
+
+##########
+# ArgChk #
+##########
+if nargin<1, N=1E7; end
+if nargin<2, mod='s'; end
+if nargin>2, usage("Pi(N,mod)"); end # let all ranks complain
+flag=0; # code much simpler
+flag=flag || ~isscalar(N) || ~isnumeric(N);
+flag=flag | fix(N)~=N | N<1;
+ mod=lower(mod); mods='sr';
+flag=flag | isempty(findstr(mod, mods)); # let them all error out
+if flag, usage("Pi( <int> N>0, <char> mod=='s|r' )"); end
+
+##################
+# Results struct #
+##################
+results.pi =0;
+results.err =0;
+results.time =0;
+
+
+############
+# PARALLEL # initialization, include MPI_Init time in measurement
+############
+ T=clock; #
+############
+ MPI_ANY_SOURCE = -1;
+ MPI_Init();
+ MPI_COMM_WORLD = MPI_Comm_Load("NEWORLD");
+ rnk = MPI_Comm_rank (MPI_COMM_WORLD); # let it abort if it fails
+ siz = MPI_Comm_size (MPI_COMM_WORLD);
+
+ SLV = logical(rnk); # handy shortcuts, master is rank 0
+ MST = ~ SLV; # slaves are all other
+
+############
+# PARALLEL # computation (depends on rank/size)
+############ # vectorized code, equivalent to
+ width=1/N; lsum=0; # for i=rnk:siz:N-1
+ i=rnk:siz:N-1; # x=(i+0.5)*width;
+ x=(i+0.5)*width; # lsum=lsum+4/(1+x^2);
+ lsum=sum(4./(1+x.^2)); # end
+
+############
+# PARALLEL # reduction and finish
+############
+switch mod
+ case 's', TAG=7; # Any tag would do
+ if SLV # All slaves send result back
+ MPI_Send(lsum, 0,TAG,MPI_COMM_WORLD);
+ else # Here at master
+ Sum =lsum; # save local result
+ for slv=1:siz-1 # collect in any order
+ lsum = MPI_Recv(MPI_ANY_SOURCE,TAG,MPI_COMM_WORLD);
+ Sum+=lsum; # and accumulate
+ end # order: slv or MPI_ANY_SOURCE
+ end
+ case 'r',
+ disp("not yet implemented");
+# Sum=0;
+# reduction master = rank 0 @ WORLD
+# MPI_Reduce(lsum,Sum, MPI_SUM, 0,MPI_COMM_WORLD);
+end
+
+MPI_Finalize();
+
+if MST
+ Sum = Sum/N ; # better at end: don't loose resolution
+################################# # stopwatch measurement
+results.time = etime(clock,T); # # but only at master after PI computed
+################################# # all them started T=clock;
+results.err = Sum-pi;
+results.pi = Sum # ;
+
+end
\ No newline at end of file