+= does not exist in scilab ^^
[scilab.git] / scilab / modules / mpi / macros / MPIPi.sci
1 function MPIPi(N,mod)
2 // Pi_seq: Compute PI (sequential) by num.integr. of arctan(x) in [0..1]
3 //
4 //      Pi_seq [ (N) ]
5 //
6 //  N   [1E7]   #subdivisions of the [0, 1] interval
7 // mod  ['s']   communication modality:  (s)end (r)educe
8 //  printed results struct contains
9 //      pi      estimated pi value
10 //      err     error
11 //      time    from argument xmit to pi computed
12 //
13 // Examples:
14 //
15 //  octave:1> Pi_seq
16 //  results =
17 //  {
18 //    pi =  3.1416
19 //    err = -6.4837e-14
20 //    time =  2.7761
21 //  }
22 //
23
24 ////////////////////
25 // ArgChk //
26 ////////////////////
27 if argn(2)<1,   N=1E7;  end
28 if argn(2)<2,  mod='s'; end
29 if argn(2)>2,   error("usage MPIPi(N,mod)"); end
30 flag=0;
31 [%H,%ierr] = evstr(string(N));
32 flag=flag | size(N,"*")<>1 | %ierr<>0;
33 flag=flag  |   fix(N)~=N   |           N<1;
34 if flag,        error("usage MPIPi( <int> N>0, <char> mod==''s|r'' )"), end
35
36 ////////////////////////////////////
37 // Results struct //
38 ////////////////////////////////////
39 results.pi  =0;
40 results.err =0;
41 results.time=0;
42
43 ////////////
44 // PARALLEL / initialization, include MPI_Init time in measurement
45 ////////////
46 //  T=clock; /
47 ////////////
48                 MPI_Init();                     // should have lambooted first
49    rnk =        MPI_Comm_rank();        // let it abort if it fails
50    siz =        MPI_Comm_size();
51
52     SLV = rnk;                  // handy shortcuts, master is rank 0
53     MST = ~ SLV;                        // slaves are all other
54
55 if MST 
56 disp("Master here")
57 else
58 disp("Slave here")
59 end
60
61 ////////////
62 // PARALLEL / computation (depends on rank/size)
63 ////////////                    // vectorized code, equivalent to
64   width=1/N; lsum=0;            // for i=rnk:siz:N-1
65   i=rnk:siz:N-1;                //   x=(i+0.5)*width;
66   x=(i+0.5)*width;              //   lsum=lsum+4/(1+x^2);
67   lsum=sum(4 ./(1+x.^2));       // end
68                                 // mem-bound, N=1E7 => size(i)=8E7/siz (80MB!!!)
69 ////////////
70 // PARALLEL / reduction and finish
71 ////////////
72 select mod
73   case 's',                     TAG=7;  // Any tag would do
74     if SLV                              // All slaves send result back
75         MPI_Send(lsum,             0,TAG);
76     else                                // Here at master
77             Sum =lsum;                  // save local result
78       for slv=1:siz-1                   // collect in any order
79 //      MPI_Recv(lsum,MPI_ANY_SOURCE,TAG,MPI_COMM_WORLD);
80                 MPI_Recv(lsum,TAG);
81             Sum = Sum + lsum;                   // and accumulate
82       end                               // order: slv or MPI_ANY_SOURCE
83     end
84   case 'r',             Sum=0;          // reduction master = rank 0 @ WORLD
85         MPI_Reduce(lsum,Sum, MPI_SUM,  0,MPI_COMM_WORLD);
86 end
87
88                 MPI_Finalize;
89
90 //////////////////////////////////////////////////////////////////
91 //results.time = etime(clock,T);        //
92 //////////////////////////////////////////////////////////////////
93 results.err  = Sum-%pi;
94 results.pi   = Sum // ;
95
96
97 //////////////////////////////
98 // FINAL CHECK //
99 //////////////////////////////
100   if abs(results.err)>5E-10
101     warning("pretty nice pi value! go fix it")
102   end
103 disp(results)
104 endfunction
105
106