MPI parallel and serial fortran and c examples that calculate PI

cpi_parallel.c

#include "mpi.h"
#include <stdio.h>
#include <math.h>

/* module load gcc */
/* module load mpich/ge/gcc/64/3.1 */
/* Compile me with 'mpicc -o cpi-parallel cpi-parallel.c' */


double f( double );
double f( double a )
{
    return (4.0 / (1.0 + a*a));
}

int main( int argc, char *argv[])
{
    int done = 0, n, myid, numprocs, i;
    double PI25DT = 3.141592653589793238462643;
    double mypi, pi, h, sum, x;
    double startwtime = 0.0, endwtime;
    int  namelen;
    char processor_name[MPI_MAX_PROCESSOR_NAME];

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

    fprintf(stderr,"Process %d on %s\n",
	    myid, processor_name);

    n = 0;
    while (!done)
    {
        if (myid == 0)
        {
/*
            printf("Enter the number of intervals: (0 quits) ");
            scanf("%d",&n);
*/
	    if (n==0) n=1024*numprocs; else n=0;

	    startwtime = MPI_Wtime();
        }
        MPI_Bcast(&n, 1, MPI_INT, 0, MPI_COMM_WORLD);
        if (n == 0)
            done = 1;
        else
        {
            h   = 1.0 / (double) n;
            sum = 0.0;
            for (i = myid + 1; i <= n; i += numprocs)
            {
                x = h * ((double)i - 0.5);
                sum += f(x);
            }
            mypi = h * sum;

            MPI_Reduce(&mypi, &pi, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);

            if (myid == 0)
	    {
                printf("pi is approximately %.16f, Error is %.16f\n",
                       pi, fabs(pi - PI25DT));
		endwtime = MPI_Wtime();
		printf("wall clock time = %f\n",
		       endwtime-startwtime);
	    }
        }
    }
    MPI_Finalize();

    return 0;
}

cpi_serial.c

#include <stdio.h>
#include <math.h>
#include <time.h>
#include <sys/time.h>
#include <stdint.h>

/* module load gcc */
/* Compile me with 'gcc -o cpi-serial cpi-serial.c' */

double f( double );
double f( double a )
{
    return (4.0 / (1.0 + a*a));
}

int main( int argc, char *argv[])
{
    int done = 0, n, i;
    double PI25DT = 3.141592653589793238462643;
    double mypi, pi, h, sum, x;
    double startwtime = 0.0, endwtime=0.0;
    struct timeval tv;

    n = 0;
    while (!done)
     {
      printf("Enter the number of intervals: (0 quits) ");
      scanf("%d",&n);
      
        if (n!=0) 
      {
      if (gettimeofday(&tv, NULL) == 0)
          startwtime = tv.tv_sec*1000000+tv.tv_usec;
      else
        startwtime = 0;
          
            h   = 1.0 / (double) n;
            sum = 0.0;
            for (i = 1; i <= n; i += 1)
            {
                x = h * ((double)i - 0.5);
                sum += f(x);
            }
            mypi = h * sum;

            printf("pi is approximately %.16f, Error is %.16f\n Number of intervals is %5d\n",
               mypi, fabs(mypi - PI25DT),n);
      
      if (gettimeofday(&tv, NULL) == 0)
          endwtime = tv.tv_sec*1000000+tv.tv_usec;
      else
        endwtime = 0;
    
        printf("wall clock time = %.0f usecs \n",
               (endwtime-startwtime));
       }
       else

    return 0;
    }
}

fpi_parallel.f

c**********************************************************************
c   From the mpich distribution
c   pi.f - compute pi by integrating f(x) = 4/(1 + x**2)     
c
c  module load intel
c  module load impi
c  Compile me with 'mpiifort -o fpi-parallel fpi-parallel.f'
c
c  (C) 2001 by Argonne National Laboratory.
c      See COPYRIGHT in top-level directory.
c     
c   Each node: 
c    1) receives the number of rectangles used in the approximation.
c    2) calculates the areas of it's rectangles.
c    3) Synchronizes for a global summation.
c   Node 0 prints the result.
c
c  Variables:
c
c    pi  the calculated result
c    n   number of points of integration.  
c    x           midpoint of each rectangle's interval
c    f           function to integrate
c    sum,pi      area of rectangles
c    tmp         temporary scratch space for global summation
c    i           do loop index
c****************************************************************************
      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, rc
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 )
      print *, "Process ", myid, " of ", numprocs, " is alive"

      sizetype   = 1
      sumtype    = 2
      
 10   if ( myid .eq. 0 ) then
         write(6,98)
 98      format('Enter the number of intervals: (0 quits)')
         read(5,99) n
 99      format(i10)
      endif
      
      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
         write(6, 97) pi, abs(pi - PI25DT)
 97      format('  pi is approximately: ', F18.16,
     +          '  Error is: ', F18.16)
      endif

      goto 10

 30   call MPI_FINALIZE(rc)
      stop
      end

fpi_serial.f

c
c Example program to calculate the value of pi by 
c integrating f(x) = 4 / (1 + x^2).
c
c  module load intel
c  Compile me with 'ifort -o fpi-serial fpi-serial.f'

c Compile me with 'g77 -o fpi-serial fpi-serial.f'
      program main
      double precision  PI25DT
      parameter        (PI25DT = 3.141592653589793238462643d0)
      double precision  mypi, pi, h, sum, x, f, a
      integer*4 num_iters
      integer rank, size, i, rc
      character junk
c     Loop until finished
      num_iters = 200000
c     Calculate the interval size
      h = 1.0d0 / num_iters
      sum  = 0.0d0
      do 10 i = 1, num_iters
         x = h * (dble(i) - 0.5d0)
         sum = sum + 4.d0 / (1.d0 + x * x)
 10   continue
      mypi = h * sum
c     All finished
      write(6, 97) num_iters, mypi, abs(mypi - PI25DT)
 97   format(i20, ' points: pi is approximately: ', F18.16,
     +     ' error is: ', F18.16)
      stop
      end

it's good