Public feed
#529145
/* Author: Cyle Riggs 3mail: cyle !!at!! riggs family web %d0t% net*/#include <iostream>#include <iomanip>#include <cstdlib>#include <cmath>#include <map>#include <vector>#include "mpi.h"#define TAG_POINTS 11#define TAG_INCIRCLE 12#define REALPI 3.14159265358979323846264338327950288419716939937510using namespace std;struct CPoint { double x; double y;};struct CRunInfo { double index; double error; double pi; double points; double wtime;};double drand( double low, double high );int main( int argc, char* argv[] ) { MPI_Status status; MPI_Init( &argc, &argv ); int size = -1; int rank = -1; MPI_Comm_size( MPI_COMM_WORLD, &size ); MPI_Comm_rank( MPI_COMM_WORLD, &rank ); srand( (rank+1) * time(NULL) + clock() ); unsigned long int done = 0; double maxCount = 0; if( rank == 0 ) { cout << "max # of points? " << flush; cin >> maxCount; } MPI_Bcast( &maxCount, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD ); //If the user chose 0 quit here if( maxCount == 0 ) { MPI_Finalize(); return 0; } vector<CRunInfo> run_results; for( double count = 10; count < maxCount; count++ ) { //Bcast 'count' to stay in sync MPI_Bcast( &count, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD ); double start_time = MPI_Wtime(); map<unsigned long int,unsigned long int> job_map; for( unsigned long int i = 0; i < count; i++ ) { job_map[i%size]++; } if( rank == 0 ) { cout << "Generating ~" << job_map[rank] << " random points per process" << endl; } vector<CPoint> mypoints; for( unsigned long int i = 0; i < job_map[rank]; i++ ) { CPoint p; p.x = drand( -1.0, 1.0 ); //Generate X for the point p.y = drand( -1.0, 1.0 ); //Generate Y for the point mypoints.push_back( p ); } double my_in_circle = 0.0; for( vector<CPoint>::iterator point = mypoints.begin(); point != mypoints.end(); point++ ) { if( (pow(point->x,2)+pow(point->y,2)) < 1 ) { my_in_circle++; } } double final_in_circle = 0.0; //Add together everybody's count of points in the circle MPI_Reduce( &my_in_circle, &final_in_circle, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD ); if( rank == 0 ) { cout << "Estimating pi..........." << flush; double est_pi = (4.0 * (final_in_circle/count)); cout << "OK" << endl; cout << "Checking accuracy......." << flush; double error = est_pi - REALPI; if( error < 0 ) { error *= -1; } cout << "OK" << endl; double end_time = MPI_Wtime(); //PI is the 4 * (# of points that fell in circle / the total number of points) cout << setprecision(50); cout << endl << endl; cout << "Runtime: " << end_time - start_time << endl; cout << "Count: " << count << endl; cout << "Pi: " << est_pi << endl; cout << "Error: " << error << endl; CRunInfo run_info; run_info.wtime = end_time - start_time; run_info.pi = est_pi; run_info.points = count; run_info.error = error; run_results.push_back( run_info ); } } if( rank == 0 ) { double best_error = 99999999; unsigned long int best_index = 0; for( unsigned long int i = 0; i < run_results.size(); i++ ) { if( run_results[i].error < best_error ) { cout << "DEBUG: " << run_results[i].error << " < " << best_error << endl; best_index = i; best_error = run_results[i].error; } } cout.setf( ios::fixed ); cout << endl << endl << endl << endl; cout << "Most accurate run was run #" << best_index + 1 << ":" << endl; cout << "Error: " << run_results[best_index].error << endl; cout << "Pi: " << run_results[best_index].pi << endl; cout << "# of points: " << run_results[best_index].points << endl; } MPI_Finalize(); return 0;}double drand( double low, double high ) { double retval = static_cast<double>( rand() ) / RAND_MAX; return low + ((high - low) * retval);}