octave_packages/openmpi_ext-1.0.2/Pi.m

   1 ## Copyright (C) 2004-2007 Javier Fernández Baldomero, Mancia Anguita López
   2 ## This code has been adjusted for octave3.2.3 and octave 3.3.50+ in
   3 ## 2009 by  Riccardo Corradini <riccardocorradini@yahoo.it>
   4 ##
   5 ## This program is free software; you can redistribute it and/or modify
   6 ## it under the terms of the GNU General Public License as published by
   7 ## the Free Software Foundation; either version 2 of the License, or
   8 ## (at your option) any later version.
   9 ##
  10 ## This program is distributed in the hope that it will be useful,
  11 ## but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13 ## GNU General Public License for more details.
  14 ##
  15 ## You should have received a copy of the GNU General Public License
  16 ## along with this program; If not, see <http://www.gnu.org/licenses/>.
  17
  18
  19 # mpirun -np 5 octave -q --eval "Pi(2E7,'s')"
  20
  21 function Pi(N,mod)
  22 # Pi:   Classic PI computation by numeric integration of arctan'(x) in [0..1]
  23 #
  24 #       Pi [ ( N [ ,mod ] ) ]
  25 #
  26 #  N    [1E7]   #subdivisions of the [0, 1] interval
  27 #  mod  ['s']   communication modality:  (s)end (r)educe
  28 #
  29 #  printed results struct contains
  30 #       pi      estimated pi value
  31 #       err     error
  32 #       time    from argument xmit to pi computed
  33 #
  34
  35
  36 ##########
  37 # ArgChk #
  38 ##########
  39 if nargin<1,    N=1E7;  end
  40 if nargin<2,  mod='s';  end
  41 if nargin>2,    usage("Pi(N,mod)"); end         # let all ranks complain
  42 flag=0;                                         # code much simpler
  43 flag=flag || ~isscalar(N) || ~isnumeric(N);
  44 flag=flag  |   fix(N)~=N   |           N<1;
  45                    mod=lower(mod); mods='sr';
  46 flag=flag  | isempty(findstr(mod,  mods));      # let them all error out
  47 if flag,        usage("Pi( <int> N>0, <char> mod=='s|r' )"); end
  48
  49 ##################
  50 # Results struct #
  51 ##################
  52 results.pi   =0;
  53 results.err  =0;
  54 results.time =0;
  55
  56
  57 ############
  58 # PARALLEL # initialization, include MPI_Init time in measurement
  59 ############
  60   T=clock; #
  61 ############
  62    MPI_ANY_SOURCE = -1;
  63    MPI_Init();
  64    MPI_COMM_WORLD = MPI_Comm_Load("NEWORLD");
  65    rnk   =      MPI_Comm_rank (MPI_COMM_WORLD); # let it abort if it fails
  66    siz   =      MPI_Comm_size (MPI_COMM_WORLD);
  67
  68     SLV = logical(rnk);                 # handy shortcuts, master is rank 0
  69     MST = ~ SLV;                        # slaves are all other
  70
  71 ############
  72 # PARALLEL # computation (depends on rank/size)
  73 ############                    # vectorized code, equivalent to
  74   width=1/N; lsum=0;            # for i=rnk:siz:N-1
  75   i=rnk:siz:N-1;                #   x=(i+0.5)*width;
  76   x=(i+0.5)*width;              #   lsum=lsum+4/(1+x^2);
  77   lsum=sum(4./(1+x.^2));        # end
  78
  79 ############
  80 # PARALLEL # reduction and finish
  81 ############
  82 switch mod
  83   case 's',                     TAG=7;  # Any tag would do
  84     if SLV                              # All slaves send result back
  85         MPI_Send(lsum,             0,TAG,MPI_COMM_WORLD);
  86     else                                # Here at master
  87             Sum =lsum;                  # save local result
  88       for slv=1:siz-1                   # collect in any order
  89             lsum = MPI_Recv(MPI_ANY_SOURCE,TAG,MPI_COMM_WORLD);
  90             Sum+=lsum;                  # and accumulate
  91       end                               # order: slv or MPI_ANY_SOURCE
  92     end
  93   case 'r',
  94         disp("not yet implemented");
  95 #       Sum=0;
  96 # reduction master = rank 0 @ WORLD
  97 #       MPI_Reduce(lsum,Sum, MPI_SUM,  0,MPI_COMM_WORLD);
  98 end
  99
 100 MPI_Finalize();
 101
 102 if MST
 103     Sum      = Sum/N ;                  # better at end: don't loose resolution
 104 #################################       # stopwatch measurement
 105 results.time = etime(clock,T);  #       # but only at master after PI computed
 106 #################################       # all them started T=clock;
 107 results.err  = Sum-pi;
 108 results.pi   = Sum # ;
 109
 110 end