2012年10月5日 星期五

OpenMP and MPICH2 hybrid example

Base on OPENMP and MPICH2 Quick Install and Examples

libgomp: 4.4.6,GCC OpenMP v3.0 shared support library
mpich2: 1.4.1p1, A high-performance implementation of MPI
libgomp and mpich2-devel should be installed both in cent145 and cent146

[C/C++ Examples]


[Example 1] 

[root@cent146 hybrid]# cat hybrid_hi.c
#include
#include "mpi.h"
#include

int main(int argc, char *argv[]) {
  int numprocs, rank, namelen;
  char processor_name[MPI_MAX_PROCESSOR_NAME];
  int iam = 0, np = 1;

  MPI_Init(&argc, &argv);
  MPI_Comm_size(MPI_COMM_WORLD, &numprocs);
  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
  MPI_Get_processor_name(processor_name, &namelen);

  #pragma omp parallel default(shared) private(iam, np)
  {
    np = omp_get_num_threads();
    iam = omp_get_thread_num();
    printf("Hello from thread %d out of %d from process %d out of %d on %s\n",
           iam, np, rank, numprocs, processor_name);
  }

  MPI_Finalize();
}

[root@cent146 hybrid]# cat make_hybrid_hi.sh
#!/bin/bash
mpicc  -o hybrid_hi hybrid_hi.c -L/usr/lib64/mpich2/lib/ -lmpl -lopa -fopenmp

[root@cent146 hybrid]# mpiexec -f hydra.hosts -n 2 ./hybrid_hi
Hello from thread 0 out of 1 from process 1 out of 2 on cent146
Hello from thread 0 out of 1 from process 0 out of 2 on cent145


[Example 2]

[root@cent146 hybrid]# cat hybrid_pi.c
#include
#include
#include
#define NBIN 100000
#define MAX_THREADS 8
void main(int argc,char **argv) {
  int nbin,myid,nproc,nthreads,tid;
  double step,sum[MAX_THREADS]={0.0},pi=0.0,pig;
  MPI_Init(&argc,&argv);
  MPI_Comm_rank(MPI_COMM_WORLD,&myid);
  MPI_Comm_size(MPI_COMM_WORLD,&nproc);
  nbin = NBIN/nproc; step = 1.0/(nbin*nproc);
#pragma omp parallel private(tid)
  {
    int i;
    double x;
    nthreads = omp_get_num_threads();
    tid = omp_get_thread_num();
    for (i=nbin*myid+tid; i
        x = (i+0.5)*step; sum[tid] += 4.0/(1.0+x*x);}
    printf("rank tid sum = %d %d %e\n",myid,tid,sum[tid]);
  }
  for(tid=0; tid
  MPI_Allreduce(&pi,&pig,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
  if (myid==0) printf("PI = %f\n",pig);
  MPI_Finalize();
}

[root@cent146 hybrid]# cat make_hybrid_pi.sh
#!/bin/bash
        mpicc  -o hybrid_pi hybrid_pi.c -L/usr/lib64/mpich2/lib/ -lmpl -lopa -fopenmp

[root@cent146 hybrid]# mpiexec -f hydra.hosts -n 2 ./hybrid_pi
rank tid sum = 1 0 1.287002e+05
rank tid sum = 0 0 1.854590e+05
PI = 3.141593

[Example 3]

[root@cent146 hybrid]# cat hybrid.c
#include
#include
#include
#include
#ifdef USE_MPI
  #include
#endif /* USE_MPI */
#ifdef _OPENMP
  #include
#endif /* _OPENMP */

int read_slab_info() {
  /* This should read info from a file or something,
     but we fake it */
  return 80;
}

double process_slab(int snum)
{
  int i, j;
  double x;
  for (i = 0; i < 10000; i++)
    for (j = 0; j < 10000; j++)
      x += sqrt((i-j)*(i-j) / (sqrt((i*i) + (j*j)) + 1));
  return x;
}

void exit_on_error(char *message)
{
  fprintf(stderr, "%s\n", message);
#ifdef USE_MPI
  MPI_Finalize();
#endif
  exit(1);
}


int main(int argc, char **argv)
{
  int i, j, p, me, nprocs, num_threads, num_slabs, spp;
  int *my_slabs, *count;
  double x, sum;
#ifdef _OPENMP
  int np;
#endif /* _OPENMP */
#ifdef USE_MPI
  int namelen;
  char processor_name[MPI_MAX_PROCESSOR_NAME];
#endif /* USE_MPI */

#ifdef USE_MPI
  MPI_Init(&argc, &argv);
  MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
  MPI_Comm_rank(MPI_COMM_WORLD, &me);
  MPI_Get_processor_name(processor_name, &namelen);
#else /* USE_MPI */
  nprocs = 1;
  me = 0;
#endif /* USE_MPI */

#ifdef _OPENMP
  np = omp_get_num_procs();
  omp_set_num_threads(np);
  num_threads = omp_get_max_threads();
#else /* _OPENMP */
  num_threads = 1;
#endif /* _OPENMP */

  printf("Process %d of %d", me, nprocs);
#ifdef USE_MPI
  printf(" running on %s", processor_name);
#endif /* USE_MPI */
#ifdef _OPENMP
  printf(" using OpenMP with %d threads",
    num_threads);
#endif /* _OPENMP */
  printf("\n");

  /* Master process reads slab data */
  if (!me) num_slabs = read_slab_info();
#ifdef USE_MPI
  if (MPI_Bcast(&num_slabs, 1, MPI_INT, 0,
      MPI_COMM_WORLD) != MPI_SUCCESS)
    exit_on_error("Error in MPI_Bcast()");
#endif /* USE_MPI */

  if (num_slabs < nprocs)
    exit_on_error("Number of slabs may not exceed \
      number of processes");
  /* maximum number of slabs per process */
  spp = (int)ceil((double)num_slabs /
  (double)nprocs);

  if (!me) printf("No more than %d slabs will \
    assigned to each process\n", spp);

  /* allocate list and count of slabs for each
    process */
  if (!(my_slabs = (int *)malloc(nprocs*spp*
  sizeof(int)))) {
  perror("my_slabs");
  exit(2);
  }

  if (!(count = (int *)malloc(nprocs*sizeof(int)))) {
    perror("count");
    exit(2);
  }

  /* initialize slab counts */
  for (p = 0; p < nprocs; p++) count[p] = 0;
  /* round robin assignment of slabs to processes
    for better potential
   * load balancing
   */
  for (i = j = p = 0; i < num_slabs; i++) {
    my_slabs[p*spp+j] = i;
    count[p]++;
    if (p == nprocs -1)
      p = 0, j++;
    else
      p++;
  }


  /* each process works on its own list of slabs,
     but OpenMP threads
   * divide up the slabs on each process because
     of OpenMP directive
   */

#pragma omp parallel for reduction(+: x)
  for (i = 0; i < count[me]; i++) {
    printf("%d: slab %d being processed", me,
      my_slabs[me*spp+i]);
#ifdef _OPENMP
    printf(" by thread %d", omp_get_thread_num());
#endif /* _OPENMP */
    printf("\n");
    x += process_slab(my_slabs[me*spp+i]);
  }


#ifdef USE_MPI
  if (MPI_Reduce(&x, &sum, 1, MPI_DOUBLE, MPI_SUM, 0,
      MPI_COMM_WORLD) != MPI_SUCCESS)
    exit_on_error("Error in MPI_Reduce()");
#else /* USE_MPI */
  sum = x;
#endif /* USE_MPI */

  if (!me) printf("Sum is %lg\n", sum);

#ifdef USE_MPI
  printf("%d: Calling MPI_Finalize()\n", me);
  MPI_Finalize();
#endif /* USE_MPI */

  exit(0);
}

[root@cent146 hybrid]# cat make_hybrid.sh
#!/bin/bash
if [ "$1" == "" ] || [ "$1" == "hybrid" ]; then
        mpicc  -o hybrid hybrid.c -DUSE_MPI -L/usr/lib64/mpich2/lib/ -lmpl -lopa -lm -fopenmp
elif  [ "$1" == "mpi" ]; then
        mpicc  -o hybrid_mpi hybrid.c -DUSE_MPI -L/usr/lib64/mpich2/lib/ -lmpl -lopa -lm
elif  [ "$1" == "omp" ]; then
        mpicc  -o hybrid_omp hybrid.c -L/usr/lib64/mpich2/lib/ -lm -fopenmp
else
        echo "Illegal make parameter, try:"
        echo "$0 < hybrid / mpi / omp >"
fi

[root@cent146 hybrid]# mpiexec -f hydra.hosts -n 2 ./hybrid
Process 1 of 2 running on cent146 using OpenMP with 1 threads
Process 0 of 2 running on cent145 using OpenMP with 1 threads
No more than 40 slabs will     assigned to each process
0: slab 0 being processed by thread 0
1: slab 1 being processed by thread 0
1: slab 3 being processed by thread 0
0: slab 2 being processed by thread 0
1: slab 5 being processed by thread 0
0: slab 4 being processed by thread 0
    :
    :
    :
[root@cent146 hybrid]# mpiexec -f hydra.hosts -n 2 ./hybrid_mpi
Process 1 of 2 running on cent146
Process 0 of 2 running on cent145
No more than 40 slabs will     assigned to each process
0: slab 0 being processed
1: slab 1 being processed
1: slab 3 being processed
0: slab 2 being processed
1: slab 5 being processed
0: slab 4 being processed
    :
    :
    :
[root@cent146 hybrid]# ./hybrid_ompProcess 0 of 1 using OpenMP with 1 threads
No more than 80 slabs will     assigned to each process
0: slab 0 being processed by thread 0
0: slab 1 being processed by thread 0
0: slab 2 being processed by thread 0
0: slab 3 being processed by thread 0

[Fortran Examples]


[Example 1]

[root@cent146 hybrid]# cat hi_f.f
        program f1 !Hello World MPI/F90 style
        implicit none
        include 'mpif.h'
        integer::myrank, numprocs,strlen, ierr, comm = MPI_COMM_WORLD
        character(80)::str

        call mpi_init(ierr)
        call mpi_comm_size(comm, numprocs, ierr)
        call mpi_comm_rank(comm, myrank, ierr)

        if (myrank .EQ. 0) print *, 'Num procs ', numprocs

        call mpi_get_processor_name(str, strlen,ierr)

        print *, 'Hello world : I am processor ', myrank, ':', str

        call mpi_finalize(ierr)

        end

[root@cent146 hybrid]# mpif77 -o hi_f hi_f.f

[root@cent146 hybrid]# mpiexec -f hydra.hosts -n 2 ./hi_f
 Hello world : I am processor            1 :cent146                                     
 Num procs            2
 Hello world : I am processor            0 :cent145                                     

[Example 2]

[root@cent146 hybrid]# cat hybrid_hi_f.f
        program hybrid_hi_f ! hello from OMP_threads
        implicit none
        include 'mpif.h'
        integer nthreads, mythrd, omp_get_num_threads
        integer::myrank, numprocs, ierr,strlen, comm=MPI_COMM_WORLD
        !integer myrank numprocs, ierr, comm
        !parameter (comm = MPI_COMM_WORLD)
        character(40)::str
        integer omp_get_thread_num

        call mpi_init(ierr)
        call mpi_comm_size(comm, numprocs, ierr)
        call mpi_comm_rank(comm, myrank, ierr)


!$OMP PARALLEL PRIVATE(mythrd) SHARED(nthreads)

        nthreads = omp_get_num_threads()
        mythrd = omp_get_thread_num()

        call mpi_get_processor_name(str, strlen, ierr)

        if (mythrd.eq.0) print *, "Num threads = ", nthreads
        print *, 'th ' , mythrd, 'of', nthreads, 'prc', myrank , ':' , str
        !print *, ' proc ', myrank , ':' , str


!$OMP END PARALLEL

        call mpi_finalize(ierr)
        end

[root@cent146 hybrid]# mpif77 -o hybrid_hi_f hybrid_hi_f.f -fopenmp

[root@cent146 hybrid]# mpiexec -f hydra.hosts -n 2 ./hybrid_hi_f
 Num threads =            1
 Num threads =            1
 th            0 of           1 prc           1 :cent146                                
 th            0 of           1 prc           0 :cent145                                

[reference]

http://www.math.ntu.edu.tw/~wwang/cola_lab/knowledge/download/Parallel_Computing/7-Hybrid%20MPI+OpenMP.ppt
http://www.linux-mag.com/id/1631/
張貼者： 於 沒有留言:
標籤： , , , ,

2012年10月3日 星期三

OPENMP and MPICH2 Quick Install and Examples

Environment: (@VM)


        host : CentOS-6.0
        kernel: 2.6.32-71.el6.x86_64 #1 SMP
        MemTotal: 7891424 kB(/proc/meminfo)
        CPU(Quad core) : Intel(R) Core(TM) i3-2130 CPU @ 3.40GHz
        HDD： 1TB＊2 （LVM ＝＝＞ 2T）

        guest: CentOS-6.0
        kernel: 2.6.32-71.el6.x86_64 #1 SMP
        MemTotal: 1G
        CPU *1
        HDD： 200GB 
172.16.173.143     172.16.173.144
+------------+     +------------+
|mpich2      |     |mpich2      |
|            |     |            |
|------------|     |------------|     
|cent143     |     |cent144     |
+-----+------+     +------+-----+            +----NAT(@VM)
      |                  |                  |
------+------------------+------------------+
      |            
+-----+------+     
|mpich2      |     
|openmp      |     
|------------|     
|cent142     |     
+------------+     
172.16.173.142     

Purpose:

1. install a openmp environment in a server and test it
2. install a mpi(with mpich2) environment in all servers and test it

OPENMP:

[Setup]


[root@cent142 ~]# yum install libgomp.x86_64 gcc-gfortran-4.4.6-4.el6.x86_64 -y


[C/C++ Code Example]

[root@cent142 ~]# mkdir -p ~/src/openmp; cd ~/src/openmp
[root@cent142 openmp]# vi ttt.c
#include
#include
int main()
{
        int x;
        x=2;
        #pragma omp parallel num_threads(2) shared(x)
        {
                if(omp_get_thread_num() == 1)
                        x = 5;
                else
                {
                        printf("1: Thread %d: x = %d\n", omp_get_thread_num(),x);
                }

                #pragma omp barrier
                if( omp_get_thread_num() == 0)
                {
                        printf("2: Thread# %d, x = %d\n", omp_get_thread_num(),x);
                }
                else
                {
                        printf("3: Thread# %d, x = %d\n", omp_get_thread_num(),x);
                }
        }
        return 0;
}
[root@cent142 openmp]# gcc -fopenmp -o ttt ttt.c

[root@cent142 openmp]# ./ttt

1: Thread 0: x = 2
2: Thread# 0, x = 5
3: Thread# 1, x = 5



[Fortran Code Example]




[root@cent142 openmp]# vi ttt.f
        PROGRAM A2
        INCLUDE "omp_lib.h"
        INTEGER X
        X = 2
!$OMP PARALLEL NUM_THREADS(2) SHARED(X)
        IF( OMP_GET_THREAD_NUM() .EQ. 0) THEN
                X=5
        ELSE
                ! PRINT 1: The following read of x has a race
                  PRINT *,"1: THREAD# ", OMP_GET_THREAD_NUM(), "X = ", X
        ENDIF

!$OMP BARRIER
        IF (OMP_GET_THREAD_NUM() .EQ. 0) THEN
                ! PRINT 2
                  PRINT *,"2: THREAD# ", OMP_GET_THREAD_NUM(), "X =", X
        ELSE
                ! PRINT 3
                  PRINT *,"3: THREAD# ", OMP_GET_THREAD_NUM(), "X =", X
        ENDIF
!$OMP END PARALLEL

        END
[root@cent142 openmp]# gfortran -fopenmp -o ttt.for ttt.f

[root@cent142 openmp]# ./ttt.for

 1: THREAD#            1 X =            5
 2: THREAD#            0 X =           5
 3: THREAD#            1 X =           5


MPICH2:

[Setup]

(1) Install Packages
[root@cent142 ~]# yum install mpich2.x86_64 mpich2-devel.x86_64
 -y

[root@cent143 ~]# yum install mpich2.x86_64 mpich2-devel.x86_64
 -y


[root@cent144 ~]# yum install mpich2.x86_64 mpich2-devel.x86_64
 -y
[root@cent142 mpich2]# mpich2version
MPICH2 Version:         1.2.1
MPICH2 Release date:    Unknown, built on Fri Nov 12 05:07:30 GMT 2010
MPICH2 Device:          ch3:nemesis
MPICH2 configure:       --build=x86_64-unknown-linux-gnu....
MPICH2 CC:      gcc ....
MPICH2 CXX:     c++ ....
MPICH2 F77:     gfortran ....
MPICH2 F90:     gfortran ....




(2)/etc/hosts for all servers(scp /etc/hosts cent14X:/etc/.)
172.16.173.142  cent142
172.16.173.143  cent143
172.16.173.144  cent144


(3)a hosts file for mpi(~/src/mpich2/mpd.hosts)
cent142
cent143
cent144

(4)mpd.conf for all servers( password when inserting a task)
## all servers must have the same password
cat > /etc/mpd.conf << "EOF"
secretword=********
EOF

chmod 600 /etc/mpd.conf

scp /etc/mpd.conf cent142:/etc/.
scp /etc/mpd.conf cent143:/etc/.
scp /etc/mpd.conf cent144:/etc/.

(5)password when client try to insert a task to server
##it must the same with server
cat > ~/.mpd.conf << "EOF"
secretword=********
EOF

(6)password-less ssh

(7)start up mpich2
cd ~/src/mpich2
mpdboot -n 3 -f ./mpd.hosts
mpdboot_cent142 (handle_mpd_output 420): from mpd on cent142, invalid port info:
no_port
(1) check the ~/.mpd.conf exist , and only named as ~/.mpd.conf
(2) permission 600 for ~/.mpd.conf and /etc/mpd.conf
(3) x86 - x64 hybrid architechure is not allow for mpich2
    openmpi or mpich2 instead

mpdtrace
cent142
cent144
cent143

[root@cent142 mpich2]# mpiexec  -n 6  /bin/hostname
cent143
cent142
cent142
cent143
cent144
cent144

[root@cent142 mpich2]# mpiexec  -n 3 date : -n 4 hostname
cent144
銝?10?? 3 18:53:19 CST 2012
cent143
銝?10?? 3 18:53:19 CST 2012
cent142
cent142
銝?10?? 3 18:53:19 CST 2012



[root@cent142 mpich2]# mpiexec -n 3 /bin/date : -n 2 -host cent143 /bin/hostname
瞻T 10瞻禱  3 15:02:00 CST 2012
cent143
cent143
瞻T 10瞻禱  3 15:02:00 CST 2012
瞻T 10瞻禱  3 15:02:00 CST 2012

[root@cent142 mpich2]# mpiexec -machinefile mpd2.hosts -n 3 /bin/date
銝?10?? 3 18:54:47 CST 2012
銝?10?? 3 18:54:47 CST 2012
銝?10?? 3 18:54:47 CST 2012

(8)close mpd service for all servers
mpdallexit

[C/C++ Code Example]

[root@cent142 ~]# mkdir -p ~/src/mpich2; cd ~/src/mpich2

[root@cent142 mpich2]# cat pi_cc.cc

#include
#include "mpi.h"
#include

int main(int argc, char *argv[])
{
    int n, rank, size, i;
    double PI25DT = 3.141592653589793238462643;
    double mypi, pi, h, sum, x;

    MPI::Init(argc, argv);
    size = MPI::COMM_WORLD.Get_size();
    rank = MPI::COMM_WORLD.Get_rank();

    while (1) {
        if (rank == 0) {
            std::cout << "Enter the number of intervals: (0 quits)"
                 << std::endl;
            std::cin >> n;
        }

        MPI::COMM_WORLD.Bcast(&n, 1, MPI::INT, 0);
        if (n==0)
            break;
        else {
            h = 1.0 / (double) n;
            sum = 0.0;
            for (i = rank + 1; i <= n; i += size) {
                x = h * ((double)i - 0.5);
                sum += (4.0 / (1.0 + x*x));
            }
            mypi = h * sum;

            MPI::COMM_WORLD.Reduce(&mypi, &pi, 1, MPI::DOUBLE,
                                   MPI::SUM, 0);
            if (rank == 0)
                std::cout << "pi is approximately " << pi
                     << ", Error is " << fabs(pi - PI25DT)
                     << std::endl;
        }
    }
    MPI::Finalize();
    return 0;
}

[root@cent142 mpich2]# mpic++  -o pi_cc pi_cc.cc

[root@cent142 mpich2]# mpdboot -n 3 -f ./mpd.hosts
[root@cent142 mpich2]# mpiexec -n 12 ./pi_cc

[root@cent142 mpich2]# ps -ef | grep pi_cc
root      4039  2350  0 16:27 pts/0    00:00:00 python2.6 /usr/bin/mpiexec -n 12 ./pi_cc
root      4045  4043 49 16:27 ?        00:01:01 ./pi_cc
root      4046  4044 49 16:27 ?        00:01:00 ./pi_cc
root      4047  4041  0 16:27 ?        00:00:00 ./pi_cc
root      4052  1357  0 16:29 pts/1    00:00:00 grep pi_cc

[root@cent143 ~]# ps -ef | grep pi_cc
root      4171  4166 20 16:27 ?        00:00:29 ./pi_cc
root      4172  4168 19 16:27 ?        00:00:27 ./pi_cc
root      4173  4169 18 16:27 ?        00:00:26 ./pi_cc
root      4174  4165 20 16:27 ?        00:00:29 ./pi_cc
root      4175  4170 20 16:27 ?        00:00:28 ./pi_cc
root      4180  2277  0 16:30 pts/0    00:00:00 grep pi_cc

[root@cent144 ~]# ps -ef | grep pi_cc
root      3696  3692 25 16:27 ?        00:00:16 ./pi_cc
root      3697  3691 22 16:27 ?        00:00:14 ./pi_cc
root      3698  3694 25 16:27 ?        00:00:16 ./pi_cc
root      3699  3695 23 16:27 ?        00:00:15 ./pi_cc
root      3703  2428  0 16:28 pts/0    00:00:00 grep pi_cc


Enter the number of intervals: (0 quits)
50
pi is approximately 3.14163, Error is 3.33333e-05
Enter the number of intervals: (0 quits)
0
[root@cent142 mpich2]# mpdallexit

[Fortran Code Example]

[root@cent142 mpich2]# cat pi_ff.f
      program main
      include "mpif.h"
      double precision  PI25DT
      parameter        (PI25DT = 3.141592653589793238462643d0)
      double precision  mypi, pi, h, sum, x, f, a
      integer n, myid, numprocs, i, ierr
c                                 function to integrate
      f(a) = 4.d0 / (1.d0 + a*a)

      call MPI_INIT(ierr)
      call MPI_COMM_RANK(MPI_COMM_WORLD, myid, ierr)
      call MPI_COMM_SIZE(MPI_COMM_WORLD, numprocs, ierr)


 10   if ( myid .eq. 0 ) then
         print *, 'Enter the number of intervals: (0 quits) '
         read(*,*) n
      endif
c                                 broadcast n
      call MPI_BCAST(n,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
c                                 check for quit signal
      if ( n .le. 0 ) goto 30
c                                 calculate the interval size
      h = 1.0d0/n
      sum  = 0.0d0
      do 20 i = myid+1, n, numprocs
         x = h * (dble(i) - 0.5d0)
         sum = sum + f(x)
 20   continue
      mypi = h * sum
c                                 collect all the partial sums
      call MPI_REDUCE(mypi,pi,1,MPI_DOUBLE_PRECISION,MPI_SUM,0,
     &                  MPI_COMM_WORLD,ierr)
c                                 node 0 prints the answer.
      if (myid .eq. 0) then
         print *, 'pi is ', pi, ' Error is', abs(pi - PI25DT)
      endif
      goto 10
 30   call MPI_FINALIZE(ierr)
      stop
      end

[root@cent142 mpich2]# cat make_pi_ff.sh

#!/bin/bash

## mpif77 -f77=gfortran -o pi_ff pi_ff.f
## mpif77 -f77=f77 -o pi_ff pi_ff.f      <-- so different with gfortran>
mpif77 -o pi_ff pi_ff.f

[root@cent142 mpich2]# ./make_pi_ff.sh

[root@cent142 mpich2]# mpdboot -n 3 -f ./mpd.hosts
[root@cent142 mpich2]# mpiexec -n 12 ./pi_ff

10

0
 Enter the number of intervals: (0 quits)
 pi is    3.1424259850010983       Error is  8.33331411305149317E-004
 Enter the number of intervals: (0 quits)


[A]comparasion between -n when computating of pi with Monte Carlo



[root@cent142 mpich2]# time mpiexec -n 3 ./monte_cc 0.00001
pi =  3.14159455782312946326
points: 735000
in: 577268, out: 157732, to exit


real    0m8.805s
user    0m0.052s
sys     0m0.015s
[root@cent142 mpich2]# time mpiexec -n 6 ./monte_cc 0.00001
pi =  3.14159455782312946326
points: 735000
in: 577268, out: 157732, to exit


real    0m5.406s
user    0m0.047s
sys     0m0.015s

[root@cent142 mpich2]# time mpiexec -n 42 ./monte_cc 0.00001
pi =  3.14158287705326033645
points: 1004500
in: 788930, out: 215570, to exit


real    0m6.829s
user    0m0.053s
sys     0m0.012s

[source code]


[root@cent142 mpich2]# cat monte_cc.c
#include
#include
#include "mpi.h"
#include "mpe.h"
#define CHUNKSIZE      1000
/* We'd like a value that gives the maximum value returned by the function
   random, but no such value is *portable*.  RAND_MAX is available on many
   systems but is not always the correct value for random (it isn't for
   Solaris).  The value ((unsigned(1)<<31 but="but" common="common" guaranteed="guaranteed" is="is" not="not" p="p">#define INT_MAX 1000000000

/* message tags */
#define REQUEST  1
#define REPLY    2
int main( int argc, char *argv[] )
{
    int iter;
    int in, out, i, iters, max, ix, iy, ranks[1], done, temp;
    double x, y, Pi, error, epsilon;
    int numprocs, myid, server, totalin, totalout, workerid;
    int rands[CHUNKSIZE], request;
    MPI_Comm world, workers;
    MPI_Group world_group, worker_group;
    MPI_Status status;

    if(argc < 2)
    {
        printf("you must give a number of tolerance of pi when approaching it!!, ex:\n mpiexec -n ./monte_cc 0.01\n");
        return -1;
    }
    MPI_Init(&argc,&argv);
    world  = MPI_COMM_WORLD;
    MPI_Comm_size(world,&numprocs);
    MPI_Comm_rank(world,&myid);
    server = numprocs-1;        /* last proc is server */
    if (myid == 0)
        sscanf( argv[1], "%lf", &epsilon );
    MPI_Bcast( &epsilon, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD );
    MPI_Comm_group( world, &world_group );
    ranks[0] = server;
    MPI_Group_excl( world_group, 1, ranks, &worker_group );
    MPI_Comm_create( world, worker_group, &workers );
    MPI_Group_free(&worker_group);
    if (myid == server) {       /* I am the rand server */
        do {
            MPI_Recv(&request, 1, MPI_INT, MPI_ANY_SOURCE, REQUEST,
                     world, &status);
            if (request) {
                for (i = 0; i < CHUNKSIZE; ) {
                        rands[i] = random();
                        if (rands[i] <= INT_MAX) i++;
                }
                MPI_Send(rands, CHUNKSIZE, MPI_INT,
                         status.MPI_SOURCE, REPLY, world);
            }
        }
        while( request>0 );
    }
    else {                      /* I am a worker process */
        request = 1;
        done = in = out = 0;
        max  = INT_MAX;         /* max int, for normalization */
        MPI_Send( &request, 1, MPI_INT, server, REQUEST, world );
        MPI_Comm_rank( workers, &workerid );
        iter = 0;
        while (!done) {
            iter++;
            request = 1;
            MPI_Recv( rands, CHUNKSIZE, MPI_INT, server, REPLY,
                     world, &status );
            for (i=0; i                x = (((double) rands[i++])/max) * 2 - 1;
                y = (((double) rands[i++])/max) * 2 - 1;
                if (x*x + y*y < 1.0)
                    in++;
                else
                    out++;
            }
            MPI_Allreduce(&in, &totalin, 1, MPI_INT, MPI_SUM,
                          workers);
            MPI_Allreduce(&out, &totalout, 1, MPI_INT, MPI_SUM,
                          workers);
            Pi = (4.0*totalin)/(totalin + totalout);
            error = fabs( Pi-3.141592653589793238462643);
            done = (error < epsilon || (totalin+totalout) > 1000000);
            request = (done) ? 0 : 1;
            if (myid == 0) {
                printf( "\rpi = %23.20f", Pi );
                MPI_Send( &request, 1, MPI_INT, server, REQUEST,
                         world );
            }
            else {
                if (request)
                    MPI_Send(&request, 1, MPI_INT, server, REQUEST,
                             world);
            }
        }
        MPI_Comm_free(&workers);
    }

    if (myid == 0) {
        printf( "\npoints: %d\nin: %d, out: %d, to exit\n",
               totalin+totalout, totalin, totalout );
        getchar();
    }
    MPI_Finalize();
}


[B]mpiexec for fortran

[root@cent142 mpich2]# mpiexec -n 3 ./pitm_ff
20
0
 Enter the number of intervals: (0 quits)
 pi is    3.1418009868930934       Error is  2.08333303300278772E-004
 time is   7.91406631469726563E-003  seconds
 Enter the number of intervals: (0 quits)


[source code]


[root@cent142 mpich2]# cat pitm_ff.f
      program main
      include "mpif.h"
      double precision  PI25DT
      parameter        (PI25DT = 3.141592653589793238462643d0)
      double precision  mypi, pi, h, sum, x, f, a
      double precision starttime, endtime
      integer n, myid, numprocs, i, ierr
c                                 function to integrate
      f(a) = 4.d0 / (1.d0 + a*a)

       call MPI_INIT(ierr)
      call MPI_COMM_RANK(MPI_COMM_WORLD, myid, ierr)
      call MPI_COMM_SIZE(MPI_COMM_WORLD, numprocs, ierr)
 10   if ( myid .eq. 0 ) then
         print *, 'Enter the number of intervals: (0 quits) '
         read(*,*) n
      endif
c                                 broadcast n
      starttime = MPI_WTIME()
      call MPI_BCAST(n,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
c                                 check for quit signal
      if ( n .le. 0 ) goto 30
c                                 calculate the interval size
      h = 1.0d0/n
      sum  = 0.0d0
      do 20 i = myid+1, n, numprocs
         x = h * (dble(i) - 0.5d0)
         sum = sum + f(x)
 20   continue
      mypi = h * sum
c                                 collect all the partial sums
      call MPI_REDUCE(mypi,pi,1,MPI_DOUBLE_PRECISION,MPI_SUM,0,
     &                  MPI_COMM_WORLD,ierr)
c                                 node 0 prints the answer.
      endtime = MPI_WTIME()
      if (myid .eq. 0) then
         print *, 'pi is ', pi, ' Error is', abs(pi - PI25DT)
         print *, 'time is ', endtime-starttime, ' seconds'
      endif
      goto 10
 30   call MPI_FINALIZE(ierr)
      stop
      end


[Reference]



張貼者： 於 4 則留言:
標籤： , , ,
訂閱： 文章 (Atom)

文章分類