Add a useful package (from Source forge) for octave

[CreaPhase.git] / octave_packages / openmpi_ext-1.0.2 / Pi.m

diff --git a/octave_packages/openmpi_ext-1.0.2/Pi.m b/octave_packages/openmpi_ext-1.0.2/Pi.m
new file mode 100644 (file)
index 0000000..2f764b8
--- /dev/null
+++ b/octave_packages/openmpi_ext-1.0.2/Pi.m
@@ -0,0 +1,110 @@
+## 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