|
| 1 | +// Author: Wes Kendall |
| 2 | +// Copyright 2012 www.mpitutorial.com |
| 3 | +// This code is provided freely with the tutorials on mpitutorial.com. Feel |
| 4 | +// free to modify it for your own use. Any distribution of the code must |
| 5 | +// either provide a link to www.mpitutorial.com or keep this header in tact. |
| 6 | +// |
| 7 | +// Program that computes the average of an array of elements in parallel using |
| 8 | +// MPI_Scatter and MPI_Gather |
| 9 | +// |
| 10 | +#include <stdio.h> |
| 11 | +#include <stdlib.h> |
| 12 | +#include <mpi.h> |
| 13 | +#include <assert.h> |
| 14 | +// #include <math.h> |
| 15 | +#include <cmath> |
| 16 | + |
| 17 | +#include <complex> |
| 18 | +#include <boost/numeric/ublas/matrix.hpp> |
| 19 | +#include <boost/numeric/ublas/io.hpp> |
| 20 | +#include <bitset> |
| 21 | +#include <typeinfo> |
| 22 | + |
| 23 | +double pi; |
| 24 | +std::complex<double> I; |
| 25 | +namespace ublas = boost::numeric::ublas; |
| 26 | + |
| 27 | + |
| 28 | +// Creates an array of random numbers. Each number has a value from 0 - 1 |
| 29 | +double *create_rand_nums(int num_elements) { |
| 30 | + double *rand_nums = (double *)malloc(sizeof(double) * num_elements); |
| 31 | + assert(rand_nums != NULL); |
| 32 | + int i; |
| 33 | + for (i = 0; i < num_elements; i++) { |
| 34 | + rand_nums[i] = (rand() / (double)RAND_MAX); |
| 35 | + } |
| 36 | + return rand_nums; |
| 37 | +} |
| 38 | + |
| 39 | + |
| 40 | +double *create_nums(int num_elements,int x) { |
| 41 | + double *nums = (double *)malloc(sizeof(double) * num_elements); |
| 42 | + assert(nums != NULL); |
| 43 | + int i; |
| 44 | + for (i = 0; i < num_elements; i++) { |
| 45 | + nums[i] = x; |
| 46 | + printf("create num:%f\n",nums[i] ); |
| 47 | + } |
| 48 | + return nums; |
| 49 | +} |
| 50 | + |
| 51 | +double *create_onelinephases(double arr[],int size) { |
| 52 | + double *nums = (double *)malloc(sizeof(double) * size); |
| 53 | + assert(nums != NULL); |
| 54 | + int i; |
| 55 | + for (i = 0; i < size; i++) { |
| 56 | + nums[i] = arr[i]; |
| 57 | + } |
| 58 | + return nums; |
| 59 | +} |
| 60 | +double *compute_qft(double *array, int n_qft_point, int n) { |
| 61 | + std::complex<double> *nums_complex = (std::complex<double> *)malloc(sizeof(std::complex<double>) * n_qft_point); |
| 62 | + |
| 63 | + int i; |
| 64 | + for (i = 0; i < n_qft_point; i++) |
| 65 | + { |
| 66 | + nums_complex[i] = array[i]+n;//exp(2*pi*array[i]); |
| 67 | + //printf("world_rank: with number %f\n",nums[i] ); |
| 68 | + } |
| 69 | + |
| 70 | + |
| 71 | + double *nums = (double *)malloc(sizeof(double) * n_qft_point); |
| 72 | + assert(nums != NULL); |
| 73 | + ublas::vector<std::complex<double> > v1(4), v2(2); |
| 74 | + |
| 75 | + for (i = 0; i < n_qft_point; i++) { |
| 76 | + |
| 77 | + nums[i] = array[i]+n; |
| 78 | + //printf("world_rank: %d with number %f\n",n,nums[i] ); |
| 79 | + } |
| 80 | + return nums; |
| 81 | +} |
| 82 | +// Computes the average of an array of numbers |
| 83 | +double compute_avg(double *array, int num_elements) { |
| 84 | + double sum = 0.f; |
| 85 | + int i; |
| 86 | + for (i = 0; i < num_elements; i++) { |
| 87 | + sum += array[i]; |
| 88 | + } |
| 89 | + return sum / num_elements; |
| 90 | +} |
| 91 | +// Computes the average of an array of numbers |
| 92 | +double *compute_neg(double *array, int num_elements,int n) { |
| 93 | + double *nums = (double *)malloc(sizeof(double) * num_elements); |
| 94 | + assert(nums != NULL); |
| 95 | + ublas::vector<std::complex<double> > v1(4), v2(2); |
| 96 | + int i; |
| 97 | + for (i = 0; i < num_elements; i++) { |
| 98 | + |
| 99 | + nums[i] = array[i]+n; |
| 100 | + //printf("world_rank: %d with number %f\n",n,nums[i] ); |
| 101 | + } |
| 102 | + return nums; |
| 103 | +} |
| 104 | + |
| 105 | + |
| 106 | +int main(int argc, char** argv) { |
| 107 | + |
| 108 | + |
| 109 | + |
| 110 | + |
| 111 | + |
| 112 | + |
| 113 | + // std::string binary = std::bitset<8>(pow(2,3)).to_string(); //to binary |
| 114 | + // std::cout<<binary[0]<<"\n"; |
| 115 | + |
| 116 | + // unsigned long decimal = std::bitset<8>(binary).to_ulong(); |
| 117 | + // std::cout<<decimal<<"\n"; |
| 118 | + |
| 119 | + |
| 120 | + |
| 121 | + |
| 122 | + |
| 123 | + |
| 124 | + std::complex<double> I(0, 1); |
| 125 | + pi = 4 * atan(1.0); |
| 126 | + // std::complex<double> e(-2*pi/2,0); |
| 127 | + // std::complex<double> z = exp(I*e); |
| 128 | + // double x = 2; |
| 129 | + // std::cout << real(z*2.0) <<'\n'; |
| 130 | + ublas::vector<std::complex<double> > v1(4), v2(2); |
| 131 | + // for (unsigned i = 0; i < 4; ++i) |
| 132 | + // v1 (i) = i; |
| 133 | + // for (unsigned i = 0; i < 2; ++i) |
| 134 | + // v2 (i) = i; |
| 135 | + |
| 136 | + // ublas::matrix<double> m(v1.size(), v2.size()); |
| 137 | + // std::cout |
| 138 | + // << v1 << '\n' |
| 139 | + // << v2 << '\n' |
| 140 | + // << outer_prod(v1, v2)<< '\n'; |
| 141 | + int num_elements_per_proc = 3;//atoi(argv[1]); |
| 142 | + |
| 143 | + // Seed the random number generator to get different results each time |
| 144 | + |
| 145 | + |
| 146 | + MPI_Init(NULL, NULL); |
| 147 | + |
| 148 | + int world_rank; |
| 149 | + MPI_Comm_rank(MPI_COMM_WORLD, &world_rank); |
| 150 | + int world_size; |
| 151 | + MPI_Comm_size(MPI_COMM_WORLD, &world_size); |
| 152 | + //printf("world_size: %d\n",world_size ); |
| 153 | + |
| 154 | + // Create a random array of elements on the root process. Its total |
| 155 | + // size will be the number of elements per process times the number |
| 156 | + // of processes |
| 157 | + |
| 158 | + int num_point_qft = 3; |
| 159 | + int num_states = pow(2,num_point_qft); |
| 160 | + double *rand_nums = NULL; |
| 161 | + if (world_rank == 0) { |
| 162 | + |
| 163 | + |
| 164 | + |
| 165 | + double final_const = 1/(pow(2,num_point_qft/2)); |
| 166 | + |
| 167 | + |
| 168 | + |
| 169 | + // int *decimal_bits = (int *)malloc(sizeof(int) * num_states); |
| 170 | + std::string binary_bits[num_states]; |
| 171 | + for (int i = 0; i < num_states; ++i) |
| 172 | + { |
| 173 | + binary_bits[i] = std::bitset<3>(i).to_string(); // thisssss <3> cannot be replaced with num_point_qft |
| 174 | + //std::cout<<binary_bits[i][0]<<"\n"; |
| 175 | + } |
| 176 | + double phases[num_states][num_point_qft]; |
| 177 | + |
| 178 | + for (int i = 0; i < num_states; ++i) |
| 179 | + { |
| 180 | + |
| 181 | + for (int j = 0; j < num_point_qft; ++j) |
| 182 | + { |
| 183 | + int start_checking_point = num_point_qft - 1 -j; |
| 184 | + for(int k=0; k < num_point_qft; ++k) |
| 185 | + { |
| 186 | + if(start_checking_point<num_point_qft && binary_bits[i][j]=='1') |
| 187 | + { |
| 188 | + phases[i][start_checking_point]+=pow(2,-(k+1)); |
| 189 | + } |
| 190 | + start_checking_point+=1; |
| 191 | + } |
| 192 | + } |
| 193 | + } |
| 194 | + |
| 195 | + double oneline_phases[num_states*num_point_qft]; |
| 196 | + int k=0; |
| 197 | + for (int i = 0; i < num_states; ++i) |
| 198 | + { |
| 199 | + for(int j=0; j< num_point_qft;++j) |
| 200 | + { |
| 201 | + //printf("%f\t",phases[i][j] ); |
| 202 | + oneline_phases[k]=phases[i][j]; |
| 203 | + k++; |
| 204 | + } |
| 205 | + //printf("\n"); |
| 206 | + } |
| 207 | + // printf("start oneline \n"); |
| 208 | + |
| 209 | + // for (int i = 0; i < num_point_qft*num_states; ++i) |
| 210 | + // { |
| 211 | + // printf("%f,",oneline_phases[i] ); |
| 212 | + // } |
| 213 | + // printf("end oneline \n"); |
| 214 | + |
| 215 | + |
| 216 | + // rand_nums = create_nums(num_elements_per_proc * world_size,world_rank); |
| 217 | + rand_nums = create_onelinephases(oneline_phases, num_point_qft*num_states); |
| 218 | + } |
| 219 | + |
| 220 | + // For each process, create a buffer that will hold a subset of the entire |
| 221 | + // array |
| 222 | + double *sub_rand_nums = (double *)malloc(sizeof(double) * num_elements_per_proc); |
| 223 | + assert(sub_rand_nums != NULL); |
| 224 | + |
| 225 | + // Scatter the random numbers from the root process to all processes in |
| 226 | + // the MPI world |
| 227 | + MPI_Scatter(rand_nums, num_elements_per_proc, MPI_DOUBLE, sub_rand_nums, |
| 228 | + num_elements_per_proc, MPI_DOUBLE, 0, MPI_COMM_WORLD); |
| 229 | + |
| 230 | + // Compute the average of your subset |
| 231 | + double *rec_indi_nums = NULL; |
| 232 | + rec_indi_nums = compute_qft(sub_rand_nums, num_elements_per_proc,world_rank); |
| 233 | + |
| 234 | + // Gather all partial averages down to the root process |
| 235 | + double *rec_gather_nums = NULL; |
| 236 | + if (world_rank == 0) { |
| 237 | + rec_gather_nums = (double *)malloc(sizeof(double) * num_elements_per_proc* world_size); |
| 238 | + assert(rec_gather_nums != NULL); |
| 239 | + } |
| 240 | + MPI_Gather(rec_indi_nums, num_elements_per_proc, MPI_DOUBLE, rec_gather_nums, num_elements_per_proc, MPI_DOUBLE, 0, MPI_COMM_WORLD); |
| 241 | + |
| 242 | + // Now that we have all of the partial averages on the root, compute the |
| 243 | + // total average of all numbers. Since we are assuming each process computed |
| 244 | + // an average across an equal amount of elements, this computation will |
| 245 | + // produce the correct answer. |
| 246 | + if (world_rank == 0) { |
| 247 | + //double avg = compute_avg(sub_avgs, world_size); |
| 248 | + // printf("Avg of all elements is %f\n", avg); |
| 249 | + // // Compute the average across the original data for comparison |
| 250 | + // double original_data_avg = |
| 251 | + // compute_avg(rand_nums, num_elements_per_proc * world_size); |
| 252 | + // printf("Avg computed across original data is %f\n", original_data_avg); |
| 253 | + // int i; |
| 254 | + // for (i = 0; i < num_elements_per_proc*world_size; i++) { |
| 255 | + // printf("%f kk", rec_gather_nums[i]); |
| 256 | + // printf( "\n" ); |
| 257 | + // } |
| 258 | + |
| 259 | + // ublas::vector<std::complex<double> > v1(4), v2(2); |
| 260 | + // for (double i = 0; i < 4; ++i) |
| 261 | + // v1 (i) = exp(I*i); |
| 262 | + // for (double i = 0; i < 2; ++i) |
| 263 | + // v2 (i) = i; |
| 264 | + |
| 265 | + // ublas::matrix<double> m(v1.size(), v2.size()); |
| 266 | + // std::cout |
| 267 | + // << v1 << '\n' |
| 268 | + // << v2 << '\n' |
| 269 | + // << outer_prod(v1, v2)<< '\n'; |
| 270 | + } |
| 271 | + |
| 272 | + |
| 273 | + |
| 274 | + |
| 275 | + |
| 276 | + |
| 277 | + |
| 278 | + |
| 279 | + |
| 280 | + // Clean up |
| 281 | + if (world_rank == 0) { |
| 282 | + free(rec_indi_nums); |
| 283 | + free(rec_gather_nums); |
| 284 | + } |
| 285 | + free(sub_rand_nums); |
| 286 | + |
| 287 | + MPI_Barrier(MPI_COMM_WORLD); |
| 288 | + MPI_Finalize(); |
| 289 | +} |
0 commit comments