Changeset 1328


Ignore:
Timestamp:
11/15/17 12:14:34 (7 years ago)
Author:
yushan
Message:

dev_omp

Location:
XIOS/dev/branch_openmp
Files:
179 edited

Legend:

Unmodified
Added
Removed
  • XIOS/dev/branch_openmp/bld.cfg

    r1287 r1328  
    3838#bld::target test_new_features.exe  
    3939#bld::target test_unstruct_complete.exe  
    40 bld::target test_omp.exe  
    41 bld::target test_complete_omp.exe  
    42 bld::target test_remap_omp.exe  
    43 bld::target test_unstruct_omp.exe 
     40#bld::target test_omp.exe  
     41#bld::target test_complete_omp.exe  
     42bld::target test_remap.exe  
     43bld::target test_remap_ref.exe  
     44#bld::target test_unstruct_omp.exe 
    4445#bld::target test_netcdf_omp.exe 
    45 #bld::target test_client.exe  
    46 #bld::target test_complete.exe 
     46bld::target test_client.exe  
     47bld::target test_complete.exe 
    4748#bld::target test_remap.exe 
    4849#bld::target test_unstruct.exe 
  • XIOS/dev/branch_openmp/extern/remap/src/clipper.cpp

    r1205 r1328  
    42744274{ 
    42754275  //The equation of a line in general form (Ax + By + C = 0) 
    4276   //given 2 points (x,y) & (x,y) is ... 
    4277   //(y - y)x + (x - x)y + (y - y)x - (x - x)y = 0 
    4278   //A = (y - y); B = (x - x); C = (y - y)x - (x - x)y 
    4279   //perpendicular distance of point (x,y) = (Ax + By + C)/Sqrt(A + B) 
     4276  //given 2 points (x¹,y¹) & (x²,y²) is ... 
     4277  //(y¹ - y²)x + (x² - x¹)y + (y² - y¹)x¹ - (x² - x¹)y¹ = 0 
     4278  //A = (y¹ - y²); B = (x² - x¹); C = (y² - y¹)x¹ - (x² - x¹)y¹ 
     4279  //perpendicular distance of point (x³,y³) = (Ax³ + By³ + C)/Sqrt(A² + B²) 
    42804280  //see http://en.wikipedia.org/wiki/Perpendicular_distance 
    42814281  double A = double(ln1.Y - ln2.Y); 
  • XIOS/dev/branch_openmp/extern/remap/src/cputime.cpp

    r694 r1328  
    11#include "mpi.hpp" 
     2using namespace ep_lib; 
    23 
    34namespace sphereRemap { 
  • XIOS/dev/branch_openmp/extern/remap/src/elt.hpp

    r1172 r1328  
    5353struct Elt : Polyg 
    5454{ 
    55   Elt() {} 
    56   Elt(const double *bounds_lon, const double *bounds_lat, int max_num_vert) 
    57   { 
    58     int k = 0; 
    59     vertex[k++] = xyz(bounds_lon[0], bounds_lat[0]); 
    60     for (int i = 1; i < max_num_vert; i++) 
    61     { 
    62       vertex[k] = xyz(bounds_lon[i], bounds_lat[i]); 
    63       /* netCDF convention: if first vertex repeats element is finished (at least three vertices == triagle) */ 
    64       if (k >= 3 && squaredist(vertex[k], vertex[0]) < EPS*EPS)  
    65         break; 
    66       /* eliminate zero edges: move to next vertex only if it is different */ 
    67       if (squaredist(vertex[k], vertex[k-1]) > EPS*EPS) 
    68         k++; 
    69       //else cout << "Removed edge " << k << " due to zero length (coinciding endpoints)." << endl ; 
    70     } 
    71     n = k; 
    72     x = barycentre(vertex, n); 
    73   } 
     55        Elt() {} 
     56        Elt(const double *bounds_lon, const double *bounds_lat, int max_num_vert) 
     57        { 
     58                int k = 0; 
     59                vertex[k++] = xyz(bounds_lon[0], bounds_lat[0]); 
     60                for (int i = 1; i < max_num_vert; i++) 
     61                { 
     62                        vertex[k] = xyz(bounds_lon[i], bounds_lat[i]); 
     63                        /* netCDF convention: if first vertex repeats element is finished (at least three vertices == triagle) */ 
     64                        if (k >= 3 && squaredist(vertex[k], vertex[0]) < EPS*EPS)  
     65                                break; 
     66                        /* eliminate zero edges: move to next vertex only if it is different */ 
     67                        if (squaredist(vertex[k], vertex[k-1]) > EPS*EPS) 
     68                                k++; 
     69                        else 
     70                                /* cout << "Removed edge " << k << " due to zero length (coinciding endpoints)." << endl */ ; 
     71                } 
     72                n = k; 
     73                x = barycentre(vertex, n); 
     74        } 
    7475 
    7576        Elt& operator=(const Elt& rhs) 
     
    9596        } 
    9697 
    97   void delete_intersections() 
    98   { 
    99     for (list<Polyg*>::iterator it = this->is.begin(); it != this->is.end(); it++) 
    100     { 
    101       Polyg* poly = *it; 
    102       delete poly; 
    103     } 
    104   } 
     98        void delete_intersections() 
     99        { 
     100                for (list<Polyg*>::iterator it = this->is.begin(); it != this->is.end(); it++) 
     101                { 
     102                        Polyg* poly = *it; 
     103                        delete poly; 
     104                } 
     105        } 
    105106 
    106107  void insert_vertex(int i, const Coord& v) 
  • XIOS/dev/branch_openmp/extern/remap/src/gridRemap.cpp

    r1155 r1328  
    1111 
    1212CRemapGrid srcGrid; 
    13 #pragma omp threadprivate(srcGrid) 
    14  
    1513CRemapGrid tgtGrid; 
    16 #pragma omp threadprivate(tgtGrid) 
    1714 
    1815} 
  • XIOS/dev/branch_openmp/extern/remap/src/gridRemap.hpp

    r1220 r1328  
    1414Coord readPole(std::istream&); 
    1515 
    16  
     16extern CRemapGrid srcGrid; 
     17extern CRemapGrid tgtGrid; 
    1718 
    1819} 
  • XIOS/dev/branch_openmp/extern/remap/src/intersect.cpp

    r1220 r1328  
    1515namespace sphereRemap { 
    1616 
    17 extern CRemapGrid srcGrid; 
    18 #pragma omp threadprivate(srcGrid) 
    19  
    20 extern CRemapGrid tgtGrid; 
    21 #pragma omp threadprivate(tgtGrid) 
    22  
    2317using namespace std; 
    2418 
     
    2721int neighbour_idx(const Elt& a, const Elt& b) 
    2822{ 
    29   for (int i = 0; i < a.n; i++) 
    30   { 
    31     for (int j = 0; j < b.n; j++) 
    32     { 
    33       assert(squaredist(a.vertex[ i       ], b.vertex[ j       ]) > EPS*EPS || 
    34              squaredist(a.vertex[(i+1)%a.n], b.vertex[(j+1)%b.n]) > EPS*EPS); 
    35       if (   squaredist(a.vertex[ i       ], b.vertex[ j           ]) < 1e-13*1e-13 && 
    36              squaredist(a.vertex[(i+1)%a.n], b.vertex[(j+b.n-1)%b.n]) < 1e-13*1e-13) 
    37       { 
    38         return i; 
    39       } 
    40     } 
    41   } 
    42   return NOT_FOUND; 
     23        for (int i = 0; i < a.n; i++) 
     24        { 
     25                for (int j = 0; j < b.n; j++) 
     26                { 
     27                        assert(squaredist(a.vertex[ i       ], b.vertex[ j       ]) > EPS*EPS || 
     28                               squaredist(a.vertex[(i+1)%a.n], b.vertex[(j+1)%b.n]) > EPS*EPS); 
     29                        if (   squaredist(a.vertex[ i       ], b.vertex[ j           ]) < 1e-13*1e-13 && 
     30                               squaredist(a.vertex[(i+1)%a.n], b.vertex[(j+b.n-1)%b.n]) < 1e-13*1e-13) 
     31                        { 
     32                                return i; 
     33                        } 
     34                } 
     35        } 
     36        return NOT_FOUND; 
    4337} 
    4438 
     
    211205bool isNeighbour(Elt& a, const Elt& b) 
    212206{ 
    213   // return neighbour_idx(a, b) != NOT_FOUND; 
     207        // return neighbour_idx(a, b) != NOT_FOUND; 
    214208  return insertNeighbour(a,b,false) ; 
    215209} 
     
    219213void intersect(Elt *a, Elt *b) 
    220214{ 
    221   int na = a->n; /* vertices of a */ 
    222   int nb = b->n; /* vertices of b */ 
    223   Coord *c   = new Coord[na+nb]; 
    224   Coord *c2  = new Coord[na+nb]; 
    225   Coord *xc  = new Coord[na+nb]; 
    226   Coord *xc2 = new Coord[na+nb]; 
    227   Coord gc, gc2; 
    228   double *d = new double[na+nb]; 
    229   double *d2 = new double[na+nb]; 
    230   double are, are2; 
    231   Ipt ipt[NMAX*NMAX]; 
    232   Ipt ipt2[NMAX*NMAX]; 
    233   ptsec(a, b, ipt); 
    234   /* make ipt2 transpose of ipt */ 
    235   for (int ii = 0; ii < na; ii++) 
    236     for (int jj = 0; jj < nb; jj++) 
    237       ipt2[jj*na+ii] = ipt[ii*nb+jj]; 
    238   list<Sgm> iscot; 
    239   recense(a, b, ipt, iscot, 0); 
    240   recense(b, a, ipt2, iscot, 1); 
    241  
    242   int nseg = iscot.size(); 
    243   int nc = 0; 
    244   int nc2 = 0; 
    245   while (iscot.size() && nc < 2) 
    246     nc = assemble(iscot, c, d, xc); 
    247   while (iscot.size() && nc2 < 2) 
    248     nc2 = assemble(iscot, c2, d2, xc2); 
    249 //  assert(nseg == nc + nc2 || nseg == 1); // unused segment 
     215        int na = a->n; /* vertices of a */ 
     216        int nb = b->n; /* vertices of b */ 
     217        Coord *c   = new Coord[na+nb]; 
     218        Coord *c2  = new Coord[na+nb]; 
     219        Coord *xc  = new Coord[na+nb]; 
     220        Coord *xc2 = new Coord[na+nb]; 
     221        Coord gc, gc2; 
     222        double *d = new double[na+nb]; 
     223        double *d2 = new double[na+nb]; 
     224        double are, are2; 
     225        Ipt ipt[NMAX*NMAX]; 
     226        Ipt ipt2[NMAX*NMAX]; 
     227        ptsec(a, b, ipt); 
     228        /* make ipt2 transpose of ipt */ 
     229        for (int ii = 0; ii < na; ii++) 
     230                for (int jj = 0; jj < nb; jj++) 
     231                        ipt2[jj*na+ii] = ipt[ii*nb+jj]; 
     232        list<Sgm> iscot; 
     233        recense(a, b, ipt, iscot, 0); 
     234        recense(b, a, ipt2, iscot, 1); 
     235 
     236        int nseg = iscot.size(); 
     237        int nc = 0; 
     238        int nc2 = 0; 
     239        while (iscot.size() && nc < 2) 
     240                nc = assemble(iscot, c, d, xc); 
     241        while (iscot.size() && nc2 < 2) 
     242                nc2 = assemble(iscot, c2, d2, xc2); 
     243//      assert(nseg == nc + nc2 || nseg == 1); // unused segment 
    250244 
    251245        if (!(nseg == nc + nc2 || nseg == 1)) 
     
    266260 
    267261//    intersect_ym(a,b) ; 
    268   if (nc == 1) nc = 0; 
    269   if (nc2 == 1) nc2 = 0; 
    270   gc = barycentre(xc, nc); 
    271   gc2 = barycentre(xc2, nc2); 
    272   orient(nc, xc, c, d, gc); 
    273  
    274   Coord pole = srcGrid.pole; 
    275   if (pole == ORIGIN) pole = tgtGrid.pole; 
    276   const double MINBASE = 1e-11; 
    277   if (nc == 2) /* nc is the number of vertices of super mesh element */ 
    278   { 
    279     double base = arcdist(xc[0], xc[1]); 
     262        if (nc == 1) nc = 0; 
     263        if (nc2 == 1) nc2 = 0; 
     264        gc = barycentre(xc, nc); 
     265        gc2 = barycentre(xc2, nc2); 
     266        orient(nc, xc, c, d, gc); 
     267 
     268        Coord pole = srcGrid.pole; 
     269        if (pole == ORIGIN) pole = tgtGrid.pole; 
     270        const double MINBASE = 1e-11; 
     271        if (nc == 2) /* nc is the number of vertices of super mesh element */ 
     272        { 
     273                double base = arcdist(xc[0], xc[1]); 
    280274cerr << "DID ARRIVE " << base << xc[0] << xc[1] << endl; 
    281     gc = midpoint(gc, midpointSC(xc[0], xc[1])); 
    282     /* intersection area `are` must be zero here unless we have one great and one small circle */ 
    283     are = alun(base, fabs(scalarprod(xc[0], pole))); 
    284   } 
    285   else 
    286   { 
    287     are = airbar(nc, xc, c, d, pole, gc); 
    288   } 
    289   if (nc2 == 2) 
    290   { 
    291     double base = arcdist(xc2[0], xc2[1]); 
     275                gc = midpoint(gc, midpointSC(xc[0], xc[1])); 
     276                /* intersection area `are` must be zero here unless we have one great and one small circle */ 
     277                are = alun(base, fabs(scalarprod(xc[0], pole))); 
     278        } 
     279        else 
     280        { 
     281                are = airbar(nc, xc, c, d, pole, gc); 
     282        } 
     283        if (nc2 == 2) 
     284        { 
     285                double base = arcdist(xc2[0], xc2[1]); 
    292286cerr << "DID ARRIVE " << base << xc2[0] << xc2[1] << endl; 
    293     assert(base > MINBASE); 
    294     gc2 = midpoint(gc2, midpointSC(xc2[0], xc2[1])); 
    295     are2 = alun(base, fabs(scalarprod(xc2[0], pole))); // 0 
    296   } 
    297   else 
    298   { 
    299     are2 = airbar(nc2, xc2, c2, d2, pole, gc2); 
    300   } 
     287                assert(base > MINBASE); 
     288                gc2 = midpoint(gc2, midpointSC(xc2[0], xc2[1])); 
     289                are2 = alun(base, fabs(scalarprod(xc2[0], pole))); // 0 
     290        } 
     291        else 
     292        { 
     293                are2 = airbar(nc2, xc2, c2, d2, pole, gc2); 
     294        } 
    301295 
    302296//  double ym_area=intersect_ym(a,b) ; 
    303297 
    304298  if (nc > 1) 
    305   { 
    306     /* create one super mesh polygon that src and dest point to */ 
    307     Polyg *is = new Polyg; 
    308     is->x = gc; 
    309     is->area = are; 
    310     is->id = b->id; 
    311     is->src_id = b->src_id; 
    312     is->n = nc; 
    313     (a->is).push_back(is); 
    314     (b->is).push_back(is); 
     299        { 
     300                /* create one super mesh polygon that src and dest point to */ 
     301                Polyg *is = new Polyg; 
     302                is->x = gc; 
     303                is->area = are; 
     304                is->id = b->id; 
     305                is->src_id = b->src_id; 
     306                is->n = nc; 
     307                (a->is).push_back(is); 
     308                (b->is).push_back(is); 
    315309/* 
    316     if (  2*fabs(are-ym_area)/(are+ym_area) > 1.1 && ym_area>1e-8) 
     310          if (  2*fabs(are-ym_area)/(are+ym_area) > 1.1 && ym_area>1e-8) 
    317311    { 
    318312      cout<<"Big area difference : "<<are<<"  "<<ym_area<<endl ; 
     
    320314    } 
    321315*/ 
    322 //    cout<<"intersection : "<<are<<" "<< ym_area<<"  diff : "<<fabs(are-ym_area)<<"  ratio : "<<fabs(are-ym_area)/(0.5*(are+ym_area))<<endl ; 
    323   } 
    324   if (nc2 > 1) 
    325   { 
    326     Polyg *is = new Polyg; 
    327     is->x = gc2; 
    328     is->area = are2; 
    329     is->id = b->id; /* intersection holds id of corresponding source element (see Elt class definition for details about id) */ 
    330     is->src_id = b->src_id; 
    331     is->n = nc2; 
    332     (a->is).push_back(is); 
    333     (b->is).push_back(is); 
     316//              cout<<"intersection : "<<are<<" "<< ym_area<<"  diff : "<<fabs(are-ym_area)<<"  ratio : "<<fabs(are-ym_area)/(0.5*(are+ym_area))<<endl ; 
     317        } 
     318        if (nc2 > 1) 
     319        { 
     320                Polyg *is = new Polyg; 
     321                is->x = gc2; 
     322                is->area = are2; 
     323                is->id = b->id; /* intersection holds id of corresponding source element (see Elt class definition for details about id) */ 
     324                is->src_id = b->src_id; 
     325                is->n = nc2; 
     326                (a->is).push_back(is); 
     327                (b->is).push_back(is); 
    334328/* 
    335     if (  2*fabs(are-ym_area)/(are+ym_area) > 1.1 && ym_area>1e-8 ) 
     329    if (        2*fabs(are-ym_area)/(are+ym_area) > 1.1 && ym_area>1e-8 ) 
    336330    { 
    337331      cout<<"Big area difference : "<<are<<"  "<<ym_area<<endl ; 
     
    339333    } 
    340334*/ 
    341 //    cout<<"intersection : "<<are2<<" "<< ym_area<<"  diff : "<<fabs(are-ym_area)<<"  ratio : "<<fabs(are-ym_area)/(0.5*(are+ym_area))<<endl ; 
    342   } 
     335//              cout<<"intersection : "<<are2<<" "<< ym_area<<"  diff : "<<fabs(are-ym_area)<<"  ratio : "<<fabs(are-ym_area)/(0.5*(are+ym_area))<<endl ; 
     336        } 
    343337/* 
    344338  if (nc<=1 && nc2<=1) 
     
    350344  } 
    351345*/ 
    352   delete [] c; 
    353   delete [] c2; 
    354   delete [] xc; 
    355   delete [] xc2; 
    356   delete [] d; 
    357   delete [] d2; 
    358 } 
    359  
    360 } 
     346        delete [] c; 
     347        delete [] c2; 
     348        delete [] xc; 
     349        delete [] xc2; 
     350        delete [] d; 
     351        delete [] d2; 
     352} 
     353 
     354} 
  • XIOS/dev/branch_openmp/extern/remap/src/intersection_ym.cpp

    r1220 r1328  
    1313 
    1414namespace sphereRemap { 
    15  
    16 extern CRemapGrid srcGrid; 
    17 #pragma omp threadprivate(srcGrid) 
    18  
    19 extern CRemapGrid tgtGrid; 
    20 #pragma omp threadprivate(tgtGrid) 
    21  
    2215 
    2316using namespace std; 
  • XIOS/dev/branch_openmp/extern/remap/src/libmapper.cpp

    r1205 r1328  
    1414#include "mapper.hpp" 
    1515#include "cputime.hpp" // cputime 
    16 #include <stdio.h> 
     16 
     17using namespace ep_lib; 
    1718 
    1819using namespace sphereRemap ; 
     
    2122   and deallocated during the second step (computing the weights) */ 
    2223Mapper *mapper; 
    23 #pragma omp threadprivate(mapper) 
     24 
    2425 
    2526/** xxx_bounds_yyy is of length n_vert_per_cell_xxx*n_cell_xxx 
     
    3334                     int order, int* n_weights) 
    3435{ 
    35   assert(src_bounds_lon); 
    36   assert(src_bounds_lat); 
    37   assert(n_vert_per_cell_src >= 3); 
    38   assert(n_cell_src >= 4); 
    39   assert(dst_bounds_lon); 
    40   assert(dst_bounds_lat); 
    41   assert(n_vert_per_cell_dst >= 3); 
    42   assert(n_cell_dst >= 4); 
    43   assert(1 <= order && order <= 2); 
     36        assert(src_bounds_lon); 
     37        assert(src_bounds_lat); 
     38        assert(n_vert_per_cell_src >= 3); 
     39        assert(n_cell_src >= 4); 
     40        assert(dst_bounds_lon); 
     41        assert(dst_bounds_lat); 
     42        assert(n_vert_per_cell_dst >= 3); 
     43        assert(n_cell_dst >= 4); 
     44        assert(1 <= order && order <= 2); 
    4445 
    4546  mapper = new Mapper(MPI_COMM_WORLD); 
     
    8081        double tic = cputime(); 
    8182        mapper = new Mapper(MPI_COMM_WORLD); 
    82         mapper->setVerbosity(PROGRESS) ; 
     83  mapper->setVerbosity(PROGRESS) ; 
    8384        mapper->buildSSTree(src_msh, dst_msh); 
    8485        double tac = cputime(); 
     
    149150        char **argv = NULL; 
    150151        MPI_Init(&argc, &argv);*/ 
    151         //MPI_Init(NULL, NULL); 
    152         int provided; 
    153         MPI_Init_thread(NULL, NULL, 3, &provided); 
    154         assert(provided >= 3); 
     152        MPI_Init(NULL, NULL); 
    155153} 
    156154 
  • XIOS/dev/branch_openmp/extern/remap/src/mapper.cpp

    r1220 r1328  
    1212 
    1313#include "mapper.hpp" 
     14using namespace ep_lib; 
    1415 
    1516namespace sphereRemap { 
    16  
    17 extern CRemapGrid srcGrid; 
    18 #pragma omp threadprivate(srcGrid) 
    19  
    20 extern CRemapGrid tgtGrid; 
    21 #pragma omp threadprivate(tgtGrid) 
    22  
    2317 
    2418/* A subdivition of an array into N sub-arays 
     
    7266void Mapper::setTargetMesh(const double* boundsLon, const double* boundsLat, int nVertex, int nbCells, const double* pole, const long int* globalId) 
    7367{ 
    74     tgtGrid.pole = Coord(pole[0], pole[1], pole[2]); 
    75  
    76     int mpiRank, mpiSize; 
    77     MPI_Comm_rank(communicator, &mpiRank); 
    78     MPI_Comm_size(communicator, &mpiSize); 
    79  
    80     targetElements.reserve(nbCells); 
    81     targetMesh.reserve(nbCells); 
    82  
    83     targetGlobalId.resize(nbCells) ; 
    84     if (globalId==NULL) 
    85     { 
    86         long int offset ; 
    87         long int nb=nbCells ; 
    88         MPI_Scan(&nb,&offset,1,MPI_LONG,MPI_SUM,communicator) ; 
    89         offset=offset-nb ; 
    90         for(int i=0;i<nbCells;i++) targetGlobalId[i]=offset+i ; 
    91     } 
    92     else targetGlobalId.assign(globalId,globalId+nbCells); 
    93  
    94     for (int i = 0; i < nbCells; i++) 
    95     { 
    96         int offs = i*nVertex; 
    97         Elt elt(boundsLon + offs, boundsLat + offs, nVertex); 
    98         targetElements.push_back(elt); 
    99         targetMesh.push_back(Node(elt.x, cptRadius(elt), &sourceElements.back())); 
    100         cptEltGeom(targetElements[i], Coord(pole[0], pole[1], pole[2])); 
    101     } 
     68  tgtGrid.pole = Coord(pole[0], pole[1], pole[2]); 
     69 
     70  int mpiRank, mpiSize; 
     71  MPI_Comm_rank(communicator, &mpiRank); 
     72  MPI_Comm_size(communicator, &mpiSize); 
     73 
     74  targetElements.reserve(nbCells); 
     75  targetMesh.reserve(nbCells); 
     76 
     77  targetGlobalId.resize(nbCells) ; 
     78  if (globalId==NULL) 
     79  { 
     80    long int offset ; 
     81    long int nb=nbCells ; 
     82    MPI_Scan(&nb,&offset,1,MPI_LONG,MPI_SUM,communicator) ; 
     83    offset=offset-nb ; 
     84    for(int i=0;i<nbCells;i++) targetGlobalId[i]=offset+i ; 
     85  } 
     86  else targetGlobalId.assign(globalId,globalId+nbCells); 
     87 
     88  for (int i = 0; i < nbCells; i++) 
     89  { 
     90    int offs = i*nVertex; 
     91    Elt elt(boundsLon + offs, boundsLat + offs, nVertex); 
     92    targetElements.push_back(elt); 
     93    targetMesh.push_back(Node(elt.x, cptRadius(elt), &sourceElements.back())); 
     94    cptEltGeom(targetElements[i], Coord(pole[0], pole[1], pole[2])); 
     95  } 
    10296 
    10397 
     
    106100void Mapper::setSourceValue(const double* val) 
    107101{ 
    108     int size=sourceElements.size() ; 
    109     for(int i=0;i<size;++i) sourceElements[i].val=val[i] ; 
     102  int size=sourceElements.size() ; 
     103  for(int i=0;i<size;++i) sourceElements[i].val=val[i] ; 
    110104} 
    111105 
    112106void Mapper::getTargetValue(double* val) 
    113107{ 
    114     int size=targetElements.size() ; 
    115     for(int i=0;i<size;++i) val[i]=targetElements[i].val ; 
     108  int size=targetElements.size() ; 
     109  for(int i=0;i<size;++i) val[i]=targetElements[i].val ; 
    116110} 
    117111 
    118112vector<double> Mapper::computeWeights(int interpOrder, bool renormalize, bool quantity) 
    119113{ 
    120     vector<double> timings; 
    121     int mpiSize, mpiRank; 
    122     MPI_Comm_size(communicator, &mpiSize); 
    123     MPI_Comm_rank(communicator, &mpiRank); 
    124  
    125     this->buildSSTree(sourceMesh, targetMesh); 
    126  
    127     if (mpiRank == 0 && verbose) cout << "Computing intersections ..." << endl; 
    128     double tic = cputime(); 
    129     computeIntersection(&targetElements[0], targetElements.size()); 
    130     timings.push_back(cputime() - tic); 
    131  
    132     tic = cputime(); 
    133     if (interpOrder == 2) { 
    134         if (mpiRank == 0 && verbose) cout << "Computing grads ..." << endl; 
    135         buildMeshTopology(); 
    136         computeGrads(); 
    137     } 
    138     timings.push_back(cputime() - tic); 
    139  
    140     /* Prepare computation of weights */ 
    141     /* compute number of intersections which for the first order case 
    142        corresponds to the number of edges in the remap matrix */ 
    143     int nIntersections = 0; 
    144     for (int j = 0; j < targetElements.size(); j++) 
    145     { 
    146         Elt &elt = targetElements[j]; 
    147         for (list<Polyg*>::iterator it = elt.is.begin(); it != elt.is.end(); it++) 
    148             nIntersections++; 
    149     } 
    150     /* overallocate for NMAX neighbours for each elements */ 
    151     remapMatrix = new double[nIntersections*NMAX]; 
    152     srcAddress = new int[nIntersections*NMAX]; 
    153     srcRank = new int[nIntersections*NMAX]; 
    154     dstAddress = new int[nIntersections*NMAX]; 
    155     sourceWeightId =new long[nIntersections*NMAX]; 
    156     targetWeightId =new long[nIntersections*NMAX]; 
    157  
    158  
    159     if (mpiRank == 0 && verbose) cout << "Remapping..." << endl; 
    160     tic = cputime(); 
    161     nWeights = remap(&targetElements[0], targetElements.size(), interpOrder, renormalize, quantity); 
    162     timings.push_back(cputime() - tic); 
    163  
    164     for (int i = 0; i < targetElements.size(); i++) targetElements[i].delete_intersections(); 
    165  
    166     return timings; 
     114  vector<double> timings; 
     115  int mpiSize, mpiRank; 
     116  MPI_Comm_size(communicator, &mpiSize); 
     117  MPI_Comm_rank(communicator, &mpiRank); 
     118 
     119 
     120  this->buildSSTree(sourceMesh, targetMesh); 
     121 
     122  if (mpiRank == 0 && verbose) cout << "Computing intersections ..." << endl; 
     123  double tic = cputime(); 
     124  computeIntersection(&targetElements[0], targetElements.size()); 
     125  timings.push_back(cputime() - tic); 
     126 
     127        tic = cputime(); 
     128        if (interpOrder == 2) { 
     129                if (mpiRank == 0 && verbose) cout << "Computing grads ..." << endl; 
     130                buildMeshTopology(); 
     131                computeGrads(); 
     132        } 
     133        timings.push_back(cputime() - tic); 
     134 
     135        /* Prepare computation of weights */ 
     136        /* compute number of intersections which for the first order case 
     137           corresponds to the number of edges in the remap matrix */ 
     138        int nIntersections = 0; 
     139        for (int j = 0; j < targetElements.size(); j++) 
     140        { 
     141                Elt &elt = targetElements[j]; 
     142                for (list<Polyg*>::iterator it = elt.is.begin(); it != elt.is.end(); it++) 
     143                        nIntersections++; 
     144        } 
     145        /* overallocate for NMAX neighbours for each elements */ 
     146        remapMatrix = new double[nIntersections*NMAX]; 
     147        srcAddress = new int[nIntersections*NMAX]; 
     148        srcRank = new int[nIntersections*NMAX]; 
     149        dstAddress = new int[nIntersections*NMAX]; 
     150  sourceWeightId =new long[nIntersections*NMAX]; 
     151  targetWeightId =new long[nIntersections*NMAX]; 
     152 
     153 
     154        if (mpiRank == 0 && verbose) cout << "Remapping..." << endl; 
     155        tic = cputime(); 
     156        nWeights = remap(&targetElements[0], targetElements.size(), interpOrder, renormalize, quantity); 
     157        timings.push_back(cputime() - tic); 
     158 
     159  for (int i = 0; i < targetElements.size(); i++) targetElements[i].delete_intersections(); 
     160 
     161        return timings; 
    167162} 
    168163 
    169164/** 
    170   @param elements are cells of the target grid that are distributed over CPUs 
    171   indepentently of the distribution of the SS-tree. 
    172   @param nbElements is the size of the elements array. 
    173   @param order is the order of interpolaton (must be 1 or 2). 
    174   */ 
     165   @param elements are cells of the target grid that are distributed over CPUs 
     166          indepentently of the distribution of the SS-tree. 
     167   @param nbElements is the size of the elements array. 
     168   @param order is the order of interpolaton (must be 1 or 2). 
     169*/ 
    175170int Mapper::remap(Elt *elements, int nbElements, int order, bool renormalize, bool quantity) 
    176171{ 
    177     int mpiSize, mpiRank; 
    178     MPI_Comm_size(communicator, &mpiSize); 
    179     MPI_Comm_rank(communicator, &mpiRank); 
    180  
    181     /* create list of intersections (super mesh elements) for each rank */ 
    182     multimap<int, Polyg *> *elementList = new multimap<int, Polyg *>[mpiSize]; 
    183     for (int j = 0; j < nbElements; j++) 
    184     { 
    185         Elt& e = elements[j]; 
    186         for (list<Polyg *>::iterator it = e.is.begin(); it != e.is.end(); it++) 
    187             elementList[(*it)->id.rank].insert(pair<int, Polyg *>((*it)->id.ind, *it)); 
    188     } 
    189  
    190     int *nbSendElement = new int[mpiSize]; 
    191     int **sendElement = new int*[mpiSize]; /* indices of elements required from other rank */ 
    192     double **recvValue = new double*[mpiSize]; 
    193     double **recvArea = new double*[mpiSize]; 
    194     Coord **recvGrad = new Coord*[mpiSize]; 
    195     GloId **recvNeighIds = new GloId*[mpiSize]; /* ids of the of the source neighbours which also contribute through gradient */ 
    196     for (int rank = 0; rank < mpiSize; rank++) 
    197     { 
    198         /* get size for allocation */ 
    199         int last = -1; /* compares unequal to any index */ 
    200         int index = -1; /* increased to starting index 0 in first iteration */ 
    201         for (multimap<int, Polyg *>::iterator it = elementList[rank].begin(); it != elementList[rank].end(); ++it) 
    202         { 
    203             if (last != it->first) 
    204                 index++; 
    205             (it->second)->id.ind = index; 
    206             last = it->first; 
    207         } 
    208         nbSendElement[rank] = index + 1; 
    209  
    210         /* if size is non-zero allocate and collect indices of elements on other ranks that we intersect */ 
    211         if (nbSendElement[rank] > 0) 
    212         { 
    213             sendElement[rank] = new int[nbSendElement[rank]]; 
    214             recvValue[rank]   = new double[nbSendElement[rank]]; 
    215             recvArea[rank]    = new double[nbSendElement[rank]]; 
    216             if (order == 2) 
    217             { 
    218                 recvNeighIds[rank] = new GloId[nbSendElement[rank]*(NMAX+1)]; 
    219                 recvGrad[rank]    = new Coord[nbSendElement[rank]*(NMAX+1)]; 
    220             } 
    221             else 
    222                 recvNeighIds[rank] = new GloId[nbSendElement[rank]]; 
    223  
    224             last = -1; 
    225             index = -1; 
    226             for (multimap<int, Polyg *>::iterator it = elementList[rank].begin(); it != elementList[rank].end(); ++it) 
    227             { 
    228                 if (last != it->first) 
    229                     index++; 
    230                 sendElement[rank][index] = it->first; 
    231                 last = it->first; 
    232             } 
    233         } 
    234     } 
    235  
    236     /* communicate sizes of source elements to be sent (index lists and later values and gradients) */ 
    237     int *nbRecvElement = new int[mpiSize]; 
    238     MPI_Alltoall(nbSendElement, 1, MPI_INT, nbRecvElement, 1, MPI_INT, communicator); 
    239  
    240     /* communicate indices of source elements on other ranks whoes value and gradient we need (since intersection) */ 
    241     int nbSendRequest = 0; 
    242     int nbRecvRequest = 0; 
    243     int **recvElement = new int*[mpiSize]; 
    244     double **sendValue = new double*[mpiSize]; 
    245     double **sendArea = new double*[mpiSize]; 
    246     Coord **sendGrad = new Coord*[mpiSize]; 
    247     GloId **sendNeighIds = new GloId*[mpiSize]; 
    248     MPI_Request *sendRequest = new MPI_Request[4*mpiSize]; 
    249     MPI_Request *recvRequest = new MPI_Request[4*mpiSize]; 
    250     for (int rank = 0; rank < mpiSize; rank++) 
    251     { 
    252         if (nbSendElement[rank] > 0) 
    253         { 
    254             MPI_Issend(sendElement[rank], nbSendElement[rank], MPI_INT, rank, 0, communicator, &sendRequest[nbSendRequest]); 
    255             nbSendRequest++; 
    256         } 
    257  
    258         if (nbRecvElement[rank] > 0) 
    259         { 
    260             recvElement[rank] = new int[nbRecvElement[rank]]; 
    261             sendValue[rank]   = new double[nbRecvElement[rank]]; 
    262             sendArea[rank]   = new double[nbRecvElement[rank]]; 
    263             if (order == 2) 
    264             { 
    265                 sendNeighIds[rank] = new GloId[nbRecvElement[rank]*(NMAX+1)]; 
    266                 sendGrad[rank]    = new Coord[nbRecvElement[rank]*(NMAX+1)]; 
    267             } 
    268             else 
    269             { 
    270                 sendNeighIds[rank] = new GloId[nbRecvElement[rank]]; 
    271             } 
    272             MPI_Irecv(recvElement[rank], nbRecvElement[rank], MPI_INT, rank, 0, communicator, &recvRequest[nbRecvRequest]); 
    273             nbRecvRequest++; 
    274         } 
    275     } 
    276  
    277     MPI_Status *status = new MPI_Status[4*mpiSize]; 
    278  
    279     MPI_Waitall(nbSendRequest, sendRequest, status); 
    280     MPI_Waitall(nbRecvRequest, recvRequest, status); 
    281  
    282     /* for all indices that have been received from requesting ranks: pack values and gradients, then send */ 
    283     nbSendRequest = 0; 
    284     nbRecvRequest = 0; 
    285     for (int rank = 0; rank < mpiSize; rank++) 
    286     { 
    287         if (nbRecvElement[rank] > 0) 
    288         { 
    289             int jj = 0; // jj == j if no weight writing 
    290             for (int j = 0; j < nbRecvElement[rank]; j++) 
    291             { 
    292                 sendValue[rank][j] = sstree.localElements[recvElement[rank][j]].val; 
    293                 sendArea[rank][j] = sstree.localElements[recvElement[rank][j]].area; 
    294                 if (order == 2) 
    295                 { 
    296                     sendGrad[rank][jj] = sstree.localElements[recvElement[rank][j]].grad; 
    297                     sendNeighIds[rank][jj] = sstree.localElements[recvElement[rank][j]].src_id; 
    298                     jj++; 
    299                     for (int i = 0; i < NMAX; i++) 
    300                     { 
    301                         sendGrad[rank][jj] = sstree.localElements[recvElement[rank][j]].gradNeigh[i]; 
    302                         sendNeighIds[rank][jj] = sstree.localElements[recvElement[rank][j]].neighId[i]; 
    303                         jj++; 
    304                     } 
    305                 } 
    306                 else 
    307                     sendNeighIds[rank][j] = sstree.localElements[recvElement[rank][j]].src_id; 
    308             } 
    309             MPI_Issend(sendValue[rank], nbRecvElement[rank], MPI_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]); 
    310             nbSendRequest++; 
    311             MPI_Issend(sendArea[rank], nbRecvElement[rank], MPI_DOUBLE, rank, 1, communicator, &sendRequest[nbSendRequest]); 
    312             nbSendRequest++; 
    313             if (order == 2) 
    314             { 
    315                 MPI_Issend(sendGrad[rank], 3*nbRecvElement[rank]*(NMAX+1), 
    316                         MPI_DOUBLE, rank, 2, communicator, &sendRequest[nbSendRequest]); 
    317                 nbSendRequest++; 
    318                 MPI_Issend(sendNeighIds[rank], 4*nbRecvElement[rank]*(NMAX+1), MPI_INT, rank, 3, communicator, &sendRequest[nbSendRequest]); 
    319                 //ym  --> attention taille GloId 
    320                 nbSendRequest++; 
    321             } 
    322             else 
    323             { 
    324                 MPI_Issend(sendNeighIds[rank], 4*nbRecvElement[rank], MPI_INT, rank, 4, communicator, &sendRequest[nbSendRequest]); 
    325                 //ym  --> attention taille GloId 
    326                 nbSendRequest++;                 
    327             } 
    328         } 
    329         if (nbSendElement[rank] > 0) 
    330         { 
    331             MPI_Irecv(recvValue[rank],  nbSendElement[rank], MPI_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]); 
    332             nbRecvRequest++; 
    333             MPI_Irecv(recvArea[rank],  nbSendElement[rank], MPI_DOUBLE, rank, 1, communicator, &recvRequest[nbRecvRequest]); 
    334             nbRecvRequest++; 
    335             if (order == 2) 
    336             { 
    337                 MPI_Irecv(recvGrad[rank], 3*nbSendElement[rank]*(NMAX+1), 
    338                         MPI_DOUBLE, rank, 2, communicator, &recvRequest[nbRecvRequest]); 
    339                 nbRecvRequest++; 
    340                 MPI_Irecv(recvNeighIds[rank], 4*nbSendElement[rank]*(NMAX+1), MPI_INT, rank, 3, communicator, &recvRequest[nbRecvRequest]); 
    341                 //ym  --> attention taille GloId 
    342                 nbRecvRequest++; 
    343             } 
    344             else 
    345             { 
    346                 MPI_Irecv(recvNeighIds[rank], 4*nbSendElement[rank], MPI_INT, rank, 4, communicator, &recvRequest[nbRecvRequest]); 
    347                 //ym  --> attention taille GloId 
    348                 nbRecvRequest++; 
    349             } 
    350         } 
    351     } 
    352      
    353     MPI_Waitall(nbRecvRequest, recvRequest, status); 
    354     MPI_Waitall(nbSendRequest, sendRequest, status);  
    355      
    356      
    357  
    358     /* now that all values and gradients are available use them to computed interpolated values on target 
    359        and also to compute weights */ 
    360     int i = 0; 
    361     for (int j = 0; j < nbElements; j++) 
    362     { 
    363         Elt& e = elements[j]; 
    364  
    365         /* since for the 2nd order case source grid elements can contribute to a destination grid element over several "paths" 
    366            (step1: gradient is computed using neighbours on same grid, step2: intersection uses several elements on other grid) 
    367            accumulate them so that there is only one final weight between two elements */ 
    368         map<GloId,double> wgt_map; 
    369  
    370         /* for destination element `e` loop over all intersetions/the corresponding source elements */ 
    371         for (list<Polyg *>::iterator it = e.is.begin(); it != e.is.end(); it++) 
    372         { 
    373             /* it is the intersection element, so it->x and it->area are barycentre and area of intersection element (super mesh) 
    374                but it->id is id of the source element that it intersects */ 
    375             int n1 = (*it)->id.ind; 
    376             int rank = (*it)->id.rank; 
    377             double fk = recvValue[rank][n1]; 
    378             double srcArea = recvArea[rank][n1]; 
    379             double w = (*it)->area; 
    380             if (quantity) w/=srcArea ; 
    381  
    382             /* first order: src value times weight (weight = supermesh area), later divide by target area */ 
    383             int kk = (order == 2) ? n1 * (NMAX + 1) : n1; 
    384             GloId neighID = recvNeighIds[rank][kk]; 
    385             wgt_map[neighID] += w; 
    386  
    387             if (order == 2) 
    388             { 
    389                 for (int k = 0; k < NMAX+1; k++) 
    390                 { 
    391                     int kk = n1 * (NMAX + 1) + k; 
    392                     GloId neighID = recvNeighIds[rank][kk]; 
    393                     if (neighID.ind != -1)  wgt_map[neighID] += w * scalarprod(recvGrad[rank][kk], (*it)->x); 
    394                 } 
    395  
    396             } 
    397         } 
    398  
    399         double renorm=0; 
    400         if (renormalize)  
    401             for (map<GloId,double>::iterator it = wgt_map.begin(); it != wgt_map.end(); it++) renorm+=it->second / e.area; 
    402         else renorm=1. ; 
    403  
    404         for (map<GloId,double>::iterator it = wgt_map.begin(); it != wgt_map.end(); it++) 
    405         { 
    406             if (quantity)  this->remapMatrix[i] = (it->second ) / renorm; 
    407             else this->remapMatrix[i] = (it->second / e.area) / renorm; 
    408             this->srcAddress[i] = it->first.ind; 
    409             this->srcRank[i] = it->first.rank; 
    410             this->dstAddress[i] = j; 
    411             this->sourceWeightId[i]= it->first.globalId ; 
    412             this->targetWeightId[i]= targetGlobalId[j] ; 
    413             i++; 
    414         } 
    415     } 
    416      
    417     /* free all memory allocated in this function */ 
    418     for (int rank = 0; rank < mpiSize; rank++) 
    419     { 
    420         if (nbSendElement[rank] > 0) 
    421         { 
    422             delete[] sendElement[rank]; 
    423             delete[] recvValue[rank]; 
    424             delete[] recvArea[rank]; 
    425             if (order == 2) 
    426             { 
    427                 delete[] recvGrad[rank]; 
    428             } 
    429             delete[] recvNeighIds[rank]; 
    430         } 
    431         if (nbRecvElement[rank] > 0) 
    432         { 
    433             delete[] recvElement[rank]; 
    434             delete[] sendValue[rank]; 
    435             delete[] sendArea[rank]; 
    436             if (order == 2) 
    437                 delete[] sendGrad[rank]; 
    438             delete[] sendNeighIds[rank]; 
    439         } 
    440     } 
    441     delete[] status; 
    442     delete[] sendRequest; 
    443     delete[] recvRequest; 
    444     delete[] elementList; 
    445     delete[] nbSendElement; 
    446     delete[] nbRecvElement; 
    447     delete[] sendElement; 
    448     delete[] recvElement; 
    449     delete[] sendValue; 
    450     delete[] recvValue; 
    451     delete[] sendGrad; 
    452     delete[] recvGrad; 
    453     delete[] sendNeighIds; 
    454     delete[] recvNeighIds; 
    455      
    456     return i; 
     172        int mpiSize, mpiRank; 
     173        MPI_Comm_size(communicator, &mpiSize); 
     174        MPI_Comm_rank(communicator, &mpiRank); 
     175 
     176        /* create list of intersections (super mesh elements) for each rank */ 
     177        multimap<int, Polyg *> *elementList = new multimap<int, Polyg *>[mpiSize]; 
     178        for (int j = 0; j < nbElements; j++) 
     179        { 
     180                Elt& e = elements[j]; 
     181                for (list<Polyg *>::iterator it = e.is.begin(); it != e.is.end(); it++) 
     182                        elementList[(*it)->id.rank].insert(pair<int, Polyg *>((*it)->id.ind, *it)); 
     183        } 
     184 
     185        int *nbSendElement = new int[mpiSize]; 
     186        int **sendElement = new int*[mpiSize]; /* indices of elements required from other rank */ 
     187        double **recvValue = new double*[mpiSize]; 
     188        double **recvArea = new double*[mpiSize]; 
     189        Coord **recvGrad = new Coord*[mpiSize]; 
     190        GloId **recvNeighIds = new GloId*[mpiSize]; /* ids of the of the source neighbours which also contribute through gradient */ 
     191        for (int rank = 0; rank < mpiSize; rank++) 
     192        { 
     193                /* get size for allocation */ 
     194                int last = -1; /* compares unequal to any index */ 
     195                int index = -1; /* increased to starting index 0 in first iteration */ 
     196                for (multimap<int, Polyg *>::iterator it = elementList[rank].begin(); it != elementList[rank].end(); ++it) 
     197                { 
     198                        if (last != it->first) 
     199                                index++; 
     200                        (it->second)->id.ind = index; 
     201                        last = it->first; 
     202                } 
     203                nbSendElement[rank] = index + 1; 
     204 
     205                /* if size is non-zero allocate and collect indices of elements on other ranks that we intersect */ 
     206                if (nbSendElement[rank] > 0) 
     207                { 
     208                        sendElement[rank] = new int[nbSendElement[rank]]; 
     209                        recvValue[rank]   = new double[nbSendElement[rank]]; 
     210                        recvArea[rank]    = new double[nbSendElement[rank]]; 
     211                        if (order == 2) 
     212                        { 
     213                                recvNeighIds[rank] = new GloId[nbSendElement[rank]*(NMAX+1)]; 
     214                                recvGrad[rank]    = new Coord[nbSendElement[rank]*(NMAX+1)]; 
     215                        } 
     216                        else 
     217                                recvNeighIds[rank] = new GloId[nbSendElement[rank]]; 
     218 
     219                        last = -1; 
     220                        index = -1; 
     221                        for (multimap<int, Polyg *>::iterator it = elementList[rank].begin(); it != elementList[rank].end(); ++it) 
     222                        { 
     223                                if (last != it->first) 
     224                                        index++; 
     225                                sendElement[rank][index] = it->first; 
     226                                last = it->first; 
     227                        } 
     228                } 
     229        } 
     230 
     231        /* communicate sizes of source elements to be sent (index lists and later values and gradients) */ 
     232        int *nbRecvElement = new int[mpiSize]; 
     233        MPI_Alltoall(nbSendElement, 1, MPI_INT, nbRecvElement, 1, MPI_INT, communicator); 
     234 
     235        /* communicate indices of source elements on other ranks whoes value and gradient we need (since intersection) */ 
     236        int nbSendRequest = 0; 
     237        int nbRecvRequest = 0; 
     238        int **recvElement = new int*[mpiSize]; 
     239        double **sendValue = new double*[mpiSize]; 
     240        double **sendArea = new double*[mpiSize]; 
     241        Coord **sendGrad = new Coord*[mpiSize]; 
     242        GloId **sendNeighIds = new GloId*[mpiSize]; 
     243        MPI_Request *sendRequest = new MPI_Request[4*mpiSize]; 
     244        MPI_Request *recvRequest = new MPI_Request[4*mpiSize]; 
     245        for (int rank = 0; rank < mpiSize; rank++) 
     246        { 
     247                if (nbSendElement[rank] > 0) 
     248                { 
     249                        MPI_Issend(sendElement[rank], nbSendElement[rank], MPI_INT, rank, 0, communicator, &sendRequest[nbSendRequest]); 
     250                        nbSendRequest++; 
     251                } 
     252 
     253                if (nbRecvElement[rank] > 0) 
     254                { 
     255                        recvElement[rank] = new int[nbRecvElement[rank]]; 
     256                        sendValue[rank]   = new double[nbRecvElement[rank]]; 
     257                        sendArea[rank]   = new double[nbRecvElement[rank]]; 
     258                        if (order == 2) 
     259                        { 
     260                                sendNeighIds[rank] = new GloId[nbRecvElement[rank]*(NMAX+1)]; 
     261                                sendGrad[rank]    = new Coord[nbRecvElement[rank]*(NMAX+1)]; 
     262                        } 
     263                        else 
     264                        { 
     265                                sendNeighIds[rank] = new GloId[nbRecvElement[rank]]; 
     266                        } 
     267                        MPI_Irecv(recvElement[rank], nbRecvElement[rank], MPI_INT, rank, 0, communicator, &recvRequest[nbRecvRequest]); 
     268                        nbRecvRequest++; 
     269                } 
     270        } 
     271 
     272  MPI_Status *status = new MPI_Status[4*mpiSize]; 
     273 
     274  MPI_Waitall(nbSendRequest, sendRequest, status); 
     275  MPI_Waitall(nbRecvRequest, recvRequest, status); 
     276 
     277        /* for all indices that have been received from requesting ranks: pack values and gradients, then send */ 
     278        nbSendRequest = 0; 
     279        nbRecvRequest = 0; 
     280        for (int rank = 0; rank < mpiSize; rank++) 
     281        { 
     282                if (nbRecvElement[rank] > 0) 
     283                { 
     284                        int jj = 0; // jj == j if no weight writing 
     285                        for (int j = 0; j < nbRecvElement[rank]; j++) 
     286                        { 
     287                                sendValue[rank][j] = sstree.localElements[recvElement[rank][j]].val; 
     288                                sendArea[rank][j] = sstree.localElements[recvElement[rank][j]].area; 
     289                                if (order == 2) 
     290                                { 
     291                                        sendGrad[rank][jj] = sstree.localElements[recvElement[rank][j]].grad; 
     292                                        sendNeighIds[rank][jj] = sstree.localElements[recvElement[rank][j]].src_id; 
     293                                        jj++; 
     294                                        for (int i = 0; i < NMAX; i++) 
     295                                        { 
     296                                                sendGrad[rank][jj] = sstree.localElements[recvElement[rank][j]].gradNeigh[i]; 
     297            sendNeighIds[rank][jj] = sstree.localElements[recvElement[rank][j]].neighId[i]; 
     298                                                jj++; 
     299                                        } 
     300                                } 
     301                                else 
     302                                        sendNeighIds[rank][j] = sstree.localElements[recvElement[rank][j]].src_id; 
     303                        } 
     304                        MPI_Issend(sendValue[rank],  nbRecvElement[rank], MPI_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]); 
     305                        nbSendRequest++; 
     306                        MPI_Issend(sendArea[rank],  nbRecvElement[rank], MPI_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]); 
     307                        nbSendRequest++; 
     308                        if (order == 2) 
     309                        { 
     310                                MPI_Issend(sendGrad[rank], 3*nbRecvElement[rank]*(NMAX+1), 
     311                                                                MPI_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]); 
     312                                nbSendRequest++; 
     313                                MPI_Issend(sendNeighIds[rank], 4*nbRecvElement[rank]*(NMAX+1), MPI_INT, rank, 0, communicator, &sendRequest[nbSendRequest]); 
     314//ym  --> attention taille GloId 
     315                                nbSendRequest++; 
     316                        } 
     317                        else 
     318                        { 
     319                                MPI_Issend(sendNeighIds[rank], 4*nbRecvElement[rank], MPI_INT, rank, 0, communicator, &sendRequest[nbSendRequest]); 
     320//ym  --> attention taille GloId 
     321                                nbSendRequest++; 
     322                        } 
     323                } 
     324                if (nbSendElement[rank] > 0) 
     325                { 
     326                        MPI_Irecv(recvValue[rank],  nbSendElement[rank], MPI_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]); 
     327                        nbRecvRequest++; 
     328                        MPI_Irecv(recvArea[rank],  nbSendElement[rank], MPI_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]); 
     329                        nbRecvRequest++; 
     330                        if (order == 2) 
     331                        { 
     332                                MPI_Irecv(recvGrad[rank], 3*nbSendElement[rank]*(NMAX+1), 
     333                                                MPI_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]); 
     334                                nbRecvRequest++; 
     335                                MPI_Irecv(recvNeighIds[rank], 4*nbSendElement[rank]*(NMAX+1), MPI_INT, rank, 0, communicator, &recvRequest[nbRecvRequest]); 
     336//ym  --> attention taille GloId 
     337                                nbRecvRequest++; 
     338                        } 
     339                        else 
     340                        { 
     341                                MPI_Irecv(recvNeighIds[rank], 4*nbSendElement[rank], MPI_INT, rank, 0, communicator, &recvRequest[nbRecvRequest]); 
     342//ym  --> attention taille GloId 
     343                                nbRecvRequest++; 
     344                        } 
     345                } 
     346        } 
     347         
     348        MPI_Waitall(nbSendRequest, sendRequest, status); 
     349        MPI_Waitall(nbRecvRequest, recvRequest, status); 
     350         
     351 
     352        /* now that all values and gradients are available use them to computed interpolated values on target 
     353           and also to compute weights */ 
     354        int i = 0; 
     355        for (int j = 0; j < nbElements; j++) 
     356        { 
     357                Elt& e = elements[j]; 
     358 
     359                /* since for the 2nd order case source grid elements can contribute to a destination grid element over several "paths" 
     360                   (step1: gradient is computed using neighbours on same grid, step2: intersection uses several elements on other grid) 
     361                   accumulate them so that there is only one final weight between two elements */ 
     362                map<GloId,double> wgt_map; 
     363 
     364                /* for destination element `e` loop over all intersetions/the corresponding source elements */ 
     365                for (list<Polyg *>::iterator it = e.is.begin(); it != e.is.end(); it++) 
     366                { 
     367                        /* it is the intersection element, so it->x and it->area are barycentre and area of intersection element (super mesh) 
     368                        but it->id is id of the source element that it intersects */ 
     369                        int n1 = (*it)->id.ind; 
     370                        int rank = (*it)->id.rank; 
     371                        double fk = recvValue[rank][n1]; 
     372                        double srcArea = recvArea[rank][n1]; 
     373                        double w = (*it)->area; 
     374      if (quantity) w/=srcArea ; 
     375 
     376                        /* first order: src value times weight (weight = supermesh area), later divide by target area */ 
     377                        int kk = (order == 2) ? n1 * (NMAX + 1) : n1; 
     378                        GloId neighID = recvNeighIds[rank][kk]; 
     379                        wgt_map[neighID] += w; 
     380 
     381                        if (order == 2) 
     382                        { 
     383                                for (int k = 0; k < NMAX+1; k++) 
     384                                { 
     385                                        int kk = n1 * (NMAX + 1) + k; 
     386                                        GloId neighID = recvNeighIds[rank][kk]; 
     387                                        if (neighID.ind != -1)  wgt_map[neighID] += w * scalarprod(recvGrad[rank][kk], (*it)->x); 
     388                                } 
     389 
     390                        } 
     391                } 
     392 
     393    double renorm=0; 
     394    if (renormalize)  
     395      for (map<GloId,double>::iterator it = wgt_map.begin(); it != wgt_map.end(); it++) renorm+=it->second / e.area; 
     396    else renorm=1. ; 
     397 
     398    for (map<GloId,double>::iterator it = wgt_map.begin(); it != wgt_map.end(); it++) 
     399                { 
     400      if (quantity)  this->remapMatrix[i] = (it->second ) / renorm; 
     401                        else this->remapMatrix[i] = (it->second / e.area) / renorm; 
     402                        this->srcAddress[i] = it->first.ind; 
     403                        this->srcRank[i] = it->first.rank; 
     404                        this->dstAddress[i] = j; 
     405      this->sourceWeightId[i]= it->first.globalId ; 
     406      this->targetWeightId[i]= targetGlobalId[j] ; 
     407                        i++; 
     408                } 
     409        } 
     410 
     411        /* free all memory allocated in this function */ 
     412        for (int rank = 0; rank < mpiSize; rank++) 
     413        { 
     414                if (nbSendElement[rank] > 0) 
     415                { 
     416                        delete[] sendElement[rank]; 
     417                        delete[] recvValue[rank]; 
     418                        delete[] recvArea[rank]; 
     419                        if (order == 2) 
     420                        { 
     421                                delete[] recvGrad[rank]; 
     422                        } 
     423                        delete[] recvNeighIds[rank]; 
     424                } 
     425                if (nbRecvElement[rank] > 0) 
     426                { 
     427                        delete[] recvElement[rank]; 
     428                        delete[] sendValue[rank]; 
     429                        delete[] sendArea[rank]; 
     430                        if (order == 2) 
     431                                delete[] sendGrad[rank]; 
     432                        delete[] sendNeighIds[rank]; 
     433                } 
     434        } 
     435        delete[] status; 
     436        delete[] sendRequest; 
     437        delete[] recvRequest; 
     438        delete[] elementList; 
     439        delete[] nbSendElement; 
     440        delete[] nbRecvElement; 
     441        delete[] sendElement; 
     442        delete[] recvElement; 
     443        delete[] sendValue; 
     444        delete[] recvValue; 
     445        delete[] sendGrad; 
     446        delete[] recvGrad; 
     447        delete[] sendNeighIds; 
     448        delete[] recvNeighIds; 
     449        return i; 
    457450} 
    458451 
    459452void Mapper::computeGrads() 
    460453{ 
    461     /* array of pointers to collect local elements and elements received from other cpu */ 
    462     vector<Elt*> globalElements(sstree.nbLocalElements + nbNeighbourElements); 
    463     int index = 0; 
    464     for (int i = 0; i < sstree.nbLocalElements; i++, index++) 
    465         globalElements[index] = &(sstree.localElements[i]); 
    466     for (int i = 0; i < nbNeighbourElements; i++, index++) 
    467         globalElements[index] = &neighbourElements[i]; 
    468  
    469     update_baryc(sstree.localElements, sstree.nbLocalElements); 
    470     computeGradients(&globalElements[0], sstree.nbLocalElements); 
     454        /* array of pointers to collect local elements and elements received from other cpu */ 
     455        vector<Elt*> globalElements(sstree.nbLocalElements + nbNeighbourElements); 
     456        int index = 0; 
     457        for (int i = 0; i < sstree.nbLocalElements; i++, index++) 
     458                globalElements[index] = &(sstree.localElements[i]); 
     459        for (int i = 0; i < nbNeighbourElements; i++, index++) 
     460                globalElements[index] = &neighbourElements[i]; 
     461 
     462        update_baryc(sstree.localElements, sstree.nbLocalElements); 
     463        computeGradients(&globalElements[0], sstree.nbLocalElements); 
    471464} 
    472465 
    473466/** for each element of the source grid, finds all the neighbouring elements that share an edge 
    474   (filling array neighbourElements). This is used later to compute gradients */ 
     467    (filling array neighbourElements). This is used later to compute gradients */ 
    475468void Mapper::buildMeshTopology() 
    476469{ 
    477     int mpiSize, mpiRank; 
    478     MPI_Comm_size(communicator, &mpiSize); 
    479     MPI_Comm_rank(communicator, &mpiRank); 
    480  
    481     vector<Node> *routingList = new vector<Node>[mpiSize]; 
    482     vector<vector<int> > routes(sstree.localTree.leafs.size()); 
    483  
    484     sstree.routeIntersections(routes, sstree.localTree.leafs); 
    485  
    486     for (int i = 0; i < routes.size(); ++i) 
    487         for (int k = 0; k < routes[i].size(); ++k) 
    488             routingList[routes[i][k]].push_back(sstree.localTree.leafs[i]); 
    489     routingList[mpiRank].clear(); 
    490  
    491  
    492     CMPIRouting mpiRoute(communicator); 
    493     mpiRoute.init(routes); 
    494     int nRecv = mpiRoute.getTotalSourceElement(); 
    495  
    496     int *nbSendNode = new int[mpiSize]; 
    497     int *nbRecvNode = new int[mpiSize]; 
    498     int *sendMessageSize = new int[mpiSize]; 
    499     int *recvMessageSize = new int[mpiSize]; 
    500  
    501     for (int rank = 0; rank < mpiSize; rank++) 
    502     { 
    503         nbSendNode[rank] = routingList[rank].size(); 
    504         sendMessageSize[rank] = 0; 
    505         for (size_t j = 0; j < routingList[rank].size(); j++) 
    506         { 
    507             Elt *elt = (Elt *) (routingList[rank][j].data); 
    508             sendMessageSize[rank] += packedPolygonSize(*elt); 
    509         } 
    510     } 
    511  
    512     MPI_Alltoall(nbSendNode, 1, MPI_INT, nbRecvNode, 1, MPI_INT, communicator); 
    513     MPI_Alltoall(sendMessageSize, 1, MPI_INT, recvMessageSize, 1, MPI_INT, communicator); 
    514  
    515     char **sendBuffer = new char*[mpiSize]; 
    516     char **recvBuffer = new char*[mpiSize]; 
    517     int *pos = new int[mpiSize]; 
    518  
    519     for (int rank = 0; rank < mpiSize; rank++) 
    520     { 
    521         if (nbSendNode[rank] > 0) sendBuffer[rank] = new char[sendMessageSize[rank]]; 
    522         if (nbRecvNode[rank] > 0) recvBuffer[rank] = new char[recvMessageSize[rank]]; 
    523     } 
    524  
    525     for (int rank = 0; rank < mpiSize; rank++) 
    526     { 
    527         pos[rank] = 0; 
    528         for (size_t j = 0; j < routingList[rank].size(); j++) 
    529         { 
    530             Elt *elt = (Elt *) (routingList[rank][j].data); 
    531             packPolygon(*elt, sendBuffer[rank], pos[rank]); 
    532         } 
    533     } 
    534     delete [] routingList; 
    535  
    536  
    537     int nbSendRequest = 0; 
    538     int nbRecvRequest = 0; 
    539     MPI_Request *sendRequest = new MPI_Request[mpiSize]; 
    540     MPI_Request *recvRequest = new MPI_Request[mpiSize]; 
    541     MPI_Status  *status      = new MPI_Status[mpiSize]; 
    542  
    543     for (int rank = 0; rank < mpiSize; rank++) 
    544     { 
    545         if (nbSendNode[rank] > 0) 
    546         { 
    547             MPI_Issend(sendBuffer[rank], sendMessageSize[rank], MPI_CHAR, rank, 0, communicator, &sendRequest[nbSendRequest]); 
    548             nbSendRequest++; 
    549         } 
    550         if (nbRecvNode[rank] > 0) 
    551         { 
    552             MPI_Irecv(recvBuffer[rank], recvMessageSize[rank], MPI_CHAR, rank, 0, communicator, &recvRequest[nbRecvRequest]); 
    553             nbRecvRequest++; 
    554         } 
    555     } 
    556      
    557     MPI_Waitall(nbRecvRequest, recvRequest, status); 
    558     MPI_Waitall(nbSendRequest, sendRequest, status); 
    559  
    560     for (int rank = 0; rank < mpiSize; rank++) 
    561         if (nbSendNode[rank] > 0) delete [] sendBuffer[rank]; 
    562     delete [] sendBuffer; 
    563  
    564     char **sendBuffer2 = new char*[mpiSize]; 
    565     char **recvBuffer2 = new char*[mpiSize]; 
    566  
    567     for (int rank = 0; rank < mpiSize; rank++) 
    568     { 
    569         nbSendNode[rank] = 0; 
    570         sendMessageSize[rank] = 0; 
    571  
    572         if (nbRecvNode[rank] > 0) 
    573         { 
    574             set<NodePtr> neighbourList; 
    575             pos[rank] = 0; 
    576             for (int j = 0; j < nbRecvNode[rank]; j++) 
    577             { 
    578                 Elt elt; 
    579                 unpackPolygon(elt, recvBuffer[rank], pos[rank]); 
    580                 Node node(elt.x, cptRadius(elt), &elt); 
    581                 findNeighbour(sstree.localTree.root, &node, neighbourList); 
    582             } 
    583             nbSendNode[rank] = neighbourList.size(); 
    584             for (set<NodePtr>::iterator it = neighbourList.begin(); it != neighbourList.end(); it++) 
    585             { 
    586                 Elt *elt = (Elt *) ((*it)->data); 
    587                 sendMessageSize[rank] += packedPolygonSize(*elt); 
    588             } 
    589  
    590             sendBuffer2[rank] = new char[sendMessageSize[rank]]; 
    591             pos[rank] = 0; 
    592  
    593             for (set<NodePtr>::iterator it = neighbourList.begin(); it != neighbourList.end(); it++) 
    594             { 
    595                 Elt *elt = (Elt *) ((*it)->data); 
    596                 packPolygon(*elt, sendBuffer2[rank], pos[rank]); 
    597             } 
    598         } 
    599     } 
    600     for (int rank = 0; rank < mpiSize; rank++) 
    601         if (nbRecvNode[rank] > 0) delete [] recvBuffer[rank]; 
    602     delete [] recvBuffer; 
    603  
    604  
    605     MPI_Barrier(communicator); 
    606     MPI_Alltoall(nbSendNode, 1, MPI_INT, nbRecvNode, 1, MPI_INT, communicator); 
    607     MPI_Alltoall(sendMessageSize, 1, MPI_INT, recvMessageSize, 1, MPI_INT, communicator); 
    608  
    609     for (int rank = 0; rank < mpiSize; rank++) 
    610         if (nbRecvNode[rank] > 0) recvBuffer2[rank] = new char[recvMessageSize[rank]]; 
    611  
    612     nbSendRequest = 0; 
    613     nbRecvRequest = 0; 
    614  
    615     for (int rank = 0; rank < mpiSize; rank++) 
    616     { 
    617         if (nbSendNode[rank] > 0) 
    618         { 
    619             MPI_Issend(sendBuffer2[rank], sendMessageSize[rank], MPI_CHAR, rank, 0, communicator, &sendRequest[nbSendRequest]); 
    620             nbSendRequest++; 
    621         } 
    622         if (nbRecvNode[rank] > 0) 
    623         { 
    624             MPI_Irecv(recvBuffer2[rank], recvMessageSize[rank], MPI_CHAR, rank, 0, communicator, &recvRequest[nbRecvRequest]); 
    625             nbRecvRequest++; 
    626         } 
    627     } 
    628  
    629     MPI_Waitall(nbRecvRequest, recvRequest, status); 
    630     MPI_Waitall(nbSendRequest, sendRequest, status); 
    631  
    632     int nbNeighbourNodes = 0; 
    633     for (int rank = 0; rank < mpiSize; rank++) 
    634         nbNeighbourNodes += nbRecvNode[rank]; 
    635  
    636     neighbourElements = new Elt[nbNeighbourNodes]; 
    637     nbNeighbourElements = nbNeighbourNodes; 
    638  
    639     int index = 0; 
    640     for (int rank = 0; rank < mpiSize; rank++) 
    641     { 
    642         pos[rank] = 0; 
    643         for (int j = 0; j < nbRecvNode[rank]; j++) 
    644         { 
    645             unpackPolygon(neighbourElements[index], recvBuffer2[rank], pos[rank]); 
    646             neighbourElements[index].id.ind = sstree.localTree.leafs.size() + index; 
    647             index++; 
    648         } 
    649     } 
    650     for (int rank = 0; rank < mpiSize; rank++) 
    651     { 
    652         if (nbRecvNode[rank] > 0) delete [] recvBuffer2[rank]; 
    653         if (nbSendNode[rank] > 0) delete [] sendBuffer2[rank]; 
    654     } 
    655     delete [] recvBuffer2; 
    656     delete [] sendBuffer2; 
    657     delete [] sendMessageSize; 
    658     delete [] recvMessageSize; 
    659     delete [] nbSendNode; 
    660     delete [] nbRecvNode; 
    661     delete [] sendRequest; 
    662     delete [] recvRequest; 
    663     delete [] status; 
    664     delete [] pos; 
    665  
    666     /* re-compute on received elements to avoid having to send this information */ 
    667     neighbourNodes.resize(nbNeighbourNodes); 
    668     setCirclesAndLinks(neighbourElements, neighbourNodes); 
    669     cptAllEltsGeom(neighbourElements, nbNeighbourNodes, srcGrid.pole); 
    670  
    671     /* the local SS tree must include nodes from other cpus if they are potential 
    672        intersector of a local node */ 
    673     sstree.localTree.insertNodes(neighbourNodes); 
    674  
    675     /* for every local element, 
    676        use the SS-tree to find all elements (including neighbourElements) 
    677        who are potential neighbours because their circles intersect, 
    678        then check all canditates for common edges to build up connectivity information 
    679        */ 
    680     for (int j = 0; j < sstree.localTree.leafs.size(); j++) 
    681     { 
    682         Node& node = sstree.localTree.leafs[j]; 
    683  
    684         /* find all leafs whoes circles that intersect node's circle and save into node->intersectors */ 
    685         node.search(sstree.localTree.root); 
    686  
    687         Elt *elt = (Elt *)(node.data); 
    688  
    689         for (int i = 0; i < elt->n; i++) elt->neighbour[i] = NOT_FOUND; 
    690  
    691         /* for element `elt` loop through all nodes in the SS-tree 
    692            whoes circles intersect with the circle around `elt` (the SS intersectors) 
    693            and check if they are neighbours in the sense that the two elements share an edge. 
    694            If they do, save this information for elt */ 
    695         for (list<NodePtr>::iterator it = (node.intersectors).begin(); it != (node.intersectors).end(); ++it) 
    696         { 
    697             Elt *elt2 = (Elt *)((*it)->data); 
    698             set_neighbour(*elt, *elt2); 
    699         } 
    700  
    701     } 
     470        int mpiSize, mpiRank; 
     471        MPI_Comm_size(communicator, &mpiSize); 
     472        MPI_Comm_rank(communicator, &mpiRank); 
     473 
     474        vector<Node> *routingList = new vector<Node>[mpiSize]; 
     475        vector<vector<int> > routes(sstree.localTree.leafs.size()); 
     476 
     477        sstree.routeIntersections(routes, sstree.localTree.leafs); 
     478 
     479        for (int i = 0; i < routes.size(); ++i) 
     480                for (int k = 0; k < routes[i].size(); ++k) 
     481                        routingList[routes[i][k]].push_back(sstree.localTree.leafs[i]); 
     482        routingList[mpiRank].clear(); 
     483 
     484 
     485        CMPIRouting mpiRoute(communicator); 
     486        mpiRoute.init(routes); 
     487        int nRecv = mpiRoute.getTotalSourceElement(); 
     488 
     489        int *nbSendNode = new int[mpiSize]; 
     490        int *nbRecvNode = new int[mpiSize]; 
     491        int *sendMessageSize = new int[mpiSize]; 
     492        int *recvMessageSize = new int[mpiSize]; 
     493 
     494        for (int rank = 0; rank < mpiSize; rank++) 
     495        { 
     496                nbSendNode[rank] = routingList[rank].size(); 
     497                sendMessageSize[rank] = 0; 
     498                for (size_t j = 0; j < routingList[rank].size(); j++) 
     499                { 
     500                        Elt *elt = (Elt *) (routingList[rank][j].data); 
     501                        sendMessageSize[rank] += packedPolygonSize(*elt); 
     502                } 
     503        } 
     504 
     505        MPI_Alltoall(nbSendNode, 1, MPI_INT, nbRecvNode, 1, MPI_INT, communicator); 
     506        MPI_Alltoall(sendMessageSize, 1, MPI_INT, recvMessageSize, 1, MPI_INT, communicator); 
     507 
     508        char **sendBuffer = new char*[mpiSize]; 
     509        char **recvBuffer = new char*[mpiSize]; 
     510        int *pos = new int[mpiSize]; 
     511 
     512        for (int rank = 0; rank < mpiSize; rank++) 
     513        { 
     514                if (nbSendNode[rank] > 0) sendBuffer[rank] = new char[sendMessageSize[rank]]; 
     515                if (nbRecvNode[rank] > 0) recvBuffer[rank] = new char[recvMessageSize[rank]]; 
     516        } 
     517 
     518        for (int rank = 0; rank < mpiSize; rank++) 
     519        { 
     520                pos[rank] = 0; 
     521                for (size_t j = 0; j < routingList[rank].size(); j++) 
     522                { 
     523                        Elt *elt = (Elt *) (routingList[rank][j].data); 
     524                        packPolygon(*elt, sendBuffer[rank], pos[rank]); 
     525                } 
     526        } 
     527        delete [] routingList; 
     528 
     529 
     530        int nbSendRequest = 0; 
     531        int nbRecvRequest = 0; 
     532        MPI_Request *sendRequest = new MPI_Request[mpiSize]; 
     533        MPI_Request *recvRequest = new MPI_Request[mpiSize]; 
     534        MPI_Status  *status      = new MPI_Status[mpiSize]; 
     535 
     536        for (int rank = 0; rank < mpiSize; rank++) 
     537        { 
     538                if (nbSendNode[rank] > 0) 
     539                { 
     540                        MPI_Issend(sendBuffer[rank], sendMessageSize[rank], MPI_CHAR, rank, 0, communicator, &sendRequest[nbSendRequest]); 
     541                        nbSendRequest++; 
     542                } 
     543                if (nbRecvNode[rank] > 0) 
     544                { 
     545                        MPI_Irecv(recvBuffer[rank], recvMessageSize[rank], MPI_CHAR, rank, 0, communicator, &recvRequest[nbRecvRequest]); 
     546                        nbRecvRequest++; 
     547                } 
     548        } 
     549 
     550  MPI_Waitall(nbRecvRequest, recvRequest, status); 
     551  MPI_Waitall(nbSendRequest, sendRequest, status); 
     552 
     553        for (int rank = 0; rank < mpiSize; rank++) 
     554                if (nbSendNode[rank] > 0) delete [] sendBuffer[rank]; 
     555        delete [] sendBuffer; 
     556 
     557        char **sendBuffer2 = new char*[mpiSize]; 
     558        char **recvBuffer2 = new char*[mpiSize]; 
     559 
     560        for (int rank = 0; rank < mpiSize; rank++) 
     561        { 
     562                nbSendNode[rank] = 0; 
     563                sendMessageSize[rank] = 0; 
     564 
     565                if (nbRecvNode[rank] > 0) 
     566                { 
     567                        set<NodePtr> neighbourList; 
     568                        pos[rank] = 0; 
     569                        for (int j = 0; j < nbRecvNode[rank]; j++) 
     570                        { 
     571                                Elt elt; 
     572                                unpackPolygon(elt, recvBuffer[rank], pos[rank]); 
     573                                Node node(elt.x, cptRadius(elt), &elt); 
     574                                findNeighbour(sstree.localTree.root, &node, neighbourList); 
     575                        } 
     576                        nbSendNode[rank] = neighbourList.size(); 
     577                        for (set<NodePtr>::iterator it = neighbourList.begin(); it != neighbourList.end(); it++) 
     578                        { 
     579                                Elt *elt = (Elt *) ((*it)->data); 
     580                                sendMessageSize[rank] += packedPolygonSize(*elt); 
     581                        } 
     582 
     583                        sendBuffer2[rank] = new char[sendMessageSize[rank]]; 
     584                        pos[rank] = 0; 
     585 
     586                        for (set<NodePtr>::iterator it = neighbourList.begin(); it != neighbourList.end(); it++) 
     587                        { 
     588                                Elt *elt = (Elt *) ((*it)->data); 
     589                                packPolygon(*elt, sendBuffer2[rank], pos[rank]); 
     590                        } 
     591                } 
     592        } 
     593        for (int rank = 0; rank < mpiSize; rank++) 
     594                if (nbRecvNode[rank] > 0) delete [] recvBuffer[rank]; 
     595        delete [] recvBuffer; 
     596 
     597 
     598        MPI_Barrier(communicator); 
     599        MPI_Alltoall(nbSendNode, 1, MPI_INT, nbRecvNode, 1, MPI_INT, communicator); 
     600        MPI_Alltoall(sendMessageSize, 1, MPI_INT, recvMessageSize, 1, MPI_INT, communicator); 
     601 
     602        for (int rank = 0; rank < mpiSize; rank++) 
     603                if (nbRecvNode[rank] > 0) recvBuffer2[rank] = new char[recvMessageSize[rank]]; 
     604 
     605        nbSendRequest = 0; 
     606        nbRecvRequest = 0; 
     607 
     608        for (int rank = 0; rank < mpiSize; rank++) 
     609        { 
     610                if (nbSendNode[rank] > 0) 
     611                { 
     612                        MPI_Issend(sendBuffer2[rank], sendMessageSize[rank], MPI_CHAR, rank, 0, communicator, &sendRequest[nbSendRequest]); 
     613                        nbSendRequest++; 
     614                } 
     615                if (nbRecvNode[rank] > 0) 
     616                { 
     617                        MPI_Irecv(recvBuffer2[rank], recvMessageSize[rank], MPI_CHAR, rank, 0, communicator, &recvRequest[nbRecvRequest]); 
     618                        nbRecvRequest++; 
     619                } 
     620        } 
     621 
     622        MPI_Waitall(nbRecvRequest, recvRequest, status); 
     623        MPI_Waitall(nbSendRequest, sendRequest, status); 
     624 
     625        int nbNeighbourNodes = 0; 
     626        for (int rank = 0; rank < mpiSize; rank++) 
     627                nbNeighbourNodes += nbRecvNode[rank]; 
     628 
     629        neighbourElements = new Elt[nbNeighbourNodes]; 
     630        nbNeighbourElements = nbNeighbourNodes; 
     631 
     632        int index = 0; 
     633        for (int rank = 0; rank < mpiSize; rank++) 
     634        { 
     635                pos[rank] = 0; 
     636                for (int j = 0; j < nbRecvNode[rank]; j++) 
     637                { 
     638                        unpackPolygon(neighbourElements[index], recvBuffer2[rank], pos[rank]); 
     639                        neighbourElements[index].id.ind = sstree.localTree.leafs.size() + index; 
     640                        index++; 
     641                } 
     642        } 
     643        for (int rank = 0; rank < mpiSize; rank++) 
     644        { 
     645                if (nbRecvNode[rank] > 0) delete [] recvBuffer2[rank]; 
     646                if (nbSendNode[rank] > 0) delete [] sendBuffer2[rank]; 
     647        } 
     648        delete [] recvBuffer2; 
     649        delete [] sendBuffer2; 
     650        delete [] sendMessageSize; 
     651        delete [] recvMessageSize; 
     652        delete [] nbSendNode; 
     653        delete [] nbRecvNode; 
     654        delete [] sendRequest; 
     655        delete [] recvRequest; 
     656        delete [] status; 
     657        delete [] pos; 
     658 
     659        /* re-compute on received elements to avoid having to send this information */ 
     660        neighbourNodes.resize(nbNeighbourNodes); 
     661        setCirclesAndLinks(neighbourElements, neighbourNodes); 
     662        cptAllEltsGeom(neighbourElements, nbNeighbourNodes, srcGrid.pole); 
     663 
     664        /* the local SS tree must include nodes from other cpus if they are potential 
     665           intersector of a local node */ 
     666        sstree.localTree.insertNodes(neighbourNodes); 
     667 
     668        /* for every local element, 
     669           use the SS-tree to find all elements (including neighbourElements) 
     670           who are potential neighbours because their circles intersect, 
     671           then check all canditates for common edges to build up connectivity information 
     672        */ 
     673        for (int j = 0; j < sstree.localTree.leafs.size(); j++) 
     674        { 
     675                Node& node = sstree.localTree.leafs[j]; 
     676 
     677                /* find all leafs whoes circles that intersect node's circle and save into node->intersectors */ 
     678                node.search(sstree.localTree.root); 
     679 
     680                Elt *elt = (Elt *)(node.data); 
     681 
     682                for (int i = 0; i < elt->n; i++) elt->neighbour[i] = NOT_FOUND; 
     683 
     684                /* for element `elt` loop through all nodes in the SS-tree 
     685                   whoes circles intersect with the circle around `elt` (the SS intersectors) 
     686                   and check if they are neighbours in the sense that the two elements share an edge. 
     687                   If they do, save this information for elt */ 
     688                for (list<NodePtr>::iterator it = (node.intersectors).begin(); it != (node.intersectors).end(); ++it) 
     689                { 
     690                        Elt *elt2 = (Elt *)((*it)->data); 
     691                        set_neighbour(*elt, *elt2); 
     692                } 
     693 
     694        } 
    702695} 
    703696 
     
    705698void Mapper::computeIntersection(Elt *elements, int nbElements) 
    706699{ 
    707     int mpiSize, mpiRank; 
    708     MPI_Comm_size(communicator, &mpiSize); 
    709     MPI_Comm_rank(communicator, &mpiRank); 
    710  
    711     MPI_Barrier(communicator); 
    712  
    713     vector<Node> *routingList = new vector<Node>[mpiSize]; 
    714  
    715     vector<Node> routeNodes;  routeNodes.reserve(nbElements); 
    716     for (int j = 0; j < nbElements; j++) 
    717     { 
    718         elements[j].id.ind = j; 
    719         elements[j].id.rank = mpiRank; 
    720         routeNodes.push_back(Node(elements[j].x, cptRadius(elements[j]), &elements[j])); 
    721     } 
    722  
    723     vector<vector<int> > routes(routeNodes.size()); 
    724     sstree.routeIntersections(routes, routeNodes); 
    725     for (int i = 0; i < routes.size(); ++i) 
    726         for (int k = 0; k < routes[i].size(); ++k) 
    727             routingList[routes[i][k]].push_back(routeNodes[i]); 
    728  
    729     if (verbose >= 2) 
    730     { 
    731         cout << " --> rank  " << mpiRank << " nbElements " << nbElements << " : "; 
    732         for (int rank = 0; rank < mpiSize; rank++) 
    733             cout << routingList[rank].size() << "   "; 
    734         cout << endl; 
    735     } 
    736     MPI_Barrier(communicator); 
    737  
    738     int *nbSendNode = new int[mpiSize]; 
    739     int *nbRecvNode = new int[mpiSize]; 
    740     int *sentMessageSize = new int[mpiSize]; 
    741     int *recvMessageSize = new int[mpiSize]; 
    742  
    743     for (int rank = 0; rank < mpiSize; rank++) 
    744     { 
    745         nbSendNode[rank] = routingList[rank].size(); 
    746         sentMessageSize[rank] = 0; 
    747         for (size_t j = 0; j < routingList[rank].size(); j++) 
    748         { 
    749             Elt *elt = (Elt *) (routingList[rank][j].data); 
    750             sentMessageSize[rank] += packedPolygonSize(*elt); 
    751         } 
    752     } 
    753  
    754     MPI_Alltoall(nbSendNode, 1, MPI_INT, nbRecvNode, 1, MPI_INT, communicator); 
    755     MPI_Alltoall(sentMessageSize, 1, MPI_INT, recvMessageSize, 1, MPI_INT, communicator); 
    756  
    757     int total = 0; 
    758  
    759     for (int rank = 0; rank < mpiSize; rank++) 
    760     { 
    761         total = total + nbRecvNode[rank]; 
    762     } 
    763  
    764     if (verbose >= 2) cout << "---> rank " << mpiRank << " : compute intersection : total received nodes  " << total << endl; 
    765     char **sendBuffer = new char*[mpiSize]; 
    766     char **recvBuffer = new char*[mpiSize]; 
    767     int *pos = new int[mpiSize]; 
    768  
    769     for (int rank = 0; rank < mpiSize; rank++) 
    770     { 
    771         if (nbSendNode[rank] > 0) sendBuffer[rank] = new char[sentMessageSize[rank]]; 
    772         if (nbRecvNode[rank] > 0) recvBuffer[rank] = new char[recvMessageSize[rank]]; 
    773     } 
    774  
    775     for (int rank = 0; rank < mpiSize; rank++) 
    776     { 
    777         pos[rank] = 0; 
    778         for (size_t j = 0; j < routingList[rank].size(); j++) 
    779         { 
    780             Elt* elt = (Elt *) (routingList[rank][j].data); 
    781             packPolygon(*elt, sendBuffer[rank], pos[rank]); 
    782         } 
    783     } 
    784     delete [] routingList; 
    785  
    786     int nbSendRequest = 0; 
    787     int nbRecvRequest = 0; 
    788     MPI_Request *sendRequest = new MPI_Request[mpiSize]; 
    789     MPI_Request *recvRequest = new MPI_Request[mpiSize]; 
    790     MPI_Status   *status = new MPI_Status[mpiSize]; 
    791  
    792     for (int rank = 0; rank < mpiSize; rank++) 
    793     { 
    794         if (nbSendNode[rank] > 0) 
    795         { 
    796             MPI_Issend(sendBuffer[rank], sentMessageSize[rank], MPI_CHAR, rank, 0, communicator, &sendRequest[nbSendRequest]); 
    797             nbSendRequest++; 
    798         } 
    799         if (nbRecvNode[rank] > 0) 
    800         { 
    801             MPI_Irecv(recvBuffer[rank], recvMessageSize[rank], MPI_CHAR, rank, 0, communicator, &recvRequest[nbRecvRequest]); 
    802             nbRecvRequest++; 
    803         } 
    804     } 
    805  
    806     MPI_Waitall(nbRecvRequest, recvRequest, status); 
    807     MPI_Waitall(nbSendRequest, sendRequest, status); 
    808  
    809     char **sendBuffer2 = new char*[mpiSize]; 
    810     char **recvBuffer2 = new char*[mpiSize]; 
    811  
    812     double tic = cputime(); 
    813     for (int rank = 0; rank < mpiSize; rank++) 
    814     { 
    815         sentMessageSize[rank] = 0; 
    816  
    817         if (nbRecvNode[rank] > 0) 
    818         { 
    819             Elt *recvElt = new Elt[nbRecvNode[rank]]; 
    820             pos[rank] = 0; 
    821             for (int j = 0; j < nbRecvNode[rank]; j++) 
    822             { 
    823                 unpackPolygon(recvElt[j], recvBuffer[rank], pos[rank]); 
    824                 cptEltGeom(recvElt[j], tgtGrid.pole); 
    825                 Node recvNode(recvElt[j].x, cptRadius(recvElt[j]), &recvElt[j]); 
    826                 recvNode.search(sstree.localTree.root); 
    827                 /* for a node holding an element of the target, loop throught candidates for intersecting source */ 
    828                 for (list<NodePtr>::iterator it = (recvNode.intersectors).begin(); it != (recvNode.intersectors).end(); ++it) 
    829                 { 
    830                     Elt *elt2 = (Elt *) ((*it)->data); 
    831                     /* recvElt is target, elt2 is source */ 
    832                     intersect(&recvElt[j], elt2); 
    833                     //intersect_ym(&recvElt[j], elt2); 
    834                 } 
    835                 if (recvElt[j].is.size() > 0) sentMessageSize[rank] += packIntersectionSize(recvElt[j]); 
    836  
    837                 // here recvNode goes out of scope 
    838             } 
    839  
    840             if (sentMessageSize[rank] > 0) 
    841             { 
    842                 sentMessageSize[rank] += sizeof(int); 
    843                 sendBuffer2[rank] = new char[sentMessageSize[rank]]; 
    844                 *((int *) sendBuffer2[rank]) = 0; 
    845                 pos[rank] = sizeof(int); 
    846                 for (int j = 0; j < nbRecvNode[rank]; j++) 
    847                 { 
    848                     packIntersection(recvElt[j], sendBuffer2[rank], pos[rank]); 
    849                     //FIXME should be deleted: recvElt[j].delete_intersections(); // intersection areas have been packed to buffer and won't be used any more 
    850                 } 
    851             } 
    852             delete [] recvElt; 
    853         } 
    854     } 
    855     delete [] pos; 
    856  
    857     for (int rank = 0; rank < mpiSize; rank++) 
    858     { 
    859         if (nbSendNode[rank] > 0) delete [] sendBuffer[rank]; 
    860         if (nbRecvNode[rank] > 0) delete [] recvBuffer[rank]; 
    861         nbSendNode[rank] = 0; 
    862     } 
    863  
    864     if (verbose >= 2) cout << "Rank " << mpiRank << "  Compute (internal) intersection " << cputime() - tic << " s" << endl; 
    865     MPI_Alltoall(sentMessageSize, 1, MPI_INT, recvMessageSize, 1, MPI_INT, communicator); 
    866  
    867     for (int rank = 0; rank < mpiSize; rank++) 
    868         if (recvMessageSize[rank] > 0) 
    869             recvBuffer2[rank] = new char[recvMessageSize[rank]]; 
    870  
    871     nbSendRequest = 0; 
    872     nbRecvRequest = 0; 
    873  
    874     for (int rank = 0; rank < mpiSize; rank++) 
    875     { 
    876         if (sentMessageSize[rank] > 0) 
    877         { 
    878             MPI_Issend(sendBuffer2[rank], sentMessageSize[rank], MPI_CHAR, rank, 0, communicator, &sendRequest[nbSendRequest]); 
    879             nbSendRequest++; 
    880         } 
    881         if (recvMessageSize[rank] > 0) 
    882         { 
    883             MPI_Irecv(recvBuffer2[rank], recvMessageSize[rank], MPI_CHAR, rank, 0, communicator, &recvRequest[nbRecvRequest]); 
    884             nbRecvRequest++; 
    885         } 
    886     } 
    887      
    888     MPI_Waitall(nbRecvRequest, recvRequest, status); 
    889     MPI_Waitall(nbSendRequest, sendRequest, status); 
    890  
    891     delete [] sendRequest; 
    892     delete [] recvRequest; 
    893     delete [] status; 
    894     for (int rank = 0; rank < mpiSize; rank++) 
    895     { 
    896         if (nbRecvNode[rank] > 0) 
    897         { 
    898             if (sentMessageSize[rank] > 0) 
    899                 delete [] sendBuffer2[rank]; 
    900         } 
    901  
    902         if (recvMessageSize[rank] > 0) 
    903         { 
    904             unpackIntersection(elements, recvBuffer2[rank]); 
    905             delete [] recvBuffer2[rank]; 
    906         } 
    907     } 
    908     delete [] sendBuffer2; 
    909     delete [] recvBuffer2; 
    910     delete [] sendBuffer; 
    911     delete [] recvBuffer; 
    912  
    913     delete [] nbSendNode; 
    914     delete [] nbRecvNode; 
    915     delete [] sentMessageSize; 
    916     delete [] recvMessageSize; 
     700        int mpiSize, mpiRank; 
     701        MPI_Comm_size(communicator, &mpiSize); 
     702        MPI_Comm_rank(communicator, &mpiRank); 
     703 
     704        MPI_Barrier(communicator); 
     705 
     706        vector<Node> *routingList = new vector<Node>[mpiSize]; 
     707 
     708        vector<Node> routeNodes;  routeNodes.reserve(nbElements); 
     709        for (int j = 0; j < nbElements; j++) 
     710        { 
     711                elements[j].id.ind = j; 
     712                elements[j].id.rank = mpiRank; 
     713                routeNodes.push_back(Node(elements[j].x, cptRadius(elements[j]), &elements[j])); 
     714        } 
     715 
     716        vector<vector<int> > routes(routeNodes.size()); 
     717        sstree.routeIntersections(routes, routeNodes); 
     718        for (int i = 0; i < routes.size(); ++i) 
     719                for (int k = 0; k < routes[i].size(); ++k) 
     720                        routingList[routes[i][k]].push_back(routeNodes[i]); 
     721 
     722        if (verbose >= 2) 
     723        { 
     724                cout << " --> rank  " << mpiRank << " nbElements " << nbElements << " : "; 
     725                for (int rank = 0; rank < mpiSize; rank++) 
     726                        cout << routingList[rank].size() << "   "; 
     727                cout << endl; 
     728        } 
     729        MPI_Barrier(communicator); 
     730 
     731        int *nbSendNode = new int[mpiSize]; 
     732        int *nbRecvNode = new int[mpiSize]; 
     733        int *sentMessageSize = new int[mpiSize]; 
     734        int *recvMessageSize = new int[mpiSize]; 
     735 
     736        for (int rank = 0; rank < mpiSize; rank++) 
     737        { 
     738                nbSendNode[rank] = routingList[rank].size(); 
     739                sentMessageSize[rank] = 0; 
     740                for (size_t j = 0; j < routingList[rank].size(); j++) 
     741                { 
     742                        Elt *elt = (Elt *) (routingList[rank][j].data); 
     743                        sentMessageSize[rank] += packedPolygonSize(*elt); 
     744                } 
     745        } 
     746 
     747        MPI_Alltoall(nbSendNode, 1, MPI_INT, nbRecvNode, 1, MPI_INT, communicator); 
     748        MPI_Alltoall(sentMessageSize, 1, MPI_INT, recvMessageSize, 1, MPI_INT, communicator); 
     749 
     750        int total = 0; 
     751 
     752        for (int rank = 0; rank < mpiSize; rank++) 
     753        { 
     754                total = total + nbRecvNode[rank]; 
     755        } 
     756 
     757        if (verbose >= 2) cout << "---> rank " << mpiRank << " : compute intersection : total received nodes  " << total << endl; 
     758        char **sendBuffer = new char*[mpiSize]; 
     759        char **recvBuffer = new char*[mpiSize]; 
     760        int *pos = new int[mpiSize]; 
     761 
     762        for (int rank = 0; rank < mpiSize; rank++) 
     763        { 
     764                if (nbSendNode[rank] > 0) sendBuffer[rank] = new char[sentMessageSize[rank]]; 
     765                if (nbRecvNode[rank] > 0) recvBuffer[rank] = new char[recvMessageSize[rank]]; 
     766        } 
     767 
     768        for (int rank = 0; rank < mpiSize; rank++) 
     769        { 
     770                pos[rank] = 0; 
     771                for (size_t j = 0; j < routingList[rank].size(); j++) 
     772                { 
     773                        Elt* elt = (Elt *) (routingList[rank][j].data); 
     774                        packPolygon(*elt, sendBuffer[rank], pos[rank]); 
     775                } 
     776        } 
     777        delete [] routingList; 
     778 
     779        int nbSendRequest = 0; 
     780        int nbRecvRequest = 0; 
     781        MPI_Request *sendRequest = new MPI_Request[mpiSize]; 
     782        MPI_Request *recvRequest = new MPI_Request[mpiSize]; 
     783        MPI_Status   *status = new MPI_Status[mpiSize]; 
     784 
     785        for (int rank = 0; rank < mpiSize; rank++) 
     786        { 
     787                if (nbSendNode[rank] > 0) 
     788                { 
     789                        MPI_Issend(sendBuffer[rank], sentMessageSize[rank], MPI_CHAR, rank, 0, communicator, &sendRequest[nbSendRequest]); 
     790                        nbSendRequest++; 
     791                } 
     792                if (nbRecvNode[rank] > 0) 
     793                { 
     794                        MPI_Irecv(recvBuffer[rank], recvMessageSize[rank], MPI_CHAR, rank, 0, communicator, &recvRequest[nbRecvRequest]); 
     795                        nbRecvRequest++; 
     796                } 
     797        } 
     798 
     799        MPI_Waitall(nbRecvRequest, recvRequest, status); 
     800        MPI_Waitall(nbSendRequest, sendRequest, status); 
     801 
     802        char **sendBuffer2 = new char*[mpiSize]; 
     803        char **recvBuffer2 = new char*[mpiSize]; 
     804 
     805        double tic = cputime(); 
     806        for (int rank = 0; rank < mpiSize; rank++) 
     807        { 
     808                sentMessageSize[rank] = 0; 
     809 
     810                if (nbRecvNode[rank] > 0) 
     811                { 
     812                        Elt *recvElt = new Elt[nbRecvNode[rank]]; 
     813                        pos[rank] = 0; 
     814                        for (int j = 0; j < nbRecvNode[rank]; j++) 
     815                        { 
     816                                unpackPolygon(recvElt[j], recvBuffer[rank], pos[rank]); 
     817                                cptEltGeom(recvElt[j], tgtGrid.pole); 
     818                                Node recvNode(recvElt[j].x, cptRadius(recvElt[j]), &recvElt[j]); 
     819                                recvNode.search(sstree.localTree.root); 
     820                                /* for a node holding an element of the target, loop throught candidates for intersecting source */ 
     821                                for (list<NodePtr>::iterator it = (recvNode.intersectors).begin(); it != (recvNode.intersectors).end(); ++it) 
     822                                { 
     823                                        Elt *elt2 = (Elt *) ((*it)->data); 
     824                                        /* recvElt is target, elt2 is source */ 
     825//                                      intersect(&recvElt[j], elt2); 
     826                                        intersect_ym(&recvElt[j], elt2); 
     827                                } 
     828                                if (recvElt[j].is.size() > 0) sentMessageSize[rank] += packIntersectionSize(recvElt[j]); 
     829 
     830                                // here recvNode goes out of scope 
     831                        } 
     832 
     833                        if (sentMessageSize[rank] > 0) 
     834                        { 
     835                                sentMessageSize[rank] += sizeof(int); 
     836                                sendBuffer2[rank] = new char[sentMessageSize[rank]]; 
     837                                *((int *) sendBuffer2[rank]) = 0; 
     838                                pos[rank] = sizeof(int); 
     839                                for (int j = 0; j < nbRecvNode[rank]; j++) 
     840                                { 
     841                                        packIntersection(recvElt[j], sendBuffer2[rank], pos[rank]); 
     842                                        //FIXME should be deleted: recvElt[j].delete_intersections(); // intersection areas have been packed to buffer and won't be used any more 
     843                                } 
     844                        } 
     845                        delete [] recvElt; 
     846                } 
     847        } 
     848        delete [] pos; 
     849 
     850        for (int rank = 0; rank < mpiSize; rank++) 
     851        { 
     852                if (nbSendNode[rank] > 0) delete [] sendBuffer[rank]; 
     853                if (nbRecvNode[rank] > 0) delete [] recvBuffer[rank]; 
     854                nbSendNode[rank] = 0; 
     855        } 
     856 
     857        if (verbose >= 2) cout << "Rank " << mpiRank << "  Compute (internal) intersection " << cputime() - tic << " s" << endl; 
     858        MPI_Alltoall(sentMessageSize, 1, MPI_INT, recvMessageSize, 1, MPI_INT, communicator); 
     859 
     860        for (int rank = 0; rank < mpiSize; rank++) 
     861                if (recvMessageSize[rank] > 0) 
     862                        recvBuffer2[rank] = new char[recvMessageSize[rank]]; 
     863 
     864        nbSendRequest = 0; 
     865        nbRecvRequest = 0; 
     866 
     867        for (int rank = 0; rank < mpiSize; rank++) 
     868        { 
     869                if (sentMessageSize[rank] > 0) 
     870                { 
     871                        MPI_Issend(sendBuffer2[rank], sentMessageSize[rank], MPI_CHAR, rank, 0, communicator, &sendRequest[nbSendRequest]); 
     872                        nbSendRequest++; 
     873                } 
     874                if (recvMessageSize[rank] > 0) 
     875                { 
     876                        MPI_Irecv(recvBuffer2[rank], recvMessageSize[rank], MPI_CHAR, rank, 0, communicator, &recvRequest[nbRecvRequest]); 
     877                        nbRecvRequest++; 
     878                } 
     879        } 
     880 
     881  MPI_Waitall(nbRecvRequest, recvRequest, status); 
     882  MPI_Waitall(nbSendRequest, sendRequest, status); 
     883 
     884        delete [] sendRequest; 
     885        delete [] recvRequest; 
     886        delete [] status; 
     887        for (int rank = 0; rank < mpiSize; rank++) 
     888        { 
     889                if (nbRecvNode[rank] > 0) 
     890                { 
     891                        if (sentMessageSize[rank] > 0) 
     892                                delete [] sendBuffer2[rank]; 
     893                } 
     894 
     895                if (recvMessageSize[rank] > 0) 
     896                { 
     897                        unpackIntersection(elements, recvBuffer2[rank]); 
     898                        delete [] recvBuffer2[rank]; 
     899                } 
     900        } 
     901        delete [] sendBuffer2; 
     902        delete [] recvBuffer2; 
     903        delete [] sendBuffer; 
     904        delete [] recvBuffer; 
     905 
     906        delete [] nbSendNode; 
     907        delete [] nbRecvNode; 
     908        delete [] sentMessageSize; 
     909        delete [] recvMessageSize; 
    917910} 
    918911 
    919912Mapper::~Mapper() 
    920913{ 
    921     delete [] remapMatrix; 
    922     delete [] srcAddress; 
    923     delete [] srcRank; 
    924     delete [] dstAddress; 
    925     if (neighbourElements) delete [] neighbourElements; 
    926 } 
    927  
    928 } 
     914        delete [] remapMatrix; 
     915        delete [] srcAddress; 
     916        delete [] srcRank; 
     917        delete [] dstAddress; 
     918        if (neighbourElements) delete [] neighbourElements; 
     919} 
     920 
     921} 
  • XIOS/dev/branch_openmp/extern/remap/src/mapper.hpp

    r1134 r1328  
    33#include "parallel_tree.hpp" 
    44#include "mpi.hpp" 
    5  
    6 #ifdef _usingEP 
    7 #include "ep_declaration.hpp" 
    8 #endif 
    95 
    106namespace sphereRemap { 
     
    2218{ 
    2319public: 
    24        Mapper(ep_lib::MPI_Comm comm=MPI_COMM_WORLD) : communicator(comm), verbose(SILENT), neighbourElements(NULL), sstree(comm) {} 
     20       Mapper(ep_lib::MPI_Comm comm) : communicator(comm), verbose(SILENT), neighbourElements(NULL), sstree(comm) {} 
    2521       ~Mapper(); 
    2622       void setVerbosity(verbosity v) {verbose=v ;} 
  • XIOS/dev/branch_openmp/extern/remap/src/mpi_cascade.cpp

    r688 r1328  
    11#include "mpi_cascade.hpp" 
    22#include <iostream> 
     3using namespace ep_lib; 
    34 
    45namespace sphereRemap { 
  • XIOS/dev/branch_openmp/extern/remap/src/mpi_cascade.hpp

    r694 r1328  
    1212{ 
    1313public: 
    14         CCascadeLevel(MPI_Comm comm) : comm(comm) 
    15         { 
    16                 MPI_Comm_size(comm, &size); 
    17                 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; } 
     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; } 
    2121 
    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; } 
     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; } 
    2525 
    26         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; 
     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; 
    3131}; 
    3232 
     
    3434{ 
    3535public: 
    36         //  
    37         CMPICascade(int nodes_per_level, MPI_Comm comm); 
     36  CMPICascade(int nodes_per_level, ep_lib::MPI_Comm comm); 
    3837 
    39         int num_levels; 
    40         std::vector<CCascadeLevel> level; 
     38  int num_levels; 
     39  std::vector<CCascadeLevel> level; 
    4140}; 
    4241 
  • XIOS/dev/branch_openmp/extern/remap/src/mpi_routing.cpp

    r1289 r1328  
    55#include "timerRemap.hpp" 
    66#include <iostream> 
    7 #ifdef _usingEP 
    8 #include "ep_declaration.hpp" 
    9 #endif 
     7using namespace ep_lib; 
    108 
    119namespace sphereRemap { 
     
    125123        CTimer::get("CMPIRouting::init(reduce_scatter)").print(); 
    126124 
    127         MPI_Info info_null; 
    128  
    129         // MPI_Alloc_mem(nbTarget *sizeof(int), info_null, &targetRank); 
    130         // MPI_Alloc_mem(nbSource *sizeof(int), info_null, &sourceRank); 
    131125        MPI_Alloc_mem(nbTarget *sizeof(int), MPI_INFO_NULL, &targetRank); 
    132126        MPI_Alloc_mem(nbSource *sizeof(int), MPI_INFO_NULL, &sourceRank); 
     
    158152        { 
    159153                #ifdef _usingEP 
    160                 MPI_Irecv(&sourceRank[i], 1, MPI_INT, -1, 0, communicator, &request[indexRequest]); 
     154                MPI_Irecv(&sourceRank[i], 1, MPI_INT, -2, 0, communicator, &request[indexRequest]); 
    161155                #else 
    162156                MPI_Irecv(&sourceRank[i], 1, MPI_INT, MPI_ANY_SOURCE, 0, communicator, &request[indexRequest]); 
     
    182176        { 
    183177                #ifdef _usingEP 
    184                 MPI_Irecv(&sourceRank[i], 1, MPI_INT, -1, 0, communicator, &request[indexRequest]); 
     178                MPI_Irecv(&sourceRank[i], 1, MPI_INT, -2, 0, communicator, &request[indexRequest]); 
    185179                #else 
    186180                MPI_Irecv(&sourceRank[i], 1, MPI_INT, MPI_ANY_SOURCE, 0, communicator, &request[indexRequest]); 
  • XIOS/dev/branch_openmp/extern/remap/src/mpi_routing.hpp

    r694 r1328  
    1111{ 
    1212 
    13         MPI_Comm communicator; 
     13        ep_lib::MPI_Comm communicator; 
    1414        int mpiRank; 
    1515        int mpiSize; 
     
    2929 
    3030public: 
    31         CMPIRouting(MPI_Comm comm); 
     31        CMPIRouting(ep_lib::MPI_Comm comm); 
    3232        ~CMPIRouting(); 
    3333        template<typename T> void init(const std::vector<T>& route, CMPICascade *cascade = NULL); 
     
    4444template <typename T> 
    4545void alltoalls_known(const std::vector<std::vector<T> >& send, std::vector<std::vector<T> >& recv, 
    46                      const std::vector<int>& ranks, MPI_Comm communicator); 
     46                     const std::vector<int>& ranks, ep_lib::MPI_Comm communicator); 
    4747 
    4848template <typename T> 
    4949void alltoalls_unknown(const std::vector<std::vector<T> >& send, std::vector<std::vector<T> >& recv, 
    50                        const std::vector<int>& ranks, MPI_Comm communicator); 
     50                       const std::vector<int>& ranks, ep_lib::MPI_Comm communicator); 
    5151} 
    5252#endif 
  • XIOS/dev/branch_openmp/extern/remap/src/node.cpp

    r1205 r1328  
    254254NodePtr insert(NodePtr thIs, NodePtr node) 
    255255{ 
    256   int la = thIs->level; // node to be inserted 
    257   int lb = node->level; // node where insertation 
    258   assert(la < lb); // node to be inserted must have lower level then parent 
    259   //if (thIs->parent) assert(find_in_tree1(thIs) == true); 
    260   NodePtr q = NULL; 
    261   NodePtr chd = NULL; 
    262   node->move(thIs); 
    263   if (la == lb - 1) 
    264   { 
     256        int la = thIs->level; // node to be inserted 
     257        int lb = node->level; // node where insertation 
     258        assert(la < lb); // node to be inserted must have lower level then parent 
     259        //if (thIs->parent) assert(find_in_tree1(thIs) == true); 
     260        NodePtr q = NULL; 
     261        NodePtr chd = NULL; 
     262        node->move(thIs); 
     263        if (la == lb - 1) 
     264        { 
    265265    node->child.push_back(thIs); 
    266     thIs->parent = node; 
    267     if (node->child.size() > MAX_NODE_SZ &&  node->tree->canSplit() ) // with us as additional child `node` is now too large :( 
    268     return (node->reinserted || node->parent == NULL) ? split(node) : reinsert(node); 
    269   } 
    270   else // la < lb - 1 
    271   { 
    272     chd = thIs->closest(node->child); 
    273     q = insert(thIs, chd); 
    274   } 
    275   if ((node->updateCount + 1) % UPDATE_EVERY == 0) 
    276     node->update(); 
    277   else 
    278   { 
    279     if (q) node->remove(q); 
    280     node->inflate(chd); 
    281   } 
     266                thIs->parent = node; 
     267                if (node->child.size() > MAX_NODE_SZ &&  node->tree->canSplit() ) // with us as additional child `node` is now too large :( 
     268                        return (node->reinserted || node->parent == NULL) ? split(node) : reinsert(node); 
     269        } 
     270        else // la < lb - 1 
     271        { 
     272                chd = thIs->closest(node->child); 
     273                q = insert(thIs, chd); 
     274        } 
     275        if ((node->updateCount + 1) % UPDATE_EVERY == 0) 
     276                node->update(); 
     277        else 
     278        { 
     279                if (q) node->remove(q); 
     280                node->inflate(chd); 
     281        } 
    282282 
    283283  return q; 
  • XIOS/dev/branch_openmp/extern/remap/src/node.hpp

    r1153 r1328  
    1515struct Circle 
    1616{ 
    17   Coord centre; 
    18   double radius; 
     17        Coord centre; 
     18        double radius; 
    1919}; 
    2020 
     
    116116struct Node 
    117117{ 
    118   int level; /* FIXME leafs are 0 and root is max level? */ 
    119   int leafCount; /* number of leafs that are descendants of this node (the elements in it's cycle) */ 
    120   Coord centre; 
    121   double radius; 
    122   NodePtr parent, ref; 
    123   std::vector<NodePtr> child; 
    124   std::list<NodePtr> intersectors; 
    125   bool reinserted; 
    126   int updateCount;  // double var; 
    127   CBasicTree* tree; 
    128   void *data; 
    129   int route; 
     118        int level; /* FIXME leafs are 0 and root is max level? */ 
     119        int leafCount; /* number of leafs that are descendants of this node (the elements in it's cycle) */ 
     120        Coord centre; 
     121        double radius; 
     122        NodePtr parent, ref; 
     123        std::vector<NodePtr> child; 
     124        std::list<NodePtr> intersectors; 
     125        bool reinserted; 
     126        int updateCount;  // double var; 
     127        CBasicTree* tree; 
     128        void *data; 
     129        int route; 
    130130  bool toDelete ; 
    131131 
    132   Node() : level(0), leafCount(1), centre(ORIGIN), radius(0), reinserted(false), updateCount(0), toDelete(false) {} 
    133   Node(const Coord& centre, double radius, void *data) 
    134     : level(0), leafCount(1), centre(centre), radius(radius), reinserted(false), updateCount(0), data(data), toDelete(false) {} 
     132        Node() : level(0), leafCount(1), centre(ORIGIN), radius(0), reinserted(false), updateCount(0), toDelete(false) {} 
     133        Node(const Coord& centre, double radius, void *data) 
     134                : level(0), leafCount(1), centre(centre), radius(radius), reinserted(false), updateCount(0), data(data), toDelete(false) {} 
    135135 
    136136//#ifdef DEBUG 
     
    178178//#endif 
    179179 
    180   void move(const NodePtr node); 
    181   void remove(const NodePtr node); 
    182   void inflate(const NodePtr node); 
    183   void update(); 
     180        void move(const NodePtr node); 
     181        void remove(const NodePtr node); 
     182        void inflate(const NodePtr node); 
     183        void update(); 
    184184  void output(std::ostream& flux, int level, int color) ; 
    185   NodePtr closest(std::vector<NodePtr>& list, int n = CLOSEST); 
    186   NodePtr farthest(std::vector<NodePtr>& list); 
    187   void findClosest(int level, NodePtr src, double& minDist, NodePtr &closest); 
    188  
    189   void search(NodePtr node); 
    190   bool centreInside(Node &node); 
    191   bool intersects(NodePtr node); 
    192   bool isInside(Node &node); 
    193   int incluCheck(); 
     185        NodePtr closest(std::vector<NodePtr>& list, int n = CLOSEST); 
     186        NodePtr farthest(std::vector<NodePtr>& list); 
     187        void findClosest(int level, NodePtr src, double& minDist, NodePtr &closest); 
     188 
     189        void search(NodePtr node); 
     190        bool centreInside(Node &node); 
     191        bool intersects(NodePtr node); 
     192        bool isInside(Node &node); 
     193        int incluCheck(); 
    194194  void checkParent(void) ; 
    195   void printChildren(); 
    196   void assignRoute(std::vector<int>::iterator& rank, int level); 
    197   void assignCircleAndPropagateUp(Coord *centres, double *radia, int level); 
    198   void printLevel(int level); 
    199   void routeNode(NodePtr node, int level); 
    200   void routingIntersecting(std::vector<Node>* routingList, NodePtr node); 
    201   void routeIntersection(std::vector<int>& routes, NodePtr node); 
     195        void printChildren(); 
     196        void assignRoute(std::vector<int>::iterator& rank, int level); 
     197        void assignCircleAndPropagateUp(Coord *centres, double *radia, int level); 
     198        void printLevel(int level); 
     199        void routeNode(NodePtr node, int level); 
     200        void routingIntersecting(std::vector<Node>* routingList, NodePtr node); 
     201        void routeIntersection(std::vector<int>& routes, NodePtr node); 
    202202  void getNodeLevel(int level,std::list<NodePtr>& NodeList) ; 
    203203  bool removeDeletedNodes(int assignLevel) ; 
  • XIOS/dev/branch_openmp/extern/remap/src/parallel_tree.cpp

    r1295 r1328  
    1212 
    1313#include "parallel_tree.hpp" 
     14using namespace ep_lib; 
    1415 
    1516namespace sphereRemap { 
    16  
    17 extern CRemapGrid srcGrid; 
    18 #pragma omp threadprivate(srcGrid) 
    19  
    20 extern CRemapGrid tgtGrid; 
    21 #pragma omp threadprivate(tgtGrid) 
    2217 
    2318static const int assignLevel = 2; 
     
    292287{ 
    293288 
    294   int assignLevel = 2; 
    295   int nbSampleNodes = 2*ipow(MAX_NODE_SZ + 1, assignLevel); 
     289        int assignLevel = 2; 
     290        int nbSampleNodes = 2*ipow(MAX_NODE_SZ + 1, assignLevel); 
    296291 
    297292 
     
    305300  MPI_Comm_size(communicator,&commSize) ; 
    306301   
    307   // make multiple of two 
    308   nbSampleNodes /= 2; 
    309   nbSampleNodes *= 2; 
    310   //assert( nbTot > nbSampleNodes*commSize) ; 
     302        // make multiple of two 
     303        nbSampleNodes /= 2; 
     304        nbSampleNodes *= 2; 
     305//  assert( nbTot > nbSampleNodes*commSize) ; 
    311306     
    312307  int nbSampleNodes1 = nbSampleNodes * (nb1*commSize)/(1.*nbTot) ; 
     
    314309   
    315310 
    316   //assert(node.size() > nbSampleNodes); 
    317   //assert(node2.size() > nbSampleNodes); 
    318   //assert(node.size() + node2.size() > nbSampleNodes); 
    319   vector<Node> sampleNodes; sampleNodes.reserve(nbSampleNodes1+nbSampleNodes2); 
    320  
    321   vector<int> randomArray1(node.size()); 
    322   randomizeArray(randomArray1); 
    323   vector<int> randomArray2(node2.size()); 
    324   randomizeArray(randomArray2); 
    325  
    326   for (int i = 0; i <nbSampleNodes1; i++) sampleNodes.push_back(Node(node[randomArray1[i%nb1]].centre,  node[randomArray1[i%nb1]].radius, NULL)); 
    327   for (int i = 0; i <nbSampleNodes2; i++) sampleNodes.push_back(Node(node2[randomArray2[i%nb2]].centre, node2[randomArray2[i%nb2]].radius, NULL)); 
     311//      assert(node.size() > nbSampleNodes); 
     312//      assert(node2.size() > nbSampleNodes); 
     313//      assert(node.size() + node2.size() > nbSampleNodes); 
     314        vector<Node> sampleNodes; sampleNodes.reserve(nbSampleNodes1+nbSampleNodes2); 
     315 
     316        vector<int> randomArray1(node.size()); 
     317        randomizeArray(randomArray1); 
     318        vector<int> randomArray2(node2.size()); 
     319        randomizeArray(randomArray2); 
     320 
     321        for (int i = 0; i <nbSampleNodes1; i++) sampleNodes.push_back(Node(node[randomArray1[i%nb1]].centre,  node[randomArray1[i%nb1]].radius, NULL)); 
     322        for (int i = 0; i <nbSampleNodes2; i++) sampleNodes.push_back(Node(node2[randomArray2[i%nb2]].centre, node2[randomArray2[i%nb2]].radius, NULL)); 
    328323 
    329324        CTimer::get("buildParallelSampleTree").resume(); 
     
    336331        CTimer::get("parallelRouteNode").resume(); 
    337332        vector<int> route(node.size()); 
     333        cout<<"node.size = "<<node.size()<<endl; 
    338334        routeNodes(route /*out*/, node); 
    339335        CTimer::get("parallelRouteNode").suspend(); 
  • XIOS/dev/branch_openmp/extern/remap/src/parallel_tree.hpp

    r1134 r1328  
    66#include "mpi_cascade.hpp" 
    77#include "mpi.hpp" 
    8 #ifdef _usingEP 
    9 #include "ep_declaration.hpp" 
    10 #endif 
    11  
    128namespace sphereRemap { 
    139 
     
    1511{ 
    1612public: 
    17         CParallelTree(ep_lib::MPI_Comm comm); 
    18         ~CParallelTree(); 
     13  CParallelTree(ep_lib::MPI_Comm comm); 
     14  ~CParallelTree(); 
    1915 
    20         void build(vector<Node>& node, vector<Node>& node2); 
     16  void build(vector<Node>& node, vector<Node>& node2); 
    2117 
    22         void routeNodes(vector<int>& route, vector<Node>& nodes, int level = 0); 
    23         void routeIntersections(vector<vector<int> >& route, vector<Node>& nodes, int level = 0); 
     18  void routeNodes(vector<int>& route, vector<Node>& nodes, int level = 0); 
     19  void routeIntersections(vector<vector<int> >& route, vector<Node>& nodes, int level = 0); 
    2420 
    25         int nbLocalElements; 
    26         Elt* localElements; 
     21  int nbLocalElements; 
     22  Elt* localElements; 
    2723 
    28         CTree localTree; 
     24  CTree localTree; 
    2925 
    3026private: 
    31         void updateCirclesForRouting(Coord rootCentre, double rootRadius, int level = 0); 
    32         void buildSampleTreeCascade(vector<Node>& sampleNodes, int level = 0); 
    33         void buildLocalTree(const vector<Node>& node, const vector<int>& route); 
    34         void buildRouteTree(); 
     27  void updateCirclesForRouting(Coord rootCentre, double rootRadius, int level = 0); 
     28  void buildSampleTreeCascade(vector<Node>& sampleNodes, int level = 0); 
     29  void buildLocalTree(const vector<Node>& node, const vector<int>& route); 
     30  void buildRouteTree(); 
    3531 
    36         //CSampleTree sampleTree; 
    37         vector<CSampleTree> treeCascade; // first for sample tree, then for routing tree 
    38         CMPICascade cascade; 
    39         ep_lib::MPI_Comm communicator ; 
     32  //CSampleTree sampleTree; 
     33  vector<CSampleTree> treeCascade; // first for sample tree, then for routing tree 
     34  CMPICascade cascade; 
     35  ep_lib::MPI_Comm communicator ; 
    4036 
    4137}; 
  • XIOS/dev/branch_openmp/extern/remap/src/polyg.cpp

    r1289 r1328  
    77 
    88#include "polyg.hpp" 
    9 #include <stdio.h> 
    109 
    1110namespace sphereRemap { 
     
    4645Coord barycentre(const Coord *x, int n) 
    4746{ 
    48  
    4947        if (n == 0) return ORIGIN; 
    5048        Coord bc = ORIGIN; 
    5149        for (int i = 0; i < n; i++) 
    52         { 
    5350                bc = bc + x[i]; 
    54         }        
    5551        /* both distances can be equal down to roundoff when norm(bc) < mashineepsilon  
    5652           which can occur when weighted with tiny area */ 
    57    
    58          
    59   //assert(squaredist(bc, proj(bc)) <= squaredist(bc, proj(bc * (-1.0))));       
     53 
     54        assert(squaredist(bc, proj(bc)) <= squaredist(bc, proj(bc * (-1.0)))); 
     55        //if (squaredist(bc, proj(bc)) > squaredist(bc, proj(bc * (-1.0)))) return proj(bc * (-1.0)); 
    6056 
    6157        return proj(bc); 
     
    165161{ 
    166162        if (N < 3) 
    167                 return 0; /* polygons with less than three vertices have zero area */ 
     163          return 0; /* polygons with less than three vertices have zero area */ 
    168164        Coord t[3]; 
    169165        t[0] = barycentre(x, N); 
     
    177173                t[1] = x[i]; 
    178174                t[2] = x[ii]; 
    179     double sc=scalarprod(crossprod(t[1] - t[0], t[2] - t[0]), t[0]) ; 
     175                double sc=scalarprod(crossprod(t[1] - t[0], t[2] - t[0]), t[0]) ; 
    180176                assert(sc >= -1e-10); // Error: tri a l'env (wrong orientation) 
    181177                double area_gc = triarea(t[0], t[1], t[2]); 
    182                 //if(area_gc<=0) printf("area_gc = %e\n", area_gc); 
    183178                double area_sc_gc_moon = 0; 
    184179                if (d[i]) /* handle small circle case */ 
     
    188183                        char sgl = (mext > 0) ? -1 : 1; 
    189184                        area_sc_gc_moon = sgl * alun(arcdist(t[1], t[2]), fabs(scalarprod(t[1], pole))); 
    190                         //if(area_sc_gc_moon<=0) printf("area_sc_gc_moon = %e\n", area_sc_gc_moon); 
    191185                        gg_exact = gg_exact + sc_gc_moon_normalintegral(t[1], t[2], pole); 
    192186                } 
    193187                area += area_gc + area_sc_gc_moon; /* for "spherical circle segment" sum triangular part (at) and "small moon" and => account for small circle */ 
    194188                g[i] = barycentre(t, 3) * (area_gc + area_sc_gc_moon); 
    195                 //printf("g[%d] = (%e,%e,%e) * (%e+%e) = (%e,%e,%e) norm = %e\n", i, barycentre(t, 3).x, barycentre(t, 3).y, barycentre(t, 3).z, area_gc,  area_sc_gc_moon, g[i].x, g[i].y, g[i].z, norm(g[i])); 
    196189        } 
    197190        gg = barycentre(g, N); 
  • XIOS/dev/branch_openmp/extern/remap/src/timerRemap.cpp

    r1146 r1328  
    44#include <map> 
    55#include <iostream> 
     6using namespace ep_lib; 
    67 
    78namespace sphereRemap { 
     
    910using namespace std; 
    1011 
    11 map<string,CTimer*> *CTimer::allTimer = 0; 
     12//map<string,CTimer*> CTimer::allTimer; 
     13map<string,CTimer*> *CTimer::allTimer_ptr = 0; 
    1214 
    1315CTimer::CTimer(const string& name_) : name(name_) 
     
    5557CTimer& CTimer::get(const string name) 
    5658{ 
    57         if(allTimer == 0) allTimer = new map<string,CTimer*>; 
    5859        map<string,CTimer*>::iterator it; 
    59         it=(*allTimer).find(name); 
    60         if (it==(*allTimer).end()) it=(*allTimer).insert(pair<string,CTimer*>(name,new CTimer(name))).first; 
     60        if(allTimer_ptr == 0) allTimer_ptr = new map<string,CTimer*>; 
     61        //it=allTimer.find(name); 
     62        it=allTimer_ptr->find(name); 
     63        //if (it==allTimer.end()) it=allTimer.insert(pair<string,CTimer*>(name,new CTimer(name))).first; 
     64        if (it==allTimer_ptr->end()) it=allTimer_ptr->insert(pair<string,CTimer*>(name,new CTimer(name))).first; 
    6165        return *(it->second); 
    6266} 
  • XIOS/dev/branch_openmp/extern/remap/src/timerRemap.hpp

    r1146 r1328  
    2727    void print(void); 
    2828    //static map<string,CTimer*> allTimer; 
    29     static map<string,CTimer*> *allTimer; 
    30     #pragma omp threadprivate(allTimer) 
    31      
     29    static map<string,CTimer*> *allTimer_ptr; 
    3230    static double getTime(void); 
    3331    static CTimer& get(string name); 
  • XIOS/dev/branch_openmp/extern/remap/src/tree.cpp

    r1172 r1328  
    2929void CBasicTree::routeNodes(vector<int>& route, vector<Node>& nodes, int assignLevel) 
    3030{ 
    31   for (int i = 0; i < nodes.size(); i++) 
    32   { 
    33     root->routeNode(&nodes[i], assignLevel); 
    34     route[i] = nodes[i].route; 
    35   } 
     31        for (int i = 0; i < nodes.size(); i++) 
     32        { 
     33                root->routeNode(&nodes[i], assignLevel); 
     34                route[i] = nodes[i].route; 
     35        } 
    3636} 
    3737 
    3838void CBasicTree::routeIntersections(vector<vector<int> >& routes, vector<Node>& nodes) 
    3939{ 
    40   for (int i = 0; i < nodes.size(); i++) 
    41     root->routeIntersection(routes[i], &nodes[i]); 
     40        for (int i = 0; i < nodes.size(); i++) 
     41                root->routeIntersection(routes[i], &nodes[i]); 
    4242} 
    4343 
    4444void CBasicTree::build(vector<Node>& nodes) 
    4545{ 
    46   newRoot(1); 
    47   insertNodes(nodes); 
     46        newRoot(1); 
     47        insertNodes(nodes); 
    4848} 
    4949 
    5050void CBasicTree::output(ostream& flux, int level) 
    5151{ 
    52   root->output(flux,level,0) ; 
     52        root->output(flux,level,0) ; 
    5353} 
    5454void CBasicTree::slim(int nbIts) 
    5555{ 
    56   for (int i = 0; i < nbIts; i++) 
    57   { 
    58     for (int level = root->level - 1; level > 0; level--) 
    59     { 
    60       slim2(root, level); 
    61       ri = 0; 
    62       emptyPool(); 
    63     } 
    64  
    65     for (int level = 2; level < root->level; level++) 
    66     { 
    67       slim2(root, level); 
    68       ri = 0; 
    69       emptyPool(); 
    70     } 
    71   } 
     56        for (int i = 0; i < nbIts; i++) 
     57        { 
     58                for (int level = root->level - 1; level > 0; level--) 
     59                { 
     60                        slim2(root, level); 
     61                        ri = 0; 
     62                        emptyPool(); 
     63                } 
     64 
     65                for (int level = 2; level < root->level; level++) 
     66                { 
     67                        slim2(root, level); 
     68                        ri = 0; 
     69                        emptyPool(); 
     70                } 
     71        } 
    7272} 
    7373 
     
    7676void CBasicTree::insertNode(NodePtr node) 
    7777{ 
    78   node->tree = this; 
    79   increaseLevelSize(0); 
    80   push_back(node); 
    81  
    82   NodePtr q; 
    83   while (pool.size()) 
    84   { 
    85     q = pool.front(); 
    86     pool.pop_front(); 
    87     q = insert(q, root); 
    88     if (ri) 
    89     { 
    90       delete q; 
    91       ri = 0; 
    92     } 
    93   } 
     78        node->tree = this; 
     79        increaseLevelSize(0); 
     80        push_back(node); 
     81 
     82        NodePtr q; 
     83        while (pool.size()) 
     84        { 
     85                q = pool.front(); 
     86                pool.pop_front(); 
     87                q = insert(q, root); 
     88                if (ri) 
     89                { 
     90                        delete q; 
     91                        ri = 0; 
     92                } 
     93        } 
    9494} 
    9595 
    9696void CBasicTree::emptyPool(void) 
    9797{ 
    98   while (pool.size()) 
    99   { 
    100     NodePtr q = pool.front(); 
    101     pool.pop_front(); 
    102     q = insert(q, root); 
    103     if (ri) 
    104     { 
    105       delete q; 
    106       ri = 0; 
    107     } 
    108   } 
     98        while (pool.size()) 
     99        { 
     100                NodePtr q = pool.front(); 
     101                pool.pop_front(); 
     102                q = insert(q, root); 
     103                if (ri) 
     104                { 
     105                        delete q; 
     106                        ri = 0; 
     107                } 
     108        } 
    109109} 
    110110 
     
    142142        root->parent = 0; 
    143143        root->leafCount = 0; 
    144         // initialize root node on the sphere 
    145         root->centre.x=1 ;  
    146         root->centre.y=0 ;  
    147         root->centre.z=0 ;  
     144// initialize root node on the sphere 
     145  root->centre.x=1 ; root->centre.y=0 ; root->centre.z=0 ;  
    148146        root->radius = 0.; 
    149147        root->reinserted = false; 
  • XIOS/dev/branch_openmp/extern/remap/src/tree.hpp

    r1172 r1328  
    1414class CBasicTree 
    1515{ 
    16   public: 
     16public: 
    1717 
    18   NodePtr root; /* The main tree is stored as Nodes which can be reached through traversal starting here */ 
    19   NodePtr ref; // FIXME this reference, set by a node is odd, try to remove 
    20   int ri; /** this is set to one by a node in case of reinsertion */ 
    21   vector<int> levelSize; /** e.g. levelSize[0] == leafs.size() */ 
    22   vector<Node> leafs; /** leafs are stored in vector for easy access and rest of the tree nodes as separate allocations, only reachable through tree traversal */ 
     18        NodePtr root; /* The main tree is stored as Nodes which can be reached through traversal starting here */ 
     19        NodePtr ref; // FIXME this reference, set by a node is odd, try to remove 
     20        int ri; /** this is set to one by a node in case of reinsertion */ 
     21        vector<int> levelSize; /** e.g. levelSize[0] == leafs.size() */ 
     22        vector<Node> leafs; /** leafs are stored in vector for easy access and rest of the tree nodes as separate allocations, only reachable through tree traversal */ 
    2323 
    24   CBasicTree() : ri(0), levelSize(MAX_LEVEL_SIZE), root(NULL), isAssignedLevel(false), okSplit(true), isActiveOkSplit(false) {}  
    25   ~CBasicTree();  
    26   void build(vector<Node>& nodes); 
    27   void slim(int nbIts = 1); 
    28   virtual void insertNodes(vector<Node>& node) = 0; 
     24        CBasicTree() : ri(0), levelSize(MAX_LEVEL_SIZE), root(NULL), isAssignedLevel(false), okSplit(true), isActiveOkSplit(false) {}  
     25        ~CBasicTree();  
     26        void build(vector<Node>& nodes); 
     27        void slim(int nbIts = 1); 
     28        virtual void insertNodes(vector<Node>& node) = 0; 
    2929 
    30   void routeNodes(vector<int>& route, vector<Node>& nodes, int assignLevel); 
    31   void routeIntersections(vector<vector<int> >& route, vector<Node>& nodes); 
     30        void routeNodes(vector<int>& route, vector<Node>& nodes, int assignLevel); 
     31        void routeIntersections(vector<vector<int> >& route, vector<Node>& nodes); 
    3232 
    33   void push_back(NodePtr node); 
    34   void push_front(NodePtr node); 
    35   void increaseLevelSize(int level); 
    36   void decreaseLevelSize(int level); 
    37   void newRoot(int level); 
    38   void insertNode(NodePtr node); 
     33        void push_back(NodePtr node); 
     34        void push_front(NodePtr node); 
     35        void increaseLevelSize(int level); 
     36        void decreaseLevelSize(int level); 
     37        void newRoot(int level); 
     38        void insertNode(NodePtr node); 
    3939  void output(ostream& flux, int level) ; 
    4040 
    41   int keepNodes; 
     41        int keepNodes; 
    4242  bool isAssignedLevel ;  
    4343  int assignLevel; 
     
    5050 
    5151   
    52   private: 
    53   deque<NodePtr > pool; 
     52private: 
     53        deque<NodePtr > pool; 
    5454         
    5555  bool okSplit ; 
    5656  
    57   protected: 
     57protected: 
    5858  void emptyPool(); 
    59   CBasicTree(int keepNodes_, int assignLevel_) : ri(0), levelSize(MAX_LEVEL_SIZE), root(NULL), keepNodes(keepNodes_), assignLevel(assignLevel_), isAssignedLevel(true),  
    60                                                  okSplit(true), isActiveOkSplit(false) {}  
     59  CBasicTree(int keepNodes_, int assignLevel_) : ri(0), levelSize(MAX_LEVEL_SIZE), root(NULL), keepNodes(keepNodes_), assignLevel(assignLevel_), isAssignedLevel(true), okSplit(true), isActiveOkSplit(false) {}  
    6160}; 
    6261 
    6362class CTree : public CBasicTree 
    6463{ 
    65   public: 
    66   void insertNodes(vector<Node>& nodes); 
     64public: 
     65        void insertNodes(vector<Node>& nodes); 
    6766}; 
    6867 
     
    7069{ 
    7170 
    72   public: 
    73   CSampleTree(int keepNodes_, int assignLevel_) : CBasicTree(keepNodes_,assignLevel_) {} 
     71public: 
     72        CSampleTree(int keepNodes_, int assignLevel_) : CBasicTree(keepNodes_,assignLevel_) {} 
    7473  void slimAssignedLevel() ; 
    7574  void removeExtraNode(void) ; 
    76   void insertNodes(vector<Node>& nodes); 
     75        void insertNodes(vector<Node>& nodes); 
    7776}; 
    7877 
  • XIOS/dev/branch_openmp/extern/src_ep_dev/ep_fortran.cpp

    r1287 r1328  
    3131      { 
    3232        fc_comm_map.insert(std::make_pair( std::make_pair( fint, omp_get_thread_num()) , comm)); 
    33         //printf("EP_Comm_c2f : MAP %p insert: %d, %d, %p\n", &fc_comm_map, fint, omp_get_thread_num(), comm.ep_comm_ptr); 
     33        printf("EP_Comm_c2f : MAP %p insert: %d, %d, %p\n", &fc_comm_map, fint, omp_get_thread_num(), comm.ep_comm_ptr); 
    3434      } 
    3535    } 
     
    5454      MPI_Comm comm_ptr; 
    5555      comm_ptr = it->second; 
    56       //printf("EP_Comm_f2c : MAP %p find: %d, %d, %p\n", &fc_comm_map, it->first.first, it->first.second, comm_ptr.ep_comm_ptr); 
     56      printf("EP_Comm_f2c : MAP %p find: %d, %d, %p\n", &fc_comm_map, it->first.first, it->first.second, comm_ptr.ep_comm_ptr); 
    5757      return  comm_ptr; 
    5858    } 
  • XIOS/dev/branch_openmp/extern/src_ep_dev/ep_intercomm.cpp

    r1287 r1328  
    4747 
    4848      MPI_Waitall(2, request, status); 
     49 
     50      //MPI_Send(&leader_ranks[0], 3, static_cast< ::MPI_Datatype>(MPI_INT), remote_leader, tag, peer_comm); 
     51      //MPI_Recv(&leader_ranks[3], 3, static_cast< ::MPI_Datatype>(MPI_INT), remote_leader, tag, peer_comm, &status[1]); 
    4952    } 
    5053 
  • XIOS/dev/branch_openmp/extern/src_ep_dev/ep_lib.cpp

    r1295 r1328  
    1010#pragma omp threadprivate(EP_PendingRequests) 
    1111 
     12 
     13 
    1214namespace ep_lib 
    1315{  
     16  bool MPI_Comm::is_null() 
     17  { 
     18    if(!this->is_intercomm) 
     19      return this->mpi_comm == MPI_COMM_NULL.mpi_comm; 
     20    else 
     21      return this->ep_comm_ptr->intercomm->mpi_inter_comm == MPI_COMM_NULL.mpi_comm; 
     22  } 
    1423 
    1524  int tag_combine(int real_tag, int src, int dest) 
  • XIOS/dev/branch_openmp/extern/src_ep_dev/ep_send.cpp

    r1295 r1328  
    1919    if(!comm.is_ep) 
    2020      return ::MPI_Send(buf, count, to_mpi_type(datatype), dest, tag, to_mpi_comm(comm.mpi_comm)); 
    21  
    22     MPI_Request request; 
    23     MPI_Status status; 
    24     MPI_Isend(buf, count, datatype, dest, tag, comm, &request); 
    25     MPI_Wait(&request, &status); 
    26  
     21    if(comm.is_intercomm) 
     22    { 
     23      MPI_Request request; 
     24      MPI_Status status; 
     25      MPI_Isend(buf, count, datatype, dest, tag, comm, &request); 
     26      MPI_Wait(&request, &status); 
     27    } 
     28    else 
     29    { 
     30      int ep_src_loc  = comm.ep_comm_ptr->size_rank_info[1].first; 
     31      int ep_dest_loc = comm.ep_comm_ptr->comm_list->rank_map->at(dest).first; 
     32      int mpi_tag     = tag_combine(tag, ep_src_loc, ep_dest_loc); 
     33      int mpi_dest    = comm.ep_comm_ptr->comm_list->rank_map->at(dest).second; 
     34 
     35      ::MPI_Send(buf, count, to_mpi_type(datatype), mpi_dest, mpi_tag, to_mpi_comm(comm.mpi_comm)); 
     36      //printf("call mpi_send for intracomm, dest = %d, tag = %d\n", dest, tag); 
     37    } 
    2738    //check_sum_send(buf, count, datatype, dest, tag, comm); 
    2839 
  • XIOS/dev/branch_openmp/extern/src_ep_dev/ep_type.hpp

    r1295 r1328  
    344344    } 
    345345 
     346    bool is_null(); 
     347 
    346348  }; 
    347349 
  • XIOS/dev/branch_openmp/inputs/REMAP/iodef.xml

    r1220 r1328  
    5151 
    5252      <file_group id="write_files" > 
    53         <file id="output_2D" name="output_2D" enabled=".TRUE."> 
     53        <file id="output_2D" name="output_2D" output_freq="1ts" enabled=".TRUE."> 
    5454          <field field_ref="src_field_2D" name="field_src" enabled=".TRUE."/> 
    5555          <field field_ref="src_field_2D_clone" name="field_src_clone" default_value="100000" enabled=".TRUE."/> 
     
    5757          <field field_ref="dst_field_2D" name="field_dst_regular_1" enabled=".TRUE." /> 
    5858          <field field_ref="dst_field_2D_regular_pole" name="field_dst_regular_2" enabled=".TRUE."/> 
    59           <field field_ref="dst_field_2D_clone" name="field_dst_regular_3" detect_missing_value=".false." default_value="100000" enabled=".TRUE."/> 
     59          <field field_ref="dst_field_2D_clone" name="field_dst_regular_3" detect_missing_value=".false." default_value="100000" enabled=".FALSE."/> 
    6060          <field field_ref="dst_field_2D_extract" name="field_dst_regular_4" enabled=".TRUE."/> 
    6161        </file>  
  • XIOS/dev/branch_openmp/src/array_new.hpp

    r1134 r1328  
    554554        TinyVector<int,N_rank> vect; 
    555555        size_t ne; 
    556  
     556         
    557557        ret =  buffer.get(numDim); 
    558558        ret &= buffer.get(vect.data(), N_rank); 
  • XIOS/dev/branch_openmp/src/attribute_enum.hpp

    r1134 r1328  
    1414namespace xios 
    1515{ 
    16     /// ////////////////////// Déclarations ////////////////////// /// 
    17     /*! 
    18       \class CAttributeEnum 
    19       This class implements the attribute representing enumeration 
     16      /// ////////////////////// Déclarations ////////////////////// /// 
     17        /*! 
     18        \class CAttributeEnum 
     19        This class implements the attribute representing enumeration 
    2020      */ 
    21     template <class T> 
    22         class CAttributeEnum : public CAttribute, public CEnum<T> 
    23     { 
     21      template <class T> 
     22         class CAttributeEnum : public CAttribute, public CEnum<T> 
     23      { 
    2424        typedef typename T::t_enum T_enum ; 
    2525        public : 
    2626 
    27         /// Constructeurs /// 
    28         explicit CAttributeEnum(const StdString & id); 
    29         CAttributeEnum(const StdString & id, 
    30                 xios_map<StdString, CAttribute*> & umap); 
    31         CAttributeEnum(const StdString & id, const T_enum & value); 
    32         CAttributeEnum(const StdString & id, const T_enum & value, 
    33                 xios_map<StdString, CAttribute*> & umap); 
     27            /// Constructeurs /// 
     28            explicit CAttributeEnum(const StdString & id); 
     29            CAttributeEnum(const StdString & id, 
     30                               xios_map<StdString, CAttribute*> & umap); 
     31            CAttributeEnum(const StdString & id, const T_enum & value); 
     32            CAttributeEnum(const StdString & id, const T_enum & value, 
     33                               xios_map<StdString, CAttribute*> & umap); 
    3434 
    35         /// Accesseur /// 
    36         T_enum getValue(void) const; 
    37         string getStringValue(void) const; 
     35            /// Accesseur /// 
     36            T_enum getValue(void) const; 
     37            string getStringValue(void) const; 
    3838 
    3939 
    40         /// Mutateurs /// 
    41         void setValue(const T_enum & value); 
     40            /// Mutateurs /// 
     41            void setValue(const T_enum & value); 
     42             
     43            void set(const CAttribute& attr) ; 
     44            void set(const CAttributeEnum& attr) ; 
     45            void reset(void); 
     46             
     47            void setInheritedValue(const CAttributeEnum& attr ); 
     48            void setInheritedValue(const CAttribute& attr ); 
     49            T_enum getInheritedValue(void)  const; 
     50            string getInheritedStringValue(void) const; 
     51            bool hasInheritedValue(void) const;           
     52           
     53            bool isEqual(const CAttributeEnum& attr ); 
     54            bool isEqual(const CAttribute& attr ); 
    4255 
    43         void set(const CAttribute& attr) ; 
    44         void set(const CAttributeEnum& attr) ; 
    45         void reset(void); 
     56            /// Destructeur /// 
     57            virtual ~CAttributeEnum(void) { } 
    4658 
    47         void setInheritedValue(const CAttributeEnum& attr ); 
    48         void setInheritedValue(const CAttribute& attr ); 
    49         T_enum getInheritedValue(void)  const; 
    50         string getInheritedStringValue(void) const; 
    51         bool hasInheritedValue(void) const;           
     59            /// Operateur /// 
     60            CAttributeEnum& operator=(const T_enum & value); 
    5261 
    53         bool isEqual(const CAttributeEnum& attr ); 
    54         bool isEqual(const CAttribute& attr ); 
     62            /// Autre /// 
     63            virtual StdString toString(void) const { return _toString();} 
     64            virtual void fromString(const StdString & str) { if (str==resetInheritanceStr) { reset(); _canInherite=false ;}  else _fromString(str);} 
    5565 
    56         /// Destructeur /// 
    57         virtual ~CAttributeEnum(void) { } 
     66            virtual bool toBuffer  (CBufferOut& buffer) const { return _toBuffer(buffer);}  
     67            virtual bool fromBuffer(CBufferIn& buffer) { return _fromBuffer(buffer); }  
     68             
     69            virtual void generateCInterface(ostream& oss,const string& className) ; 
     70            virtual void generateFortran2003Interface(ostream& oss,const string& className) ; 
     71            virtual void generateFortranInterfaceDeclaration_(ostream& oss,const string& className) ; 
     72            virtual void generateFortranInterfaceBody_(ostream& oss,const string& className) ; 
     73            virtual void generateFortranInterfaceDeclaration(ostream& oss,const string& className) ; 
     74            virtual void generateFortranInterfaceGetDeclaration_(ostream& oss,const string& className) ; 
     75            virtual void generateFortranInterfaceGetBody_(ostream& oss,const string& className) ; 
     76            virtual void generateFortranInterfaceGetDeclaration(ostream& oss,const string& className) ;       
    5877 
    59         /// Operateur /// 
    60         CAttributeEnum& operator=(const T_enum & value); 
    61  
    62         /// Autre /// 
    63         virtual StdString toString(void) const { return _toString();} 
    64         virtual void fromString(const StdString & str) { if (str==resetInheritanceStr) { reset(); _canInherite=false ;}  else _fromString(str);} 
    65  
    66         virtual bool toBuffer  (CBufferOut& buffer) const { return _toBuffer(buffer);}  
    67         virtual bool fromBuffer(CBufferIn& buffer) { return _fromBuffer(buffer); }  
    68  
    69         virtual void generateCInterface(ostream& oss,const string& className) ; 
    70         virtual void generateFortran2003Interface(ostream& oss,const string& className) ; 
    71         virtual void generateFortranInterfaceDeclaration_(ostream& oss,const string& className) ; 
    72         virtual void generateFortranInterfaceBody_(ostream& oss,const string& className) ; 
    73         virtual void generateFortranInterfaceDeclaration(ostream& oss,const string& className) ; 
    74         virtual void generateFortranInterfaceGetDeclaration_(ostream& oss,const string& className) ; 
    75         virtual void generateFortranInterfaceGetBody_(ostream& oss,const string& className) ; 
    76         virtual void generateFortranInterfaceGetDeclaration(ostream& oss,const string& className) ;       
    77  
    78         private : 
    79         StdString _toString(void) const; 
    80         void _fromString(const StdString & str); 
    81         bool _toBuffer  (CBufferOut& buffer) const; 
    82         bool _fromBuffer(CBufferIn& buffer) ; 
    83         CEnum<T> inheritedValue ; 
    84     }; // class CAttributeEnum     
    85  
     78         private : 
     79          StdString _toString(void) const; 
     80          void _fromString(const StdString & str); 
     81          bool _toBuffer  (CBufferOut& buffer) const; 
     82          bool _fromBuffer(CBufferIn& buffer) ; 
     83          CEnum<T> inheritedValue ; 
     84      }; // class CAttributeEnum     
     85    
    8686} // namespace xios 
    8787 
    8888#endif // __XIOS_ATTRIBUTE_ENUM__ 
    89  
  • XIOS/dev/branch_openmp/src/attribute_enum_impl.hpp

    r1134 r1328  
    1010namespace xios 
    1111{ 
    12   /// ////////////////////// Définitions ////////////////////// /// 
     12  /// ////////////////////// Définitions ////////////////////// /// 
    1313  template <class T> 
    1414  CAttributeEnum<T>::CAttributeEnum(const StdString & id) 
     
    3030     umap.insert(umap.end(), std::make_pair(id, this)); 
    3131  } 
    32   
     32 
    3333  template <class T> 
    3434  CAttributeEnum<T>::CAttributeEnum 
     
    4040     umap.insert(umap.end(), std::make_pair(id, this)); 
    4141  } 
    42   
     42 
    4343  ///-------------------------------------------------------------- 
    4444  template <class T> 
     
    5454     return CEnum<T>::get(); 
    5555  } 
    56   
     56 
    5757  template <class T> 
    5858  string CAttributeEnum<T>::getStringValue(void) const 
    5959  { 
    60     return CEnum<T>::toString(); 
    61   } 
    62  
     60     return CEnum<T>::toString(); 
     61  } 
    6362 
    6463  template <class T> 
     
    7170  void CAttributeEnum<T>::set(const CAttribute& attr) 
    7271  { 
    73      this->set(dynamic_cast<const CAttributeEnum<T>& >(attr)); 
    74   } 
    75     
    76   template <class T> 
     72    this->set(dynamic_cast<const CAttributeEnum<T>& >(attr)); 
     73  } 
     74 
     75 template <class T> 
    7776  void CAttributeEnum<T>::set(const CAttributeEnum& attr) 
    7877  { 
    79      CEnum<T>::set(attr); 
    80   } 
    81   
     78    CEnum<T>::set(attr); 
     79  } 
     80 
    8281  template <class T> 
    8382  void CAttributeEnum<T>::setInheritedValue(const CAttribute& attr) 
    8483  { 
    85      this->setInheritedValue(dynamic_cast<const CAttributeEnum<T>& >(attr)); 
    86   } 
    87    
     84    this->setInheritedValue(dynamic_cast<const CAttributeEnum<T>& >(attr)); 
     85  } 
     86 
    8887  template <class T> 
    8988  void CAttributeEnum<T>::setInheritedValue(const CAttributeEnum& attr) 
    9089  { 
    91      if (this->isEmpty() && _canInherite && attr.hasInheritedValue()) inheritedValue.set(attr.getInheritedValue()); 
    92   } 
    93    
     90    if (this->isEmpty() && _canInherite && attr.hasInheritedValue()) inheritedValue.set(attr.getInheritedValue()); 
     91  } 
     92 
    9493  template <class T> 
    9594  typename T::t_enum CAttributeEnum<T>::getInheritedValue(void) const 
    9695  { 
    97      if (this->isEmpty()) return inheritedValue.get(); 
    98      else return getValue(); 
    99   } 
    100  
    101   template <class T> 
    102       string CAttributeEnum<T>::getInheritedStringValue(void) const 
    103       { 
    104           if (this->isEmpty()) return inheritedValue.toString(); 
    105           else return CEnum<T>::toString();; 
    106       } 
    107  
    108   template <class T> 
    109       bool CAttributeEnum<T>::hasInheritedValue(void) const 
    110       { 
    111           return !this->isEmpty() || !inheritedValue.isEmpty(); 
    112       } 
    113  
    114   template <class T> 
    115       bool CAttributeEnum<T>::isEqual(const CAttribute& attr) 
    116       { 
    117           return (this->isEqual(dynamic_cast<const CAttributeEnum<T>& >(attr))); 
    118       } 
    119  
    120   template <class T> 
    121       bool CAttributeEnum<T>::isEqual(const CAttributeEnum& attr) 
    122       { 
    123           return ((dynamic_cast<const CEnum<T>& >(*this)) == (dynamic_cast<const CEnum<T>& >(attr))); 
    124       } 
     96    if (this->isEmpty()) return inheritedValue.get(); 
     97    else return getValue(); 
     98  } 
     99 
     100  template <class T> 
     101  string CAttributeEnum<T>::getInheritedStringValue(void) const 
     102  { 
     103     if (this->isEmpty()) return inheritedValue.toString(); 
     104     else return CEnum<T>::toString();; 
     105  } 
     106 
     107  template <class T> 
     108  bool CAttributeEnum<T>::hasInheritedValue(void) const 
     109  { 
     110    return !this->isEmpty() || !inheritedValue.isEmpty(); 
     111  } 
     112 
     113  template <class T> 
     114  bool CAttributeEnum<T>::isEqual(const CAttribute& attr) 
     115  { 
     116    return (this->isEqual(dynamic_cast<const CAttributeEnum<T>& >(attr))); 
     117  } 
     118 
     119  template <class T> 
     120  bool CAttributeEnum<T>::isEqual(const CAttributeEnum& attr) 
     121  { 
     122    return ((dynamic_cast<const CEnum<T>& >(*this)) == (dynamic_cast<const CEnum<T>& >(attr))); 
     123  } 
    125124 
    126125  //--------------------------------------------------------------- 
    127126 
    128127  template <class T> 
    129       CAttributeEnum<T>& CAttributeEnum<T>::operator=(const T_enum & value) 
    130       { 
    131           this->setValue(value); 
    132           return *this; 
    133       } 
     128  CAttributeEnum<T>& CAttributeEnum<T>::operator=(const T_enum & value) 
     129  { 
     130     this->setValue(value); 
     131     return *this; 
     132  } 
    134133 
    135134  //--------------------------------------------------------------- 
    136135 
    137136  template <class T> 
    138       StdString CAttributeEnum<T>::_toString(void) const 
    139       { 
    140           StdOStringStream oss; 
    141           if (!CEnum<T>::isEmpty() && this->hasId()) 
    142               oss << this->getName() << "=\"" << CEnum<T>::toString() << "\""; 
    143           return (oss.str()); 
    144       } 
    145  
    146   template <class T> 
    147       void CAttributeEnum<T>::_fromString(const StdString & str) 
    148       { 
    149           CEnum<T>::fromString(str); 
    150       } 
    151  
    152   template <class T> 
    153       bool CAttributeEnum<T>::_toBuffer (CBufferOut& buffer) const 
    154       { 
    155           return CEnum<T>::toBuffer(buffer); 
    156       } 
    157  
    158   template <class T> 
    159       bool CAttributeEnum<T>::_fromBuffer(CBufferIn& buffer) 
    160       { 
    161           return CEnum<T>::fromBuffer(buffer); 
    162       } 
    163  
    164   template <typename T> 
    165       void CAttributeEnum<T>::generateCInterface(ostream& oss,const string& className) 
    166       { 
    167           CInterface::AttributeCInterface<CEnumBase>(oss, className, this->getName()); 
    168       } 
    169  
    170   template <typename T> 
    171       void CAttributeEnum<T>::generateFortran2003Interface(ostream& oss,const string& className) 
    172       { 
    173           CInterface::AttributeFortran2003Interface<string>(oss, className, this->getName()); 
    174       } 
    175  
    176   template <typename T> 
    177       void CAttributeEnum<T>::generateFortranInterfaceDeclaration_(ostream& oss,const string& className) 
    178       { 
    179           CInterface::AttributeFortranInterfaceDeclaration<string>(oss, className, this->getName()+"_"); 
    180       } 
    181  
    182   template <typename T> 
    183       void CAttributeEnum<T>::generateFortranInterfaceBody_(ostream& oss,const string& className) 
    184       { 
    185           CInterface::AttributeFortranInterfaceBody<string>(oss, className, this->getName()); 
    186       } 
    187  
    188   template <typename T> 
    189       void CAttributeEnum<T>::generateFortranInterfaceDeclaration(ostream& oss,const string& className) 
    190       { 
    191           CInterface::AttributeFortranInterfaceDeclaration<string>(oss, className, this->getName()); 
    192       } 
    193  
    194   template <typename T> 
    195       void CAttributeEnum<T>::generateFortranInterfaceGetDeclaration_(ostream& oss,const string& className) 
    196       { 
    197           CInterface::AttributeFortranInterfaceGetDeclaration<string>(oss, className, this->getName()+"_"); 
    198       } 
    199  
    200   template <typename T> 
    201       void CAttributeEnum<T>::generateFortranInterfaceGetBody_(ostream& oss,const string& className) 
    202       { 
    203           CInterface::AttributeFortranInterfaceGetBody<string>(oss, className, this->getName()); 
    204       } 
    205  
    206   template <typename T> 
    207       void CAttributeEnum<T>::generateFortranInterfaceGetDeclaration(ostream& oss,const string& className) 
    208       { 
    209           CInterface::AttributeFortranInterfaceGetDeclaration<string>(oss, className, this->getName()); 
    210       } 
     137  StdString CAttributeEnum<T>::_toString(void) const 
     138  { 
     139     StdOStringStream oss; 
     140     if (!CEnum<T>::isEmpty() && this->hasId()) 
     141        oss << this->getName() << "=\"" << CEnum<T>::toString() << "\""; 
     142     return (oss.str()); 
     143  } 
     144 
     145  template <class T> 
     146  void CAttributeEnum<T>::_fromString(const StdString & str) 
     147  { 
     148    CEnum<T>::fromString(str); 
     149  } 
     150 
     151  template <class T> 
     152  bool CAttributeEnum<T>::_toBuffer (CBufferOut& buffer) const 
     153  { 
     154     return CEnum<T>::toBuffer(buffer); 
     155  } 
     156 
     157  template <class T> 
     158  bool CAttributeEnum<T>::_fromBuffer(CBufferIn& buffer) 
     159  { 
     160    return CEnum<T>::fromBuffer(buffer); 
     161  } 
     162 
     163  template <typename T> 
     164  void CAttributeEnum<T>::generateCInterface(ostream& oss,const string& className) 
     165  { 
     166    CInterface::AttributeCInterface<CEnumBase>(oss, className, this->getName()); 
     167  } 
     168 
     169  template <typename T> 
     170  void CAttributeEnum<T>::generateFortran2003Interface(ostream& oss,const string& className) 
     171  { 
     172    CInterface::AttributeFortran2003Interface<string>(oss, className, this->getName()); 
     173  } 
     174 
     175  template <typename T> 
     176  void CAttributeEnum<T>::generateFortranInterfaceDeclaration_(ostream& oss,const string& className) 
     177  { 
     178    CInterface::AttributeFortranInterfaceDeclaration<string>(oss, className, this->getName()+"_"); 
     179  } 
     180 
     181  template <typename T> 
     182  void CAttributeEnum<T>::generateFortranInterfaceBody_(ostream& oss,const string& className) 
     183  { 
     184    CInterface::AttributeFortranInterfaceBody<string>(oss, className, this->getName()); 
     185  } 
     186 
     187  template <typename T> 
     188  void CAttributeEnum<T>::generateFortranInterfaceDeclaration(ostream& oss,const string& className) 
     189  { 
     190    CInterface::AttributeFortranInterfaceDeclaration<string>(oss, className, this->getName()); 
     191  } 
     192 
     193  template <typename T> 
     194  void CAttributeEnum<T>::generateFortranInterfaceGetDeclaration_(ostream& oss,const string& className) 
     195  { 
     196    CInterface::AttributeFortranInterfaceGetDeclaration<string>(oss, className, this->getName()+"_"); 
     197  } 
     198 
     199  template <typename T> 
     200  void CAttributeEnum<T>::generateFortranInterfaceGetBody_(ostream& oss,const string& className) 
     201  { 
     202    CInterface::AttributeFortranInterfaceGetBody<string>(oss, className, this->getName()); 
     203  } 
     204 
     205  template <typename T> 
     206  void CAttributeEnum<T>::generateFortranInterfaceGetDeclaration(ostream& oss,const string& className) 
     207  { 
     208    CInterface::AttributeFortranInterfaceGetDeclaration<string>(oss, className, this->getName()); 
     209  } 
    211210} // namespace xios 
    212211 
    213212#endif // __XIOS_ATTRIBUTE_ENUM_IMPL_HPP__ 
    214  
  • XIOS/dev/branch_openmp/src/attribute_map.hpp

    r1134 r1328  
    7676            /// Propriété statique /// 
    7777            static CAttributeMap * Current; 
    78             #pragma omp threadprivate (Current) 
    7978 
    8079      };  // class CAttributeMap 
  • XIOS/dev/branch_openmp/src/buffer_client.cpp

    r1205 r1328  
    77#include "mpi.hpp" 
    88#include "tracer.hpp" 
     9 
     10 
     11using namespace ep_lib; 
    912 
    1013namespace xios 
     
    2730    buffer[1] = new char[bufferSize]; 
    2831    retBuffer = new CBufferOut(buffer[current], bufferSize); 
    29     #pragma omp critical (_output) 
    3032    info(10) << "CClientBuffer: allocated 2 x " << bufferSize << " bytes for server " << serverRank << " with a maximum of " << maxBufferedEvents << " buffered events" << endl; 
    3133  } 
  • XIOS/dev/branch_openmp/src/buffer_client.hpp

    r1205 r1328  
    66#include "mpi.hpp" 
    77#include "cxios.hpp" 
    8 #ifdef _usingEP 
    9 #include "ep_declaration.hpp" 
    10 #endif 
    118 
    129namespace xios 
     
    1613    public: 
    1714      static size_t maxRequestSize; 
    18       #pragma omp threadprivate(maxRequestSize) 
    1915 
    20       CClientBuffer(MPI_Comm intercomm, int serverRank, StdSize bufferSize, StdSize estimatedMaxEventSize, StdSize maxBufferedEvents); 
     16      CClientBuffer(ep_lib::MPI_Comm intercomm, int serverRank, StdSize bufferSize, StdSize estimatedMaxEventSize, StdSize maxBufferedEvents); 
    2117      ~CClientBuffer(); 
    2218 
     
    4036      bool pending; 
    4137 
    42       MPI_Request request; 
     38      ep_lib::MPI_Request request; 
    4339 
    4440      CBufferOut* retBuffer; 
    45       const MPI_Comm interComm; 
     41      const ep_lib::MPI_Comm interComm; 
    4642  }; 
    4743} 
  • XIOS/dev/branch_openmp/src/buffer_server.hpp

    r1134 r1328  
    44#include "xios_spl.hpp" 
    55#include "buffer.hpp" 
    6 #include "mpi_std.hpp" 
     6#include "mpi.hpp" 
    77#include "cxios.hpp" 
    88 
  • XIOS/dev/branch_openmp/src/calendar.cpp

    r1134 r1328  
    117117      const CDate& CCalendar::update(int step) 
    118118      { 
    119         #pragma omp critical (_output) 
    120         info(80)<< "update step : " << step << " timestep " << this->timestep << std::endl; 
     119        info(20) << "update step : " << step << " timestep " << this->timestep << std::endl; 
    121120        return (this->currentDate = this->getInitDate() + step * this->timestep); 
    122121      } 
  • XIOS/dev/branch_openmp/src/client.cpp

    r1205 r1328  
    1111#include "timer.hpp" 
    1212#include "buffer_client.hpp" 
    13 #include "log.hpp" 
    14  
     13using namespace ep_lib; 
    1514 
    1615namespace xios 
    1716{ 
    18     extern int test_omp_rank; 
    19     #pragma omp threadprivate(test_omp_rank) 
    2017 
    2118    MPI_Comm CClient::intraComm ; 
    2219    MPI_Comm CClient::interComm ; 
     20    //std::list<MPI_Comm> CClient::contextInterComms; 
    2321    std::list<MPI_Comm> *CClient::contextInterComms_ptr = 0; 
    2422    int CClient::serverLeader ; 
     
    2826    StdOFStream CClient::m_errorStream; 
    2927 
    30     StdOFStream CClient::array_infoStream[16]; 
    31  
    32     void CClient::initialize(const string& codeId, MPI_Comm& localComm, MPI_Comm& returnComm) 
     28    void CClient::initialize(const string& codeId,MPI_Comm& localComm,MPI_Comm& returnComm) 
    3329    { 
    3430      int initialized ; 
     
    4137      { 
    4238// localComm doesn't given 
    43  
    4439        if (localComm == MPI_COMM_NULL) 
    4540        { 
    4641          if (!is_MPI_Initialized) 
    4742          { 
    48             //MPI_Init(NULL, NULL); 
    49             int return_level; 
    50             MPI_Init_thread(NULL, NULL, 3, &return_level); 
    51             assert(return_level == 3); 
     43            MPI_Init(NULL, NULL); 
    5244          } 
    5345          CTimer::get("XIOS").resume() ; 
     
    6153          int myColor ; 
    6254          int i,c ; 
    63  
    64           MPI_Comm_size(CXios::globalComm,&size); 
     55          MPI_Comm newComm ; 
     56 
     57          MPI_Comm_size(CXios::globalComm,&size) ; 
    6558          MPI_Comm_rank(CXios::globalComm,&rank); 
    66         
    6759 
    6860          hashAll=new unsigned long[size] ; 
     
    10698            MPI_Comm_size(intraComm,&intraCommSize) ; 
    10799            MPI_Comm_rank(intraComm,&intraCommRank) ; 
    108              
    109             #pragma omp critical(_output) 
    110             { 
    111               info(10)<<"intercommCreate::client "<<test_omp_rank<< " "<< &test_omp_rank <<" intraCommSize : "<<intraCommSize 
    112                  <<" intraCommRank :"<<intraCommRank<<"  serverLeader "<< serverLeader 
    113                  <<" globalComm : "<< &(CXios::globalComm) << endl ;   
    114             } 
    115  
    116              
    117             //test_sendrecv(CXios::globalComm); 
     100            info(50)<<"intercommCreate::client "<<rank<<" intraCommSize : "<<intraCommSize 
     101                 <<" intraCommRank :"<<intraCommRank<<"  clientLeader "<< serverLeader<<endl ; 
    118102            MPI_Intercomm_create(intraComm,0,CXios::globalComm,serverLeader,0,&interComm) ; 
    119  
    120103          } 
    121104          else 
     
    140123      } 
    141124      // using OASIS 
    142       else 
     125/*      else 
    143126      { 
    144127        // localComm doesn't given 
     
    165148        else MPI_Comm_dup(intraComm,&interComm) ; 
    166149      } 
    167  
     150*/ 
    168151      MPI_Comm_dup(intraComm,&returnComm) ; 
    169  
    170152    } 
    171153 
     
    174156    { 
    175157      CContext::setCurrent(id) ; 
    176       CContext* context = CContext::create(id); 
    177  
    178       int tmp_rank; 
    179       MPI_Comm_rank(contextComm,&tmp_rank) ; 
    180        
     158      CContext* context=CContext::create(id); 
    181159      StdString idServer(id); 
    182160      idServer += "_server"; 
     
    185163      { 
    186164        int size,rank,globalRank ; 
    187         //size_t message_size ; 
    188         //int leaderRank ; 
     165        size_t message_size ; 
     166        int leaderRank ; 
    189167        MPI_Comm contextInterComm ; 
    190168 
     
    197175        CMessage msg ; 
    198176        msg<<idServer<<size<<globalRank ; 
    199  
     177//        msg<<id<<size<<globalRank ; 
    200178 
    201179        int messageSize=msg.size() ; 
     
    208186 
    209187        MPI_Intercomm_create(contextComm,0,CXios::globalComm,serverLeader,10+globalRank,&contextInterComm) ; 
    210          
    211         #pragma omp critical(_output) 
    212         info(10)<<" RANK "<< tmp_rank<<" Register new Context : "<<id<<endl ; 
    213  
     188        info(10)<<"Register new Context : "<<id<<endl ; 
    214189 
    215190        MPI_Comm inter ; 
     
    217192        MPI_Barrier(inter) ; 
    218193 
    219          
    220194        context->initClient(contextComm,contextInterComm) ; 
    221195 
    222          
     196        //contextInterComms.push_back(contextInterComm); 
    223197        if(contextInterComms_ptr == NULL) contextInterComms_ptr = new std::list<MPI_Comm>; 
    224198        contextInterComms_ptr->push_back(contextInterComm); 
    225          
    226199        MPI_Comm_free(&inter); 
    227200      } 
     
    240213        // Finally, we should return current context to context client 
    241214        CContext::setCurrent(id); 
    242          
     215 
     216        //contextInterComms.push_back(contextInterComm); 
    243217        if(contextInterComms_ptr == NULL) contextInterComms_ptr = new std::list<MPI_Comm>; 
    244218        contextInterComms_ptr->push_back(contextInterComm); 
    245  
    246219      } 
    247220    } 
     
    253226 
    254227      MPI_Comm_rank(intraComm,&rank) ; 
    255  
     228  
    256229      if (!CXios::isServer) 
    257230      { 
     
    263236      } 
    264237 
    265       for (std::list<MPI_Comm>::iterator it = contextInterComms_ptr->begin(); it != contextInterComms_ptr->end(); ++it) 
     238      //for (std::list<MPI_Comm>::iterator it = contextInterComms.begin(); it != contextInterComms.end(); it++) 
     239      for (std::list<MPI_Comm>::iterator it = contextInterComms_ptr->begin(); it != contextInterComms_ptr->end(); it++) 
    266240        MPI_Comm_free(&(*it)); 
    267        
    268241      MPI_Comm_free(&interComm); 
    269242      MPI_Comm_free(&intraComm); 
     
    274247      if (!is_MPI_Initialized) 
    275248      { 
    276         if (CXios::usingOasis) oasis_finalize(); 
    277         else MPI_Finalize(); 
     249        //if (CXios::usingOasis) oasis_finalize(); 
     250        //else 
     251        MPI_Finalize() ; 
    278252      } 
    279253       
    280       #pragma omp critical (_output) 
    281       info(20) << "Client "<<rank<<" : Client side context is finalized "<< endl ; 
    282  
    283    /*#pragma omp critical (_output) 
    284    { 
     254      info(20) << "Client side context is finalized"<<endl ; 
    285255      report(0) <<" Performance report : Whole time from XIOS init and finalize: "<< CTimer::get("XIOS init/finalize").getCumulatedTime()<<" s"<<endl ; 
    286256      report(0) <<" Performance report : total time spent for XIOS : "<< CTimer::get("XIOS").getCumulatedTime()<<" s"<<endl ; 
     
    292262      report(0)<< " Memory report : increasing it by a factor will increase performance, depending of the volume of data wrote in file at each time step of the file"<<endl ; 
    293263      report(100)<<CTimer::getAllCumulatedTime()<<endl ; 
    294    }*/ 
    295     
    296264   } 
    297265 
     
    322290 
    323291      fileNameClient << fileName << "_" << std::setfill('0') << std::setw(numDigit) << getRank() << ext; 
    324        
    325292      fb->open(fileNameClient.str().c_str(), std::ios::out); 
    326293      if (!fb->is_open()) 
    327294        ERROR("void CClient::openStream(const StdString& fileName, const StdString& ext, std::filebuf* fb)", 
    328             << std::endl << "Can not open <" << fileNameClient << "> file to write the client log(s)."); 
     295              << std::endl << "Can not open <" << fileNameClient << "> file to write the client log(s)."); 
    329296    } 
    330297 
     
    337304    void CClient::openInfoStream(const StdString& fileName) 
    338305    { 
    339       //std::filebuf* fb = m_infoStream.rdbuf(); 
    340  
    341       info_FB[omp_get_thread_num()] = array_infoStream[omp_get_thread_num()].rdbuf(); 
    342            
    343       openStream(fileName, ".out", info_FB[omp_get_thread_num()]); 
    344  
    345       info.write2File(info_FB[omp_get_thread_num()]); 
    346       report.write2File(info_FB[omp_get_thread_num()]); 
    347        
     306      std::filebuf* fb = m_infoStream.rdbuf(); 
     307      openStream(fileName, ".out", fb); 
     308 
     309      info.write2File(fb); 
     310      report.write2File(fb); 
    348311    } 
    349312 
  • XIOS/dev/branch_openmp/src/client.hpp

    r1164 r1328  
    1010    { 
    1111      public: 
    12         static void initialize(const string& codeId, MPI_Comm& localComm, MPI_Comm& returnComm); 
     12        static void initialize(const string& codeId, ep_lib::MPI_Comm& localComm, ep_lib::MPI_Comm& returnComm); 
    1313        static void finalize(void); 
    1414        static void registerContext(const string& id, ep_lib::MPI_Comm contextComm); 
    1515 
    16         static MPI_Comm intraComm; 
    17         #pragma omp threadprivate(intraComm) 
    18  
    19         static MPI_Comm interComm; 
    20         #pragma omp threadprivate(interComm) 
    21  
     16        static ep_lib::MPI_Comm intraComm; 
     17        static ep_lib::MPI_Comm interComm; 
    2218        //static std::list<MPI_Comm> contextInterComms; 
    23          
    24         static std::list<MPI_Comm> * contextInterComms_ptr; 
    25         #pragma omp threadprivate(contextInterComms_ptr) 
    26  
     19        static std::list<ep_lib::MPI_Comm> *contextInterComms_ptr; 
    2720        static int serverLeader; 
    28         #pragma omp threadprivate(serverLeader) 
    29  
    3021        static bool is_MPI_Initialized ; 
    31         #pragma omp threadprivate(is_MPI_Initialized) 
    3222 
    3323        //! Get rank of the current process 
     
    5040      protected: 
    5141        static int rank; 
    52         #pragma omp threadprivate(rank) 
    53  
    5442        static StdOFStream m_infoStream; 
    55         #pragma omp threadprivate(m_infoStream)  
    56  
    5743        static StdOFStream m_errorStream; 
    58         #pragma omp threadprivate(m_errorStream) 
    59  
    60         static StdOFStream array_infoStream[16]; 
    6144 
    6245        static void openStream(const StdString& fileName, const StdString& ext, std::filebuf* fb); 
  • XIOS/dev/branch_openmp/src/client_client_dht_template.hpp

    r1134 r1328  
    1313#include "xios_spl.hpp" 
    1414#include "array_new.hpp" 
    15 #include "mpi_std.hpp" 
     15#include "mpi.hpp" 
    1616#include "policy.hpp" 
    1717#include <boost/unordered_map.hpp> 
     
    8787                           const ep_lib::MPI_Comm& clientIntraComm, 
    8888                           std::vector<ep_lib::MPI_Request>& requestSendInfo); 
     89    void sendInfoToClients(int clientDestRank, unsigned char* info, int infoSize, 
     90                           const ep_lib::MPI_Comm& clientIntraComm, 
     91                           ep_lib::MPI_Request* requestSendInfo); 
    8992 
    9093    void recvInfoFromClients(int clientSrcRank, unsigned char* info, int infoSize, 
    9194                            const ep_lib::MPI_Comm& clientIntraComm, 
    9295                            std::vector<ep_lib::MPI_Request>& requestRecvInfo); 
     96    void recvInfoFromClients(int clientSrcRank, unsigned char* info, int infoSize, 
     97                             const ep_lib::MPI_Comm& clientIntraComm, 
     98                             ep_lib::MPI_Request* requestRecvInfo); 
     99                                                         
    93100 
    94101    // Send global index to clients 
     
    96103                            const ep_lib::MPI_Comm& clientIntraComm, 
    97104                            std::vector<ep_lib::MPI_Request>& requestSendIndexGlobal); 
     105    void sendIndexToClients(int clientDestRank, size_t* indices, size_t indiceSize, 
     106                            const ep_lib::MPI_Comm& clientIntraComm, 
     107                            ep_lib::MPI_Request* requestSendIndexGlobal); 
    98108 
    99109    void recvIndexFromClients(int clientSrcRank, size_t* indices, size_t indiceSize, 
    100110                             const ep_lib::MPI_Comm& clientIntraComm, 
    101111                             std::vector<ep_lib::MPI_Request>& requestRecvIndex); 
     112    void recvIndexFromClients(int clientSrcRank, size_t* indices, size_t indiceSize, 
     113                              const ep_lib::MPI_Comm& clientIntraComm, 
     114                              ep_lib::MPI_Request* requestRecvIndex); 
    102115 
    103116    void sendRecvOnReturn(const std::vector<int>& sendNbRank, std::vector<int>& sendNbElements, 
  • XIOS/dev/branch_openmp/src/client_client_dht_template_impl.hpp

    r1209 r1328  
    1010#include "utils.hpp" 
    1111#include "mpi_tag.hpp" 
    12 #ifdef _usingEP 
    13 #include "ep_declaration.hpp" 
    14 #include "ep_lib.hpp" 
    15 #endif 
    16  
    1712 
    1813namespace xios 
     
    2217  : H(clientIntraComm), index2InfoMapping_(), indexToInfoMappingLevel_(), nbClient_(0) 
    2318{ 
    24   MPI_Comm_size(clientIntraComm, &nbClient_); 
     19  ep_lib::MPI_Comm_size(clientIntraComm, &nbClient_); 
    2520  this->computeMPICommLevel(); 
    2621  int nbLvl = this->getNbLevel(); 
     
    4237  : H(clientIntraComm), index2InfoMapping_(), indexToInfoMappingLevel_(), nbClient_(0) 
    4338{ 
    44   MPI_Comm_size(clientIntraComm, &nbClient_); 
     39  ep_lib::MPI_Comm_size(clientIntraComm, &nbClient_); 
    4540  this->computeMPICommLevel(); 
    4641  int nbLvl = this->getNbLevel(); 
     
    6762  : H(clientIntraComm), index2InfoMapping_(), indexToInfoMappingLevel_(), nbClient_(0) 
    6863{ 
    69   MPI_Comm_size(clientIntraComm, &nbClient_); 
     64  ep_lib::MPI_Comm_size(clientIntraComm, &nbClient_); 
    7065  this->computeMPICommLevel(); 
    7166  int nbLvl = this->getNbLevel(); 
     
    10499{ 
    105100  int clientRank; 
    106   MPI_Comm_rank(commLevel,&clientRank); 
    107   //ep_lib::MPI_Barrier(commLevel); 
     101  ep_lib::MPI_Comm_rank(commLevel,&clientRank); 
    108102  int groupRankBegin = this->getGroupBegin()[level]; 
    109103  int nbClient = this->getNbInGroup()[level]; 
     
    175169    recvIndexBuff = new unsigned long[recvNbIndexCount]; 
    176170 
    177   int request_size = 0; 
    178  
    179   int currentIndex = 0; 
    180   int nbRecvClient = recvRankClient.size(); 
    181  
    182   int position = 0; 
    183  
    184   for (int idx = 0; idx < nbRecvClient; ++idx) 
     171int request_size = 0; 
     172  for (int idx = 0; idx < recvRankClient.size(); ++idx) 
    185173  { 
    186174    if (0 != recvNbIndexClientCount[idx]) 
    187     { 
    188       request_size++; 
    189     } 
     175      request_size ++; 
    190176  } 
    191177 
    192178  request_size += client2ClientIndex.size(); 
    193179 
    194  
    195180  std::vector<ep_lib::MPI_Request> request(request_size); 
    196  
     181   
    197182  std::vector<int>::iterator itbRecvIndex = recvRankClient.begin(), itRecvIndex, 
    198183                             iteRecvIndex = recvRankClient.end(), 
    199184                           itbRecvNbIndex = recvNbIndexClientCount.begin(), 
    200185                           itRecvNbIndex; 
    201    
    202    
     186  int currentIndex = 0; 
     187  int nbRecvClient = recvRankClient.size(); 
     188  int request_position = 0; 
     189  for (int idx = 0; idx < nbRecvClient; ++idx) 
     190  { 
     191    if (0 != recvNbIndexClientCount[idx]) 
     192      recvIndexFromClients(recvRankClient[idx], recvIndexBuff+currentIndex, recvNbIndexClientCount[idx], commLevel, &request[request_position++]); 
     193      //recvIndexFromClients(recvRankClient[idx], recvIndexBuff+currentIndex, recvNbIndexClientCount[idx], commLevel, request); 
     194    currentIndex += recvNbIndexClientCount[idx]; 
     195  } 
     196 
    203197  boost::unordered_map<int, size_t* >::iterator itbIndex = client2ClientIndex.begin(), itIndex, 
    204198                                                iteIndex = client2ClientIndex.end(); 
    205199  for (itIndex = itbIndex; itIndex != iteIndex; ++itIndex) 
    206   { 
    207     MPI_Isend(itIndex->second, sendNbIndexBuff[itIndex->first-groupRankBegin], MPI_UNSIGNED_LONG, itIndex->first, MPI_DHT_INDEX, commLevel, &request[position]); 
    208     position++; 
     200    sendIndexToClients(itIndex->first, (itIndex->second), sendNbIndexBuff[itIndex->first-groupRankBegin], commLevel, &request[request_position++]); 
    209201    //sendIndexToClients(itIndex->first, (itIndex->second), sendNbIndexBuff[itIndex->first-groupRankBegin], commLevel, request); 
    210   } 
    211      
    212   for (int idx = 0; idx < nbRecvClient; ++idx) 
    213   { 
    214     if (0 != recvNbIndexClientCount[idx]) 
    215     { 
    216       MPI_Irecv(recvIndexBuff+currentIndex, recvNbIndexClientCount[idx], MPI_UNSIGNED_LONG, 
    217             recvRankClient[idx], MPI_DHT_INDEX, commLevel, &request[position]); 
    218       position++; 
    219       //recvIndexFromClients(recvRankClient[idx], recvIndexBuff+currentIndex, recvNbIndexClientCount[idx], commLevel, request); 
    220     } 
    221     currentIndex += recvNbIndexClientCount[idx]; 
    222   } 
    223  
    224    
    225   std::vector<ep_lib::MPI_Status> status(request_size); 
    226   MPI_Waitall(request.size(), &request[0], &status[0]); 
    227  
     202 
     203  std::vector<ep_lib::MPI_Status> status(request.size()); 
     204  ep_lib::MPI_Waitall(request.size(), &request[0], &status[0]); 
    228205 
    229206  CArray<size_t,1>* tmpGlobalIndex; 
     
    238215    --level; 
    239216    computeIndexInfoMappingLevel(*tmpGlobalIndex, this->internalComm_, level); 
    240        
    241217  } 
    242218  else // Now, we are in the last level where necessary mappings are. 
     
    279255  } 
    280256 
    281   request_size = 0; 
     257int requestOnReturn_size=0; 
    282258  for (int idx = 0; idx < recvRankOnReturn.size(); ++idx) 
    283259  { 
    284260    if (0 != recvNbIndexOnReturn[idx]) 
    285261    { 
    286       request_size += 2; 
     262      requestOnReturn_size += 2; 
    287263    } 
    288264  } 
     
    292268    if (0 != sendNbIndexOnReturn[idx]) 
    293269    { 
    294       request_size += 2; 
    295     } 
    296   } 
    297  
    298   std::vector<ep_lib::MPI_Request> requestOnReturn(request_size); 
     270      requestOnReturn_size += 2; 
     271    } 
     272  } 
     273 
     274  int requestOnReturn_position=0; 
     275 
     276  std::vector<ep_lib::MPI_Request> requestOnReturn(requestOnReturn_size); 
    299277  currentIndex = 0; 
    300   position = 0; 
    301278  for (int idx = 0; idx < recvRankOnReturn.size(); ++idx) 
    302279  { 
     
    304281    { 
    305282      //recvIndexFromClients(recvRankOnReturn[idx], recvIndexBuffOnReturn+currentIndex, recvNbIndexOnReturn[idx], commLevel, requestOnReturn); 
    306       MPI_Irecv(recvIndexBuffOnReturn+currentIndex, recvNbIndexOnReturn[idx], MPI_UNSIGNED_LONG, 
    307             recvRankOnReturn[idx], MPI_DHT_INDEX, commLevel, &requestOnReturn[position]); 
    308       position++; 
    309283      //recvInfoFromClients(recvRankOnReturn[idx], 
    310284      //                    recvInfoBuffOnReturn+currentIndex*ProcessDHTElement<InfoType>::typeSize(), 
    311285      //                    recvNbIndexOnReturn[idx]*ProcessDHTElement<InfoType>::typeSize(), 
    312286      //                    commLevel, requestOnReturn); 
    313       MPI_Irecv(recvInfoBuffOnReturn+currentIndex*ProcessDHTElement<InfoType>::typeSize(),  
    314                 recvNbIndexOnReturn[idx]*ProcessDHTElement<InfoType>::typeSize(), MPI_CHAR, 
    315                 recvRankOnReturn[idx], MPI_DHT_INFO, commLevel, &requestOnReturn[position]); 
    316       position++; 
     287      recvIndexFromClients(recvRankOnReturn[idx], recvIndexBuffOnReturn+currentIndex, recvNbIndexOnReturn[idx], commLevel, &requestOnReturn[requestOnReturn_position++]); 
     288      recvInfoFromClients(recvRankOnReturn[idx], 
     289                          recvInfoBuffOnReturn+currentIndex*ProcessDHTElement<InfoType>::typeSize(), 
     290                          recvNbIndexOnReturn[idx]*ProcessDHTElement<InfoType>::typeSize(), 
     291                          commLevel, &requestOnReturn[requestOnReturn_position++]); 
    317292    } 
    318293    currentIndex += recvNbIndexOnReturn[idx]; 
     
    349324      //sendIndexToClients(rank, client2ClientIndexOnReturn[rank], 
    350325      //                   sendNbIndexOnReturn[idx], commLevel, requestOnReturn); 
    351       MPI_Isend(client2ClientIndexOnReturn[rank], sendNbIndexOnReturn[idx], MPI_UNSIGNED_LONG, 
    352             rank, MPI_DHT_INDEX, commLevel, &requestOnReturn[position]); 
    353       position++; 
    354326      //sendInfoToClients(rank, client2ClientInfoOnReturn[rank], 
    355327      //                  sendNbIndexOnReturn[idx]*ProcessDHTElement<InfoType>::typeSize(), commLevel, requestOnReturn); 
    356       MPI_Isend(client2ClientInfoOnReturn[rank], sendNbIndexOnReturn[idx]*ProcessDHTElement<InfoType>::typeSize(), MPI_CHAR, 
    357             rank, MPI_DHT_INFO, commLevel, &requestOnReturn[position]); 
    358       position++; 
     328      sendIndexToClients(rank, client2ClientIndexOnReturn[rank], 
     329                         sendNbIndexOnReturn[idx], commLevel, &requestOnReturn[requestOnReturn_position++]); 
     330      sendInfoToClients(rank, client2ClientInfoOnReturn[rank], 
     331                        sendNbIndexOnReturn[idx]*ProcessDHTElement<InfoType>::typeSize(), commLevel, &requestOnReturn[requestOnReturn_position++]); 
     332       
    359333    } 
    360334    currentIndex += recvNbIndexClientCount[idx]; 
     
    362336 
    363337  std::vector<ep_lib::MPI_Status> statusOnReturn(requestOnReturn.size()); 
    364   MPI_Waitall(requestOnReturn.size(), &requestOnReturn[0], &statusOnReturn[0]); 
     338  ep_lib::MPI_Waitall(requestOnReturn.size(), &requestOnReturn[0], &statusOnReturn[0]); 
    365339 
    366340  Index2VectorInfoTypeMap indexToInfoMapping; 
     
    432406{ 
    433407  int clientRank; 
    434   MPI_Comm_rank(commLevel,&clientRank); 
     408  ep_lib::MPI_Comm_rank(commLevel,&clientRank); 
    435409  computeSendRecvRank(level, clientRank); 
    436   //ep_lib::MPI_Barrier(commLevel); 
    437410 
    438411  int groupRankBegin = this->getGroupBegin()[level]; 
     
    481454    { 
    482455      client2ClientIndex[indexClient + groupRankBegin][sendNbIndexBuff[indexClient]] = it->first;; 
     456  //          ProcessDHTElement<InfoType>::packElement(it->second, client2ClientInfo[indexClient + groupRankBegin], sendNbInfo[indexClient]); 
    483457      ProcessDHTElement<InfoType>::packElement(infoTmp[idx], client2ClientInfo[indexClient + groupRankBegin], sendNbInfo[indexClient]); 
    484458      ++sendNbIndexBuff[indexClient]; 
     
    494468  int recvNbIndexCount = 0; 
    495469  for (int idx = 0; idx < recvNbIndexClientCount.size(); ++idx) 
    496   {  
    497470    recvNbIndexCount += recvNbIndexClientCount[idx]; 
    498   } 
    499471 
    500472  unsigned long* recvIndexBuff; 
     
    509481  // it will send a message to the correct clients. 
    510482  // Contents of the message are index and its corresponding informatioin 
    511   int request_size = 0;   
     483  int request_size = 0; 
     484   for (int idx = 0; idx < recvRankClient.size(); ++idx) 
     485   { 
     486     if (0 != recvNbIndexClientCount[idx]) 
     487     { 
     488       request_size += 2; 
     489     } 
     490   } 
     491  
     492   request_size += client2ClientIndex.size(); 
     493   request_size += client2ClientInfo.size(); 
     494  
     495   std::vector<ep_lib::MPI_Request> request(request_size); 
    512496  int currentIndex = 0; 
    513497  int nbRecvClient = recvRankClient.size(); 
    514   int current_pos = 0;  
    515  
     498  int request_position=0; 
    516499  for (int idx = 0; idx < nbRecvClient; ++idx) 
    517500  { 
    518501    if (0 != recvNbIndexClientCount[idx]) 
    519502    { 
    520       request_size += 2; 
    521     } 
    522     //currentIndex += recvNbIndexClientCount[idx]; 
    523   } 
    524  
    525   request_size += client2ClientIndex.size(); 
    526   request_size += client2ClientInfo.size(); 
    527  
    528  
    529  
    530   std::vector<ep_lib::MPI_Request> request(request_size); 
    531    
    532   //unsigned long* tmp_send_buf_long[client2ClientIndex.size()]; 
    533   //unsigned char* tmp_send_buf_char[client2ClientInfo.size()]; 
    534    
    535   int info_position = 0; 
    536   int index_position = 0; 
    537  
     503      //recvIndexFromClients(recvRankClient[idx], recvIndexBuff+currentIndex, recvNbIndexClientCount[idx], commLevel, request); 
     504      //recvInfoFromClients(recvRankClient[idx], 
     505      //                    recvInfoBuff+currentIndex*ProcessDHTElement<InfoType>::typeSize(), 
     506      //                    recvNbIndexClientCount[idx]*ProcessDHTElement<InfoType>::typeSize(), 
     507      //                    commLevel, request); 
     508        recvIndexFromClients(recvRankClient[idx], recvIndexBuff+currentIndex, recvNbIndexClientCount[idx], commLevel, &request[request_position++]); 
     509        recvInfoFromClients(recvRankClient[idx], 
     510                            recvInfoBuff+currentIndex*ProcessDHTElement<InfoType>::typeSize(), 
     511                            recvNbIndexClientCount[idx]*ProcessDHTElement<InfoType>::typeSize(), 
     512                            commLevel, &request[request_position++]); 
     513    } 
     514    currentIndex += recvNbIndexClientCount[idx]; 
     515  } 
    538516 
    539517  boost::unordered_map<int, size_t* >::iterator itbIndex = client2ClientIndex.begin(), itIndex, 
    540518                                                iteIndex = client2ClientIndex.end(); 
    541519  for (itIndex = itbIndex; itIndex != iteIndex; ++itIndex) 
    542   { 
     520    sendIndexToClients(itIndex->first, itIndex->second, sendNbIndexBuff[itIndex->first-groupRankBegin], commLevel, &request[request_position++]); 
    543521    //sendIndexToClients(itIndex->first, itIndex->second, sendNbIndexBuff[itIndex->first-groupRankBegin], commLevel, request); 
    544  
    545     //tmp_send_buf_long[index_position] = new unsigned long[sendNbIndexBuff[itIndex->first-groupRankBegin]]; 
    546     //for(int i=0; i<sendNbIndexBuff[itIndex->first-groupRankBegin]; i++) 
    547     //{ 
    548     //  tmp_send_buf_long[index_position][i] = (static_cast<unsigned long * >(itIndex->second))[i]; 
    549     //} 
    550     //MPI_Isend(tmp_send_buf_long[current_pos], sendNbIndexBuff[itIndex->first-groupRankBegin], MPI_UNSIGNED_LONG, 
    551     MPI_Isend(itIndex->second, sendNbIndexBuff[itIndex->first-groupRankBegin], MPI_UNSIGNED_LONG, 
    552               itIndex->first, MPI_DHT_INDEX, commLevel, &request[current_pos]); 
    553     current_pos++;  
    554     index_position++; 
    555  
    556   } 
    557  
    558522  boost::unordered_map<int, unsigned char*>::iterator itbInfo = client2ClientInfo.begin(), itInfo, 
    559523                                                      iteInfo = client2ClientInfo.end(); 
    560524  for (itInfo = itbInfo; itInfo != iteInfo; ++itInfo) 
    561   { 
     525    sendInfoToClients(itInfo->first, itInfo->second, sendNbInfo[itInfo->first-groupRankBegin], commLevel, &request[request_position++]); 
    562526    //sendInfoToClients(itInfo->first, itInfo->second, sendNbInfo[itInfo->first-groupRankBegin], commLevel, request); 
    563527 
    564     //tmp_send_buf_char[info_position] = new unsigned char[sendNbInfo[itInfo->first-groupRankBegin]]; 
    565     //for(int i=0; i<sendNbInfo[itInfo->first-groupRankBegin]; i++) 
    566     //{ 
    567     //  tmp_send_buf_char[info_position][i] = (static_cast<unsigned char * >(itInfo->second))[i]; 
    568     //} 
    569  
    570     MPI_Isend(itInfo->second, sendNbInfo[itInfo->first-groupRankBegin], MPI_CHAR, 
    571               itInfo->first, MPI_DHT_INFO, commLevel, &request[current_pos]); 
    572  
    573     current_pos++; 
    574     info_position++; 
    575   } 
    576    
    577   for (int idx = 0; idx < nbRecvClient; ++idx) 
    578   { 
    579     if (0 != recvNbIndexClientCount[idx]) 
    580     { 
    581       //recvIndexFromClients(recvRankClient[idx], recvIndexBuff+currentIndex, recvNbIndexClientCount[idx], commLevel, request); 
    582       MPI_Irecv(recvIndexBuff+currentIndex, recvNbIndexClientCount[idx], MPI_UNSIGNED_LONG, 
    583                 recvRankClient[idx], MPI_DHT_INDEX, commLevel, &request[current_pos]); 
    584       current_pos++; 
    585        
    586        
    587       MPI_Irecv(recvInfoBuff+currentIndex*ProcessDHTElement<InfoType>::typeSize(),  
    588                 recvNbIndexClientCount[idx]*ProcessDHTElement<InfoType>::typeSize(),  
    589                 MPI_CHAR, recvRankClient[idx], MPI_DHT_INFO, commLevel, &request[current_pos]); 
    590        
    591       current_pos++; 
    592        
    593  
    594  
    595       // recvInfoFromClients(recvRankClient[idx], 
    596       //                     recvInfoBuff+currentIndex*ProcessDHTElement<InfoType>::typeSize(), 
    597       //                     recvNbIndexClientCount[idx]*ProcessDHTElement<InfoType>::typeSize(), 
    598       //                     commLevel, request); 
    599     } 
    600     currentIndex += recvNbIndexClientCount[idx]; 
    601   } 
    602  
    603528  std::vector<ep_lib::MPI_Status> status(request.size()); 
    604    
    605   MPI_Waitall(request.size(), &request[0], &status[0]); 
    606    
    607   
    608   //for(int i=0; i<client2ClientInfo.size(); i++) 
    609   //  delete[] tmp_send_buf_char[i]; 
    610  
    611    
    612  
    613   //for(int i=0; i<client2ClientIndex.size(); i++) 
    614   //  delete[] tmp_send_buf_long[i]; 
    615  
     529  ep_lib::MPI_Waitall(request.size(), &request[0], &status[0]); 
    616530 
    617531  Index2VectorInfoTypeMap indexToInfoMapping; 
     
    654568  else 
    655569    index2InfoMapping_.swap(indexToInfoMapping); 
    656    
    657570} 
    658571 
     
    670583                                                       std::vector<ep_lib::MPI_Request>& requestSendIndex) 
    671584{ 
    672   printf("should not call this function sendIndexToClients");   
    673585  ep_lib::MPI_Request request; 
    674586  requestSendIndex.push_back(request); 
    675   MPI_Isend(indices, indiceSize, MPI_UNSIGNED_LONG, 
     587  ep_lib::MPI_Isend(indices, indiceSize, MPI_UNSIGNED_LONG, 
    676588            clientDestRank, MPI_DHT_INDEX, clientIntraComm, &(requestSendIndex.back())); 
     589} 
     590 
     591/*! 
     592  Send message containing index to clients 
     593  \param [in] clientDestRank rank of destination client 
     594  \param [in] indices index to send 
     595  \param [in] indiceSize size of index array to send 
     596  \param [in] clientIntraComm communication group of client 
     597  \param [in] requestSendIndex sending request 
     598*/ 
     599template<typename T, typename H> 
     600void CClientClientDHTTemplate<T,H>::sendIndexToClients(int clientDestRank, size_t* indices, size_t indiceSize, 
     601                                                       const ep_lib::MPI_Comm& clientIntraComm, 
     602                                                       ep_lib::MPI_Request* requestSendIndex) 
     603{ 
     604  ep_lib::MPI_Isend(indices, indiceSize, MPI_UNSIGNED_LONG, 
     605            clientDestRank, MPI_DHT_INDEX, clientIntraComm, requestSendIndex); 
    677606} 
    678607 
     
    689618                                                         std::vector<ep_lib::MPI_Request>& requestRecvIndex) 
    690619{ 
    691   printf("should not call this function recvIndexFromClients"); 
    692620  ep_lib::MPI_Request request; 
    693621  requestRecvIndex.push_back(request); 
    694   MPI_Irecv(indices, indiceSize, MPI_UNSIGNED_LONG, 
     622  ep_lib::MPI_Irecv(indices, indiceSize, MPI_UNSIGNED_LONG, 
    695623            clientSrcRank, MPI_DHT_INDEX, clientIntraComm, &(requestRecvIndex.back())); 
     624} 
     625 
     626/*! 
     627  Receive message containing index to clients 
     628  \param [in] clientDestRank rank of destination client 
     629  \param [in] indices index to send 
     630  \param [in] clientIntraComm communication group of client 
     631  \param [in] requestRecvIndex receiving request 
     632*/ 
     633template<typename T, typename H> 
     634void CClientClientDHTTemplate<T,H>::recvIndexFromClients(int clientSrcRank, size_t* indices, size_t indiceSize, 
     635                                                         const ep_lib::MPI_Comm& clientIntraComm, 
     636                                                         ep_lib::MPI_Request *requestRecvIndex) 
     637{ 
     638  ep_lib::MPI_Irecv(indices, indiceSize, MPI_UNSIGNED_LONG, 
     639            clientSrcRank, MPI_DHT_INDEX, clientIntraComm, requestRecvIndex); 
    696640} 
    697641 
     
    709653                                                      std::vector<ep_lib::MPI_Request>& requestSendInfo) 
    710654{ 
    711   printf("should not call this function sendInfoToClients"); 
    712655  ep_lib::MPI_Request request; 
    713656  requestSendInfo.push_back(request); 
    714   MPI_Isend(info, infoSize, MPI_CHAR, 
     657 
     658  ep_lib::MPI_Isend(info, infoSize, MPI_CHAR, 
    715659            clientDestRank, MPI_DHT_INFO, clientIntraComm, &(requestSendInfo.back())); 
     660} 
     661 
     662/*! 
     663  Send message containing information to clients 
     664  \param [in] clientDestRank rank of destination client 
     665  \param [in] info info array to send 
     666  \param [in] infoSize info array size to send 
     667  \param [in] clientIntraComm communication group of client 
     668  \param [in] requestSendInfo sending request 
     669*/ 
     670template<typename T, typename H> 
     671void CClientClientDHTTemplate<T,H>::sendInfoToClients(int clientDestRank, unsigned char* info, int infoSize, 
     672                                                      const ep_lib::MPI_Comm& clientIntraComm, 
     673                                                      ep_lib::MPI_Request *requestSendInfo) 
     674{ 
     675  ep_lib::MPI_Isend(info, infoSize, MPI_CHAR, 
     676            clientDestRank, MPI_DHT_INFO, clientIntraComm, requestSendInfo); 
    716677} 
    717678 
     
    729690                                                        std::vector<ep_lib::MPI_Request>& requestRecvInfo) 
    730691{ 
    731   printf("should not call this function recvInfoFromClients\n"); 
    732692  ep_lib::MPI_Request request; 
    733693  requestRecvInfo.push_back(request); 
    734694 
    735   MPI_Irecv(info, infoSize, MPI_CHAR, 
     695  ep_lib::MPI_Irecv(info, infoSize, MPI_CHAR, 
    736696            clientSrcRank, MPI_DHT_INFO, clientIntraComm, &(requestRecvInfo.back())); 
     697} 
     698 
     699/*! 
     700  Receive message containing information from other clients 
     701  \param [in] clientDestRank rank of destination client 
     702  \param [in] info info array to receive 
     703  \param [in] infoSize info array size to receive 
     704  \param [in] clientIntraComm communication group of client 
     705  \param [in] requestRecvInfo list of receiving request 
     706*/ 
     707template<typename T, typename H> 
     708void CClientClientDHTTemplate<T,H>::recvInfoFromClients(int clientSrcRank, unsigned char* info, int infoSize, 
     709                                                        const ep_lib::MPI_Comm& clientIntraComm, 
     710                                                        ep_lib::MPI_Request* requestRecvInfo) 
     711{ 
     712  ep_lib::MPI_Irecv(info, infoSize, MPI_CHAR, 
     713            clientSrcRank, MPI_DHT_INFO, clientIntraComm, requestRecvInfo); 
    737714} 
    738715 
     
    807784 
    808785  int nRequest = 0; 
    809    
     786  for (int idx = 0; idx < recvNbRank.size(); ++idx) 
     787  { 
     788    ep_lib::MPI_Irecv(&recvNbElements[0]+idx, 1, MPI_INT, 
     789              recvNbRank[idx], MPI_DHT_INDEX_1, this->internalComm_, &request[nRequest]); 
     790    ++nRequest; 
     791  } 
    810792 
    811793  for (int idx = 0; idx < sendNbRank.size(); ++idx) 
    812794  { 
    813     MPI_Isend(&sendNbElements[0]+idx, 1, MPI_INT, 
     795    ep_lib::MPI_Isend(&sendNbElements[0]+idx, 1, MPI_INT, 
    814796              sendNbRank[idx], MPI_DHT_INDEX_1, this->internalComm_, &request[nRequest]); 
    815797    ++nRequest; 
    816798  } 
    817    
    818   for (int idx = 0; idx < recvNbRank.size(); ++idx) 
    819   { 
    820     MPI_Irecv(&recvNbElements[0]+idx, 1, MPI_INT, 
    821               recvNbRank[idx], MPI_DHT_INDEX_1, this->internalComm_, &request[nRequest]); 
    822     ++nRequest; 
    823   } 
    824  
    825   MPI_Waitall(sendNbRank.size()+recvNbRank.size(), &request[0], &requestStatus[0]); 
     799 
     800  ep_lib::MPI_Waitall(sendNbRank.size()+recvNbRank.size(), &request[0], &requestStatus[0]); 
    826801} 
    827802 
     
    852827  std::vector<ep_lib::MPI_Request> request(sendBuffSize+recvBuffSize); 
    853828  std::vector<ep_lib::MPI_Status> requestStatus(sendBuffSize+recvBuffSize); 
    854   //ep_lib::MPI_Request request[sendBuffSize+recvBuffSize]; 
    855   //ep_lib::MPI_Status requestStatus[sendBuffSize+recvBuffSize]; 
    856    
    857   int my_rank; 
    858   MPI_Comm_rank(this->internalComm_, &my_rank); 
    859    
     829 
    860830  int nRequest = 0; 
    861831  for (int idx = 0; idx < recvBuffSize; ++idx) 
    862832  { 
    863     MPI_Irecv(&recvBuff[2*idx], 2, MPI_INT, 
     833    ep_lib::MPI_Irecv(&recvBuff[0]+2*idx, 2, MPI_INT, 
    864834              recvRank[idx], MPI_DHT_INDEX_0, this->internalComm_, &request[nRequest]); 
    865835    ++nRequest; 
    866836  } 
    867    
    868837 
    869838  for (int idx = 0; idx < sendBuffSize; ++idx) 
     
    873842    sendBuff[idx*2+1] = sendNbElements[offSet]; 
    874843  } 
    875    
    876    
    877844 
    878845  for (int idx = 0; idx < sendBuffSize; ++idx) 
    879846  { 
    880     MPI_Isend(&sendBuff[idx*2], 2, MPI_INT, 
     847    ep_lib::MPI_Isend(&sendBuff[idx*2], 2, MPI_INT, 
    881848              sendRank[idx], MPI_DHT_INDEX_0, this->internalComm_, &request[nRequest]); 
    882849    ++nRequest; 
    883850  } 
    884    
    885    
    886  
    887   //MPI_Barrier(this->internalComm_); 
    888  
    889   MPI_Waitall(sendBuffSize+recvBuffSize, &request[0], &requestStatus[0]); 
    890   //MPI_Waitall(sendBuffSize+recvBuffSize, request, requestStatus); 
    891  
    892    
     851 
     852  ep_lib::MPI_Waitall(sendBuffSize+recvBuffSize, &request[0], &requestStatus[0]); 
    893853  int nbRecvRank = 0, nbRecvElements = 0; 
    894854  recvNbRank.clear(); 
     
    902862    } 
    903863  } 
    904  
    905  
    906    
    907    
    908 } 
    909  
    910 } 
    911  
     864} 
     865 
     866} 
  • XIOS/dev/branch_openmp/src/client_server_mapping.cpp

    r1287 r1328  
    88 */ 
    99#include "client_server_mapping.hpp" 
     10using namespace ep_lib; 
    1011 
    1112namespace xios { 
     
    6465  MPI_Allgather(&nbConnectedServer,1,MPI_INT,recvCount,1,MPI_INT,clientIntraComm) ; 
    6566 
    66    
    67   // for(int i=0; i<nbClient; i++) 
    68   //   printf("MPI_Allgather : recvCount[%d] = %d\n", i, recvCount[i]); 
    69  
    7067  displ[0]=0 ; 
    7168  for(int n=1;n<nbClient;n++) displ[n]=displ[n-1]+recvCount[n-1] ; 
     
    7572 
    7673  MPI_Allgatherv(sendBuff,nbConnectedServer,MPI_INT,recvBuff,recvCount,displ,MPI_INT,clientIntraComm) ; 
    77  
    78   // for(int i=0; i<recvSize; i++) 
    79   //   printf("MPI_Allgatherv : recvBuff[%d] = %d\n", i, recvBuff[i]); 
    80  
    81  
    8274  for(int n=0;n<recvSize;n++) clientRes[recvBuff[n]]++ ; 
    8375 
  • XIOS/dev/branch_openmp/src/client_server_mapping.hpp

    r1134 r1328  
    1414#include "mpi.hpp" 
    1515#include <boost/unordered_map.hpp> 
    16 #ifdef _usingEP 
    17 #include "ep_declaration.hpp" 
    18 #endif 
    19  
    2016 
    2117namespace xios { 
  • XIOS/dev/branch_openmp/src/client_server_mapping_distributed.cpp

    r907 r1328  
    1515#include "context.hpp" 
    1616#include "context_client.hpp" 
     17using namespace ep_lib; 
    1718 
    1819namespace xios 
  • XIOS/dev/branch_openmp/src/context_client.cpp

    r1205 r1328  
    1111#include "timer.hpp" 
    1212#include "cxios.hpp" 
     13using namespace ep_lib; 
    1314 
    1415namespace xios 
     
    2021    \cxtSer [in] cxtSer Pointer to context of server side. (It is only used on case of attached mode) 
    2122    */ 
    22     CContextClient::CContextClient(CContext* parent, ep_lib::MPI_Comm intraComm_, ep_lib::MPI_Comm interComm_, CContext* cxtSer) 
     23    CContextClient::CContextClient(CContext* parent, MPI_Comm intraComm_, MPI_Comm interComm_, CContext* cxtSer) 
    2324     : mapBufferSize_(), parentServer(cxtSer), maxBufferedEvents(4) 
    2425    { 
     
    293294       if (ratio < minBufferSizeEventSizeRatio) minBufferSizeEventSizeRatio = ratio; 
    294295     } 
    295      #ifdef _usingMPI 
    296      MPI_Allreduce(MPI_IN_PLACE, &minBufferSizeEventSizeRatio, 1, MPI_DOUBLE, MPI_MIN, intraComm); 
    297      #elif _usingEP 
     296     //MPI_Allreduce(MPI_IN_PLACE, &minBufferSizeEventSizeRatio, 1, MPI_DOUBLE, MPI_MIN, intraComm); 
    298297     MPI_Allreduce(&minBufferSizeEventSizeRatio, &minBufferSizeEventSizeRatio, 1, MPI_DOUBLE, MPI_MIN, intraComm); 
    299      #endif 
    300       
     298 
    301299     if (minBufferSizeEventSizeRatio < 1.0) 
    302300     { 
     
    402400     for (itMap = itbMap; itMap != iteMap; ++itMap) 
    403401     { 
    404        //report(10) << " Memory report : Context <" << context->getId() << "> : client side : memory used for buffer of each connection to server" << endl 
    405        //           << "  +) To server with rank " << itMap->first << " : " << itMap->second << " bytes " << endl; 
     402       report(10) << " Memory report : Context <" << context->getId() << "> : client side : memory used for buffer of each connection to server" << endl 
     403                  << "  +) To server with rank " << itMap->first << " : " << itMap->second << " bytes " << endl; 
    406404       totalBuf += itMap->second; 
    407405     } 
    408      //report(0) << " Memory report : Context <" << context->getId() << "> : client side : total memory used for buffer " << totalBuf << " bytes" << endl; 
     406     report(0) << " Memory report : Context <" << context->getId() << "> : client side : total memory used for buffer " << totalBuf << " bytes" << endl; 
    409407 
    410408     releaseBuffers(); 
  • XIOS/dev/branch_openmp/src/context_server.cpp

    r1179 r1328  
    1010#include "file.hpp" 
    1111#include "grid.hpp" 
    12 #include "mpi_std.hpp" 
     12#include "mpi.hpp" 
    1313#include "tracer.hpp" 
    1414#include "timer.hpp" 
     
    1818#include <boost/functional/hash.hpp> 
    1919 
    20  
     20using namespace ep_lib; 
    2121 
    2222namespace xios 
    2323{ 
    2424 
    25   CContextServer::CContextServer(CContext* parent, ep_lib::MPI_Comm intraComm_, ep_lib::MPI_Comm interComm_) 
     25  CContextServer::CContextServer(CContext* parent,MPI_Comm intraComm_,MPI_Comm interComm_) 
    2626  { 
    2727    context=parent; 
     
    7272    int count; 
    7373    char * addr; 
    74     ep_lib::MPI_Status status; 
     74    MPI_Status status; 
    7575    map<int,CServerBuffer*>::iterator it; 
    76  
    77     for(rank=0;rank<commSize;rank++) 
    78     { 
     76    bool okLoop; 
     77 
     78    traceOff(); 
     79    MPI_Iprobe(-2, 20,interComm,&flag,&status); 
     80    traceOn(); 
     81 
     82    if (flag==true) 
     83    { 
     84      #ifdef _usingMPI 
     85      rank=status.MPI_SOURCE ; 
     86      #elif _usingEP 
     87      rank=status.ep_src ; 
     88      #endif 
     89      okLoop = true; 
    7990      if (pendingRequest.find(rank)==pendingRequest.end()) 
    80       { 
    81         traceOff(); 
    82         ep_lib::MPI_Iprobe(rank,20,interComm,&flag,&status); 
    83         traceOn(); 
    84         if (flag) 
     91        okLoop = !listenPendingRequest(status) ; 
     92      if (okLoop) 
     93      { 
     94        for(rank=0;rank<commSize;rank++) 
    8595        { 
    86           it=buffers.find(rank); 
    87           if (it==buffers.end()) // Receive the buffer size and allocate the buffer 
     96          if (pendingRequest.find(rank)==pendingRequest.end()) 
    8897          { 
    89             StdSize buffSize = 0; 
    90             ep_lib::MPI_Recv(&buffSize, 1, MPI_LONG, rank, 20, interComm, &status); 
    91             mapBufferSize_.insert(std::make_pair(rank, buffSize)); 
    92             it=(buffers.insert(pair<int,CServerBuffer*>(rank,new CServerBuffer(buffSize)))).first; 
    93           } 
    94           else 
    95           { 
    96              
    97             ep_lib::MPI_Get_count(&status,MPI_CHAR,&count); 
    98             if (it->second->isBufferFree(count)) 
    99             { 
    100               addr=(char*)it->second->getBuffer(count); 
    101               ep_lib::MPI_Irecv(addr,count,MPI_CHAR,rank,20,interComm,&pendingRequest[rank]); 
    102               bufferRequest[rank]=addr; 
    103             } 
     98 
     99            traceOff(); 
     100            MPI_Iprobe(rank, 20,interComm,&flag,&status); 
     101            traceOn(); 
     102            if (flag==true) listenPendingRequest(status) ; 
    104103          } 
    105104        } 
     
    108107  } 
    109108 
     109  bool CContextServer::listenPendingRequest(MPI_Status& status) 
     110  { 
     111    int count; 
     112    char * addr; 
     113    map<int,CServerBuffer*>::iterator it; 
     114    #ifdef _usingMPI 
     115    int rank=status.MPI_SOURCE ; 
     116    #elif _usingEP 
     117    int rank=status.ep_src; 
     118    #endif 
     119 
     120    it=buffers.find(rank); 
     121    if (it==buffers.end()) // Receive the buffer size and allocate the buffer 
     122    { 
     123       StdSize buffSize = 0; 
     124       MPI_Recv(&buffSize, 1, MPI_LONG, rank, 20, interComm, &status); 
     125       mapBufferSize_.insert(std::make_pair(rank, buffSize)); 
     126       it=(buffers.insert(pair<int,CServerBuffer*>(rank,new CServerBuffer(buffSize)))).first; 
     127       return true; 
     128    } 
     129    else 
     130    { 
     131      MPI_Get_count(&status,MPI_CHAR,&count); 
     132      if (it->second->isBufferFree(count)) 
     133      { 
     134         addr=(char*)it->second->getBuffer(count); 
     135         MPI_Irecv(addr,count,MPI_CHAR,rank,20,interComm,&pendingRequest[rank]); 
     136         bufferRequest[rank]=addr; 
     137         return true; 
     138      } 
     139      else 
     140        return false; 
     141    } 
     142  } 
     143 
    110144  void CContextServer::checkPendingRequest(void) 
    111145  { 
    112     map<int,ep_lib::MPI_Request>::iterator it; 
     146    map<int,MPI_Request>::iterator it; 
    113147    list<int> recvRequest; 
    114148    list<int>::iterator itRecv; 
     
    116150    int flag; 
    117151    int count; 
    118     ep_lib::MPI_Status status; 
    119  
    120     for(it=pendingRequest.begin();it!=pendingRequest.end();++it) 
     152    MPI_Status status; 
     153 
     154    for(it=pendingRequest.begin();it!=pendingRequest.end();it++) 
    121155    { 
    122156      rank=it->first; 
    123157      traceOff(); 
    124       ep_lib::MPI_Test(& it->second, &flag, &status); 
     158      MPI_Test(& it->second, &flag, &status); 
    125159      traceOn(); 
    126160      if (flag==true) 
    127161      { 
    128162        recvRequest.push_back(rank); 
    129         ep_lib::MPI_Get_count(&status,MPI_CHAR,&count); 
     163        MPI_Get_count(&status,MPI_CHAR,&count); 
    130164        processRequest(rank,bufferRequest[rank],count); 
    131165      } 
     
    220254    { 
    221255      finished=true; 
    222       #pragma omp critical (_output) 
    223256      info(20)<<"Server Side context <"<<context->getId()<<"> finalized"<<endl; 
    224257      std::map<int, StdSize>::const_iterator itbMap = mapBufferSize_.begin(), 
     
    227260      for (itMap = itbMap; itMap != iteMap; ++itMap) 
    228261      { 
    229         //report(10)<< " Memory report : Context <"<<context->getId()<<"> : server side : memory used for buffer of each connection to client" << endl 
    230         //          << "  +) With client of rank " << itMap->first << " : " << itMap->second << " bytes " << endl; 
     262        report(10)<< " Memory report : Context <"<<context->getId()<<"> : server side : memory used for buffer of each connection to client" << endl 
     263                  << "  +) With client of rank " << itMap->first << " : " << itMap->second << " bytes " << endl; 
    231264        totalBuf += itMap->second; 
    232265      } 
    233266      context->finalize(); 
    234       //report(0)<< " Memory report : Context <"<<context->getId()<<"> : server side : total memory used for buffer "<<totalBuf<<" bytes"<<endl; 
     267      report(0)<< " Memory report : Context <"<<context->getId()<<"> : server side : total memory used for buffer "<<totalBuf<<" bytes"<<endl; 
    235268    } 
    236269    else if (event.classId==CContext::GetType()) CContext::dispatchEvent(event); 
  • XIOS/dev/branch_openmp/src/context_server.hpp

    r1134 r1328  
    1717    bool eventLoop(bool enableEventsProcessing = true); 
    1818    void listen(void) ; 
     19    bool listenPendingRequest(ep_lib::MPI_Status& status); 
    1920    void checkPendingRequest(void) ; 
    2021    void processRequest(int rank, char* buff,int count) ; 
  • XIOS/dev/branch_openmp/src/cxios.cpp

    r1205 r1328  
    1111#include "memtrack.hpp" 
    1212#include "registry.hpp" 
     13using namespace ep_lib; 
    1314 
    1415namespace xios 
    1516{ 
    16  
    17   extern int test_omp_rank; 
    18   #pragma omp threadprivate(test_omp_rank) 
    19  
    20   const string CXios::rootFile="./iodef.xml" ; 
    21   const string CXios::xiosCodeId="xios.x" ; 
    22   const string CXios::clientFile="./xios_client"; 
    23   const string CXios::serverFile="./xios_server"; 
    24  
     17  string CXios::rootFile="./iodef.xml" ; 
     18  string CXios::xiosCodeId="xios.x" ; 
     19  string CXios::clientFile="./xios_client"; 
     20  string CXios::serverFile="./xios_server"; 
    2521 
    2622  bool CXios::isClient ; 
    2723  bool CXios::isServer ; 
    28  
    29  
    3024  MPI_Comm CXios::globalComm ; 
    31  
    32    
    3325  bool CXios::usingOasis ; 
    3426  bool CXios::usingServer = false; 
    35  
    36  
    3727  double CXios::bufferSizeFactor = 1.0; 
    3828  const double CXios::defaultBufferSizeFactor = 1.0; 
    3929  StdSize CXios::minBufferSize = 1024 * sizeof(double); 
    40  
    41  
    4230  bool CXios::printLogs2Files; 
    4331  bool CXios::isOptPerformance = true; 
     
    4937  { 
    5038    set_new_handler(noMemory); 
    51      
    52      
    53     #pragma omp critical 
    54     { 
    55       parseFile(rootFile);   
    56     } 
    57     #pragma omp barrier 
     39    parseFile(rootFile); 
    5840    parseXiosConfig(); 
    5941  } 
     
    8769      ERROR("CXios::parseXiosConfig()", "recv_field_timeout cannot be negative."); 
    8870 
    89   
     71    //globalComm=MPI_COMM_WORLD ; 
    9072    int num_ep; 
    9173    if(isClient)   
     
    9375      num_ep = omp_get_num_threads(); 
    9476    } 
    95      
     77         
    9678    if(isServer)  
    9779    {  
    9880      num_ep = omp_get_num_threads(); 
    9981    } 
    100      
     82         
    10183    MPI_Info info; 
    10284    #pragma omp master 
     
    10688      passage = ep_comm;   
    10789    } 
    108      
     90         
    10991    #pragma omp barrier 
    110  
    111        
     92     
     93           
    11294    CXios::globalComm = passage[omp_get_thread_num()]; 
    113  
    114     int tmp_rank; 
    115     MPI_Comm_rank(CXios::globalComm, &tmp_rank); 
    116  
    117      
    118     test_omp_rank = tmp_rank; 
    119      
    12095  } 
    12196 
     
    133108 
    134109    CClient::initialize(codeId,localComm,returnComm) ; 
    135  
    136110    if (CClient::getRank()==0) globalRegistry = new CRegistry(returnComm) ; 
    137111 
     
    142116    if (printLogs2Files) 
    143117    { 
    144       #pragma omp critical 
    145118      CClient::openInfoStream(clientFile); 
    146119      CClient::openErrorStream(clientFile); 
     
    158131     if (CClient::getRank()==0) 
    159132     { 
    160        #pragma omp critical (_output) 
    161133       info(80)<<"Write data base Registry"<<endl<<globalRegistry->toString()<<endl ; 
    162134       globalRegistry->toFile("xios_registry.bin") ; 
     
    184156  void CXios::initServer() 
    185157  { 
    186     int initialized; 
    187     MPI_Initialized(&initialized); 
    188     if (initialized) CServer::is_MPI_Initialized=true ; 
    189     else CServer::is_MPI_Initialized=false ; 
    190        
    191   
    192     if(!CServer::is_MPI_Initialized) 
    193     { 
    194       MPI_Init(NULL, NULL); 
    195     } 
    196        
    197158    set_new_handler(noMemory); 
    198159    std::set<StdString> parseList; 
     
    210171     
    211172    initServer(); 
    212      
    213      
     173 
    214174    // Initialize all aspects MPI 
    215175    CServer::initialize(); 
  • XIOS/dev/branch_openmp/src/cxios.hpp

    r1134 r1328  
    55#include "mpi.hpp" 
    66#include "registry.hpp" 
    7 #include "log.hpp" 
    87 
    98namespace xios 
     
    1514  { 
    1615    public: 
    17       static void initialize(void) ; 
    18       static void initClientSide(const string & codeId, ep_lib::MPI_Comm& localComm, ep_lib::MPI_Comm& returnComm) ; 
    19       static void initServerSide(void) ; 
    20       static void clientFinalize(void) ; 
    21       static void parseFile(const string& filename) ; 
     16     static void initialize(void) ; 
     17     static void initClientSide(const string & codeId, ep_lib::MPI_Comm& localComm, ep_lib::MPI_Comm& returnComm) ; 
     18     static void initServerSide(void) ; 
     19     static void clientFinalize(void) ; 
     20     static void parseFile(const string& filename) ; 
    2221 
    23       template <typename T> 
    24       static T getin(const string& id,const T& defaultValue) ; 
     22     template <typename T> 
     23     static T getin(const string& id,const T& defaultValue) ; 
    2524 
    26       template <typename T> 
    27       static T getin(const string& id) ; 
     25     template <typename T> 
     26     static T getin(const string& id) ; 
    2827 
    2928    public: 
    30       static const string rootFile; //!< Configuration filename 
    31       static const string xiosCodeId ; //!< Identity for XIOS 
    32       static const string clientFile; //!< Filename template for client 
    33       static const string serverFile; //!< Filename template for server 
     29     static string rootFile ; //!< Configuration filename 
     30     static string xiosCodeId ; //!< Identity for XIOS 
     31     static string clientFile; //!< Filename template for client 
     32     static string serverFile; //!< Filename template for server 
    3433 
    35       static bool isClient ; //!< Check if xios is client 
    36       static bool isServer ; //!< Check if xios is server 
    37       #pragma omp threadprivate(isClient, isServer) 
     34     static bool isClient ; //!< Check if xios is client 
     35     static bool isServer ; //!< Check if xios is server 
    3836 
    39       static MPI_Comm globalComm ; //!< Global communicator 
    40       #pragma omp threadprivate(globalComm) 
     37     static ep_lib::MPI_Comm globalComm ; //!< Global communicator 
    4138 
    42       static bool printLogs2Files; //!< Printing out logs into files 
    43       static bool usingOasis ; //!< Using Oasis 
    44       static bool usingServer ; //!< Using server (server mode) 
    45       static double bufferSizeFactor; //!< Factor used to tune the buffer size 
    46       static const double defaultBufferSizeFactor; //!< Default factor value 
    47       static StdSize minBufferSize; //!< Minimum buffer size 
    48       static bool isOptPerformance; //!< Check if buffer size is for performance (as large as possible) 
    49       #pragma omp threadprivate(printLogs2Files, usingOasis, usingServer, bufferSizeFactor, minBufferSize, isOptPerformance) 
    50        
    51       static CRegistry* globalRegistry ; //!< global registry which is wrote by the root process of the servers 
    52       static double recvFieldTimeout; 
    53       #pragma omp threadprivate(recvFieldTimeout) 
    54        
    55        
     39     static bool printLogs2Files; //!< Printing out logs into files 
     40     static bool usingOasis ; //!< Using Oasis 
     41     static bool usingServer ; //!< Using server (server mode) 
     42     static double bufferSizeFactor; //!< Factor used to tune the buffer size 
     43     static const double defaultBufferSizeFactor; //!< Default factor value 
     44     static StdSize minBufferSize; //!< Minimum buffer size 
     45     static bool isOptPerformance; //!< Check if buffer size is for performance (as large as possible) 
     46     static CRegistry* globalRegistry ; //!< global registry which is wrote by the root process of the servers 
     47     static double recvFieldTimeout; //!< Time to wait for data before issuing an error when receiving a field 
     48 
    5649    public: 
    5750     //! Setting xios to use server mode 
  • XIOS/dev/branch_openmp/src/data_output.cpp

    r1205 r1328  
    44#include "group_template.hpp" 
    55#include "context.hpp" 
    6 //mpi.hpp 
     6 
    77namespace xios 
    88{ 
  • XIOS/dev/branch_openmp/src/data_output.hpp

    r1287 r1328