Ignore:
Timestamp:
05/18/17 17:40:03 (7 years ago)
Author:
yushan
Message:

test_remap back to work. No thread for now

Location:
XIOS/dev/branch_yushan_merged/extern
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • XIOS/dev/branch_yushan_merged/extern/remap/src/mapper.cpp

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

  • XIOS/dev/branch_yushan_merged/extern/src_ep_dev/ep_gatherv.cpp

    r1134 r1138  
    484484    mpi_size = comm.ep_comm_ptr->size_rank_info[2].second; 
    485485     
    486  
    487     assert(accumulate(recvcounts, recvcounts+ep_size-1, 0) == displs[ep_size-1]); // Only for continuous gather. 
     486    if(ep_size == mpi_size)  
     487      return ::MPI_Allgatherv(sendbuf, sendcount, static_cast< ::MPI_Datatype>(datatype), recvbuf, recvcounts, displs, 
     488                              static_cast< ::MPI_Datatype>(datatype), static_cast< ::MPI_Comm>(comm.mpi_comm)); 
     489     
     490 
     491    assert(accumulate(recvcounts, recvcounts+ep_size-1, 0) >= displs[ep_size-1]); // Only for continuous gather. 
    488492 
    489493 

  • XIOS/dev/branch_yushan_merged/extern/src_ep_dev/ep_memory.cpp

    r1134 r1138  
    1010  int MPI_Alloc_mem(MPI_Aint size, MPI_Info info, void *baseptr) 
    1111  { 
    12     ::MPI_Alloc_mem(size.mpi_aint, static_cast< ::MPI_Info>(info.mpi_info), baseptr); 
     12    //::MPI_Alloc_mem(size.mpi_aint, static_cast< ::MPI_Info>(info.mpi_info), baseptr); 
     13    ::MPI_Alloc_mem(size.mpi_aint, MPI_INFO_NULL_STD, baseptr); 
    1314    return 0; 
    1415   } 
    … …  
    1617  int MPI_Alloc_mem(unsigned long size, MPI_Info info, void *baseptr) 
    1718  { 
    18     ::MPI_Alloc_mem(size, static_cast< ::MPI_Info>(info.mpi_info), baseptr); 
     19    //::MPI_Alloc_mem(size, static_cast< ::MPI_Info>(info.mpi_info), baseptr); 
     20    ::MPI_Alloc_mem(size, MPI_INFO_NULL_STD, baseptr); 
    1921    return 0; 
    2022  } 

  • XIOS/dev/branch_yushan_merged/extern/src_ep_dev/ep_wait.cpp

    r1134 r1138  
    11/*! 
    2    \file ep_wait.cpp 
    3    \since 2 may 2016 
     2  \file ep_wait.cpp 
     3  \since 2 may 2016 
    44 
    5    \brief Definitions of MPI wait function: MPI_Wait, MPI_Waitall, MPI_Waitsome, MPI_Waitany 
    6  */ 
     5  \brief Definitions of MPI wait function: MPI_Wait, MPI_Waitall, MPI_Waitsome, MPI_Waitany 
     6  */ 
    77 
    88#include "ep_lib.hpp" 
    … …  
    1616namespace ep_lib {       
    1717 
    18         int MPI_Wait(MPI_Request *request, MPI_Status *status) 
    19         { 
     18    int MPI_Wait(MPI_Request *request, MPI_Status *status) 
     19    { 
    2020 
    21     if(request->type == 1) 
    22     { 
    23       ::MPI_Request mpi_request = static_cast< ::MPI_Request >(request->mpi_request); 
    24       ::MPI_Status mpi_status; 
    25       ::MPI_Wait(&mpi_request, &mpi_status); 
    26  
    27  
    28        
    29       status->mpi_status = &mpi_status; 
    30       status->ep_src = request->ep_src; 
    31       status->ep_tag = request->ep_tag; 
    32       status->ep_datatype = request->ep_datatype; 
    33  
    34       return 0; 
    35     } 
    36  
    37     if(request->type == 2) 
    38     { 
    39       int flag = false; 
    40       MPI_Message message; 
    41  
    42       while(!flag) 
    43       { 
    44         Message_Check(request->comm); 
    45         #pragma omp flush 
    46         MPI_Improbe(request->ep_src, request->ep_tag, request->comm, &flag, &message, status); 
    47       } 
    48  
    49       int count; 
    50       MPI_Get_count(status, request->ep_datatype, &count); 
    51       MPI_Mrecv(request->buf, count, request->ep_datatype, &message, status); 
    52       status->ep_datatype = request->ep_datatype; 
    53  
    54       //check_sum_recv(request->buf, count, request->ep_datatype, request->ep_src, request->ep_tag, request->comm, 2); 
    55  
    56       return 0; 
    57     } 
    58  
    59     if(request->type == 3) 
    60     { 
    61       ::MPI_Request mpi_request = static_cast< ::MPI_Request >(request->mpi_request); 
    62       ::MPI_Status mpi_status; 
    63       ::MPI_Wait(&mpi_request, &mpi_status); 
    64        
    65       status->mpi_status = new ::MPI_Status(mpi_status); 
    66       status->ep_src = request->ep_src; 
    67       status->ep_tag = request->ep_tag; 
    68       status->ep_datatype = request->ep_datatype; 
    69  
    70       //int count; 
    71       //MPI_Get_count(status, request->ep_datatype, &count); 
    72       //check_sum_recv(request->buf, count, request->ep_datatype, request->ep_src, request->ep_tag, request->comm, 2); 
    73     } 
    74           return MPI_SUCCESS; 
    75         } 
     21        if(request->type == 1) 
     22        { 
     23            ::MPI_Request mpi_request = static_cast< ::MPI_Request >(request->mpi_request); 
     24            ::MPI_Status mpi_status; 
     25            ::MPI_Wait(&mpi_request, &mpi_status); 
    7626 
    7727 
    7828 
    79          
     29            status->mpi_status = &mpi_status; 
     30            status->ep_src = request->ep_src; 
     31            status->ep_tag = request->ep_tag; 
     32            status->ep_datatype = request->ep_datatype; 
     33 
     34            return 0; 
     35        } 
     36 
     37        if(request->type == 2) 
     38        { 
     39            int flag = false; 
     40            MPI_Message message; 
     41 
     42            while(!flag) 
     43            { 
     44                Message_Check(request->comm); 
     45#pragma omp flush 
     46                MPI_Improbe(request->ep_src, request->ep_tag, request->comm, &flag, &message, status); 
     47            } 
     48 
     49            int count; 
     50            MPI_Get_count(status, request->ep_datatype, &count); 
     51            MPI_Mrecv(request->buf, count, request->ep_datatype, &message, status); 
     52            status->ep_datatype = request->ep_datatype; 
     53 
     54            //check_sum_recv(request->buf, count, request->ep_datatype, request->ep_src, request->ep_tag, request->comm, 2); 
     55 
     56            return 0; 
     57        } 
     58 
     59        if(request->type == 3) 
     60        { 
     61            ::MPI_Request mpi_request = static_cast< ::MPI_Request >(request->mpi_request); 
     62            ::MPI_Status mpi_status; 
     63            ::MPI_Wait(&mpi_request, &mpi_status); 
     64 
     65            status->mpi_status = new ::MPI_Status(mpi_status); 
     66            status->ep_src = request->ep_src; 
     67            status->ep_tag = request->ep_tag; 
     68            status->ep_datatype = request->ep_datatype; 
     69 
     70            //int count; 
     71            //MPI_Get_count(status, request->ep_datatype, &count); 
     72            //check_sum_recv(request->buf, count, request->ep_datatype, request->ep_src, request->ep_tag, request->comm, 2); 
     73        } 
     74        return MPI_SUCCESS; 
     75    } 
    8076 
    8177 
    82         int MPI_Waitall(int count, MPI_Request *array_of_requests, MPI_Status *array_of_statuses) 
    83         { 
    84     int dest_rank; 
    85     MPI_Comm_rank(MPI_COMM_WORLD, &dest_rank); 
    86      
    87     int finished = 0; 
    88     bool finished_index[count]; 
    8978 
    90           for(int i=0; i<count; i++) 
     79 
     80 
     81 
     82    int MPI_Waitall(int count, MPI_Request *array_of_requests, MPI_Status *array_of_statuses) 
    9183    { 
    92       finished_index[i] = false; 
     84        int dest_rank; 
     85        MPI_Comm_rank(MPI_COMM_WORLD, &dest_rank); 
     86 
     87        int finished = 0; 
     88        bool finished_index[count]; 
     89 
     90        for(int i=0; i<count; i++) 
     91        { 
     92            finished_index[i] = false; 
     93        } 
     94 
     95        while(finished < count) 
     96        { 
     97            for(int i=0; i<count; i++) 
     98            { 
     99                if(finished_index[i] == false) // this request has not been tested. 
     100                { 
     101                    if(array_of_requests[i].type != 2) // isend or imrecv 
     102                    {       
     103                        MPI_Wait(&array_of_requests[i], &array_of_statuses[i]); 
     104                        if(array_of_requests[i].type == 3) 
     105                        { 
     106                            //int check_count; 
     107                            //MPI_Get_count(&array_of_statuses[i], array_of_requests[i].ep_datatype, &check_count); 
     108                            //check_sum_recv(array_of_requests[i].buf, count, array_of_requests[i].ep_datatype, array_of_requests[i].ep_src, array_of_requests[i].ep_tag, array_of_requests[i].comm, 2); 
     109                        } 
     110                        finished++; 
     111                        finished_index[i] = true; 
     112                    } 
     113                    else // irecv 
     114                    { 
     115                        int flag = false; 
     116                        MPI_Message message; 
     117 
     118                        MPI_Improbe(array_of_requests[i].ep_src, array_of_requests[i].ep_tag, array_of_requests[i].comm, &flag, &message, &array_of_statuses[i]); 
     119 
     120                        if(flag) 
     121                        { 
     122                            //printf("dest_rank = %d, Waiting one message with src = %d, tag = %d, buf = %p\n", dest_rank, array_of_requests[i].ep_src, array_of_requests[i].ep_tag, array_of_requests[i].buf); 
     123                            int recv_count; 
     124                            MPI_Get_count(&array_of_statuses[i], array_of_requests[i].ep_datatype, &recv_count); 
     125                            MPI_Mrecv(array_of_requests[i].buf, recv_count, array_of_requests[i].ep_datatype, &message, &array_of_statuses[i]); 
     126                            //check_sum_recv(array_of_requests[i].buf, recv_count, array_of_requests[i].ep_datatype, array_of_requests[i].ep_src, array_of_requests[i].ep_tag, array_of_requests[i].comm, 2); 
     127 
     128                            finished++; 
     129                            finished_index[i] = true; 
     130                        } 
     131                    } 
     132                } 
     133            }     
     134        } 
     135        return MPI_SUCCESS; 
    93136    } 
    94  
    95     while(finished < count) 
    96     { 
    97       for(int i=0; i<count; i++) 
    98       { 
    99         if(finished_index[i] == false) // this request has not been tested. 
    100         { 
    101           if(array_of_requests[i].type != 2) // isend or imrecv 
    102           {       
    103             MPI_Wait(&array_of_requests[i], &array_of_statuses[i]); 
    104             if(array_of_requests[i].type == 3) 
    105             { 
    106               //int check_count; 
    107               //MPI_Get_count(&array_of_statuses[i], array_of_requests[i].ep_datatype, &check_count); 
    108               //check_sum_recv(array_of_requests[i].buf, count, array_of_requests[i].ep_datatype, array_of_requests[i].ep_src, array_of_requests[i].ep_tag, array_of_requests[i].comm, 2); 
    109             } 
    110             finished++; 
    111             finished_index[i] = true; 
    112           } 
    113           else // irecv 
    114           { 
    115             int flag = false; 
    116             MPI_Message message; 
    117              
    118             MPI_Improbe(array_of_requests[i].ep_src, array_of_requests[i].ep_tag, array_of_requests[i].comm, &flag, &message, &array_of_statuses[i]); 
    119              
    120             if(flag) 
    121             { 
    122               //printf("dest_rank = %d, Waiting one message with src = %d, tag = %d, buf = %p\n", dest_rank, array_of_requests[i].ep_src, array_of_requests[i].ep_tag, array_of_requests[i].buf); 
    123               int recv_count; 
    124               MPI_Get_count(&array_of_statuses[i], array_of_requests[i].ep_datatype, &recv_count); 
    125               MPI_Mrecv(array_of_requests[i].buf, recv_count, array_of_requests[i].ep_datatype, &message, &array_of_statuses[i]); 
    126               //check_sum_recv(array_of_requests[i].buf, recv_count, array_of_requests[i].ep_datatype, array_of_requests[i].ep_src, array_of_requests[i].ep_tag, array_of_requests[i].comm, 2); 
    127  
    128               finished++; 
    129               finished_index[i] = true; 
    130             } 
    131           } 
    132         } 
    133       }     
    134     } 
    135           return MPI_SUCCESS; 
    136         } 
    137137 
    138138 
Note: See TracChangeset for help on using the changeset viewer.