고른 표본 추출을 통한 빠른 정렬

고른 표본 추출을 통한 병렬 정렬(Parallel Sorting by Regular Sampling: PSRS)은 n개의 프로세스에서 골고루 표본을 추출하여 이것을 주축으로 하여 값들을 정렬하는 병렬 정렬 알고리즘의 일종입니다. 균형이 잘 잡히고, 프로세스의 수가 2의 n제곱 꼴이 되지 않아도 된다는 등의 장점이 있습니다.

MPI를 배우면서 이것을 한번 C+MPI로 구현해 보았습니다. p개의 노드에서 n/p개만큼의 [0..1) 범위의 float형 난수를 발생한 다음에 정렬을 하는 과정입니다. 잘 작성된 코드는 아니지만, 참고할 수 있는 코드일 수도 있어서 실어 봅니다.

  1 /**
  2  * PSRS implementation using MPI.
  3  *
  4  * Date: 5th of December, 2006
  5  * Author: Jay
  6  *
  7  * Compile options:
  8  *     PRINT_MSG: print some messages to stdout.
  9  *     OUTPUT_FILE: write 'glist_nn.txt', and 'slist_nn.txt'.
 10  *                  'glist_nn' contains sorted generated numbers
 11  *                  of each node. 'slist_nn' is final result
 12  *                  of each node.
 13  */
 14 #include <stdlib.h>
 15 #include <stdio.h>
 16 #include <mpi.h>
 17 #include <limits.h>
 18 #include <time.h>
 19 #include <stddef.h>
 20 #include <string.h>
 21 
 22 /* Upper bound of generated floating point */
 23 #define UPPER_BOUND 1.0
 24 
 25 /* float comparision function */
 26 int float_comp(const void* aa, const void* bb)
 27 {
 28     const float* a = aa;
 29     const float* b = bb;
 30     if (*a == *b) return 0;
 31     if (*a < *b) return -1;
 32     return 1;
 33 }
 34 
 35 /* Generates random sequence into array A with size */
 36 void generate_random_sequence(float* A, size_t size)
 37 {
 38     size_t i;
 39     for (i=0; i<size; i++)
 40         A[i] = (float)rand()/RAND_MAX;
 41 }
 42 
 43 /* Returns pointer to the lower_bound,
 44  * which means the lowest position
 45  * that the element can be inserted
 46  * with sorted order. */
 47 float* lower_bound(float* first, float* last, float val)
 48 {
 49     ptrdiff_t len = last - first;
 50     while (len > 0)
 51     {
 52         ptrdiff_t half = len / 2;
 53         float* middle = first + half;
 54         if (*middle < val)
 55         {
 56             first = middle + 1;
 57             len = len - half - 1;
 58         }
 59         else
 60             len = half;
 61     }
 62     return first;
 63 }
 64 
 65 /* Writes each elements in float array to the file */
 66 void output(const char* filename, float* A, size_t size)
 67 {
 68     FILE* f = fopen(filename, "w");
 69     int i;
 70     for (i=0; i<size; i++)
 71     {
 72         fprintf(f, "%1.5f\n", A[i]);
 73     }
 74     fclose(f);
 75 }
 76 
 77 /*
 78  * First command-line argument: n
 79  */
 80 int main(int argc, char* argv[])
 81 {
 82     int n_number = 1000000;
 83     if (argc > 1) n_number = atoi(argv[1]);
 84     int id, p;
 85 
 86     MPI_Init(&argc, &argv);
 87 
 88     MPI_Barrier(MPI_COMM_WORLD);
 89     double elapsed_time = -MPI_Wtime();
 90 
 91     MPI_Comm_rank(MPI_COMM_WORLD, &id);
 92     MPI_Comm_size(MPI_COMM_WORLD, &p);
 93 
 94     /* Some nodes should control 1 more
 95      * element if n % p != 0. */
 96     float A[(n_number+(p-1))/p];
 97     float* single_list = A;
 98     int n_control = n_number/p;
 99     int n_sorted = n_control;
100     if (n_number%p > id) n_control++;
101     int i;
102 
103 #ifdef PRINT_MSG    
104     fprintf(stdout, "Process %d generates %d of %d elements.\n",
105             id, n_control, n_number);
106     fflush(stdout);
107 #endif
108 
109     /* For unique random number, add (id*1000) */
110     srand( time(NULL) + id * 1000);
111     generate_random_sequence(A, n_control);
112 
113     /* Phase 1 */
114     /* Each process quicksort their own list and each one picks samples */
115     qsort(A, n_control, sizeof(float), float_comp);
116 
117     if (p > 1)
118     {
119 
120         float samples[p];
121         for (i=0; i<p; i++)
122             samples[i] = A[i*n_control/p];
123 
124         /* Phase 2 */
125         /* Node 0 gathers samples, sorts them, and picks pivots. */
126         float all_samples[p*p];
127         MPI_Gather(samples, p, MPI_FLOAT, all_samples, p,
128                 MPI_FLOAT, 0, MPI_COMM_WORLD);
129         float pivots[p-1];
130         if (!id)
131         {
132             qsort(all_samples, p*p, sizeof(float), float_comp);
133 
134             for (i=0; i<p-1; i++)
135                 pivots[i] = all_samples[(i+1)*p+p/2-1];
136         }
137         /* Node 0 broadcasts pivots and each process 
138          * partitions its own list. */
139         MPI_Bcast(pivots, p-1, MPI_FLOAT, 0, MPI_COMM_WORLD);
140         int send_cnts[p], send_disp[p];
141         send_disp[0] = 0;
142         for (i=1; i<p; i++)
143         {
144             send_disp[i] =
145                 (float*)(lower_bound(A, A+n_control, pivots[i-1]))-A;
146             send_cnts[i-1] = send_disp[i] - send_disp[i-1];
147         }
148         send_cnts[p-1] = n_control - send_disp[p-1];
149 
150         /* Phase 3 */
151         /* First, exchanges the number of elements that 
152          * each one is going to exchange. */
153         int recv_cnts[p], recv_disp[p+1];
154         MPI_Alltoall(send_cnts, 1, MPI_FLOAT, recv_cnts, 1,
155                 MPI_FLOAT, MPI_COMM_WORLD);
156         recv_disp[0] = 0;
157         for (i=1; i<p; i++)
158             recv_disp[i] = recv_disp[i-1] + recv_cnts[i-1];
159         recv_disp[p] = recv_disp[p-1]+recv_cnts[p-1];
160         float partitions[recv_disp[p]];
161         /* Exchanges elements to appropriate nodes. */
162         MPI_Alltoallv(A, send_cnts, send_disp, MPI_FLOAT, partitions,
163                 recv_cnts, recv_disp, MPI_FLOAT, MPI_COMM_WORLD);
164 
165         /* Phase 4 */
166         /* Each node merges its own partitions into a single list. */
167         int j;
168         int merge_disp[p];
169         n_sorted = recv_disp[p];
170         single_list = malloc(n_sorted*sizeof(float));
171         memcpy(merge_disp, recv_disp, p*sizeof(int));
172         for (i=0; i<n_sorted; i++)
173         {
174             float min = UPPER_BOUND;
175             int min_pos = 0;
176             for (j=0; j<p; j++)
177                 if (merge_disp[j] < recv_disp[j+1]
178                         && min > partitions[merge_disp[j]])
179                 {
180                     min = partitions[merge_disp[j]];
181                     min_pos = j;
182                 }
183             single_list[i] = min;
184             merge_disp[min_pos]++;
185         }
186 
187     }
188     /* Synchronizes for checking maximum elapsed time among nodes. */
189     MPI_Barrier(MPI_COMM_WORLD);
190     elapsed_time += MPI_Wtime();
191 
192 #ifdef PRINT_MSG    
193     fprintf(stdout, "Process %d now has sorted the list that contains \
194 %d of %d elements.\n", id, n_sorted, n_number);
195     fflush(stdout);
196 #endif
197 
198     if (!id)
199         printf("Elapsed Time with %d processes: %10.6f\n",
200                 p, elapsed_time);
201 
202     /* Output (elapsed_time doesn't count for file output!) */
203 #ifdef OUTPUT_FILE      
204     char filename[100];
205     strcpy(filename, "glistxx.txt");
206     filename[5] = '0' + id / 10;
207     filename[6] = '0' + id % 10;
208     output(filename, A, n_control);
209     filename[0] = 's';
210     output(filename, single_list, n_sorted);
211 #endif
212     if (single_list != A) free(single_list);
213     MPI_Finalize();
214 
215     return 0;
216 }
217