Changeset 1642
- Timestamp:
- 01/23/19 10:31:44 (6 years ago)
- Location:
- XIOS/dev/branch_openmp
- Files:
-
- 257 edited
Legend:
- Unmodified
- Added
- Removed
-
XIOS/dev/branch_openmp/extern/remap/src/cputime.cpp
r1328 r1642 1 1 #include "mpi.hpp" 2 using namespace ep_lib;3 2 4 3 namespace sphereRemap { -
XIOS/dev/branch_openmp/extern/remap/src/elt.hpp
r1328 r1642 48 48 int n; /* number of vertices */ 49 49 double area; 50 double given_area ; 50 51 Coord x; /* barycentre */ 51 52 }; … … 80 81 n = rhs.n; 81 82 area = rhs.area; 83 given_area = rhs.given_area; 82 84 x = rhs.x; 83 85 val = rhs.val; -
XIOS/dev/branch_openmp/extern/remap/src/gridRemap.cpp
r1335 r1642 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
r1460 r1642 14 14 Coord readPole(std::istream&); 15 15 16 //extern CRemapGrid srcGrid;17 //extern CRemapGrid tgtGrid;16 extern CRemapGrid srcGrid; 17 extern CRemapGrid tgtGrid; 18 18 19 19 } -
XIOS/dev/branch_openmp/extern/remap/src/intersect.cpp
r1335 r1642 14 14 15 15 namespace sphereRemap { 16 17 extern CRemapGrid srcGrid;18 #pragma omp threadprivate(srcGrid)19 20 extern CRemapGrid tgtGrid;21 #pragma omp threadprivate(tgtGrid)22 16 23 17 using namespace std; -
XIOS/dev/branch_openmp/extern/remap/src/intersection_ym.cpp
r1335 r1642 8 8 #include <stdlib.h> 9 9 #include <limits> 10 #include <array> 11 #include <cstdint> 12 #include "earcut.hpp" 13 #include <fstream> 14 10 15 11 16 #define epsilon 1e-3 // epsilon distance ratio over side lenght for approximate small circle by great circle … … 14 19 namespace sphereRemap { 15 20 16 extern CRemapGrid srcGrid;17 #pragma omp threadprivate(srcGrid)18 19 extern CRemapGrid tgtGrid;20 #pragma omp threadprivate(tgtGrid)21 22 23 21 using namespace std; 24 22 using namespace ClipperLib ; 23 25 24 26 25 double intersect_ym(Elt *a, Elt *b) 27 26 { 27 28 using N = uint32_t; 29 using Point = array<double, 2>; 30 vector<Point> vect_points; 31 vector< vector<Point> > polyline; 28 32 29 33 // transform small circle into piece of great circle if necessary … … 52 56 double cos_alpha; 53 57 58 /// vector<p2t::Point*> polyline; 54 59 for(int n=0; n<na;n++) 55 60 { … … 58 63 a_gno[n].y=scalarprod(dstPolygon[n],Oy)/cos_alpha ; 59 64 a_gno[n].z=scalarprod(dstPolygon[n],Oz)/cos_alpha ; // must be equal to 1 60 } 65 66 vect_points.push_back( array<double, 2>() ); 67 vect_points[n][0] = a_gno[n].x; 68 vect_points[n][1] = a_gno[n].y; 69 70 } 71 72 polyline.push_back(vect_points); 73 vector<N> indices_a_gno = mapbox::earcut<N>(polyline); 74 75 double area_a_gno=0 ; 76 for(int i=0;i<indices_a_gno.size()/3;++i) 77 { 78 Coord x0 = Ox * polyline[0][indices_a_gno[3*i]][0] + Oy* polyline[0][indices_a_gno[3*i]][1] + Oz ; 79 Coord x1 = Ox * polyline[0][indices_a_gno[3*i+1]][0] + Oy* polyline[0][indices_a_gno[3*i+1]][1] + Oz ; 80 Coord x2 = Ox * polyline[0][indices_a_gno[3*i+2]][0] + Oy* polyline[0][indices_a_gno[3*i+2]][1] + Oz ; 81 area_a_gno+=triarea(x0 * (1./norm(x0)),x1* (1./norm(x1)), x2* (1./norm(x2))) ; 82 } 83 84 vect_points.clear(); 85 polyline.clear(); 86 indices_a_gno.clear(); 87 88 61 89 62 90 for(int n=0; n<nb;n++) … … 66 94 b_gno[n].y=scalarprod(srcPolygon[n],Oy)/cos_alpha ; 67 95 b_gno[n].z=scalarprod(srcPolygon[n],Oz)/cos_alpha ; // must be equal to 1 68 } 69 96 97 vect_points.push_back( array<double, 2>() ); 98 vect_points[n][0] = b_gno[n].x; 99 vect_points[n][1] = b_gno[n].y; 100 } 101 102 103 polyline.push_back(vect_points); 104 vector<N> indices_b_gno = mapbox::earcut<N>(polyline); 105 106 double area_b_gno=0 ; 107 for(int i=0;i<indices_b_gno.size()/3;++i) 108 { 109 Coord x0 = Ox * polyline[0][indices_b_gno[3*i]][0] + Oy* polyline[0][indices_b_gno[3*i]][1] + Oz ; 110 Coord x1 = Ox * polyline[0][indices_b_gno[3*i+1]][0] + Oy* polyline[0][indices_b_gno[3*i+1]][1] + Oz ; 111 Coord x2 = Ox * polyline[0][indices_b_gno[3*i+2]][0] + Oy* polyline[0][indices_b_gno[3*i+2]][1] + Oz ; 112 area_b_gno+=triarea(x0 * (1./norm(x0)),x1* (1./norm(x1)), x2* (1./norm(x2))) ; 113 } 114 115 vect_points.clear(); 116 polyline.clear(); 117 indices_b_gno.clear(); 70 118 71 119 … … 116 164 clip.AddPaths(dst, ptClip, true); 117 165 clip.Execute(ctIntersection, intersection); 118 166 119 167 double area=0 ; 120 168 121 169 for(int ni=0;ni<intersection.size(); ni++) 122 170 { 123 // go back into real coordinate on the sphere124 171 Coord* intersectPolygon=new Coord[intersection[ni].size()] ; 125 172 for(int n=0; n < intersection[ni].size(); n++) 126 173 { 127 double x=intersection[ni][n].X/xscale+xoffset ; 128 double y=intersection[ni][n].Y/yscale+yoffset ; 129 130 intersectPolygon[n]=Ox*x+Oy*y+Oz ; 131 intersectPolygon[n]=intersectPolygon[n]*(1./norm(intersectPolygon[n])) ; 132 } 133 134 // remove redondants vertex 135 int nv=0 ; 174 intersectPolygon[n].x=intersection[ni][n].X/xscale+xoffset ; 175 intersectPolygon[n].y=intersection[ni][n].Y/yscale+yoffset ; 176 } 177 178 179 int nv=0; 180 136 181 for(int n=0; n < intersection[ni].size(); n++) 137 182 { 138 if (norm(intersectPolygon[n]-intersectPolygon[(n+1)%intersection[ni].size()])>fusion_vertex) 183 double dx=intersectPolygon[n].x-intersectPolygon[(n+1)%intersection[ni].size()].x ; 184 double dy=intersectPolygon[n].y-intersectPolygon[(n+1)%intersection[ni].size()].y ; 185 186 if (dx*dx+dy*dy>fusion_vertex*fusion_vertex) 187 { 188 intersectPolygon[nv]=intersectPolygon[n] ; 189 190 vect_points.push_back( array<double, 2>() ); 191 vect_points[nv][0] = intersectPolygon[n].x; 192 vect_points[nv][1] = intersectPolygon[n].y; 193 194 nv++ ; 195 } 196 197 198 } 199 200 polyline.push_back(vect_points); 201 vect_points.clear(); 202 203 if (nv>2) 204 { 205 206 vector<N> indices = mapbox::earcut<N>(polyline); 207 208 double area2=0 ; 209 for(int i=0;i<indices.size()/3;++i) 139 210 { 140 intersectPolygon[nv]=intersectPolygon[n] ; 141 nv++ ; 211 Coord x0 = Ox * polyline[0][indices[3*i]][0] + Oy* polyline[0][indices[3*i]][1] + Oz ; 212 Coord x1 = Ox * polyline[0][indices[3*i+1]][0] + Oy* polyline[0][indices[3*i+1]][1] + Oz ; 213 Coord x2 = Ox * polyline[0][indices[3*i+2]][0] + Oy* polyline[0][indices[3*i+2]][1] + Oz ; 214 area2+=triarea(x0 * (1./norm(x0)),x1* (1./norm(x1)), x2* (1./norm(x2))) ; 142 215 } 143 } 144 145 146 if (nv>2) 147 { 216 217 polyline.clear(); 218 219 for(int n=0; n < nv; n++) 220 { 221 intersectPolygon[n] = Ox*intersectPolygon[n].x+Oy*intersectPolygon[n].y+Oz; 222 intersectPolygon[n] = intersectPolygon[n]*(1./norm(intersectPolygon[n])) ; 223 } 224 225 148 226 // assign intersection to source and destination polygons 149 227 Polyg *is = new Polyg; 150 228 is->x = exact_barycentre(intersectPolygon,nv); 151 is->area = polygonarea(intersectPolygon,nv) ; 229 // is->area = polygonarea(intersectPolygon,nv) ; 230 is->area = area2 ; 231 152 232 // if (is->area < 1e-12) cout<<"Small intersection : "<<is->area<<endl ; 153 233 if (is->area==0.) delete is ; … … 168 248 delete[] b_gno ; 169 249 return area ; 250 170 251 } 252 253 171 254 172 255 void createGreatCirclePolygon(const Elt& element, const Coord& pole, vector<Coord>& coordinates) -
XIOS/dev/branch_openmp/extern/remap/src/libmapper.cpp
r1335 r1642 15 15 #include "cputime.hpp" // cputime 16 16 17 using namespace ep_lib;18 19 17 using namespace sphereRemap ; 20 18 … … 22 20 and deallocated during the second step (computing the weights) */ 23 21 Mapper *mapper; 24 #pragma omp threadprivate(mapper) 22 25 23 26 24 /** xxx_bounds_yyy is of length n_vert_per_cell_xxx*n_cell_xxx … … 43 41 assert(n_cell_dst >= 4); 44 42 assert(1 <= order && order <= 2); 45 46 mapper = new Mapper(MPI_COMM_WORLD); 43 double* src_area=NULL ; 44 double* dst_area=NULL ; 45 mapper = new Mapper(EP_COMM_WORLD); 47 46 mapper->setVerbosity(PROGRESS) ; 48 mapper->setSourceMesh(src_bounds_lon, src_bounds_lat, n_vert_per_cell_src, n_cell_src, src_pole ) ;49 mapper->setTargetMesh(dst_bounds_lon, dst_bounds_lat, n_vert_per_cell_dst, n_cell_dst, dst_pole ) ;47 mapper->setSourceMesh(src_bounds_lon, src_bounds_lat, src_area, n_vert_per_cell_src, n_cell_src, src_pole ) ; 48 mapper->setTargetMesh(dst_bounds_lon, dst_bounds_lat, dst_area, n_vert_per_cell_dst, n_cell_dst, dst_pole ) ; 50 49 51 50 /* -
XIOS/dev/branch_openmp/extern/remap/src/mapper.cpp
r1545 r1642 12 12 13 13 #include "mapper.hpp" 14 using namespace ep_lib;15 14 16 15 namespace sphereRemap { 17 18 extern CRemapGrid srcGrid;19 #pragma omp threadprivate(srcGrid)20 21 extern CRemapGrid tgtGrid;22 #pragma omp threadprivate(tgtGrid)23 24 16 25 17 /* A subdivition of an array into N sub-arays … … 35 27 36 28 37 void Mapper::setSourceMesh(const double* boundsLon, const double* boundsLat, int nVertex, int nbCells, const double* pole, const long int* globalId)29 void Mapper::setSourceMesh(const double* boundsLon, const double* boundsLat, const double* area, int nVertex, int nbCells, const double* pole, const long int* globalId) 38 30 { 39 31 srcGrid.pole = Coord(pole[0], pole[1], pole[2]); 40 32 41 33 int mpiRank, mpiSize; 42 MPI_Comm_rank(communicator, &mpiRank);43 MPI_Comm_size(communicator, &mpiSize);34 ep_lib::MPI_Comm_rank(communicator, &mpiRank); 35 ep_lib::MPI_Comm_size(communicator, &mpiSize); 44 36 45 37 sourceElements.reserve(nbCells); … … 51 43 long int offset ; 52 44 long int nb=nbCells ; 53 MPI_Scan(&nb,&offset,1,MPI_LONG,MPI_SUM,communicator) ;45 ep_lib::MPI_Scan(&nb,&offset,1,EP_LONG,EP_SUM,communicator) ; 54 46 offset=offset-nb ; 55 47 for(int i=0;i<nbCells;i++) sourceGlobalId[i]=offset+i ; … … 67 59 sourceMesh.push_back(Node(elt.x, cptRadius(elt), &sourceElements.back())); 68 60 cptEltGeom(sourceElements[i], Coord(pole[0], pole[1], pole[2])); 69 } 70 71 } 72 73 void Mapper::setTargetMesh(const double* boundsLon, const double* boundsLat, int nVertex, int nbCells, const double* pole, const long int* globalId) 61 if (area!=NULL) sourceElements[i].given_area=area[i] ; 62 else sourceElements[i].given_area=sourceElements[i].area ; 63 } 64 65 } 66 67 void Mapper::setTargetMesh(const double* boundsLon, const double* boundsLat, const double* area, int nVertex, int nbCells, const double* pole, const long int* globalId) 74 68 { 75 69 tgtGrid.pole = Coord(pole[0], pole[1], pole[2]); 76 70 77 71 int mpiRank, mpiSize; 78 MPI_Comm_rank(communicator, &mpiRank);79 MPI_Comm_size(communicator, &mpiSize);72 ep_lib::MPI_Comm_rank(communicator, &mpiRank); 73 ep_lib::MPI_Comm_size(communicator, &mpiSize); 80 74 81 75 targetElements.reserve(nbCells); … … 87 81 long int offset ; 88 82 long int nb=nbCells ; 89 MPI_Scan(&nb,&offset,1,MPI_LONG,MPI_SUM,communicator) ;83 ep_lib::MPI_Scan(&nb,&offset,1,EP_LONG,EP_SUM,communicator) ; 90 84 offset=offset-nb ; 91 85 for(int i=0;i<nbCells;i++) targetGlobalId[i]=offset+i ; … … 100 94 targetMesh.push_back(Node(elt.x, cptRadius(elt), &sourceElements.back())); 101 95 cptEltGeom(targetElements[i], Coord(pole[0], pole[1], pole[2])); 96 if (area!=NULL) targetElements[i].given_area=area[i] ; 97 else targetElements[i].given_area=targetElements[i].area ; 102 98 } 103 99 … … 121 117 vector<double> timings; 122 118 int mpiSize, mpiRank; 123 MPI_Comm_size(communicator, &mpiSize);124 MPI_Comm_rank(communicator, &mpiRank);119 ep_lib::MPI_Comm_size(communicator, &mpiSize); 120 ep_lib::MPI_Comm_rank(communicator, &mpiRank); 125 121 126 122 this->buildSSTree(sourceMesh, targetMesh); … … 132 128 133 129 tic = cputime(); 134 if (interpOrder == 2) 135 { 130 if (interpOrder == 2) { 136 131 if (mpiRank == 0 && verbose) cout << "Computing grads ..." << endl; 137 132 buildMeshTopology(); … … 142 137 /* Prepare computation of weights */ 143 138 /* compute number of intersections which for the first order case 144 corresponds to the number of edges in the remap matrix */139 corresponds to the number of edges in the remap matrix */ 145 140 int nIntersections = 0; 146 141 for (int j = 0; j < targetElements.size(); j++) … … 178 173 { 179 174 int mpiSize, mpiRank; 180 MPI_Comm_size(communicator, &mpiSize);181 MPI_Comm_rank(communicator, &mpiRank);175 ep_lib::MPI_Comm_size(communicator, &mpiSize); 176 ep_lib::MPI_Comm_rank(communicator, &mpiRank); 182 177 183 178 /* create list of intersections (super mesh elements) for each rank */ … … 187 182 Elt& e = elements[j]; 188 183 for (list<Polyg *>::iterator it = e.is.begin(); it != e.is.end(); it++) 189 {190 184 elementList[(*it)->id.rank].insert(pair<int, Polyg *>((*it)->id.ind, *it)); 191 }192 185 } 193 186 … … 196 189 double **recvValue = new double*[mpiSize]; 197 190 double **recvArea = new double*[mpiSize]; 191 double **recvGivenArea = new double*[mpiSize]; 198 192 Coord **recvGrad = new Coord*[mpiSize]; 199 193 GloId **recvNeighIds = new GloId*[mpiSize]; /* ids of the of the source neighbours which also contribute through gradient */ … … 218 212 recvValue[rank] = new double[nbSendElement[rank]]; 219 213 recvArea[rank] = new double[nbSendElement[rank]]; 214 recvGivenArea[rank] = new double[nbSendElement[rank]]; 220 215 if (order == 2) 221 216 { … … 240 235 /* communicate sizes of source elements to be sent (index lists and later values and gradients) */ 241 236 int *nbRecvElement = new int[mpiSize]; 242 MPI_Alltoall(nbSendElement, 1, MPI_INT, nbRecvElement, 1, MPI_INT, communicator); 243 237 ep_lib::MPI_Alltoall(nbSendElement, 1, EP_INT, nbRecvElement, 1, EP_INT, communicator); 244 238 245 239 /* communicate indices of source elements on other ranks whoes value and gradient we need (since intersection) */ … … 249 243 double **sendValue = new double*[mpiSize]; 250 244 double **sendArea = new double*[mpiSize]; 245 double **sendGivenArea = new double*[mpiSize]; 251 246 Coord **sendGrad = new Coord*[mpiSize]; 252 247 GloId **sendNeighIds = new GloId*[mpiSize]; 253 MPI_Request *sendRequest = new MPI_Request[4*mpiSize];254 MPI_Request *recvRequest = new MPI_Request[4*mpiSize];248 ep_lib::MPI_Request *sendRequest = new ep_lib::MPI_Request[5*mpiSize]; 249 ep_lib::MPI_Request *recvRequest = new ep_lib::MPI_Request[5*mpiSize]; 255 250 for (int rank = 0; rank < mpiSize; rank++) 256 251 { 257 252 if (nbSendElement[rank] > 0) 258 253 { 259 MPI_Issend(sendElement[rank], nbSendElement[rank], MPI_INT, rank, 0, communicator, &sendRequest[nbSendRequest]);254 ep_lib::MPI_Issend(sendElement[rank], nbSendElement[rank], EP_INT, rank, 0, communicator, &sendRequest[nbSendRequest]); 260 255 nbSendRequest++; 261 256 } … … 266 261 sendValue[rank] = new double[nbRecvElement[rank]]; 267 262 sendArea[rank] = new double[nbRecvElement[rank]]; 263 sendGivenArea[rank] = new double[nbRecvElement[rank]]; 268 264 if (order == 2) 269 265 { … … 275 271 sendNeighIds[rank] = new GloId[nbRecvElement[rank]]; 276 272 } 277 MPI_Irecv(recvElement[rank], nbRecvElement[rank], MPI_INT, rank, 0, communicator, &recvRequest[nbRecvRequest]);273 ep_lib::MPI_Irecv(recvElement[rank], nbRecvElement[rank], EP_INT, rank, 0, communicator, &recvRequest[nbRecvRequest]); 278 274 nbRecvRequest++; 279 275 } 280 276 } 281 MPI_Status *status = new MPI_Status[4*mpiSize];282 283 MPI_Waitall(nbSendRequest, sendRequest, status);284 MPI_Waitall(nbRecvRequest, recvRequest, status);277 ep_lib::MPI_Status *status = new ep_lib::MPI_Status[5*mpiSize]; 278 279 ep_lib::MPI_Waitall(nbSendRequest, sendRequest, status); 280 ep_lib::MPI_Waitall(nbRecvRequest, recvRequest, status); 285 281 286 282 /* for all indices that have been received from requesting ranks: pack values and gradients, then send */ … … 296 292 sendValue[rank][j] = sstree.localElements[recvElement[rank][j]].val; 297 293 sendArea[rank][j] = sstree.localElements[recvElement[rank][j]].area; 294 sendGivenArea[rank][j] = sstree.localElements[recvElement[rank][j]].given_area; 298 295 if (order == 2) 299 296 { 300 297 sendGrad[rank][jj] = sstree.localElements[recvElement[rank][j]].grad; 298 // cout<<"grad "<<jj<<" "<<recvElement[rank][j]<<" "<<sendGrad[rank][jj]<<" "<<sstree.localElements[recvElement[rank][j]].grad<<endl ; 301 299 sendNeighIds[rank][jj] = sstree.localElements[recvElement[rank][j]].src_id; 302 300 jj++; … … 304 302 { 305 303 sendGrad[rank][jj] = sstree.localElements[recvElement[rank][j]].gradNeigh[i]; 304 // cout<<"grad "<<jj<<" "<<sendGrad[rank][jj]<<" "<<sstree.localElements[recvElement[rank][j]].grad<<endl ; 306 305 sendNeighIds[rank][jj] = sstree.localElements[recvElement[rank][j]].neighId[i]; 307 jj++; 308 } 309 } 310 else 311 sendNeighIds[rank][j] = sstree.localElements[recvElement[rank][j]].src_id; 312 } 313 MPI_Issend(sendValue[rank], nbRecvElement[rank], MPI_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]); 314 nbSendRequest++; 315 MPI_Issend(sendArea[rank], nbRecvElement[rank], MPI_DOUBLE, rank, 1, communicator, &sendRequest[nbSendRequest]); 316 nbSendRequest++; 317 if (order == 2) 318 { 319 MPI_Issend(sendGrad[rank], 3*nbRecvElement[rank]*(NMAX+1), 320 MPI_DOUBLE, rank, 2, communicator, &sendRequest[nbSendRequest]); 321 nbSendRequest++; 322 MPI_Issend(sendNeighIds[rank], 4*nbRecvElement[rank]*(NMAX+1), MPI_INT, rank, 3, communicator, &sendRequest[nbSendRequest]); 306 jj++; 307 } 308 } 309 else 310 sendNeighIds[rank][j] = sstree.localElements[recvElement[rank][j]].src_id; 311 } 312 ep_lib::MPI_Issend(sendValue[rank], nbRecvElement[rank], EP_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]); 313 nbSendRequest++; 314 ep_lib::MPI_Issend(sendArea[rank], nbRecvElement[rank], EP_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]); 315 nbSendRequest++; 316 ep_lib::MPI_Issend(sendGivenArea[rank], nbRecvElement[rank], EP_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]); 317 nbSendRequest++; 318 if (order == 2) 319 { 320 ep_lib::MPI_Issend(sendGrad[rank], 3*nbRecvElement[rank]*(NMAX+1), EP_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]); 321 nbSendRequest++; 322 ep_lib::MPI_Issend(sendNeighIds[rank], 4*nbRecvElement[rank]*(NMAX+1), EP_INT, rank, 0, communicator, &sendRequest[nbSendRequest]); 323 323 //ym --> attention taille GloId 324 325 326 327 328 MPI_Issend(sendNeighIds[rank], 4*nbRecvElement[rank], MPI_INT, rank, 4, communicator, &sendRequest[nbSendRequest]);324 nbSendRequest++; 325 } 326 else 327 { 328 ep_lib::MPI_Issend(sendNeighIds[rank], 4*nbRecvElement[rank], EP_INT, rank, 0, communicator, &sendRequest[nbSendRequest]); 329 329 //ym --> attention taille GloId 330 nbSendRequest++; 331 } 332 } 333 if (nbSendElement[rank] > 0) 334 { 335 MPI_Irecv(recvValue[rank], nbSendElement[rank], MPI_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]); 336 nbRecvRequest++; 337 MPI_Irecv(recvArea[rank], nbSendElement[rank], MPI_DOUBLE, rank, 1, communicator, &recvRequest[nbRecvRequest]); 338 nbRecvRequest++; 339 if (order == 2) 340 { 341 MPI_Irecv(recvGrad[rank], 3*nbSendElement[rank]*(NMAX+1), 342 MPI_DOUBLE, rank, 2, communicator, &recvRequest[nbRecvRequest]); 343 nbRecvRequest++; 344 MPI_Irecv(recvNeighIds[rank], 4*nbSendElement[rank]*(NMAX+1), MPI_INT, rank, 3, communicator, &recvRequest[nbRecvRequest]); 345 //ym --> attention taille GloId 346 nbRecvRequest++; 347 } 348 else 349 { 350 MPI_Irecv(recvNeighIds[rank], 4*nbSendElement[rank], MPI_INT, rank, 4, communicator, &recvRequest[nbRecvRequest]); 330 nbSendRequest++; 331 } 332 } 333 if (nbSendElement[rank] > 0) 334 { 335 ep_lib::MPI_Irecv(recvValue[rank], nbSendElement[rank], EP_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]); 336 nbRecvRequest++; 337 ep_lib::MPI_Irecv(recvArea[rank], nbSendElement[rank], EP_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]); 338 nbRecvRequest++; 339 ep_lib::MPI_Irecv(recvGivenArea[rank], nbSendElement[rank], EP_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]); 340 nbRecvRequest++; 341 if (order == 2) 342 { 343 ep_lib::MPI_Irecv(recvGrad[rank], 3*nbSendElement[rank]*(NMAX+1), 344 EP_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]); 345 nbRecvRequest++; 346 ep_lib::MPI_Irecv(recvNeighIds[rank], 4*nbSendElement[rank]*(NMAX+1), EP_INT, rank, 0, communicator, &recvRequest[nbRecvRequest]); 351 347 //ym --> attention taille GloId 352 348 nbRecvRequest++; 353 349 } 350 else 351 { 352 ep_lib::MPI_Irecv(recvNeighIds[rank], 4*nbSendElement[rank], EP_INT, rank, 0, communicator, &recvRequest[nbRecvRequest]); 353 //ym --> attention taille GloId 354 nbRecvRequest++; 355 } 354 356 } 355 357 } 356 358 357 MPI_Waitall(nbSendRequest, sendRequest, status);358 MPI_Waitall(nbRecvRequest, recvRequest, status);359 359 ep_lib::MPI_Waitall(nbSendRequest, sendRequest, status); 360 ep_lib::MPI_Waitall(nbRecvRequest, recvRequest, status); 361 360 362 361 363 /* now that all values and gradients are available use them to computed interpolated values on target 362 and also to compute weights */364 and also to compute weights */ 363 365 int i = 0; 364 366 for (int j = 0; j < nbElements; j++) … … 367 369 368 370 /* since for the 2nd order case source grid elements can contribute to a destination grid element over several "paths" 369 (step1: gradient is computed using neighbours on same grid, step2: intersection uses several elements on other grid)370 accumulate them so that there is only one final weight between two elements */371 (step1: gradient is computed using neighbours on same grid, step2: intersection uses several elements on other grid) 372 accumulate them so that there is only one final weight between two elements */ 371 373 map<GloId,double> wgt_map; 372 374 … … 380 382 double fk = recvValue[rank][n1]; 381 383 double srcArea = recvArea[rank][n1]; 384 double srcGivenArea = recvGivenArea[rank][n1]; 382 385 double w = (*it)->area; 383 386 if (quantity) w/=srcArea ; 387 else w=w*srcGivenArea/srcArea*e.area/e.given_area ; 384 388 385 389 /* first order: src value times weight (weight = supermesh area), later divide by target area */ … … 402 406 double renorm=0; 403 407 if (renormalize) 404 for (map<GloId,double>::iterator it = wgt_map.begin(); it != wgt_map.end(); it++) renorm+=it->second / e.area; 408 { 409 if (quantity) for (map<GloId,double>::iterator it = wgt_map.begin(); it != wgt_map.end(); it++) renorm+=it->second ; 410 else for (map<GloId,double>::iterator it = wgt_map.begin(); it != wgt_map.end(); it++) renorm+=it->second / e.area; 411 } 405 412 else renorm=1. ; 406 413 … … 417 424 } 418 425 } 419 420 //MPI_Barrier(communicator);421 426 422 427 /* free all memory allocated in this function */ … … 428 433 delete[] recvValue[rank]; 429 434 delete[] recvArea[rank]; 435 delete[] recvGivenArea[rank]; 430 436 if (order == 2) 431 437 { … … 439 445 delete[] sendValue[rank]; 440 446 delete[] sendArea[rank]; 447 delete[] sendGivenArea[rank]; 441 448 if (order == 2) 442 449 delete[] sendGrad[rank]; … … 463 470 void Mapper::computeGrads() 464 471 { 465 466 467 468 469 470 471 472 473 474 472 /* array of pointers to collect local elements and elements received from other cpu */ 473 vector<Elt*> globalElements(sstree.nbLocalElements + nbNeighbourElements); 474 int index = 0; 475 for (int i = 0; i < sstree.nbLocalElements; i++, index++) 476 globalElements[index] = &(sstree.localElements[i]); 477 for (int i = 0; i < nbNeighbourElements; i++, index++) 478 globalElements[index] = &neighbourElements[i]; 479 480 update_baryc(sstree.localElements, sstree.nbLocalElements); 481 computeGradients(&globalElements[0], sstree.nbLocalElements); 475 482 } 476 483 … … 479 486 void Mapper::buildMeshTopology() 480 487 { 481 int mpiSize, mpiRank; 482 MPI_Comm_size(communicator, &mpiSize); 483 MPI_Comm_rank(communicator, &mpiRank); 484 485 vector<Node> *routingList = new vector<Node>[mpiSize]; 486 vector<vector<int> > routes(sstree.localTree.leafs.size()); 487 488 sstree.routeIntersections(routes, sstree.localTree.leafs); 489 490 for (int i = 0; i < routes.size(); ++i) 491 for (int k = 0; k < routes[i].size(); ++k) 492 routingList[routes[i][k]].push_back(sstree.localTree.leafs[i]); 493 routingList[mpiRank].clear(); 494 495 496 CMPIRouting mpiRoute(communicator); 497 mpiRoute.init(routes); 498 int nRecv = mpiRoute.getTotalSourceElement(); 499 500 int *nbSendNode = new int[mpiSize]; 501 int *nbRecvNode = new int[mpiSize]; 502 int *sendMessageSize = new int[mpiSize]; 503 int *recvMessageSize = new int[mpiSize]; 504 505 for (int rank = 0; rank < mpiSize; rank++) 506 { 507 nbSendNode[rank] = routingList[rank].size(); 508 sendMessageSize[rank] = 0; 509 for (size_t j = 0; j < routingList[rank].size(); j++) 510 { 511 Elt *elt = (Elt *) (routingList[rank][j].data); 512 sendMessageSize[rank] += packedPolygonSize(*elt); 513 } 514 } 515 516 MPI_Alltoall(nbSendNode, 1, MPI_INT, nbRecvNode, 1, MPI_INT, communicator); 517 MPI_Alltoall(sendMessageSize, 1, MPI_INT, recvMessageSize, 1, MPI_INT, communicator); 518 519 char **sendBuffer = new char*[mpiSize]; 520 char **recvBuffer = new char*[mpiSize]; 521 int *pos = new int[mpiSize]; 522 523 for (int rank = 0; rank < mpiSize; rank++) 524 { 525 if (nbSendNode[rank] > 0) sendBuffer[rank] = new char[sendMessageSize[rank]]; 526 if (nbRecvNode[rank] > 0) recvBuffer[rank] = new char[recvMessageSize[rank]]; 527 } 528 529 for (int rank = 0; rank < mpiSize; rank++) 530 { 531 pos[rank] = 0; 532 for (size_t j = 0; j < routingList[rank].size(); j++) 533 { 534 Elt *elt = (Elt *) (routingList[rank][j].data); 535 packPolygon(*elt, sendBuffer[rank], pos[rank]); 536 } 537 } 538 delete [] routingList; 539 540 541 int nbSendRequest = 0; 542 int nbRecvRequest = 0; 543 MPI_Request *sendRequest = new MPI_Request[mpiSize]; 544 MPI_Request *recvRequest = new MPI_Request[mpiSize]; 545 MPI_Status *status = new MPI_Status[mpiSize]; 546 547 for (int rank = 0; rank < mpiSize; rank++) 548 { 549 if (nbSendNode[rank] > 0) 550 { 551 MPI_Issend(sendBuffer[rank], sendMessageSize[rank], MPI_CHAR, rank, 0, communicator, &sendRequest[nbSendRequest]); 552 nbSendRequest++; 553 } 554 if (nbRecvNode[rank] > 0) 555 { 556 MPI_Irecv(recvBuffer[rank], recvMessageSize[rank], MPI_CHAR, rank, 0, communicator, &recvRequest[nbRecvRequest]); 557 nbRecvRequest++; 558 } 559 } 560 561 MPI_Waitall(nbRecvRequest, recvRequest, status); 562 MPI_Waitall(nbSendRequest, sendRequest, status); 563 564 for (int rank = 0; rank < mpiSize; rank++) 565 if (nbSendNode[rank] > 0) delete [] sendBuffer[rank]; 566 delete [] sendBuffer; 567 568 char **sendBuffer2 = new char*[mpiSize]; 569 char **recvBuffer2 = new char*[mpiSize]; 570 571 for (int rank = 0; rank < mpiSize; rank++) 572 { 573 nbSendNode[rank] = 0; 574 sendMessageSize[rank] = 0; 575 576 if (nbRecvNode[rank] > 0) 577 { 578 set<NodePtr> neighbourList; 579 pos[rank] = 0; 580 for (int j = 0; j < nbRecvNode[rank]; j++) 581 { 582 Elt elt; 583 unpackPolygon(elt, recvBuffer[rank], pos[rank]); 584 Node node(elt.x, cptRadius(elt), &elt); 585 findNeighbour(sstree.localTree.root, &node, neighbourList); 586 } 587 nbSendNode[rank] = neighbourList.size(); 588 for (set<NodePtr>::iterator it = neighbourList.begin(); it != neighbourList.end(); it++) 589 { 590 Elt *elt = (Elt *) ((*it)->data); 591 sendMessageSize[rank] += packedPolygonSize(*elt); 592 } 593 594 sendBuffer2[rank] = new char[sendMessageSize[rank]]; 595 pos[rank] = 0; 596 597 for (set<NodePtr>::iterator it = neighbourList.begin(); it != neighbourList.end(); it++) 598 { 599 Elt *elt = (Elt *) ((*it)->data); 600 packPolygon(*elt, sendBuffer2[rank], pos[rank]); 601 } 602 } 603 } 604 for (int rank = 0; rank < mpiSize; rank++) 605 if (nbRecvNode[rank] > 0) delete [] recvBuffer[rank]; 606 delete [] recvBuffer; 607 608 609 MPI_Barrier(communicator); 610 MPI_Alltoall(nbSendNode, 1, MPI_INT, nbRecvNode, 1, MPI_INT, communicator); 611 MPI_Alltoall(sendMessageSize, 1, MPI_INT, recvMessageSize, 1, MPI_INT, communicator); 612 613 for (int rank = 0; rank < mpiSize; rank++) 614 if (nbRecvNode[rank] > 0) recvBuffer2[rank] = new char[recvMessageSize[rank]]; 615 616 nbSendRequest = 0; 617 nbRecvRequest = 0; 618 619 for (int rank = 0; rank < mpiSize; rank++) 620 { 621 if (nbSendNode[rank] > 0) 622 { 623 MPI_Issend(sendBuffer2[rank], sendMessageSize[rank], MPI_CHAR, rank, 0, communicator, &sendRequest[nbSendRequest]); 624 nbSendRequest++; 625 } 626 if (nbRecvNode[rank] > 0) 627 { 628 MPI_Irecv(recvBuffer2[rank], recvMessageSize[rank], MPI_CHAR, rank, 0, communicator, &recvRequest[nbRecvRequest]); 629 nbRecvRequest++; 630 } 631 } 632 633 MPI_Waitall(nbRecvRequest, recvRequest, status); 634 MPI_Waitall(nbSendRequest, sendRequest, status); 635 636 int nbNeighbourNodes = 0; 637 for (int rank = 0; rank < mpiSize; rank++) 638 nbNeighbourNodes += nbRecvNode[rank]; 639 640 neighbourElements = new Elt[nbNeighbourNodes]; 641 nbNeighbourElements = nbNeighbourNodes; 642 643 int index = 0; 644 for (int rank = 0; rank < mpiSize; rank++) 645 { 646 pos[rank] = 0; 647 for (int j = 0; j < nbRecvNode[rank]; j++) 648 { 649 unpackPolygon(neighbourElements[index], recvBuffer2[rank], pos[rank]); 650 neighbourElements[index].id.ind = sstree.localTree.leafs.size() + index; 651 index++; 652 } 653 } 654 for (int rank = 0; rank < mpiSize; rank++) 655 { 656 if (nbRecvNode[rank] > 0) delete [] recvBuffer2[rank]; 657 if (nbSendNode[rank] > 0) delete [] sendBuffer2[rank]; 658 } 659 delete [] recvBuffer2; 660 delete [] sendBuffer2; 661 delete [] sendMessageSize; 662 delete [] recvMessageSize; 663 delete [] nbSendNode; 664 delete [] nbRecvNode; 665 delete [] sendRequest; 666 delete [] recvRequest; 667 delete [] status; 668 delete [] pos; 669 670 /* re-compute on received elements to avoid having to send this information */ 671 neighbourNodes.resize(nbNeighbourNodes); 672 setCirclesAndLinks(neighbourElements, neighbourNodes); 673 cptAllEltsGeom(neighbourElements, nbNeighbourNodes, srcGrid.pole); 674 675 /* the local SS tree must include nodes from other cpus if they are potential 488 int mpiSize, mpiRank; 489 ep_lib::MPI_Comm_size(communicator, &mpiSize); 490 ep_lib::MPI_Comm_rank(communicator, &mpiRank); 491 492 vector<Node> *routingList = new vector<Node>[mpiSize]; 493 vector<vector<int> > routes(sstree.localTree.leafs.size()); 494 495 sstree.routeIntersections(routes, sstree.localTree.leafs); 496 497 for (int i = 0; i < routes.size(); ++i) 498 for (int k = 0; k < routes[i].size(); ++k) 499 routingList[routes[i][k]].push_back(sstree.localTree.leafs[i]); 500 routingList[mpiRank].clear(); 501 502 503 CMPIRouting mpiRoute(communicator); 504 mpiRoute.init(routes); 505 int nRecv = mpiRoute.getTotalSourceElement(); 506 // cout << mpiRank << " NRECV " << nRecv << "(" << routes.size() << ")"<< endl; 507 508 int *nbSendNode = new int[mpiSize]; 509 int *nbRecvNode = new int[mpiSize]; 510 int *sendMessageSize = new int[mpiSize]; 511 int *recvMessageSize = new int[mpiSize]; 512 513 for (int rank = 0; rank < mpiSize; rank++) 514 { 515 nbSendNode[rank] = routingList[rank].size(); 516 sendMessageSize[rank] = 0; 517 for (size_t j = 0; j < routingList[rank].size(); j++) 518 { 519 Elt *elt = (Elt *) (routingList[rank][j].data); 520 sendMessageSize[rank] += packedPolygonSize(*elt); 521 } 522 } 523 524 ep_lib::MPI_Alltoall(nbSendNode, 1, EP_INT, nbRecvNode, 1, EP_INT, communicator); 525 ep_lib::MPI_Alltoall(sendMessageSize, 1, EP_INT, recvMessageSize, 1, EP_INT, communicator); 526 527 char **sendBuffer = new char*[mpiSize]; 528 char **recvBuffer = new char*[mpiSize]; 529 int *pos = new int[mpiSize]; 530 531 for (int rank = 0; rank < mpiSize; rank++) 532 { 533 if (nbSendNode[rank] > 0) sendBuffer[rank] = new char[sendMessageSize[rank]]; 534 if (nbRecvNode[rank] > 0) recvBuffer[rank] = new char[recvMessageSize[rank]]; 535 } 536 537 for (int rank = 0; rank < mpiSize; rank++) 538 { 539 pos[rank] = 0; 540 for (size_t j = 0; j < routingList[rank].size(); j++) 541 { 542 Elt *elt = (Elt *) (routingList[rank][j].data); 543 packPolygon(*elt, sendBuffer[rank], pos[rank]); 544 } 545 } 546 delete [] routingList; 547 548 549 int nbSendRequest = 0; 550 int nbRecvRequest = 0; 551 ep_lib::MPI_Request *sendRequest = new ep_lib::MPI_Request[mpiSize]; 552 ep_lib::MPI_Request *recvRequest = new ep_lib::MPI_Request[mpiSize]; 553 ep_lib::MPI_Status *status = new ep_lib::MPI_Status[mpiSize]; 554 555 for (int rank = 0; rank < mpiSize; rank++) 556 { 557 if (nbSendNode[rank] > 0) 558 { 559 ep_lib::MPI_Issend(sendBuffer[rank], sendMessageSize[rank], EP_CHAR, rank, 0, communicator, &sendRequest[nbSendRequest]); 560 nbSendRequest++; 561 } 562 if (nbRecvNode[rank] > 0) 563 { 564 ep_lib::MPI_Irecv(recvBuffer[rank], recvMessageSize[rank], EP_CHAR, rank, 0, communicator, &recvRequest[nbRecvRequest]); 565 nbRecvRequest++; 566 } 567 } 568 569 ep_lib::MPI_Waitall(nbRecvRequest, recvRequest, status); 570 ep_lib::MPI_Waitall(nbSendRequest, sendRequest, status); 571 572 for (int rank = 0; rank < mpiSize; rank++) 573 if (nbSendNode[rank] > 0) delete [] sendBuffer[rank]; 574 delete [] sendBuffer; 575 576 char **sendBuffer2 = new char*[mpiSize]; 577 char **recvBuffer2 = new char*[mpiSize]; 578 579 for (int rank = 0; rank < mpiSize; rank++) 580 { 581 nbSendNode[rank] = 0; 582 sendMessageSize[rank] = 0; 583 584 if (nbRecvNode[rank] > 0) 585 { 586 set<NodePtr> neighbourList; 587 pos[rank] = 0; 588 for (int j = 0; j < nbRecvNode[rank]; j++) 589 { 590 Elt elt; 591 unpackPolygon(elt, recvBuffer[rank], pos[rank]); 592 Node node(elt.x, cptRadius(elt), &elt); 593 findNeighbour(sstree.localTree.root, &node, neighbourList); 594 } 595 nbSendNode[rank] = neighbourList.size(); 596 for (set<NodePtr>::iterator it = neighbourList.begin(); it != neighbourList.end(); it++) 597 { 598 Elt *elt = (Elt *) ((*it)->data); 599 sendMessageSize[rank] += packedPolygonSize(*elt); 600 } 601 602 sendBuffer2[rank] = new char[sendMessageSize[rank]]; 603 pos[rank] = 0; 604 605 for (set<NodePtr>::iterator it = neighbourList.begin(); it != neighbourList.end(); it++) 606 { 607 Elt *elt = (Elt *) ((*it)->data); 608 packPolygon(*elt, sendBuffer2[rank], pos[rank]); 609 } 610 } 611 } 612 for (int rank = 0; rank < mpiSize; rank++) 613 if (nbRecvNode[rank] > 0) delete [] recvBuffer[rank]; 614 delete [] recvBuffer; 615 616 617 ep_lib::MPI_Barrier(communicator); 618 ep_lib::MPI_Alltoall(nbSendNode, 1, EP_INT, nbRecvNode, 1, EP_INT, communicator); 619 ep_lib::MPI_Alltoall(sendMessageSize, 1, EP_INT, recvMessageSize, 1, EP_INT, communicator); 620 621 for (int rank = 0; rank < mpiSize; rank++) 622 if (nbRecvNode[rank] > 0) recvBuffer2[rank] = new char[recvMessageSize[rank]]; 623 624 nbSendRequest = 0; 625 nbRecvRequest = 0; 626 627 for (int rank = 0; rank < mpiSize; rank++) 628 { 629 if (nbSendNode[rank] > 0) 630 { 631 ep_lib::MPI_Issend(sendBuffer2[rank], sendMessageSize[rank], EP_CHAR, rank, 0, communicator, &sendRequest[nbSendRequest]); 632 nbSendRequest++; 633 } 634 if (nbRecvNode[rank] > 0) 635 { 636 ep_lib::MPI_Irecv(recvBuffer2[rank], recvMessageSize[rank], EP_CHAR, rank, 0, communicator, &recvRequest[nbRecvRequest]); 637 nbRecvRequest++; 638 } 639 } 640 641 ep_lib::MPI_Waitall(nbRecvRequest, recvRequest, status); 642 ep_lib::MPI_Waitall(nbSendRequest, sendRequest, status); 643 644 int nbNeighbourNodes = 0; 645 for (int rank = 0; rank < mpiSize; rank++) 646 nbNeighbourNodes += nbRecvNode[rank]; 647 648 neighbourElements = new Elt[nbNeighbourNodes]; 649 nbNeighbourElements = nbNeighbourNodes; 650 651 int index = 0; 652 for (int rank = 0; rank < mpiSize; rank++) 653 { 654 pos[rank] = 0; 655 for (int j = 0; j < nbRecvNode[rank]; j++) 656 { 657 unpackPolygon(neighbourElements[index], recvBuffer2[rank], pos[rank]); 658 neighbourElements[index].id.ind = sstree.localTree.leafs.size() + index; 659 index++; 660 } 661 } 662 for (int rank = 0; rank < mpiSize; rank++) 663 { 664 if (nbRecvNode[rank] > 0) delete [] recvBuffer2[rank]; 665 if (nbSendNode[rank] > 0) delete [] sendBuffer2[rank]; 666 } 667 delete [] recvBuffer2; 668 delete [] sendBuffer2; 669 delete [] sendMessageSize; 670 delete [] recvMessageSize; 671 delete [] nbSendNode; 672 delete [] nbRecvNode; 673 delete [] sendRequest; 674 delete [] recvRequest; 675 delete [] status; 676 delete [] pos; 677 678 /* re-compute on received elements to avoid having to send this information */ 679 neighbourNodes.resize(nbNeighbourNodes); 680 setCirclesAndLinks(neighbourElements, neighbourNodes); 681 cptAllEltsGeom(neighbourElements, nbNeighbourNodes, srcGrid.pole); 682 683 /* the local SS tree must include nodes from other cpus if they are potential 676 684 intersector of a local node */ 677 678 679 685 sstree.localTree.insertNodes(neighbourNodes); 686 687 /* for every local element, 680 688 use the SS-tree to find all elements (including neighbourElements) 681 689 who are potential neighbours because their circles intersect, 682 683 684 685 686 687 688 689 690 691 692 693 694 695 690 then check all canditates for common edges to build up connectivity information 691 */ 692 for (int j = 0; j < sstree.localTree.leafs.size(); j++) 693 { 694 Node& node = sstree.localTree.leafs[j]; 695 696 /* find all leafs whoes circles that intersect node's circle and save into node->intersectors */ 697 node.search(sstree.localTree.root); 698 699 Elt *elt = (Elt *)(node.data); 700 701 for (int i = 0; i < elt->n; i++) elt->neighbour[i] = NOT_FOUND; 702 703 /* for element `elt` loop through all nodes in the SS-tree 696 704 whoes circles intersect with the circle around `elt` (the SS intersectors) 697 705 and check if they are neighbours in the sense that the two elements share an edge. 698 706 If they do, save this information for elt */ 699 700 701 702 703 707 for (list<NodePtr>::iterator it = (node.intersectors).begin(); it != (node.intersectors).end(); ++it) 708 { 709 Elt *elt2 = (Elt *)((*it)->data); 710 set_neighbour(*elt, *elt2); 711 } 704 712 705 713 /* 706 707 708 709 710 714 for (int i = 0; i < elt->n; i++) 715 { 716 if (elt->neighbour[i] == NOT_FOUND) 717 error_exit("neighbour not found"); 718 } 711 719 */ 712 720 } 713 721 } 714 722 … … 716 724 void Mapper::computeIntersection(Elt *elements, int nbElements) 717 725 { 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 MPI_Alltoall(nbSendNode, 1, MPI_INT, nbRecvNode, 1, MPI_INT, communicator);766 MPI_Alltoall(sentMessageSize, 1, MPI_INT, recvMessageSize, 1, MPI_INT, communicator);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 MPI_Request *sendRequest = newMPI_Request[mpiSize];800 MPI_Request *recvRequest = newMPI_Request[mpiSize];801 MPI_Status *status = newMPI_Status[mpiSize];802 803 804 805 806 807 MPI_Issend(sendBuffer[rank], sentMessageSize[rank], MPI_CHAR, rank, 0, communicator, &sendRequest[nbSendRequest]);808 809 810 811 812 MPI_Irecv(recvBuffer[rank], recvMessageSize[rank], MPI_CHAR, rank, 0, communicator, &recvRequest[nbRecvRequest]);813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 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 MPI_Alltoall(sentMessageSize, 1, MPI_INT, recvMessageSize, 1, MPI_INT, communicator);878 879 880 881 882 883 884 885 886 887 888 889 890 MPI_Issend(sendBuffer2[rank], sentMessageSize[rank], MPI_CHAR, rank, 0, communicator, &sendRequest[nbSendRequest]);891 892 893 894 895 MPI_Irecv(recvBuffer2[rank], recvMessageSize[rank], MPI_CHAR, rank, 0, communicator, &recvRequest[nbRecvRequest]);896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 916 917 918 919 920 921 922 923 924 925 926 927 928 726 int mpiSize, mpiRank; 727 ep_lib::MPI_Comm_size(communicator, &mpiSize); 728 ep_lib::MPI_Comm_rank(communicator, &mpiRank); 729 730 ep_lib::MPI_Barrier(communicator); 731 732 vector<Node> *routingList = new vector<Node>[mpiSize]; 733 734 vector<Node> routeNodes; routeNodes.reserve(nbElements); 735 for (int j = 0; j < nbElements; j++) 736 { 737 elements[j].id.ind = j; 738 elements[j].id.rank = mpiRank; 739 routeNodes.push_back(Node(elements[j].x, cptRadius(elements[j]), &elements[j])); 740 } 741 742 vector<vector<int> > routes(routeNodes.size()); 743 sstree.routeIntersections(routes, routeNodes); 744 for (int i = 0; i < routes.size(); ++i) 745 for (int k = 0; k < routes[i].size(); ++k) 746 routingList[routes[i][k]].push_back(routeNodes[i]); 747 748 if (verbose >= 2) 749 { 750 cout << " --> rank " << mpiRank << " nbElements " << nbElements << " : "; 751 for (int rank = 0; rank < mpiSize; rank++) 752 cout << routingList[rank].size() << " "; 753 cout << endl; 754 } 755 ep_lib::MPI_Barrier(communicator); 756 757 int *nbSendNode = new int[mpiSize]; 758 int *nbRecvNode = new int[mpiSize]; 759 int *sentMessageSize = new int[mpiSize]; 760 int *recvMessageSize = new int[mpiSize]; 761 762 for (int rank = 0; rank < mpiSize; rank++) 763 { 764 nbSendNode[rank] = routingList[rank].size(); 765 sentMessageSize[rank] = 0; 766 for (size_t j = 0; j < routingList[rank].size(); j++) 767 { 768 Elt *elt = (Elt *) (routingList[rank][j].data); 769 sentMessageSize[rank] += packedPolygonSize(*elt); 770 } 771 } 772 773 ep_lib::MPI_Alltoall(nbSendNode, 1, EP_INT, nbRecvNode, 1, EP_INT, communicator); 774 ep_lib::MPI_Alltoall(sentMessageSize, 1, EP_INT, recvMessageSize, 1, EP_INT, communicator); 775 776 int total = 0; 777 778 for (int rank = 0; rank < mpiSize; rank++) 779 { 780 total = total + nbRecvNode[rank]; 781 } 782 783 if (verbose >= 2) cout << "---> rank " << mpiRank << " : compute intersection : total received nodes " << total << endl; 784 char **sendBuffer = new char*[mpiSize]; 785 char **recvBuffer = new char*[mpiSize]; 786 int *pos = new int[mpiSize]; 787 788 for (int rank = 0; rank < mpiSize; rank++) 789 { 790 if (nbSendNode[rank] > 0) sendBuffer[rank] = new char[sentMessageSize[rank]]; 791 if (nbRecvNode[rank] > 0) recvBuffer[rank] = new char[recvMessageSize[rank]]; 792 } 793 794 for (int rank = 0; rank < mpiSize; rank++) 795 { 796 pos[rank] = 0; 797 for (size_t j = 0; j < routingList[rank].size(); j++) 798 { 799 Elt* elt = (Elt *) (routingList[rank][j].data); 800 packPolygon(*elt, sendBuffer[rank], pos[rank]); 801 } 802 } 803 delete [] routingList; 804 805 int nbSendRequest = 0; 806 int nbRecvRequest = 0; 807 ep_lib::MPI_Request *sendRequest = new ep_lib::MPI_Request[mpiSize]; 808 ep_lib::MPI_Request *recvRequest = new ep_lib::MPI_Request[mpiSize]; 809 ep_lib::MPI_Status *status = new ep_lib::MPI_Status[mpiSize]; 810 811 for (int rank = 0; rank < mpiSize; rank++) 812 { 813 if (nbSendNode[rank] > 0) 814 { 815 ep_lib::MPI_Issend(sendBuffer[rank], sentMessageSize[rank], EP_CHAR, rank, 0, communicator, &sendRequest[nbSendRequest]); 816 nbSendRequest++; 817 } 818 if (nbRecvNode[rank] > 0) 819 { 820 ep_lib::MPI_Irecv(recvBuffer[rank], recvMessageSize[rank], EP_CHAR, rank, 0, communicator, &recvRequest[nbRecvRequest]); 821 nbRecvRequest++; 822 } 823 } 824 825 ep_lib::MPI_Waitall(nbRecvRequest, recvRequest, status); 826 ep_lib::MPI_Waitall(nbSendRequest, sendRequest, status); 827 char **sendBuffer2 = new char*[mpiSize]; 828 char **recvBuffer2 = new char*[mpiSize]; 829 830 double tic = cputime(); 831 for (int rank = 0; rank < mpiSize; rank++) 832 { 833 sentMessageSize[rank] = 0; 834 835 if (nbRecvNode[rank] > 0) 836 { 837 Elt *recvElt = new Elt[nbRecvNode[rank]]; 838 pos[rank] = 0; 839 for (int j = 0; j < nbRecvNode[rank]; j++) 840 { 841 unpackPolygon(recvElt[j], recvBuffer[rank], pos[rank]); 842 cptEltGeom(recvElt[j], tgtGrid.pole); 843 Node recvNode(recvElt[j].x, cptRadius(recvElt[j]), &recvElt[j]); 844 recvNode.search(sstree.localTree.root); 845 /* for a node holding an element of the target, loop throught candidates for intersecting source */ 846 for (list<NodePtr>::iterator it = (recvNode.intersectors).begin(); it != (recvNode.intersectors).end(); ++it) 847 { 848 Elt *elt2 = (Elt *) ((*it)->data); 849 /* recvElt is target, elt2 is source */ 850 // intersect(&recvElt[j], elt2); 851 intersect_ym(&recvElt[j], elt2); 852 } 853 854 if (recvElt[j].is.size() > 0) sentMessageSize[rank] += packIntersectionSize(recvElt[j]); 855 856 // here recvNode goes out of scope 857 } 858 859 if (sentMessageSize[rank] > 0) 860 { 861 sentMessageSize[rank] += sizeof(int); 862 sendBuffer2[rank] = new char[sentMessageSize[rank]]; 863 *((int *) sendBuffer2[rank]) = 0; 864 pos[rank] = sizeof(int); 865 for (int j = 0; j < nbRecvNode[rank]; j++) 866 { 867 packIntersection(recvElt[j], sendBuffer2[rank], pos[rank]); 868 //FIXME should be deleted: recvElt[j].delete_intersections(); // intersection areas have been packed to buffer and won't be used any more 869 } 870 } 871 delete [] recvElt; 872 873 } 874 } 875 delete [] pos; 876 877 for (int rank = 0; rank < mpiSize; rank++) 878 { 879 if (nbSendNode[rank] > 0) delete [] sendBuffer[rank]; 880 if (nbRecvNode[rank] > 0) delete [] recvBuffer[rank]; 881 nbSendNode[rank] = 0; 882 } 883 884 if (verbose >= 2) cout << "Rank " << mpiRank << " Compute (internal) intersection " << cputime() - tic << " s" << endl; 885 ep_lib::MPI_Alltoall(sentMessageSize, 1, EP_INT, recvMessageSize, 1, EP_INT, communicator); 886 887 for (int rank = 0; rank < mpiSize; rank++) 888 if (recvMessageSize[rank] > 0) 889 recvBuffer2[rank] = new char[recvMessageSize[rank]]; 890 891 nbSendRequest = 0; 892 nbRecvRequest = 0; 893 894 for (int rank = 0; rank < mpiSize; rank++) 895 { 896 if (sentMessageSize[rank] > 0) 897 { 898 ep_lib::MPI_Issend(sendBuffer2[rank], sentMessageSize[rank], EP_CHAR, rank, 0, communicator, &sendRequest[nbSendRequest]); 899 nbSendRequest++; 900 } 901 if (recvMessageSize[rank] > 0) 902 { 903 ep_lib::MPI_Irecv(recvBuffer2[rank], recvMessageSize[rank], EP_CHAR, rank, 0, communicator, &recvRequest[nbRecvRequest]); 904 nbRecvRequest++; 905 } 906 } 907 908 ep_lib::MPI_Waitall(nbRecvRequest, recvRequest, status); 909 ep_lib::MPI_Waitall(nbSendRequest, sendRequest, status); 910 911 delete [] sendRequest; 912 delete [] recvRequest; 913 delete [] status; 914 for (int rank = 0; rank < mpiSize; rank++) 915 { 916 if (nbRecvNode[rank] > 0) 917 { 918 if (sentMessageSize[rank] > 0) 919 delete [] sendBuffer2[rank]; 920 } 921 922 if (recvMessageSize[rank] > 0) 923 { 924 unpackIntersection(elements, recvBuffer2[rank]); 925 delete [] recvBuffer2[rank]; 926 } 927 } 928 delete [] sendBuffer2; 929 delete [] recvBuffer2; 930 delete [] sendBuffer; 931 delete [] recvBuffer; 932 933 delete [] nbSendNode; 934 delete [] nbRecvNode; 935 delete [] sentMessageSize; 936 delete [] recvMessageSize; 929 937 } 930 938 931 939 Mapper::~Mapper() 932 940 { 933 934 935 936 937 938 } 939 940 } 941 delete [] remapMatrix; 942 delete [] srcAddress; 943 delete [] srcRank; 944 delete [] dstAddress; 945 if (neighbourElements) delete [] neighbourElements; 946 } 947 948 } -
XIOS/dev/branch_openmp/extern/remap/src/mapper.hpp
r1538 r1642 18 18 { 19 19 public: 20 Mapper(ep_lib::MPI_Comm comm) : communicator(comm), verbose(SILENT), neighbourElements(NULL), sstree(comm) {} 21 20 Mapper(ep_lib::MPI_Comm comm=EP_COMM_WORLD) : communicator(comm), verbose(SILENT), neighbourElements(NULL), sstree(comm) {} 22 21 ~Mapper(); 23 22 void setVerbosity(verbosity v) {verbose=v ;} 24 23 25 void setSourceMesh(const double* boundsLon, const double* boundsLat, int nVertex, int nbCells, const double* pole, const long int* globalId=NULL) ;26 void setTargetMesh(const double* boundsLon, const double* boundsLat, int nVertex, int nbCells, const double* pole, const long int* globalId=NULL) ;24 void setSourceMesh(const double* boundsLon, const double* boundsLat, const double* area, int nVertex, int nbCells, const double* pole, const long int* globalId=NULL) ; 25 void setTargetMesh(const double* boundsLon, const double* boundsLat, const double* area, int nVertex, int nbCells, const double* pole, const long int* globalId=NULL) ; 27 26 void setSourceValue(const double* val) ; 28 27 void getTargetValue(double* val) ; -
XIOS/dev/branch_openmp/extern/remap/src/meshutil.cpp
r877 r1642 2 2 #include "elt.hpp" 3 3 #include "polyg.hpp" 4 #include "intersection_ym.hpp" 5 #include "earcut.hpp" 6 #include <vector> 4 7 5 8 namespace sphereRemap { 6 9 7 10 using namespace std; 11 12 double computePolygoneArea(Elt& a, const Coord &pole) 13 { 14 using N = uint32_t; 15 using Point = array<double, 2>; 16 vector<Point> vect_points; 17 vector< vector<Point> > polyline; 18 19 vector<Coord> dstPolygon ; 20 createGreatCirclePolygon(a, pole, dstPolygon) ; 21 22 int na=dstPolygon.size() ; 23 Coord *a_gno = new Coord[na]; 24 25 Coord OC=barycentre(a.vertex,a.n) ; 26 Coord Oz=OC ; 27 Coord Ox=crossprod(Coord(0,0,1),Oz) ; 28 // choose Ox not too small to avoid rounding error 29 if (norm(Ox)< 0.1) Ox=crossprod(Coord(0,1,0),Oz) ; 30 Ox=Ox*(1./norm(Ox)) ; 31 Coord Oy=crossprod(Oz,Ox) ; 32 double cos_alpha; 33 34 for(int n=0; n<na;n++) 35 { 36 cos_alpha=scalarprod(OC,dstPolygon[n]) ; 37 a_gno[n].x=scalarprod(dstPolygon[n],Ox)/cos_alpha ; 38 a_gno[n].y=scalarprod(dstPolygon[n],Oy)/cos_alpha ; 39 a_gno[n].z=scalarprod(dstPolygon[n],Oz)/cos_alpha ; // must be equal to 1 40 41 vect_points.push_back( array<double, 2>() ); 42 vect_points[n][0] = a_gno[n].x; 43 vect_points[n][1] = a_gno[n].y; 44 45 } 46 47 polyline.push_back(vect_points); 48 vector<N> indices_a_gno = mapbox::earcut<N>(polyline); 49 50 double area_a_gno=0 ; 51 for(int i=0;i<indices_a_gno.size()/3;++i) 52 { 53 Coord x0 = Ox * polyline[0][indices_a_gno[3*i]][0] + Oy* polyline[0][indices_a_gno[3*i]][1] + Oz ; 54 Coord x1 = Ox * polyline[0][indices_a_gno[3*i+1]][0] + Oy* polyline[0][indices_a_gno[3*i+1]][1] + Oz ; 55 Coord x2 = Ox * polyline[0][indices_a_gno[3*i+2]][0] + Oy* polyline[0][indices_a_gno[3*i+2]][1] + Oz ; 56 area_a_gno+=triarea(x0 * (1./norm(x0)),x1* (1./norm(x1)), x2* (1./norm(x2))) ; 57 } 58 59 vect_points.clear(); 60 polyline.clear(); 61 indices_a_gno.clear(); 62 return area_a_gno ; 63 } 64 8 65 9 66 void cptEltGeom(Elt& elt, const Coord &pole) … … 14 71 elt.area = airbar(elt.n, elt.vertex, elt.edge, elt.d, pole, gg); 15 72 elt.x = gg; 16 } 73 // overwrite area computation 74 75 elt.area = computePolygoneArea(elt, pole) ; 76 } 77 17 78 18 79 void cptAllEltsGeom(Elt *elt, int N, const Coord &pole) … … 133 194 { 134 195 for (int i = 0; i < elts[j]->n; i++) 135 if ( elts[j]->neighbour[i]== NOT_FOUND) neighbours[i]=NULL ; 196 if ( elts[j]->neighbour[i]== NOT_FOUND) neighbours[i]=NULL ; // no neighbour 197 else if (elts[elts[j]->neighbour[i]]->is.size() == 0) neighbours[i]=NULL ; // neighbour has none supermesh cell 136 198 else neighbours[i] = elts[elts[j]->neighbour[i]]; 137 199 -
XIOS/dev/branch_openmp/extern/remap/src/mpi_cascade.cpp
r1328 r1642 1 1 #include "mpi_cascade.hpp" 2 2 #include <iostream> 3 using namespace ep_lib;4 3 5 4 namespace sphereRemap { 6 5 7 CMPICascade::CMPICascade(int nodes_per_level, MPI_Comm comm)6 CMPICascade::CMPICascade(int nodes_per_level, ep_lib::MPI_Comm comm) 8 7 { 9 8 int remaining_levels; 10 MPI_Comm intraComm;9 ep_lib::MPI_Comm intraComm; 11 10 int l = 0; // current level 12 11 do { … … 16 15 level[l].p_grp_size = level[l].size/level[l].group_size; 17 16 18 MPI_Comm_split(comm, level[l].colour(), level[l].key(), &intraComm);19 MPI_Comm_split(comm, level[l].p_colour(), level[l].p_key(), &(level[l].pg_comm));17 ep_lib::MPI_Comm_split(comm, level[l].colour(), level[l].key(), &intraComm); 18 ep_lib::MPI_Comm_split(comm, level[l].p_colour(), level[l].p_key(), &(level[l].pg_comm)); 20 19 comm = intraComm; 21 20 l++; -
XIOS/dev/branch_openmp/extern/remap/src/mpi_cascade.hpp
r1538 r1642 12 12 { 13 13 public: 14 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 CMPICascade(int nodes_per_level, ep_lib::MPI_Comm comm); 36 // 37 CMPICascade(int nodes_per_level, ep_lib::MPI_Comm comm); 37 38 38 39 39 int num_levels; 40 std::vector<CCascadeLevel> level; 40 41 }; 41 42 -
XIOS/dev/branch_openmp/extern/remap/src/mpi_routing.cpp
r1362 r1642 5 5 #include "timerRemap.hpp" 6 6 #include <iostream> 7 using namespace ep_lib;8 7 9 8 namespace sphereRemap { … … 11 10 const int verbose = 0; 12 11 13 CMPIRouting::CMPIRouting( MPI_Comm comm) : communicator(comm)14 { 15 MPI_Comm_rank(comm, &mpiRank);16 MPI_Comm_size(comm, &mpiSize);12 CMPIRouting::CMPIRouting(ep_lib::MPI_Comm comm) : communicator(comm) 13 { 14 ep_lib::MPI_Comm_rank(comm, &mpiRank); 15 ep_lib::MPI_Comm_size(comm, &mpiSize); 17 16 } 18 17 … … 20 19 but message lengths are *known* to receiver */ 21 20 template <typename T> 22 void alltoalls_known(const vector<vector<T> >& send, vector<vector<T> >& recv, const vector<int>& ranks, MPI_Comm communicator)23 { 24 vector< MPI_Request> request(ranks.size() * 2);25 vector< MPI_Status> status(ranks.size() * 2);21 void alltoalls_known(const vector<vector<T> >& send, vector<vector<T> >& recv, const vector<int>& ranks, ep_lib::MPI_Comm communicator) 22 { 23 vector<ep_lib::MPI_Request> request(ranks.size() * 2); 24 vector<ep_lib::MPI_Status> status(ranks.size() * 2); 26 25 27 26 // communicate data … … 29 28 for (int i = 0; i < ranks.size(); i++) 30 29 if (recv[i].size()) 31 MPI_Irecv(&recv[i][0], recv[i].size()*sizeof(T), MPI_CHAR, ranks[i], 0, communicator, &request[nbRequest++]);30 ep_lib::MPI_Irecv(&recv[i][0], recv[i].size()*sizeof(T), EP_CHAR, ranks[i], 0, communicator, &request[nbRequest++]); 32 31 for (int i = 0; i < ranks.size(); i++) 33 32 if (send[i].size()) 34 MPI_Isend((void *) &send[i][0], send[i].size()*sizeof(T), MPI_CHAR, ranks[i], 0, communicator, &request[nbRequest++]);35 MPI_Waitall(nbRequest, &request[0], &status[0]);33 ep_lib::MPI_Isend((void *) &send[i][0], send[i].size()*sizeof(T), EP_CHAR, ranks[i], 0, communicator, &request[nbRequest++]); 34 ep_lib::MPI_Waitall(nbRequest, &request[0], &status[0]); 36 35 } 37 36 … … 39 38 but message lengths are *unknown* to receiver */ 40 39 template <typename T> 41 void alltoalls_unknown(const vector<vector<T> >& send, vector<vector<T> >& recv, const vector<int>& ranks, MPI_Comm communicator)42 { 43 vector< MPI_Request> request(ranks.size() * 2);44 vector< MPI_Status> status(ranks.size() * 2);40 void alltoalls_unknown(const vector<vector<T> >& send, vector<vector<T> >& recv, const vector<int>& ranks, ep_lib::MPI_Comm communicator) 41 { 42 vector<ep_lib::MPI_Request> request(ranks.size() * 2); 43 vector<ep_lib::MPI_Status> status(ranks.size() * 2); 45 44 46 45 // communicate sizes … … 51 50 sendSizes[i] = send[i].size(); 52 51 for (int i = 0; i < ranks.size(); i++) 53 MPI_Irecv(&recvSizes[i], 1, MPI_INT, ranks[i], 0, communicator, &request[nbRequest++]);52 ep_lib::MPI_Irecv(&recvSizes[i], 1, EP_INT, ranks[i], 0, communicator, &request[nbRequest++]); 54 53 for (int i = 0; i < ranks.size(); i++) 55 MPI_Isend(&sendSizes[i], 1, MPI_INT, ranks[i], 0, communicator, &request[nbRequest++]);56 MPI_Waitall(nbRequest, &request[0], &status[0]);54 ep_lib::MPI_Isend(&sendSizes[i], 1, EP_INT, ranks[i], 0, communicator, &request[nbRequest++]); 55 ep_lib::MPI_Waitall(nbRequest, &request[0], &status[0]); 57 56 58 57 // allocate … … 119 118 CTimer::get("CMPIRouting::init(reduce_scatter)").reset(); 120 119 CTimer::get("CMPIRouting::init(reduce_scatter)").resume(); 121 MPI_Reduce_scatter(toSend, &nbSource, recvCount, MPI_INT, MPI_SUM, communicator);120 ep_lib::MPI_Reduce_scatter(toSend, &nbSource, recvCount, EP_INT, EP_SUM, communicator); 122 121 CTimer::get("CMPIRouting::init(reduce_scatter)").suspend(); 123 122 CTimer::get("CMPIRouting::init(reduce_scatter)").print(); 124 123 125 MPI_Alloc_mem(nbTarget *sizeof(int), MPI_INFO_NULL, &targetRank);126 MPI_Alloc_mem(nbSource *sizeof(int), MPI_INFO_NULL, &sourceRank);124 ep_lib::MPI_Alloc_mem(nbTarget *sizeof(int), EP_INFO_NULL, &targetRank); 125 ep_lib::MPI_Alloc_mem(nbSource *sizeof(int), EP_INFO_NULL, &sourceRank); 127 126 128 127 targetRankToIndex = new int[mpiSize]; … … 138 137 } 139 138 140 MPI_Barrier(communicator);139 ep_lib::MPI_Barrier(communicator); 141 140 CTimer::get("CMPIRouting::init(get_source)").reset(); 142 141 CTimer::get("CMPIRouting::init(get_source)").resume(); 143 142 144 MPI_Request *request = newMPI_Request[nbSource + nbTarget];145 MPI_Status *status = newMPI_Status[nbSource + nbTarget];143 ep_lib::MPI_Request *request = new ep_lib::MPI_Request[nbSource + nbTarget]; 144 ep_lib::MPI_Status *status = new ep_lib::MPI_Status[nbSource + nbTarget]; 146 145 147 146 int indexRequest = 0; … … 151 150 for (int i = 0; i < nbSource; i++) 152 151 { 153 #ifdef _using EP154 MPI_Irecv(&sourceRank[i], 1, MPI_INT, -2, 0, communicator, &request[indexRequest++]);155 #el se156 MPI_Irecv(&sourceRank[i], 1, MPI_INT, MPI_ANY_SOURCE, 0, communicator, &request[indexRequest++]);152 #ifdef _usingMPI 153 ep_lib::MPI_Irecv(&sourceRank[i], 1, EP_INT, MPI_ANY_SOURCE, 0, communicator, &request[indexRequest]); 154 #elif _usingEP 155 ep_lib::MPI_Irecv(&sourceRank[i], 1, EP_INT, -2, 0, communicator, &request[indexRequest]); 157 156 #endif 157 indexRequest++; 158 158 } 159 159 MPI_Barrier(communicator); 160 160 for (int i = 0; i < nbTarget; i++) 161 161 { 162 MPI_Isend(&mpiRank, 1, MPI_INT, targetRank[i], 0, communicator, &request[indexRequest++]); 162 ep_lib::MPI_Isend(&mpiRank, 1, EP_INT, targetRank[i], 0, communicator, &request[indexRequest]); 163 indexRequest++; 163 164 } 164 165 MPI_Waitall(indexRequest, request, status); … … 173 174 for (int i = 0; i < nbSource; i++) 174 175 { 175 #ifdef _using EP176 MPI_Irecv(&sourceRank[i], 1, MPI_INT, -2, 0, communicator, &request[indexRequest]);177 #el se178 MPI_Irecv(&sourceRank[i], 1, MPI_INT, MPI_ANY_SOURCE, 0, communicator, &request[indexRequest]);176 #ifdef _usingMPI 177 ep_lib::MPI_Irecv(&sourceRank[i], 1, EP_INT, MPI_ANY_SOURCE, 0, communicator, &request[indexRequest]); 178 #elif _usingEP 179 ep_lib::MPI_Irecv(&sourceRank[i], 1, EP_INT, -2, 0, communicator, &request[indexRequest]); 179 180 #endif 180 181 } 182 183 for (int i = 0; i < nbTarget; i++) 184 { 185 MPI_Isend(&mpiRank, 1, MPI_INT, targetRank[i], 0, communicator, &request[indexRequest]);181 indexRequest++; 182 } 183 184 for (int i = 0; i < nbTarget; i++) 185 { 186 ep_lib::MPI_Isend(&mpiRank, 1, EP_INT, targetRank[i], 0, communicator, &request[indexRequest]); 186 187 indexRequest++; 187 188 } … … 208 209 for (int i = 0; i < nbSource; i++) 209 210 { 210 MPI_Irecv(&nbSourceElement[i], 1, MPI_INT, sourceRank[i], 0, communicator, &request[indexRequest]);211 ep_lib::MPI_Irecv(&nbSourceElement[i], 1, EP_INT, sourceRank[i], 0, communicator, &request[indexRequest]); 211 212 indexRequest++; 212 213 } … … 215 216 { 216 217 totalTargetElement += nbTargetElement[i]; 217 MPI_Isend(&nbTargetElement[i], 1, MPI_INT, targetRank[i], 0, communicator, &request[indexRequest]);218 ep_lib::MPI_Isend(&nbTargetElement[i], 1, EP_INT, targetRank[i], 0, communicator, &request[indexRequest]); 218 219 indexRequest++; 219 220 } … … 283 284 284 285 285 MPI_Request* request=newMPI_Request[nbSource+nbTarget];286 MPI_Status* status=newMPI_Status[nbSource+nbTarget];286 ep_lib::MPI_Request* request=new ep_lib::MPI_Request[nbSource+nbTarget]; 287 ep_lib::MPI_Status* status=new ep_lib::MPI_Status[nbSource+nbTarget]; 287 288 int indexRequest=0; 288 289 289 MPI_Barrier(communicator);290 ep_lib::MPI_Barrier(communicator); 290 291 CTimer::get("CMPIRouting::transferToTarget").reset(); 291 292 CTimer::get("CMPIRouting::transferToTarget").resume(); … … 293 294 for(int i=0; i<nbSource; i++) 294 295 { 295 MPI_Irecv(sourceBuffer[i],nbSourceElement[i]*sizeof(T),MPI_CHAR, sourceRank[i], 0, communicator, &request[indexRequest]);296 ep_lib::MPI_Irecv(sourceBuffer[i],nbSourceElement[i]*sizeof(T),EP_CHAR, sourceRank[i], 0, communicator, &request[indexRequest]); 296 297 indexRequest++; 297 298 } … … 299 300 for(int i=0;i<nbTarget; i++) 300 301 { 301 MPI_Isend(targetBuffer[i],nbTargetElement[i]*sizeof(T), MPI_CHAR, targetRank[i], 0, communicator, &request[indexRequest]);302 indexRequest++; 303 } 304 305 MPI_Waitall(indexRequest,request,status);302 ep_lib::MPI_Isend(targetBuffer[i],nbTargetElement[i]*sizeof(T), EP_CHAR, targetRank[i], 0, communicator, &request[indexRequest]); 303 indexRequest++; 304 } 305 306 ep_lib::MPI_Waitall(indexRequest,request,status); 306 307 307 308 CTimer::get("CMPIRouting::transferToTarget").suspend(); 308 309 CTimer::get("CMPIRouting::transferToTarget").print(); 309 MPI_Barrier(communicator);310 ep_lib::MPI_Barrier(communicator); 310 311 311 312 // unpack the data … … 347 348 } 348 349 349 MPI_Request *request = newMPI_Request[nbSource + nbTarget];350 MPI_Status *status = newMPI_Status[nbSource + nbTarget];350 ep_lib::MPI_Request *request = new ep_lib::MPI_Request[nbSource + nbTarget]; 351 ep_lib::MPI_Status *status = new ep_lib::MPI_Status[nbSource + nbTarget]; 351 352 int indexRequest = 0; 352 353 353 MPI_Barrier(communicator);354 ep_lib::MPI_Barrier(communicator); 354 355 CTimer::get("CMPIRouting::transferToTarget(messageSize)").reset(); 355 356 CTimer::get("CMPIRouting::transferToTarget(messageSize)").resume(); … … 357 358 for(int i=0; i<nbSource; i++) 358 359 { 359 MPI_Irecv(&sourceMessageSize[i],1,MPI_INT, sourceRank[i], 0, communicator, &request[indexRequest]);360 ep_lib::MPI_Irecv(&sourceMessageSize[i],1,EP_INT, sourceRank[i], 0, communicator, &request[indexRequest]); 360 361 indexRequest++; 361 362 } … … 363 364 for(int i=0; i<nbTarget; i++) 364 365 { 365 MPI_Isend(&targetMessageSize[i],1, MPI_INT, targetRank[i], 0, communicator, &request[indexRequest]);366 indexRequest++; 367 } 368 369 MPI_Waitall(indexRequest,request,status);370 371 MPI_Barrier(communicator);366 ep_lib::MPI_Isend(&targetMessageSize[i],1, EP_INT, targetRank[i], 0, communicator, &request[indexRequest]); 367 indexRequest++; 368 } 369 370 ep_lib::MPI_Waitall(indexRequest,request,status); 371 372 ep_lib::MPI_Barrier(communicator); 372 373 CTimer::get("CMPIRouting::transferToTarget(messageSize)").suspend(); 373 374 CTimer::get("CMPIRouting::transferToTarget(messageSize)").print(); … … 402 403 for(int i=0; i<nbSource; i++) 403 404 { 404 MPI_Irecv(sourceBuffer[i],sourceMessageSize[i],MPI_CHAR, sourceRank[i], 0, communicator, &request[indexRequest]);405 ep_lib::MPI_Irecv(sourceBuffer[i],sourceMessageSize[i],EP_CHAR, sourceRank[i], 0, communicator, &request[indexRequest]); 405 406 indexRequest++; 406 407 } … … 408 409 for(int i=0;i<nbTarget; i++) 409 410 { 410 MPI_Isend(targetBuffer[i],targetMessageSize[i], MPI_CHAR, targetRank[i], 0, communicator, &request[indexRequest]);411 ep_lib::MPI_Isend(targetBuffer[i],targetMessageSize[i], EP_CHAR, targetRank[i], 0, communicator, &request[indexRequest]); 411 412 indexRequest++; 412 413 } … … 467 468 } 468 469 469 MPI_Request* request=newMPI_Request[nbSource+nbTarget];470 MPI_Status* status=newMPI_Status[nbSource+nbTarget];470 ep_lib::MPI_Request* request=new ep_lib::MPI_Request[nbSource+nbTarget]; 471 ep_lib::MPI_Status* status=new ep_lib::MPI_Status[nbSource+nbTarget]; 471 472 int indexRequest=0; 472 473 473 474 for(int i=0; i<nbSource; i++) 474 475 { 475 MPI_Isend(sourceBuffer[i],nbSourceElement[i]*sizeof(T),MPI_CHAR, sourceRank[i], 0, communicator, &request[indexRequest]);476 ep_lib::MPI_Isend(sourceBuffer[i],nbSourceElement[i]*sizeof(T),EP_CHAR, sourceRank[i], 0, communicator, &request[indexRequest]); 476 477 indexRequest++; 477 478 } … … 479 480 for(int i=0;i<nbTarget; i++) 480 481 { 481 MPI_Irecv(targetBuffer[i],nbTargetElement[i]*sizeof(T), MPI_CHAR, targetRank[i], 0, communicator, &request[indexRequest]);482 indexRequest++; 483 } 484 485 MPI_Waitall(indexRequest,request,status);482 ep_lib::MPI_Irecv(targetBuffer[i],nbTargetElement[i]*sizeof(T), EP_CHAR, targetRank[i], 0, communicator, &request[indexRequest]); 483 indexRequest++; 484 } 485 486 ep_lib::MPI_Waitall(indexRequest,request,status); 486 487 487 488 // unpack the data … … 523 524 } 524 525 525 MPI_Request *request = newMPI_Request[nbSource + nbTarget];526 MPI_Status *status = newMPI_Status[nbSource + nbTarget];526 ep_lib::MPI_Request *request = new ep_lib::MPI_Request[nbSource + nbTarget]; 527 ep_lib::MPI_Status *status = new ep_lib::MPI_Status[nbSource + nbTarget]; 527 528 int indexRequest = 0; 528 529 for (int i = 0; i < nbSource; i++) 529 530 { 530 MPI_Isend(&sourceMessageSize[i], 1, MPI_INT, sourceRank[i], 0, communicator, &request[indexRequest]);531 indexRequest++; 532 } 533 for (int i = 0; i < nbTarget; i++) 534 { 535 MPI_Irecv(&targetMessageSize[i], 1, MPI_INT, targetRank[i], 0, communicator, &request[indexRequest]);536 indexRequest++; 537 } 538 MPI_Waitall(indexRequest, request, status);531 ep_lib::MPI_Isend(&sourceMessageSize[i], 1, EP_INT, sourceRank[i], 0, communicator, &request[indexRequest]); 532 indexRequest++; 533 } 534 for (int i = 0; i < nbTarget; i++) 535 { 536 ep_lib::MPI_Irecv(&targetMessageSize[i], 1, EP_INT, targetRank[i], 0, communicator, &request[indexRequest]); 537 indexRequest++; 538 } 539 ep_lib::MPI_Waitall(indexRequest, request, status); 539 540 540 541 for (int i = 0; i < nbTarget; i++) … … 564 565 for (int i = 0; i < nbSource; i++) 565 566 { 566 MPI_Isend(sourceBuffer[i], sourceMessageSize[i], MPI_CHAR, sourceRank[i], 0, communicator, &request[indexRequest]);567 indexRequest++; 568 } 569 for (int i = 0; i < nbTarget; i++) 570 { 571 MPI_Irecv(targetBuffer[i], targetMessageSize[i], MPI_CHAR, targetRank[i], 0, communicator, &request[indexRequest]);572 indexRequest++; 573 } 574 MPI_Waitall(indexRequest, request, status);567 ep_lib::MPI_Isend(sourceBuffer[i], sourceMessageSize[i], EP_CHAR, sourceRank[i], 0, communicator, &request[indexRequest]); 568 indexRequest++; 569 } 570 for (int i = 0; i < nbTarget; i++) 571 { 572 ep_lib::MPI_Irecv(targetBuffer[i], targetMessageSize[i], EP_CHAR, targetRank[i], 0, communicator, &request[indexRequest]); 573 indexRequest++; 574 } 575 ep_lib::MPI_Waitall(indexRequest, request, status); 575 576 576 577 // unpack the data … … 612 613 613 614 template void alltoalls_unknown(const std::vector<std::vector<NES> >& send, std::vector<std::vector<NES> >& recv, 614 const std::vector<int>& ranks, MPI_Comm communicator);615 const std::vector<int>& ranks, ep_lib::MPI_Comm communicator); 615 616 616 617 template void alltoalls_known(const std::vector<std::vector<int> >& send, std::vector<std::vector<int> >& recv, 617 const std::vector<int>& ranks, MPI_Comm communicator);618 619 } 618 const std::vector<int>& ranks, ep_lib::MPI_Comm communicator); 619 620 } -
XIOS/dev/branch_openmp/extern/remap/src/parallel_tree.cpp
r1545 r1642 12 12 13 13 #include "parallel_tree.hpp" 14 using namespace ep_lib;15 14 16 15 namespace sphereRemap { 17 18 extern CRemapGrid srcGrid;19 #pragma omp threadprivate(srcGrid)20 21 extern CRemapGrid tgtGrid;22 #pragma omp threadprivate(tgtGrid)23 16 24 17 static const int assignLevel = 2; … … 121 114 } 122 115 123 CParallelTree::CParallelTree(MPI_Comm comm) : communicator(comm), cascade(MIN_NODE_SZ*MIN_NODE_SZ*2, comm) 124 { 125 treeCascade.reserve(cascade.num_levels); 126 for (int lev = 0; lev < cascade.num_levels; lev++) 127 treeCascade.push_back(CSampleTree(cascade.level[lev].group_size, assignLevel)); 116 //CParallelTree::CParallelTree(MPI_Comm comm) : communicator(comm), cascade(MIN_NODE_SZ*MIN_NODE_SZ, comm) 117 CParallelTree::CParallelTree(ep_lib::MPI_Comm comm) : communicator(comm), cascade(MAX_NODE_SZ*MAX_NODE_SZ*2, comm) 118 { 119 treeCascade.reserve(cascade.num_levels); 120 for (int lev = 0; lev < cascade.num_levels; lev++) 121 treeCascade.push_back(CSampleTree(cascade.level[lev].group_size, assignLevel)); 128 122 } 129 123 … … 157 151 158 152 int nrecv; // global number of samples THIS WILL BE THE NUMBER OF LEAFS IN THE SAMPLE TREE 159 MPI_Allreduce(&n, &nrecv, 1, MPI_INT, MPI_SUM, comm.comm); // => size of sample tree does not depend on keepNodes! 160 153 ep_lib::MPI_Allreduce(&n, &nrecv, 1, EP_INT, EP_SUM, comm.comm); // => size of sample tree does not depend on keepNodes! 161 154 double ratio = blocSize / (1.0 * nrecv); 162 155 int nsend = ratio * n + 1; // nsend = n_local_samples / n_global_samples * blocksize + 1 = blocksize/comm.size … … 164 157 165 158 int *counts = new int[comm.size]; 166 MPI_Allgather(&nsend, 1, MPI_INT, counts, 1, MPI_INT, comm.comm);159 ep_lib::MPI_Allgather(&nsend, 1, EP_INT, counts, 1, EP_INT, comm.comm); 167 160 168 161 nrecv = 0; … … 190 183 /* each process needs the sample elements from all processes */ 191 184 double *recvBuffer = new double[nrecv*4]; 192 MPI_Allgatherv(sendBuffer, 4 * nsend, MPI_DOUBLE, recvBuffer, counts, displs, MPI_DOUBLE, comm.comm);185 ep_lib::MPI_Allgatherv(sendBuffer, 4 * nsend, EP_DOUBLE, recvBuffer, counts, displs, EP_DOUBLE, comm.comm); 193 186 delete[] sendBuffer; 194 187 delete[] counts; … … 248 241 << " node size : "<<node.size()<<" bloc size : "<<blocSize<<" total number of leaf : "<<tree.leafs.size()<<endl ; 249 242 /* 250 MPI_Allreduce(&ok, &allok, 1, MPI_INT, MPI_PROD, communicator);243 MPI_Allreduce(&ok, &allok, 1, EP_INT, MPI_PROD, communicator); 251 244 if (!allok) { 252 245 MPI_Finalize(); … … 254 247 } 255 248 */ 256 MPI_Abort(MPI_COMM_WORLD,-1) ;249 ep_lib::MPI_Abort(EP_COMM_WORLD,-1) ; 257 250 } 258 251 /* … … 272 265 { 273 266 CMPIRouting MPIRoute(communicator); 274 MPI_Barrier(communicator);267 ep_lib::MPI_Barrier(communicator); 275 268 CTimer::get("buildLocalTree(initRoute)").resume(); 276 269 MPIRoute.init(route); … … 297 290 298 291 int mpiRank; 299 MPI_Comm_rank(communicator, &mpiRank);292 ep_lib::MPI_Comm_rank(communicator, &mpiRank); 300 293 localTree.leafs.reserve(nbLocalElements); 301 294 for (int i = 0; i < nbLocalElements; i++) … … 323 316 nb1=node.size() ; nb2=node2.size() ; 324 317 nb=nb1+nb2 ; 325 MPI_Allreduce(&nb, &nbTot, 1, MPI_LONG, MPI_SUM, communicator) ;318 ep_lib::MPI_Allreduce(&nb, &nbTot, 1, EP_LONG, EP_SUM, communicator) ; 326 319 int commSize ; 327 MPI_Comm_size(communicator,&commSize) ;320 ep_lib::MPI_Comm_size(communicator,&commSize) ; 328 321 329 322 // make multiple of two … … 508 501 // gather circles on this level of the cascade 509 502 int pg_size; 510 MPI_Comm_size(cascade.level[level].pg_comm, &pg_size);503 ep_lib::MPI_Comm_size(cascade.level[level].pg_comm, &pg_size); 511 504 vector<Coord> allRootCentres(pg_size); 512 505 vector<double> allRootRadia(pg_size); 513 MPI_Allgather(&rootCentre, 3, MPI_DOUBLE, &allRootCentres[0], 3, MPI_DOUBLE, cascade.level[level].pg_comm);514 MPI_Allgather(&rootRadius, 1, MPI_DOUBLE, &allRootRadia[0], 1, MPI_DOUBLE, cascade.level[level].pg_comm);506 ep_lib::MPI_Allgather(&rootCentre, 3, EP_DOUBLE, &allRootCentres[0], 3, EP_DOUBLE, cascade.level[level].pg_comm); 507 ep_lib::MPI_Allgather(&rootRadius, 1, EP_DOUBLE, &allRootRadia[0], 1, EP_DOUBLE, cascade.level[level].pg_comm); 515 508 516 509 // now allRootsRadia and allRootCentres must be inserted into second levels of us and propagated to root -
XIOS/dev/branch_openmp/extern/remap/src/parallel_tree.hpp
r1460 r1642 12 12 { 13 13 public: 14 15 14 CParallelTree(ep_lib::MPI_Comm comm); 15 ~CParallelTree(); 16 16 17 17 void build(vector<Node>& node, vector<Node>& node2); 18 18 19 20 19 void routeNodes(vector<int>& route, vector<Node>& nodes, int level = 0); 20 void routeIntersections(vector<vector<int> >& route, vector<Node>& nodes, int level = 0); 21 21 22 23 22 int nbLocalElements; 23 Elt* localElements; 24 24 25 25 CTree localTree; 26 26 27 27 private: 28 29 30 31 28 void updateCirclesForRouting(Coord rootCentre, double rootRadius, int level = 0); 29 void buildSampleTreeCascade(vector<Node>& sampleNodes, int level = 0); 30 void buildLocalTree(const vector<Node>& node, const vector<int>& route); 31 void buildRouteTree(); 32 32 33 34 35 33 //CSampleTree sampleTree; 34 vector<CSampleTree> treeCascade; // first for sample tree, then for routing tree 35 CMPICascade cascade; 36 36 ep_lib::MPI_Comm communicator ; 37 37 -
XIOS/dev/branch_openmp/extern/remap/src/polyg.cpp
r1328 r1642 20 20 void orient(int N, Coord *vertex, Coord *edge, double *d, const Coord &g) 21 21 { 22 23 24 25 26 27 28 29 30 31 32 33 34 35 22 Coord ga = vertex[0] - g; 23 Coord gb = vertex[1] - g; 24 Coord vertical = crossprod(ga, gb); 25 if (N > 2 && scalarprod(g, vertical) < 0) // (GAxGB).G 26 { 27 for (int i = 0; i < N/2; i++) 28 swap(vertex[i], vertex[N-1-i]); 29 30 for (int i = 0; i < (N-1)/2; i++) 31 { 32 swap(edge[N-2-i], edge[i]); 33 swap(d[i], d[N-2-i]); 34 } 35 } 36 36 } 37 37 … … 39 39 void normals(Coord *x, int n, Coord *a) 40 40 { 41 42 41 for (int i = 0; i<n; i++) 42 a[i] = crossprod(x[(i+1)%n], x[i]); 43 43 } 44 44 45 45 Coord barycentre(const Coord *x, int n) 46 46 { 47 48 49 50 51 52 53 54 55 56 57 47 if (n == 0) return ORIGIN; 48 Coord bc = ORIGIN; 49 for (int i = 0; i < n; i++) 50 bc = bc + x[i]; 51 /* both distances can be equal down to roundoff when norm(bc) < mashineepsilon 52 which can occur when weighted with tiny area */ 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)); 56 57 return proj(bc); 58 58 } 59 59 … … 62 62 static Coord tetrah_side_diff_centre(Coord a, Coord b) 63 63 { 64 64 Coord n = crossprod(a,b); 65 65 double sinc2 = n.x*n.x + n.y*n.y + n.z*n.z; 66 66 assert(sinc2 < 1.0 + EPS); 67 67 68 68 // exact: u = asin(sinc)/sinc - 1; asin(sinc) = geodesic length of arc ab 69 69 // approx: 70 70 // double u = sinc2/6. + (3./40.)*sinc2*sinc2; 71 71 72 73 74 75 72 // exact 73 if (sinc2 > 1.0 - EPS) /* if round-off leads to sinc > 1 asin produces NaN */ 74 return n * (M_PI_2 - 1); 75 double sinc = sqrt(sinc2); 76 76 double u = asin(sinc)/sinc - 1; 77 77 … … 82 82 Coord gc_normalintegral(const Coord *x, int n) 83 83 { 84 85 86 87 88 84 Coord m = barycentre(x, n); 85 Coord bc = crossprod(x[n-1]-m, x[0]-m) + tetrah_side_diff_centre(x[n-1], x[0]); 86 for (int i = 1; i < n; i++) 87 bc = bc + crossprod(x[i-1]-m, x[i]-m) + tetrah_side_diff_centre(x[i-1], x[i]); 88 return bc*0.5; 89 89 } 90 90 91 91 Coord exact_barycentre(const Coord *x, int n) 92 92 { 93 94 95 96 97 98 99 93 if (n >= 3) 94 { 95 return proj(gc_normalintegral(x, n)); 96 } 97 else if (n == 0) return ORIGIN; 98 else if (n == 2) return midpoint(x[0], x[1]); 99 else if (n == 1) return x[0]; 100 100 } 101 101 102 102 Coord sc_gc_moon_normalintegral(Coord a, Coord b, Coord pole) 103 103 { 104 double hemisphere = (a.z > 0) ? 1: -1; 105 106 double lat = hemisphere * (M_PI_2 - acos(a.z)); 107 double lon1 = atan2(a.y, a.x); 108 double lon2 = atan2(b.y, b.x); 109 double lon_diff = lon2 - lon1; 110 111 // wraparound at lon=-pi=pi 112 if (lon_diff < -M_PI) lon_diff += 2.0*M_PI; 113 else if (lon_diff > M_PI) lon_diff -= 2.0*M_PI; 114 115 // integral of the normal over the surface bound by great arcs a-pole and b-pole and small arc a-b 116 Coord sc_normalintegral = Coord(0.5*(sin(lon2)-sin(lon1))*(M_PI_2 - lat - 0.5*sin(2.0*lat)), 117 0.5*(cos(lon1)-cos(lon2))*(M_PI_2 - lat - 0.5*sin(2.0*lat)), 118 hemisphere * lon_diff * 0.25 * (cos(2.0*lat) + 1.0)); 119 Coord p = Coord(0,0,hemisphere); // TODO assumes north pole is (0,0,1) 120 Coord t[] = {a, b, p}; 121 if (hemisphere < 0) swap(t[0], t[1]); 122 return (sc_normalintegral - gc_normalintegral(t, 3)) * hemisphere; 123 } 124 125 126 double triarea(Coord& A, Coord& B, Coord& C) 127 { 128 double a = ds(B, C); 129 double b = ds(C, A); 130 double c = ds(A, B); 131 double s = 0.5 * (a + b + c); 132 double t = tan(0.5*s) * tan(0.5*(s - a)) * tan(0.5*(s - b)) * tan(0.5*(s - c)); 133 // assert(t >= 0); 134 if (t<1e-20) return 0. ; 135 return 4 * atan(sqrt(t)); 104 double hemisphere = (a.z > 0) ? 1: -1; 105 106 double lat = hemisphere * (M_PI_2 - acos(a.z)); 107 double lon1 = atan2(a.y, a.x); 108 double lon2 = atan2(b.y, b.x); 109 double lon_diff = lon2 - lon1; 110 111 // wraparound at lon=-pi=pi 112 if (lon_diff < -M_PI) lon_diff += 2.0*M_PI; 113 else if (lon_diff > M_PI) lon_diff -= 2.0*M_PI; 114 115 // integral of the normal over the surface bound by great arcs a-pole and b-pole and small arc a-b 116 Coord sc_normalintegral = Coord(0.5*(sin(lon2)-sin(lon1))*(M_PI_2 - lat - 0.5*sin(2.0*lat)), 117 0.5*(cos(lon1)-cos(lon2))*(M_PI_2 - lat - 0.5*sin(2.0*lat)), 118 hemisphere * lon_diff * 0.25 * (cos(2.0*lat) + 1.0)); 119 Coord p = Coord(0,0,hemisphere); // TODO assumes north pole is (0,0,1) 120 Coord t[] = {a, b, p}; 121 if (hemisphere < 0) swap(t[0], t[1]); 122 return (sc_normalintegral - gc_normalintegral(t, 3)) * hemisphere; 123 } 124 125 126 double triarea(const Coord& A, const Coord& B, const Coord& C) 127 { 128 double a = ds(B, C); 129 double b = ds(C, A); 130 double c = ds(A, B); 131 double tmp ; 132 133 if (a<b) { tmp=a ; a=b ; b=tmp ; } 134 if (c > a) { tmp=a ; a=c ; c=b, b=tmp; } 135 else if (c > b) { tmp=c ; c=b ; b=tmp ; } 136 137 double s = 0.5 * (a + b + c); 138 double t = tan(0.25*(a+(b+c))) * tan(0.25*(c-(a-b))) * tan(0.25*( c + (a-b) )) * tan(0.25*( a + (b - c))); 139 if (t>0) return 4 * atan(sqrt(t)); 140 else 141 { 142 std::cout<<"double triarea(const Coord& A, const Coord& B, const Coord& C) : t < 0 !!! t="<<t<<endl ; 143 return 0 ; 144 } 136 145 } 137 146 … … 142 151 double alun(double b, double d) 143 152 { 144 145 146 147 148 149 150 153 double a = acos(d); 154 assert(b <= 2 * a); 155 double s = a + 0.5 * b; 156 double t = tan(0.5 * s) * tan(0.5 * (s - a)) * tan(0.5 * (s - a)) * tan(0.5 * (s - b)); 157 double r = sqrt(1 - d*d); 158 double p = 2 * asin(sin(0.5*b) / r); 159 return p*(1 - d) - 4*atan(sqrt(t)); 151 160 } 152 161 … … 160 169 double airbar(int N, const Coord *x, const Coord *c, double *d, const Coord& pole, Coord& gg) 161 170 { 162 163 return 0; /* polygons with less than three vertices have zero area */164 165 166 167 168 169 170 171 172 173 174 171 if (N < 3) 172 return 0; /* polygons with less then three vertices have zero area */ 173 Coord t[3]; 174 t[0] = barycentre(x, N); 175 Coord *g = new Coord[N]; 176 double area = 0; 177 Coord gg_exact = gc_normalintegral(x, N); 178 for (int i = 0; i < N; i++) 179 { 180 /* for "spherical circle segment" sum triangular part and "small moon" and => account for small circle */ 181 int ii = (i + 1) % N; 182 t[1] = x[i]; 183 t[2] = x[ii]; 175 184 double sc=scalarprod(crossprod(t[1] - t[0], t[2] - t[0]), t[0]) ; 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 185 assert(sc >= -1e-10); // Error: tri a l'env (wrong orientation) 186 double area_gc = triarea(t[0], t[1], t[2]); 187 double area_sc_gc_moon = 0; 188 if (d[i]) /* handle small circle case */ 189 { 190 Coord m = midpoint(t[1], t[2]); 191 double mext = scalarprod(m, c[i]) - d[i]; 192 char sgl = (mext > 0) ? -1 : 1; 193 area_sc_gc_moon = sgl * alun(arcdist(t[1], t[2]), fabs(scalarprod(t[1], pole))); 194 gg_exact = gg_exact + sc_gc_moon_normalintegral(t[1], t[2], pole); 195 } 196 area += area_gc + area_sc_gc_moon; /* for "spherical circle segment" sum triangular part (at) and "small moon" and => account for small circle */ 197 g[i] = barycentre(t, 3) * (area_gc + area_sc_gc_moon); 198 } 199 gg = barycentre(g, N); 200 gg_exact = proj(gg_exact); 201 delete[] g; 202 gg = gg_exact; 203 return area; 195 204 } 196 205 197 206 double polygonarea(Coord *vertices, int N) 198 207 { 199 200 201 202 203 204 205 206 208 assert(N >= 3); 209 210 /* compute polygon area as sum of triangles */ 211 Coord centre = barycentre(vertices, N); 212 double area = 0; 213 for (int i = 0; i < N; i++) 214 area += triarea(centre, vertices[i], vertices[(i+1)%N]); 215 return area; 207 216 } 208 217 209 218 int packedPolygonSize(const Elt& e) 210 219 { 211 return sizeof(e.id) + sizeof(e.src_id) + sizeof(e.x) + sizeof(e.val)+212 220 return sizeof(e.id) + sizeof(e.src_id) + sizeof(e.x) + sizeof(e.val) + sizeof(e.given_area)+ 221 sizeof(e.n) + e.n*(sizeof(double)+sizeof(Coord)); 213 222 } 214 223 215 224 void packPolygon(const Elt& e, char *buffer, int& pos) 216 225 { 217 *((GloId *) &(buffer[pos])) = e.id; 218 pos += sizeof(e.id); 219 *((GloId *) &(buffer[pos])) = e.src_id; 220 pos += sizeof(e.src_id); 221 222 *((Coord *) &(buffer[pos])) = e.x; 223 pos += sizeof(e.x); 224 225 *((double*) &(buffer[pos])) = e.val; 226 pos += sizeof(e.val); 227 228 *((int *) & (buffer[pos])) = e.n; 229 pos += sizeof(e.n); 230 231 for (int i = 0; i < e.n; i++) 232 { 233 *((double *) & (buffer[pos])) = e.d[i]; 234 pos += sizeof(e.d[i]); 235 236 *((Coord *) &(buffer[pos])) = e.vertex[i]; 237 pos += sizeof(e.vertex[i]); 238 } 226 *((GloId *) &(buffer[pos])) = e.id; 227 pos += sizeof(e.id); 228 *((GloId *) &(buffer[pos])) = e.src_id; 229 pos += sizeof(e.src_id); 230 231 *((Coord *) &(buffer[pos])) = e.x; 232 pos += sizeof(e.x); 233 234 *((double*) &(buffer[pos])) = e.val; 235 pos += sizeof(e.val); 236 237 *((double*) &(buffer[pos])) = e.given_area; 238 pos += sizeof(e.val); 239 240 *((int *) & (buffer[pos])) = e.n; 241 pos += sizeof(e.n); 242 243 for (int i = 0; i < e.n; i++) 244 { 245 *((double *) & (buffer[pos])) = e.d[i]; 246 pos += sizeof(e.d[i]); 247 248 *((Coord *) &(buffer[pos])) = e.vertex[i]; 249 pos += sizeof(e.vertex[i]); 250 } 239 251 240 252 } … … 242 254 void unpackPolygon(Elt& e, const char *buffer, int& pos) 243 255 { 244 e.id = *((GloId *) &(buffer[pos])); 245 pos += sizeof(e.id); 246 e.src_id = *((GloId *) &(buffer[pos])); 247 pos += sizeof(e.src_id); 248 249 e.x = *((Coord *) & (buffer[pos]) ); 250 pos += sizeof(e.x); 251 252 e.val = *((double *) & (buffer[pos])); 253 pos += sizeof(double); 254 255 e.n = *((int *) & (buffer[pos])); 256 pos += sizeof(int); 257 258 for (int i = 0; i < e.n; i++) 259 { 260 e.d[i] = *((double *) & (buffer[pos])); 261 pos += sizeof(double); 262 263 e.vertex[i] = *((Coord *) & (buffer[pos])); 264 pos += sizeof(Coord); 265 } 256 e.id = *((GloId *) &(buffer[pos])); 257 pos += sizeof(e.id); 258 e.src_id = *((GloId *) &(buffer[pos])); 259 pos += sizeof(e.src_id); 260 261 e.x = *((Coord *) & (buffer[pos]) ); 262 pos += sizeof(e.x); 263 264 e.val = *((double *) & (buffer[pos])); 265 pos += sizeof(double); 266 267 e.given_area = *((double *) & (buffer[pos])); 268 pos += sizeof(double); 269 270 e.n = *((int *) & (buffer[pos])); 271 pos += sizeof(int); 272 273 for (int i = 0; i < e.n; i++) 274 { 275 e.d[i] = *((double *) & (buffer[pos])); 276 pos += sizeof(double); 277 278 e.vertex[i] = *((Coord *) & (buffer[pos])); 279 pos += sizeof(Coord); 280 } 266 281 } 267 282 268 283 int packIntersectionSize(const Elt& elt) 269 284 { 270 285 return elt.is.size() * (2*sizeof(int)+ sizeof(GloId) + 5*sizeof(double)); 271 286 } 272 287 … … 274 289 { 275 290 for (list<Polyg *>::const_iterator it = e.is.begin(); it != e.is.end(); ++it) 276 277 278 279 280 291 { 292 *((int *) &(buffer[0])) += 1; 293 294 *((int *) &(buffer[pos])) = e.id.ind; 295 pos += sizeof(int); 281 296 282 297 *((double *) &(buffer[pos])) = e.area; 283 298 pos += sizeof(double); 284 299 285 *((GloId *) &(buffer[pos])) = (*it)->id; 286 pos += sizeof(GloId); 300 301 *((GloId *) &(buffer[pos])) = (*it)->id; 302 pos += sizeof(GloId); 287 303 288 289 290 291 292 293 294 295 304 *((int *) &(buffer[pos])) = (*it)->n; 305 pos += sizeof(int); 306 *((double *) &(buffer[pos])) = (*it)->area; 307 pos += sizeof(double); 308 309 *((Coord *) &(buffer[pos])) = (*it)->x; 310 pos += sizeof(Coord); 311 } 296 312 } 297 313 298 314 void unpackIntersection(Elt* e, const char* buffer) 299 315 { 300 301 316 int ind; 317 int pos = 0; 302 318 303 304 305 306 307 308 309 310 319 int n = *((int *) & (buffer[pos])); 320 pos += sizeof(int); 321 for (int i = 0; i < n; i++) 322 { 323 ind = *((int*) &(buffer[pos])); 324 pos+=sizeof(int); 325 326 Elt& elt= e[ind]; 311 327 312 328 elt.area=*((double *) & (buffer[pos])); 313 pos += sizeof(double); 314 315 Polyg *polygon = new Polyg; 316 317 polygon->id = *((GloId *) & (buffer[pos])); 318 pos += sizeof(GloId); 319 320 polygon->n = *((int *) & (buffer[pos])); 321 pos += sizeof(int); 322 323 polygon->area = *((double *) & (buffer[pos])); 324 pos += sizeof(double); 325 326 polygon->x = *((Coord *) & (buffer[pos])); 327 pos += sizeof(Coord); 328 329 elt.is.push_back(polygon); 330 } 331 } 332 333 } 329 pos += sizeof(double); 330 331 332 Polyg *polygon = new Polyg; 333 334 polygon->id = *((GloId *) & (buffer[pos])); 335 pos += sizeof(GloId); 336 337 polygon->n = *((int *) & (buffer[pos])); 338 pos += sizeof(int); 339 340 polygon->area = *((double *) & (buffer[pos])); 341 pos += sizeof(double); 342 343 polygon->x = *((Coord *) & (buffer[pos])); 344 pos += sizeof(Coord); 345 346 elt.is.push_back(polygon); 347 } 348 } 349 350 } -
XIOS/dev/branch_openmp/extern/remap/src/polyg.hpp
r688 r1642 23 23 void unpackIntersection(Elt *e, const char *buffer); 24 24 int packIntersectionSize(const Elt& e); 25 double triarea( const Coord& A, const Coord& B, const Coord& C) ; 25 26 26 27 } -
XIOS/dev/branch_openmp/extern/remap/src/timerRemap.cpp
r1328 r1642 4 4 #include <map> 5 5 #include <iostream> 6 using namespace ep_lib;7 6 8 7 namespace sphereRemap { … … 10 9 using namespace std; 11 10 12 //map<string,CTimer*> CTimer::allTimer; 13 map<string,CTimer*> *CTimer::allTimer_ptr = 0; 11 map<string,CTimer*> CTimer::allTimer; 14 12 15 13 CTimer::CTimer(const string& name_) : name(name_) … … 58 56 { 59 57 map<string,CTimer*>::iterator it; 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; 58 it=allTimer.find(name); 59 if (it==allTimer.end()) it=allTimer.insert(pair<string,CTimer*>(name,new CTimer(name))).first; 65 60 return *(it->second); 66 61 } -
XIOS/dev/branch_openmp/extern/remap/src/timerRemap.hpp
r1334 r1642 26 26 double getCumulatedTime(void); 27 27 void print(void); 28 //static map<string,CTimer*> allTimer; 29 static map<string,CTimer*> *allTimer_ptr; 30 #pragma omp threadprivate(allTimer_ptr) 28 static map<string,CTimer*> allTimer; 31 29 static double getTime(void); 32 30 static CTimer& get(string name); -
XIOS/dev/branch_openmp/extern/src_ep_dev/ep_allgather.cpp
r1539 r1642 92 92 int local_sendcount = num_ep * count; 93 93 94 ::MPI_Allgather(&local_sendcount, 1, to_mpi_type( MPI_INT), mpi_recvcounts, 1, to_mpi_type(MPI_INT), to_mpi_comm(comm->mpi_comm));94 ::MPI_Allgather(&local_sendcount, 1, to_mpi_type(EP_INT), mpi_recvcounts, 1, to_mpi_type(EP_INT), to_mpi_comm(comm->mpi_comm)); 95 95 96 96 mpi_displs[0] = 0; -
XIOS/dev/branch_openmp/extern/src_ep_dev/ep_allgatherv.cpp
r1539 r1642 56 56 vector<int>local_displs(num_ep, 0); 57 57 58 MPI_Gather_local(&sendcount, 1, MPI_INT, local_recvcounts.data(), 0, comm);58 MPI_Gather_local(&sendcount, 1, EP_INT, local_recvcounts.data(), 0, comm); 59 59 for(int i=1; i<num_ep; i++) local_displs[i] = local_displs[i-1] + local_recvcounts[i-1]; 60 60 … … 75 75 76 76 int local_sendcount = std::accumulate(local_recvcounts.begin(), local_recvcounts.end(), 0); 77 ::MPI_Allgather(&local_sendcount, 1, to_mpi_type( MPI_INT), mpi_recvcounts.data(), 1, to_mpi_type(MPI_INT), to_mpi_comm(comm->mpi_comm));77 ::MPI_Allgather(&local_sendcount, 1, to_mpi_type(EP_INT), mpi_recvcounts.data(), 1, to_mpi_type(EP_INT), to_mpi_comm(comm->mpi_comm)); 78 78 79 79 for(int i=1; i<mpi_size; i++) -
XIOS/dev/branch_openmp/extern/src_ep_dev/ep_allocate.cpp
r1539 r1642 33 33 34 34 int num_ep_max; 35 MPI_Allreduce(&num_ep, &num_ep_max, 1, MPI_INT, MPI_MAX, comm);35 MPI_Allreduce(&num_ep, &num_ep_max, 1, EP_INT, MPI_MAX, comm); 36 36 37 37 assert(num_ep_max > 1); -
XIOS/dev/branch_openmp/extern/src_ep_dev/ep_create.cpp
r1539 r1642 73 73 } 74 74 75 ::MPI_Allgather(&num_ep, 1, to_mpi_type( MPI_INT), &recv_num_ep[0], 1, to_mpi_type(MPI_INT), mpi_base_comm);75 ::MPI_Allgather(&num_ep, 1, to_mpi_type(EP_INT), &recv_num_ep[0], 1, to_mpi_type(EP_INT), mpi_base_comm); 76 76 77 77 … … 124 124 for(int j=0; j<recv_num_ep[i]; j++) 125 125 { 126 out_comm_hdls[0]->ep_rank_map->insert(std::pair< int, std::pair<int,int> >(ind, j, i)); 126 //out_comm_hdls[0]->ep_rank_map->insert(std::pair< int, std::pair<int,int> >(ind, j, i)); 127 (*(out_comm_hdls[0]->ep_rank_map))[ind] = std::make_pair(j,i);//->insert(std::pair< int, std::pair<int,int> >(ind, j, i)); 127 128 ind++; 128 129 } -
XIOS/dev/branch_openmp/extern/src_ep_dev/ep_declaration.cpp
r1520 r1642 4 4 5 5 ::MPI_Comm MPI_COMM_WORLD_STD = MPI_COMM_WORLD; 6 #undef MPI_COMM_WORLD6 //#undef MPI_COMM_WORLD 7 7 8 8 9 9 ::MPI_Comm MPI_COMM_NULL_STD = MPI_COMM_NULL; 10 #undef MPI_COMM_NULL10 //#undef MPI_COMM_NULL 11 11 12 12 13 13 ::MPI_Request MPI_REQUEST_NULL_STD = MPI_REQUEST_NULL; 14 #undef MPI_REQUEST_NULL14 //#undef MPI_REQUEST_NULL 15 15 16 16 ::MPI_Info MPI_INFO_NULL_STD = MPI_INFO_NULL; 17 #undef MPI_INFO_NULL17 //#undef MPI_INFO_NULL 18 18 19 19 ::MPI_Datatype MPI_INT_STD = MPI_INT; … … 26 26 ::MPI_Datatype MPI_UINT64_T_STD = MPI_UINT64_T; 27 27 ::MPI_Datatype MPI_LONG_LONG_INT_STD = MPI_LONG_LONG_INT; 28 28 ::MPI_Datatype MPI_LONG_LONG_STD = MPI_LONG_LONG; 29 /* 29 30 #undef MPI_INT 30 31 #undef MPI_FLOAT … … 36 37 #undef MPI_UINT64_T 37 38 #undef MPI_LONG_LONG_INT 38 39 */ 39 40 40 41 ::MPI_Op MPI_SUM_STD = MPI_SUM; … … 43 44 ::MPI_Op MPI_LOR_STD = MPI_LOR; 44 45 ::MPI_Op MPI_REPLACE_STD = MPI_REPLACE; 45 46 /* 46 47 #undef MPI_SUM 47 48 #undef MPI_MAX … … 49 50 #undef MPI_LOR 50 51 #undef MPI_REPLACE 51 52 */ 52 53 53 54 // _STD defined in ep_type.cpp … … 62 63 extern ::MPI_Datatype MPI_UINT64_T_STD; 63 64 extern ::MPI_Datatype MPI_LONG_LONG_INT_STD; 65 extern ::MPI_Datatype MPI_LONG_LONG_STD; 66 64 67 65 68 … … 77 80 extern ::MPI_Info MPI_INFO_NULL_STD; 78 81 79 ep_lib::MPI_Datatype MPI_INT = &MPI_INT_STD; 80 ep_lib::MPI_Datatype MPI_FLOAT = &MPI_FLOAT_STD; 81 ep_lib::MPI_Datatype MPI_DOUBLE = &MPI_DOUBLE_STD; 82 ep_lib::MPI_Datatype MPI_CHAR = &MPI_CHAR_STD; 83 ep_lib::MPI_Datatype MPI_LONG = &MPI_LONG_STD; 84 ep_lib::MPI_Datatype MPI_UNSIGNED_LONG = &MPI_UNSIGNED_LONG_STD; 85 ep_lib::MPI_Datatype MPI_UNSIGNED_CHAR = &MPI_UNSIGNED_CHAR_STD; 86 ep_lib::MPI_Datatype MPI_UINT64_T = &MPI_UINT64_T_STD; 87 ep_lib::MPI_Datatype MPI_LONG_LONG_INT = &MPI_LONG_LONG_INT_STD; 88 89 90 ep_lib::MPI_Op MPI_SUM = &MPI_SUM_STD; 91 ep_lib::MPI_Op MPI_MAX = &MPI_MAX_STD; 92 ep_lib::MPI_Op MPI_MIN = &MPI_MIN_STD; 93 ep_lib::MPI_Op MPI_LOR = &MPI_LOR_STD; 94 ep_lib::MPI_Op MPI_REPLACE = &MPI_REPLACE_STD; 95 96 ep_lib::ep_comm EP_COMM_WORLD(&MPI_COMM_WORLD_STD); 97 ep_lib::ep_comm EP_COMM_NULL(&MPI_COMM_NULL_STD); 98 99 ep_lib::MPI_Comm MPI_COMM_WORLD = &EP_COMM_WORLD; 100 ep_lib::MPI_Comm MPI_COMM_NULL = &EP_COMM_NULL; 101 102 //ep_lib::ep_status EP_STATUS_IGNORE(&MPI_STATUS_IGNORE_STD); 103 ep_lib::ep_request EP_REQUEST_NULL(&MPI_REQUEST_NULL_STD); 104 ep_lib::ep_info EP_INFO_NULL(&MPI_INFO_NULL_STD); 105 106 //ep_lib::MPI_Status MPI_STATUS_IGNORE = &EP_STATUS_IGNORE; 107 ep_lib::MPI_Request MPI_REQUEST_NULL = &EP_REQUEST_NULL; 108 ep_lib::MPI_Info MPI_INFO_NULL = &EP_INFO_NULL; 82 ep_lib::MPI_Datatype EP_INT = &MPI_INT_STD; 83 ep_lib::MPI_Datatype EP_FLOAT = &MPI_FLOAT_STD; 84 ep_lib::MPI_Datatype EP_DOUBLE = &MPI_DOUBLE_STD; 85 ep_lib::MPI_Datatype EP_CHAR = &MPI_CHAR_STD; 86 ep_lib::MPI_Datatype EP_LONG = &MPI_LONG_STD; 87 ep_lib::MPI_Datatype EP_UNSIGNED_LONG = &MPI_UNSIGNED_LONG_STD; 88 ep_lib::MPI_Datatype EP_UNSIGNED_CHAR = &MPI_UNSIGNED_CHAR_STD; 89 ep_lib::MPI_Datatype EP_UINT64_T = &MPI_UINT64_T_STD; 90 ep_lib::MPI_Datatype EP_LONG_LONG_INT = &MPI_LONG_LONG_INT_STD; 91 ep_lib::MPI_Datatype EP_LONG_LONG = &MPI_LONG_LONG_STD; 109 92 110 93 111 94 95 ep_lib::MPI_Op EP_SUM = &MPI_SUM_STD; 96 ep_lib::MPI_Op EP_MAX = &MPI_MAX_STD; 97 ep_lib::MPI_Op EP_MIN = &MPI_MIN_STD; 98 ep_lib::MPI_Op EP_LOR = &MPI_LOR_STD; 99 ep_lib::MPI_Op EP_REPLACE = &MPI_REPLACE_STD; 100 101 ep_lib::ep_comm EP_COMM_WORLD_t(&MPI_COMM_WORLD_STD); 102 ep_lib::ep_comm EP_COMM_NULL_t(&MPI_COMM_NULL_STD); 103 104 ep_lib::MPI_Comm EP_COMM_WORLD = &EP_COMM_WORLD_t; 105 ep_lib::MPI_Comm EP_COMM_NULL = &EP_COMM_NULL_t; 106 107 //ep_lib::ep_status EP_STATUS_IGNORE(&MPI_STATUS_IGNORE_STD); 108 ep_lib::ep_request EP_REQUEST_NULL_t(&MPI_REQUEST_NULL_STD); 109 ep_lib::ep_info EP_INFO_NULL_t(&MPI_INFO_NULL_STD); 110 111 //ep_lib::MPI_Status MPI_STATUS_IGNORE = &EP_STATUS_IGNORE; 112 ep_lib::MPI_Request EP_REQUEST_NULL = &EP_REQUEST_NULL_t; 113 ep_lib::MPI_Info EP_INFO_NULL = &EP_INFO_NULL_t; 114 115 116 -
XIOS/dev/branch_openmp/extern/src_ep_dev/ep_declaration.hpp
r1520 r1642 1 1 #ifndef EP_DECLARATION_HPP_INCLUDED 2 2 #define EP_DECLARATION_HPP_INCLUDED 3 3 /* 4 4 #undef MPI_INT 5 5 #undef MPI_FLOAT … … 18 18 #undef MPI_REPLACE 19 19 20 #undef MPI_COMM_WORLD20 //#undef MPI_COMM_WORLD 21 21 #undef MPI_COMM_NULL 22 22 … … 24 24 #undef MPI_STATUS_IGNORE 25 25 #undef MPI_INFO_NULL 26 */ 27 extern ep_lib::MPI_Datatype EP_INT; 28 extern ep_lib::MPI_Datatype EP_FLOAT; 29 extern ep_lib::MPI_Datatype EP_DOUBLE; 30 extern ep_lib::MPI_Datatype EP_CHAR; 31 extern ep_lib::MPI_Datatype EP_LONG; 32 extern ep_lib::MPI_Datatype EP_UNSIGNED_LONG; 33 extern ep_lib::MPI_Datatype EP_UNSIGNED_CHAR; 34 extern ep_lib::MPI_Datatype EP_UINT64_T; 35 extern ep_lib::MPI_Datatype EP_LONG_LONG_INT; 36 extern ep_lib::MPI_Datatype EP_LONG_LONG; 26 37 27 extern ep_lib::MPI_Datatype MPI_INT;28 extern ep_lib::MPI_Datatype MPI_FLOAT;29 extern ep_lib::MPI_Datatype MPI_DOUBLE;30 extern ep_lib::MPI_Datatype MPI_CHAR;31 extern ep_lib::MPI_Datatype MPI_LONG;32 extern ep_lib::MPI_Datatype MPI_UNSIGNED_LONG;33 extern ep_lib::MPI_Datatype MPI_UNSIGNED_CHAR;34 extern ep_lib::MPI_Datatype MPI_UINT64_T;35 extern ep_lib::MPI_Datatype MPI_LONG_LONG_INT;36 38 37 extern ep_lib::MPI_Op MPI_SUM;38 extern ep_lib::MPI_Op MPI_MAX;39 extern ep_lib::MPI_Op MPI_MIN;40 extern ep_lib::MPI_Op MPI_LOR;41 extern ep_lib::MPI_Op MPI_REPLACE;39 extern ep_lib::MPI_Op EP_SUM; 40 extern ep_lib::MPI_Op EP_MAX; 41 extern ep_lib::MPI_Op EP_MIN; 42 extern ep_lib::MPI_Op EP_LOR; 43 extern ep_lib::MPI_Op EP_REPLACE; 42 44 43 extern ep_lib::MPI_Comm MPI_COMM_WORLD;44 extern ep_lib::MPI_Comm MPI_COMM_NULL;45 extern ep_lib::MPI_Comm EP_COMM_WORLD; 46 extern ep_lib::MPI_Comm EP_COMM_NULL; 45 47 46 extern ep_lib::MPI_Status MPI_STATUS_IGNORE; 47 extern ep_lib::MPI_Request MPI_REQUEST_NULL; 48 extern ep_lib::MPI_Info MPI_INFO_NULL; 49 48 extern ep_lib::MPI_Status EP_STATUS_IGNORE; 49 extern ep_lib::MPI_Request EP_REQUEST_NULL; 50 extern ep_lib::MPI_Info EP_INFO_NULL; 50 51 51 52 #endif // EP_DECLARATION_HPP_INCLUDED -
XIOS/dev/branch_openmp/extern/src_ep_dev/ep_exscan.cpp
r1540 r1642 78 78 MPI_Barrier_local(comm); 79 79 80 if(op == MPI_SUM)81 { 82 if(datatype == MPI_INT )80 if(op == EP_SUM) 81 { 82 if(datatype == EP_INT ) 83 83 { 84 84 assert(datasize == sizeof(int)); … … 87 87 } 88 88 89 else if(datatype == MPI_FLOAT )89 else if(datatype == EP_FLOAT ) 90 90 { 91 91 assert(datasize == sizeof(float)); … … 95 95 96 96 97 else if(datatype == MPI_DOUBLE )97 else if(datatype == EP_DOUBLE ) 98 98 { 99 99 assert(datasize == sizeof(double)); … … 102 102 } 103 103 104 else if(datatype == MPI_CHAR )104 else if(datatype == EP_CHAR ) 105 105 { 106 106 assert(datasize == sizeof(char)); … … 109 109 } 110 110 111 else if(datatype == MPI_LONG )111 else if(datatype == EP_LONG ) 112 112 { 113 113 assert(datasize == sizeof(long)); … … 116 116 } 117 117 118 else if(datatype == MPI_UNSIGNED_LONG )118 else if(datatype == EP_UNSIGNED_LONG ) 119 119 { 120 120 assert(datasize == sizeof(unsigned long)); … … 123 123 } 124 124 125 else if(datatype == MPI_LONG_LONG_INT )125 else if(datatype == EP_LONG_LONG_INT ) 126 126 { 127 127 assert(datasize == sizeof(long long int)); … … 139 139 } 140 140 141 else if(op == MPI_MAX)142 { 143 if(datatype == MPI_INT )141 else if(op == EP_MAX) 142 { 143 if(datatype == EP_INT ) 144 144 { 145 145 assert(datasize == sizeof(int)); … … 148 148 } 149 149 150 else if(datatype == MPI_FLOAT )150 else if(datatype == EP_FLOAT ) 151 151 { 152 152 assert(datasize == sizeof(float)); … … 155 155 } 156 156 157 else if(datatype == MPI_DOUBLE )157 else if(datatype == EP_DOUBLE ) 158 158 { 159 159 assert(datasize == sizeof(double)); … … 162 162 } 163 163 164 else if(datatype == MPI_CHAR )164 else if(datatype == EP_CHAR ) 165 165 { 166 166 assert(datasize == sizeof(char)); … … 169 169 } 170 170 171 else if(datatype == MPI_LONG )171 else if(datatype == EP_LONG ) 172 172 { 173 173 assert(datasize == sizeof(long)); … … 176 176 } 177 177 178 else if(datatype == MPI_UNSIGNED_LONG )178 else if(datatype == EP_UNSIGNED_LONG ) 179 179 { 180 180 assert(datasize == sizeof(unsigned long)); … … 183 183 } 184 184 185 else if(datatype == MPI_LONG_LONG_INT )185 else if(datatype == EP_LONG_LONG_INT ) 186 186 { 187 187 assert(datasize == sizeof(long long int)); … … 197 197 } 198 198 199 else if(op == MPI_MIN)200 { 201 if(datatype == MPI_INT )199 else if(op == EP_MIN) 200 { 201 if(datatype == EP_INT ) 202 202 { 203 203 assert(datasize == sizeof(int)); … … 206 206 } 207 207 208 else if(datatype == MPI_FLOAT )208 else if(datatype == EP_FLOAT ) 209 209 { 210 210 assert(datasize == sizeof(float)); … … 213 213 } 214 214 215 else if(datatype == MPI_DOUBLE )215 else if(datatype == EP_DOUBLE ) 216 216 { 217 217 assert(datasize == sizeof(double)); … … 220 220 } 221 221 222 else if(datatype == MPI_CHAR )222 else if(datatype == EP_CHAR ) 223 223 { 224 224 assert(datasize == sizeof(char)); … … 227 227 } 228 228 229 else if(datatype == MPI_LONG )229 else if(datatype == EP_LONG ) 230 230 { 231 231 assert(datasize == sizeof(long)); … … 234 234 } 235 235 236 else if(datatype == MPI_UNSIGNED_LONG )236 else if(datatype == EP_UNSIGNED_LONG ) 237 237 { 238 238 assert(datasize == sizeof(unsigned long)); … … 241 241 } 242 242 243 else if(datatype == MPI_LONG_LONG_INT )243 else if(datatype == EP_LONG_LONG_INT ) 244 244 { 245 245 assert(datasize == sizeof(long long int)); -
XIOS/dev/branch_openmp/extern/src_ep_dev/ep_finalize.cpp
r1539 r1642 9 9 int MPI_Finalize() 10 10 { 11 #pragma omp master11 //pragma omp master 12 12 { 13 13 printf("calling EP Finalize\n"); -
XIOS/dev/branch_openmp/extern/src_ep_dev/ep_fortran.cpp
r1520 r1642 8 8 9 9 10 namespace ep_lib11 {10 //namespace ep_lib 11 //{ 12 12 13 void* EP_Comm_c2f( MPI_Comm comm)13 void* EP_Comm_c2f(ep_lib::MPI_Comm comm) 14 14 { 15 15 Debug("MPI_Comm_c2f"); 16 void* fint = new ::MPI_Fint;16 void* fint = new MPI_Fint; 17 17 #ifdef _intelmpi 18 *static_cast< ::MPI_Fint*>(fint) = (::MPI_Fint)(to_mpi_comm(comm->mpi_comm)); 18 //*static_cast< ::MPI_Fint*>(fint) = (::MPI_Fint)(to_mpi_comm(comm->mpi_comm)); 19 *static_cast< MPI_Fint*>(fint) = MPI_Comm_c2f(to_mpi_comm(comm->mpi_comm)); 19 20 #elif _openmpi 20 *static_cast< ::MPI_Fint*>(fint) = MPI_Comm_c2f(to_mpi_comm(comm->mpi_comm));21 *static_cast< MPI_Fint*>(fint) = MPI_Comm_c2f(to_mpi_comm(comm->mpi_comm)); 21 22 #endif 22 23 23 std::map<std::pair< ::MPI_Fint, int>,MPI_Comm > ::iterator it;24 std::map<std::pair< MPI_Fint, int>, ep_lib::MPI_Comm > ::iterator it; 24 25 25 26 #pragma omp critical (fc_comm_map) 26 27 { 27 it = fc_comm_map.find(std::make_pair(*static_cast< ::MPI_Fint*>(fint), omp_get_thread_num()));28 if(it == fc_comm_map.end())28 it = ep_lib::fc_comm_map.find(std::make_pair(*static_cast< MPI_Fint*>(fint), omp_get_thread_num())); 29 if(it == ep_lib::fc_comm_map.end()) 29 30 { 30 fc_comm_map.insert(std::make_pair( std::make_pair( *static_cast< ::MPI_Fint*>(fint), omp_get_thread_num()) , comm));31 printf("EP_Comm_c2f : MAP %p insert: %d, %d, %p\n", & fc_comm_map, *static_cast< ::MPI_Fint*>(fint), omp_get_thread_num(), comm->ep_comm_ptr);31 ep_lib::fc_comm_map.insert(std::make_pair( std::make_pair( *static_cast< MPI_Fint*>(fint), omp_get_thread_num()) , comm)); 32 printf("EP_Comm_c2f : MAP %p insert: %d, %d, %p\n", &(ep_lib::fc_comm_map), *static_cast< MPI_Fint*>(fint), omp_get_thread_num(), comm->ep_comm_ptr); 32 33 } 33 34 } … … 38 39 } 39 40 40 MPI_Comm EP_Comm_f2c(void* comm)41 ep_lib::MPI_Comm EP_Comm_f2c(void* comm) 41 42 { 42 43 Debug("MPI_Comm_f2c"); 43 44 44 45 45 std::map<std::pair< ::MPI_Fint, int>,MPI_Comm > ::iterator it;46 std::map<std::pair< MPI_Fint, int>, ep_lib::MPI_Comm > ::iterator it; 46 47 47 48 #pragma omp critical (fc_comm_map) 48 it = fc_comm_map.find(std::make_pair(*static_cast< ::MPI_Fint*>(comm), omp_get_thread_num()));49 it = ep_lib::fc_comm_map.find(std::make_pair(*static_cast< MPI_Fint*>(comm), omp_get_thread_num())); 49 50 50 if(it != fc_comm_map.end())51 if(it != ep_lib::fc_comm_map.end()) 51 52 { 52 MPI_Comm comm_ptr;53 ep_lib::MPI_Comm comm_ptr; 53 54 comm_ptr = it->second; 54 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);55 printf("EP_Comm_f2c : MAP %p find: %d, %d, %p\n", &(ep_lib::fc_comm_map), it->first.first, it->first.second, comm_ptr->ep_comm_ptr); 55 56 return comm_ptr; 56 57 } … … 58 59 59 60 61 MPI_Comm *base_comm = new MPI_Comm; 60 62 #ifdef _intelmpi 61 ::MPI_Comm *base_comm = new ::MPI_Comm; 62 *base_comm = (::MPI_Comm)(*static_cast< ::MPI_Fint*>(comm)); 63 *base_comm = (MPI_Comm)(*static_cast< MPI_Fint*>(comm)); 63 64 #elif _openmpi 64 ::MPI_Comm *base_comm = ::MPI_Comm_f2c(*static_cast< ::MPI_Fint*>(comm));65 *base_comm = MPI_Comm_f2c(*static_cast< MPI_Fint*>(comm)); 65 66 #endif 66 67 67 if(*base_comm != to_mpi_comm( MPI_COMM_NULL->mpi_comm))68 if(*base_comm != to_mpi_comm(EP_COMM_NULL->mpi_comm)) 68 69 { 69 70 if(omp_get_thread_num() == 0) 70 71 { 71 72 int num_ep = omp_get_num_threads(); 72 MPI_Comm *new_comm;73 MPI_Info info;74 MPI_Comm_create_endpoints(base_comm, num_ep, info, new_comm);75 passage = new_comm;73 ep_lib::MPI_Comm *new_comm; 74 ep_lib::MPI_Info info; 75 ep_lib::MPI_Comm_create_endpoints(base_comm, num_ep, info, new_comm); 76 ep_lib::passage = new_comm; 76 77 } 77 78 #pragma omp barrier 78 79 79 MPI_Comm return_comm =passage[omp_get_thread_num()];80 ep_lib::MPI_Comm return_comm = ep_lib::passage[omp_get_thread_num()]; 80 81 return return_comm; 81 82 82 83 } 83 return MPI_COMM_NULL;84 return EP_COMM_NULL; 84 85 85 86 } 86 87 87 }88 //} 88 89 89 90 -
XIOS/dev/branch_openmp/extern/src_ep_dev/ep_gatherv.cpp
r1539 r1642 87 87 } 88 88 89 MPI_Bcast(recvcounts, ep_size, MPI_INT, root, comm);90 MPI_Bcast(displs, ep_size, MPI_INT, root, comm);89 MPI_Bcast(recvcounts, ep_size, EP_INT, root, comm); 90 MPI_Bcast(displs, ep_size, EP_INT, root, comm); 91 91 92 if(mpi_rank == root_mpi_rank) MPI_Gather_local(&sendcount, 1, MPI_INT, local_recvcounts.data(), root_ep_loc, comm);93 else MPI_Gather_local(&sendcount, 1, MPI_INT, local_recvcounts.data(), 0, comm);92 if(mpi_rank == root_mpi_rank) MPI_Gather_local(&sendcount, 1, EP_INT, local_recvcounts.data(), root_ep_loc, comm); 93 else MPI_Gather_local(&sendcount, 1, EP_INT, local_recvcounts.data(), 0, comm); 94 94 95 95 -
XIOS/dev/branch_openmp/extern/src_ep_dev/ep_intercomm.cpp
r1556 r1642 33 33 } 34 34 35 MPI_Bcast(mpi_rank_of_leader, 2, MPI_INT, local_leader, local_comm);35 MPI_Bcast(mpi_rank_of_leader, 2, EP_INT, local_leader, local_comm); 36 36 37 37 if(mpi_rank_of_leader[0] != mpi_rank_of_leader[1]) … … 86 86 MPI_Comm_rank(peer_comm, &local_leader_rank_in_peer); 87 87 ::MPI_Comm_rank(to_mpi_comm(peer_comm->mpi_comm), &local_leader_rank_in_peer_mpi); 88 ::MPI_Comm_rank(to_mpi_comm( MPI_COMM_WORLD->mpi_comm), &local_leader_rank_in_world);88 ::MPI_Comm_rank(to_mpi_comm(EP_COMM_WORLD->mpi_comm), &local_leader_rank_in_world); 89 89 90 90 send_quadruple[0] = ep_size; … … 99 99 if(remote_leader > local_leader_rank_in_peer) 100 100 { 101 MPI_Isend(send_quadruple, 4, MPI_INT, remote_leader, tag, peer_comm, &request);101 MPI_Isend(send_quadruple, 4, EP_INT, remote_leader, tag, peer_comm, &request); 102 102 MPI_Wait(&request, &status); 103 103 104 MPI_Irecv(recv_quadruple, 4, MPI_INT, remote_leader, tag, peer_comm, &request);104 MPI_Irecv(recv_quadruple, 4, EP_INT, remote_leader, tag, peer_comm, &request); 105 105 MPI_Wait(&request, &status); 106 106 } 107 107 else 108 108 { 109 MPI_Irecv(recv_quadruple, 4, MPI_INT, remote_leader, tag, peer_comm, &request);109 MPI_Irecv(recv_quadruple, 4, EP_INT, remote_leader, tag, peer_comm, &request); 110 110 MPI_Wait(&request, &status); 111 111 112 MPI_Isend(send_quadruple, 4, MPI_INT, remote_leader, tag, peer_comm, &request);112 MPI_Isend(send_quadruple, 4, EP_INT, remote_leader, tag, peer_comm, &request); 113 113 MPI_Wait(&request, &status); 114 114 } … … 123 123 } 124 124 125 MPI_Bcast(send_quadruple, 4, MPI_INT, local_leader, local_comm);126 MPI_Bcast(recv_quadruple, 4, MPI_INT, local_leader, local_comm);125 MPI_Bcast(send_quadruple, 4, EP_INT, local_leader, local_comm); 126 MPI_Bcast(recv_quadruple, 4, EP_INT, local_leader, local_comm); 127 127 128 128 if(!is_local_leader) … … 152 152 153 153 int rank_in_world; 154 ::MPI_Comm_rank(to_mpi_comm( MPI_COMM_WORLD->mpi_comm), &rank_in_world);154 ::MPI_Comm_rank(to_mpi_comm(EP_COMM_WORLD->mpi_comm), &rank_in_world); 155 155 156 156 int *ranks_in_world_local = new int[ep_size]; 157 157 int *ranks_in_world_remote = new int[remote_ep_size]; 158 158 159 MPI_Allgather(&rank_in_world, 1, MPI_INT, ranks_in_world_local, 1, MPI_INT, local_comm);159 MPI_Allgather(&rank_in_world, 1, EP_INT, ranks_in_world_local, 1, EP_INT, local_comm); 160 160 161 161 if(is_local_leader) … … 166 166 if(remote_leader > local_leader_rank_in_peer) 167 167 { 168 MPI_Isend(ranks_in_world_local, ep_size, MPI_INT, remote_leader, tag, peer_comm, &request);168 MPI_Isend(ranks_in_world_local, ep_size, EP_INT, remote_leader, tag, peer_comm, &request); 169 169 MPI_Wait(&request, &status); 170 170 171 MPI_Irecv(ranks_in_world_remote, remote_ep_size, MPI_INT, remote_leader, tag, peer_comm, &request);171 MPI_Irecv(ranks_in_world_remote, remote_ep_size, EP_INT, remote_leader, tag, peer_comm, &request); 172 172 MPI_Wait(&request, &status); 173 173 } 174 174 else 175 175 { 176 MPI_Irecv(ranks_in_world_remote, remote_ep_size, MPI_INT, remote_leader, tag, peer_comm, &request);176 MPI_Irecv(ranks_in_world_remote, remote_ep_size, EP_INT, remote_leader, tag, peer_comm, &request); 177 177 MPI_Wait(&request, &status); 178 178 179 MPI_Isend(ranks_in_world_local, ep_size, MPI_INT, remote_leader, tag, peer_comm, &request);179 MPI_Isend(ranks_in_world_local, ep_size, EP_INT, remote_leader, tag, peer_comm, &request); 180 180 MPI_Wait(&request, &status); 181 181 } … … 185 185 } 186 186 187 MPI_Bcast(ranks_in_world_remote, remote_ep_size, MPI_INT, local_leader, local_comm);187 MPI_Bcast(ranks_in_world_remote, remote_ep_size, EP_INT, local_leader, local_comm); 188 188 189 189 #ifdef _showinfo … … 313 313 int *mpi_rank_list = new int[mpi_size]; 314 314 315 ::MPI_Allgather(&ownership, 1, to_mpi_type( MPI_INT), ownership_list, 1, to_mpi_type(MPI_INT), to_mpi_comm(local_comm->mpi_comm));316 ::MPI_Allgather(&mpi_rank, 1, to_mpi_type( MPI_INT), mpi_rank_list, 1, to_mpi_type(MPI_INT), to_mpi_comm(local_comm->mpi_comm));315 ::MPI_Allgather(&ownership, 1, to_mpi_type(EP_INT), ownership_list, 1, to_mpi_type(EP_INT), to_mpi_comm(local_comm->mpi_comm)); 316 ::MPI_Allgather(&mpi_rank, 1, to_mpi_type(EP_INT), mpi_rank_list, 1, to_mpi_type(EP_INT), to_mpi_comm(local_comm->mpi_comm)); 317 317 318 318 … … 347 347 } 348 348 349 ::MPI_Bcast(&local_leader_rank_in_extracted_comm, 1, to_mpi_type( MPI_INT), local_comm->ep_rank_map->at(local_leader).second, to_mpi_comm(local_comm->mpi_comm));349 ::MPI_Bcast(&local_leader_rank_in_extracted_comm, 1, to_mpi_type(EP_INT), local_comm->ep_rank_map->at(local_leader).second, to_mpi_comm(local_comm->mpi_comm)); 350 350 351 351 ::MPI_Comm *intracomm = new ::MPI_Comm; 352 bool is_real_involved = ownership && extracted_comm != to_mpi_comm( MPI_COMM_NULL->mpi_comm);352 bool is_real_involved = ownership && extracted_comm != to_mpi_comm(EP_COMM_NULL->mpi_comm); 353 353 354 354 if(is_real_involved) … … 489 489 int *remote_rank_info = new int[2*remote_ep_size]; 490 490 491 MPI_Allgather(rank_info, 2, MPI_INT, local_rank_info, 2, MPI_INT, local_comm);491 MPI_Allgather(rank_info, 2, EP_INT, local_rank_info, 2, EP_INT, local_comm); 492 492 493 493 if(is_local_leader) … … 498 498 if(priority) 499 499 { 500 MPI_Isend(local_rank_info, 2*ep_size, MPI_INT, remote_leader, tag, peer_comm, &request);500 MPI_Isend(local_rank_info, 2*ep_size, EP_INT, remote_leader, tag, peer_comm, &request); 501 501 MPI_Wait(&request, &status); 502 502 503 MPI_Irecv(remote_rank_info, 2*remote_ep_size, MPI_INT, remote_leader, tag, peer_comm, &request);503 MPI_Irecv(remote_rank_info, 2*remote_ep_size, EP_INT, remote_leader, tag, peer_comm, &request); 504 504 MPI_Wait(&request, &status); 505 505 } 506 506 else 507 507 { 508 MPI_Irecv(remote_rank_info, 2*remote_ep_size, MPI_INT, remote_leader, tag, peer_comm, &request);509 MPI_Wait(&request, &status); 510 511 MPI_Isend(local_rank_info, 2*ep_size, MPI_INT, remote_leader, tag, peer_comm, &request);512 MPI_Wait(&request, &status); 513 } 514 } 515 516 MPI_Bcast(remote_rank_info, 2*remote_ep_size, MPI_INT, local_leader, local_comm);508 MPI_Irecv(remote_rank_info, 2*remote_ep_size, EP_INT, remote_leader, tag, peer_comm, &request); 509 MPI_Wait(&request, &status); 510 511 MPI_Isend(local_rank_info, 2*ep_size, EP_INT, remote_leader, tag, peer_comm, &request); 512 MPI_Wait(&request, &status); 513 } 514 } 515 516 MPI_Bcast(remote_rank_info, 2*remote_ep_size, EP_INT, local_leader, local_comm); 517 517 518 518 for(int i=0; i<remote_ep_size; i++) … … 537 537 } 538 538 */ 539 539 //MPI_Barrier(*newintercomm); 540 //MPI_Barrier(*newintercomm); 541 //MPI_Barrier(*newintercomm); 540 542 541 543 } -
XIOS/dev/branch_openmp/extern/src_ep_dev/ep_lib.cpp
r1540 r1642 131 131 bool valid_type(MPI_Datatype datatype) 132 132 { 133 if( datatype == MPI_INT134 || datatype == MPI_FLOAT135 || datatype == MPI_DOUBLE136 || datatype == MPI_CHAR137 || datatype == MPI_LONG138 || datatype == MPI_UNSIGNED_LONG139 || datatype == MPI_LONG_LONG)133 if( datatype == EP_INT 134 || datatype == EP_FLOAT 135 || datatype == EP_DOUBLE 136 || datatype == EP_CHAR 137 || datatype == EP_LONG 138 || datatype == EP_UNSIGNED_LONG 139 || datatype == EP_LONG_LONG) 140 140 { 141 141 return true; … … 148 148 bool valid_op(MPI_Op op) 149 149 { 150 if( op == MPI_MAX151 || op == MPI_MIN152 || op == MPI_SUM153 || op == MPI_LOR)150 if( op == EP_MAX 151 || op == EP_MIN 152 || op == EP_SUM 153 || op == EP_LOR) 154 154 { 155 155 return true; -
XIOS/dev/branch_openmp/extern/src_ep_dev/ep_lib.hpp
r1539 r1642 51 51 int MPI_Imrecv(void *buf, int count, MPI_Datatype datatype, MPI_Message *message, MPI_Request *request); 52 52 53 int my_recv (void *buf, int count, MPI_Datatype datatype, int src, int tag, MPI_Comm comm, MPI_Status *status); 54 int my_irecv (void *buf, int count, MPI_Datatype datatype, int src, int tag, MPI_Comm comm, MPI_Request *request); 53 55 54 56 int MPI_Test(MPI_Request *request, int *flag, MPI_Status *status); 55 57 int MPI_Testall(int count, MPI_Request *array_of_requests, int *flag, MPI_Status *array_of_statuses); 56 58 59 int my_test(MPI_Request *request, int *flag, MPI_Status *status); 60 int my_testall(int count, MPI_Request *array_of_requests, int *flag, MPI_Status *array_of_statuses); 61 57 62 int MPI_Wait(MPI_Request *request, MPI_Status *status); 58 63 int MPI_Waitall(int count, MPI_Request *array_of_requests, MPI_Status *array_of_statuses); 59 64 65 int my_wait(MPI_Request *request, MPI_Status *status); 66 int my_waitall(int count, MPI_Request *array_of_requests, MPI_Status *array_of_statuses); 60 67 61 68 int MPI_Iprobe(int source, int tag, MPI_Comm comm, int *flag, MPI_Status *status); -
XIOS/dev/branch_openmp/extern/src_ep_dev/ep_lib_fortran.hpp
r1369 r1642 4 4 #include "ep_type.hpp" 5 5 6 namespace ep_lib7 {6 //namespace ep_lib 7 //{ 8 8 9 void* EP_Comm_c2f( MPI_Comm comm);10 MPI_Comm EP_Comm_f2c(void* comm);11 }9 void* EP_Comm_c2f(ep_lib::MPI_Comm comm); 10 ep_lib::MPI_Comm EP_Comm_f2c(void* comm); 11 //} 12 12 13 13 -
XIOS/dev/branch_openmp/extern/src_ep_dev/ep_merge.cpp
r1539 r1642 36 36 int sum_high; 37 37 38 MPI_Allreduce(&int_high, &sum_high, 1, MPI_INT, MPI_SUM, *newintracomm);38 MPI_Allreduce(&int_high, &sum_high, 1, EP_INT, EP_SUM, *newintracomm); 39 39 40 40 if(sum_high==0 || sum_high==ep_size+remote_ep_size) … … 55 55 56 56 57 MPI_Allgather(my_triple, 3, MPI_INT, my_triple_list, 3, MPI_INT, *newintracomm);57 MPI_Allgather(my_triple, 3, EP_INT, my_triple_list, 3, EP_INT, *newintracomm); 58 58 59 59 -
XIOS/dev/branch_openmp/extern/src_ep_dev/ep_message.cpp
r1539 r1642 147 147 { 148 148 Debug("Message probing for intracomm\n"); 149 149 /* 150 150 #pragma omp critical (_mpi_call) 151 151 { … … 158 158 } 159 159 160 160 */ 161 ::MPI_Improbe(MPI_ANY_SOURCE, MPI_ANY_TAG, to_mpi_comm(comm->mpi_comm), &flag, &message, &status); 162 161 163 if(flag) 162 164 { -
XIOS/dev/branch_openmp/extern/src_ep_dev/ep_mpi.hpp
r1520 r1642 1 1 #ifndef EP_MPI_HPP_INCLUDED 2 2 #define EP_MPI_HPP_INCLUDED 3 3 #ifdef _usingEP 4 4 #include "ep_type.hpp" 5 5 #endif 6 6 MPI_Datatype to_mpi_type(ep_lib::MPI_Datatype type); 7 7 MPI_Op to_mpi_op(ep_lib::MPI_Op op); -
XIOS/dev/branch_openmp/extern/src_ep_dev/ep_probe.cpp
r1539 r1642 68 68 *flag = false; 69 69 70 Message_Check(comm);71 72 #pragma omp flush73 74 70 #pragma omp critical (_query) 75 71 for(Message_list::iterator it = comm->ep_comm_ptr->message_queue->begin(); it!= comm->ep_comm_ptr->message_queue->end(); ++it) … … 98 94 } 99 95 } 96 if(*flag) return 0; 97 98 Message_Check(comm); 99 100 #pragma omp flush 101 102 #pragma omp critical (_query) 103 for(Message_list::iterator it = comm->ep_comm_ptr->message_queue->begin(); it!= comm->ep_comm_ptr->message_queue->end(); ++it) 104 { 105 bool src_matched = src<0? true: (*it)->ep_src == src; 106 bool tag_matched = tag<0? true: (*it)->ep_tag == tag; 107 108 if(src_matched && tag_matched) 109 { 110 Debug("find message\n"); 111 112 status->mpi_status = new ::MPI_Status(*static_cast< ::MPI_Status*>((*it)->mpi_status)); 113 status->ep_src = (*it)->ep_src; 114 status->ep_tag = (*it)->ep_tag; 115 116 if(comm->is_intercomm) 117 { 118 for(INTER_RANK_MAP::iterator iter = comm->inter_rank_map->begin(); iter != comm->inter_rank_map->end(); iter++) 119 { 120 if(iter->second == (*it)->ep_src) status->ep_src=iter->first; 121 } 122 } 123 124 *flag = true; 125 break; 126 } 127 } 128 if(*flag) return 0; 100 129 } 101 130 … … 129 158 *flag = false; 130 159 131 Message_Check(comm);132 133 #pragma omp flush134 135 160 #pragma omp critical (_query) 136 161 if(! comm->ep_comm_ptr->message_queue->empty()) … … 179 204 } 180 205 } 206 207 if(*flag) return 0; 208 209 Message_Check(comm); 210 211 #pragma omp flush 212 213 #pragma omp critical (_query) 214 if(! comm->ep_comm_ptr->message_queue->empty()) 215 { 216 for(Message_list::iterator it = comm->ep_comm_ptr->message_queue->begin(); it!= comm->ep_comm_ptr->message_queue->end(); ++it) 217 { 218 219 bool src_matched = src<0? true: (*it)->ep_src == src; 220 bool tag_matched = tag<0? true: (*it)->ep_tag == tag; 221 222 if(src_matched && tag_matched) 223 { 224 *flag = true; 225 226 status->mpi_status = new ::MPI_Status(*static_cast< ::MPI_Status*>((*it)->mpi_status)); 227 memcheck("new "<< status->mpi_status << " : in ep_lib::MPI_Improbe, status->mpi_status = new ::MPI_Status"); 228 status->ep_src = (*it)->ep_src; 229 status->ep_tag = (*it)->ep_tag; 230 231 (*message)->mpi_message = new ::MPI_Message(*static_cast< ::MPI_Message*>((*it)->mpi_message)); 232 memcheck("new "<< (*message)->mpi_message <<" : in ep_lib::MPI_Improbe, (*message)->mpi_message = new ::MPI_Message"); 233 (*message)->ep_src = (*it)->ep_src; 234 (*message)->ep_tag = (*it)->ep_tag; 235 236 237 #pragma omp critical (_query2) 238 { 239 memcheck("delete "<< (*it)->mpi_message <<" : in ep_lib::Message_Check, delete (*it)->mpi_message"); 240 memcheck("delete "<< (*it)->mpi_status <<" : in ep_lib::Message_Check, delete (*it)->mpi_status"); 241 memcheck("delete "<< (*it) <<" : in ep_lib::Message_Check, delete (*it)"); 242 243 244 delete (*it)->mpi_message; 245 delete (*it)->mpi_status; 246 delete *it; 247 248 249 comm->ep_comm_ptr->message_queue->erase(it); 250 memcheck("message_queue["<<mpi_rank<<","<<ep_rank_loc<<"]->size = "<<comm->ep_comm_ptr->message_queue->size()); 251 #pragma omp flush 252 } 253 254 break; 255 } 256 257 } 258 } 259 260 if(*flag) return 0; 261 181 262 } 182 263 -
XIOS/dev/branch_openmp/extern/src_ep_dev/ep_recv.cpp
r1539 r1642 22 22 int MPI_Recv(void *buf, int count, MPI_Datatype datatype, int src, int tag, MPI_Comm comm, MPI_Status *status) 23 23 { 24 25 24 if(!comm->is_ep) return MPI_Recv_mpi(buf, count, datatype, src, tag, comm, status); 26 25 -
XIOS/dev/branch_openmp/extern/src_ep_dev/ep_reduce.cpp
r1540 r1642 80 80 memcpy(recvbuf, comm->my_buffer->void_buffer[0], datasize * count); 81 81 82 if(op == MPI_MAX)83 { 84 if(datatype == MPI_INT)82 if(op == EP_MAX) 83 { 84 if(datatype == EP_INT) 85 85 { 86 86 assert(datasize == sizeof(int)); … … 96 96 } 97 97 98 else if(datatype == MPI_DOUBLE)98 else if(datatype == EP_DOUBLE) 99 99 { 100 100 assert(datasize == sizeof(double)); … … 103 103 } 104 104 105 else if(datatype == MPI_CHAR)105 else if(datatype == EP_CHAR) 106 106 { 107 107 assert(datasize == sizeof(char)); … … 110 110 } 111 111 112 else if(datatype == MPI_LONG)112 else if(datatype == EP_LONG) 113 113 { 114 114 assert(datasize == sizeof(long)); … … 117 117 } 118 118 119 else if(datatype == MPI_UNSIGNED_LONG)119 else if(datatype == EP_UNSIGNED_LONG) 120 120 { 121 121 assert(datasize == sizeof(unsigned long)); … … 124 124 } 125 125 126 else if(datatype == MPI_LONG_LONG_INT)126 else if(datatype == EP_LONG_LONG_INT) 127 127 { 128 128 assert(datasize == sizeof(long long)); … … 139 139 } 140 140 141 else if(op == MPI_MIN)142 { 143 if(datatype == MPI_INT)141 else if(op == EP_MIN) 142 { 143 if(datatype ==EP_INT) 144 144 { 145 145 assert(datasize == sizeof(int)); … … 155 155 } 156 156 157 else if(datatype == MPI_DOUBLE)157 else if(datatype == EP_DOUBLE) 158 158 { 159 159 assert(datasize == sizeof(double)); … … 162 162 } 163 163 164 else if(datatype == MPI_CHAR)164 else if(datatype == EP_CHAR) 165 165 { 166 166 assert(datasize == sizeof(char)); … … 169 169 } 170 170 171 else if(datatype == MPI_LONG)171 else if(datatype == EP_LONG) 172 172 { 173 173 assert(datasize == sizeof(long)); … … 176 176 } 177 177 178 else if(datatype == MPI_UNSIGNED_LONG)178 else if(datatype == EP_UNSIGNED_LONG) 179 179 { 180 180 assert(datasize == sizeof(unsigned long)); … … 183 183 } 184 184 185 else if(datatype == MPI_LONG_LONG_INT)185 else if(datatype == EP_LONG_LONG_INT) 186 186 { 187 187 assert(datasize == sizeof(long long)); … … 199 199 200 200 201 else if(op == MPI_SUM)202 { 203 if(datatype== MPI_INT)201 else if(op == EP_SUM) 202 { 203 if(datatype==EP_INT) 204 204 { 205 205 assert(datasize == sizeof(int)); … … 215 215 } 216 216 217 else if(datatype == MPI_DOUBLE)217 else if(datatype == EP_DOUBLE) 218 218 { 219 219 assert(datasize == sizeof(double)); … … 222 222 } 223 223 224 else if(datatype == MPI_CHAR)224 else if(datatype == EP_CHAR) 225 225 { 226 226 assert(datasize == sizeof(char)); … … 229 229 } 230 230 231 else if(datatype == MPI_LONG)231 else if(datatype == EP_LONG) 232 232 { 233 233 assert(datasize == sizeof(long)); … … 236 236 } 237 237 238 else if(datatype == MPI_UNSIGNED_LONG)238 else if(datatype ==EP_UNSIGNED_LONG) 239 239 { 240 240 assert(datasize == sizeof(unsigned long)); … … 243 243 } 244 244 245 else if(datatype == MPI_LONG_LONG_INT)245 else if(datatype ==EP_LONG_LONG_INT) 246 246 { 247 247 assert(datasize == sizeof(long long)); … … 258 258 } 259 259 260 else if(op == MPI_LOR)261 { 262 if(datatype != MPI_INT)260 else if(op == EP_LOR) 261 { 262 if(datatype != EP_INT) 263 263 printf("datatype Error, must be MPI_INT\n"); 264 264 else -
XIOS/dev/branch_openmp/extern/src_ep_dev/ep_reduce_scatter.cpp
r1539 r1642 61 61 62 62 int my_recvcount = recvcounts[ep_rank]; 63 MPI_Gather_local(&my_recvcount, 1, MPI_INT, local_recvcounts.data(), 0, comm);64 MPI_Bcast_local(local_recvcounts.data(), num_ep, MPI_INT, 0, comm);63 MPI_Gather_local(&my_recvcount, 1, EP_INT, local_recvcounts.data(), 0, comm); 64 MPI_Bcast_local(local_recvcounts.data(), num_ep, EP_INT, 0, comm); 65 65 66 66 int my_displs = std::accumulate(recvcounts, recvcounts+ep_rank, 0); 67 MPI_Gather_local(&my_displs, 1, MPI_INT, local_displs.data(), 0, comm);68 MPI_Bcast_local(local_displs.data(), num_ep, MPI_INT, 0, comm);67 MPI_Gather_local(&my_displs, 1, EP_INT, local_displs.data(), 0, comm); 68 MPI_Bcast_local(local_displs.data(), num_ep, EP_INT, 0, comm); 69 69 70 70 -
XIOS/dev/branch_openmp/extern/src_ep_dev/ep_scan.cpp
r1540 r1642 60 60 if(ep_rank_loc == 0 && mpi_rank != 0) 61 61 { 62 if(op == MPI_SUM)63 { 64 if(datatype == MPI_INT)62 if(op == EP_SUM) 63 { 64 if(datatype == EP_INT) 65 65 { 66 66 assert(datasize == sizeof(int)); … … 68 68 } 69 69 70 else if(datatype == MPI_FLOAT)70 else if(datatype == EP_FLOAT) 71 71 { 72 72 assert( datasize == sizeof(float)); … … 74 74 } 75 75 76 else if(datatype == MPI_DOUBLE )76 else if(datatype == EP_DOUBLE ) 77 77 { 78 78 assert( datasize == sizeof(double)); … … 80 80 } 81 81 82 else if(datatype == MPI_CHAR)82 else if(datatype == EP_CHAR) 83 83 { 84 84 assert( datasize == sizeof(char)); … … 86 86 } 87 87 88 else if(datatype == MPI_LONG)88 else if(datatype == EP_LONG) 89 89 { 90 90 assert( datasize == sizeof(long)); … … 93 93 94 94 95 else if(datatype == MPI_UNSIGNED_LONG)95 else if(datatype == EP_UNSIGNED_LONG) 96 96 { 97 97 assert(datasize == sizeof(unsigned long)); … … 99 99 } 100 100 101 else if(datatype == MPI_LONG_LONG_INT)101 else if(datatype == EP_LONG_LONG_INT) 102 102 { 103 103 assert(datasize == sizeof(long long int)); … … 112 112 } 113 113 114 else if(op == MPI_MAX)115 { 116 if(datatype == MPI_INT)114 else if(op == EP_MAX) 115 { 116 if(datatype == EP_INT) 117 117 { 118 118 assert( datasize == sizeof(int)); … … 120 120 } 121 121 122 else if(datatype == MPI_FLOAT )122 else if(datatype == EP_FLOAT ) 123 123 { 124 124 assert( datasize == sizeof(float)); … … 126 126 } 127 127 128 else if(datatype == MPI_DOUBLE )128 else if(datatype == EP_DOUBLE ) 129 129 { 130 130 assert( datasize == sizeof(double)); … … 132 132 } 133 133 134 else if(datatype == MPI_CHAR )134 else if(datatype == EP_CHAR ) 135 135 { 136 136 assert(datasize == sizeof(char)); … … 138 138 } 139 139 140 else if(datatype == MPI_LONG)140 else if(datatype == EP_LONG) 141 141 { 142 142 assert( datasize == sizeof(long)); … … 144 144 } 145 145 146 else if(datatype == MPI_UNSIGNED_LONG)146 else if(datatype == EP_UNSIGNED_LONG) 147 147 { 148 148 assert( datasize == sizeof(unsigned long)); … … 150 150 } 151 151 152 else if(datatype == MPI_LONG_LONG_INT)152 else if(datatype == EP_LONG_LONG_INT) 153 153 { 154 154 assert(datasize == sizeof(long long int)); … … 163 163 } 164 164 165 else if(op == MPI_MIN)166 { 167 if(datatype == MPI_INT )165 else if(op == EP_MIN) 166 { 167 if(datatype == EP_INT ) 168 168 { 169 169 assert (datasize == sizeof(int)); … … 171 171 } 172 172 173 else if(datatype == MPI_FLOAT )173 else if(datatype == EP_FLOAT ) 174 174 { 175 175 assert( datasize == sizeof(float)); … … 177 177 } 178 178 179 else if(datatype == MPI_DOUBLE )179 else if(datatype == EP_DOUBLE ) 180 180 { 181 181 assert( datasize == sizeof(double)); … … 183 183 } 184 184 185 else if(datatype == MPI_CHAR )185 else if(datatype == EP_CHAR ) 186 186 { 187 187 assert( datasize == sizeof(char)); … … 189 189 } 190 190 191 else if(datatype == MPI_LONG )191 else if(datatype == EP_LONG ) 192 192 { 193 193 assert( datasize == sizeof(long)); … … 195 195 } 196 196 197 else if(datatype == MPI_UNSIGNED_LONG )197 else if(datatype == EP_UNSIGNED_LONG ) 198 198 { 199 199 assert( datasize == sizeof(unsigned long)); … … 201 201 } 202 202 203 else if(datatype == MPI_LONG_LONG_INT)203 else if(datatype == EP_LONG_LONG_INT) 204 204 { 205 205 assert(datasize == sizeof(long long int)); … … 235 235 236 236 237 if(op == MPI_SUM)238 { 239 if(datatype == MPI_INT )237 if(op == EP_SUM) 238 { 239 if(datatype == EP_INT ) 240 240 { 241 241 assert (datasize == sizeof(int)); … … 244 244 } 245 245 246 else if(datatype == MPI_FLOAT )246 else if(datatype == EP_FLOAT ) 247 247 { 248 248 assert(datasize == sizeof(float)); … … 252 252 253 253 254 else if(datatype == MPI_DOUBLE )254 else if(datatype == EP_DOUBLE ) 255 255 { 256 256 assert(datasize == sizeof(double)); … … 259 259 } 260 260 261 else if(datatype == MPI_CHAR )261 else if(datatype == EP_CHAR ) 262 262 { 263 263 assert(datasize == sizeof(char)); … … 266 266 } 267 267 268 else if(datatype == MPI_LONG )268 else if(datatype == EP_LONG ) 269 269 { 270 270 assert(datasize == sizeof(long)); … … 273 273 } 274 274 275 else if(datatype == MPI_UNSIGNED_LONG )275 else if(datatype == EP_UNSIGNED_LONG ) 276 276 { 277 277 assert(datasize == sizeof(unsigned long)); … … 280 280 } 281 281 282 else if(datatype == MPI_LONG_LONG_INT )282 else if(datatype == EP_LONG_LONG_INT ) 283 283 { 284 284 assert(datasize == sizeof(long long int)); … … 296 296 } 297 297 298 else if(op == MPI_MAX)299 { 300 if(datatype == MPI_INT)298 else if(op == EP_MAX) 299 { 300 if(datatype == EP_INT) 301 301 { 302 302 assert(datasize == sizeof(int)); … … 305 305 } 306 306 307 else if(datatype == MPI_FLOAT )307 else if(datatype == EP_FLOAT ) 308 308 { 309 309 assert(datasize == sizeof(float)); … … 312 312 } 313 313 314 else if(datatype == MPI_DOUBLE )314 else if(datatype == EP_DOUBLE ) 315 315 { 316 316 assert(datasize == sizeof(double)); … … 319 319 } 320 320 321 else if(datatype == MPI_CHAR )321 else if(datatype == EP_CHAR ) 322 322 { 323 323 assert(datasize == sizeof(char)); … … 326 326 } 327 327 328 else if(datatype == MPI_LONG )328 else if(datatype == EP_LONG ) 329 329 { 330 330 assert(datasize == sizeof(long)); … … 333 333 } 334 334 335 else if(datatype == MPI_UNSIGNED_LONG )335 else if(datatype == EP_UNSIGNED_LONG ) 336 336 { 337 337 assert(datasize == sizeof(unsigned long)); … … 340 340 } 341 341 342 else if(datatype == MPI_LONG_LONG_INT )342 else if(datatype == EP_LONG_LONG_INT ) 343 343 { 344 344 assert(datasize == sizeof(long long int)); … … 355 355 } 356 356 357 else if(op == MPI_MIN)358 { 359 if(datatype == MPI_INT )357 else if(op == EP_MIN) 358 { 359 if(datatype == EP_INT ) 360 360 { 361 361 assert(datasize == sizeof(int)); … … 364 364 } 365 365 366 else if(datatype == MPI_FLOAT )366 else if(datatype == EP_FLOAT ) 367 367 { 368 368 assert(datasize == sizeof(float)); … … 371 371 } 372 372 373 else if(datatype == MPI_DOUBLE )373 else if(datatype == EP_DOUBLE ) 374 374 { 375 375 assert(datasize == sizeof(double)); … … 378 378 } 379 379 380 else if(datatype == MPI_CHAR )380 else if(datatype == EP_CHAR ) 381 381 { 382 382 assert(datasize == sizeof(char)); … … 385 385 } 386 386 387 else if(datatype == MPI_LONG )387 else if(datatype == EP_LONG ) 388 388 { 389 389 assert(datasize == sizeof(long)); … … 392 392 } 393 393 394 else if(datatype == MPI_UNSIGNED_LONG )394 else if(datatype == EP_UNSIGNED_LONG ) 395 395 { 396 396 assert(datasize == sizeof(unsigned long)); … … 399 399 } 400 400 401 else if(datatype == MPI_LONG_LONG_INT )401 else if(datatype == EP_LONG_LONG_INT ) 402 402 { 403 403 assert(datasize == sizeof(long long int)); -
XIOS/dev/branch_openmp/extern/src_ep_dev/ep_scatter.cpp
r1539 r1642 73 73 std::vector<int>ranks(ep_size); 74 74 75 if(mpi_rank == root_mpi_rank) MPI_Gather_local(&ep_rank, 1, MPI_INT, local_ranks.data(), root_ep_loc, comm);76 else MPI_Gather_local(&ep_rank, 1, MPI_INT, local_ranks.data(), 0, comm);75 if(mpi_rank == root_mpi_rank) MPI_Gather_local(&ep_rank, 1, EP_INT, local_ranks.data(), root_ep_loc, comm); 76 else MPI_Gather_local(&ep_rank, 1, EP_INT, local_ranks.data(), 0, comm); 77 77 78 78 … … 91 91 displs[i] = displs[i-1] + recvcounts[i-1]; 92 92 93 ::MPI_Gatherv(local_ranks.data(), num_ep, to_mpi_type( MPI_INT), ranks.data(), recvcounts.data(), displs.data(), to_mpi_type(MPI_INT), root_mpi_rank, to_mpi_comm(comm->mpi_comm));93 ::MPI_Gatherv(local_ranks.data(), num_ep, to_mpi_type(EP_INT), ranks.data(), recvcounts.data(), displs.data(), to_mpi_type(EP_INT), root_mpi_rank, to_mpi_comm(comm->mpi_comm)); 94 94 } 95 95 … … 109 109 { 110 110 int local_sendcount = num_ep * count; 111 ::MPI_Gather(&local_sendcount, 1, to_mpi_type( MPI_INT), recvcounts.data(), 1, to_mpi_type(MPI_INT), root_mpi_rank, to_mpi_comm(comm->mpi_comm));111 ::MPI_Gather(&local_sendcount, 1, to_mpi_type(EP_INT), recvcounts.data(), 1, to_mpi_type(EP_INT), root_mpi_rank, to_mpi_comm(comm->mpi_comm)); 112 112 113 113 if(is_root) for(int i=1; i<mpi_size; i++) displs[i] = displs[i-1] + recvcounts[i-1]; -
XIOS/dev/branch_openmp/extern/src_ep_dev/ep_scatterv.cpp
r1539 r1642 76 76 std::vector<int>ranks(ep_size); 77 77 78 if(mpi_rank == root_mpi_rank) MPI_Gather_local(&ep_rank, 1, MPI_INT, local_ranks.data(), root_ep_loc, comm);79 else MPI_Gather_local(&ep_rank, 1, MPI_INT, local_ranks.data(), 0, comm);78 if(mpi_rank == root_mpi_rank) MPI_Gather_local(&ep_rank, 1, EP_INT, local_ranks.data(), root_ep_loc, comm); 79 else MPI_Gather_local(&ep_rank, 1, EP_INT, local_ranks.data(), 0, comm); 80 80 81 81 … … 94 94 my_displs[i] = my_displs[i-1] + recvcounts[i-1]; 95 95 96 ::MPI_Gatherv(local_ranks.data(), num_ep, to_mpi_type( MPI_INT), ranks.data(), recvcounts.data(), my_displs.data(), to_mpi_type(MPI_INT), root_mpi_rank, to_mpi_comm(comm->mpi_comm));96 ::MPI_Gatherv(local_ranks.data(), num_ep, to_mpi_type(EP_INT), ranks.data(), recvcounts.data(), my_displs.data(), to_mpi_type(EP_INT), root_mpi_rank, to_mpi_comm(comm->mpi_comm)); 97 97 } 98 98 … … 112 112 void* local_sendbuf; 113 113 int local_sendcount; 114 if(mpi_rank == root_mpi_rank) MPI_Reduce_local(&recvcount, &local_sendcount, 1, MPI_INT, MPI_SUM, root_ep_loc, comm);115 else MPI_Reduce_local(&recvcount, &local_sendcount, 1, MPI_INT, MPI_SUM, 0, comm);114 if(mpi_rank == root_mpi_rank) MPI_Reduce_local(&recvcount, &local_sendcount, 1, EP_INT, EP_SUM, root_ep_loc, comm); 115 else MPI_Reduce_local(&recvcount, &local_sendcount, 1, EP_INT, EP_SUM, 0, comm); 116 116 117 117 if(is_master) … … 119 119 local_sendbuf = new void*[datasize * local_sendcount]; 120 120 121 ::MPI_Gather(&local_sendcount, 1, to_mpi_type( MPI_INT), recvcounts.data(), 1, to_mpi_type(MPI_INT), root_mpi_rank, to_mpi_comm(comm->mpi_comm));121 ::MPI_Gather(&local_sendcount, 1, to_mpi_type(EP_INT), recvcounts.data(), 1, to_mpi_type(EP_INT), root_mpi_rank, to_mpi_comm(comm->mpi_comm)); 122 122 123 123 if(is_root) for(int i=1; i<mpi_size; i++) my_displs[i] = my_displs[i-1] + recvcounts[i-1]; … … 129 129 std::vector<int>local_displs(num_ep, 0); 130 130 131 MPI_Gather_local(&recvcount, 1, MPI_INT, local_sendcounts.data(), 0, comm);132 MPI_Bcast_local(local_sendcounts.data(), num_ep, MPI_INT, 0, comm);131 MPI_Gather_local(&recvcount, 1, EP_INT, local_sendcounts.data(), 0, comm); 132 MPI_Bcast_local(local_sendcounts.data(), num_ep, EP_INT, 0, comm); 133 133 for(int i=1; i<num_ep; i++) 134 134 local_displs[i] = local_displs[i-1] + local_sendcounts[i-1]; -
XIOS/dev/branch_openmp/extern/src_ep_dev/ep_split.cpp
r1539 r1642 54 54 vector<int> all_color_loc(num_ep); 55 55 56 MPI_Allgather(&color, 1, MPI_INT, all_color.data(), 1, MPI_INT, comm);57 MPI_Allgather_local(&color, 1, MPI_INT, all_color_loc.data(), comm);56 MPI_Allgather(&color, 1, EP_INT, all_color.data(), 1, EP_INT, comm); 57 MPI_Allgather_local(&color, 1, EP_INT, all_color_loc.data(), comm); 58 58 59 59 list<int> color_list(all_color.begin(), all_color.end()); … … 106 106 107 107 108 MPI_Allgather(&key, 1, MPI_INT, all_key.data(),1, MPI_INT, comm);109 MPI_Allgather_local(&key, 1, MPI_INT, all_key_loc.data(), comm);108 MPI_Allgather(&key, 1, EP_INT, all_key.data(),1, EP_INT, comm); 109 MPI_Allgather_local(&key, 1, EP_INT, all_key_loc.data(), comm); 110 110 111 111 vector<int> all_ep_rank(ep_size); … … 115 115 vector<int> colored_ep_rank_loc[num_color]; 116 116 117 MPI_Allgather(&ep_rank, 1, MPI_INT, all_ep_rank.data(),1, MPI_INT, comm);118 MPI_Allgather_local(&ep_rank, 1, MPI_INT, all_ep_rank_loc.data(), comm);117 MPI_Allgather(&ep_rank, 1, EP_INT, all_ep_rank.data(),1, EP_INT, comm); 118 MPI_Allgather_local(&ep_rank, 1, EP_INT, all_ep_rank_loc.data(), comm); 119 119 120 120 for(int i=0; i<num_ep; i++) … … 270 270 if(new_ep_rank_loc == 0) my_triple_vector_recv.resize(3*new_ep_size); 271 271 272 MPI_Gather_local(my_triple, 3, MPI_INT, my_triple_vector.data(), 0, *newcomm);272 MPI_Gather_local(my_triple, 3, EP_INT, my_triple_vector.data(), 0, *newcomm); 273 273 274 274 if(new_ep_rank_loc == 0) … … 277 277 int *displs = new int[new_mpi_size]; 278 278 int new_num_epx3 = new_num_ep * 3; 279 ::MPI_Allgather(&new_num_epx3, 1, to_mpi_type( MPI_INT), recvcounts, 1, to_mpi_type(MPI_INT), to_mpi_comm((*newcomm)->mpi_comm));279 ::MPI_Allgather(&new_num_epx3, 1, to_mpi_type(EP_INT), recvcounts, 1, to_mpi_type(EP_INT), to_mpi_comm((*newcomm)->mpi_comm)); 280 280 displs[0]=0; 281 281 for(int i=1; i<new_mpi_size; i++) 282 282 displs[i] = displs[i-1] + recvcounts[i-1]; 283 283 284 ::MPI_Allgatherv(my_triple_vector.data(), 3*new_num_ep, to_mpi_type( MPI_INT), my_triple_vector_recv.data(), recvcounts, displs, to_mpi_type(MPI_INT), to_mpi_comm((*newcomm)->mpi_comm));284 ::MPI_Allgatherv(my_triple_vector.data(), 3*new_num_ep, to_mpi_type(EP_INT), my_triple_vector_recv.data(), recvcounts, displs, to_mpi_type(EP_INT), to_mpi_comm((*newcomm)->mpi_comm)); 285 285 286 286 for(int i=0; i<new_ep_size; i++) 287 287 { 288 (*newcomm)->ep_comm_ptr->comm_list[0]->ep_rank_map->insert(std::pair< int, std::pair<int,int> >(my_triple_vector_recv[3*i], my_triple_vector_recv[3*i+1], my_triple_vector_recv[3*i+2])); 288 //(*newcomm)->ep_comm_ptr->comm_list[0]->ep_rank_map->insert(std::pair< int, std::pair<int,int> >(my_triple_vector_recv[3*i], my_triple_vector_recv[3*i+1], my_triple_vector_recv[3*i+2])); 289 (*((*newcomm)->ep_comm_ptr->comm_list[0]->ep_rank_map)) [ my_triple_vector_recv[3*i] ] = std::make_pair(my_triple_vector_recv[3*i+1], my_triple_vector_recv[3*i+2]); 289 290 } 290 291 -
XIOS/dev/branch_openmp/extern/src_ep_dev/ep_type.hpp
r1539 r1642 68 68 69 69 MPI_Fint() {} 70 MPI_Fint(void* fint) ;70 MPI_Fint(void* fint) {mpi_fint = fint;}; 71 71 72 72 }; -
XIOS/dev/branch_openmp/extern/src_ep_dev/ep_win.cpp
r1521 r1642 21 21 22 22 int num_ep_max; 23 MPI_Allreduce(&num_ep, &num_ep_max, 1, MPI_INT, MPI_MAX, comm);23 MPI_Allreduce(&num_ep, &num_ep_max, 1, EP_INT, EP_MAX, comm); 24 24 25 25 assert(num_ep_max > 1); … … 88 88 89 89 int num_ep_max; 90 MPI_Allreduce(&num_ep, &num_ep_max, 1, MPI_INT, MPI_MAX, (*win)->comm);90 MPI_Allreduce(&num_ep, &num_ep_max, 1, EP_INT, EP_MAX, (*win)->comm); 91 91 92 92 //printf("rank_loc = %d, thread_num = %d, num_ep_max = %d\n", rank_loc, omp_get_thread_num(), num_ep_max); … … 139 139 140 140 int num_ep_max; 141 MPI_Allreduce(&num_ep, &num_ep_max, 1, MPI_INT, MPI_MAX, win->comm);141 MPI_Allreduce(&num_ep, &num_ep_max, 1, EP_INT, EP_MAX, win->comm); 142 142 143 143 if(num_ep == 1) -
XIOS/dev/branch_openmp/src/array_new.hpp
r1460 r1642 532 532 virtual void fromString(const string& str) { istringstream iss(str); iss >> *this; initialized = true; } 533 533 virtual string toString(void) const { ostringstream oss; oss << *this; return oss.str(); } 534 535 virtual string dump(void) const 536 { 537 ostringstream oss; 538 oss << this->shape()<<" "; 539 if (this->shape().numElements() == 1 && this->shape().dataFirst()[0] == 1) 540 oss << this->dataFirst()[0]; 541 else 542 oss << this->dataFirst()[0] <<" ... "<< this->dataFirst()[this->numElements()-1]; 543 return oss.str(); 544 } 545 534 546 virtual void reset(void) { this->free(); initialized = false; } 535 547 virtual bool isEmpty(void) const { return !initialized; } -
XIOS/dev/branch_openmp/src/attribute.hpp
r1460 r1642 41 41 virtual StdString toString(void) const = 0; 42 42 virtual void fromString(const StdString & str) = 0; 43 virtual StdString dump(void) const = 0; 43 44 virtual bool isEqual(const CAttribute& ) = 0; 44 45 -
XIOS/dev/branch_openmp/src/attribute_array.hpp
r1460 r1642 55 55 virtual bool toBuffer (CBufferOut& buffer) const { return _toBuffer(buffer);} 56 56 virtual bool fromBuffer(CBufferIn& buffer) { return _fromBuffer(buffer); } 57 virtual string dump(void) const { return _dump();} 57 58 58 59 virtual void generateCInterface(ostream& oss,const string& className) ; … … 69 70 CArray<T_numtype, N_rank> inheritedValue ; 70 71 StdString _toString(void) const; 72 StdString _dump(void) const; 71 73 void _fromString(const StdString & str); 72 74 bool _toBuffer (CBufferOut& buffer) const; -
XIOS/dev/branch_openmp/src/attribute_array_impl.hpp
r1460 r1642 129 129 } 130 130 131 template <typename T_numtype, int N_rank> 132 StdString CAttributeArray<T_numtype,N_rank>::_dump(void) const 133 { 134 StdOStringStream oss; 135 if (! isEmpty() && this->hasId() && (this->numElements()!=0)) 136 oss << this->getName() << "=\"" << CArray<T_numtype, N_rank>::dump() << "\""; 137 return (oss.str()); 138 } 139 140 131 141 template <typename T_numtype, int N_rank> 132 142 void CAttributeArray<T_numtype, N_rank>::_fromString(const StdString & str) -
XIOS/dev/branch_openmp/src/attribute_enum.hpp
r1460 r1642 14 14 namespace xios 15 15 { 16 /// ////////////////////// D éclarations ////////////////////// ///16 /// ////////////////////// Declarations ////////////////////// /// 17 17 /*! 18 18 \class CAttributeEnum … … 62 62 virtual StdString toString(void) const { return _toString();} 63 63 virtual void fromString(const StdString & str) { if (str==resetInheritanceStr) { reset(); _canInherite=false ;} else _fromString(str);} 64 virtual StdString dump(void) const { return _toString();} 64 65 65 66 virtual bool toBuffer (CBufferOut& buffer) const { return _toBuffer(buffer);} -
XIOS/dev/branch_openmp/src/attribute_map.cpp
r1460 r1642 28 28 att.second->reset(); 29 29 } 30 } 31 32 ///-------------------------------------------------------------- 33 /*! 34 Dump of all non-empty attributes of an object 35 */ 36 StdString CAttributeMap::dumpXiosAttributes(void) const 37 { 38 int maxNbChar = 250; 39 StdString str; 40 typedef std::pair<StdString, CAttribute*> StdStrAttPair; 41 auto it = SuperClassMap::begin(), end = SuperClassMap::end(); 42 for (; it != end; it++) 43 { 44 const StdStrAttPair& att = *it; 45 if (!att.second->isEmpty()) 46 { 47 if (str.length() < maxNbChar) 48 { 49 str.append(att.second->dump()); 50 str.append(" "); 51 } 52 else if (str.length() == maxNbChar) 53 { 54 str.append("..."); 55 } 56 } 57 } 58 return str; 30 59 } 31 60 -
XIOS/dev/branch_openmp/src/attribute_map.hpp
r1331 r1642 38 38 void duplicateAttributes(const CAttributeMap* const _parent); 39 39 void clearAllAttributes(void); 40 StdString dumpXiosAttributes(void) const; 40 41 41 42 void clearAttribute(const StdString& key); … … 76 77 /// Propriété statique /// 77 78 static CAttributeMap * Current; 78 #pragma omp threadprivate(Current)79 79 80 80 }; // class CAttributeMap -
XIOS/dev/branch_openmp/src/attribute_template.hpp
r1482 r1642 54 54 void checkEmpty(void) const; 55 55 56 56 57 void setInheritedValue(const CAttributeTemplate& attr ); 57 58 void setInheritedValue(const CAttribute& attr ); … … 73 74 // virtual void toBinary (StdOStream & os) const; 74 75 // virtual void fromBinary(StdIStream & is); 76 virtual StdString dump(void) const { return _dump();} 75 77 76 78 virtual bool toBuffer (CBufferOut& buffer) const { return _toBuffer(buffer);} … … 97 99 bool isEqual_(const CAttributeTemplate& attr); 98 100 StdString _toString(void) const; 101 StdString _dump(void) const; 99 102 void _fromString(const StdString & str); 100 103 bool _toBuffer (CBufferOut& buffer) const; -
XIOS/dev/branch_openmp/src/attribute_template_impl.hpp
r1545 r1642 80 80 } 81 81 82 template <class T> 83 T CAttributeTemplate<T>::getValue(void) const 82 83 template <class T> 84 T CAttributeTemplate<T>::getValue(void) const 84 85 { 85 86 return CType<T>::get() ; 86 } 87 88 89 template <class T> 90 void CAttributeTemplate<T>::setValue(const T & value) 87 88 /* 89 if (SuperClass::isEmpty()) 90 { 91 ERROR("T CAttributeTemplate<T>::getValue(void) const", 92 << "[ id = " << this->getId() << "]" 93 << " L'attribut est requis mais n'est pas défini !"); 94 } 95 return (SuperClass::getValue<T>()); 96 */ 97 } 98 99 /* 100 template <class T> 101 T* CAttributeTemplate<T>::getRef(void) 102 { 103 if (SuperClass::isEmpty()) 104 { 105 ERROR("T CAttributeTemplate<T>::getValue(void) const", 106 << "[ id = " << this->getId() << "]" 107 << " L'attribut est requis mais