examples/charm++: Always use CkWallTimer, not CmiWallTimer (tested)
[charm.git] / examples / charm++ / topology / matmul3d / matmul3d.C
1 /** \file matmul3d.C
2  *  Author: Abhinav S Bhatele
3  *  Date Created: March 13th, 2008
4  *
5  */
6
7 #include "matmul3d.decl.h"
8 #include "matmul3d.h"
9 #include "TopoManager.h"
10 #if CMK_BLUEGENEP || CMK_VERSION_BLUEGENE
11 #include "essl.h"
12 #endif
13
14 Main::Main(CkArgMsg* m) {
15   if ( (m->argc != 3) && (m->argc != 7) && (m->argc != 10) ) {
16     CkPrintf("%s [array_size] [block_size]\n", m->argv[0]);
17     CkPrintf("OR %s [array_size_X] [array_size_Y] [array_size_Z] [block_size_X] [block_size_Y] [block_size_Z]\n", m->argv[0]);
18     CkPrintf("OR %s [array_size_X] [array_size_Y] [array_size_Z] [block_size_X] [block_size_Y] [block_size_Z] [torus_dim_X] [torus_dim_Y] [torus_dim_Z] \n", m->argv[0]);
19     CkAbort("Abort");
20   }
21
22   // store the main proxy
23   mainProxy = thisProxy;
24
25   // get the size of the global array, size of each chare and size of the torus [optional]
26   if(m->argc == 3) {
27     arrayDimX = arrayDimY = arrayDimZ = atoi(m->argv[1]);
28     blockDimX = blockDimY = blockDimZ = atoi(m->argv[2]);
29   }
30   else if (m->argc == 7) {
31     arrayDimX = atoi(m->argv[1]);
32     arrayDimY = atoi(m->argv[2]);
33     arrayDimZ = atoi(m->argv[3]);
34     blockDimX = atoi(m->argv[4]);
35     blockDimY = atoi(m->argv[5]);
36     blockDimZ = atoi(m->argv[6]);
37   } else {
38     arrayDimX = atoi(m->argv[1]);
39     arrayDimY = atoi(m->argv[2]);
40     arrayDimZ = atoi(m->argv[3]);
41     blockDimX = atoi(m->argv[4]);
42     blockDimY = atoi(m->argv[5]);
43     blockDimZ = atoi(m->argv[6]);
44     torusDimX = atoi(m->argv[7]);
45     torusDimY = atoi(m->argv[8]);
46     torusDimZ = atoi(m->argv[9]);
47   }
48
49   if (arrayDimX < blockDimX || arrayDimX % blockDimX != 0)
50     CkAbort("array_size_X % block_size_X != 0!");
51   if (arrayDimY < blockDimY || arrayDimY % blockDimY != 0)
52     CkAbort("array_size_Y % block_size_Y != 0!");
53   if (arrayDimZ < blockDimZ || arrayDimZ % blockDimZ != 0)
54     CkAbort("array_size_Z % block_size_Z != 0!");
55
56   num_chare_x = arrayDimX / blockDimX;
57   num_chare_y = arrayDimY / blockDimY;
58   num_chare_z = arrayDimZ / blockDimZ;
59   num_chares = num_chare_x * num_chare_y * num_chare_z;
60
61   subBlockDimXz = blockDimX/num_chare_z;
62   subBlockDimYx = blockDimY/num_chare_x;
63   subBlockDimXy = blockDimX/num_chare_y;
64
65   // print info
66   CkPrintf("Running Matrix Multiplication on %d processors with (%d, %d, %d) chares\n", CkNumPes(), num_chare_x, num_chare_y, num_chare_z);
67   CkPrintf("Array Dimensions: %d %d %d\n", arrayDimX, arrayDimY, arrayDimZ);
68   CkPrintf("Block Dimensions: %d %d %d\n", blockDimX, blockDimY, blockDimZ);
69
70   // Create new array of worker chares
71 #if USE_TOPOMAP
72   CkPrintf("Topology Mapping is being done ...\n");
73   CProxy_ComputeMap map = CProxy_ComputeMap::ckNew(num_chare_x, num_chare_y, num_chare_z, 1, 1, 1);
74   CkArrayOptions opts(num_chare_x, num_chare_y, num_chare_z);
75   opts.setMap(map);
76   compute = CProxy_Compute::ckNew(opts);
77 #elif USE_BLOCKMAP
78   CkPrintf("Block Mapping is being done ...\n");
79   CProxy_ComputeMap map = CProxy_ComputeMap::ckNew(num_chare_x, num_chare_y, num_chare_z, torusDimX, torusDimY, torusDimZ);
80   CkArrayOptions opts(num_chare_x, num_chare_y, num_chare_z);
81   opts.setMap(map);
82   compute = CProxy_Compute::ckNew(opts);
83 #else
84   compute = CProxy_Compute::ckNew(num_chare_x, num_chare_y, num_chare_z);
85 #endif
86
87   // CkPrintf("Total Hops: %d\n", hops);
88
89   // Start the computation
90   numIterations = 0;
91   startTime = CkWallTimer();
92   compute.beginCopying();
93 }
94
95 void Main::done() {
96   numIterations++;
97   if(numIterations == 1) {
98     firstTime = CkWallTimer();
99     CkPrintf("FIRST ITER TIME %f secs\n", firstTime - startTime);
100   }
101
102   if(numIterations == NUM_ITER) {
103     endTime = CkWallTimer();
104     CkPrintf("AVG TIME %f secs\n", (endTime - firstTime)/(NUM_ITER-1));
105     CkExit();
106   } else {
107     compute.resetArrays();
108   }
109 }
110
111 // Constructor, initialize values
112 Compute::Compute() {
113   // each compute will only hold blockDimX/num_chare_z no. of rows to begin with
114   A = new float[blockDimX*blockDimY];
115   B = new float[blockDimY*blockDimZ];
116   C = new float[blockDimX*blockDimZ];
117
118   int indexX = thisIndex.x;
119   int indexY = thisIndex.y;
120   int indexZ = thisIndex.z;
121
122   float tmp;
123
124   for(int i=indexZ*subBlockDimXz; i<(indexZ+1)*subBlockDimXz; i++)
125     for(int j=0; j<blockDimY; j++) {
126       tmp = (float)drand48();
127       while(tmp > MAX_LIMIT || tmp < (-1)*MAX_LIMIT)
128         tmp = (float)drand48();
129
130       A[i*blockDimY + j] = tmp;
131   }
132
133   for(int j=indexX*subBlockDimYx; j<(indexX+1)*subBlockDimYx; j++)
134     for(int k=0; k<blockDimZ; k++) {
135       tmp = (float)drand48();
136       while(tmp > MAX_LIMIT || tmp < (-1)*MAX_LIMIT)
137         tmp = (float)drand48();
138
139       B[j*blockDimZ + k] = tmp;
140   }
141
142   for(int i=0; i<blockDimX; i++)
143     for(int k=0; k<blockDimZ; k++)
144       C[i*blockDimZ + k] = 0.0;
145
146   // counters to keep track of how many messages have been received
147   countA = 0;
148   countB = 0;
149   countC = 0;
150 }
151
152 Compute::Compute(CkMigrateMessage* m) {
153
154 }
155
156 Compute::~Compute() {
157   delete [] A;
158   delete [] B;
159   delete [] C;
160 }
161
162 void Compute::beginCopying() {
163   sendA();
164   sendB();
165 }
166
167 void Compute::resetArrays() {
168   int indexX = thisIndex.x;
169   int indexY = thisIndex.y;
170   int indexZ = thisIndex.z;
171
172   float tmp;
173
174   for(int i=indexZ*subBlockDimXz; i<(indexZ+1)*subBlockDimXz; i++)
175     for(int j=0; j<blockDimY; j++) {
176       tmp = (float)drand48();
177       while(tmp > MAX_LIMIT || tmp < (-1)*MAX_LIMIT)
178         tmp = (float)drand48();
179
180       A[i*blockDimY + j] = tmp;
181   }
182
183   for(int j=indexX*subBlockDimYx; j<(indexX+1)*subBlockDimYx; j++)
184     for(int k=0; k<blockDimZ; k++) {
185       tmp = (float)drand48();
186       while(tmp > MAX_LIMIT || tmp < (-1)*MAX_LIMIT)
187         tmp = (float)drand48();
188
189       B[j*blockDimZ + k] = tmp;
190   }
191
192   for(int i=0; i<blockDimX; i++)
193     for(int k=0; k<blockDimZ; k++) {
194       C[i*blockDimZ + k] = 0.0;
195 #if USE_CKDIRECT
196       tmpC[i*blockDimZ + k] = 0.0;
197 #endif
198     }
199
200   sendA();
201   sendB();
202 }
203
204 void Compute::sendA() {
205   int indexZ = thisIndex.z;
206
207   for(int k=0; k<num_chare_z; k++)
208     if(k != indexZ) {
209       // use a local pointer for chares on the same processor
210       Compute* c = compute(thisIndex.x, thisIndex.y, k).ckLocal();
211       if(c != NULL)
212         c->receiveA(indexZ, &A[indexZ*subBlockDimXz*blockDimY], subBlockDimXz * blockDimY);
213       else
214         compute(thisIndex.x, thisIndex.y, k).receiveA(indexZ, &A[indexZ*subBlockDimXz*blockDimY], subBlockDimXz * blockDimY);
215     }
216 }
217
218 void Compute::sendB() {
219   int indexX = thisIndex.x;
220
221   for(int i=0; i<num_chare_x; i++)
222     if(i != indexX) {
223       // use a local pointer for chares on the same processor
224       Compute* c = compute(i, thisIndex.y, thisIndex.z).ckLocal();
225       if(c != NULL)
226         c->receiveB(indexX, &B[indexX*subBlockDimYx*blockDimZ], subBlockDimYx * blockDimZ);
227       else
228         compute(i, thisIndex.y, thisIndex.z).receiveB(indexX, &B[indexX*subBlockDimYx*blockDimZ], subBlockDimYx * blockDimZ);
229     }
230 }
231
232 void Compute::sendC() {
233   int indexY = thisIndex.y;
234   for(int j=0; j<num_chare_y; j++) {
235     if(j != indexY) {
236       // use a local pointer for chares on the same processor
237       Compute *c = compute(thisIndex.x, j, thisIndex.z).ckLocal();
238       if(c != NULL)
239         c->receiveC(&C[j*subBlockDimXy*blockDimZ], subBlockDimXy * blockDimZ, 1);
240       else
241         compute(thisIndex.x, j, thisIndex.z).receiveC(&C[j*subBlockDimXy*blockDimZ], subBlockDimXy * blockDimZ, 1);
242     }
243   }
244 }
245
246 void Compute::receiveA(int indexZ, float *data, int size) {
247   for(int i=0; i<subBlockDimXz; i++)
248     for(int j=0; j<blockDimY; j++)
249       A[indexZ*subBlockDimXz*blockDimY + i*blockDimY + j] = data[i*blockDimY + j];
250   countA++;
251   if(countA == num_chare_z-1)
252     doWork();
253 }
254
255 void Compute::receiveB(int indexX, float *data, int size) {
256   for(int j=0; j<subBlockDimYx; j++)
257     for(int k=0; k<blockDimZ; k++)
258       B[indexX*subBlockDimYx*blockDimZ + j*blockDimZ + k] = data[j*blockDimZ + k];
259   countB++;
260   if(countB == num_chare_x-1)
261     doWork();
262 }
263
264 void Compute::receiveC(float *data, int size, int who) {
265   int indexY = thisIndex.y;
266   if(who) {
267     for(int i=0; i<subBlockDimXy; i++)
268       for(int k=0; k<blockDimZ; k++)
269         C[indexY*subBlockDimXy*blockDimZ + i*blockDimZ + k] += data[i*blockDimZ + k];
270   }
271   countC++;
272   if(countC == num_chare_y) {
273     /*char name[30];
274     sprintf(name, "%s_%d_%d_%d", "C", thisIndex.x, thisIndex.y, thisIndex.z);
275     FILE *fp = fopen(name, "w");
276     for(int i=0; i<subBlockDimXy; i++) {
277       for(int k=0; k<blockDimZ; k++)
278         fprintf(fp, "%f ", C[indexY*subBlockDimXy*blockDimZ + i*blockDimZ + k]);
279       fprintf(fp, "\n");
280     }
281     fclose(fp);*/
282
283     // counters to keep track of how many messages have been received
284     countA = 0;
285     countB = 0;
286     countC = 0;
287
288     contribute(0, 0, CkReduction::concat, CkCallback(CkIndex_Main::done(), mainProxy));
289     // mainProxy.done();
290   }
291 }
292
293 void Compute::doWork() {
294   if(countA == num_chare_z-1 && countB == num_chare_x-1) {
295
296 #if CMK_BLUEGENEP || CMK_VERSION_BLUEGENE
297     const char trans = 'N';
298     const double alpha = 1.0;
299     const double beta = 0.0;
300
301     sgemm(&trans, &trans, blockDimX, blockDimZ, blockDimY, alpha, A, blockDimX, B, blockDimY, beta, C, blockDimX);
302 #else
303     for(int i=0; i<blockDimX; i++)
304       for(int j=0; j<blockDimY; j++)
305         for(int k=0; k<blockDimZ; k++)
306           C[i*blockDimZ+k] += A[i*blockDimY+j] * B[j*blockDimZ+k];
307 #endif
308     
309     receiveC(&C[(thisIndex.y)*subBlockDimXy*blockDimZ], subBlockDimXy*blockDimZ, 0);
310     sendC();
311   }
312 }
313
314 ComputeMap::ComputeMap(int x, int y, int z, int tx, int ty, int tz) {
315   X = x;
316   Y = y;
317   Z = z;
318   mapping = new int[X*Y*Z];
319
320   TopoManager tmgr;
321   int dimX, dimY, dimZ, dimT;
322
323 #if USE_TOPOMAP
324   dimX = tmgr.getDimNX();
325   dimY = tmgr.getDimNY();
326   dimZ = tmgr.getDimNZ();
327   dimT = tmgr.getDimNT();
328 #elif USE_BLOCKMAP
329   dimX = tx;
330   dimY = ty;
331   dimZ = tz;
332   dimT = 1;
333 #endif
334
335   // we are assuming that the no. of chares in each dimension is a 
336   // multiple of the torus dimension
337   int numCharesPerPe = X*Y*Z/CkNumPes();
338
339   int numCharesPerPeX = X / dimX;
340   int numCharesPerPeY = Y / dimY;
341   int numCharesPerPeZ = Z / dimZ;
342
343   if(dimT < 2) {    // one core per node
344     if(CkMyPe()==0) CkPrintf("DATA: %d %d %d %d : %d %d %d\n", dimX, dimY, dimZ, dimT, numCharesPerPeX, numCharesPerPeY, numCharesPerPeZ);
345     for(int i=0; i<dimX; i++)
346       for(int j=0; j<dimY; j++)
347         for(int k=0; k<dimZ; k++)
348           for(int ci=i*numCharesPerPeX; ci<(i+1)*numCharesPerPeX; ci++)
349             for(int cj=j*numCharesPerPeY; cj<(j+1)*numCharesPerPeY; cj++)
350               for(int ck=k*numCharesPerPeZ; ck<(k+1)*numCharesPerPeZ; ck++) {
351 #if USE_TOPOMAP
352                 mapping[ci*Y*Z + cj*Z + ck] = tmgr.coordinatesToRank(i, j, k);
353 #elif USE_BLOCKMAP
354                 mapping[ci*Y*Z + cj*Z + ck] = i + j*dimX + k*dimX*dimY;
355 #endif
356               }
357   } else {          // multiple cores per node
358     // In this case, we split the chares in the X dimension among the
359     // cores on the same node.
360     numCharesPerPeX /= dimT;
361       if(CkMyPe()==0) CkPrintf("%d %d %d : %d %d %d %d : %d %d %d \n", x, y, z, dimX, dimY, dimZ, dimT, numCharesPerPeX, numCharesPerPeY, numCharesPerPeZ);
362       for(int i=0; i<dimX; i++)
363         for(int j=0; j<dimY; j++)
364           for(int k=0; k<dimZ; k++)
365             for(int l=0; l<dimT; l++)
366               for(int ci=(dimT*i+l)*numCharesPerPeX; ci<(dimT*i+l+1)*numCharesPerPeX; ci++)
367                 for(int cj=j*numCharesPerPeY; cj<(j+1)*numCharesPerPeY; cj++)
368                   for(int ck=k*numCharesPerPeZ; ck<(k+1)*numCharesPerPeZ; ck++) {
369                     mapping[ci*Y*Z + cj*Z + ck] = tmgr.coordinatesToRank(i, j, k, l);
370                   }
371   }
372 }
373
374 ComputeMap::~ComputeMap() {
375   delete [] mapping;
376 }
377
378 int ComputeMap::procNum(int, const CkArrayIndex &idx) {
379   int *index = (int *)idx.data();
380   return mapping[index[0]*Y*Z + index[1]*Z + index[2]];
381 }
382
383