Changes to the advancedlb text and adding the location of an example
[charm.git] / examples / charm++ / topology / matmul2d / matmul2d.C
1 #include "matmul2d.decl.h"
2 #include "matmul2d.h"
3 #include "TopoManager.h"
4 /*
5 #include <cstdio>
6 #include <iostream>
7 #include <fstream>
8 #include <sstream>
9 */
10 #include "stdio.h"
11 #include "stdlib.h"
12
13 Main::Main(CkArgMsg* m) {
14   if (m->argc != 4){ 
15     CkPrintf("%s [N] [K] [num_chares_per_dim]\n", m->argv[0]);
16     CkAbort("Abort");
17   }
18   else {
19     N = atoi(m->argv[1]);
20     K = atoi(m->argv[2]);
21     num_chares_per_dim = atoi(m->argv[3]);
22     T = N/num_chares_per_dim;
23   }
24
25   // store the main proxy
26   mainProxy = thisProxy;
27
28   // get the size of the global array, size of each chare and size of the torus [optional]
29   /*
30   else if (m->argc == 5) {
31     arrayDimX = atoi(m->argv[1]);
32     arrayDimY = atoi(m->argv[2]);
33     blockDimX = atoi(m->argv[3]);
34     blockDimY = atoi(m->argv[4]);
35   } else if (m->argc == 7){
36     arrayDimX = atoi(m->argv[1]);
37     arrayDimY = atoi(m->argv[2]);
38     blockDimX = atoi(m->argv[3]);
39     blockDimY = atoi(m->argv[4]);
40     torusDimX = atoi(m->argv[5]);
41     torusDimY = atoi(m->argv[6]);
42   }
43   */
44
45   if (N < T || N % T != 0)
46     CkAbort("N % T != 0!");
47   if (K < T || K % T != 0)
48     CkAbort("K % T != 0!");
49
50   // print info
51   CkPrintf("Running Matrix Multiplication on %d processors with (%d, %d) chares\n", CkNumPes(), num_chares_per_dim, num_chares_per_dim);
52   CkPrintf("Array Dimensions: %d %d\n", N, K);
53   CkPrintf("Block Dimensions: %dx%d, %dx%d\n", T, K/num_chares_per_dim, K/num_chares_per_dim, T);
54   CkPrintf("Chare-array Dimensions: %d %d\n", num_chares_per_dim, num_chares_per_dim);
55
56   // Create new array of worker chares
57   compute = CProxy_Compute::ckNew(num_chares_per_dim, num_chares_per_dim);
58
59   // Start the computation
60   startTime = CmiWallTimer();
61   if(num_chares_per_dim == 1){
62     compute(0,0).compute();
63   }
64   else{
65     compute.compute();
66     //compute.start();
67   }
68 }
69
70 void Main::done(){
71   endTime = CmiWallTimer();
72   CkPrintf("Fin: %f sec\n", endTime-startTime);
73   CkExit();
74 }
75
76 // Constructor, initialize values
77 Compute::Compute() {
78
79   int s1 = (K/num_chares_per_dim)*T;
80   int s2 = T*T;
81   for(int i = 0; i < 2; i++){
82     A[i] = new float[s1];
83     B[i] = new float[s1];
84   }
85   C = new float[s2];
86
87   for(int i = 0; i < s1; i++){
88     A[0][i] = MAGIC_A;
89     B[0][i] = MAGIC_B;
90   }
91   for(int i = 0; i < s2; i++){
92     C[i] = 0;
93   }
94     
95   step = 0;  
96   row = thisIndex.y;
97   col = thisIndex.x;
98
99   whichLocal = 0;
100   remaining = 2;
101   iteration = 0;
102
103   //comps = 0;
104 }
105
106 Compute::Compute(CkMigrateMessage* m) {
107 }
108
109 Compute::~Compute() {
110   delete [] A[0];
111   delete [] A[1];
112   delete [] B[0];
113   delete [] B[1];
114   delete [] C;
115 }
116
117 void Compute::start(){
118   int newBuf = 1 - whichLocal;
119   /*
120   //1. send A
121   Compute *c = thisProxy((col-row+num_chare_x)%num_chare_x, row).ckLocal();
122   if(c == 0){
123     thisProxy((col-row+num_chare_x)%num_chare_x, row).recvBlockA(A[whichLocal], blockDimX*blockDimY, newBuf);
124   }
125   else{
126     c->recvBlockA(A[whichLocal], blockDimX*blockDimY, newBuf);
127   }
128   
129   //2. send B
130   c = thisProxy(col, (row-col+num_chare_y)%num_chare_y).ckLocal();
131   if(c == 0){
132     thisProxy(col, (row-col+num_chare_y)%num_chare_y).recvBlockB(B[whichLocal], blockDimX*blockDimY, newBuf);
133   }
134   else{
135     c->recvBlockB(B[whichLocal], blockDimX*blockDimY, newBuf);
136   }
137   */
138   whichLocal = newBuf;
139 }
140
141 void Compute::compute(){
142   //int count = 0;
143 #ifdef USE_OPT_ROUTINES
144   const char trans = 'N';
145   const double alpha = 1.0;
146   const double beta = 0.0;
147
148   sgemm(&trans, &trans, blockDimX, blockDimZ, blockDimY, alpha, A, blockDimX, B, blockDimY, beta, C, blockDimX);
149 #else
150   int i, j, k;
151
152   float *thisa = A[whichLocal];
153   float *thisb = B[whichLocal];
154
155   for(i = 0; i < T; i++){
156     for(j = 0; j < T; j++){
157       float sum = 0.0;
158       for(k = 0; k < (K/num_chares_per_dim); k++){
159         sum += thisa[i*(K/num_chares_per_dim)+k]*thisb[k*(T)+j];
160         //comps++;
161       }
162       C[i*T+j] += sum;
163     }
164   }
165 #endif
166
167   remaining = 2;
168   iteration++;
169   if(iteration == num_chares_per_dim){
170 #ifdef MATMUL2D_WRITE_FILE
171     // create a file
172     //ostringstream oss;
173     char buf[128];
174     
175     sprintf(buf, "mat.%d.%d", row, col);
176     //oss << "mat." << row << "." << col;
177     //ofstream ofs(oss.str().c_str());
178     FILE *fp = fopen(buf, "w");
179
180     for(int i = 0; i < blockDimX; i++){
181       for(int j = 0; j < blockDimX; j++){
182         //ofs << C[i*blockDimX+j];
183         fprintf(fp, "%f ", C[i*blockDimX+j]);
184       }
185       //ofs << endl;
186       fprintf(fp, "\n");
187
188     }
189     fclose(fp);
190 #endif
191     contribute(0,0,CkReduction::concat, CkCallback(CkIndex_Main::done(), mainProxy));
192     CkPrintf("[%d,%d] comps: %d iter: %d\n", thisIndex.x, thisIndex.y, -1, iteration);
193   }
194   else{
195     contribute(0,0,CkReduction::concat,CkCallback(CkIndex_Compute::resumeFromBarrier(), thisProxy));
196   }
197 }
198
199 void Compute::resumeFromBarrier(){
200   // At this point, everyone has used their A and B buffers
201   int newBuf = 1-whichLocal;
202   // We must put our own 
203   // 1. First put A
204
205   //if(num_chare_x == 0 || num_chare_y ==0)
206   //  CkPrintf("(%d,%d): 0 divisor\n", thisIndex.y, thisIndex.x);
207 #ifdef MATMUL2D_DEBUG
208   CkPrintf("(%d,%d): A nbr: (%d,%d), iteration: %d\n", thisIndex.y, thisIndex.x, row, (col-1+num_chare_x)%num_chare_x, iteration);
209 #endif
210   /*
211   CkPrintf("(%d,%d): B nbr: (%d,%d)\n", thisIndex.y, thisIndex.x, (row-1+num_chare_y)%num_chare_y, col);
212   */
213
214   int size = (K/num_chares_per_dim)*T;
215
216   //Compute *c = thisProxy((col-1+num_chares_per_dim)%num_chares_per_dim, row).ckLocal();
217   //if(c == 0){
218   thisProxy((col-1+num_chares_per_dim)%num_chares_per_dim, row).recvBlockA(A[whichLocal], size, newBuf);
219   //}
220   //else{
221   //  c->recvBlockA(A[whichLocal], blockDimX*blockDimY, newBuf);
222   //}
223   // 2. Then put B
224   //c =  thisProxy(col, (row-1+num_chare_y)%num_chare_y).ckLocal();
225   //if(c == 0){
226   thisProxy(col, (row-1+num_chares_per_dim)%num_chares_per_dim).recvBlockB(B[whichLocal], size, newBuf);
227   //}
228   //else{
229   //  c->recvBlockB(B[whichLocal], blockDimX*blockDimY, newBuf);
230   //}
231   // toggle between local buffers
232   whichLocal = newBuf;
233 }
234
235 void Compute::recvBlockA(float *block, int size, int whichBuf){
236   memcpy(A[whichBuf], block, sizeof(float)*size);
237   remaining--;
238   if(remaining == 0){
239     compute();
240   }
241 }
242
243 void Compute::recvBlockB(float *block, int size, int whichBuf){
244   memcpy(B[whichBuf], block, sizeof(float)*size);
245   remaining--;
246   if(remaining == 0){
247     compute();
248   }
249 }
250
251 ComputeMap::ComputeMap(int _adx, int _ady){
252   arrayDimX = _adx;
253   arrayDimY = _ady;
254
255   map = new int[_adx*_ady];
256   
257 }
258
259 int ComputeMap::procNum(int arrayHdl, const CkArrayIndex &idx){
260   int *index = (int *)(idx.data());
261   int row = index[1];
262   int col = index[0];
263
264   return map[row*arrayDimX + col];
265 }
266 #include "matmul2d.def.h"