Changeset 1328
- Timestamp:
- 11/15/17 12:14:34 (7 years ago)
- Location:
- XIOS/dev/branch_openmp
- Files:
-
- 179 edited
Legend:
- Unmodified
- Added
- Removed
-
XIOS/dev/branch_openmp/bld.cfg
r1287 r1328 38 38 #bld::target test_new_features.exe 39 39 #bld::target test_unstruct_complete.exe 40 bld::target test_omp.exe 41 bld::target test_complete_omp.exe 42 bld::target test_remap_omp.exe 43 bld::target test_unstruct_omp.exe 40 #bld::target test_omp.exe 41 #bld::target test_complete_omp.exe 42 bld::target test_remap.exe 43 bld::target test_remap_ref.exe 44 #bld::target test_unstruct_omp.exe 44 45 #bld::target test_netcdf_omp.exe 45 #bld::target test_client.exe46 #bld::target test_complete.exe46 bld::target test_client.exe 47 bld::target test_complete.exe 47 48 #bld::target test_remap.exe 48 49 #bld::target test_unstruct.exe -
XIOS/dev/branch_openmp/extern/remap/src/clipper.cpp
r1205 r1328 4274 4274 { 4275 4275 //The equation of a line in general form (Ax + By + C = 0) 4276 //given 2 points (x ,y) & (x,y) is ...4277 //(y - y)x + (x - x)y + (y - y)x - (x - x)y= 04278 //A = (y - y); B = (x - x); C = (y - y)x - (x - x)y4279 //perpendicular distance of point (x ,y) = (Ax + By + C)/Sqrt(A + B)4276 //given 2 points (x¹,y¹) & (x²,y²) is ... 4277 //(y¹ - y²)x + (x² - x¹)y + (y² - y¹)x¹ - (x² - x¹)y¹ = 0 4278 //A = (y¹ - y²); B = (x² - x¹); C = (y² - y¹)x¹ - (x² - x¹)y¹ 4279 //perpendicular distance of point (x³,y³) = (Ax³ + By³ + C)/Sqrt(A² + B²) 4280 4280 //see http://en.wikipedia.org/wiki/Perpendicular_distance 4281 4281 double A = double(ln1.Y - ln2.Y); -
XIOS/dev/branch_openmp/extern/remap/src/cputime.cpp
r694 r1328 1 1 #include "mpi.hpp" 2 using namespace ep_lib; 2 3 3 4 namespace sphereRemap { -
XIOS/dev/branch_openmp/extern/remap/src/elt.hpp
r1172 r1328 53 53 struct Elt : Polyg 54 54 { 55 Elt() {} 56 Elt(const double *bounds_lon, const double *bounds_lat, int max_num_vert) 57 { 58 int k = 0; 59 vertex[k++] = xyz(bounds_lon[0], bounds_lat[0]); 60 for (int i = 1; i < max_num_vert; i++) 61 { 62 vertex[k] = xyz(bounds_lon[i], bounds_lat[i]); 63 /* netCDF convention: if first vertex repeats element is finished (at least three vertices == triagle) */ 64 if (k >= 3 && squaredist(vertex[k], vertex[0]) < EPS*EPS) 65 break; 66 /* eliminate zero edges: move to next vertex only if it is different */ 67 if (squaredist(vertex[k], vertex[k-1]) > EPS*EPS) 68 k++; 69 //else cout << "Removed edge " << k << " due to zero length (coinciding endpoints)." << endl ; 70 } 71 n = k; 72 x = barycentre(vertex, n); 73 } 55 Elt() {} 56 Elt(const double *bounds_lon, const double *bounds_lat, int max_num_vert) 57 { 58 int k = 0; 59 vertex[k++] = xyz(bounds_lon[0], bounds_lat[0]); 60 for (int i = 1; i < max_num_vert; i++) 61 { 62 vertex[k] = xyz(bounds_lon[i], bounds_lat[i]); 63 /* netCDF convention: if first vertex repeats element is finished (at least three vertices == triagle) */ 64 if (k >= 3 && squaredist(vertex[k], vertex[0]) < EPS*EPS) 65 break; 66 /* eliminate zero edges: move to next vertex only if it is different */ 67 if (squaredist(vertex[k], vertex[k-1]) > EPS*EPS) 68 k++; 69 else 70 /* cout << "Removed edge " << k << " due to zero length (coinciding endpoints)." << endl */ ; 71 } 72 n = k; 73 x = barycentre(vertex, n); 74 } 74 75 75 76 Elt& operator=(const Elt& rhs) … … 95 96 } 96 97 97 98 99 100 101 102 103 104 98 void delete_intersections() 99 { 100 for (list<Polyg*>::iterator it = this->is.begin(); it != this->is.end(); it++) 101 { 102 Polyg* poly = *it; 103 delete poly; 104 } 105 } 105 106 106 107 void insert_vertex(int i, const Coord& v) -
XIOS/dev/branch_openmp/extern/remap/src/gridRemap.cpp
r1155 r1328 11 11 12 12 CRemapGrid srcGrid; 13 #pragma omp threadprivate(srcGrid)14 15 13 CRemapGrid tgtGrid; 16 #pragma omp threadprivate(tgtGrid)17 14 18 15 } -
XIOS/dev/branch_openmp/extern/remap/src/gridRemap.hpp
r1220 r1328 14 14 Coord readPole(std::istream&); 15 15 16 16 extern CRemapGrid srcGrid; 17 extern CRemapGrid tgtGrid; 17 18 18 19 } -
XIOS/dev/branch_openmp/extern/remap/src/intersect.cpp
r1220 r1328 15 15 namespace sphereRemap { 16 16 17 extern CRemapGrid srcGrid;18 #pragma omp threadprivate(srcGrid)19 20 extern CRemapGrid tgtGrid;21 #pragma omp threadprivate(tgtGrid)22 23 17 using namespace std; 24 18 … … 27 21 int neighbour_idx(const Elt& a, const Elt& b) 28 22 { 29 30 31 32 33 34 35 36 37 38 39 40 41 42 23 for (int i = 0; i < a.n; i++) 24 { 25 for (int j = 0; j < b.n; j++) 26 { 27 assert(squaredist(a.vertex[ i ], b.vertex[ j ]) > EPS*EPS || 28 squaredist(a.vertex[(i+1)%a.n], b.vertex[(j+1)%b.n]) > EPS*EPS); 29 if ( squaredist(a.vertex[ i ], b.vertex[ j ]) < 1e-13*1e-13 && 30 squaredist(a.vertex[(i+1)%a.n], b.vertex[(j+b.n-1)%b.n]) < 1e-13*1e-13) 31 { 32 return i; 33 } 34 } 35 } 36 return NOT_FOUND; 43 37 } 44 38 … … 211 205 bool isNeighbour(Elt& a, const Elt& b) 212 206 { 213 207 // return neighbour_idx(a, b) != NOT_FOUND; 214 208 return insertNeighbour(a,b,false) ; 215 209 } … … 219 213 void intersect(Elt *a, Elt *b) 220 214 { 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 // 215 int na = a->n; /* vertices of a */ 216 int nb = b->n; /* vertices of b */ 217 Coord *c = new Coord[na+nb]; 218 Coord *c2 = new Coord[na+nb]; 219 Coord *xc = new Coord[na+nb]; 220 Coord *xc2 = new Coord[na+nb]; 221 Coord gc, gc2; 222 double *d = new double[na+nb]; 223 double *d2 = new double[na+nb]; 224 double are, are2; 225 Ipt ipt[NMAX*NMAX]; 226 Ipt ipt2[NMAX*NMAX]; 227 ptsec(a, b, ipt); 228 /* make ipt2 transpose of ipt */ 229 for (int ii = 0; ii < na; ii++) 230 for (int jj = 0; jj < nb; jj++) 231 ipt2[jj*na+ii] = ipt[ii*nb+jj]; 232 list<Sgm> iscot; 233 recense(a, b, ipt, iscot, 0); 234 recense(b, a, ipt2, iscot, 1); 235 236 int nseg = iscot.size(); 237 int nc = 0; 238 int nc2 = 0; 239 while (iscot.size() && nc < 2) 240 nc = assemble(iscot, c, d, xc); 241 while (iscot.size() && nc2 < 2) 242 nc2 = assemble(iscot, c2, d2, xc2); 243 // assert(nseg == nc + nc2 || nseg == 1); // unused segment 250 244 251 245 if (!(nseg == nc + nc2 || nseg == 1)) … … 266 260 267 261 // intersect_ym(a,b) ; 268 269 270 271 272 273 274 275 276 277 278 279 262 if (nc == 1) nc = 0; 263 if (nc2 == 1) nc2 = 0; 264 gc = barycentre(xc, nc); 265 gc2 = barycentre(xc2, nc2); 266 orient(nc, xc, c, d, gc); 267 268 Coord pole = srcGrid.pole; 269 if (pole == ORIGIN) pole = tgtGrid.pole; 270 const double MINBASE = 1e-11; 271 if (nc == 2) /* nc is the number of vertices of super mesh element */ 272 { 273 double base = arcdist(xc[0], xc[1]); 280 274 cerr << "DID ARRIVE " << base << xc[0] << xc[1] << endl; 281 282 283 284 285 286 287 288 289 290 291 275 gc = midpoint(gc, midpointSC(xc[0], xc[1])); 276 /* intersection area `are` must be zero here unless we have one great and one small circle */ 277 are = alun(base, fabs(scalarprod(xc[0], pole))); 278 } 279 else 280 { 281 are = airbar(nc, xc, c, d, pole, gc); 282 } 283 if (nc2 == 2) 284 { 285 double base = arcdist(xc2[0], xc2[1]); 292 286 cerr << "DID ARRIVE " << base << xc2[0] << xc2[1] << endl; 293 294 295 296 297 298 299 300 287 assert(base > MINBASE); 288 gc2 = midpoint(gc2, midpointSC(xc2[0], xc2[1])); 289 are2 = alun(base, fabs(scalarprod(xc2[0], pole))); // 0 290 } 291 else 292 { 293 are2 = airbar(nc2, xc2, c2, d2, pole, gc2); 294 } 301 295 302 296 // double ym_area=intersect_ym(a,b) ; 303 297 304 298 if (nc > 1) 305 306 307 308 309 310 311 312 313 314 299 { 300 /* create one super mesh polygon that src and dest point to */ 301 Polyg *is = new Polyg; 302 is->x = gc; 303 is->area = are; 304 is->id = b->id; 305 is->src_id = b->src_id; 306 is->n = nc; 307 (a->is).push_back(is); 308 (b->is).push_back(is); 315 309 /* 316 if (2*fabs(are-ym_area)/(are+ym_area) > 1.1 && ym_area>1e-8)310 if ( 2*fabs(are-ym_area)/(are+ym_area) > 1.1 && ym_area>1e-8) 317 311 { 318 312 cout<<"Big area difference : "<<are<<" "<<ym_area<<endl ; … … 320 314 } 321 315 */ 322 // 323 324 325 326 327 328 329 330 331 332 333 316 // cout<<"intersection : "<<are<<" "<< ym_area<<" diff : "<<fabs(are-ym_area)<<" ratio : "<<fabs(are-ym_area)/(0.5*(are+ym_area))<<endl ; 317 } 318 if (nc2 > 1) 319 { 320 Polyg *is = new Polyg; 321 is->x = gc2; 322 is->area = are2; 323 is->id = b->id; /* intersection holds id of corresponding source element (see Elt class definition for details about id) */ 324 is->src_id = b->src_id; 325 is->n = nc2; 326 (a->is).push_back(is); 327 (b->is).push_back(is); 334 328 /* 335 if ( 329 if ( 2*fabs(are-ym_area)/(are+ym_area) > 1.1 && ym_area>1e-8 ) 336 330 { 337 331 cout<<"Big area difference : "<<are<<" "<<ym_area<<endl ; … … 339 333 } 340 334 */ 341 // 342 335 // cout<<"intersection : "<<are2<<" "<< ym_area<<" diff : "<<fabs(are-ym_area)<<" ratio : "<<fabs(are-ym_area)/(0.5*(are+ym_area))<<endl ; 336 } 343 337 /* 344 338 if (nc<=1 && nc2<=1) … … 350 344 } 351 345 */ 352 353 354 355 356 357 358 } 359 360 } 346 delete [] c; 347 delete [] c2; 348 delete [] xc; 349 delete [] xc2; 350 delete [] d; 351 delete [] d2; 352 } 353 354 } -
XIOS/dev/branch_openmp/extern/remap/src/intersection_ym.cpp
r1220 r1328 13 13 14 14 namespace sphereRemap { 15 16 extern CRemapGrid srcGrid;17 #pragma omp threadprivate(srcGrid)18 19 extern CRemapGrid tgtGrid;20 #pragma omp threadprivate(tgtGrid)21 22 15 23 16 using namespace std; -
XIOS/dev/branch_openmp/extern/remap/src/libmapper.cpp
r1205 r1328 14 14 #include "mapper.hpp" 15 15 #include "cputime.hpp" // cputime 16 #include <stdio.h> 16 17 using namespace ep_lib; 17 18 18 19 using namespace sphereRemap ; … … 21 22 and deallocated during the second step (computing the weights) */ 22 23 Mapper *mapper; 23 #pragma omp threadprivate(mapper) 24 24 25 25 26 /** xxx_bounds_yyy is of length n_vert_per_cell_xxx*n_cell_xxx … … 33 34 int order, int* n_weights) 34 35 { 35 36 37 38 39 40 41 42 43 36 assert(src_bounds_lon); 37 assert(src_bounds_lat); 38 assert(n_vert_per_cell_src >= 3); 39 assert(n_cell_src >= 4); 40 assert(dst_bounds_lon); 41 assert(dst_bounds_lat); 42 assert(n_vert_per_cell_dst >= 3); 43 assert(n_cell_dst >= 4); 44 assert(1 <= order && order <= 2); 44 45 45 46 mapper = new Mapper(MPI_COMM_WORLD); … … 80 81 double tic = cputime(); 81 82 mapper = new Mapper(MPI_COMM_WORLD); 82 83 mapper->setVerbosity(PROGRESS) ; 83 84 mapper->buildSSTree(src_msh, dst_msh); 84 85 double tac = cputime(); … … 149 150 char **argv = NULL; 150 151 MPI_Init(&argc, &argv);*/ 151 //MPI_Init(NULL, NULL); 152 int provided; 153 MPI_Init_thread(NULL, NULL, 3, &provided); 154 assert(provided >= 3); 152 MPI_Init(NULL, NULL); 155 153 } 156 154 -
XIOS/dev/branch_openmp/extern/remap/src/mapper.cpp
r1220 r1328 12 12 13 13 #include "mapper.hpp" 14 using namespace ep_lib; 14 15 15 16 namespace sphereRemap { 16 17 extern CRemapGrid srcGrid;18 #pragma omp threadprivate(srcGrid)19 20 extern CRemapGrid tgtGrid;21 #pragma omp threadprivate(tgtGrid)22 23 17 24 18 /* A subdivition of an array into N sub-arays … … 72 66 void Mapper::setTargetMesh(const double* boundsLon, const double* boundsLat, int nVertex, int nbCells, const double* pole, const long int* globalId) 73 67 { 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 68 tgtGrid.pole = Coord(pole[0], pole[1], pole[2]); 69 70 int mpiRank, mpiSize; 71 MPI_Comm_rank(communicator, &mpiRank); 72 MPI_Comm_size(communicator, &mpiSize); 73 74 targetElements.reserve(nbCells); 75 targetMesh.reserve(nbCells); 76 77 targetGlobalId.resize(nbCells) ; 78 if (globalId==NULL) 79 { 80 long int offset ; 81 long int nb=nbCells ; 82 MPI_Scan(&nb,&offset,1,MPI_LONG,MPI_SUM,communicator) ; 83 offset=offset-nb ; 84 for(int i=0;i<nbCells;i++) targetGlobalId[i]=offset+i ; 85 } 86 else targetGlobalId.assign(globalId,globalId+nbCells); 87 88 for (int i = 0; i < nbCells; i++) 89 { 90 int offs = i*nVertex; 91 Elt elt(boundsLon + offs, boundsLat + offs, nVertex); 92 targetElements.push_back(elt); 93 targetMesh.push_back(Node(elt.x, cptRadius(elt), &sourceElements.back())); 94 cptEltGeom(targetElements[i], Coord(pole[0], pole[1], pole[2])); 95 } 102 96 103 97 … … 106 100 void Mapper::setSourceValue(const double* val) 107 101 { 108 109 102 int size=sourceElements.size() ; 103 for(int i=0;i<size;++i) sourceElements[i].val=val[i] ; 110 104 } 111 105 112 106 void Mapper::getTargetValue(double* val) 113 107 { 114 115 108 int size=targetElements.size() ; 109 for(int i=0;i<size;++i) val[i]=targetElements[i].val ; 116 110 } 117 111 118 112 vector<double> Mapper::computeWeights(int interpOrder, bool renormalize, bool quantity) 119 113 { 120 vector<double> timings; 121 int mpiSize, mpiRank; 122 MPI_Comm_size(communicator, &mpiSize); 123 MPI_Comm_rank(communicator, &mpiRank); 124 125 this->buildSSTree(sourceMesh, targetMesh); 126 127 if (mpiRank == 0 && verbose) cout << "Computing intersections ..." << endl; 128 double tic = cputime(); 129 computeIntersection(&targetElements[0], targetElements.size()); 130 timings.push_back(cputime() - tic); 131 132 tic = cputime(); 133 if (interpOrder == 2) { 134 if (mpiRank == 0 && verbose) cout << "Computing grads ..." << endl; 135 buildMeshTopology(); 136 computeGrads(); 137 } 138 timings.push_back(cputime() - tic); 139 140 /* Prepare computation of weights */ 141 /* compute number of intersections which for the first order case 142 corresponds to the number of edges in the remap matrix */ 143 int nIntersections = 0; 144 for (int j = 0; j < targetElements.size(); j++) 145 { 146 Elt &elt = targetElements[j]; 147 for (list<Polyg*>::iterator it = elt.is.begin(); it != elt.is.end(); it++) 148 nIntersections++; 149 } 150 /* overallocate for NMAX neighbours for each elements */ 151 remapMatrix = new double[nIntersections*NMAX]; 152 srcAddress = new int[nIntersections*NMAX]; 153 srcRank = new int[nIntersections*NMAX]; 154 dstAddress = new int[nIntersections*NMAX]; 155 sourceWeightId =new long[nIntersections*NMAX]; 156 targetWeightId =new long[nIntersections*NMAX]; 157 158 159 if (mpiRank == 0 && verbose) cout << "Remapping..." << endl; 160 tic = cputime(); 161 nWeights = remap(&targetElements[0], targetElements.size(), interpOrder, renormalize, quantity); 162 timings.push_back(cputime() - tic); 163 164 for (int i = 0; i < targetElements.size(); i++) targetElements[i].delete_intersections(); 165 166 return timings; 114 vector<double> timings; 115 int mpiSize, mpiRank; 116 MPI_Comm_size(communicator, &mpiSize); 117 MPI_Comm_rank(communicator, &mpiRank); 118 119 120 this->buildSSTree(sourceMesh, targetMesh); 121 122 if (mpiRank == 0 && verbose) cout << "Computing intersections ..." << endl; 123 double tic = cputime(); 124 computeIntersection(&targetElements[0], targetElements.size()); 125 timings.push_back(cputime() - tic); 126 127 tic = cputime(); 128 if (interpOrder == 2) { 129 if (mpiRank == 0 && verbose) cout << "Computing grads ..." << endl; 130 buildMeshTopology(); 131 computeGrads(); 132 } 133 timings.push_back(cputime() - tic); 134 135 /* Prepare computation of weights */ 136 /* compute number of intersections which for the first order case 137 corresponds to the number of edges in the remap matrix */ 138 int nIntersections = 0; 139 for (int j = 0; j < targetElements.size(); j++) 140 { 141 Elt &elt = targetElements[j]; 142 for (list<Polyg*>::iterator it = elt.is.begin(); it != elt.is.end(); it++) 143 nIntersections++; 144 } 145 /* overallocate for NMAX neighbours for each elements */ 146 remapMatrix = new double[nIntersections*NMAX]; 147 srcAddress = new int[nIntersections*NMAX]; 148 srcRank = new int[nIntersections*NMAX]; 149 dstAddress = new int[nIntersections*NMAX]; 150 sourceWeightId =new long[nIntersections*NMAX]; 151 targetWeightId =new long[nIntersections*NMAX]; 152 153 154 if (mpiRank == 0 && verbose) cout << "Remapping..." << endl; 155 tic = cputime(); 156 nWeights = remap(&targetElements[0], targetElements.size(), interpOrder, renormalize, quantity); 157 timings.push_back(cputime() - tic); 158 159 for (int i = 0; i < targetElements.size(); i++) targetElements[i].delete_intersections(); 160 161 return timings; 167 162 } 168 163 169 164 /** 170 @param elements are cells of the target grid that are distributed over CPUs171 indepentently of the distribution of the SS-tree.172 @param nbElements is the size of the elements array.173 @param order is the order of interpolaton (must be 1 or 2).174 165 @param elements are cells of the target grid that are distributed over CPUs 166 indepentently of the distribution of the SS-tree. 167 @param nbElements is the size of the elements array. 168 @param order is the order of interpolaton (must be 1 or 2). 169 */ 175 170 int Mapper::remap(Elt *elements, int nbElements, int order, bool renormalize, bool quantity) 176 171 { 177 int mpiSize, mpiRank; 178 MPI_Comm_size(communicator, &mpiSize); 179 MPI_Comm_rank(communicator, &mpiRank); 180 181 /* create list of intersections (super mesh elements) for each rank */ 182 multimap<int, Polyg *> *elementList = new multimap<int, Polyg *>[mpiSize]; 183 for (int j = 0; j < nbElements; j++) 184 { 185 Elt& e = elements[j]; 186 for (list<Polyg *>::iterator it = e.is.begin(); it != e.is.end(); it++) 187 elementList[(*it)->id.rank].insert(pair<int, Polyg *>((*it)->id.ind, *it)); 188 } 189 190 int *nbSendElement = new int[mpiSize]; 191 int **sendElement = new int*[mpiSize]; /* indices of elements required from other rank */ 192 double **recvValue = new double*[mpiSize]; 193 double **recvArea = new double*[mpiSize]; 194 Coord **recvGrad = new Coord*[mpiSize]; 195 GloId **recvNeighIds = new GloId*[mpiSize]; /* ids of the of the source neighbours which also contribute through gradient */ 196 for (int rank = 0; rank < mpiSize; rank++) 197 { 198 /* get size for allocation */ 199 int last = -1; /* compares unequal to any index */ 200 int index = -1; /* increased to starting index 0 in first iteration */ 201 for (multimap<int, Polyg *>::iterator it = elementList[rank].begin(); it != elementList[rank].end(); ++it) 202 { 203 if (last != it->first) 204 index++; 205 (it->second)->id.ind = index; 206 last = it->first; 207 } 208 nbSendElement[rank] = index + 1; 209 210 /* if size is non-zero allocate and collect indices of elements on other ranks that we intersect */ 211 if (nbSendElement[rank] > 0) 212 { 213 sendElement[rank] = new int[nbSendElement[rank]]; 214 recvValue[rank] = new double[nbSendElement[rank]]; 215 recvArea[rank] = new double[nbSendElement[rank]]; 216 if (order == 2) 217 { 218 recvNeighIds[rank] = new GloId[nbSendElement[rank]*(NMAX+1)]; 219 recvGrad[rank] = new Coord[nbSendElement[rank]*(NMAX+1)]; 220 } 221 else 222 recvNeighIds[rank] = new GloId[nbSendElement[rank]]; 223 224 last = -1; 225 index = -1; 226 for (multimap<int, Polyg *>::iterator it = elementList[rank].begin(); it != elementList[rank].end(); ++it) 227 { 228 if (last != it->first) 229 index++; 230 sendElement[rank][index] = it->first; 231 last = it->first; 232 } 233 } 234 } 235 236 /* communicate sizes of source elements to be sent (index lists and later values and gradients) */ 237 int *nbRecvElement = new int[mpiSize]; 238 MPI_Alltoall(nbSendElement, 1, MPI_INT, nbRecvElement, 1, MPI_INT, communicator); 239 240 /* communicate indices of source elements on other ranks whoes value and gradient we need (since intersection) */ 241 int nbSendRequest = 0; 242 int nbRecvRequest = 0; 243 int **recvElement = new int*[mpiSize]; 244 double **sendValue = new double*[mpiSize]; 245 double **sendArea = new double*[mpiSize]; 246 Coord **sendGrad = new Coord*[mpiSize]; 247 GloId **sendNeighIds = new GloId*[mpiSize]; 248 MPI_Request *sendRequest = new MPI_Request[4*mpiSize]; 249 MPI_Request *recvRequest = new MPI_Request[4*mpiSize]; 250 for (int rank = 0; rank < mpiSize; rank++) 251 { 252 if (nbSendElement[rank] > 0) 253 { 254 MPI_Issend(sendElement[rank], nbSendElement[rank], MPI_INT, rank, 0, communicator, &sendRequest[nbSendRequest]); 255 nbSendRequest++; 256 } 257 258 if (nbRecvElement[rank] > 0) 259 { 260 recvElement[rank] = new int[nbRecvElement[rank]]; 261 sendValue[rank] = new double[nbRecvElement[rank]]; 262 sendArea[rank] = new double[nbRecvElement[rank]]; 263 if (order == 2) 264 { 265 sendNeighIds[rank] = new GloId[nbRecvElement[rank]*(NMAX+1)]; 266 sendGrad[rank] = new Coord[nbRecvElement[rank]*(NMAX+1)]; 267 } 268 else 269 { 270 sendNeighIds[rank] = new GloId[nbRecvElement[rank]]; 271 } 272 MPI_Irecv(recvElement[rank], nbRecvElement[rank], MPI_INT, rank, 0, communicator, &recvRequest[nbRecvRequest]); 273 nbRecvRequest++; 274 } 275 } 276 277 MPI_Status *status = new MPI_Status[4*mpiSize]; 278 279 MPI_Waitall(nbSendRequest, sendRequest, status); 280 MPI_Waitall(nbRecvRequest, recvRequest, status); 281 282 /* for all indices that have been received from requesting ranks: pack values and gradients, then send */ 283 nbSendRequest = 0; 284 nbRecvRequest = 0; 285 for (int rank = 0; rank < mpiSize; rank++) 286 { 287 if (nbRecvElement[rank] > 0) 288 { 289 int jj = 0; // jj == j if no weight writing 290 for (int j = 0; j < nbRecvElement[rank]; j++) 291 { 292 sendValue[rank][j] = sstree.localElements[recvElement[rank][j]].val; 293 sendArea[rank][j] = sstree.localElements[recvElement[rank][j]].area; 294 if (order == 2) 295 { 296 sendGrad[rank][jj] = sstree.localElements[recvElement[rank][j]].grad; 297 sendNeighIds[rank][jj] = sstree.localElements[recvElement[rank][j]].src_id; 298 jj++; 299 for (int i = 0; i < NMAX; i++) 300 { 301 sendGrad[rank][jj] = sstree.localElements[recvElement[rank][j]].gradNeigh[i]; 302 sendNeighIds[rank][jj] = sstree.localElements[recvElement[rank][j]].neighId[i]; 303 jj++; 304 } 305 } 306 else 307 sendNeighIds[rank][j] = sstree.localElements[recvElement[rank][j]].src_id; 308 } 309 MPI_Issend(sendValue[rank], nbRecvElement[rank], MPI_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]); 310 nbSendRequest++; 311 MPI_Issend(sendArea[rank], nbRecvElement[rank], MPI_DOUBLE, rank, 1, communicator, &sendRequest[nbSendRequest]); 312 nbSendRequest++; 313 if (order == 2) 314 { 315 MPI_Issend(sendGrad[rank], 3*nbRecvElement[rank]*(NMAX+1), 316 MPI_DOUBLE, rank, 2, communicator, &sendRequest[nbSendRequest]); 317 nbSendRequest++; 318 MPI_Issend(sendNeighIds[rank], 4*nbRecvElement[rank]*(NMAX+1), MPI_INT, rank, 3, communicator, &sendRequest[nbSendRequest]); 319 //ym --> attention taille GloId 320 nbSendRequest++; 321 } 322 else 323 { 324 MPI_Issend(sendNeighIds[rank], 4*nbRecvElement[rank], MPI_INT, rank, 4, communicator, &sendRequest[nbSendRequest]); 325 //ym --> attention taille GloId 326 nbSendRequest++; 327 } 328 } 329 if (nbSendElement[rank] > 0) 330 { 331 MPI_Irecv(recvValue[rank], nbSendElement[rank], MPI_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]); 332 nbRecvRequest++; 333 MPI_Irecv(recvArea[rank], nbSendElement[rank], MPI_DOUBLE, rank, 1, communicator, &recvRequest[nbRecvRequest]); 334 nbRecvRequest++; 335 if (order == 2) 336 { 337 MPI_Irecv(recvGrad[rank], 3*nbSendElement[rank]*(NMAX+1), 338 MPI_DOUBLE, rank, 2, communicator, &recvRequest[nbRecvRequest]); 339 nbRecvRequest++; 340 MPI_Irecv(recvNeighIds[rank], 4*nbSendElement[rank]*(NMAX+1), MPI_INT, rank, 3, communicator, &recvRequest[nbRecvRequest]); 341 //ym --> attention taille GloId 342 nbRecvRequest++; 343 } 344 else 345 { 346 MPI_Irecv(recvNeighIds[rank], 4*nbSendElement[rank], MPI_INT, rank, 4, communicator, &recvRequest[nbRecvRequest]); 347 //ym --> attention taille GloId 348 nbRecvRequest++; 349 } 350 } 351 } 352 353 MPI_Waitall(nbRecvRequest, recvRequest, status); 354 MPI_Waitall(nbSendRequest, sendRequest, status); 355 356 357 358 /* now that all values and gradients are available use them to computed interpolated values on target 359 and also to compute weights */ 360 int i = 0; 361 for (int j = 0; j < nbElements; j++) 362 { 363 Elt& e = elements[j]; 364 365 /* since for the 2nd order case source grid elements can contribute to a destination grid element over several "paths" 366 (step1: gradient is computed using neighbours on same grid, step2: intersection uses several elements on other grid) 367 accumulate them so that there is only one final weight between two elements */ 368 map<GloId,double> wgt_map; 369 370 /* for destination element `e` loop over all intersetions/the corresponding source elements */ 371 for (list<Polyg *>::iterator it = e.is.begin(); it != e.is.end(); it++) 372 { 373 /* it is the intersection element, so it->x and it->area are barycentre and area of intersection element (super mesh) 374 but it->id is id of the source element that it intersects */ 375 int n1 = (*it)->id.ind; 376 int rank = (*it)->id.rank; 377 double fk = recvValue[rank][n1]; 378 double srcArea = recvArea[rank][n1]; 379 double w = (*it)->area; 380 if (quantity) w/=srcArea ; 381 382 /* first order: src value times weight (weight = supermesh area), later divide by target area */ 383 int kk = (order == 2) ? n1 * (NMAX + 1) : n1; 384 GloId neighID = recvNeighIds[rank][kk]; 385 wgt_map[neighID] += w; 386 387 if (order == 2) 388 { 389 for (int k = 0; k < NMAX+1; k++) 390 { 391 int kk = n1 * (NMAX + 1) + k; 392 GloId neighID = recvNeighIds[rank][kk]; 393 if (neighID.ind != -1) wgt_map[neighID] += w * scalarprod(recvGrad[rank][kk], (*it)->x); 394 } 395 396 } 397 } 398 399 double renorm=0; 400 if (renormalize) 401 for (map<GloId,double>::iterator it = wgt_map.begin(); it != wgt_map.end(); it++) renorm+=it->second / e.area; 402 else renorm=1. ; 403 404 for (map<GloId,double>::iterator it = wgt_map.begin(); it != wgt_map.end(); it++) 405 { 406 if (quantity) this->remapMatrix[i] = (it->second ) / renorm; 407 else this->remapMatrix[i] = (it->second / e.area) / renorm; 408 this->srcAddress[i] = it->first.ind; 409 this->srcRank[i] = it->first.rank; 410 this->dstAddress[i] = j; 411 this->sourceWeightId[i]= it->first.globalId ; 412 this->targetWeightId[i]= targetGlobalId[j] ; 413 i++; 414 } 415 } 416 417 /* free all memory allocated in this function */ 418 for (int rank = 0; rank < mpiSize; rank++) 419 { 420 if (nbSendElement[rank] > 0) 421 { 422 delete[] sendElement[rank]; 423 delete[] recvValue[rank]; 424 delete[] recvArea[rank]; 425 if (order == 2) 426 { 427 delete[] recvGrad[rank]; 428 } 429 delete[] recvNeighIds[rank]; 430 } 431 if (nbRecvElement[rank] > 0) 432 { 433 delete[] recvElement[rank]; 434 delete[] sendValue[rank]; 435 delete[] sendArea[rank]; 436 if (order == 2) 437 delete[] sendGrad[rank]; 438 delete[] sendNeighIds[rank]; 439 } 440 } 441 delete[] status; 442 delete[] sendRequest; 443 delete[] recvRequest; 444 delete[] elementList; 445 delete[] nbSendElement; 446 delete[] nbRecvElement; 447 delete[] sendElement; 448 delete[] recvElement; 449 delete[] sendValue; 450 delete[] recvValue; 451 delete[] sendGrad; 452 delete[] recvGrad; 453 delete[] sendNeighIds; 454 delete[] recvNeighIds; 455 456 return i; 172 int mpiSize, mpiRank; 173 MPI_Comm_size(communicator, &mpiSize); 174 MPI_Comm_rank(communicator, &mpiRank); 175 176 /* create list of intersections (super mesh elements) for each rank */ 177 multimap<int, Polyg *> *elementList = new multimap<int, Polyg *>[mpiSize]; 178 for (int j = 0; j < nbElements; j++) 179 { 180 Elt& e = elements[j]; 181 for (list<Polyg *>::iterator it = e.is.begin(); it != e.is.end(); it++) 182 elementList[(*it)->id.rank].insert(pair<int, Polyg *>((*it)->id.ind, *it)); 183 } 184 185 int *nbSendElement = new int[mpiSize]; 186 int **sendElement = new int*[mpiSize]; /* indices of elements required from other rank */ 187 double **recvValue = new double*[mpiSize]; 188 double **recvArea = new double*[mpiSize]; 189 Coord **recvGrad = new Coord*[mpiSize]; 190 GloId **recvNeighIds = new GloId*[mpiSize]; /* ids of the of the source neighbours which also contribute through gradient */ 191 for (int rank = 0; rank < mpiSize; rank++) 192 { 193 /* get size for allocation */ 194 int last = -1; /* compares unequal to any index */ 195 int index = -1; /* increased to starting index 0 in first iteration */ 196 for (multimap<int, Polyg *>::iterator it = elementList[rank].begin(); it != elementList[rank].end(); ++it) 197 { 198 if (last != it->first) 199 index++; 200 (it->second)->id.ind = index; 201 last = it->first; 202 } 203 nbSendElement[rank] = index + 1; 204 205 /* if size is non-zero allocate and collect indices of elements on other ranks that we intersect */ 206 if (nbSendElement[rank] > 0) 207 { 208 sendElement[rank] = new int[nbSendElement[rank]]; 209 recvValue[rank] = new double[nbSendElement[rank]]; 210 recvArea[rank] = new double[nbSendElement[rank]]; 211 if (order == 2) 212 { 213 recvNeighIds[rank] = new GloId[nbSendElement[rank]*(NMAX+1)]; 214 recvGrad[rank] = new Coord[nbSendElement[rank]*(NMAX+1)]; 215 } 216 else 217 recvNeighIds[rank] = new GloId[nbSendElement[rank]]; 218 219 last = -1; 220 index = -1; 221 for (multimap<int, Polyg *>::iterator it = elementList[rank].begin(); it != elementList[rank].end(); ++it) 222 { 223 if (last != it->first) 224 index++; 225 sendElement[rank][index] = it->first; 226 last = it->first; 227 } 228 } 229 } 230 231 /* communicate sizes of source elements to be sent (index lists and later values and gradients) */ 232 int *nbRecvElement = new int[mpiSize]; 233 MPI_Alltoall(nbSendElement, 1, MPI_INT, nbRecvElement, 1, MPI_INT, communicator); 234 235 /* communicate indices of source elements on other ranks whoes value and gradient we need (since intersection) */ 236 int nbSendRequest = 0; 237 int nbRecvRequest = 0; 238 int **recvElement = new int*[mpiSize]; 239 double **sendValue = new double*[mpiSize]; 240 double **sendArea = new double*[mpiSize]; 241 Coord **sendGrad = new Coord*[mpiSize]; 242 GloId **sendNeighIds = new GloId*[mpiSize]; 243 MPI_Request *sendRequest = new MPI_Request[4*mpiSize]; 244 MPI_Request *recvRequest = new MPI_Request[4*mpiSize]; 245 for (int rank = 0; rank < mpiSize; rank++) 246 { 247 if (nbSendElement[rank] > 0) 248 { 249 MPI_Issend(sendElement[rank], nbSendElement[rank], MPI_INT, rank, 0, communicator, &sendRequest[nbSendRequest]); 250 nbSendRequest++; 251 } 252 253 if (nbRecvElement[rank] > 0) 254 { 255 recvElement[rank] = new int[nbRecvElement[rank]]; 256 sendValue[rank] = new double[nbRecvElement[rank]]; 257 sendArea[rank] = new double[nbRecvElement[rank]]; 258 if (order == 2) 259 { 260 sendNeighIds[rank] = new GloId[nbRecvElement[rank]*(NMAX+1)]; 261 sendGrad[rank] = new Coord[nbRecvElement[rank]*(NMAX+1)]; 262 } 263 else 264 { 265 sendNeighIds[rank] = new GloId[nbRecvElement[rank]]; 266 } 267 MPI_Irecv(recvElement[rank], nbRecvElement[rank], MPI_INT, rank, 0, communicator, &recvRequest[nbRecvRequest]); 268 nbRecvRequest++; 269 } 270 } 271 272 MPI_Status *status = new MPI_Status[4*mpiSize]; 273 274 MPI_Waitall(nbSendRequest, sendRequest, status); 275 MPI_Waitall(nbRecvRequest, recvRequest, status); 276 277 /* for all indices that have been received from requesting ranks: pack values and gradients, then send */ 278 nbSendRequest = 0; 279 nbRecvRequest = 0; 280 for (int rank = 0; rank < mpiSize; rank++) 281 { 282 if (nbRecvElement[rank] > 0) 283 { 284 int jj = 0; // jj == j if no weight writing 285 for (int j = 0; j < nbRecvElement[rank]; j++) 286 { 287 sendValue[rank][j] = sstree.localElements[recvElement[rank][j]].val; 288 sendArea[rank][j] = sstree.localElements[recvElement[rank][j]].area; 289 if (order == 2) 290 { 291 sendGrad[rank][jj] = sstree.localElements[recvElement[rank][j]].grad; 292 sendNeighIds[rank][jj] = sstree.localElements[recvElement[rank][j]].src_id; 293 jj++; 294 for (int i = 0; i < NMAX; i++) 295 { 296 sendGrad[rank][jj] = sstree.localElements[recvElement[rank][j]].gradNeigh[i]; 297 sendNeighIds[rank][jj] = sstree.localElements[recvElement[rank][j]].neighId[i]; 298 jj++; 299 } 300 } 301 else 302 sendNeighIds[rank][j] = sstree.localElements[recvElement[rank][j]].src_id; 303 } 304 MPI_Issend(sendValue[rank], nbRecvElement[rank], MPI_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]); 305 nbSendRequest++; 306 MPI_Issend(sendArea[rank], nbRecvElement[rank], MPI_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]); 307 nbSendRequest++; 308 if (order == 2) 309 { 310 MPI_Issend(sendGrad[rank], 3*nbRecvElement[rank]*(NMAX+1), 311 MPI_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]); 312 nbSendRequest++; 313 MPI_Issend(sendNeighIds[rank], 4*nbRecvElement[rank]*(NMAX+1), MPI_INT, rank, 0, communicator, &sendRequest[nbSendRequest]); 314 //ym --> attention taille GloId 315 nbSendRequest++; 316 } 317 else 318 { 319 MPI_Issend(sendNeighIds[rank], 4*nbRecvElement[rank], MPI_INT, rank, 0, communicator, &sendRequest[nbSendRequest]); 320 //ym --> attention taille GloId 321 nbSendRequest++; 322 } 323 } 324 if (nbSendElement[rank] > 0) 325 { 326 MPI_Irecv(recvValue[rank], nbSendElement[rank], MPI_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]); 327 nbRecvRequest++; 328 MPI_Irecv(recvArea[rank], nbSendElement[rank], MPI_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]); 329 nbRecvRequest++; 330 if (order == 2) 331 { 332 MPI_Irecv(recvGrad[rank], 3*nbSendElement[rank]*(NMAX+1), 333 MPI_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]); 334 nbRecvRequest++; 335 MPI_Irecv(recvNeighIds[rank], 4*nbSendElement[rank]*(NMAX+1), MPI_INT, rank, 0, communicator, &recvRequest[nbRecvRequest]); 336 //ym --> attention taille GloId 337 nbRecvRequest++; 338 } 339 else 340 { 341 MPI_Irecv(recvNeighIds[rank], 4*nbSendElement[rank], MPI_INT, rank, 0, communicator, &recvRequest[nbRecvRequest]); 342 //ym --> attention taille GloId 343 nbRecvRequest++; 344 } 345 } 346 } 347 348 MPI_Waitall(nbSendRequest, sendRequest, status); 349 MPI_Waitall(nbRecvRequest, recvRequest, status); 350 351 352 /* now that all values and gradients are available use them to computed interpolated values on target 353 and also to compute weights */ 354 int i = 0; 355 for (int j = 0; j < nbElements; j++) 356 { 357 Elt& e = elements[j]; 358 359 /* since for the 2nd order case source grid elements can contribute to a destination grid element over several "paths" 360 (step1: gradient is computed using neighbours on same grid, step2: intersection uses several elements on other grid) 361 accumulate them so that there is only one final weight between two elements */ 362 map<GloId,double> wgt_map; 363 364 /* for destination element `e` loop over all intersetions/the corresponding source elements */ 365 for (list<Polyg *>::iterator it = e.is.begin(); it != e.is.end(); it++) 366 { 367 /* it is the intersection element, so it->x and it->area are barycentre and area of intersection element (super mesh) 368 but it->id is id of the source element that it intersects */ 369 int n1 = (*it)->id.ind; 370 int rank = (*it)->id.rank; 371 double fk = recvValue[rank][n1]; 372 double srcArea = recvArea[rank][n1]; 373 double w = (*it)->area; 374 if (quantity) w/=srcArea ; 375 376 /* first order: src value times weight (weight = supermesh area), later divide by target area */ 377 int kk = (order == 2) ? n1 * (NMAX + 1) : n1; 378 GloId neighID = recvNeighIds[rank][kk]; 379 wgt_map[neighID] += w; 380 381 if (order == 2) 382 { 383 for (int k = 0; k < NMAX+1; k++) 384 { 385 int kk = n1 * (NMAX + 1) + k; 386 GloId neighID = recvNeighIds[rank][kk]; 387 if (neighID.ind != -1) wgt_map[neighID] += w * scalarprod(recvGrad[rank][kk], (*it)->x); 388 } 389 390 } 391 } 392 393 double renorm=0; 394 if (renormalize) 395 for (map<GloId,double>::iterator it = wgt_map.begin(); it != wgt_map.end(); it++) renorm+=it->second / e.area; 396 else renorm=1. ; 397 398 for (map<GloId,double>::iterator it = wgt_map.begin(); it != wgt_map.end(); it++) 399 { 400 if (quantity) this->remapMatrix[i] = (it->second ) / renorm; 401 else this->remapMatrix[i] = (it->second / e.area) / renorm; 402 this->srcAddress[i] = it->first.ind; 403 this->srcRank[i] = it->first.rank; 404 this->dstAddress[i] = j; 405 this->sourceWeightId[i]= it->first.globalId ; 406 this->targetWeightId[i]= targetGlobalId[j] ; 407 i++; 408 } 409 } 410 411 /* free all memory allocated in this function */ 412 for (int rank = 0; rank < mpiSize; rank++) 413 { 414 if (nbSendElement[rank] > 0) 415 { 416 delete[] sendElement[rank]; 417 delete[] recvValue[rank]; 418 delete[] recvArea[rank]; 419 if (order == 2) 420 { 421 delete[] recvGrad[rank]; 422 } 423 delete[] recvNeighIds[rank]; 424 } 425 if (nbRecvElement[rank] > 0) 426 { 427 delete[] recvElement[rank]; 428 delete[] sendValue[rank]; 429 delete[] sendArea[rank]; 430 if (order == 2) 431 delete[] sendGrad[rank]; 432 delete[] sendNeighIds[rank]; 433 } 434 } 435 delete[] status; 436 delete[] sendRequest; 437 delete[] recvRequest; 438 delete[] elementList; 439 delete[] nbSendElement; 440 delete[] nbRecvElement; 441 delete[] sendElement; 442 delete[] recvElement; 443 delete[] sendValue; 444 delete[] recvValue; 445 delete[] sendGrad; 446 delete[] recvGrad; 447 delete[] sendNeighIds; 448 delete[] recvNeighIds; 449 return i; 457 450 } 458 451 459 452 void Mapper::computeGrads() 460 453 { 461 462 463 464 465 466 467 468 469 470 454 /* array of pointers to collect local elements and elements received from other cpu */ 455 vector<Elt*> globalElements(sstree.nbLocalElements + nbNeighbourElements); 456 int index = 0; 457 for (int i = 0; i < sstree.nbLocalElements; i++, index++) 458 globalElements[index] = &(sstree.localElements[i]); 459 for (int i = 0; i < nbNeighbourElements; i++, index++) 460 globalElements[index] = &neighbourElements[i]; 461 462 update_baryc(sstree.localElements, sstree.nbLocalElements); 463 computeGradients(&globalElements[0], sstree.nbLocalElements); 471 464 } 472 465 473 466 /** for each element of the source grid, finds all the neighbouring elements that share an edge 474 (filling array neighbourElements). This is used later to compute gradients */467 (filling array neighbourElements). This is used later to compute gradients */ 475 468 void Mapper::buildMeshTopology() 476 469 { 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 intersector of a local node */673 674 675 676 use the SS-tree to find all elements (including neighbourElements)677 who are potential neighbours because their circles intersect,678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 whoes circles intersect with the circle around `elt` (the SS intersectors)693 and check if they are neighbours in the sense that the two elements share an edge.694 If they do, save this information for elt */695 696 697 698 699 700 701 470 int mpiSize, mpiRank; 471 MPI_Comm_size(communicator, &mpiSize); 472 MPI_Comm_rank(communicator, &mpiRank); 473 474 vector<Node> *routingList = new vector<Node>[mpiSize]; 475 vector<vector<int> > routes(sstree.localTree.leafs.size()); 476 477 sstree.routeIntersections(routes, sstree.localTree.leafs); 478 479 for (int i = 0; i < routes.size(); ++i) 480 for (int k = 0; k < routes[i].size(); ++k) 481 routingList[routes[i][k]].push_back(sstree.localTree.leafs[i]); 482 routingList[mpiRank].clear(); 483 484 485 CMPIRouting mpiRoute(communicator); 486 mpiRoute.init(routes); 487 int nRecv = mpiRoute.getTotalSourceElement(); 488 489 int *nbSendNode = new int[mpiSize]; 490 int *nbRecvNode = new int[mpiSize]; 491 int *sendMessageSize = new int[mpiSize]; 492 int *recvMessageSize = new int[mpiSize]; 493 494 for (int rank = 0; rank < mpiSize; rank++) 495 { 496 nbSendNode[rank] = routingList[rank].size(); 497 sendMessageSize[rank] = 0; 498 for (size_t j = 0; j < routingList[rank].size(); j++) 499 { 500 Elt *elt = (Elt *) (routingList[rank][j].data); 501 sendMessageSize[rank] += packedPolygonSize(*elt); 502 } 503 } 504 505 MPI_Alltoall(nbSendNode, 1, MPI_INT, nbRecvNode, 1, MPI_INT, communicator); 506 MPI_Alltoall(sendMessageSize, 1, MPI_INT, recvMessageSize, 1, MPI_INT, communicator); 507 508 char **sendBuffer = new char*[mpiSize]; 509 char **recvBuffer = new char*[mpiSize]; 510 int *pos = new int[mpiSize]; 511 512 for (int rank = 0; rank < mpiSize; rank++) 513 { 514 if (nbSendNode[rank] > 0) sendBuffer[rank] = new char[sendMessageSize[rank]]; 515 if (nbRecvNode[rank] > 0) recvBuffer[rank] = new char[recvMessageSize[rank]]; 516 } 517 518 for (int rank = 0; rank < mpiSize; rank++) 519 { 520 pos[rank] = 0; 521 for (size_t j = 0; j < routingList[rank].size(); j++) 522 { 523 Elt *elt = (Elt *) (routingList[rank][j].data); 524 packPolygon(*elt, sendBuffer[rank], pos[rank]); 525 } 526 } 527 delete [] routingList; 528 529 530 int nbSendRequest = 0; 531 int nbRecvRequest = 0; 532 MPI_Request *sendRequest = new MPI_Request[mpiSize]; 533 MPI_Request *recvRequest = new MPI_Request[mpiSize]; 534 MPI_Status *status = new MPI_Status[mpiSize]; 535 536 for (int rank = 0; rank < mpiSize; rank++) 537 { 538 if (nbSendNode[rank] > 0) 539 { 540 MPI_Issend(sendBuffer[rank], sendMessageSize[rank], MPI_CHAR, rank, 0, communicator, &sendRequest[nbSendRequest]); 541 nbSendRequest++; 542 } 543 if (nbRecvNode[rank] > 0) 544 { 545 MPI_Irecv(recvBuffer[rank], recvMessageSize[rank], MPI_CHAR, rank, 0, communicator, &recvRequest[nbRecvRequest]); 546 nbRecvRequest++; 547 } 548 } 549 550 MPI_Waitall(nbRecvRequest, recvRequest, status); 551 MPI_Waitall(nbSendRequest, sendRequest, status); 552 553 for (int rank = 0; rank < mpiSize; rank++) 554 if (nbSendNode[rank] > 0) delete [] sendBuffer[rank]; 555 delete [] sendBuffer; 556 557 char **sendBuffer2 = new char*[mpiSize]; 558 char **recvBuffer2 = new char*[mpiSize]; 559 560 for (int rank = 0; rank < mpiSize; rank++) 561 { 562 nbSendNode[rank] = 0; 563 sendMessageSize[rank] = 0; 564 565 if (nbRecvNode[rank] > 0) 566 { 567 set<NodePtr> neighbourList; 568 pos[rank] = 0; 569 for (int j = 0; j < nbRecvNode[rank]; j++) 570 { 571 Elt elt; 572 unpackPolygon(elt, recvBuffer[rank], pos[rank]); 573 Node node(elt.x, cptRadius(elt), &elt); 574 findNeighbour(sstree.localTree.root, &node, neighbourList); 575 } 576 nbSendNode[rank] = neighbourList.size(); 577 for (set<NodePtr>::iterator it = neighbourList.begin(); it != neighbourList.end(); it++) 578 { 579 Elt *elt = (Elt *) ((*it)->data); 580 sendMessageSize[rank] += packedPolygonSize(*elt); 581 } 582 583 sendBuffer2[rank] = new char[sendMessageSize[rank]]; 584 pos[rank] = 0; 585 586 for (set<NodePtr>::iterator it = neighbourList.begin(); it != neighbourList.end(); it++) 587 { 588 Elt *elt = (Elt *) ((*it)->data); 589 packPolygon(*elt, sendBuffer2[rank], pos[rank]); 590 } 591 } 592 } 593 for (int rank = 0; rank < mpiSize; rank++) 594 if (nbRecvNode[rank] > 0) delete [] recvBuffer[rank]; 595 delete [] recvBuffer; 596 597 598 MPI_Barrier(communicator); 599 MPI_Alltoall(nbSendNode, 1, MPI_INT, nbRecvNode, 1, MPI_INT, communicator); 600 MPI_Alltoall(sendMessageSize, 1, MPI_INT, recvMessageSize, 1, MPI_INT, communicator); 601 602 for (int rank = 0; rank < mpiSize; rank++) 603 if (nbRecvNode[rank] > 0) recvBuffer2[rank] = new char[recvMessageSize[rank]]; 604 605 nbSendRequest = 0; 606 nbRecvRequest = 0; 607 608 for (int rank = 0; rank < mpiSize; rank++) 609 { 610 if (nbSendNode[rank] > 0) 611 { 612 MPI_Issend(sendBuffer2[rank], sendMessageSize[rank], MPI_CHAR, rank, 0, communicator, &sendRequest[nbSendRequest]); 613 nbSendRequest++; 614 } 615 if (nbRecvNode[rank] > 0) 616 { 617 MPI_Irecv(recvBuffer2[rank], recvMessageSize[rank], MPI_CHAR, rank, 0, communicator, &recvRequest[nbRecvRequest]); 618 nbRecvRequest++; 619 } 620 } 621 622 MPI_Waitall(nbRecvRequest, recvRequest, status); 623 MPI_Waitall(nbSendRequest, sendRequest, status); 624 625 int nbNeighbourNodes = 0; 626 for (int rank = 0; rank < mpiSize; rank++) 627 nbNeighbourNodes += nbRecvNode[rank]; 628 629 neighbourElements = new Elt[nbNeighbourNodes]; 630 nbNeighbourElements = nbNeighbourNodes; 631 632 int index = 0; 633 for (int rank = 0; rank < mpiSize; rank++) 634 { 635 pos[rank] = 0; 636 for (int j = 0; j < nbRecvNode[rank]; j++) 637 { 638 unpackPolygon(neighbourElements[index], recvBuffer2[rank], pos[rank]); 639 neighbourElements[index].id.ind = sstree.localTree.leafs.size() + index; 640 index++; 641 } 642 } 643 for (int rank = 0; rank < mpiSize; rank++) 644 { 645 if (nbRecvNode[rank] > 0) delete [] recvBuffer2[rank]; 646 if (nbSendNode[rank] > 0) delete [] sendBuffer2[rank]; 647 } 648 delete [] recvBuffer2; 649 delete [] sendBuffer2; 650 delete [] sendMessageSize; 651 delete [] recvMessageSize; 652 delete [] nbSendNode; 653 delete [] nbRecvNode; 654 delete [] sendRequest; 655 delete [] recvRequest; 656 delete [] status; 657 delete [] pos; 658 659 /* re-compute on received elements to avoid having to send this information */ 660 neighbourNodes.resize(nbNeighbourNodes); 661 setCirclesAndLinks(neighbourElements, neighbourNodes); 662 cptAllEltsGeom(neighbourElements, nbNeighbourNodes, srcGrid.pole); 663 664 /* the local SS tree must include nodes from other cpus if they are potential 665 intersector of a local node */ 666 sstree.localTree.insertNodes(neighbourNodes); 667 668 /* for every local element, 669 use the SS-tree to find all elements (including neighbourElements) 670 who are potential neighbours because their circles intersect, 671 then check all canditates for common edges to build up connectivity information 672 */ 673 for (int j = 0; j < sstree.localTree.leafs.size(); j++) 674 { 675 Node& node = sstree.localTree.leafs[j]; 676 677 /* find all leafs whoes circles that intersect node's circle and save into node->intersectors */ 678 node.search(sstree.localTree.root); 679 680 Elt *elt = (Elt *)(node.data); 681 682 for (int i = 0; i < elt->n; i++) elt->neighbour[i] = NOT_FOUND; 683 684 /* for element `elt` loop through all nodes in the SS-tree 685 whoes circles intersect with the circle around `elt` (the SS intersectors) 686 and check if they are neighbours in the sense that the two elements share an edge. 687 If they do, save this information for elt */ 688 for (list<NodePtr>::iterator it = (node.intersectors).begin(); it != (node.intersectors).end(); ++it) 689 { 690 Elt *elt2 = (Elt *)((*it)->data); 691 set_neighbour(*elt, *elt2); 692 } 693 694 } 702 695 } 703 696 … … 705 698 void Mapper::computeIntersection(Elt *elements, int nbElements) 706 699 { 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 //intersect_ym(&recvElt[j], elt2);834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870 871 872 873 874 875 876 877 878 879 880 881 882 883 884 885 886 887 888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 916 700 int mpiSize, mpiRank; 701 MPI_Comm_size(communicator, &mpiSize); 702 MPI_Comm_rank(communicator, &mpiRank); 703 704 MPI_Barrier(communicator); 705 706 vector<Node> *routingList = new vector<Node>[mpiSize]; 707 708 vector<Node> routeNodes; routeNodes.reserve(nbElements); 709 for (int j = 0; j < nbElements; j++) 710 { 711 elements[j].id.ind = j; 712 elements[j].id.rank = mpiRank; 713 routeNodes.push_back(Node(elements[j].x, cptRadius(elements[j]), &elements[j])); 714 } 715 716 vector<vector<int> > routes(routeNodes.size()); 717 sstree.routeIntersections(routes, routeNodes); 718 for (int i = 0; i < routes.size(); ++i) 719 for (int k = 0; k < routes[i].size(); ++k) 720 routingList[routes[i][k]].push_back(routeNodes[i]); 721 722 if (verbose >= 2) 723 { 724 cout << " --> rank " << mpiRank << " nbElements " << nbElements << " : "; 725 for (int rank = 0; rank < mpiSize; rank++) 726 cout << routingList[rank].size() << " "; 727 cout << endl; 728 } 729 MPI_Barrier(communicator); 730 731 int *nbSendNode = new int[mpiSize]; 732 int *nbRecvNode = new int[mpiSize]; 733 int *sentMessageSize = new int[mpiSize]; 734 int *recvMessageSize = new int[mpiSize]; 735 736 for (int rank = 0; rank < mpiSize; rank++) 737 { 738 nbSendNode[rank] = routingList[rank].size(); 739 sentMessageSize[rank] = 0; 740 for (size_t j = 0; j < routingList[rank].size(); j++) 741 { 742 Elt *elt = (Elt *) (routingList[rank][j].data); 743 sentMessageSize[rank] += packedPolygonSize(*elt); 744 } 745 } 746 747 MPI_Alltoall(nbSendNode, 1, MPI_INT, nbRecvNode, 1, MPI_INT, communicator); 748 MPI_Alltoall(sentMessageSize, 1, MPI_INT, recvMessageSize, 1, MPI_INT, communicator); 749 750 int total = 0; 751 752 for (int rank = 0; rank < mpiSize; rank++) 753 { 754 total = total + nbRecvNode[rank]; 755 } 756 757 if (verbose >= 2) cout << "---> rank " << mpiRank << " : compute intersection : total received nodes " << total << endl; 758 char **sendBuffer = new char*[mpiSize]; 759 char **recvBuffer = new char*[mpiSize]; 760 int *pos = new int[mpiSize]; 761 762 for (int rank = 0; rank < mpiSize; rank++) 763 { 764 if (nbSendNode[rank] > 0) sendBuffer[rank] = new char[sentMessageSize[rank]]; 765 if (nbRecvNode[rank] > 0) recvBuffer[rank] = new char[recvMessageSize[rank]]; 766 } 767 768 for (int rank = 0; rank < mpiSize; rank++) 769 { 770 pos[rank] = 0; 771 for (size_t j = 0; j < routingList[rank].size(); j++) 772 { 773 Elt* elt = (Elt *) (routingList[rank][j].data); 774 packPolygon(*elt, sendBuffer[rank], pos[rank]); 775 } 776 } 777 delete [] routingList; 778 779 int nbSendRequest = 0; 780 int nbRecvRequest = 0; 781 MPI_Request *sendRequest = new MPI_Request[mpiSize]; 782 MPI_Request *recvRequest = new MPI_Request[mpiSize]; 783 MPI_Status *status = new MPI_Status[mpiSize]; 784 785 for (int rank = 0; rank < mpiSize; rank++) 786 { 787 if (nbSendNode[rank] > 0) 788 { 789 MPI_Issend(sendBuffer[rank], sentMessageSize[rank], MPI_CHAR, rank, 0, communicator, &sendRequest[nbSendRequest]); 790 nbSendRequest++; 791 } 792 if (nbRecvNode[rank] > 0) 793 { 794 MPI_Irecv(recvBuffer[rank], recvMessageSize[rank], MPI_CHAR, rank, 0, communicator, &recvRequest[nbRecvRequest]); 795 nbRecvRequest++; 796 } 797 } 798 799 MPI_Waitall(nbRecvRequest, recvRequest, status); 800 MPI_Waitall(nbSendRequest, sendRequest, status); 801 802 char **sendBuffer2 = new char*[mpiSize]; 803 char **recvBuffer2 = new char*[mpiSize]; 804 805 double tic = cputime(); 806 for (int rank = 0; rank < mpiSize; rank++) 807 { 808 sentMessageSize[rank] = 0; 809 810 if (nbRecvNode[rank] > 0) 811 { 812 Elt *recvElt = new Elt[nbRecvNode[rank]]; 813 pos[rank] = 0; 814 for (int j = 0; j < nbRecvNode[rank]; j++) 815 { 816 unpackPolygon(recvElt[j], recvBuffer[rank], pos[rank]); 817 cptEltGeom(recvElt[j], tgtGrid.pole); 818 Node recvNode(recvElt[j].x, cptRadius(recvElt[j]), &recvElt[j]); 819 recvNode.search(sstree.localTree.root); 820 /* for a node holding an element of the target, loop throught candidates for intersecting source */ 821 for (list<NodePtr>::iterator it = (recvNode.intersectors).begin(); it != (recvNode.intersectors).end(); ++it) 822 { 823 Elt *elt2 = (Elt *) ((*it)->data); 824 /* recvElt is target, elt2 is source */ 825 // intersect(&recvElt[j], elt2); 826 intersect_ym(&recvElt[j], elt2); 827 } 828 if (recvElt[j].is.size() > 0) sentMessageSize[rank] += packIntersectionSize(recvElt[j]); 829 830 // here recvNode goes out of scope 831 } 832 833 if (sentMessageSize[rank] > 0) 834 { 835 sentMessageSize[rank] += sizeof(int); 836 sendBuffer2[rank] = new char[sentMessageSize[rank]]; 837 *((int *) sendBuffer2[rank]) = 0; 838 pos[rank] = sizeof(int); 839 for (int j = 0; j < nbRecvNode[rank]; j++) 840 { 841 packIntersection(recvElt[j], sendBuffer2[rank], pos[rank]); 842 //FIXME should be deleted: recvElt[j].delete_intersections(); // intersection areas have been packed to buffer and won't be used any more 843 } 844 } 845 delete [] recvElt; 846 } 847 } 848 delete [] pos; 849 850 for (int rank = 0; rank < mpiSize; rank++) 851 { 852 if (nbSendNode[rank] > 0) delete [] sendBuffer[rank]; 853 if (nbRecvNode[rank] > 0) delete [] recvBuffer[rank]; 854 nbSendNode[rank] = 0; 855 } 856 857 if (verbose >= 2) cout << "Rank " << mpiRank << " Compute (internal) intersection " << cputime() - tic << " s" << endl; 858 MPI_Alltoall(sentMessageSize, 1, MPI_INT, recvMessageSize, 1, MPI_INT, communicator); 859 860 for (int rank = 0; rank < mpiSize; rank++) 861 if (recvMessageSize[rank] > 0) 862 recvBuffer2[rank] = new char[recvMessageSize[rank]]; 863 864 nbSendRequest = 0; 865 nbRecvRequest = 0; 866 867 for (int rank = 0; rank < mpiSize; rank++) 868 { 869 if (sentMessageSize[rank] > 0) 870 { 871 MPI_Issend(sendBuffer2[rank], sentMessageSize[rank], MPI_CHAR, rank, 0, communicator, &sendRequest[nbSendRequest]); 872 nbSendRequest++; 873 } 874 if (recvMessageSize[rank] > 0) 875 { 876 MPI_Irecv(recvBuffer2[rank], recvMessageSize[rank], MPI_CHAR, rank, 0, communicator, &recvRequest[nbRecvRequest]); 877 nbRecvRequest++; 878 } 879 } 880 881 MPI_Waitall(nbRecvRequest, recvRequest, status); 882 MPI_Waitall(nbSendRequest, sendRequest, status); 883 884 delete [] sendRequest; 885 delete [] recvRequest; 886 delete [] status; 887 for (int rank = 0; rank < mpiSize; rank++) 888 { 889 if (nbRecvNode[rank] > 0) 890 { 891 if (sentMessageSize[rank] > 0) 892 delete [] sendBuffer2[rank]; 893 } 894 895 if (recvMessageSize[rank] > 0) 896 { 897 unpackIntersection(elements, recvBuffer2[rank]); 898 delete [] recvBuffer2[rank]; 899 } 900 } 901 delete [] sendBuffer2; 902 delete [] recvBuffer2; 903 delete [] sendBuffer; 904 delete [] recvBuffer; 905 906 delete [] nbSendNode; 907 delete [] nbRecvNode; 908 delete [] sentMessageSize; 909 delete [] recvMessageSize; 917 910 } 918 911 919 912 Mapper::~Mapper() 920 913 { 921 922 923 924 925 926 } 927 928 } 914 delete [] remapMatrix; 915 delete [] srcAddress; 916 delete [] srcRank; 917 delete [] dstAddress; 918 if (neighbourElements) delete [] neighbourElements; 919 } 920 921 } -
XIOS/dev/branch_openmp/extern/remap/src/mapper.hpp
r1134 r1328 3 3 #include "parallel_tree.hpp" 4 4 #include "mpi.hpp" 5 6 #ifdef _usingEP7 #include "ep_declaration.hpp"8 #endif9 5 10 6 namespace sphereRemap { … … 22 18 { 23 19 public: 24 Mapper(ep_lib::MPI_Comm comm =MPI_COMM_WORLD) : communicator(comm), verbose(SILENT), neighbourElements(NULL), sstree(comm) {}20 Mapper(ep_lib::MPI_Comm comm) : communicator(comm), verbose(SILENT), neighbourElements(NULL), sstree(comm) {} 25 21 ~Mapper(); 26 22 void setVerbosity(verbosity v) {verbose=v ;} -
XIOS/dev/branch_openmp/extern/remap/src/mpi_cascade.cpp
r688 r1328 1 1 #include "mpi_cascade.hpp" 2 2 #include <iostream> 3 using namespace ep_lib; 3 4 4 5 namespace sphereRemap { -
XIOS/dev/branch_openmp/extern/remap/src/mpi_cascade.hpp
r694 r1328 12 12 { 13 13 public: 14 CCascadeLevel(MPI_Comm comm) : comm(comm)15 16 17 18 19 20 14 CCascadeLevel(ep_lib::MPI_Comm comm) : comm(comm) 15 { 16 ep_lib::MPI_Comm_size(comm, &size); 17 ep_lib::MPI_Comm_rank(comm, &rank); 18 } 19 int colour() const { return rank % group_size; }; 20 int key() const { return p_colour() + rank/(p_grp_size*group_size)*p_grp_size; } 21 21 22 23 24 22 // perpendicular group 23 int p_colour() const { return (rank%group_size + rank/group_size) % p_grp_size; } 24 int p_key() const { return colour() + rank/(p_grp_size*group_size)*group_size; } 25 25 26 27 28 29 30 26 ep_lib::MPI_Comm comm, pg_comm; 27 int rank; 28 int size; 29 int group_size; // group_size and p_grp_size are interchanged?? FIXME 30 int p_grp_size; 31 31 }; 32 32 … … 34 34 { 35 35 public: 36 // 37 CMPICascade(int nodes_per_level, MPI_Comm comm); 36 CMPICascade(int nodes_per_level, ep_lib::MPI_Comm comm); 38 37 39 40 38 int num_levels; 39 std::vector<CCascadeLevel> level; 41 40 }; 42 41 -
XIOS/dev/branch_openmp/extern/remap/src/mpi_routing.cpp
r1289 r1328 5 5 #include "timerRemap.hpp" 6 6 #include <iostream> 7 #ifdef _usingEP 8 #include "ep_declaration.hpp" 9 #endif 7 using namespace ep_lib; 10 8 11 9 namespace sphereRemap { … … 125 123 CTimer::get("CMPIRouting::init(reduce_scatter)").print(); 126 124 127 MPI_Info info_null;128 129 // MPI_Alloc_mem(nbTarget *sizeof(int), info_null, &targetRank);130 // MPI_Alloc_mem(nbSource *sizeof(int), info_null, &sourceRank);131 125 MPI_Alloc_mem(nbTarget *sizeof(int), MPI_INFO_NULL, &targetRank); 132 126 MPI_Alloc_mem(nbSource *sizeof(int), MPI_INFO_NULL, &sourceRank); … … 158 152 { 159 153 #ifdef _usingEP 160 MPI_Irecv(&sourceRank[i], 1, MPI_INT, - 1, 0, communicator, &request[indexRequest]);154 MPI_Irecv(&sourceRank[i], 1, MPI_INT, -2, 0, communicator, &request[indexRequest]); 161 155 #else 162 156 MPI_Irecv(&sourceRank[i], 1, MPI_INT, MPI_ANY_SOURCE, 0, communicator, &request[indexRequest]); … … 182 176 { 183 177 #ifdef _usingEP 184 MPI_Irecv(&sourceRank[i], 1, MPI_INT, - 1, 0, communicator, &request[indexRequest]);178 MPI_Irecv(&sourceRank[i], 1, MPI_INT, -2, 0, communicator, &request[indexRequest]); 185 179 #else 186 180 MPI_Irecv(&sourceRank[i], 1, MPI_INT, MPI_ANY_SOURCE, 0, communicator, &request[indexRequest]); -
XIOS/dev/branch_openmp/extern/remap/src/mpi_routing.hpp
r694 r1328 11 11 { 12 12 13 MPI_Comm communicator;13 ep_lib::MPI_Comm communicator; 14 14 int mpiRank; 15 15 int mpiSize; … … 29 29 30 30 public: 31 CMPIRouting( MPI_Comm comm);31 CMPIRouting(ep_lib::MPI_Comm comm); 32 32 ~CMPIRouting(); 33 33 template<typename T> void init(const std::vector<T>& route, CMPICascade *cascade = NULL); … … 44 44 template <typename T> 45 45 void alltoalls_known(const std::vector<std::vector<T> >& send, std::vector<std::vector<T> >& recv, 46 const std::vector<int>& ranks, MPI_Comm communicator);46 const std::vector<int>& ranks, ep_lib::MPI_Comm communicator); 47 47 48 48 template <typename T> 49 49 void alltoalls_unknown(const std::vector<std::vector<T> >& send, std::vector<std::vector<T> >& recv, 50 const std::vector<int>& ranks, MPI_Comm communicator);50 const std::vector<int>& ranks, ep_lib::MPI_Comm communicator); 51 51 } 52 52 #endif -
XIOS/dev/branch_openmp/extern/remap/src/node.cpp
r1205 r1328 254 254 NodePtr insert(NodePtr thIs, NodePtr node) 255 255 { 256 257 258 259 260 261 262 263 264 256 int la = thIs->level; // node to be inserted 257 int lb = node->level; // node where insertation 258 assert(la < lb); // node to be inserted must have lower level then parent 259 //if (thIs->parent) assert(find_in_tree1(thIs) == true); 260 NodePtr q = NULL; 261 NodePtr chd = NULL; 262 node->move(thIs); 263 if (la == lb - 1) 264 { 265 265 node->child.push_back(thIs); 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 266 thIs->parent = node; 267 if (node->child.size() > MAX_NODE_SZ && node->tree->canSplit() ) // with us as additional child `node` is now too large :( 268 return (node->reinserted || node->parent == NULL) ? split(node) : reinsert(node); 269 } 270 else // la < lb - 1 271 { 272 chd = thIs->closest(node->child); 273 q = insert(thIs, chd); 274 } 275 if ((node->updateCount + 1) % UPDATE_EVERY == 0) 276 node->update(); 277 else 278 { 279 if (q) node->remove(q); 280 node->inflate(chd); 281 } 282 282 283 283 return q; -
XIOS/dev/branch_openmp/extern/remap/src/node.hpp
r1153 r1328 15 15 struct Circle 16 16 { 17 18 17 Coord centre; 18 double radius; 19 19 }; 20 20 … … 116 116 struct Node 117 117 { 118 119 120 121 122 123 124 125 126 127 128 129 118 int level; /* FIXME leafs are 0 and root is max level? */ 119 int leafCount; /* number of leafs that are descendants of this node (the elements in it's cycle) */ 120 Coord centre; 121 double radius; 122 NodePtr parent, ref; 123 std::vector<NodePtr> child; 124 std::list<NodePtr> intersectors; 125 bool reinserted; 126 int updateCount; // double var; 127 CBasicTree* tree; 128 void *data; 129 int route; 130 130 bool toDelete ; 131 131 132 133 134 132 Node() : level(0), leafCount(1), centre(ORIGIN), radius(0), reinserted(false), updateCount(0), toDelete(false) {} 133 Node(const Coord& centre, double radius, void *data) 134 : level(0), leafCount(1), centre(centre), radius(radius), reinserted(false), updateCount(0), data(data), toDelete(false) {} 135 135 136 136 //#ifdef DEBUG … … 178 178 //#endif 179 179 180 181 182 183 180 void move(const NodePtr node); 181 void remove(const NodePtr node); 182 void inflate(const NodePtr node); 183 void update(); 184 184 void output(std::ostream& flux, int level, int color) ; 185 186 187 188 189 190 191 192 193 185 NodePtr closest(std::vector<NodePtr>& list, int n = CLOSEST); 186 NodePtr farthest(std::vector<NodePtr>& list); 187 void findClosest(int level, NodePtr src, double& minDist, NodePtr &closest); 188 189 void search(NodePtr node); 190 bool centreInside(Node &node); 191 bool intersects(NodePtr node); 192 bool isInside(Node &node); 193 int incluCheck(); 194 194 void checkParent(void) ; 195 196 197 198 199 200 201 195 void printChildren(); 196 void assignRoute(std::vector<int>::iterator& rank, int level); 197 void assignCircleAndPropagateUp(Coord *centres, double *radia, int level); 198 void printLevel(int level); 199 void routeNode(NodePtr node, int level); 200 void routingIntersecting(std::vector<Node>* routingList, NodePtr node); 201 void routeIntersection(std::vector<int>& routes, NodePtr node); 202 202 void getNodeLevel(int level,std::list<NodePtr>& NodeList) ; 203 203 bool removeDeletedNodes(int assignLevel) ; -
XIOS/dev/branch_openmp/extern/remap/src/parallel_tree.cpp
r1295 r1328 12 12 13 13 #include "parallel_tree.hpp" 14 using namespace ep_lib; 14 15 15 16 namespace sphereRemap { 16 17 extern CRemapGrid srcGrid;18 #pragma omp threadprivate(srcGrid)19 20 extern CRemapGrid tgtGrid;21 #pragma omp threadprivate(tgtGrid)22 17 23 18 static const int assignLevel = 2; … … 292 287 { 293 288 294 295 289 int assignLevel = 2; 290 int nbSampleNodes = 2*ipow(MAX_NODE_SZ + 1, assignLevel); 296 291 297 292 … … 305 300 MPI_Comm_size(communicator,&commSize) ; 306 301 307 308 309 310 //assert( nbTot > nbSampleNodes*commSize) ;302 // make multiple of two 303 nbSampleNodes /= 2; 304 nbSampleNodes *= 2; 305 // assert( nbTot > nbSampleNodes*commSize) ; 311 306 312 307 int nbSampleNodes1 = nbSampleNodes * (nb1*commSize)/(1.*nbTot) ; … … 314 309 315 310 316 //assert(node.size() > nbSampleNodes);317 //assert(node2.size() > nbSampleNodes);318 //assert(node.size() + node2.size() > nbSampleNodes);319 320 321 322 323 324 325 326 for (int i = 0; i <nbSampleNodes1; i++) sampleNodes.push_back(Node(node[randomArray1[i%nb1]].centre, node[randomArray1[i%nb1]].radius, NULL));327 for (int i = 0; i <nbSampleNodes2; i++) sampleNodes.push_back(Node(node2[randomArray2[i%nb2]].centre, node2[randomArray2[i%nb2]].radius, NULL));311 // assert(node.size() > nbSampleNodes); 312 // assert(node2.size() > nbSampleNodes); 313 // assert(node.size() + node2.size() > nbSampleNodes); 314 vector<Node> sampleNodes; sampleNodes.reserve(nbSampleNodes1+nbSampleNodes2); 315 316 vector<int> randomArray1(node.size()); 317 randomizeArray(randomArray1); 318 vector<int> randomArray2(node2.size()); 319 randomizeArray(randomArray2); 320 321 for (int i = 0; i <nbSampleNodes1; i++) sampleNodes.push_back(Node(node[randomArray1[i%nb1]].centre, node[randomArray1[i%nb1]].radius, NULL)); 322 for (int i = 0; i <nbSampleNodes2; i++) sampleNodes.push_back(Node(node2[randomArray2[i%nb2]].centre, node2[randomArray2[i%nb2]].radius, NULL)); 328 323 329 324 CTimer::get("buildParallelSampleTree").resume(); … … 336 331 CTimer::get("parallelRouteNode").resume(); 337 332 vector<int> route(node.size()); 333 cout<<"node.size = "<<node.size()<<endl; 338 334 routeNodes(route /*out*/, node); 339 335 CTimer::get("parallelRouteNode").suspend(); -
XIOS/dev/branch_openmp/extern/remap/src/parallel_tree.hpp
r1134 r1328 6 6 #include "mpi_cascade.hpp" 7 7 #include "mpi.hpp" 8 #ifdef _usingEP9 #include "ep_declaration.hpp"10 #endif11 12 8 namespace sphereRemap { 13 9 … … 15 11 { 16 12 public: 17 18 13 CParallelTree(ep_lib::MPI_Comm comm); 14 ~CParallelTree(); 19 15 20 16 void build(vector<Node>& node, vector<Node>& node2); 21 17 22 23 18 void routeNodes(vector<int>& route, vector<Node>& nodes, int level = 0); 19 void routeIntersections(vector<vector<int> >& route, vector<Node>& nodes, int level = 0); 24 20 25 26 21 int nbLocalElements; 22 Elt* localElements; 27 23 28 24 CTree localTree; 29 25 30 26 private: 31 32 33 34 27 void updateCirclesForRouting(Coord rootCentre, double rootRadius, int level = 0); 28 void buildSampleTreeCascade(vector<Node>& sampleNodes, int level = 0); 29 void buildLocalTree(const vector<Node>& node, const vector<int>& route); 30 void buildRouteTree(); 35 31 36 37 38 39 32 //CSampleTree sampleTree; 33 vector<CSampleTree> treeCascade; // first for sample tree, then for routing tree 34 CMPICascade cascade; 35 ep_lib::MPI_Comm communicator ; 40 36 41 37 }; -
XIOS/dev/branch_openmp/extern/remap/src/polyg.cpp
r1289 r1328 7 7 8 8 #include "polyg.hpp" 9 #include <stdio.h>10 9 11 10 namespace sphereRemap { … … 46 45 Coord barycentre(const Coord *x, int n) 47 46 { 48 49 47 if (n == 0) return ORIGIN; 50 48 Coord bc = ORIGIN; 51 49 for (int i = 0; i < n; i++) 52 {53 50 bc = bc + x[i]; 54 }55 51 /* both distances can be equal down to roundoff when norm(bc) < mashineepsilon 56 52 which can occur when weighted with tiny area */ 57 58 59 //assert(squaredist(bc, proj(bc)) <= squaredist(bc, proj(bc * (-1.0)))); 53 54 assert(squaredist(bc, proj(bc)) <= squaredist(bc, proj(bc * (-1.0)))); 55 //if (squaredist(bc, proj(bc)) > squaredist(bc, proj(bc * (-1.0)))) return proj(bc * (-1.0)); 60 56 61 57 return proj(bc); … … 165 161 { 166 162 if (N < 3) 167 163 return 0; /* polygons with less than three vertices have zero area */ 168 164 Coord t[3]; 169 165 t[0] = barycentre(x, N); … … 177 173 t[1] = x[i]; 178 174 t[2] = x[ii]; 179 double sc=scalarprod(crossprod(t[1] - t[0], t[2] - t[0]), t[0]) ;175 double sc=scalarprod(crossprod(t[1] - t[0], t[2] - t[0]), t[0]) ; 180 176 assert(sc >= -1e-10); // Error: tri a l'env (wrong orientation) 181 177 double area_gc = triarea(t[0], t[1], t[2]); 182 //if(area_gc<=0) printf("area_gc = %e\n", area_gc);183 178 double area_sc_gc_moon = 0; 184 179 if (d[i]) /* handle small circle case */ … … 188 183 char sgl = (mext > 0) ? -1 : 1; 189 184 area_sc_gc_moon = sgl * alun(arcdist(t[1], t[2]), fabs(scalarprod(t[1], pole))); 190 //if(area_sc_gc_moon<=0) printf("area_sc_gc_moon = %e\n", area_sc_gc_moon);191 185 gg_exact = gg_exact + sc_gc_moon_normalintegral(t[1], t[2], pole); 192 186 } 193 187 area += area_gc + area_sc_gc_moon; /* for "spherical circle segment" sum triangular part (at) and "small moon" and => account for small circle */ 194 188 g[i] = barycentre(t, 3) * (area_gc + area_sc_gc_moon); 195 //printf("g[%d] = (%e,%e,%e) * (%e+%e) = (%e,%e,%e) norm = %e\n", i, barycentre(t, 3).x, barycentre(t, 3).y, barycentre(t, 3).z, area_gc, area_sc_gc_moon, g[i].x, g[i].y, g[i].z, norm(g[i]));196 189 } 197 190 gg = barycentre(g, N); -
XIOS/dev/branch_openmp/extern/remap/src/timerRemap.cpp
r1146 r1328 4 4 #include <map> 5 5 #include <iostream> 6 using namespace ep_lib; 6 7 7 8 namespace sphereRemap { … … 9 10 using namespace std; 10 11 11 map<string,CTimer*> *CTimer::allTimer = 0; 12 //map<string,CTimer*> CTimer::allTimer; 13 map<string,CTimer*> *CTimer::allTimer_ptr = 0; 12 14 13 15 CTimer::CTimer(const string& name_) : name(name_) … … 55 57 CTimer& CTimer::get(const string name) 56 58 { 57 if(allTimer == 0) allTimer = new map<string,CTimer*>;58 59 map<string,CTimer*>::iterator it; 59 it=(*allTimer).find(name); 60 if (it==(*allTimer).end()) it=(*allTimer).insert(pair<string,CTimer*>(name,new CTimer(name))).first; 60 if(allTimer_ptr == 0) allTimer_ptr = new map<string,CTimer*>; 61 //it=allTimer.find(name); 62 it=allTimer_ptr->find(name); 63 //if (it==allTimer.end()) it=allTimer.insert(pair<string,CTimer*>(name,new CTimer(name))).first; 64 if (it==allTimer_ptr->end()) it=allTimer_ptr->insert(pair<string,CTimer*>(name,new CTimer(name))).first; 61 65 return *(it->second); 62 66 } -
XIOS/dev/branch_openmp/extern/remap/src/timerRemap.hpp
r1146 r1328 27 27 void print(void); 28 28 //static map<string,CTimer*> allTimer; 29 static map<string,CTimer*> *allTimer; 30 #pragma omp threadprivate(allTimer) 31 29 static map<string,CTimer*> *allTimer_ptr; 32 30 static double getTime(void); 33 31 static CTimer& get(string name); -
XIOS/dev/branch_openmp/extern/remap/src/tree.cpp
r1172 r1328 29 29 void CBasicTree::routeNodes(vector<int>& route, vector<Node>& nodes, int assignLevel) 30 30 { 31 32 33 34 35 31 for (int i = 0; i < nodes.size(); i++) 32 { 33 root->routeNode(&nodes[i], assignLevel); 34 route[i] = nodes[i].route; 35 } 36 36 } 37 37 38 38 void CBasicTree::routeIntersections(vector<vector<int> >& routes, vector<Node>& nodes) 39 39 { 40 41 40 for (int i = 0; i < nodes.size(); i++) 41 root->routeIntersection(routes[i], &nodes[i]); 42 42 } 43 43 44 44 void CBasicTree::build(vector<Node>& nodes) 45 45 { 46 47 46 newRoot(1); 47 insertNodes(nodes); 48 48 } 49 49 50 50 void CBasicTree::output(ostream& flux, int level) 51 51 { 52 52 root->output(flux,level,0) ; 53 53 } 54 54 void CBasicTree::slim(int nbIts) 55 55 { 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 56 for (int i = 0; i < nbIts; i++) 57 { 58 for (int level = root->level - 1; level > 0; level--) 59 { 60 slim2(root, level); 61 ri = 0; 62 emptyPool(); 63 } 64 65 for (int level = 2; level < root->level; level++) 66 { 67 slim2(root, level); 68 ri = 0; 69 emptyPool(); 70 } 71 } 72 72 } 73 73 … … 76 76 void CBasicTree::insertNode(NodePtr node) 77 77 { 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 78 node->tree = this; 79 increaseLevelSize(0); 80 push_back(node); 81 82 NodePtr q; 83 while (pool.size()) 84 { 85 q = pool.front(); 86 pool.pop_front(); 87 q = insert(q, root); 88 if (ri) 89 { 90 delete q; 91 ri = 0; 92 } 93 } 94 94 } 95 95 96 96 void CBasicTree::emptyPool(void) 97 97 { 98 99 100 101 102 103 104 105 106 107 108 98 while (pool.size()) 99 { 100 NodePtr q = pool.front(); 101 pool.pop_front(); 102 q = insert(q, root); 103 if (ri) 104 { 105 delete q; 106 ri = 0; 107 } 108 } 109 109 } 110 110 … … 142 142 root->parent = 0; 143 143 root->leafCount = 0; 144 // initialize root node on the sphere 145 root->centre.x=1 ; 146 root->centre.y=0 ; 147 root->centre.z=0 ; 144 // initialize root node on the sphere 145 root->centre.x=1 ; root->centre.y=0 ; root->centre.z=0 ; 148 146 root->radius = 0.; 149 147 root->reinserted = false; -
XIOS/dev/branch_openmp/extern/remap/src/tree.hpp
r1172 r1328 14 14 class CBasicTree 15 15 { 16 16 public: 17 17 18 19 20 21 22 18 NodePtr root; /* The main tree is stored as Nodes which can be reached through traversal starting here */ 19 NodePtr ref; // FIXME this reference, set by a node is odd, try to remove 20 int ri; /** this is set to one by a node in case of reinsertion */ 21 vector<int> levelSize; /** e.g. levelSize[0] == leafs.size() */ 22 vector<Node> leafs; /** leafs are stored in vector for easy access and rest of the tree nodes as separate allocations, only reachable through tree traversal */ 23 23 24 25 26 27 28 24 CBasicTree() : ri(0), levelSize(MAX_LEVEL_SIZE), root(NULL), isAssignedLevel(false), okSplit(true), isActiveOkSplit(false) {} 25 ~CBasicTree(); 26 void build(vector<Node>& nodes); 27 void slim(int nbIts = 1); 28 virtual void insertNodes(vector<Node>& node) = 0; 29 29 30 31 30 void routeNodes(vector<int>& route, vector<Node>& nodes, int assignLevel); 31 void routeIntersections(vector<vector<int> >& route, vector<Node>& nodes); 32 32 33 34 35 36 37 38 33 void push_back(NodePtr node); 34 void push_front(NodePtr node); 35 void increaseLevelSize(int level); 36 void decreaseLevelSize(int level); 37 void newRoot(int level); 38 void insertNode(NodePtr node); 39 39 void output(ostream& flux, int level) ; 40 40 41 41 int keepNodes; 42 42 bool isAssignedLevel ; 43 43 int assignLevel; … … 50 50 51 51 52 53 52 private: 53 deque<NodePtr > pool; 54 54 55 55 bool okSplit ; 56 56 57 57 protected: 58 58 void emptyPool(); 59 CBasicTree(int keepNodes_, int assignLevel_) : ri(0), levelSize(MAX_LEVEL_SIZE), root(NULL), keepNodes(keepNodes_), assignLevel(assignLevel_), isAssignedLevel(true), 60 okSplit(true), isActiveOkSplit(false) {} 59 CBasicTree(int keepNodes_, int assignLevel_) : ri(0), levelSize(MAX_LEVEL_SIZE), root(NULL), keepNodes(keepNodes_), assignLevel(assignLevel_), isAssignedLevel(true), okSplit(true), isActiveOkSplit(false) {} 61 60 }; 62 61 63 62 class CTree : public CBasicTree 64 63 { 65 66 64 public: 65 void insertNodes(vector<Node>& nodes); 67 66 }; 68 67 … … 70 69 { 71 70 72 73 71 public: 72 CSampleTree(int keepNodes_, int assignLevel_) : CBasicTree(keepNodes_,assignLevel_) {} 74 73 void slimAssignedLevel() ; 75 74 void removeExtraNode(void) ; 76 75 void insertNodes(vector<Node>& nodes); 77 76 }; 78 77 -
XIOS/dev/branch_openmp/extern/src_ep_dev/ep_fortran.cpp
r1287 r1328 31 31 { 32 32 fc_comm_map.insert(std::make_pair( std::make_pair( fint, omp_get_thread_num()) , comm)); 33 //printf("EP_Comm_c2f : MAP %p insert: %d, %d, %p\n", &fc_comm_map, fint, omp_get_thread_num(), comm.ep_comm_ptr);33 printf("EP_Comm_c2f : MAP %p insert: %d, %d, %p\n", &fc_comm_map, fint, omp_get_thread_num(), comm.ep_comm_ptr); 34 34 } 35 35 } … … 54 54 MPI_Comm comm_ptr; 55 55 comm_ptr = it->second; 56 //printf("EP_Comm_f2c : MAP %p find: %d, %d, %p\n", &fc_comm_map, it->first.first, it->first.second, comm_ptr.ep_comm_ptr);56 printf("EP_Comm_f2c : MAP %p find: %d, %d, %p\n", &fc_comm_map, it->first.first, it->first.second, comm_ptr.ep_comm_ptr); 57 57 return comm_ptr; 58 58 } -
XIOS/dev/branch_openmp/extern/src_ep_dev/ep_intercomm.cpp
r1287 r1328 47 47 48 48 MPI_Waitall(2, request, status); 49 50 //MPI_Send(&leader_ranks[0], 3, static_cast< ::MPI_Datatype>(MPI_INT), remote_leader, tag, peer_comm); 51 //MPI_Recv(&leader_ranks[3], 3, static_cast< ::MPI_Datatype>(MPI_INT), remote_leader, tag, peer_comm, &status[1]); 49 52 } 50 53 -
XIOS/dev/branch_openmp/extern/src_ep_dev/ep_lib.cpp
r1295 r1328 10 10 #pragma omp threadprivate(EP_PendingRequests) 11 11 12 13 12 14 namespace ep_lib 13 15 { 16 bool MPI_Comm::is_null() 17 { 18 if(!this->is_intercomm) 19 return this->mpi_comm == MPI_COMM_NULL.mpi_comm; 20 else 21 return this->ep_comm_ptr->intercomm->mpi_inter_comm == MPI_COMM_NULL.mpi_comm; 22 } 14 23 15 24 int tag_combine(int real_tag, int src, int dest) -
XIOS/dev/branch_openmp/extern/src_ep_dev/ep_send.cpp
r1295 r1328 19 19 if(!comm.is_ep) 20 20 return ::MPI_Send(buf, count, to_mpi_type(datatype), dest, tag, to_mpi_comm(comm.mpi_comm)); 21 22 MPI_Request request; 23 MPI_Status status; 24 MPI_Isend(buf, count, datatype, dest, tag, comm, &request); 25 MPI_Wait(&request, &status); 26 21 if(comm.is_intercomm) 22 { 23 MPI_Request request; 24 MPI_Status status; 25 MPI_Isend(buf, count, datatype, dest, tag, comm, &request); 26 MPI_Wait(&request, &status); 27 } 28 else 29 { 30 int ep_src_loc = comm.ep_comm_ptr->size_rank_info[1].first; 31 int ep_dest_loc = comm.ep_comm_ptr->comm_list->rank_map->at(dest).first; 32 int mpi_tag = tag_combine(tag, ep_src_loc, ep_dest_loc); 33 int mpi_dest = comm.ep_comm_ptr->comm_list->rank_map->at(dest).second; 34 35 ::MPI_Send(buf, count, to_mpi_type(datatype), mpi_dest, mpi_tag, to_mpi_comm(comm.mpi_comm)); 36 //printf("call mpi_send for intracomm, dest = %d, tag = %d\n", dest, tag); 37 } 27 38 //check_sum_send(buf, count, datatype, dest, tag, comm); 28 39 -
XIOS/dev/branch_openmp/extern/src_ep_dev/ep_type.hpp
r1295 r1328 344 344 } 345 345 346 bool is_null(); 347 346 348 }; 347 349 -
XIOS/dev/branch_openmp/inputs/REMAP/iodef.xml
r1220 r1328 51 51 52 52 <file_group id="write_files" > 53 <file id="output_2D" name="output_2D" enabled=".TRUE.">53 <file id="output_2D" name="output_2D" output_freq="1ts" enabled=".TRUE."> 54 54 <field field_ref="src_field_2D" name="field_src" enabled=".TRUE."/> 55 55 <field field_ref="src_field_2D_clone" name="field_src_clone" default_value="100000" enabled=".TRUE."/> … … 57 57 <field field_ref="dst_field_2D" name="field_dst_regular_1" enabled=".TRUE." /> 58 58 <field field_ref="dst_field_2D_regular_pole" name="field_dst_regular_2" enabled=".TRUE."/> 59 <field field_ref="dst_field_2D_clone" name="field_dst_regular_3" detect_missing_value=".false." default_value="100000" enabled=". TRUE."/>59 <field field_ref="dst_field_2D_clone" name="field_dst_regular_3" detect_missing_value=".false." default_value="100000" enabled=".FALSE."/> 60 60 <field field_ref="dst_field_2D_extract" name="field_dst_regular_4" enabled=".TRUE."/> 61 61 </file> -
XIOS/dev/branch_openmp/src/array_new.hpp
r1134 r1328 554 554 TinyVector<int,N_rank> vect; 555 555 size_t ne; 556 556 557 557 ret = buffer.get(numDim); 558 558 ret &= buffer.get(vect.data(), N_rank); -
XIOS/dev/branch_openmp/src/attribute_enum.hpp
r1134 r1328 14 14 namespace xios 15 15 { 16 /// ////////////////////// Déclarations ////////////////////// ///17 /*!18 \class CAttributeEnum19 This class implements the attribute representing enumeration16 /// ////////////////////// Déclarations ////////////////////// /// 17 /*! 18 \class CAttributeEnum 19 This class implements the attribute representing enumeration 20 20 */ 21 template <class T>22 class CAttributeEnum : public CAttribute, public CEnum<T>23 {21 template <class T> 22 class CAttributeEnum : public CAttribute, public CEnum<T> 23 { 24 24 typedef typename T::t_enum T_enum ; 25 25 public : 26 26 27 /// Constructeurs ///28 explicit CAttributeEnum(const StdString & id);29 CAttributeEnum(const StdString & id,30 xios_map<StdString, CAttribute*> & umap);31 CAttributeEnum(const StdString & id, const T_enum & value);32 CAttributeEnum(const StdString & id, const T_enum & value,33 xios_map<StdString, CAttribute*> & umap);27 /// Constructeurs /// 28 explicit CAttributeEnum(const StdString & id); 29 CAttributeEnum(const StdString & id, 30 xios_map<StdString, CAttribute*> & umap); 31 CAttributeEnum(const StdString & id, const T_enum & value); 32 CAttributeEnum(const StdString & id, const T_enum & value, 33 xios_map<StdString, CAttribute*> & umap); 34 34 35 /// Accesseur ///36 T_enum getValue(void) const;37 string getStringValue(void) const;35 /// Accesseur /// 36 T_enum getValue(void) const; 37 string getStringValue(void) const; 38 38 39 39 40 /// Mutateurs /// 41 void setValue(const T_enum & value); 40 /// Mutateurs /// 41 void setValue(const T_enum & value); 42 43 void set(const CAttribute& attr) ; 44 void set(const CAttributeEnum& attr) ; 45 void reset(void); 46 47 void setInheritedValue(const CAttributeEnum& attr ); 48 void setInheritedValue(const CAttribute& attr ); 49 T_enum getInheritedValue(void) const; 50 string getInheritedStringValue(void) const; 51 bool hasInheritedValue(void) const; 52 53 bool isEqual(const CAttributeEnum& attr ); 54 bool isEqual(const CAttribute& attr ); 42 55 43 void set(const CAttribute& attr) ; 44 void set(const CAttributeEnum& attr) ; 45 void reset(void); 56 /// Destructeur /// 57 virtual ~CAttributeEnum(void) { } 46 58 47 void setInheritedValue(const CAttributeEnum& attr ); 48 void setInheritedValue(const CAttribute& attr ); 49 T_enum getInheritedValue(void) const; 50 string getInheritedStringValue(void) const; 51 bool hasInheritedValue(void) const; 59 /// Operateur /// 60 CAttributeEnum& operator=(const T_enum & value); 52 61 53 bool isEqual(const CAttributeEnum& attr ); 54 bool isEqual(const CAttribute& attr ); 62 /// Autre /// 63 virtual StdString toString(void) const { return _toString();} 64 virtual void fromString(const StdString & str) { if (str==resetInheritanceStr) { reset(); _canInherite=false ;} else _fromString(str);} 55 65 56 /// Destructeur /// 57 virtual ~CAttributeEnum(void) { } 66 virtual bool toBuffer (CBufferOut& buffer) const { return _toBuffer(buffer);} 67 virtual bool fromBuffer(CBufferIn& buffer) { return _fromBuffer(buffer); } 68 69 virtual void generateCInterface(ostream& oss,const string& className) ; 70 virtual void generateFortran2003Interface(ostream& oss,const string& className) ; 71 virtual void generateFortranInterfaceDeclaration_(ostream& oss,const string& className) ; 72 virtual void generateFortranInterfaceBody_(ostream& oss,const string& className) ; 73 virtual void generateFortranInterfaceDeclaration(ostream& oss,const string& className) ; 74 virtual void generateFortranInterfaceGetDeclaration_(ostream& oss,const string& className) ; 75 virtual void generateFortranInterfaceGetBody_(ostream& oss,const string& className) ; 76 virtual void generateFortranInterfaceGetDeclaration(ostream& oss,const string& className) ; 58 77 59 /// Operateur /// 60 CAttributeEnum& operator=(const T_enum & value); 61 62 /// Autre /// 63 virtual StdString toString(void) const { return _toString();} 64 virtual void fromString(const StdString & str) { if (str==resetInheritanceStr) { reset(); _canInherite=false ;} else _fromString(str);} 65 66 virtual bool toBuffer (CBufferOut& buffer) const { return _toBuffer(buffer);} 67 virtual bool fromBuffer(CBufferIn& buffer) { return _fromBuffer(buffer); } 68 69 virtual void generateCInterface(ostream& oss,const string& className) ; 70 virtual void generateFortran2003Interface(ostream& oss,const string& className) ; 71 virtual void generateFortranInterfaceDeclaration_(ostream& oss,const string& className) ; 72 virtual void generateFortranInterfaceBody_(ostream& oss,const string& className) ; 73 virtual void generateFortranInterfaceDeclaration(ostream& oss,const string& className) ; 74 virtual void generateFortranInterfaceGetDeclaration_(ostream& oss,const string& className) ; 75 virtual void generateFortranInterfaceGetBody_(ostream& oss,const string& className) ; 76 virtual void generateFortranInterfaceGetDeclaration(ostream& oss,const string& className) ; 77 78 private : 79 StdString _toString(void) const; 80 void _fromString(const StdString & str); 81 bool _toBuffer (CBufferOut& buffer) const; 82 bool _fromBuffer(CBufferIn& buffer) ; 83 CEnum<T> inheritedValue ; 84 }; // class CAttributeEnum 85 78 private : 79 StdString _toString(void) const; 80 void _fromString(const StdString & str); 81 bool _toBuffer (CBufferOut& buffer) const; 82 bool _fromBuffer(CBufferIn& buffer) ; 83 CEnum<T> inheritedValue ; 84 }; // class CAttributeEnum 85 86 86 } // namespace xios 87 87 88 88 #endif // __XIOS_ATTRIBUTE_ENUM__ 89 -
XIOS/dev/branch_openmp/src/attribute_enum_impl.hpp
r1134 r1328 10 10 namespace xios 11 11 { 12 /// ////////////////////// D éfinitions ////////////////////// ///12 /// ////////////////////// Définitions ////////////////////// /// 13 13 template <class T> 14 14 CAttributeEnum<T>::CAttributeEnum(const StdString & id) … … 30 30 umap.insert(umap.end(), std::make_pair(id, this)); 31 31 } 32 32 33 33 template <class T> 34 34 CAttributeEnum<T>::CAttributeEnum … … 40 40 umap.insert(umap.end(), std::make_pair(id, this)); 41 41 } 42 42 43 43 ///-------------------------------------------------------------- 44 44 template <class T> … … 54 54 return CEnum<T>::get(); 55 55 } 56 56 57 57 template <class T> 58 58 string CAttributeEnum<T>::getStringValue(void) const 59 59 { 60 return CEnum<T>::toString(); 61 } 62 60 return CEnum<T>::toString(); 61 } 63 62 64 63 template <class T> … … 71 70 void CAttributeEnum<T>::set(const CAttribute& attr) 72 71 { 73 74 } 75 76 72 this->set(dynamic_cast<const CAttributeEnum<T>& >(attr)); 73 } 74 75 template <class T> 77 76 void CAttributeEnum<T>::set(const CAttributeEnum& attr) 78 77 { 79 80 } 81 78 CEnum<T>::set(attr); 79 } 80 82 81 template <class T> 83 82 void CAttributeEnum<T>::setInheritedValue(const CAttribute& attr) 84 83 { 85 86 } 87 84 this->setInheritedValue(dynamic_cast<const CAttributeEnum<T>& >(attr)); 85 } 86 88 87 template <class T> 89 88 void CAttributeEnum<T>::setInheritedValue(const CAttributeEnum& attr) 90 89 { 91 92 } 93 90 if (this->isEmpty() && _canInherite && attr.hasInheritedValue()) inheritedValue.set(attr.getInheritedValue()); 91 } 92 94 93 template <class T> 95 94 typename T::t_enum CAttributeEnum<T>::getInheritedValue(void) const 96 95 { 97 98 99 } 100 101 template <class T> 102 103 104 105 106 107 108 template <class T> 109 110 111 112 113 114 template <class T> 115 116 117 118 119 120 template <class T> 121 122 123 124 96 if (this->isEmpty()) return inheritedValue.get(); 97 else return getValue(); 98 } 99 100 template <class T> 101 string CAttributeEnum<T>::getInheritedStringValue(void) const 102 { 103 if (this->isEmpty()) return inheritedValue.toString(); 104 else return CEnum<T>::toString();; 105 } 106 107 template <class T> 108 bool CAttributeEnum<T>::hasInheritedValue(void) const 109 { 110 return !this->isEmpty() || !inheritedValue.isEmpty(); 111 } 112 113 template <class T> 114 bool CAttributeEnum<T>::isEqual(const CAttribute& attr) 115 { 116 return (this->isEqual(dynamic_cast<const CAttributeEnum<T>& >(attr))); 117 } 118 119 template <class T> 120 bool CAttributeEnum<T>::isEqual(const CAttributeEnum& attr) 121 { 122 return ((dynamic_cast<const CEnum<T>& >(*this)) == (dynamic_cast<const CEnum<T>& >(attr))); 123 } 125 124 126 125 //--------------------------------------------------------------- 127 126 128 127 template <class T> 129 130 131 132 133 128 CAttributeEnum<T>& CAttributeEnum<T>::operator=(const T_enum & value) 129 { 130 this->setValue(value); 131 return *this; 132 } 134 133 135 134 //--------------------------------------------------------------- 136 135 137 136 template <class T> 138 139 140 141 142 143 144 145 146 template <class T> 147 148 149 150 151 152 template <class T> 153 154 155 156 157 158 template <class T> 159 160 161 162 163 164 template <typename T> 165 166 167 168 169 170 template <typename T> 171 172 173 174 175 176 template <typename T> 177 178 179 180 181 182 template <typename T> 183 184 185 186 187 188 template <typename T> 189 190 191 192 193 194 template <typename T> 195 196 197 198 199 200 template <typename T> 201 202 203 204 205 206 template <typename T> 207 208 209 210 137 StdString CAttributeEnum<T>::_toString(void) const 138 { 139 StdOStringStream oss; 140 if (!CEnum<T>::isEmpty() && this->hasId()) 141 oss << this->getName() << "=\"" << CEnum<T>::toString() << "\""; 142 return (oss.str()); 143 } 144 145 template <class T> 146 void CAttributeEnum<T>::_fromString(const StdString & str) 147 { 148 CEnum<T>::fromString(str); 149 } 150 151 template <class T> 152 bool CAttributeEnum<T>::_toBuffer (CBufferOut& buffer) const 153 { 154 return CEnum<T>::toBuffer(buffer); 155 } 156 157 template <class T> 158 bool CAttributeEnum<T>::_fromBuffer(CBufferIn& buffer) 159 { 160 return CEnum<T>::fromBuffer(buffer); 161 } 162 163 template <typename T> 164 void CAttributeEnum<T>::generateCInterface(ostream& oss,const string& className) 165 { 166 CInterface::AttributeCInterface<CEnumBase>(oss, className, this->getName()); 167 } 168 169 template <typename T> 170 void CAttributeEnum<T>::generateFortran2003Interface(ostream& oss,const string& className) 171 { 172 CInterface::AttributeFortran2003Interface<string>(oss, className, this->getName()); 173 } 174 175 template <typename T> 176 void CAttributeEnum<T>::generateFortranInterfaceDeclaration_(ostream& oss,const string& className) 177 { 178 CInterface::AttributeFortranInterfaceDeclaration<string>(oss, className, this->getName()+"_"); 179 } 180 181 template <typename T> 182 void CAttributeEnum<T>::generateFortranInterfaceBody_(ostream& oss,const string& className) 183 { 184 CInterface::AttributeFortranInterfaceBody<string>(oss, className, this->getName()); 185 } 186 187 template <typename T> 188 void CAttributeEnum<T>::generateFortranInterfaceDeclaration(ostream& oss,const string& className) 189 { 190 CInterface::AttributeFortranInterfaceDeclaration<string>(oss, className, this->getName()); 191 } 192 193 template <typename T> 194 void CAttributeEnum<T>::generateFortranInterfaceGetDeclaration_(ostream& oss,const string& className) 195 { 196 CInterface::AttributeFortranInterfaceGetDeclaration<string>(oss, className, this->getName()+"_"); 197 } 198 199 template <typename T> 200 void CAttributeEnum<T>::generateFortranInterfaceGetBody_(ostream& oss,const string& className) 201 { 202 CInterface::AttributeFortranInterfaceGetBody<string>(oss, className, this->getName()); 203 } 204 205 template <typename T> 206 void CAttributeEnum<T>::generateFortranInterfaceGetDeclaration(ostream& oss,const string& className) 207 { 208 CInterface::AttributeFortranInterfaceGetDeclaration<string>(oss, className, this->getName()); 209 } 211 210 } // namespace xios 212 211 213 212 #endif // __XIOS_ATTRIBUTE_ENUM_IMPL_HPP__ 214 -
XIOS/dev/branch_openmp/src/attribute_map.hpp
r1134 r1328 76 76 /// Propriété statique /// 77 77 static CAttributeMap * Current; 78 #pragma omp threadprivate (Current)79 78 80 79 }; // class CAttributeMap -
XIOS/dev/branch_openmp/src/buffer_client.cpp
r1205 r1328 7 7 #include "mpi.hpp" 8 8 #include "tracer.hpp" 9 10 11 using namespace ep_lib; 9 12 10 13 namespace xios … … 27 30 buffer[1] = new char[bufferSize]; 28 31 retBuffer = new CBufferOut(buffer[current], bufferSize); 29 #pragma omp critical (_output)30 32 info(10) << "CClientBuffer: allocated 2 x " << bufferSize << " bytes for server " << serverRank << " with a maximum of " << maxBufferedEvents << " buffered events" << endl; 31 33 } -
XIOS/dev/branch_openmp/src/buffer_client.hpp
r1205 r1328 6 6 #include "mpi.hpp" 7 7 #include "cxios.hpp" 8 #ifdef _usingEP9 #include "ep_declaration.hpp"10 #endif11 8 12 9 namespace xios … … 16 13 public: 17 14 static size_t maxRequestSize; 18 #pragma omp threadprivate(maxRequestSize)19 15 20 CClientBuffer( MPI_Comm intercomm, int serverRank, StdSize bufferSize, StdSize estimatedMaxEventSize, StdSize maxBufferedEvents);16 CClientBuffer(ep_lib::MPI_Comm intercomm, int serverRank, StdSize bufferSize, StdSize estimatedMaxEventSize, StdSize maxBufferedEvents); 21 17 ~CClientBuffer(); 22 18 … … 40 36 bool pending; 41 37 42 MPI_Request request;38 ep_lib::MPI_Request request; 43 39 44 40 CBufferOut* retBuffer; 45 const MPI_Comm interComm;41 const ep_lib::MPI_Comm interComm; 46 42 }; 47 43 } -
XIOS/dev/branch_openmp/src/buffer_server.hpp
r1134 r1328 4 4 #include "xios_spl.hpp" 5 5 #include "buffer.hpp" 6 #include "mpi _std.hpp"6 #include "mpi.hpp" 7 7 #include "cxios.hpp" 8 8 -
XIOS/dev/branch_openmp/src/calendar.cpp
r1134 r1328 117 117 const CDate& CCalendar::update(int step) 118 118 { 119 #pragma omp critical (_output) 120 info(80)<< "update step : " << step << " timestep " << this->timestep << std::endl; 119 info(20) << "update step : " << step << " timestep " << this->timestep << std::endl; 121 120 return (this->currentDate = this->getInitDate() + step * this->timestep); 122 121 } -
XIOS/dev/branch_openmp/src/client.cpp
r1205 r1328 11 11 #include "timer.hpp" 12 12 #include "buffer_client.hpp" 13 #include "log.hpp" 14 13 using namespace ep_lib; 15 14 16 15 namespace xios 17 16 { 18 extern int test_omp_rank;19 #pragma omp threadprivate(test_omp_rank)20 17 21 18 MPI_Comm CClient::intraComm ; 22 19 MPI_Comm CClient::interComm ; 20 //std::list<MPI_Comm> CClient::contextInterComms; 23 21 std::list<MPI_Comm> *CClient::contextInterComms_ptr = 0; 24 22 int CClient::serverLeader ; … … 28 26 StdOFStream CClient::m_errorStream; 29 27 30 StdOFStream CClient::array_infoStream[16]; 31 32 void CClient::initialize(const string& codeId, MPI_Comm& localComm, MPI_Comm& returnComm) 28 void CClient::initialize(const string& codeId,MPI_Comm& localComm,MPI_Comm& returnComm) 33 29 { 34 30 int initialized ; … … 41 37 { 42 38 // localComm doesn't given 43 44 39 if (localComm == MPI_COMM_NULL) 45 40 { 46 41 if (!is_MPI_Initialized) 47 42 { 48 //MPI_Init(NULL, NULL); 49 int return_level; 50 MPI_Init_thread(NULL, NULL, 3, &return_level); 51 assert(return_level == 3); 43 MPI_Init(NULL, NULL); 52 44 } 53 45 CTimer::get("XIOS").resume() ; … … 61 53 int myColor ; 62 54 int i,c ; 63 64 MPI_Comm_size(CXios::globalComm,&size); 55 MPI_Comm newComm ; 56 57 MPI_Comm_size(CXios::globalComm,&size) ; 65 58 MPI_Comm_rank(CXios::globalComm,&rank); 66 67 59 68 60 hashAll=new unsigned long[size] ; … … 106 98 MPI_Comm_size(intraComm,&intraCommSize) ; 107 99 MPI_Comm_rank(intraComm,&intraCommRank) ; 108 109 #pragma omp critical(_output) 110 { 111 info(10)<<"intercommCreate::client "<<test_omp_rank<< " "<< &test_omp_rank <<" intraCommSize : "<<intraCommSize 112 <<" intraCommRank :"<<intraCommRank<<" serverLeader "<< serverLeader 113 <<" globalComm : "<< &(CXios::globalComm) << endl ; 114 } 115 116 117 //test_sendrecv(CXios::globalComm); 100 info(50)<<"intercommCreate::client "<<rank<<" intraCommSize : "<<intraCommSize 101 <<" intraCommRank :"<<intraCommRank<<" clientLeader "<< serverLeader<<endl ; 118 102 MPI_Intercomm_create(intraComm,0,CXios::globalComm,serverLeader,0,&interComm) ; 119 120 103 } 121 104 else … … 140 123 } 141 124 // using OASIS 142 else125 /* else 143 126 { 144 127 // localComm doesn't given … … 165 148 else MPI_Comm_dup(intraComm,&interComm) ; 166 149 } 167 150 */ 168 151 MPI_Comm_dup(intraComm,&returnComm) ; 169 170 152 } 171 153 … … 174 156 { 175 157 CContext::setCurrent(id) ; 176 CContext* context = CContext::create(id); 177 178 int tmp_rank; 179 MPI_Comm_rank(contextComm,&tmp_rank) ; 180 158 CContext* context=CContext::create(id); 181 159 StdString idServer(id); 182 160 idServer += "_server"; … … 185 163 { 186 164 int size,rank,globalRank ; 187 //size_t message_size ;188 //int leaderRank ;165 size_t message_size ; 166 int leaderRank ; 189 167 MPI_Comm contextInterComm ; 190 168 … … 197 175 CMessage msg ; 198 176 msg<<idServer<<size<<globalRank ; 199 177 // msg<<id<<size<<globalRank ; 200 178 201 179 int messageSize=msg.size() ; … … 208 186 209 187 MPI_Intercomm_create(contextComm,0,CXios::globalComm,serverLeader,10+globalRank,&contextInterComm) ; 210 211 #pragma omp critical(_output) 212 info(10)<<" RANK "<< tmp_rank<<" Register new Context : "<<id<<endl ; 213 188 info(10)<<"Register new Context : "<<id<<endl ; 214 189 215 190 MPI_Comm inter ; … … 217 192 MPI_Barrier(inter) ; 218 193 219 220 194 context->initClient(contextComm,contextInterComm) ; 221 195 222 196 //contextInterComms.push_back(contextInterComm); 223 197 if(contextInterComms_ptr == NULL) contextInterComms_ptr = new std::list<MPI_Comm>; 224 198 contextInterComms_ptr->push_back(contextInterComm); 225 226 199 MPI_Comm_free(&inter); 227 200 } … … 240 213 // Finally, we should return current context to context client 241 214 CContext::setCurrent(id); 242 215 216 //contextInterComms.push_back(contextInterComm); 243 217 if(contextInterComms_ptr == NULL) contextInterComms_ptr = new std::list<MPI_Comm>; 244 218 contextInterComms_ptr->push_back(contextInterComm); 245 246 219 } 247 220 } … … 253 226 254 227 MPI_Comm_rank(intraComm,&rank) ; 255 228 256 229 if (!CXios::isServer) 257 230 { … … 263 236 } 264 237 265 for (std::list<MPI_Comm>::iterator it = contextInterComms_ptr->begin(); it != contextInterComms_ptr->end(); ++it) 238 //for (std::list<MPI_Comm>::iterator it = contextInterComms.begin(); it != contextInterComms.end(); it++) 239 for (std::list<MPI_Comm>::iterator it = contextInterComms_ptr->begin(); it != contextInterComms_ptr->end(); it++) 266 240 MPI_Comm_free(&(*it)); 267 268 241 MPI_Comm_free(&interComm); 269 242 MPI_Comm_free(&intraComm); … … 274 247 if (!is_MPI_Initialized) 275 248 { 276 if (CXios::usingOasis) oasis_finalize(); 277 else MPI_Finalize(); 249 //if (CXios::usingOasis) oasis_finalize(); 250 //else 251 MPI_Finalize() ; 278 252 } 279 253 280 #pragma omp critical (_output) 281 info(20) << "Client "<<rank<<" : Client side context is finalized "<< endl ; 282 283 /*#pragma omp critical (_output) 284 { 254 info(20) << "Client side context is finalized"<<endl ; 285 255 report(0) <<" Performance report : Whole time from XIOS init and finalize: "<< CTimer::get("XIOS init/finalize").getCumulatedTime()<<" s"<<endl ; 286 256 report(0) <<" Performance report : total time spent for XIOS : "<< CTimer::get("XIOS").getCumulatedTime()<<" s"<<endl ; … … 292 262 report(0)<< " Memory report : increasing it by a factor will increase performance, depending of the volume of data wrote in file at each time step of the file"<<endl ; 293 263 report(100)<<CTimer::getAllCumulatedTime()<<endl ; 294 }*/295 296 264 } 297 265 … … 322 290 323 291 fileNameClient << fileName << "_" << std::setfill('0') << std::setw(numDigit) << getRank() << ext; 324 325 292 fb->open(fileNameClient.str().c_str(), std::ios::out); 326 293 if (!fb->is_open()) 327 294 ERROR("void CClient::openStream(const StdString& fileName, const StdString& ext, std::filebuf* fb)", 328 << std::endl << "Can not open <" << fileNameClient << "> file to write the client log(s).");295 << std::endl << "Can not open <" << fileNameClient << "> file to write the client log(s)."); 329 296 } 330 297 … … 337 304 void CClient::openInfoStream(const StdString& fileName) 338 305 { 339 //std::filebuf* fb = m_infoStream.rdbuf(); 340 341 info_FB[omp_get_thread_num()] = array_infoStream[omp_get_thread_num()].rdbuf(); 342 343 openStream(fileName, ".out", info_FB[omp_get_thread_num()]); 344 345 info.write2File(info_FB[omp_get_thread_num()]); 346 report.write2File(info_FB[omp_get_thread_num()]); 347 306 std::filebuf* fb = m_infoStream.rdbuf(); 307 openStream(fileName, ".out", fb); 308 309 info.write2File(fb); 310 report.write2File(fb); 348 311 } 349 312 -
XIOS/dev/branch_openmp/src/client.hpp
r1164 r1328 10 10 { 11 11 public: 12 static void initialize(const string& codeId, MPI_Comm& localComm,MPI_Comm& returnComm);12 static void initialize(const string& codeId, ep_lib::MPI_Comm& localComm, ep_lib::MPI_Comm& returnComm); 13 13 static void finalize(void); 14 14 static void registerContext(const string& id, ep_lib::MPI_Comm contextComm); 15 15 16 static MPI_Comm intraComm; 17 #pragma omp threadprivate(intraComm) 18 19 static MPI_Comm interComm; 20 #pragma omp threadprivate(interComm) 21 16 static ep_lib::MPI_Comm intraComm; 17 static ep_lib::MPI_Comm interComm; 22 18 //static std::list<MPI_Comm> contextInterComms; 23 24 static std::list<MPI_Comm> * contextInterComms_ptr; 25 #pragma omp threadprivate(contextInterComms_ptr) 26 19 static std::list<ep_lib::MPI_Comm> *contextInterComms_ptr; 27 20 static int serverLeader; 28 #pragma omp threadprivate(serverLeader)29 30 21 static bool is_MPI_Initialized ; 31 #pragma omp threadprivate(is_MPI_Initialized)32 22 33 23 //! Get rank of the current process … … 50 40 protected: 51 41 static int rank; 52 #pragma omp threadprivate(rank)53 54 42 static StdOFStream m_infoStream; 55 #pragma omp threadprivate(m_infoStream)56 57 43 static StdOFStream m_errorStream; 58 #pragma omp threadprivate(m_errorStream)59 60 static StdOFStream array_infoStream[16];61 44 62 45 static void openStream(const StdString& fileName, const StdString& ext, std::filebuf* fb); -
XIOS/dev/branch_openmp/src/client_client_dht_template.hpp
r1134 r1328 13 13 #include "xios_spl.hpp" 14 14 #include "array_new.hpp" 15 #include "mpi _std.hpp"15 #include "mpi.hpp" 16 16 #include "policy.hpp" 17 17 #include <boost/unordered_map.hpp> … … 87 87 const ep_lib::MPI_Comm& clientIntraComm, 88 88 std::vector<ep_lib::MPI_Request>& requestSendInfo); 89 void sendInfoToClients(int clientDestRank, unsigned char* info, int infoSize, 90 const ep_lib::MPI_Comm& clientIntraComm, 91 ep_lib::MPI_Request* requestSendInfo); 89 92 90 93 void recvInfoFromClients(int clientSrcRank, unsigned char* info, int infoSize, 91 94 const ep_lib::MPI_Comm& clientIntraComm, 92 95 std::vector<ep_lib::MPI_Request>& requestRecvInfo); 96 void recvInfoFromClients(int clientSrcRank, unsigned char* info, int infoSize, 97 const ep_lib::MPI_Comm& clientIntraComm, 98 ep_lib::MPI_Request* requestRecvInfo); 99 93 100 94 101 // Send global index to clients … … 96 103 const ep_lib::MPI_Comm& clientIntraComm, 97 104 std::vector<ep_lib::MPI_Request>& requestSendIndexGlobal); 105 void sendIndexToClients(int clientDestRank, size_t* indices, size_t indiceSize, 106 const ep_lib::MPI_Comm& clientIntraComm, 107 ep_lib::MPI_Request* requestSendIndexGlobal); 98 108 99 109 void recvIndexFromClients(int clientSrcRank, size_t* indices, size_t indiceSize, 100 110 const ep_lib::MPI_Comm& clientIntraComm, 101 111 std::vector<ep_lib::MPI_Request>& requestRecvIndex); 112 void recvIndexFromClients(int clientSrcRank, size_t* indices, size_t indiceSize, 113 const ep_lib::MPI_Comm& clientIntraComm, 114 ep_lib::MPI_Request* requestRecvIndex); 102 115 103 116 void sendRecvOnReturn(const std::vector<int>& sendNbRank, std::vector<int>& sendNbElements, -
XIOS/dev/branch_openmp/src/client_client_dht_template_impl.hpp
r1209 r1328 10 10 #include "utils.hpp" 11 11 #include "mpi_tag.hpp" 12 #ifdef _usingEP13 #include "ep_declaration.hpp"14 #include "ep_lib.hpp"15 #endif16 17 12 18 13 namespace xios … … 22 17 : H(clientIntraComm), index2InfoMapping_(), indexToInfoMappingLevel_(), nbClient_(0) 23 18 { 24 MPI_Comm_size(clientIntraComm, &nbClient_);19 ep_lib::MPI_Comm_size(clientIntraComm, &nbClient_); 25 20 this->computeMPICommLevel(); 26 21 int nbLvl = this->getNbLevel(); … … 42 37 : H(clientIntraComm), index2InfoMapping_(), indexToInfoMappingLevel_(), nbClient_(0) 43 38 { 44 MPI_Comm_size(clientIntraComm, &nbClient_);39 ep_lib::MPI_Comm_size(clientIntraComm, &nbClient_); 45 40 this->computeMPICommLevel(); 46 41 int nbLvl = this->getNbLevel(); … … 67 62 : H(clientIntraComm), index2InfoMapping_(), indexToInfoMappingLevel_(), nbClient_(0) 68 63 { 69 MPI_Comm_size(clientIntraComm, &nbClient_);64 ep_lib::MPI_Comm_size(clientIntraComm, &nbClient_); 70 65 this->computeMPICommLevel(); 71 66 int nbLvl = this->getNbLevel(); … … 104 99 { 105 100 int clientRank; 106 MPI_Comm_rank(commLevel,&clientRank); 107 //ep_lib::MPI_Barrier(commLevel); 101 ep_lib::MPI_Comm_rank(commLevel,&clientRank); 108 102 int groupRankBegin = this->getGroupBegin()[level]; 109 103 int nbClient = this->getNbInGroup()[level]; … … 175 169 recvIndexBuff = new unsigned long[recvNbIndexCount]; 176 170 177 int request_size = 0; 178 179 int currentIndex = 0; 180 int nbRecvClient = recvRankClient.size(); 181 182 int position = 0; 183 184 for (int idx = 0; idx < nbRecvClient; ++idx) 171 int request_size = 0; 172 for (int idx = 0; idx < recvRankClient.size(); ++idx) 185 173 { 186 174 if (0 != recvNbIndexClientCount[idx]) 187 { 188 request_size++; 189 } 175 request_size ++; 190 176 } 191 177 192 178 request_size += client2ClientIndex.size(); 193 179 194 195 180 std::vector<ep_lib::MPI_Request> request(request_size); 196 181 197 182 std::vector<int>::iterator itbRecvIndex = recvRankClient.begin(), itRecvIndex, 198 183 iteRecvIndex = recvRankClient.end(), 199 184 itbRecvNbIndex = recvNbIndexClientCount.begin(), 200 185 itRecvNbIndex; 201 202 186 int currentIndex = 0; 187 int nbRecvClient = recvRankClient.size(); 188 int request_position = 0; 189 for (int idx = 0; idx < nbRecvClient; ++idx) 190 { 191 if (0 != recvNbIndexClientCount[idx]) 192 recvIndexFromClients(recvRankClient[idx], recvIndexBuff+currentIndex, recvNbIndexClientCount[idx], commLevel, &request[request_position++]); 193 //recvIndexFromClients(recvRankClient[idx], recvIndexBuff+currentIndex, recvNbIndexClientCount[idx], commLevel, request); 194 currentIndex += recvNbIndexClientCount[idx]; 195 } 196 203 197 boost::unordered_map<int, size_t* >::iterator itbIndex = client2ClientIndex.begin(), itIndex, 204 198 iteIndex = client2ClientIndex.end(); 205 199 for (itIndex = itbIndex; itIndex != iteIndex; ++itIndex) 206 { 207 MPI_Isend(itIndex->second, sendNbIndexBuff[itIndex->first-groupRankBegin], MPI_UNSIGNED_LONG, itIndex->first, MPI_DHT_INDEX, commLevel, &request[position]); 208 position++; 200 sendIndexToClients(itIndex->first, (itIndex->second), sendNbIndexBuff[itIndex->first-groupRankBegin], commLevel, &request[request_position++]); 209 201 //sendIndexToClients(itIndex->first, (itIndex->second), sendNbIndexBuff[itIndex->first-groupRankBegin], commLevel, request); 210 } 211 212 for (int idx = 0; idx < nbRecvClient; ++idx) 213 { 214 if (0 != recvNbIndexClientCount[idx]) 215 { 216 MPI_Irecv(recvIndexBuff+currentIndex, recvNbIndexClientCount[idx], MPI_UNSIGNED_LONG, 217 recvRankClient[idx], MPI_DHT_INDEX, commLevel, &request[position]); 218 position++; 219 //recvIndexFromClients(recvRankClient[idx], recvIndexBuff+currentIndex, recvNbIndexClientCount[idx], commLevel, request); 220 } 221 currentIndex += recvNbIndexClientCount[idx]; 222 } 223 224 225 std::vector<ep_lib::MPI_Status> status(request_size); 226 MPI_Waitall(request.size(), &request[0], &status[0]); 227 202 203 std::vector<ep_lib::MPI_Status> status(request.size()); 204 ep_lib::MPI_Waitall(request.size(), &request[0], &status[0]); 228 205 229 206 CArray<size_t,1>* tmpGlobalIndex; … … 238 215 --level; 239 216 computeIndexInfoMappingLevel(*tmpGlobalIndex, this->internalComm_, level); 240 241 217 } 242 218 else // Now, we are in the last level where necessary mappings are. … … 279 255 } 280 256 281 request_size =0;257 int requestOnReturn_size=0; 282 258 for (int idx = 0; idx < recvRankOnReturn.size(); ++idx) 283 259 { 284 260 if (0 != recvNbIndexOnReturn[idx]) 285 261 { 286 request _size += 2;262 requestOnReturn_size += 2; 287 263 } 288 264 } … … 292 268 if (0 != sendNbIndexOnReturn[idx]) 293 269 { 294 request_size += 2; 295 } 296 } 297 298 std::vector<ep_lib::MPI_Request> requestOnReturn(request_size); 270 requestOnReturn_size += 2; 271 } 272 } 273 274 int requestOnReturn_position=0; 275 276 std::vector<ep_lib::MPI_Request> requestOnReturn(requestOnReturn_size); 299 277 currentIndex = 0; 300 position = 0;301 278 for (int idx = 0; idx < recvRankOnReturn.size(); ++idx) 302 279 { … … 304 281 { 305 282 //recvIndexFromClients(recvRankOnReturn[idx], recvIndexBuffOnReturn+currentIndex, recvNbIndexOnReturn[idx], commLevel, requestOnReturn); 306 MPI_Irecv(recvIndexBuffOnReturn+currentIndex, recvNbIndexOnReturn[idx], MPI_UNSIGNED_LONG,307 recvRankOnReturn[idx], MPI_DHT_INDEX, commLevel, &requestOnReturn[position]);308 position++;309 283 //recvInfoFromClients(recvRankOnReturn[idx], 310 284 // recvInfoBuffOnReturn+currentIndex*ProcessDHTElement<InfoType>::typeSize(), 311 285 // recvNbIndexOnReturn[idx]*ProcessDHTElement<InfoType>::typeSize(), 312 286 // commLevel, requestOnReturn); 313 MPI_Irecv(recvInfoBuffOnReturn+currentIndex*ProcessDHTElement<InfoType>::typeSize(), 314 recvNbIndexOnReturn[idx]*ProcessDHTElement<InfoType>::typeSize(), MPI_CHAR, 315 recvRankOnReturn[idx], MPI_DHT_INFO, commLevel, &requestOnReturn[position]); 316 position++; 287 recvIndexFromClients(recvRankOnReturn[idx], recvIndexBuffOnReturn+currentIndex, recvNbIndexOnReturn[idx], commLevel, &requestOnReturn[requestOnReturn_position++]); 288 recvInfoFromClients(recvRankOnReturn[idx], 289 recvInfoBuffOnReturn+currentIndex*ProcessDHTElement<InfoType>::typeSize(), 290 recvNbIndexOnReturn[idx]*ProcessDHTElement<InfoType>::typeSize(), 291 commLevel, &requestOnReturn[requestOnReturn_position++]); 317 292 } 318 293 currentIndex += recvNbIndexOnReturn[idx]; … … 349 324 //sendIndexToClients(rank, client2ClientIndexOnReturn[rank], 350 325 // sendNbIndexOnReturn[idx], commLevel, requestOnReturn); 351 MPI_Isend(client2ClientIndexOnReturn[rank], sendNbIndexOnReturn[idx], MPI_UNSIGNED_LONG,352 rank, MPI_DHT_INDEX, commLevel, &requestOnReturn[position]);353 position++;354 326 //sendInfoToClients(rank, client2ClientInfoOnReturn[rank], 355 327 // sendNbIndexOnReturn[idx]*ProcessDHTElement<InfoType>::typeSize(), commLevel, requestOnReturn); 356 MPI_Isend(client2ClientInfoOnReturn[rank], sendNbIndexOnReturn[idx]*ProcessDHTElement<InfoType>::typeSize(), MPI_CHAR, 357 rank, MPI_DHT_INFO, commLevel, &requestOnReturn[position]); 358 position++; 328 sendIndexToClients(rank, client2ClientIndexOnReturn[rank], 329 sendNbIndexOnReturn[idx], commLevel, &requestOnReturn[requestOnReturn_position++]); 330 sendInfoToClients(rank, client2ClientInfoOnReturn[rank], 331 sendNbIndexOnReturn[idx]*ProcessDHTElement<InfoType>::typeSize(), commLevel, &requestOnReturn[requestOnReturn_position++]); 332 359 333 } 360 334 currentIndex += recvNbIndexClientCount[idx]; … … 362 336 363 337 std::vector<ep_lib::MPI_Status> statusOnReturn(requestOnReturn.size()); 364 MPI_Waitall(requestOnReturn.size(), &requestOnReturn[0], &statusOnReturn[0]);338 ep_lib::MPI_Waitall(requestOnReturn.size(), &requestOnReturn[0], &statusOnReturn[0]); 365 339 366 340 Index2VectorInfoTypeMap indexToInfoMapping; … … 432 406 { 433 407 int clientRank; 434 MPI_Comm_rank(commLevel,&clientRank);408 ep_lib::MPI_Comm_rank(commLevel,&clientRank); 435 409 computeSendRecvRank(level, clientRank); 436 //ep_lib::MPI_Barrier(commLevel);437 410 438 411 int groupRankBegin = this->getGroupBegin()[level]; … … 481 454 { 482 455 client2ClientIndex[indexClient + groupRankBegin][sendNbIndexBuff[indexClient]] = it->first;; 456 // ProcessDHTElement<InfoType>::packElement(it->second, client2ClientInfo[indexClient + groupRankBegin], sendNbInfo[indexClient]); 483 457 ProcessDHTElement<InfoType>::packElement(infoTmp[idx], client2ClientInfo[indexClient + groupRankBegin], sendNbInfo[indexClient]); 484 458 ++sendNbIndexBuff[indexClient]; … … 494 468 int recvNbIndexCount = 0; 495 469 for (int idx = 0; idx < recvNbIndexClientCount.size(); ++idx) 496 {497 470 recvNbIndexCount += recvNbIndexClientCount[idx]; 498 }499 471 500 472 unsigned long* recvIndexBuff; … … 509 481 // it will send a message to the correct clients. 510 482 // Contents of the message are index and its corresponding informatioin 511 int request_size = 0; 483 int request_size = 0; 484 for (int idx = 0; idx < recvRankClient.size(); ++idx) 485 { 486 if (0 != recvNbIndexClientCount[idx]) 487 { 488 request_size += 2; 489 } 490 } 491 492 request_size += client2ClientIndex.size(); 493 request_size += client2ClientInfo.size(); 494 495 std::vector<ep_lib::MPI_Request> request(request_size); 512 496 int currentIndex = 0; 513 497 int nbRecvClient = recvRankClient.size(); 514 int current_pos = 0; 515 498 int request_position=0; 516 499 for (int idx = 0; idx < nbRecvClient; ++idx) 517 500 { 518 501 if (0 != recvNbIndexClientCount[idx]) 519 502 { 520 request_size += 2; 521 } 522 //currentIndex += recvNbIndexClientCount[idx]; 523 } 524 525 request_size += client2ClientIndex.size(); 526 request_size += client2ClientInfo.size(); 527 528 529 530 std::vector<ep_lib::MPI_Request> request(request_size); 531 532 //unsigned long* tmp_send_buf_long[client2ClientIndex.size()]; 533 //unsigned char* tmp_send_buf_char[client2ClientInfo.size()]; 534 535 int info_position = 0; 536 int index_position = 0; 537 503 //recvIndexFromClients(recvRankClient[idx], recvIndexBuff+currentIndex, recvNbIndexClientCount[idx], commLevel, request); 504 //recvInfoFromClients(recvRankClient[idx], 505 // recvInfoBuff+currentIndex*ProcessDHTElement<InfoType>::typeSize(), 506 // recvNbIndexClientCount[idx]*ProcessDHTElement<InfoType>::typeSize(), 507 // commLevel, request); 508 recvIndexFromClients(recvRankClient[idx], recvIndexBuff+currentIndex, recvNbIndexClientCount[idx], commLevel, &request[request_position++]); 509 recvInfoFromClients(recvRankClient[idx], 510 recvInfoBuff+currentIndex*ProcessDHTElement<InfoType>::typeSize(), 511 recvNbIndexClientCount[idx]*ProcessDHTElement<InfoType>::typeSize(), 512 commLevel, &request[request_position++]); 513 } 514 currentIndex += recvNbIndexClientCount[idx]; 515 } 538 516 539 517 boost::unordered_map<int, size_t* >::iterator itbIndex = client2ClientIndex.begin(), itIndex, 540 518 iteIndex = client2ClientIndex.end(); 541 519 for (itIndex = itbIndex; itIndex != iteIndex; ++itIndex) 542 {520 sendIndexToClients(itIndex->first, itIndex->second, sendNbIndexBuff[itIndex->first-groupRankBegin], commLevel, &request[request_position++]); 543 521 //sendIndexToClients(itIndex->first, itIndex->second, sendNbIndexBuff[itIndex->first-groupRankBegin], commLevel, request); 544 545 //tmp_send_buf_long[index_position] = new unsigned long[sendNbIndexBuff[itIndex->first-groupRankBegin]];546 //for(int i=0; i<sendNbIndexBuff[itIndex->first-groupRankBegin]; i++)547 //{548 // tmp_send_buf_long[index_position][i] = (static_cast<unsigned long * >(itIndex->second))[i];549 //}550 //MPI_Isend(tmp_send_buf_long[current_pos], sendNbIndexBuff[itIndex->first-groupRankBegin], MPI_UNSIGNED_LONG,551 MPI_Isend(itIndex->second, sendNbIndexBuff[itIndex->first-groupRankBegin], MPI_UNSIGNED_LONG,552 itIndex->first, MPI_DHT_INDEX, commLevel, &request[current_pos]);553 current_pos++;554 index_position++;555 556 }557 558 522 boost::unordered_map<int, unsigned char*>::iterator itbInfo = client2ClientInfo.begin(), itInfo, 559 523 iteInfo = client2ClientInfo.end(); 560 524 for (itInfo = itbInfo; itInfo != iteInfo; ++itInfo) 561 {525 sendInfoToClients(itInfo->first, itInfo->second, sendNbInfo[itInfo->first-groupRankBegin], commLevel, &request[request_position++]); 562 526 //sendInfoToClients(itInfo->first, itInfo->second, sendNbInfo[itInfo->first-groupRankBegin], commLevel, request); 563 527 564 //tmp_send_buf_char[info_position] = new unsigned char[sendNbInfo[itInfo->first-groupRankBegin]];565 //for(int i=0; i<sendNbInfo[itInfo->first-groupRankBegin]; i++)566 //{567 // tmp_send_buf_char[info_position][i] = (static_cast<unsigned char * >(itInfo->second))[i];568 //}569 570 MPI_Isend(itInfo->second, sendNbInfo[itInfo->first-groupRankBegin], MPI_CHAR,571 itInfo->first, MPI_DHT_INFO, commLevel, &request[current_pos]);572 573 current_pos++;574 info_position++;575 }576 577 for (int idx = 0; idx < nbRecvClient; ++idx)578 {579 if (0 != recvNbIndexClientCount[idx])580 {581 //recvIndexFromClients(recvRankClient[idx], recvIndexBuff+currentIndex, recvNbIndexClientCount[idx], commLevel, request);582 MPI_Irecv(recvIndexBuff+currentIndex, recvNbIndexClientCount[idx], MPI_UNSIGNED_LONG,583 recvRankClient[idx], MPI_DHT_INDEX, commLevel, &request[current_pos]);584 current_pos++;585 586 587 MPI_Irecv(recvInfoBuff+currentIndex*ProcessDHTElement<InfoType>::typeSize(),588 recvNbIndexClientCount[idx]*ProcessDHTElement<InfoType>::typeSize(),589 MPI_CHAR, recvRankClient[idx], MPI_DHT_INFO, commLevel, &request[current_pos]);590 591 current_pos++;592 593 594 595 // recvInfoFromClients(recvRankClient[idx],596 // recvInfoBuff+currentIndex*ProcessDHTElement<InfoType>::typeSize(),597 // recvNbIndexClientCount[idx]*ProcessDHTElement<InfoType>::typeSize(),598 // commLevel, request);599 }600 currentIndex += recvNbIndexClientCount[idx];601 }602 603 528 std::vector<ep_lib::MPI_Status> status(request.size()); 604 605 MPI_Waitall(request.size(), &request[0], &status[0]); 606 607 608 //for(int i=0; i<client2ClientInfo.size(); i++) 609 // delete[] tmp_send_buf_char[i]; 610 611 612 613 //for(int i=0; i<client2ClientIndex.size(); i++) 614 // delete[] tmp_send_buf_long[i]; 615 529 ep_lib::MPI_Waitall(request.size(), &request[0], &status[0]); 616 530 617 531 Index2VectorInfoTypeMap indexToInfoMapping; … … 654 568 else 655 569 index2InfoMapping_.swap(indexToInfoMapping); 656 657 570 } 658 571 … … 670 583 std::vector<ep_lib::MPI_Request>& requestSendIndex) 671 584 { 672 printf("should not call this function sendIndexToClients");673 585 ep_lib::MPI_Request request; 674 586 requestSendIndex.push_back(request); 675 MPI_Isend(indices, indiceSize, MPI_UNSIGNED_LONG,587 ep_lib::MPI_Isend(indices, indiceSize, MPI_UNSIGNED_LONG, 676 588 clientDestRank, MPI_DHT_INDEX, clientIntraComm, &(requestSendIndex.back())); 589 } 590 591 /*! 592 Send message containing index to clients 593 \param [in] clientDestRank rank of destination client 594 \param [in] indices index to send 595 \param [in] indiceSize size of index array to send 596 \param [in] clientIntraComm communication group of client 597 \param [in] requestSendIndex sending request 598 */ 599 template<typename T, typename H> 600 void CClientClientDHTTemplate<T,H>::sendIndexToClients(int clientDestRank, size_t* indices, size_t indiceSize, 601 const ep_lib::MPI_Comm& clientIntraComm, 602 ep_lib::MPI_Request* requestSendIndex) 603 { 604 ep_lib::MPI_Isend(indices, indiceSize, MPI_UNSIGNED_LONG, 605 clientDestRank, MPI_DHT_INDEX, clientIntraComm, requestSendIndex); 677 606 } 678 607 … … 689 618 std::vector<ep_lib::MPI_Request>& requestRecvIndex) 690 619 { 691 printf("should not call this function recvIndexFromClients");692 620 ep_lib::MPI_Request request; 693 621 requestRecvIndex.push_back(request); 694 MPI_Irecv(indices, indiceSize, MPI_UNSIGNED_LONG,622 ep_lib::MPI_Irecv(indices, indiceSize, MPI_UNSIGNED_LONG, 695 623 clientSrcRank, MPI_DHT_INDEX, clientIntraComm, &(requestRecvIndex.back())); 624 } 625 626 /*! 627 Receive message containing index to clients 628 \param [in] clientDestRank rank of destination client 629 \param [in] indices index to send 630 \param [in] clientIntraComm communication group of client 631 \param [in] requestRecvIndex receiving request 632 */ 633 template<typename T, typename H> 634 void CClientClientDHTTemplate<T,H>::recvIndexFromClients(int clientSrcRank, size_t* indices, size_t indiceSize, 635 const ep_lib::MPI_Comm& clientIntraComm, 636 ep_lib::MPI_Request *requestRecvIndex) 637 { 638 ep_lib::MPI_Irecv(indices, indiceSize, MPI_UNSIGNED_LONG, 639 clientSrcRank, MPI_DHT_INDEX, clientIntraComm, requestRecvIndex); 696 640 } 697 641 … … 709 653 std::vector<ep_lib::MPI_Request>& requestSendInfo) 710 654 { 711 printf("should not call this function sendInfoToClients");712 655 ep_lib::MPI_Request request; 713 656 requestSendInfo.push_back(request); 714 MPI_Isend(info, infoSize, MPI_CHAR, 657 658 ep_lib::MPI_Isend(info, infoSize, MPI_CHAR, 715 659 clientDestRank, MPI_DHT_INFO, clientIntraComm, &(requestSendInfo.back())); 660 } 661 662 /*! 663 Send message containing information to clients 664 \param [in] clientDestRank rank of destination client 665 \param [in] info info array to send 666 \param [in] infoSize info array size to send 667 \param [in] clientIntraComm communication group of client 668 \param [in] requestSendInfo sending request 669 */ 670 template<typename T, typename H> 671 void CClientClientDHTTemplate<T,H>::sendInfoToClients(int clientDestRank, unsigned char* info, int infoSize, 672 const ep_lib::MPI_Comm& clientIntraComm, 673 ep_lib::MPI_Request *requestSendInfo) 674 { 675 ep_lib::MPI_Isend(info, infoSize, MPI_CHAR, 676 clientDestRank, MPI_DHT_INFO, clientIntraComm, requestSendInfo); 716 677 } 717 678 … … 729 690 std::vector<ep_lib::MPI_Request>& requestRecvInfo) 730 691 { 731 printf("should not call this function recvInfoFromClients\n");732 692 ep_lib::MPI_Request request; 733 693 requestRecvInfo.push_back(request); 734 694 735 MPI_Irecv(info, infoSize, MPI_CHAR,695 ep_lib::MPI_Irecv(info, infoSize, MPI_CHAR, 736 696 clientSrcRank, MPI_DHT_INFO, clientIntraComm, &(requestRecvInfo.back())); 697 } 698 699 /*! 700 Receive message containing information from other clients 701 \param [in] clientDestRank rank of destination client 702 \param [in] info info array to receive 703 \param [in] infoSize info array size to receive 704 \param [in] clientIntraComm communication group of client 705 \param [in] requestRecvInfo list of receiving request 706 */ 707 template<typename T, typename H> 708 void CClientClientDHTTemplate<T,H>::recvInfoFromClients(int clientSrcRank, unsigned char* info, int infoSize, 709 const ep_lib::MPI_Comm& clientIntraComm, 710 ep_lib::MPI_Request* requestRecvInfo) 711 { 712 ep_lib::MPI_Irecv(info, infoSize, MPI_CHAR, 713 clientSrcRank, MPI_DHT_INFO, clientIntraComm, requestRecvInfo); 737 714 } 738 715 … … 807 784 808 785 int nRequest = 0; 809 786 for (int idx = 0; idx < recvNbRank.size(); ++idx) 787 { 788 ep_lib::MPI_Irecv(&recvNbElements[0]+idx, 1, MPI_INT, 789 recvNbRank[idx], MPI_DHT_INDEX_1, this->internalComm_, &request[nRequest]); 790 ++nRequest; 791 } 810 792 811 793 for (int idx = 0; idx < sendNbRank.size(); ++idx) 812 794 { 813 MPI_Isend(&sendNbElements[0]+idx, 1, MPI_INT,795 ep_lib::MPI_Isend(&sendNbElements[0]+idx, 1, MPI_INT, 814 796 sendNbRank[idx], MPI_DHT_INDEX_1, this->internalComm_, &request[nRequest]); 815 797 ++nRequest; 816 798 } 817 818 for (int idx = 0; idx < recvNbRank.size(); ++idx) 819 { 820 MPI_Irecv(&recvNbElements[0]+idx, 1, MPI_INT, 821 recvNbRank[idx], MPI_DHT_INDEX_1, this->internalComm_, &request[nRequest]); 822 ++nRequest; 823 } 824 825 MPI_Waitall(sendNbRank.size()+recvNbRank.size(), &request[0], &requestStatus[0]); 799 800 ep_lib::MPI_Waitall(sendNbRank.size()+recvNbRank.size(), &request[0], &requestStatus[0]); 826 801 } 827 802 … … 852 827 std::vector<ep_lib::MPI_Request> request(sendBuffSize+recvBuffSize); 853 828 std::vector<ep_lib::MPI_Status> requestStatus(sendBuffSize+recvBuffSize); 854 //ep_lib::MPI_Request request[sendBuffSize+recvBuffSize]; 855 //ep_lib::MPI_Status requestStatus[sendBuffSize+recvBuffSize]; 856 857 int my_rank; 858 MPI_Comm_rank(this->internalComm_, &my_rank); 859 829 860 830 int nRequest = 0; 861 831 for (int idx = 0; idx < recvBuffSize; ++idx) 862 832 { 863 MPI_Irecv(&recvBuff[2*idx], 2, MPI_INT,833 ep_lib::MPI_Irecv(&recvBuff[0]+2*idx, 2, MPI_INT, 864 834 recvRank[idx], MPI_DHT_INDEX_0, this->internalComm_, &request[nRequest]); 865 835 ++nRequest; 866 836 } 867 868 837 869 838 for (int idx = 0; idx < sendBuffSize; ++idx) … … 873 842 sendBuff[idx*2+1] = sendNbElements[offSet]; 874 843 } 875 876 877 844 878 845 for (int idx = 0; idx < sendBuffSize; ++idx) 879 846 { 880 MPI_Isend(&sendBuff[idx*2], 2, MPI_INT,847 ep_lib::MPI_Isend(&sendBuff[idx*2], 2, MPI_INT, 881 848 sendRank[idx], MPI_DHT_INDEX_0, this->internalComm_, &request[nRequest]); 882 849 ++nRequest; 883 850 } 884 885 886 887 //MPI_Barrier(this->internalComm_); 888 889 MPI_Waitall(sendBuffSize+recvBuffSize, &request[0], &requestStatus[0]); 890 //MPI_Waitall(sendBuffSize+recvBuffSize, request, requestStatus); 891 892 851 852 ep_lib::MPI_Waitall(sendBuffSize+recvBuffSize, &request[0], &requestStatus[0]); 893 853 int nbRecvRank = 0, nbRecvElements = 0; 894 854 recvNbRank.clear(); … … 902 862 } 903 863 } 904 905 906 907 908 } 909 910 } 911 864 } 865 866 } -
XIOS/dev/branch_openmp/src/client_server_mapping.cpp
r1287 r1328 8 8 */ 9 9 #include "client_server_mapping.hpp" 10 using namespace ep_lib; 10 11 11 12 namespace xios { … … 64 65 MPI_Allgather(&nbConnectedServer,1,MPI_INT,recvCount,1,MPI_INT,clientIntraComm) ; 65 66 66 67 // for(int i=0; i<nbClient; i++)68 // printf("MPI_Allgather : recvCount[%d] = %d\n", i, recvCount[i]);69 70 67 displ[0]=0 ; 71 68 for(int n=1;n<nbClient;n++) displ[n]=displ[n-1]+recvCount[n-1] ; … … 75 72 76 73 MPI_Allgatherv(sendBuff,nbConnectedServer,MPI_INT,recvBuff,recvCount,displ,MPI_INT,clientIntraComm) ; 77 78 // for(int i=0; i<recvSize; i++)79 // printf("MPI_Allgatherv : recvBuff[%d] = %d\n", i, recvBuff[i]);80 81 82 74 for(int n=0;n<recvSize;n++) clientRes[recvBuff[n]]++ ; 83 75 -
XIOS/dev/branch_openmp/src/client_server_mapping.hpp
r1134 r1328 14 14 #include "mpi.hpp" 15 15 #include <boost/unordered_map.hpp> 16 #ifdef _usingEP17 #include "ep_declaration.hpp"18 #endif19 20 16 21 17 namespace xios { -
XIOS/dev/branch_openmp/src/client_server_mapping_distributed.cpp
r907 r1328 15 15 #include "context.hpp" 16 16 #include "context_client.hpp" 17 using namespace ep_lib; 17 18 18 19 namespace xios -
XIOS/dev/branch_openmp/src/context_client.cpp
r1205 r1328 11 11 #include "timer.hpp" 12 12 #include "cxios.hpp" 13 using namespace ep_lib; 13 14 14 15 namespace xios … … 20 21 \cxtSer [in] cxtSer Pointer to context of server side. (It is only used on case of attached mode) 21 22 */ 22 CContextClient::CContextClient(CContext* parent, ep_lib::MPI_Comm intraComm_, ep_lib::MPI_Comm interComm_, CContext* cxtSer)23 CContextClient::CContextClient(CContext* parent, MPI_Comm intraComm_, MPI_Comm interComm_, CContext* cxtSer) 23 24 : mapBufferSize_(), parentServer(cxtSer), maxBufferedEvents(4) 24 25 { … … 293 294 if (ratio < minBufferSizeEventSizeRatio) minBufferSizeEventSizeRatio = ratio; 294 295 } 295 #ifdef _usingMPI 296 MPI_Allreduce(MPI_IN_PLACE, &minBufferSizeEventSizeRatio, 1, MPI_DOUBLE, MPI_MIN, intraComm); 297 #elif _usingEP 296 //MPI_Allreduce(MPI_IN_PLACE, &minBufferSizeEventSizeRatio, 1, MPI_DOUBLE, MPI_MIN, intraComm); 298 297 MPI_Allreduce(&minBufferSizeEventSizeRatio, &minBufferSizeEventSizeRatio, 1, MPI_DOUBLE, MPI_MIN, intraComm); 299 #endif 300 298 301 299 if (minBufferSizeEventSizeRatio < 1.0) 302 300 { … … 402 400 for (itMap = itbMap; itMap != iteMap; ++itMap) 403 401 { 404 //report(10) << " Memory report : Context <" << context->getId() << "> : client side : memory used for buffer of each connection to server" << endl405 //<< " +) To server with rank " << itMap->first << " : " << itMap->second << " bytes " << endl;402 report(10) << " Memory report : Context <" << context->getId() << "> : client side : memory used for buffer of each connection to server" << endl 403 << " +) To server with rank " << itMap->first << " : " << itMap->second << " bytes " << endl; 406 404 totalBuf += itMap->second; 407 405 } 408 //report(0) << " Memory report : Context <" << context->getId() << "> : client side : total memory used for buffer " << totalBuf << " bytes" << endl;406 report(0) << " Memory report : Context <" << context->getId() << "> : client side : total memory used for buffer " << totalBuf << " bytes" << endl; 409 407 410 408 releaseBuffers(); -
XIOS/dev/branch_openmp/src/context_server.cpp
r1179 r1328 10 10 #include "file.hpp" 11 11 #include "grid.hpp" 12 #include "mpi _std.hpp"12 #include "mpi.hpp" 13 13 #include "tracer.hpp" 14 14 #include "timer.hpp" … … 18 18 #include <boost/functional/hash.hpp> 19 19 20 20 using namespace ep_lib; 21 21 22 22 namespace xios 23 23 { 24 24 25 CContextServer::CContextServer(CContext* parent, ep_lib::MPI_Comm intraComm_, ep_lib::MPI_Comm interComm_)25 CContextServer::CContextServer(CContext* parent,MPI_Comm intraComm_,MPI_Comm interComm_) 26 26 { 27 27 context=parent; … … 72 72 int count; 73 73 char * addr; 74 ep_lib::MPI_Status status;74 MPI_Status status; 75 75 map<int,CServerBuffer*>::iterator it; 76 77 for(rank=0;rank<commSize;rank++) 78 { 76 bool okLoop; 77 78 traceOff(); 79 MPI_Iprobe(-2, 20,interComm,&flag,&status); 80 traceOn(); 81 82 if (flag==true) 83 { 84 #ifdef _usingMPI 85 rank=status.MPI_SOURCE ; 86 #elif _usingEP 87 rank=status.ep_src ; 88 #endif 89 okLoop = true; 79 90 if (pendingRequest.find(rank)==pendingRequest.end()) 80 { 81 traceOff(); 82 ep_lib::MPI_Iprobe(rank,20,interComm,&flag,&status); 83 traceOn(); 84 if (flag) 91 okLoop = !listenPendingRequest(status) ; 92 if (okLoop) 93 { 94 for(rank=0;rank<commSize;rank++) 85 95 { 86 it=buffers.find(rank); 87 if (it==buffers.end()) // Receive the buffer size and allocate the buffer 96 if (pendingRequest.find(rank)==pendingRequest.end()) 88 97 { 89 StdSize buffSize = 0; 90 ep_lib::MPI_Recv(&buffSize, 1, MPI_LONG, rank, 20, interComm, &status); 91 mapBufferSize_.insert(std::make_pair(rank, buffSize)); 92 it=(buffers.insert(pair<int,CServerBuffer*>(rank,new CServerBuffer(buffSize)))).first; 93 } 94 else 95 { 96 97 ep_lib::MPI_Get_count(&status,MPI_CHAR,&count); 98 if (it->second->isBufferFree(count)) 99 { 100 addr=(char*)it->second->getBuffer(count); 101 ep_lib::MPI_Irecv(addr,count,MPI_CHAR,rank,20,interComm,&pendingRequest[rank]); 102 bufferRequest[rank]=addr; 103 } 98 99 traceOff(); 100 MPI_Iprobe(rank, 20,interComm,&flag,&status); 101 traceOn(); 102 if (flag==true) listenPendingRequest(status) ; 104 103 } 105 104 } … … 108 107 } 109 108 109 bool CContextServer::listenPendingRequest(MPI_Status& status) 110 { 111 int count; 112 char * addr; 113 map<int,CServerBuffer*>::iterator it; 114 #ifdef _usingMPI 115 int rank=status.MPI_SOURCE ; 116 #elif _usingEP 117 int rank=status.ep_src; 118 #endif 119 120 it=buffers.find(rank); 121 if (it==buffers.end()) // Receive the buffer size and allocate the buffer 122 { 123 StdSize buffSize = 0; 124 MPI_Recv(&buffSize, 1, MPI_LONG, rank, 20, interComm, &status); 125 mapBufferSize_.insert(std::make_pair(rank, buffSize)); 126 it=(buffers.insert(pair<int,CServerBuffer*>(rank,new CServerBuffer(buffSize)))).first; 127 return true; 128 } 129 else 130 { 131 MPI_Get_count(&status,MPI_CHAR,&count); 132 if (it->second->isBufferFree(count)) 133 { 134 addr=(char*)it->second->getBuffer(count); 135 MPI_Irecv(addr,count,MPI_CHAR,rank,20,interComm,&pendingRequest[rank]); 136 bufferRequest[rank]=addr; 137 return true; 138 } 139 else 140 return false; 141 } 142 } 143 110 144 void CContextServer::checkPendingRequest(void) 111 145 { 112 map<int, ep_lib::MPI_Request>::iterator it;146 map<int,MPI_Request>::iterator it; 113 147 list<int> recvRequest; 114 148 list<int>::iterator itRecv; … … 116 150 int flag; 117 151 int count; 118 ep_lib::MPI_Status status;119 120 for(it=pendingRequest.begin();it!=pendingRequest.end(); ++it)152 MPI_Status status; 153 154 for(it=pendingRequest.begin();it!=pendingRequest.end();it++) 121 155 { 122 156 rank=it->first; 123 157 traceOff(); 124 ep_lib::MPI_Test(& it->second, &flag, &status);158 MPI_Test(& it->second, &flag, &status); 125 159 traceOn(); 126 160 if (flag==true) 127 161 { 128 162 recvRequest.push_back(rank); 129 ep_lib::MPI_Get_count(&status,MPI_CHAR,&count);163 MPI_Get_count(&status,MPI_CHAR,&count); 130 164 processRequest(rank,bufferRequest[rank],count); 131 165 } … … 220 254 { 221 255 finished=true; 222 #pragma omp critical (_output)223 256 info(20)<<"Server Side context <"<<context->getId()<<"> finalized"<<endl; 224 257 std::map<int, StdSize>::const_iterator itbMap = mapBufferSize_.begin(), … … 227 260 for (itMap = itbMap; itMap != iteMap; ++itMap) 228 261 { 229 //report(10)<< " Memory report : Context <"<<context->getId()<<"> : server side : memory used for buffer of each connection to client" << endl230 //<< " +) With client of rank " << itMap->first << " : " << itMap->second << " bytes " << endl;262 report(10)<< " Memory report : Context <"<<context->getId()<<"> : server side : memory used for buffer of each connection to client" << endl 263 << " +) With client of rank " << itMap->first << " : " << itMap->second << " bytes " << endl; 231 264 totalBuf += itMap->second; 232 265 } 233 266 context->finalize(); 234 //report(0)<< " Memory report : Context <"<<context->getId()<<"> : server side : total memory used for buffer "<<totalBuf<<" bytes"<<endl;267 report(0)<< " Memory report : Context <"<<context->getId()<<"> : server side : total memory used for buffer "<<totalBuf<<" bytes"<<endl; 235 268 } 236 269 else if (event.classId==CContext::GetType()) CContext::dispatchEvent(event); -
XIOS/dev/branch_openmp/src/context_server.hpp
r1134 r1328 17 17 bool eventLoop(bool enableEventsProcessing = true); 18 18 void listen(void) ; 19 bool listenPendingRequest(ep_lib::MPI_Status& status); 19 20 void checkPendingRequest(void) ; 20 21 void processRequest(int rank, char* buff,int count) ; -
XIOS/dev/branch_openmp/src/cxios.cpp
r1205 r1328 11 11 #include "memtrack.hpp" 12 12 #include "registry.hpp" 13 using namespace ep_lib; 13 14 14 15 namespace xios 15 16 { 16 17 extern int test_omp_rank; 18 #pragma omp threadprivate(test_omp_rank) 19 20 const string CXios::rootFile="./iodef.xml" ; 21 const string CXios::xiosCodeId="xios.x" ; 22 const string CXios::clientFile="./xios_client"; 23 const string CXios::serverFile="./xios_server"; 24 17 string CXios::rootFile="./iodef.xml" ; 18 string CXios::xiosCodeId="xios.x" ; 19 string CXios::clientFile="./xios_client"; 20 string CXios::serverFile="./xios_server"; 25 21 26 22 bool CXios::isClient ; 27 23 bool CXios::isServer ; 28 29 30 24 MPI_Comm CXios::globalComm ; 31 32 33 25 bool CXios::usingOasis ; 34 26 bool CXios::usingServer = false; 35 36 37 27 double CXios::bufferSizeFactor = 1.0; 38 28 const double CXios::defaultBufferSizeFactor = 1.0; 39 29 StdSize CXios::minBufferSize = 1024 * sizeof(double); 40 41 42 30 bool CXios::printLogs2Files; 43 31 bool CXios::isOptPerformance = true; … … 49 37 { 50 38 set_new_handler(noMemory); 51 52 53 #pragma omp critical 54 { 55 parseFile(rootFile); 56 } 57 #pragma omp barrier 39 parseFile(rootFile); 58 40 parseXiosConfig(); 59 41 } … … 87 69 ERROR("CXios::parseXiosConfig()", "recv_field_timeout cannot be negative."); 88 70 89 71 //globalComm=MPI_COMM_WORLD ; 90 72 int num_ep; 91 73 if(isClient) … … 93 75 num_ep = omp_get_num_threads(); 94 76 } 95 77 96 78 if(isServer) 97 79 { 98 80 num_ep = omp_get_num_threads(); 99 81 } 100 82 101 83 MPI_Info info; 102 84 #pragma omp master … … 106 88 passage = ep_comm; 107 89 } 108 90 109 91 #pragma omp barrier 110 111 92 93 112 94 CXios::globalComm = passage[omp_get_thread_num()]; 113 114 int tmp_rank;115 MPI_Comm_rank(CXios::globalComm, &tmp_rank);116 117 118 test_omp_rank = tmp_rank;119 120 95 } 121 96 … … 133 108 134 109 CClient::initialize(codeId,localComm,returnComm) ; 135 136 110 if (CClient::getRank()==0) globalRegistry = new CRegistry(returnComm) ; 137 111 … … 142 116 if (printLogs2Files) 143 117 { 144 #pragma omp critical145 118 CClient::openInfoStream(clientFile); 146 119 CClient::openErrorStream(clientFile); … … 158 131 if (CClient::getRank()==0) 159 132 { 160 #pragma omp critical (_output)161 133 info(80)<<"Write data base Registry"<<endl<<globalRegistry->toString()<<endl ; 162 134 globalRegistry->toFile("xios_registry.bin") ; … … 184 156 void CXios::initServer() 185 157 { 186 int initialized;187 MPI_Initialized(&initialized);188 if (initialized) CServer::is_MPI_Initialized=true ;189 else CServer::is_MPI_Initialized=false ;190 191 192 if(!CServer::is_MPI_Initialized)193 {194 MPI_Init(NULL, NULL);195 }196 197 158 set_new_handler(noMemory); 198 159 std::set<StdString> parseList; … … 210 171 211 172 initServer(); 212 213 173 214 174 // Initialize all aspects MPI 215 175 CServer::initialize(); -
XIOS/dev/branch_openmp/src/cxios.hpp
r1134 r1328 5 5 #include "mpi.hpp" 6 6 #include "registry.hpp" 7 #include "log.hpp"8 7 9 8 namespace xios … … 15 14 { 16 15 public: 17 18 19 20 21 16 static void initialize(void) ; 17 static void initClientSide(const string & codeId, ep_lib::MPI_Comm& localComm, ep_lib::MPI_Comm& returnComm) ; 18 static void initServerSide(void) ; 19 static void clientFinalize(void) ; 20 static void parseFile(const string& filename) ; 22 21 23 24 22 template <typename T> 23 static T getin(const string& id,const T& defaultValue) ; 25 24 26 27 25 template <typename T> 26 static T getin(const string& id) ; 28 27 29 28 public: 30 static const string rootFile; //!< Configuration filename31 static conststring xiosCodeId ; //!< Identity for XIOS32 static conststring clientFile; //!< Filename template for client33 static conststring serverFile; //!< Filename template for server29 static string rootFile ; //!< Configuration filename 30 static string xiosCodeId ; //!< Identity for XIOS 31 static string clientFile; //!< Filename template for client 32 static string serverFile; //!< Filename template for server 34 33 35 static bool isClient ; //!< Check if xios is client 36 static bool isServer ; //!< Check if xios is server 37 #pragma omp threadprivate(isClient, isServer) 34 static bool isClient ; //!< Check if xios is client 35 static bool isServer ; //!< Check if xios is server 38 36 39 static MPI_Comm globalComm ; //!< Global communicator 40 #pragma omp threadprivate(globalComm) 37 static ep_lib::MPI_Comm globalComm ; //!< Global communicator 41 38 42 static bool printLogs2Files; //!< Printing out logs into files 43 static bool usingOasis ; //!< Using Oasis 44 static bool usingServer ; //!< Using server (server mode) 45 static double bufferSizeFactor; //!< Factor used to tune the buffer size 46 static const double defaultBufferSizeFactor; //!< Default factor value 47 static StdSize minBufferSize; //!< Minimum buffer size 48 static bool isOptPerformance; //!< Check if buffer size is for performance (as large as possible) 49 #pragma omp threadprivate(printLogs2Files, usingOasis, usingServer, bufferSizeFactor, minBufferSize, isOptPerformance) 50 51 static CRegistry* globalRegistry ; //!< global registry which is wrote by the root process of the servers 52 static double recvFieldTimeout; 53 #pragma omp threadprivate(recvFieldTimeout) 54 55 39 static bool printLogs2Files; //!< Printing out logs into files 40 static bool usingOasis ; //!< Using Oasis 41 static bool usingServer ; //!< Using server (server mode) 42 static double bufferSizeFactor; //!< Factor used to tune the buffer size 43 static const double defaultBufferSizeFactor; //!< Default factor value 44 static StdSize minBufferSize; //!< Minimum buffer size 45 static bool isOptPerformance; //!< Check if buffer size is for performance (as large as possible) 46 static CRegistry* globalRegistry ; //!< global registry which is wrote by the root process of the servers 47 static double recvFieldTimeout; //!< Time to wait for data before issuing an error when receiving a field 48 56 49 public: 57 50 //! Setting xios to use server mode -
XIOS/dev/branch_openmp/src/data_output.cpp
r1205 r1328 4 4 #include "group_template.hpp" 5 5 #include "context.hpp" 6 //mpi.hpp 6 7 7 namespace xios 8 8 { -
XIOS/dev/branch_openmp/src/data_output.hpp
r1287 r1328