ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/BioCocoa/BioSwarm/trunk/Simulation/BioSwarmMPIController.m
Revision: 350
Committed: Wed Aug 19 16:45:07 2015 UTC (4 years, 9 months ago) by schristley
File size: 54704 byte(s)
Log Message:
MPI timing
Line File contents
1 /*
2 BioSwarmMPIController.m
3
4 The top level controller for an MPI simulation run.
5
6 Copyright (C) 2014 Scott Christley
7 Author: Scott Christley <schristley@mac.com>
8 Date: December 2014
9
10 This file is part of the BioSwarm Framework.
11
12 This library is free software; you can redistribute it and/or
13 modify it under the terms of the GNU Lesser General Public
14 License as published by the Free Software Foundation; either
15 version 2 of the License, or (at your option) any later version.
16
17 This library is distributed in the hope that it will be useful,
18 but WITHOUT ANY WARRANTY; without even the implied warranty of
19 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 Library General Public License for more details.
21
22 You should have received a copy of the GNU Lesser General Public
23 License along with this library; if not, write to the Free
24 Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
25 Boston, MA 02111 USA.
26 */
27
28 #import "BioSwarmMPIController.h"
29 #import "SubcellularElementModel.h"
30
31 #include <mpi.h>
32
33 #import "internal.h"
34
35 @implementation BioSwarmMPIController
36
37 - initWithModel: (NSDictionary *)aModel andManager: anObj
38 {
39 id param;
40 int i, j;
41 BOOL good = YES;
42
43 [super initWithModel: aModel andManager: anObj];
44
45 // assume models are in subspace
46 [spatialModel setSubspace: YES];
47 [spatialModel setSubspaceTranslator: self];
48
49 // look for MPI settings
50 NSDictionary *mpi = [aModel objectForKey: @"MPI"];
51 CHECK_PARAMF(mpi, "MPI", good);
52 if (!good) return nil;
53
54 [model setObject: @"YES" forKey: @"isMPI"];
55
56 NSArray *a = [mpi objectForKey: @"topologySize"];
57 CHECK_PARAM(a, "topologySize");
58 if ([a count] < 1) {
59 printf("ERROR: Size of topologySize is not valid: %d\n", [a count]);
60 return nil;
61 }
62
63 int dims = [spatialModel dimensions];
64 if ([a count] != dims) {
65 printf("ERROR: Dimensions of toplogySize does not match spatial dimensions, %ld != %d\n", (unsigned long)[a count], dims);
66 exit(1);
67 }
68
69 topologySize = malloc(dims * sizeof(int));
70 periodic = malloc(dims * sizeof(int));
71 coordinates = malloc(dims * sizeof(int));
72 for (i = 0; i < dims; ++i) {
73 topologySize[i] = [[a objectAtIndex: i] intValue];
74 periodic[i] = 0;
75 }
76
77 // initialize MPI
78 int res = MPI_Init(NULL, NULL);
79
80 // Create topology, allow reorder of ranks
81 MPI_Cart_create(MPI_COMM_WORLD, dims, topologySize, periodic, 1, &topologyMPIComm);
82
83 // get MPI rank
84 res = MPI_Comm_rank(MPI_COMM_WORLD, &rank);
85 printf("my rank = %d, res = %d\n", rank, res);
86
87 // coordinates
88 MPI_Cart_coords(topologyMPIComm, rank, dims, coordinates);
89 if (dims == 2) printf("rank = %d, coordinates = %d %d\n", rank, coordinates[0], coordinates[1]);
90 else printf("rank = %d, coordinates = %d %d %d\n", rank, coordinates[0], coordinates[1], coordinates[2]);
91
92 int src, dest;
93 MPI_Cart_shift(topologyMPIComm, 0, 1, &src, &dest);
94 printf("rank = %d, down, src = %d, dest = %d\n", rank, src, dest);
95 MPI_Cart_shift(topologyMPIComm, 0, -1, &src, &dest);
96 printf("rank = %d, up, src = %d, dest = %d\n", rank, src, dest);
97 MPI_Cart_shift(topologyMPIComm, 1, 1, &src, &dest);
98 printf("rank = %d, right, src = %d, dest = %d\n", rank, src, dest);
99 MPI_Cart_shift(topologyMPIComm, 1, -1, &src, &dest);
100 printf("rank = %d, left, src = %d, dest = %d\n", rank, src, dest);
101
102 // cache ranks for neighbor MPI nodes
103 for (i = S_CURRENT; i < S_XYZ_TOT; ++i) neighborRanks[i] = -1;
104 neighborRanks[S_CURRENT] = rank;
105
106 int maxNeighbors;
107 if (dims == 2) maxNeighbors = S_XY_TOT;
108 else maxNeighbors = S_XYZ_TOT;
109 for (i = S_START; i < maxNeighbors; ++i) {
110 int c[3];
111 int nRank = MPI_PROC_NULL;
112 bc_mpi_sector_coordinate_shift(coordinates[0], coordinates[1], coordinates[2], i, &(c[0]), &(c[1]), &(c[2]));
113 if (c[0] < 0) { neighborRanks[i] = MPI_PROC_NULL; continue; }
114 if (c[1] < 0) { neighborRanks[i] = MPI_PROC_NULL; continue; }
115 if (c[2] < 0) { neighborRanks[i] = MPI_PROC_NULL; continue; }
116 if (c[0] >= topologySize[0]) { neighborRanks[i] = MPI_PROC_NULL; continue; }
117 if (c[1] >= topologySize[1]) { neighborRanks[i] = MPI_PROC_NULL; continue; }
118 if (c[2] >= topologySize[2]) { neighborRanks[i] = MPI_PROC_NULL; continue; }
119
120 MPI_Cart_rank(topologyMPIComm, c, &nRank);
121 neighborRanks[i] = nRank;
122 printf("rank = %d (%d, %d, %d), %s neighbor = %d (%d, %d, %d)\n", rank, coordinates[0], coordinates[1], coordinates[2],
123 bc_sector_direction_string(i), nRank, c[0], c[1], c[2]);
124 }
125
126 if (TIME_PROFILE) bc_record_time([[NSString stringWithFormat: @"initialize MPI rank = %d", rank] UTF8String]);
127
128 return self;
129 }
130
131 - (void)dealloc
132 {
133 if (topologySize) free(topologySize);
134
135 [super dealloc];
136 }
137
138 - (int)rank { return rank; }
139 - (int *)neighborRanks { return neighborRanks; }
140 - (MPI_Comm)topologyMPICommunicator { return topologyMPIComm; }
141
142 - (NSString *)outputSuffix { return [NSString stringWithFormat: @"rank%d_%d", rank, outputRun]; }
143
144 - (BOOL)saveModelSpecification:(NSDictionary *)aModel forRun:(int)aNum
145 {
146 // only one MPI process needs to save
147 if (rank == 0) return [super saveModelSpecification: aModel forRun: aNum];
148 return YES;
149 }
150
151 //
152 // Override core BioSwarmController methods to add MPI
153 //
154
155 - (void)allocateHierarchy: (int)modelIndex withEncode: (NSString *)anEncode
156 {
157 // default
158 [super allocateHierarchy: modelIndex withEncode: anEncode];
159
160 // MPI
161 id modelObject = [modelInstances objectAtIndex: modelIndex];
162 [modelObject allocateMPIHierarchyWithEncode: anEncode];
163 }
164
165 - (void)initializeHierarchy
166 {
167 printf("initializeHierarchy\n");
168 if (TIME_PROFILE) bc_record_time([[NSString stringWithFormat: @"start MPI init simulation, rank = %d", rank] UTF8String]);
169 [super initializeHierarchy];
170 if (TIME_PROFILE) bc_record_time([[NSString stringWithFormat: @"end MPI init simulation, rank = %d", rank] UTF8String]);
171
172 // if initialized from file, allow specific MPI initialization
173 if (initializeFile) {
174 printf("initializeHierarchy saveIF\n");
175 NSInvocation *anInv = [NSInvocation invocationWithMethodSignature: [BioSwarmModel instanceMethodSignatureForSelector: @selector(initializeMPIFileCompleted)]];
176 [anInv setSelector: @selector(initializeMPIFileCompleted)];
177 [self performHierarchyInvocation: anInv];
178 }
179 }
180
181 - (void)initializeHierarchyForModel:(int)modelIndex
182 {
183 // default
184 [super initializeHierarchyForModel: modelIndex];
185
186 // MPI
187 id modelObject = [modelInstances objectAtIndex: modelIndex];
188 [modelObject initializeMPIHierarchy];
189 }
190
191 //
192 // Translate coordinates from/to world and subspace
193 //
194 // MPI uses a flipped coordinate system where (y, x) == (0, 0) is in the
195 // upper left corner for a 2D system. (y, x, z) == (0, 0, 0) is in the
196 // upper left front corner for a 3D system.
197 //
198 // BioSwarm has coordinate system where (x, y, z) == (0, 0, 0) is in the
199 // lower left corner.
200 //
201 // Translating back/forth requires using both the flipped coordinate system
202 // and the proper coordinate index.
203 //
204
205 // tweak float precision
206 // if precision is lost in addition, out != in + a
207 // adjust by scaled epsilon to maintain correct inequality
208 static float tweak_precision(float in, double a)
209 {
210 float out = in + a;
211 int i;
212
213 // out should be less than a
214 if ((in < 0) && !(out < a)) {
215 int scale = determine_scale_float(out);
216 float tweak = 0.0000001;
217 for (i = 0; i < scale; ++i) tweak *= 10;
218 printf("WARNING: tweak float precision, %e + %lf !< %f, tweak %f - %f = %f \n", in, a, out, out, tweak, out - tweak);
219 out = out - tweak;
220 return out;
221 }
222
223 // out should be greater than a
224 if ((in > 0) && !(out > a)) {
225 int scale = determine_scale_float(out);
226 float tweak = 0.0000001;
227 for (i = 0; i < scale; ++i) tweak *= 10;
228 printf("WARNING: tweak float precision, %e + %lf !> %f, tweak %f + %f = %f \n", in, a, out, out, tweak, out + tweak);
229 out = out + tweak;
230 return out;
231 }
232
233 return out;
234 }
235
236 static float tweak_precision_direction(float in, double a, BOOL roundUp)
237 {
238 float out = in + a;
239 float temp = out - a;
240 int i;
241
242 if (roundUp) {
243 if (temp < in) { // error rounded down, we want up
244 int scale = determine_scale_float(out);
245 float tweak = 0.0000001;
246 for (i = 0; i < scale; ++i) tweak *= 10;
247 //printf("WARNING: tweak float precision (roundUp = %d), %e + %e !> %e, tweak %e + %e = %e \n", roundUp, in, a, out, out, tweak, out + tweak);
248 out = out + tweak;
249 return out;
250 }
251 } else {
252 if (temp > in) { // error rounded up, we want down
253 int scale = determine_scale_float(out);
254 float tweak = 0.0000001;
255 for (i = 0; i < scale; ++i) tweak *= 10;
256 //printf("WARNING: tweak float precision (roundUp = %d), %e + %e !< %e, tweak %e - %e = %e \n", roundUp, in, a, out, out, tweak, out - tweak);
257 out = out - tweak;
258 return out;
259 }
260 }
261
262 return out;
263 }
264
265 - (void)translateIntegerCoordinatesToWorld: (int *)localCoordinates forModel: aModel
266 {
267 int *gridSize = [spatialModel gridSize];
268 int i, dims = [spatialModel dimensions];
269
270 localCoordinates[0] = coordinates[1] * gridSize[0] + localCoordinates[0];
271 localCoordinates[1] = (topologySize[1] - coordinates[0] - 1) * gridSize[0] + localCoordinates[1];
272 if (dims == 3) localCoordinates[2] = coordinates[2] * gridSize[2] + localCoordinates[2];
273 }
274
275 - (void)translateFloatCoordinatesToWorld: (float *)localCoordinates forModel: aModel
276 {
277 double *domainSize = [spatialModel domainSize];
278 int i, dims = [spatialModel dimensions];
279
280 // single precision can be problematic so tweak with EPS if necessary
281 localCoordinates[0] = tweak_precision(localCoordinates[0], coordinates[1] * domainSize[0]);
282 localCoordinates[1] = tweak_precision(localCoordinates[1], (topologySize[1] - coordinates[0] - 1) * domainSize[0]);
283 if (dims == 3) localCoordinates[2] = tweak_precision(localCoordinates[2], coordinates[2] * domainSize[2]);
284 }
285
286 - (void)translateFloatCoordinatesToWorld: (float *)localCoordinates forModel: aModel direction: (int)direction
287 {
288 double *domainSize = [spatialModel domainSize];
289 int i, dims = [spatialModel dimensions];
290
291 // single precision can be problematic so tweak with EPS if necessary
292 localCoordinates[0] = tweak_precision_direction(localCoordinates[0], coordinates[1] * domainSize[0], direction_round(direction, 0));
293 localCoordinates[1] = tweak_precision_direction(localCoordinates[1], (topologySize[1] - coordinates[0] - 1) * domainSize[0], direction_round(direction, 1));
294 if (dims == 3) localCoordinates[2] = tweak_precision_direction(localCoordinates[2], coordinates[2] * domainSize[2], direction_round(direction, 2));
295 }
296
297 - (void)translateDoubleCoordinatesToWorld: (double *)localCoordinates forModel: aModel
298 {
299 double *domainSize = [spatialModel domainSize];
300 int i, dims = [spatialModel dimensions];
301
302 localCoordinates[0] = coordinates[1] * domainSize[0] + localCoordinates[0];
303 localCoordinates[1] = (topologySize[1] - coordinates[0] - 1) * domainSize[0] + localCoordinates[1];
304 if (dims == 3) localCoordinates[2] = coordinates[2] * domainSize[2] + localCoordinates[2];
305 }
306
307 - (void)translateIntegerCoordinatesToSubspace: (int *)worldCoordinates forModel: aModel
308 {
309 int *gridSize = [spatialModel gridSize];
310 int i, dims = [spatialModel dimensions];
311
312 worldCoordinates[0] = worldCoordinates[0] - coordinates[1] * gridSize[0];
313 worldCoordinates[1] = worldCoordinates[1] - (topologySize[1] - coordinates[0] - 1) * gridSize[1];
314 if (dims == 3) worldCoordinates[2] = worldCoordinates[2] - coordinates[2] * gridSize[2];
315 }
316
317 - (void)translateFloatCoordinatesToSubspace: (float *)worldCoordinates forModel: aModel
318 {
319 double *domainSize = [spatialModel domainSize];
320 int i, dims = [spatialModel dimensions];
321
322 worldCoordinates[0] = worldCoordinates[0] - coordinates[1] * domainSize[0];
323 worldCoordinates[1] = worldCoordinates[1] - (topologySize[1] - coordinates[0] - 1) * domainSize[1];
324 if (dims == 3) worldCoordinates[2] = worldCoordinates[2] - coordinates[2] * domainSize[2];
325 }
326
327 - (void)translateDoubleCoordinatesToSubspace: (double *)worldCoordinates forModel: aModel
328 {
329 double *domainSize = [spatialModel domainSize];
330 int i, dims = [spatialModel dimensions];
331
332 worldCoordinates[0] = worldCoordinates[0] - coordinates[1] * domainSize[0];
333 worldCoordinates[1] = worldCoordinates[1] - (topologySize[1] - coordinates[0] - 1) * domainSize[1];
334 if (dims == 3) worldCoordinates[2] = worldCoordinates[2] - coordinates[2] * domainSize[2];
335 }
336
337 @end
338
339 //
340 // GPU methods
341 //
342
343 @implementation BioSwarmMPIController (GPURun)
344
345 - (void)allocateHierarchyGPUData: (int)numInstances
346 {
347 // allocate standard GPU memory
348 [super allocateHierarchyGPUData: numInstances];
349
350 // allocate GPU-MPI memory
351 id aModel = [modelInstances objectAtIndex: 0];
352 int numSubModels = [aModel numModelsInHierarchy];
353 [aModel allocateHierarchyGPUMPIData: gpuData modelInstances: numInstances];
354 }
355
356 @end
357
358 //
359 // MPI coordinates system
360 // fortran matrix style, x coordinate is rows (up, down), y coordinate is columns (left, right)
361 // origin in upper, left, front corner
362 //
363 void bc_mpi_sector_coordinate_shift(int xc, int yc, int zc, int direction, int *X, int *Y, int *Z)
364 {
365 *X = xc;
366 *Y = yc;
367 *Z = zc;
368 switch (direction) {
369 case S_LEFT:
370 *Y = *Y - 1; break;
371 case S_RIGHT:
372 *Y = *Y + 1; break;
373 case S_UP:
374 *X = *X - 1; break;
375 case S_DOWN:
376 *X = *X + 1; break;
377 case S_LEFT_UP:
378 *X = *X - 1; *Y = *Y - 1; break;
379 case S_RIGHT_UP:
380 *X = *X - 1; *Y = *Y + 1; break;
381 case S_LEFT_DOWN:
382 *X = *X + 1; *Y = *Y - 1; break;
383 case S_RIGHT_DOWN:
384 *X = *X + 1; *Y = *Y + 1; break;
385 case S_FRONT:
386 *Z = *Z - 1; break;
387 case S_BACK:
388 *Z = *Z + 1; break;
389 case S_LEFT_FRONT:
390 *Y = *Y - 1; *Z = *Z - 1; break;
391 case S_RIGHT_FRONT:
392 *Y = *Y + 1; *Z = *Z - 1; break;
393 case S_LEFT_BACK:
394 *Y = *Y - 1; *Z = *Z + 1; break;
395 case S_RIGHT_BACK:
396 *Y = *Y + 1; *Z = *Z + 1; break;
397 case S_UP_FRONT:
398 *X = *X - 1; *Z = *Z - 1; break;
399 case S_DOWN_FRONT:
400 *X = *X + 1; *Z = *Z - 1; break;
401 case S_UP_BACK:
402 *X = *X - 1; *Z = *Z + 1; break;
403 case S_DOWN_BACK:
404 *X = *X + 1; *Z = *Z + 1; break;
405 case S_LEFT_UP_FRONT:
406 *X = *X - 1; *Y = *Y - 1; *Z = *Z - 1; break;
407 case S_RIGHT_UP_FRONT:
408 *X = *X - 1; *Y = *Y + 1; *Z = *Z - 1; break;
409 case S_LEFT_DOWN_FRONT:
410 *X = *X + 1; *Y = *Y - 1; *Z = *Z - 1; break;
411 case S_RIGHT_DOWN_FRONT:
412 *X = *X + 1; *Y = *Y + 1; *Z = *Z - 1; break;
413 case S_LEFT_UP_BACK:
414 *X = *X - 1; *Y = *Y - 1; *Z = *Z + 1; break;
415 case S_RIGHT_UP_BACK:
416 *X = *X - 1; *Y = *Y + 1; *Z = *Z + 1; break;
417 case S_LEFT_DOWN_BACK:
418 *X = *X + 1; *Y = *Y - 1; *Z = *Z + 1; break;
419 case S_RIGHT_DOWN_BACK:
420 *X = *X + 1; *Y = *Y + 1; *Z = *Z + 1; break;
421 }
422 }
423
424 //
425 // GPU-MPI data transfer interface functions
426 //
427
428 void BC_MPI_initiate_transfer(void *model, int mode, double currentTime, double nextTime)
429 {
430 id modelObject = (id)model;
431 int modelNumber = [modelObject modelNumber];
432 id modelController = [modelObject modelController];
433 ModelGPUData **gpuData = [modelController gpuData];
434
435 printf("initiate transfer MPI data: %lf\n", currentTime);
436 [modelObject initiateGPUMPIData: gpuData[modelNumber] mode: mode];
437 }
438
439 void BC_MPI_finalize_transfer(void *model, int mode, double currentTime, double nextTime)
440 {
441 id modelObject = (id)model;
442
443 printf("finalize transfer MPI data: %lf\n", currentTime);
444 [modelObject finalizeGPUMPIData: NULL mode: mode];
445 }
446
447 void BC_MPI_record_time(void *model, const char *msg)
448 {
449 id modelObject = (id)model;
450 id modelController = [modelObject modelController];
451
452 if (TIME_PROFILE) bc_record_time([[NSString stringWithFormat: @"%s, rank = %d", msg, [modelController rank]] UTF8String]);
453 }
454
455 //
456 // MPI enhancements to standard BioSwarmModel
457 //
458
459 @implementation BioSwarmModel (MPI)
460
461 - (void)allocateMPIHierarchyWithEncode:(NSString *)anEncode
462 {
463 int i;
464
465 printf("allocateMPIHierarchyWithEncode:\n");
466
467 if ([anEncode isEqual: [BioSwarmModel intEncode]]) {
468 encodeTag = BCENCODE_INT;
469 } else if ([anEncode isEqual: [BioSwarmModel longEncode]]) {
470 encodeTag = BCENCODE_LONG;
471 } else if ([anEncode isEqual: [BioSwarmModel floatEncode]]) {
472 encodeTag = BCENCODE_FLOAT;
473 } else if ([anEncode isEqual: [BioSwarmModel doubleEncode]]) {
474 encodeTag = BCENCODE_DOUBLE;
475 } else {
476 // throw exception or something
477 NSLog(@"ERROR: BioSwarmModel unsupported encoding %@\n", anEncode);
478 return;
479 }
480 dataEncode = [anEncode retain];
481
482 if (numOfSpecies) {
483 switch (encodeTag) {
484 case BCENCODE_INT:
485 speciesData = malloc(sizeof(int) * numOfSpecies);
486 for (i = 0; i < numOfSpecies; ++i) ((int *)speciesData)[i] = 0;
487 break;
488 case BCENCODE_LONG:
489 speciesData = malloc(sizeof(long) * numOfSpecies);
490 for (i = 0; i < numOfSpecies; ++i) ((long *)speciesData)[i] = 0;
491 break;
492 case BCENCODE_FLOAT:
493 speciesData = malloc(sizeof(float) * numOfSpecies);
494 for (i = 0; i < numOfSpecies; ++i) ((float *)speciesData)[i] = 0;
495 break;
496 case BCENCODE_DOUBLE:
497 speciesData = malloc(sizeof(double) * numOfSpecies);
498 for (i = 0; i < numOfSpecies; ++i) ((double *)speciesData)[i] = 0;
499 break;
500 }
501 }
502
503 [self allocateMPISimulationWithEncode: anEncode];
504 for (i = 0; i < [childModels count]; ++i)
505 [[childModels objectAtIndex: i] allocateMPIHierarchyWithEncode: anEncode];
506 }
507
508 - (void)initializeMPIHierarchy
509 {
510 int i;
511 printf("initializeMPIHierarchy\n");
512
513 [self initializeMPISimulation];
514 for (i = 0; i < [childModels count]; ++i)
515 [[childModels objectAtIndex: i] initializeMPIHierarchy];
516 }
517
518 // these are implemented by subclasses
519 - (void)allocateMPISimulationWithEncode: (NSString *)anEncode
520 {
521 }
522
523 - (void)initializeMPISimulation
524 {
525 }
526
527 - (void)initializeMPIFileCompleted
528 {
529 }
530
531 @end
532
533 //
534 // SubcellularElementMethod
535 //
536 @implementation SubcellularElementModel (GPUMPIRun)
537
538 - (void)allocGPUMPIKernel: (void *)data
539 {
540 ModelGPUData *gpuData = (ModelGPUData *)data;
541 SEMfunctions *gpuFunctions = (SEMfunctions *)gpuData->gpuFunctions;
542 SpatialModel *spatialModel = [modelController spatialModel];
543 int dims = [spatialModel dimensions];
544 int neighborSectorSize;
545 if (dims == 2) neighborSectorSize = sectorsInDomain;
546 else neighborSectorSize = sectorsInDomain*sectorsInDomain;
547
548 gpuData->mpiData = (gpuFunctions->allocGPUMPIKernel)(self, gpuData->gpuPtrs, neighborSectorSize, sectorsInDomain, incomingSectors, outgoingSectors);
549 printf("allocGPUMPIKernel: %d %d %p %p %p\n", neighborSectorSize, sectorsInDomain, gpuData, gpuFunctions, gpuData->mpiData);
550 }
551
552 - (void)transferGPUMPIKernel: (void *)data toGPU: (BOOL)aFlag
553 {
554 printf("transferGPUMPIKernel:\n");
555 ModelGPUData *gpuData = (ModelGPUData *)data;
556 SEMfunctions *gpuFunctions = (SEMfunctions *)gpuData->gpuFunctions;
557
558 printf("transferGPUMPIKernel: %d %p %p %p\n", aFlag, gpuData, gpuFunctions, gpuData->mpiData);
559 (gpuFunctions->initGPUMPIKernel)(self, gpuData->mpiData, aFlag, gpuData->gpuPtrs, incomingElements, outgoingElements, totalElements);
560 }
561
562 - (void)initiateGPUMPIData: (void *)data mode: (int)aMode
563 {
564 SpatialModel *spatialModel = [modelController spatialModel];
565 int dims = [spatialModel dimensions];
566 int *neighborRanks = [modelController neighborRanks];
567 MPI_Comm comm = [modelController topologyMPICommunicator];
568 int direction;
569 int i, j, k;
570
571 int neighborSize, neighborSectorSize;
572 if (dims == 2) { neighborSize = S_XY_TOT; neighborSectorSize = sectorsInDomain; }
573 else { neighborSize = S_XYZ_TOT; neighborSectorSize = sectorsInDomain*sectorsInDomain; }
574
575 float (*oElements)[neighborSize][4 * neighborSectorSize * maxElementsPerSector] = (void *)outgoingElements;
576 int (*oSectors)[neighborSectorSize][neighborSize] = (void *)outgoingSectors;
577 float (*iElements)[neighborSize][4 * neighborSectorSize * maxElementsPerSector] = (void *)incomingElements;
578 int (*iSectors)[neighborSectorSize][neighborSize] = (void *)incomingSectors;
579
580 switch (aMode) {
581
582 case BC_MPI_SEM_MODE_INIT: {
583 printf("BC_MPI_SEM_MODE_INIT\n");
584 if (TIME_PROFILE) bc_record_time([[NSString stringWithFormat: @"start MPI_INIT rank = %d", [modelController rank]] UTF8String]);
585
586 // element map for each MPI neighbor
587 elementMapping = (id *)malloc(neighborSize * sizeof(id));
588 markedElements = (id *)malloc(neighborSize * sizeof(id));
589 for (direction = S_START; direction < neighborSize; ++direction) {
590 elementMapping[direction] = [NSMutableDictionary new];
591 markedElements[direction] = [NSMutableDictionary new];
592 }
593
594 if (receiveRequests) free(receiveRequests);
595 receiveRequests = malloc(neighborSize * sizeof(MPI_Request));
596 if (sendRequests) free(sendRequests);
597 sendRequests = malloc(neighborSize * sizeof(MPI_Request));
598 MPI_Request (*recvRequest)[neighborSize] = receiveRequests;
599 MPI_Request (*sendRequest)[neighborSize] = sendRequests;
600
601 // initiate the data receive
602 // max size plus 1 int for end of data flag
603 void *receiveElements = malloc(sizeof(int) * 4 * neighborSectorSize * maxElementsPerSector * neighborSize + sizeof(int));
604 float (*rElements)[neighborSize][4 * neighborSectorSize * maxElementsPerSector] = (void *)receiveElements;
605
606 for (direction = S_START; direction < neighborSize; ++direction) {
607 int dataSize = 0;
608 switch (direction) {
609 case S_LEFT:
610 case S_RIGHT:
611 case S_UP:
612 case S_DOWN:
613 case S_FRONT:
614 case S_BACK:
615 dataSize = 4 * neighborSectorSize * maxElementsPerSector;
616 break;
617
618 case S_LEFT_UP:
619 case S_RIGHT_UP:
620 case S_LEFT_DOWN:
621 case S_RIGHT_DOWN:
622 case S_LEFT_FRONT:
623 case S_RIGHT_FRONT:
624 case S_LEFT_BACK:
625 case S_RIGHT_BACK:
626 case S_UP_FRONT:
627 case S_DOWN_FRONT:
628 case S_UP_BACK:
629 case S_DOWN_BACK:
630 if (dims == 2) dataSize = 4 * maxElementsPerSector;
631 else dataSize = 4 * sectorsInDomain * maxElementsPerSector;
632 break;
633
634 case S_LEFT_UP_FRONT:
635 case S_RIGHT_UP_FRONT:
636 case S_LEFT_DOWN_FRONT:
637 case S_RIGHT_DOWN_FRONT:
638 case S_LEFT_UP_BACK:
639 case S_RIGHT_UP_BACK:
640 case S_LEFT_DOWN_BACK:
641 case S_RIGHT_DOWN_BACK:
642 dataSize = 4 * maxElementsPerSector;
643 break;
644 }
645 ++dataSize; // end flag
646 //printf("receive: rank = %d, %s neighbor rank = %d, size = %d\n", [modelController rank], bc_sector_direction_string(direction), neighborRanks[direction], dataSize);
647 MPI_Irecv(&((*rElements)[direction][0]), dataSize, MPI_FLOAT, neighborRanks[direction], 1, comm, &(*recvRequest)[direction]);
648 }
649
650 // pack the elements
651 // No GPU data yet so we use CPU data
652 void *packElements = malloc(sizeof(int) * 4 * neighborSectorSize * maxElementsPerSector * neighborSize + sizeof(int));
653 float (*pElements)[neighborSize][4 * neighborSectorSize * maxElementsPerSector] = packElements;
654 int currentPos[neighborSize];
655 int (*seTable)[maxElementsPerSector][maxSectors] = (void *)sectorElementTable;
656
657 // TODO: need to handle float/double switch
658 Grid2DArray *X = [[elementData gridWithName: @"X"] objectForKey: @"matrix"];
659 Grid2DArray *Y = [[elementData gridWithName: @"Y"] objectForKey: @"matrix"];
660 Grid2DArray *Z = [[elementData gridWithName: @"Z"] objectForKey: @"matrix"];
661 float (*xGrid)[maxElements] = [X getMatrix];
662 float (*yGrid)[maxElements] = [Y getMatrix];
663 float (*zGrid)[maxElements] = [Z getMatrix];
664 float coords[3];
665
666 for (direction = S_START; direction < neighborSize; ++direction) {
667 currentPos[direction] = 0;
668 for (i = 0; i < neighborSectorSize; ++i) {
669 int sectorNum = (*oSectors)[i][direction];
670 if (sectorNum == -1) continue;
671
672 for (j = 0; j < maxElementsPerSector; ++j) {
673 int elemNum = (*seTable)[j][sectorNum];
674 if (elemNum == -1) continue;
675
676 // convert coordinates to world
677 coords[0] = (*xGrid)[elemNum]; coords[1] = (*yGrid)[elemNum]; coords[2] = (*zGrid)[elemNum];
678 if ([spatialModel isSubspace]) [spatialModel translateFloatCoordinatesToWorld: coords forModel: self direction: direction];
679
680 (*pElements)[direction][currentPos[direction]] = (float)elemNum;
681 (*pElements)[direction][currentPos[direction]+1] = coords[0];
682 (*pElements)[direction][currentPos[direction]+2] = coords[1];
683 (*pElements)[direction][currentPos[direction]+3] = coords[2];
684 currentPos[direction] += 4;
685 }
686 }
687 printf("rank = %d, direction = %d, pack = %d\n", [modelController rank], direction, currentPos[direction] / 4);
688 // add end of data flag
689 (*pElements)[direction][currentPos[direction]] = -1.0;
690 ++currentPos[direction];
691 }
692
693 // send the data
694 for (direction = S_START; direction < neighborSize; ++direction) {
695 MPI_Isend(&((*pElements)[direction][0]), currentPos[direction], MPI_FLOAT, neighborRanks[direction], 1, comm, &(*sendRequest)[direction]);
696 }
697
698 printf("rank = %d, wait on send\n", [modelController rank]);
699 MPI_Waitall(neighborSize - 1, &(*sendRequest)[S_START], MPI_STATUSES_IGNORE);
700 printf("rank = %d, wait on receive\n", [modelController rank]);
701 MPI_Waitall(neighborSize - 1, &(*recvRequest)[S_START], MPI_STATUSES_IGNORE);
702
703 // We assume the elements and sectors do not exist
704
705 for (direction = S_START; direction < neighborSize; ++direction) {
706 for (i = 0; i < neighborSectorSize * maxElementsPerSector; ++i) {
707 if (neighborRanks[direction] == MPI_PROC_NULL) break;
708
709 int elemNum = (*rElements)[direction][4*i];
710 if (elemNum == -1) break;
711
712 // convert coordinates to local subspace
713 coords[0] = (*rElements)[direction][4*i+1];
714 coords[1] = (*rElements)[direction][4*i+2];
715 coords[2] = (*rElements)[direction][4*i+3];
716 if ([spatialModel isSubspace]) [spatialModel translateFloatCoordinatesToSubspace: coords forModel: self];
717 //if ([modelController rank] == 2) printf("remote (%d) = (%f %f %f), (%f %f %f)\n", elemNum, (*rElements)[direction][4*i+1], (*rElements)[direction][4*i+2], (*rElements)[direction][4*i+3], coords[0], coords[1], coords[2]);
718
719 if (totalElements >= maxElements) {
720 printf("ERROR: exceeded maximum number of elements: %d >= %d\n", totalElements, maxElements);
721 abort();
722 }
723
724 // add element
725 (*xGrid)[totalElements] = coords[0]; (*yGrid)[totalElements] = coords[1]; (*zGrid)[totalElements] = coords[2];
726 [elementMapping[direction] setObject: [NSNumber numberWithInt: totalElements] forKey: [NSNumber numberWithInt: elemNum]];
727 ++totalElements;
728 }
729 printf("rank = %d, direction = %d, elements = %d %f\n", [modelController rank], direction, i, (*rElements)[direction][0]);
730 }
731
732 free(receiveElements);
733 free(packElements);
734
735 if (TIME_PROFILE) bc_record_time([[NSString stringWithFormat: @"end MPI_INIT rank = %d", [modelController rank]] UTF8String]);
736 break;
737 }
738
739 case BC_MPI_SEM_MODE_F1: {
740 printf("BC_MPI_SEM_MODE_F1\n");
741 if (TIME_PROFILE) bc_record_time([[NSString stringWithFormat: @"start MPI_F1 rank = %d", [modelController rank]] UTF8String]);
742
743 if (receiveRequests) free(receiveRequests);
744 receiveRequests = malloc(neighborSize * sizeof(MPI_Request));
745 if (sendRequests) free(sendRequests);
746 sendRequests = malloc(neighborSize * sizeof(MPI_Request));
747 MPI_Request (*recvRequest)[neighborSize] = receiveRequests;
748 MPI_Request (*sendRequest)[neighborSize] = sendRequests;
749
750 //printf("pointers: %p %p %p %p\n", outgoingElements, outgoingSectors, incomingElements, incomingSectors);
751
752 // receive data from MPI neighbors
753 // we use non-blocking communication
754
755 // max size plus 1 int for end of data flag
756 void *receiveElements = malloc(sizeof(int) * 4 * neighborSectorSize * maxElementsPerSector * neighborSize + sizeof(int));
757 float (*rElements)[neighborSize][4 * neighborSectorSize * maxElementsPerSector] = (void *)receiveElements;
758
759 for (direction = S_START; direction < neighborSize; ++direction) {
760 int dataSize = 0;
761 switch (direction) {
762 case S_LEFT:
763 case S_RIGHT:
764 case S_UP:
765 case S_DOWN:
766 case S_FRONT:
767 case S_BACK:
768 dataSize = 4 * neighborSectorSize * maxElementsPerSector;
769 break;
770
771 case S_LEFT_UP:
772 case S_RIGHT_UP:
773 case S_LEFT_DOWN:
774 case S_RIGHT_DOWN:
775 case S_LEFT_FRONT:
776 case S_RIGHT_FRONT:
777 case S_LEFT_BACK:
778 case S_RIGHT_BACK:
779 case S_UP_FRONT:
780 case S_DOWN_FRONT:
781 case S_UP_BACK:
782 case S_DOWN_BACK:
783 if (dims == 2) dataSize = 4 * maxElementsPerSector;
784 else dataSize = 4 * sectorsInDomain * maxElementsPerSector;
785 break;
786
787 case S_LEFT_UP_FRONT:
788 case S_RIGHT_UP_FRONT:
789 case S_LEFT_DOWN_FRONT:
790 case S_RIGHT_DOWN_FRONT:
791 case S_LEFT_UP_BACK:
792 case S_RIGHT_UP_BACK:
793 case S_LEFT_DOWN_BACK:
794 case S_RIGHT_DOWN_BACK:
795 dataSize = 4 * maxElementsPerSector;
796 break;
797 }
798 ++dataSize; // end flag
799 printf("receive: rank = %d, %s neighbor rank = %d, size = %d\n", [modelController rank], bc_sector_direction_string(direction), neighborRanks[direction], dataSize);
800 MPI_Irecv(&((*rElements)[direction][0]), dataSize, MPI_FLOAT, neighborRanks[direction], 1, comm, &(*recvRequest)[direction]);
801 }
802
803 // gather data from GPU
804 [self transferGPUMPIKernel: data toGPU: NO];
805
806 // max size plus one int for end marker
807 int currentPos[neighborSize];
808 void *packElements = malloc(sizeof(int) * 4 * neighborSectorSize * maxElementsPerSector * neighborSize + sizeof(int));
809 float (*pElements)[neighborSize][4 * neighborSectorSize * maxElementsPerSector] = packElements;
810 NSMutableArray *openElements = [NSMutableArray array];
811
812 for (direction = S_START; direction < neighborSize; ++direction) {
813 currentPos[direction] = 0;
814 for (i = 0; i < neighborSectorSize; ++i) {
815 int sectorNum = (*oSectors)[i][direction];
816 if (sectorNum == -1) continue;
817
818 int idx = 4 * i * maxElementsPerSector;
819 for (j = 0; j < maxElementsPerSector; ++j) {
820 int elemNum = (int)(*oElements)[direction][idx+4*j];
821 if (elemNum == -1) continue;
822
823 // found element
824 [openElements addObject: [NSNumber numberWithInt: elemNum]];
825 (*pElements)[direction][currentPos[direction]] = (float)elemNum;
826 (*pElements)[direction][currentPos[direction]+1] = (*oElements)[direction][idx+4*j+1];
827 (*pElements)[direction][currentPos[direction]+2] = (*oElements)[direction][idx+4*j+2];
828 (*pElements)[direction][currentPos[direction]+3] = (*oElements)[direction][idx+4*j+3];
829 currentPos[direction] += 4;
830 }
831 }
832 printf("rank = %d, direction = %d, pack = %d\n", [modelController rank], direction, currentPos[direction] / 4);
833 (*pElements)[direction][currentPos[direction]] = -1.0;
834 ++currentPos[direction];
835 }
836 printf("rank = %d, total pack = %d\n", [modelController rank], [openElements count]);
837
838 // send data to MPI neighbors
839 for (direction = S_START; direction < neighborSize; ++direction) {
840 MPI_Isend(&((*pElements)[direction][0]), currentPos[direction], MPI_FLOAT, neighborRanks[direction], 1, comm, &(*sendRequest)[direction]);
841 }
842
843 printf("rank = %d, wait on send\n", [modelController rank]);
844 MPI_Waitall(neighborSize - 1, &(*sendRequest)[S_START], MPI_STATUSES_IGNORE);
845 printf("rank = %d, wait on receive\n", [modelController rank]);
846 MPI_Waitall(neighborSize - 1, &(*recvRequest)[S_START], MPI_STATUSES_IGNORE);
847
848 // for (direction = S_START; direction < neighborSize; ++direction) {
849 // printf("rank = %d, receive (%d, %d) = %f\n", [modelController rank], direction, neighborRanks[direction], (*iElements)[direction][0]);
850 //}
851
852 // Extract elements and copy data for GPU transfer
853 // This data is not in any sector order, but we do not care because
854 // we just update data based upon the element number.
855 // The sector numbers do not change for F1.
856 for (direction = S_START; direction < neighborSize; ++direction) {
857 currentPos[direction] = 0;
858 for (i = 0; i < 4 * neighborSectorSize * maxElementsPerSector; ++i) (*iElements)[direction][i] = -1;
859
860 for (i = 0; i < neighborSectorSize * maxElementsPerSector; ++i) {
861 if (neighborRanks[direction] == MPI_PROC_NULL) break;
862
863 int remoteElem = (*rElements)[direction][4*i];
864 if (remoteElem == -1) break;
865
866 // search for element in map
867 NSNumber *n = [elementMapping[direction] objectForKey: [NSNumber numberWithInt: remoteElem]];
868 if (!n) printf("MPI_ERROR: could not find matching element in map (%d, %f, %f, %f)\n", remoteElem, (*rElements)[direction][4*i+1], (*rElements)[direction][4*i+2], (*rElements)[direction][4*i+3]);
869 else {
870 // copy data for GPU transfer
871 int elemNum = [n intValue];
872 (*iElements)[direction][currentPos[direction]] = (float)elemNum;
873 (*iElements)[direction][currentPos[direction]+1] = (*rElements)[direction][4*i+1];
874 (*iElements)[direction][currentPos[direction]+2] = (*rElements)[direction][4*i+2];
875 (*iElements)[direction][currentPos[direction]+3] = (*rElements)[direction][4*i+3];
876 currentPos[direction] += 4;
877
878 if (elemNum == 308) printf("rank = %d, incoming element map %d = %d, %f, %f, %f\n", [modelController rank], elemNum, remoteElem, (*rElements)[direction][4*i+1], (*rElements)[direction][4*i+2], (*rElements)[direction][4*i+3]);
879 }
880 }
881 //printf("rank = %d, direction = %d\n", [modelController rank], direction);
882 }
883
884 // transfer data to GPU
885 [self transferGPUMPIKernel: data toGPU: YES];
886
887 free(receiveElements);
888 free(packElements);
889
890 if (TIME_PROFILE) bc_record_time([[NSString stringWithFormat: @"end MPI_F1 rank = %d", [modelController rank]] UTF8String]);
891 break;
892 }
893
894 case BC_MPI_SEM_MODE_MOVE_READ_ONLY: {
895 printf("BC_MPI_SEM_MODE_MOVE_READ_ONLY\n");
896 if (TIME_PROFILE) bc_record_time([[NSString stringWithFormat: @"start MPI_MOVE_READ_ONLY rank = %d", [modelController rank]] UTF8String]);
897
898 if (receiveRequests) free(receiveRequests);
899 receiveRequests = malloc(neighborSize * sizeof(MPI_Request));
900 if (sendRequests) free(sendRequests);
901 sendRequests = malloc(neighborSize * sizeof(MPI_Request));
902 MPI_Request (*recvRequest)[neighborSize] = receiveRequests;
903 MPI_Request (*sendRequest)[neighborSize] = sendRequests;
904
905 // max size plus 1 int for end of data flag
906 void *receiveElements = malloc(sizeof(int) * 4 * neighborSectorSize * maxElementsPerSector * neighborSize + sizeof(int));
907 float (*rElements)[neighborSize][4 * neighborSectorSize * maxElementsPerSector] = (void *)receiveElements;
908
909 for (direction = S_START; direction < neighborSize; ++direction) {
910 //printf("receive: rank = %d, %s neighbor rank = %d, %p\n", [modelController rank], bc_sector_direction_string(direction), neighborRanks[direction], &((*iElements)[direction][0]));
911 int dataSize = 0;
912 switch (direction) {
913 case S_LEFT:
914 case S_RIGHT:
915 case S_UP:
916 case S_DOWN:
917 case S_FRONT:
918 case S_BACK:
919 dataSize = 4 * neighborSectorSize * maxElementsPerSector;
920 break;
921
922 case S_LEFT_UP:
923 case S_RIGHT_UP:
924 case S_LEFT_DOWN:
925 case S_RIGHT_DOWN:
926 case S_LEFT_FRONT:
927 case S_RIGHT_FRONT:
928 case S_LEFT_BACK:
929 case S_RIGHT_BACK:
930 case S_UP_FRONT:
931 case S_DOWN_FRONT:
932 case S_UP_BACK:
933 case S_DOWN_BACK:
934 if (dims == 2) dataSize = 4 * maxElementsPerSector;
935 else dataSize = 4 * sectorsInDomain * maxElementsPerSector;
936 break;
937
938 case S_LEFT_UP_FRONT:
939 case S_RIGHT_UP_FRONT:
940 case S_LEFT_DOWN_FRONT:
941 case S_RIGHT_DOWN_FRONT:
942 case S_LEFT_UP_BACK:
943 case S_RIGHT_UP_BACK:
944 case S_LEFT_DOWN_BACK:
945 case S_RIGHT_DOWN_BACK:
946 dataSize = 4 * maxElementsPerSector;
947 break;
948 }
949 ++dataSize; // end flag
950 MPI_Irecv(&((*rElements)[direction][0]), dataSize, MPI_FLOAT, neighborRanks[direction], 1, comm, &(*recvRequest)[direction]);
951 }
952
953 // gather data from GPU
954 [self transferGPUMPIKernel: data toGPU: NO];
955
956 // pack data
957 // max size plus one int for end marker
958 int currentPos[neighborSize];
959 void *packElements = malloc(sizeof(int) * 4 * neighborSectorSize * maxElementsPerSector * neighborSize + sizeof(int));
960 float (*pElements)[neighborSize][4 * neighborSectorSize * maxElementsPerSector] = packElements;
961 float coords[3];
962
963 for (direction = S_START; direction < neighborSize; ++direction) {
964 currentPos[direction] = 0;
965 for (i = 0; i < neighborSectorSize; ++i) {
966 int sectorNum = (*oSectors)[i][direction];
967 if (sectorNum == -1) continue;
968
969 int idx = 4 * i * maxElementsPerSector;
970 for (j = 0; j < maxElementsPerSector; ++j) {
971 int elemNum = (int)(*oElements)[direction][idx+4*j];
972 if (elemNum == -1) continue;
973
974 // convert coordinates to world
975 coords[0] = (*oElements)[direction][idx+4*j+1]; coords[1] = (*oElements)[direction][idx+4*j+2]; coords[2] = (*oElements)[direction][idx+4*j+3];
976 if ([spatialModel isSubspace]) [spatialModel translateFloatCoordinatesToWorld: coords forModel: self direction: direction];
977
978 printf("MOVE_READ_ONLY: rank = %d, direction = %d, send element %d, %.15f, %.15f, %.15f (%.15f, %.15f, %.15f)\n", [modelController rank], direction, elemNum, coords[0], coords[1], coords[2], (*oElements)[direction][idx+4*j+1], (*oElements)[direction][idx+4*j+2], (*oElements)[direction][idx+4*j+3]);
979
980 // found element
981 (*pElements)[direction][currentPos[direction]] = (float)elemNum;
982 (*pElements)[direction][currentPos[direction]+1] = coords[0];
983 (*pElements)[direction][currentPos[direction]+2] = coords[1];
984 (*pElements)[direction][currentPos[direction]+3] = coords[2];
985 currentPos[direction] += 4;
986 }
987 }
988 printf("MOVE_READ_ONLY: rank = %d, direction = %d, pack = %d\n", [modelController rank], direction, currentPos[direction] / 4);
989 (*pElements)[direction][currentPos[direction]] = -1.0;
990 ++currentPos[direction];
991 }
992
993 // send data to MPI neighbors
994
995 for (direction = S_START; direction < neighborSize; ++direction) {
996 //printf("send: rank = %d, %s neighbor rank = %d\n", [modelController rank], bc_sector_direction_string(direction), neighborRanks[direction]);
997 MPI_Isend(&((*pElements)[direction][0]), currentPos[direction], MPI_FLOAT, neighborRanks[direction], 1, comm, &(*sendRequest)[direction]);
998 }
999
1000 printf("MOVE_READ_ONLY: rank = %d, wait on send\n", [modelController rank]);
1001 MPI_Waitall(neighborSize - 1, &(*sendRequest)[S_START], MPI_STATUSES_IGNORE);
1002 printf("MOVE_READ_ONLY: rank = %d, wait on receive\n", [modelController rank]);
1003 MPI_Waitall(neighborSize - 1, &(*recvRequest)[S_START], MPI_STATUSES_IGNORE);
1004
1005 //
1006 // Given a set of elements that have moved from a read-only border sector
1007 // into a local sector. Each element has essentially migrated from one
1008 // MPI node to another, it was under our neighbor's control before but
1009 // now we control it. As an element in a read-only border sector, it
1010 // has an entry in the element map and it has a local element number.
1011 // We remove it from the element map, let it keep its local element
1012 // number, and update the sector tables to include that element.
1013 //
1014 for (direction = S_START; direction < neighborSize; ++direction)
1015 [markedElements[direction] removeAllObjects];
1016
1017 for (direction = S_START; direction < neighborSize; ++direction) {
1018 currentPos[direction] = 0;
1019 for (i = 0; i < 4 * neighborSectorSize * maxElementsPerSector; ++i) (*iElements)[direction][i] = -1;
1020
1021 for (i = 0; i < neighborSectorSize * maxElementsPerSector; ++i) {
1022 if (neighborRanks[direction] == MPI_PROC_NULL) break;
1023
1024 int remoteElem = (*rElements)[direction][4*i];
1025 if (remoteElem == -1) break;
1026
1027 // search for element in map
1028 NSNumber *rn = [NSNumber numberWithInt: remoteElem];
1029 NSNumber *n = [elementMapping[direction] objectForKey: rn];
1030 if (!n) printf("MPI_ERROR: could not find matching element in map (%d, %f, %f, %f)\n", remoteElem, (*rElements)[direction][4*i+1], (*rElements)[direction][4*i+2], (*rElements)[direction][4*i+3]);
1031 else {
1032 // convert coordinates to local subspace
1033 coords[0] = (*rElements)[direction][4*i+1];
1034 coords[1] = (*rElements)[direction][4*i+2];
1035 coords[2] = (*rElements)[direction][4*i+3];
1036 if ([spatialModel isSubspace]) [spatialModel translateFloatCoordinatesToSubspace: coords forModel: self];
1037
1038 // copy data for GPU transfer
1039 int elemNum = [n intValue];
1040 (*iElements)[direction][currentPos[direction]] = (float)elemNum;
1041 (*iElements)[direction][currentPos[direction]+1] = coords[0];
1042 (*iElements)[direction][currentPos[direction]+2] = coords[1];
1043 (*iElements)[direction][currentPos[direction]+3] = coords[2];
1044 currentPos[direction] += 4;
1045
1046 [markedElements[direction] setObject: n forKey: rn];
1047 [elementMapping[direction] removeObjectForKey: rn];
1048 printf("MOVE_READ_ONLY: rank = %d, direction = %d, add incoming element %d = %d, %f, %f, %f\n", [modelController rank], direction, elemNum, remoteElem, (*rElements)[direction][4*i+1], (*rElements)[direction][4*i+2], (*rElements)[direction][4*i+3]);
1049 }
1050 }
1051 //printf("rank = %d, direction = %d\n", [modelController rank], direction);
1052 }
1053
1054 //
1055 // On the sending side, the MPI node has lost control of the element
1056 // so that local element number is available, however that element is
1057 // not really deleted because we expect is to show up in our neighbor
1058 // read-only sector border list. But we do not know the new element
1059 // number that will be assigned in our MPI neighbor. We need to
1060 // let our neighbors know their new element numbers so they can add to their map.
1061 //
1062
1063 for (direction = S_START; direction < neighborSize; ++direction) {
1064 //printf("receive: rank = %d, %s neighbor rank = %d, %p\n", [modelController rank], bc_sector_direction_string(direction), neighborRanks[direction], &((*iElements)[direction][0]));
1065 int dataSize = 0;
1066 switch (direction) {
1067 case S_LEFT:
1068 case S_RIGHT:
1069 case S_UP:
1070 case S_DOWN:
1071 case S_FRONT:
1072 case S_BACK:
1073 dataSize = 4 * neighborSectorSize * maxElementsPerSector;
1074 break;
1075
1076 case S_LEFT_UP:
1077 case S_RIGHT_UP:
1078 case S_LEFT_DOWN:
1079 case S_RIGHT_DOWN:
1080 case S_LEFT_FRONT:
1081 case S_RIGHT_FRONT:
1082 case S_LEFT_BACK:
1083 case S_RIGHT_BACK:
1084 case S_UP_FRONT:
1085 case S_DOWN_FRONT:
1086 case S_UP_BACK:
1087 case S_DOWN_BACK:
1088 if (dims == 2) dataSize = 4 * maxElementsPerSector;
1089 else dataSize = 4 * sectorsInDomain * maxElementsPerSector;
1090 break;
1091
1092 case S_LEFT_UP_FRONT:
1093 case S_RIGHT_UP_FRONT:
1094 case S_LEFT_DOWN_FRONT:
1095 case S_RIGHT_DOWN_FRONT:
1096 case S_LEFT_UP_BACK:
1097 case S_RIGHT_UP_BACK:
1098 case S_LEFT_DOWN_BACK:
1099 case S_RIGHT_DOWN_BACK:
1100 dataSize = 4 * maxElementsPerSector;
1101 break;
1102 }
1103 ++dataSize; // end flag
1104 MPI_Irecv(&((*rElements)[direction][0]), dataSize, MPI_FLOAT, neighborRanks[direction], 1, comm, &(*recvRequest)[direction]);
1105 }
1106
1107 // pack data
1108 for (direction = S_START; direction < neighborSize; ++direction) {
1109 currentPos[direction] = 0;
1110
1111 NSArray *allKeys = [markedElements[direction] allKeys];
1112 printf("MOVE_READ_ONLY: rank = %d, direction = %d, send elements for map = %d\n", [modelController rank], direction, [allKeys count]);
1113 for (i = 0; i < [allKeys count]; ++i) {
1114 NSNumber *rn = [allKeys objectAtIndex: i];
1115 (*pElements)[direction][currentPos[direction]] = [rn floatValue];
1116 (*pElements)[direction][currentPos[direction]+1] = [[markedElements[direction] objectForKey: rn] floatValue];
1117 currentPos[direction] += 2;
1118 }
1119 (*pElements)[direction][currentPos[direction]] = -1.0;
1120 ++currentPos[direction];
1121 }
1122
1123 // send data to MPI neighbors
1124
1125 for (direction = S_START; direction < neighborSize; ++direction) {
1126 MPI_Isend(&((*pElements)[direction][0]), currentPos[direction], MPI_FLOAT, neighborRanks[direction], 1, comm, &(*sendRequest)[direction]);
1127 }
1128
1129 printf("rank = %d, wait on send\n", [modelController rank]);
1130 MPI_Waitall(neighborSize - 1, &(*sendRequest)[S_START], MPI_STATUSES_IGNORE);
1131 printf("rank = %d, wait on receive\n", [modelController rank]);
1132 MPI_Waitall(neighborSize - 1, &(*recvRequest)[S_START], MPI_STATUSES_IGNORE);
1133
1134 // update element maps
1135 for (direction = S_START; direction < neighborSize; ++direction) {
1136 for (i = 0; i < neighborSectorSize * maxElementsPerSector; ++i) {
1137 if (neighborRanks[direction] == MPI_PROC_NULL) break;
1138
1139 int elemNum = (*rElements)[direction][2*i];
1140 if (elemNum == -1) break;
1141 int remoteElem = (*rElements)[direction][2*i+1];
1142 [elementMapping[direction] setObject: [NSNumber numberWithInt: elemNum] forKey: [NSNumber numberWithInt: remoteElem]];
1143 printf("MOVE_READ_ONLY: rank = %d, direction = %d, add element map %d = %d\n", [modelController rank], direction, remoteElem, elemNum);
1144 }
1145 }
1146
1147 // transfer data to GPU
1148 [self transferGPUMPIKernel: data toGPU: YES];
1149
1150 free(receiveElements);
1151 free(packElements);
1152
1153 if (TIME_PROFILE) bc_record_time([[NSString stringWithFormat: @"end MPI_MOVE_READ_ONLY rank = %d", [modelController rank]] UTF8String]);
1154 break;
1155 }
1156
1157 case BC_MPI_SEM_MODE_MOVE: {
1158 printf("BC_MPI_SEM_MODE_MOVE\n");
1159 if (TIME_PROFILE) bc_record_time([[NSString stringWithFormat: @"start MPI_MOVE rank = %d", [modelController rank]] UTF8String]);
1160
1161 if (receiveRequests) free(receiveRequests);
1162 receiveRequests = malloc(neighborSize * sizeof(MPI_Request));
1163 if (sendRequests) free(sendRequests);
1164 sendRequests = malloc(neighborSize * sizeof(MPI_Request));
1165 MPI_Request (*recvRequest)[neighborSize] = receiveRequests;
1166 MPI_Request (*sendRequest)[neighborSize] = sendRequests;
1167
1168 // max size plus 1 int for end of data flag
1169 void *receiveElements = malloc(sizeof(int) * 4 * neighborSectorSize * maxElementsPerSector * neighborSize + sizeof(int));
1170 float (*rElements)[neighborSize][4 * neighborSectorSize * maxElementsPerSector] = (void *)receiveElements;
1171
1172 for (direction = S_START; direction < neighborSize; ++direction) {
1173 //printf("receive: rank = %d, %s neighbor rank = %d, %p\n", [modelController rank], bc_sector_direction_string(direction), neighborRanks[direction], &((*iElements)[direction][0]));
1174 int dataSize = 0;
1175 switch (direction) {
1176 case S_LEFT:
1177 case S_RIGHT:
1178 case S_UP:
1179 case S_DOWN:
1180 case S_FRONT:
1181 case S_BACK:
1182 dataSize = 4 * neighborSectorSize * maxElementsPerSector;
1183 break;
1184
1185 case S_LEFT_UP:
1186 case S_RIGHT_UP:
1187 case S_LEFT_DOWN:
1188 case S_RIGHT_DOWN:
1189 case S_LEFT_FRONT:
1190 case S_RIGHT_FRONT:
1191 case S_LEFT_BACK:
1192 case S_RIGHT_BACK:
1193 case S_UP_FRONT:
1194 case S_DOWN_FRONT:
1195 case S_UP_BACK:
1196 case S_DOWN_BACK:
1197 if (dims == 2) dataSize = 4 * maxElementsPerSector;
1198 else dataSize = 4 * sectorsInDomain * maxElementsPerSector;
1199 break;
1200
1201 case S_LEFT_UP_FRONT:
1202 case S_RIGHT_UP_FRONT:
1203 case S_LEFT_DOWN_FRONT:
1204 case S_RIGHT_DOWN_FRONT:
1205 case S_LEFT_UP_BACK:
1206 case S_RIGHT_UP_BACK:
1207 case S_LEFT_DOWN_BACK:
1208 case S_RIGHT_DOWN_BACK:
1209 dataSize = 4 * maxElementsPerSector;
1210 break;
1211 }
1212 ++dataSize; // end flag
1213 MPI_Irecv(&((*rElements)[direction][0]), dataSize, MPI_FLOAT, neighborRanks[direction], 1, comm, &(*recvRequest)[direction]);
1214 }
1215
1216 // transfer data from GPU
1217 [self transferGPUMPIKernel: data toGPU: NO];
1218
1219 // pack data
1220 // max size plus one int for end marker
1221 int currentPos[neighborSize];
1222 void *packElements = malloc(sizeof(int) * 4 * neighborSectorSize * maxElementsPerSector * neighborSize + sizeof(int));
1223 float (*pElements)[neighborSize][4 * neighborSectorSize * maxElementsPerSector] = packElements;
1224 float coords[3];
1225
1226 for (direction = S_START; direction < neighborSize; ++direction) {
1227 currentPos[direction] = 0;
1228 for (i = 0; i < neighborSectorSize; ++i) {
1229 int sectorNum = (*oSectors)[i][direction];
1230 if (sectorNum == -1) continue;
1231
1232 int idx = 4 * i * maxElementsPerSector;
1233 for (j = 0; j < maxElementsPerSector; ++j) {
1234 int elemNum = (int)(*oElements)[direction][idx+4*j];
1235 if (elemNum == -1) continue;
1236
1237 // convert coordinates to world
1238 coords[0] = (*oElements)[direction][idx+4*j+1]; coords[1] = (*oElements)[direction][idx+4*j+2]; coords[2] = (*oElements)[direction][idx+4*j+3];
1239 if ([spatialModel isSubspace]) [spatialModel translateFloatCoordinatesToWorld: coords forModel: self direction: direction];
1240
1241 // found element
1242 (*pElements)[direction][currentPos[direction]] = (float)elemNum;
1243 (*pElements)[direction][currentPos[direction]+1] = coords[0];
1244 (*pElements)[direction][currentPos[direction]+2] = coords[1];
1245 (*pElements)[direction][currentPos[direction]+3] = coords[2];
1246 currentPos[direction] += 4;
1247 }
1248 }
1249 printf("BC_MPI_SEM_MODE_MOVE: rank = %d, direction = %d, pack = %d\n", [modelController rank], direction, currentPos[direction] / 4);
1250 (*pElements)[direction][currentPos[direction]] = -1.0;
1251 ++currentPos[direction];
1252 }
1253
1254 // send data to MPI neighbors
1255
1256 for (direction = S_START; direction < neighborSize; ++direction) {
1257 //printf("send: rank = %d, %s neighbor rank = %d, size = %d\n", [modelController rank], bc_sector_direction_string(direction), neighborRanks[direction], currentPos[direction]);
1258 MPI_Isend(&((*pElements)[direction][0]), currentPos[direction], MPI_FLOAT, neighborRanks[direction], 1, comm, &(*sendRequest)[direction]);
1259 }
1260
1261 printf("rank = %d, wait on send\n", [modelController rank]);
1262 MPI_Waitall(neighborSize - 1, &(*sendRequest)[S_START], MPI_STATUSES_IGNORE);
1263 printf("rank = %d, wait on receive\n", [modelController rank]);
1264 MPI_Waitall(neighborSize - 1, &(*recvRequest)[S_START], MPI_STATUSES_IGNORE);
1265
1266 //
1267 // Given the set of elements for the read-only border sectors, need to
1268 // determine elements in four categories.
1269 // (1) elements that stayed in sector
1270 // (2) elements moved from local sector to read-only border
1271 // (3) elements moved out of read-only border (removed)
1272 // (4) new elements in read-only border (insert)
1273 //
1274 // to determine these categores
1275 // (1) exists in element map
1276 // (2) exists in element map, put in map by MOVE_READ_ONLY abve
1277 // (3) exists in element map, not in given element set
1278 // (4) does not exist in element map
1279 //
1280
1281 NSMutableArray *insertElements[neighborSize];
1282 for (direction = S_START; direction < neighborSize; ++direction) {
1283 NSMutableDictionary *d = [NSMutableDictionary dictionaryWithDictionary: elementMapping[direction]];
1284 insertElements[direction] = [NSMutableArray array];
1285
1286 currentPos[direction] = 0;
1287 for (i = 0; i < 4 * neighborSectorSize * maxElementsPerSector; ++i) (*iElements)[direction][i] = -1;
1288
1289 for (i = 0; i < neighborSectorSize * maxElementsPerSector; ++i) {
1290 if (neighborRanks[direction] == MPI_PROC_NULL) break;
1291
1292 int remoteElem = (*rElements)[direction][4*i];
1293 if (remoteElem == -1) break;
1294
1295 // search for element in map
1296 NSNumber *rn = [NSNumber numberWithInt: remoteElem];
1297 NSNumber *n = [d objectForKey: rn];
1298 if (n) {
1299 // existing element
1300 [d removeObjectForKey: rn];
1301 //printf("rank = %d, existing incoming element %d = %d\n", [modelController rank], [n intValue], [rn intValue]);
1302 } else {
1303 // new element, save for later
1304 [insertElements[direction] addObject: rn];
1305 }
1306 }
1307
1308 // any remaining elements in map are category (3) and removed
1309 // these element numbers are available for reassignment
1310 NSArray *allKeys = [d allKeys];
1311 for (i = 0; i < [allKeys count]; ++i) {
1312 NSNumber *rn = [allKeys objectAtIndex: i];
1313 NSNumber *n = [d objectForKey: rn];
1314 [availableElements addObject: n];
1315 [elementMapping[direction] removeObjectForKey: rn];
1316 printf("MOVE: rank = %d, remove incoming element %d = %d\n", [modelController rank], [n intValue], [rn intValue]);
1317 }
1318
1319 // assign new elements a local element number
1320 for (i = 0; i < [insertElements[direction] count]; ++i) {
1321 NSNumber *rn = [insertElements[direction] objectAtIndex: i];
1322 if ([availableElements count] > 0) {
1323 // grab an available elements
1324 NSNumber *n = [availableElements lastObject];
1325 [availableElements removeLastObject];
1326 [elementMapping[direction] setObject: n forKey: rn];
1327 } else {
1328 // no available elements so add to end
1329 [elementMapping[direction] setObject: [NSNumber numberWithInt: totalElements] forKey: rn];
1330 ++totalElements;
1331 }
1332 printf("MOVE: rank = %d, direction = %d, add incoming element %d = %d\n", [modelController rank], direction, [[elementMapping[direction] objectForKey: rn] intValue], [rn intValue]);
1333 }
1334
1335 // now that we have the element numbers resolved
1336 // copy data for transfer down to GPU
1337 for (i = 0; i < neighborSectorSize * maxElementsPerSector; ++i) {
1338 if (neighborRanks[direction] == MPI_PROC_NULL) break;
1339
1340 int remoteElem = (*rElements)[direction][4*i];
1341 if (remoteElem == -1) break;
1342
1343 // search for element in map
1344 NSNumber *rn = [NSNumber numberWithInt: remoteElem];
1345 NSNumber *n = [elementMapping[direction] objectForKey: rn];
1346 if (!n) printf("MPI_ERROR: could not find matching element in map (%d, %f, %f, %f)\n", remoteElem, (*rElements)[direction][4*i+1], (*rElements)[direction][4*i+2], (*rElements)[direction][4*i+3]);
1347 else {
1348 // convert coordinates to local subspace
1349 coords[0] = (*rElements)[direction][4*i+1];
1350 coords[1] = (*rElements)[direction][4*i+2];
1351 coords[2] = (*rElements)[direction][4*i+3];
1352 if ([spatialModel isSubspace]) [spatialModel translateFloatCoordinatesToSubspace: coords forModel: self];
1353
1354 // copy data for GPU transfer
1355 int elemNum = [n intValue];
1356 (*iElements)[direction][currentPos[direction]] = (float)elemNum;
1357 (*iElements)[direction][currentPos[direction]+1] = coords[0];
1358 (*iElements)[direction][currentPos[direction]+2] = coords[1];
1359 (*iElements)[direction][currentPos[direction]+3] = coords[2];
1360 if ([insertElements[direction] containsObject: rn])
1361 printf("MOVE: rank = %d, direction = %d, incoming element %d, %.15f %.15f %.15f (%.15f %.15f %.15f)\n", [modelController rank], direction, elemNum, coords[0], coords[1], coords[2],
1362 (*rElements)[direction][4*i+1], (*rElements)[direction][4*i+2], (*rElements)[direction][4*i+3]);
1363
1364 currentPos[direction] += 4;
1365 }
1366 }
1367 //printf("rank = %d, direction = %d\n", [modelController rank], direction);
1368 }
1369
1370 // transfer data to GPU
1371 [self transferGPUMPIKernel: data toGPU: YES];
1372
1373 free(receiveElements);
1374 free(packElements);
1375
1376 if (TIME_PROFILE) bc_record_time([[NSString stringWithFormat: @"end MPI_MOVE rank = %d", [modelController rank]] UTF8String]);
1377 break;
1378 }
1379
1380 default:
1381 printf("ERROR: Unknown SubcellularElementModel MPI mode (%d)\n", aMode);
1382 }
1383
1384 }
1385
1386 @end