2015-11-14 132 views
-1

我正在嘗試使用montecarlo方法和使用並行C代碼來查找PI的值。我已經寫了serail代碼並且工作正常。但是,並行代碼給我一些時間0或PI錯誤的值負值使用蒙特卡洛方法多線程計算pi值

我的代碼

#include <pthread.h> 
#include <stdio.h> 
#include <stdlib.h> 
#include <math.h> 
#define NUM_THREADS 4 //number of threads 
#define TOT_COUNT 10000055 //total number of iterations 



void *doCalcs(void *threadid) 
{ 
long longTid; 
longTid = (long)threadid; 
int tid = (int)longTid; //obtain the integer value of thread id 
//using malloc for the return variable in order make 
//sure that it is not destroyed once the thread call is finished 
float *in_count = (float *)malloc(sizeof(float)); 
*in_count=0; 
unsigned int rand_state = rand(); 
//get the total number of iterations for a thread 
float tot_iterations= TOT_COUNT/NUM_THREADS; 
int counter=0; 
//calculation 
for(counter=0;counter<tot_iterations;counter++){ 
    //float x = (double)random()/RAND_MAX; 
    //float y = (double)random()/RAND_MAX; 
    //float result = sqrt((x*x) + (y*y)); 
    double x = rand_r(&rand_state)/((double)RAND_MAX + 1) * 2.0 - 1.0; 
    double y = rand_r(&rand_state)/((double)RAND_MAX + 1) * 2.0 - 1.0; 
    float result = sqrt((x*x) + (y*y)); 
    if(result<1){ 
     *in_count+=1; //check if the generated value is inside a unit circle 
    } 
} 
//get the remaining iterations calculated by thread 0 
if(tid==0){ 
     float remainder = TOT_COUNT%NUM_THREADS; 
     for(counter=0;counter<remainder;counter++){ 
      float x = (double)random()/RAND_MAX; 
     float y = (double)random()/RAND_MAX; 
      float result = sqrt((x*x) + (y*y)); 
     if(result<1){ 
      *in_count+=1; //check if the generated value is inside a unit circle 
     } 
    } 
} 
} 


int main(int argc, char *argv[]) 
    { 
    pthread_t threads[NUM_THREADS]; 
    int rc; 
    long t; 
    void *status; 
    float tot_in=0; 
    for(t=0;t<NUM_THREADS;t++){ 
     rc = pthread_create(&threads[t], NULL, doCalcs, (void *)t); 
     if (rc){ 
    printf("ERROR; return code from pthread_create() is %d\n", rc); 
    exit(-1); 
    } 
} 
//join the threads 
for(t=0;t<NUM_THREADS;t++){ 
pthread_join(threads[t], &status); 
//printf("Return from thread %ld is : %f\n",t, *(float*)status); 
tot_in+=*(float*)status; //keep track of the total in count 
} 
printf("Value for PI is %f \n",1, 4*(tot_in/TOT_COUNT)); 
/* Last thing that main() should do */ 
pthread_exit(NULL); 
    }    
+0

該代碼中有許多可疑的事情,比如爲什麼你爲'in_count'動態地分配一個浮點數?你順便把它看作一個整數。哦,而且你並沒有在任何地方釋放導致內存泄漏的配置。 –

+0

我建議使用'double'。如果你想'浮動',爲什麼不計算'22/7'或'355/113'? –

+0

@WeatherVane他想嘗試蒙特卡洛方法,我認爲教育。 – vladon

回答

0

你的代碼是不是C++,這是不好的,非常不好的老式C.

也就是說C++:

#include <cmath> 
#include <iostream> 
#include <numeric> 
#include <random> 
#include <thread> 
#include <vector> 

constexpr auto num_threads = 4; //number of threads 
constexpr auto total_count = 10000055; //total number of iterations 

void doCalcs(int total_iterations, int & in_count_result) 
{ 
    auto seed = std::random_device{}(); 
    auto gen = std::mt19937{ seed }; 
    auto dist = std::uniform_real_distribution<>{0, 1}; 

    auto in_count{ 0 }; 

    //calculation 
    for (auto counter = 0; counter < total_iterations; ++counter) { 
     auto x = dist(gen); 
     auto y = dist(gen); 
     auto result = std::sqrt(std::pow(x, 2) + std::pow(y, 2)); 
     if (result < 1) { 
      ++in_count; //check if the generated value is inside a unit circle 
     } 
    } 

    in_count_result = in_count; 
} 

void main() 
{ 
    std::vector<std::thread> threads(num_threads); 
    std::vector<int> in_count(num_threads); 
    in_count.resize(num_threads); 

    for (size_t i = 0; i < num_threads; ++i) { 
     int total_iterations = total_count/num_threads; 
     if (i == 0) { 
      total_iterations += total_count % num_threads; // get the remaining iterations calculated by thread 0 
     } 

     threads.emplace_back(doCalcs, total_iterations, std::ref(in_count[i])); 
    } 

    for (auto & thread : threads) { 
     if (thread.joinable()) { 
      thread.join(); 
     } 
    } 

    double pi_value = 4.0 * static_cast<double>(std::accumulate(in_count.begin(), in_count.end(), 0))/static_cast<double>(total_count); 
    std::cout << "Value of PI is: " << pi_value << std::endl; 
} 

PS這也不是那麼好,請閱讀future s,promise s和std::async