/* ****************************************************************************
 *
 *   ompipi.c -- OpenMP/MPI hybrid demo program
 *
 *   Computes pi by integration -- simple Riemann sum with midpoint heights.
 *   Accepts sub-intervals on command line, or sets an appropriate default.
 *
 *   djames (at) tacc.utexas.edu
 *   10 Apr 2014
 *   OpenMP hybrid version by Steve Lantz, Cornell CAC, Apr 2017 - Sept 2024
 *
 * ****************************************************************************
 */

#include <mpi.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <omp.h>
#define N_TRIALS 10

const int    DEFAULT_SUBINTS = 1000000;            // default sub-intervals
const int    ROOT_RANK       = 0;                  // rank of root process
const double ACTUAL_PI  = 3.14159265358979323846;  // for error calculations

double f( double x );
int    getSubints( int argc, char* argv[] );
void   reportResults( double computedPi, int totalTasks, int nThreads,
                      int totalSubints, double firstTime, double bestTime,
                      double tickTime );

// ----------------------------------------------------------------------------

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

    int totalTasks;     // total MPI tasks (processes)
    int myRank;         // rank number (process ID) of current process
    int totalSubints;   // total number of sub-intervals in Riemann sum
    int statusCode;     // status code returned by certain functions
    int j;              // loop index; counts timing trials
    int i;              // loop index; counts through sub-intervals
    double xVal;        // x coordinate at which we're evaluating function
    double myPortion;   // portion of integral computed by this process
    double computedPi;  // integral (area) summed over all processes
    int nThreads;       // number of threads per task
    int rqd, pvd;       // threading levels required of and provided by MPI

    double startTime, stopTime, elapsedTime[N_TRIALS], tickTime;
    double minTime = 9999.0;
    
    rqd = MPI_THREAD_FUNNELED;
    MPI_Init_thread( &argc, &argv, rqd, &pvd );
    MPI_Comm_size( MPI_COMM_WORLD, &totalTasks );
    MPI_Comm_rank( MPI_COMM_WORLD, &myRank );

    if ( myRank == ROOT_RANK ) {  // root process
        startTime = MPI_Wtime();
        totalSubints = getSubints( argc, argv );
    }  // end if root process

    // Divide the labor across tasks (and yes -- the root process
    // does its share of the computation).  Here we're using stride
    // to distribute the tasks...

    statusCode   // return value ignored
      = MPI_Bcast( &totalSubints, 1, MPI_INTEGER, ROOT_RANK, MPI_COMM_WORLD );

    double width = 1.0 / ( (double)( totalSubints ) );
    double halfWidth = 0.5*width;

    for ( j=0; j<N_TRIALS; j++ ) {
        myPortion = 0.0;

        #ifdef _OPENMP
        #pragma omp parallel for private( xVal ) reduction( +:myPortion )
        #endif
        for ( i=myRank; i<totalSubints; i+=totalTasks ) {
            // DON'T just increment xVal, since roundoff error accumulates
            xVal = (double)(i)*width + halfWidth;
            myPortion += f( xVal );
        }  // for each Riemann rectangle i

        if ( myRank == ROOT_RANK ) {  // root process
            stopTime = MPI_Wtime();
            elapsedTime[j] = stopTime - startTime;
            if ( minTime > elapsedTime[j] ) {
                minTime = elapsedTime[j];
            }  // end if minTime
            startTime = stopTime;
        }  // end if root process

        myPortion *= width;

        // Assemble the results into a final result...

        computedPi = 0.0;
        statusCode  // return value ignored
          = MPI_Reduce( &myPortion, &computedPi, 1,
                         MPI_DOUBLE, MPI_SUM, ROOT_RANK, MPI_COMM_WORLD );

    }  // for each trial j

    // Report the results...

    if ( myRank == ROOT_RANK ) {  // root process
        nThreads = 1;
        #ifdef _OPENMP
        nThreads = omp_get_max_threads();
        #endif
        tickTime = MPI_Wtick();
        reportResults( computedPi, totalTasks, nThreads,
                       totalSubints, elapsedTime[0], minTime,
                       tickTime );
    }  // end if root process

    MPI_Finalize();
    return 0;

}  // end main function

// ----------------------------------------------------------------------------

int getSubints( int argc, char* argv[] ) {

    int totalSubints = DEFAULT_SUBINTS;
      // the total number of sub-intervals as specified by user

    if ( argc > 1 )  {
        totalSubints = atoi( argv[1] );  // no input validation
    }  // end if no args else clause

    return totalSubints;

}  // end getSubints function

// ----------------------------------------------------------------------------

void reportResults( double computedPi, int totalTasks, int nThreads,
                    int totalSubints, double firstTime, double bestTime,
                    double tickTime ) {

    double absError = fabs(  computedPi - ACTUAL_PI );
    double relError = absError / ACTUAL_PI;

    printf( "\n\t**********************************************\n"    );
    printf(   "\t         Computed value: %17.15f\n",    computedPi   );
    printf(   "\t         Absolute error: %7.2e\n",      absError     );
    printf(   "\t         Relative error: %7.2e\n",      relError     );
    printf(   "\t    Total MPI processes: %d\n",         totalTasks   );
    printf(   "\t    Threads per process: %d\n",         nThreads     );
    printf(   "\t    Total sub-intervals: %d\n",         totalSubints );
    printf(   "\t             First time: %7.2e secs\n", firstTime    );
    printf(   "\t              Best time: %7.2e secs\n", bestTime     );
    printf(   "\t        Time resolution: %7.2e secs\n", tickTime     );
    printf(   "\t**********************************************\n\n"  );

}  // end reportResults function

// ----------------------------------------------------------------------------

double f( double x ) {
    double value = 4.0 / ( 1.0 + x*x );
    return value;
}  // end function f