# include <cstdlib>
# include <iostream>
# include <fstream>
# include <cmath>

using namespace std;

# include "mpi.h"

int main ( int argc, char *argv[] );
double boundary_condition ( double x, double time );
double initial_condition ( double x, double time );
double rhs ( double x, double time );
void update ( int id, int p );

//****************************************************************************80

int main ( int argc, char *argv[] )

//****************************************************************************80
//
//  Purpose:
//
//    MAIN is the main program for HEAT_MPI.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 April 2008
//
//  Author:
// 
//    John Burkardt
//
//  Reference:
//
//    William Gropp, Ewing Lusk, Anthony Skjellum,
//    Using MPI: Portable Parallel Programming with the
//    Message-Passing Interface,
//    Second Edition,
//    MIT Press, 1999,
//    ISBN: 0262571323,
//    LC: QA76.642.G76.
//
//    Marc Snir, Steve Otto, Steven Huss-Lederman, David Walker, 
//    Jack Dongarra,
//    MPI: The Complete Reference,
//    Volume I: The MPI Core,
//    Second Edition,
//    MIT Press, 1998,
//    ISBN: 0-262-69216-3,
//     LC: QA76.642.M65.
//
{
  int id;
  int p;
  double wtime;

  cout << "\n";
  cout << "HEAT_MPI:\n";
  cout << "  C++/MPI version\n";
  cout << "\n";
  cout << "  Solve the 1D time-dependent heat equation.\n";

  MPI::Init ( argc, argv );

  id = MPI::COMM_WORLD.Get_rank ( );

  p = MPI::COMM_WORLD.Get_size ( );
//
//  Record the starting time.
//
  if ( id == 0 ) 
  {
    wtime = MPI::Wtime ( );
  }

  update ( id, p );
//
//  Record the final time.
//
  if ( id == 0 )
  {
    wtime = MPI::Wtime ( ) - wtime;

    cout << "\n";       
    cout << "  Wall clock elapsed seconds = " << wtime << "\n";      
  }

  MPI::Finalize ( );

  cout << "\n";
  cout << "HEAT_MPI:\n";
  cout << "  Normal end of execution.\n";

  return 0;
}
//****************************************************************************80

void update ( int id, int p )

//****************************************************************************80
//
//  Purpose:
//
//    UPDATE computes the solution of the heat equation.
//
//  Discussion:
//
//    If there is only one processor ( P == 1 ), then the program writes the
//    values of X and H to files.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 April 2008
//
//  Author:
// 
//    John Burkardt
//
//  Parameters:
//
//    Input, int ID, the id of this processor.
//
//    Input, int P, the number of processors.
//
{
  double cfl;
  double *h;
  ofstream h_file;
  double *h_new;
  int i;
  int j;
  int j_min = 0;
  int j_max = 100;
  double k = 0.002;
  int n = 11;
  MPI::Status status;
  int tag;
  double time;
  double time_delta;
  double time_max = 10.0;
  double time_min = 0.0;
  double time_new;
  double *x;
  double x_delta;
  ofstream x_file;
  double x_max = 1.0;
  double x_min = 0.0;
//
//  Have process 0 print out some information.
//
  if ( id == 0 )
  {
    cout << "\n";
    cout << "  Compute an approximate solution to the time dependent\n";
    cout << "  one dimensional heat equation:\n";
    cout << "\n";
    cout << "    dH/dt - K * d2H/dx2 = f(x,t)\n";
    cout << "\n";
    cout << "  for " << x_min << " = x_min < x < x_max = " << x_max << "\n";
    cout << "\n";
    cout << "  and " << time_min << " = time_min < t <= t_max = " << time_max << "\n";
    cout << "\n";
    cout << "  Boundary conditions are specified at x_min and x_max.\n";
    cout << "  Initial conditions are specified at time_min.\n";
    cout << "\n";
    cout << "  The finite difference method is used to discretize the\n";
    cout << "  differential equation.\n";
    cout << "\n";
    cout << "  This uses " << p * n << " equally spaced points in X\n";
    cout << "  and " << j_max << " equally spaced points in time.\n";
    cout << "\n";
    cout << "  Parallel execution is done using " << p << " processors.\n";
    cout << "  Domain decomposition is used.\n";
    cout << "  Each processor works on " << n << " nodes, \n";
    cout << "  and shares some information with its immediate neighbors.\n";
  }
//
//  Set the X coordinates of the N nodes.
//  We don't actually need ghost values of X but we'll throw them in
//  as X[0] and X[N+1].
//
  x = new double[n+2];

  for ( i = 0; i <= n + 1; i++ )
  {
    x[i] = ( ( double ) (         id * n + i - 1 ) * x_max
           + ( double ) ( p * n - id * n - i     ) * x_min )
           / ( double ) ( p * n              - 1 );
  }
//
//  In single processor mode, write out the X coordinates for display.
//
  if ( p == 1 )
  {
    x_file.open ( "x_data.txt" );
    for ( i = 1; i <= n; i++ )
    {
      x_file << "  %f", x[i];
    }
    x_file << "\n";

    x_file.close ( );
  }
//
//  Set the values of H at the initial time.
//
  time = time_min;
  h = new double[n+2];
  h_new = new double[n+2];
  h[0] = 0.0;
  for ( i = 1; i <= n; i++ )
  {
    h[i] = initial_condition ( x[i], time );
  }
  h[n+1] = 0.0;
  
  time_delta = ( time_max - time_min ) / ( double ) ( j_max - j_min );
  x_delta = ( x_max - x_min ) / ( double ) ( p * n - 1 );
//
//  Check the CFL condition, have processor 0 print out its value,
//  and quit if it is too large.
//
  cfl = k * time_delta / x_delta / x_delta;

  if ( id == 0 ) 
  {
    cout << "\n";
    cout << "UPDATE\n";
    cout << "  CFL stability criterion value = " << cfl << "\n";;
  }

  if ( 0.5 <= cfl ) 
  {
    if ( id == 0 )
    {
      cout << "\n";
      cout << "UPDATE - Warning!\n";
      cout << "  Computation cancelled!\n";
      cout << "  CFL condition failed.\n";
      cout << "  0.5 <= K * dT / dX / dX = " << cfl << "\n";
    }
    return;
  }
//
//  In single processor mode, write out the values of H.
//
  if ( p == 1 )
  {
    h_file.open ( "h_data.txt" );

    for ( i = 1; i <= n; i++ )
    {
      h_file << "  " << h[i];
    }
    h_file << "\n";
  }
//
//  Compute the values of H at the next time, based on current data.
//
  for ( j = 1; j <= j_max; j++ )
  {

    time_new = ( ( double ) (         j - j_min ) * time_max
               + ( double ) ( j_max - j         ) * time_min )
               / ( double ) ( j_max     - j_min );
//
//  Send H[1] to ID-1.
//
    if ( 0 < id )
    {
      tag = 1;
      MPI::COMM_WORLD.Send ( &h[1], 1, MPI::DOUBLE, id-1, tag );
    }
//
//  Receive H[N+1] from ID+1.
//
    if ( id < p-1 )
    {
      tag = 1;
      MPI::COMM_WORLD.Recv ( &h[n+1], 1,  MPI::DOUBLE, id+1, tag, status );
    }
//
//  Send H[N] to ID+1.
//
    if ( id < p-1 )
    {
      tag = 2;
      MPI::COMM_WORLD.Send ( &h[n], 1, MPI::DOUBLE, id+1, tag );
    }
//
//  Receive H[0] from ID-1.
//
    if ( 0 < id )
    {
      tag = 2;
      MPI::COMM_WORLD.Recv ( &h[0], 1, MPI::DOUBLE, id-1, tag, status );
    }
//
//  Update the temperature based on the four point stencil.
//
    for ( i = 1; i <= n; i++ )
    {
      h_new[i] = h[i] 
      + ( time_delta * k / x_delta / x_delta ) * ( h[i-1] - 2.0 * h[i] + h[i+1] ) 
      + time_delta * rhs ( x[i], time );
    }
//
//  H at the extreme left and right boundaries was incorrectly computed
//  using the differential equation.  Replace that calculation by
//  the boundary conditions.
//
    if ( 0 == id )
    {
      h_new[1] = boundary_condition ( x[1], time_new );
    }
    if ( id == p - 1 )
    {
      h_new[n] = boundary_condition ( x[n], time_new );
    }
//
//  Update time and temperature.
//
    time = time_new;

    for ( i = 1; i <= n; i++ )
    {
      h[i] = h_new[i];
    }
//
//  In single processor mode, add current solution data to output file.
//
    if ( p == 1 )
    {
      for ( i = 1; i <= n; i++ )
      {
        h_file << "  " << h[i];
      }
      h_file << "\n";
    }
  }

  if ( p == 1 )
  {
    h_file.close ( );
  }

  delete [] h;
  delete [] h_new;
  delete [] x;

  return;
}
//****************************************************************************80

double boundary_condition ( double x, double time )

//****************************************************************************80
//
//  Purpose:
//
//    BOUNDARY_CONDITION evaluates the boundary condition of the differential equation.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 April 2008
//
//  Author:
// 
//    John Burkardt
//
//  Parameters:
//
//    Input, double X, TIME, the position and time.
//
//    Output, double BOUNDARY_CONDITION, the value of the boundary condition.
//
{
  double value;
//
//  Left condition:
//
  if ( x < 0.5 )
  {
    value = 100.0 + 10.0 * sin ( time );
  }
  else
  {
    value = 75.0;
  }
  return value;
}
//****************************************************************************80

double initial_condition ( double x, double time )

//****************************************************************************80
//
//  Purpose:
//
//    INITIAL_CONDITION evaluates the initial condition of the differential equation.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 April 2008
//
//  Author:
// 
//    John Burkardt
//
//  Parameters:
//
//    Input, double X, TIME, the position and time.
//
//    Output, double INITIAL_CONDITION, the value of the initial condition.
//
{
  double value;

  value = 95.0;

  return value;
}
//****************************************************************************80

double rhs ( double x, double time )

//****************************************************************************80
//
//  Purpose:
//
//    RHS evaluates the right hand side of the differential equation.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 April 2008
//
//  Author:
// 
//    John Burkardt
//
//  Parameters:
//
//    Input, double X, TIME, the position and time.
//
//    Output, double RHS, the value of the right hand side function.
//
{
  double value;

  value = 0.0;

  return value;
}

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

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=100; 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;
}

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

int main(int argc, char *argv[])
{
   int rank, numprocs;

   MPI_Init(&argc, &argv);
   MPI_Comm_size(MPI_COMM_WORLD,&numprocs);
   MPI_Comm_rank(MPI_COMM_WORLD,&rank);
   printf("Hello World from proces %d of %d.\n",rank,numprocs);
   MPI_Finalize();
}

program hworld
USE MPI

character*(MPI_MAX_PROCESSOR_NAME) host
integer ncores,rank,ierr,hostname_len

call MPI_INIT(ierr)
call MPI_COMM_SIZE(MPI_COMM_WORLD,ncores,ierr)
call MPI_COMM_RANK(MPI_COMM_WORLD,rank,ierr)
call MPI_GET_PROCESSOR_NAME(host,hostname_len,ierr)
write (6,'(''  Hello from mpif90 on host '' , a10, '', core '',i5, '' of '',i5,'' cores'')') &
        host,rank,ncores
call MPI_FINALIZE(ierr)

stop
end

c**********************************************************************
c     This program calculates the value of pi, using numerical integration
c     with parallel processing.  The user selects the number of points of
c     integration.  By selecting more points you get more accurate results
c     at the expense of additional computation
c
c     This version is written using p4 calls to handle message passing
c     It should run without changes on most workstation clusters and MPPs.
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     Constants:
c
c     SIZETYPE    initial message to the cube
c     ALLNODES    used to load all nodes in cube with a node process
c     INTSIZ      four bytes for an integer
c     DBLSIZ      eight bytes for double precision
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)

      integer   INTSIZ , DBLSIZ,  ALLNODES,   ANYNODE
      parameter(INTSIZ=4,DBLSIZ=8,ALLNODES=-1,ANYNODE=-1)

      double precision  pi, h, sum, x, f, a, temp
      integer n, myid, numnodes, i, rc
      integer sumtype, sizetype, masternode
      integer status(3)

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, numnodes, ierr )
c     print *, "Process ", myid, " of ", numnodes, " is alive"

      sizetype   = 10
      sumtype    = 17
      masternode =  0

 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)

         do i=1,numnodes-1
            call MPI_SEND(n,1,MPI_INTEGER,i,sizetype,MPI_COMM_WORLD,rc)
         enddo

      else

         call MPI_RECV(n,1,MPI_INTEGER,masternode,sizetype,
     +        MPI_COMM_WORLD,status,rc)

      endif

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, numnodes
         x = h * (dble(i) - 0.5d0)
         sum = sum + f(x)
 20         continue
      pi = h * sum

      if (myid .ne. 0) then

         call MPI_SEND(pi,1,MPI_DOUBLE_PRECISION,masternode,sumtype,
     +        MPI_COMM_WORLD,rc)

      else

         do i=1,numnodes-1
            call MPI_RECV(temp,1,MPI_DOUBLE_PRECISION,i,sumtype,
     +           MPI_COMM_WORLD,status,rc)
            pi = pi + temp
         enddo
      endif

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)
      end

 pi.f