1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 3 /* TODO PetscArrayExchangeBegin/End */ 4 /* TODO blocksize */ 5 /* TODO move to API ? */ 6 static PetscErrorCode ExchangeArrayByRank_Private(PetscObject obj, MPI_Datatype dt, PetscInt nsranks, const PetscMPIInt sranks[], PetscInt ssize[], const void *sarr[], PetscInt nrranks, const PetscMPIInt rranks[], PetscInt *rsize_out[], void **rarr_out[]) 7 { 8 PetscInt r; 9 PetscInt *rsize; 10 void **rarr; 11 MPI_Request *sreq, *rreq; 12 PetscMPIInt tag, unitsize; 13 MPI_Comm comm; 14 PetscErrorCode ierr; 15 16 PetscFunctionBegin; 17 ierr = MPI_Type_size(dt, &unitsize);CHKERRQ(ierr); 18 ierr = PetscObjectGetComm(obj, &comm);CHKERRQ(ierr); 19 ierr = PetscMalloc2(nrranks, &rsize, nrranks, &rarr);CHKERRQ(ierr); 20 ierr = PetscMalloc2(nrranks, &rreq, nsranks, &sreq);CHKERRQ(ierr); 21 /* exchange array size */ 22 ierr = PetscObjectGetNewTag(obj,&tag);CHKERRQ(ierr); 23 for (r=0; r<nrranks; r++) { 24 ierr = MPI_Irecv(&rsize[r], 1, MPIU_INT, rranks[r], tag, comm, &rreq[r]);CHKERRQ(ierr); 25 } 26 for (r=0; r<nsranks; r++) { 27 ierr = MPI_Isend(&ssize[r], 1, MPIU_INT, sranks[r], tag, comm, &sreq[r]);CHKERRQ(ierr); 28 } 29 ierr = MPI_Waitall(nrranks, rreq, MPI_STATUSES_IGNORE);CHKERRQ(ierr); 30 /* exchange array */ 31 ierr = PetscObjectGetNewTag(obj,&tag);CHKERRQ(ierr); 32 for (r=0; r<nrranks; r++) { 33 ierr = PetscMalloc(rsize[r]*unitsize, &rarr[r]);CHKERRQ(ierr); 34 ierr = MPI_Irecv(rarr[r], rsize[r], dt, rranks[r], tag, comm, &rreq[r]);CHKERRQ(ierr); 35 } 36 for (r=0; r<nsranks; r++) { 37 ierr = MPI_Isend(sarr[r], ssize[r], dt, sranks[r], tag, comm, &sreq[r]);CHKERRQ(ierr); 38 } 39 ierr = MPI_Waitall(nrranks, rreq, MPI_STATUSES_IGNORE);CHKERRQ(ierr); 40 ierr = MPI_Waitall(nsranks, sreq, MPI_STATUSES_IGNORE);CHKERRQ(ierr); 41 ierr = PetscFree2(rreq, sreq);CHKERRQ(ierr); 42 *rsize_out = rsize; 43 *rarr_out = rarr; 44 PetscFunctionReturn(0); 45 } 46 47 /* TODO VecExchangeBegin/End */ 48 /* TODO move to API ? */ 49 static PetscErrorCode ExchangeVecByRank_Private(PetscObject obj, PetscInt nsranks, const PetscMPIInt sranks[], Vec svecs[], PetscInt nrranks, const PetscMPIInt rranks[], Vec *rvecs[]) 50 { 51 PetscInt r; 52 PetscInt *ssize, *rsize; 53 PetscScalar **rarr; 54 const PetscScalar **sarr; 55 Vec *rvecs_; 56 MPI_Request *sreq, *rreq; 57 PetscErrorCode ierr; 58 59 PetscFunctionBegin; 60 ierr = PetscMalloc4(nsranks, &ssize, nsranks, &sarr, nrranks, &rreq, nsranks, &sreq);CHKERRQ(ierr); 61 for (r=0; r<nsranks; r++) { 62 ierr = VecGetLocalSize(svecs[r], &ssize[r]);CHKERRQ(ierr); 63 ierr = VecGetArrayRead(svecs[r], &sarr[r]);CHKERRQ(ierr); 64 } 65 ierr = ExchangeArrayByRank_Private(obj, MPIU_SCALAR, nsranks, sranks, ssize, (const void**)sarr, nrranks, rranks, &rsize, (void***)&rarr);CHKERRQ(ierr); 66 ierr = PetscMalloc1(nrranks, &rvecs_);CHKERRQ(ierr); 67 for (r=0; r<nrranks; r++) { 68 /* set array in two steps to mimic PETSC_OWN_POINTER */ 69 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF, 1, rsize[r], NULL, &rvecs_[r]);CHKERRQ(ierr); 70 ierr = VecReplaceArray(rvecs_[r], rarr[r]);CHKERRQ(ierr); 71 } 72 for (r=0; r<nsranks; r++) { 73 ierr = VecRestoreArrayRead(svecs[r], &sarr[r]);CHKERRQ(ierr); 74 } 75 ierr = PetscFree2(rsize, rarr);CHKERRQ(ierr); 76 ierr = PetscFree4(ssize, sarr, rreq, sreq);CHKERRQ(ierr); 77 *rvecs = rvecs_; 78 PetscFunctionReturn(0); 79 } 80 81 static PetscErrorCode SortByRemote_Private(PetscSF sf, PetscInt *rmine1[], PetscInt *rremote1[]) 82 { 83 PetscInt nleaves; 84 PetscInt nranks; 85 const PetscMPIInt *ranks; 86 const PetscInt *roffset, *rmine, *rremote; 87 PetscInt n, o, r; 88 PetscErrorCode ierr; 89 90 PetscFunctionBegin; 91 ierr = PetscSFGetRanks(sf, &nranks, &ranks, &roffset, &rmine, &rremote);CHKERRQ(ierr); 92 nleaves = roffset[nranks]; 93 ierr = PetscMalloc2(nleaves, rmine1, nleaves, rremote1);CHKERRQ(ierr); 94 for (r=0; r<nranks; r++) { 95 /* simultaneously sort rank-wise portions of rmine & rremote by values in rremote 96 - to unify order with the other side */ 97 o = roffset[r]; 98 n = roffset[r+1] - o; 99 ierr = PetscMemcpy(&(*rmine1)[o], &rmine[o], n*sizeof(PetscInt));CHKERRQ(ierr); 100 ierr = PetscMemcpy(&(*rremote1)[o], &rremote[o], n*sizeof(PetscInt));CHKERRQ(ierr); 101 ierr = PetscSortIntWithArray(n, &(*rremote1)[o], &(*rmine1)[o]);CHKERRQ(ierr); 102 } 103 PetscFunctionReturn(0); 104 } 105 106 static PetscErrorCode GetRecursiveConeCoordinatesPerRank_Private(DM dm, PetscSF sf, PetscInt rmine[], Vec *coordinatesPerRank[]) 107 { 108 IS pointsPerRank, conesPerRank; 109 PetscInt nranks; 110 const PetscMPIInt *ranks; 111 const PetscInt *roffset; 112 PetscInt n, o, r; 113 PetscErrorCode ierr; 114 115 PetscFunctionBegin; 116 ierr = DMGetCoordinatesLocalSetUp(dm);CHKERRQ(ierr); 117 ierr = PetscSFGetRanks(sf, &nranks, &ranks, &roffset, NULL, NULL);CHKERRQ(ierr); 118 ierr = PetscMalloc1(nranks, coordinatesPerRank);CHKERRQ(ierr); 119 for (r=0; r<nranks; r++) { 120 o = roffset[r]; 121 n = roffset[r+1] - o; 122 ierr = ISCreateGeneral(PETSC_COMM_SELF, n, &rmine[o], PETSC_USE_POINTER, &pointsPerRank);CHKERRQ(ierr); 123 ierr = DMPlexGetConeRecursive(dm, pointsPerRank, &conesPerRank);CHKERRQ(ierr); 124 ierr = DMGetCoordinatesLocalTuple(dm, conesPerRank, NULL, &(*coordinatesPerRank)[r]);CHKERRQ(ierr); 125 ierr = ISDestroy(&pointsPerRank);CHKERRQ(ierr); 126 ierr = ISDestroy(&conesPerRank);CHKERRQ(ierr); 127 } 128 PetscFunctionReturn(0); 129 } 130 131 static PetscErrorCode PetscSFComputeMultiRootOriginalNumberingByRank_Private(PetscSF sf, PetscSF imsf, PetscInt *irmine1[]) 132 { 133 PetscInt *mRootsOrigNumbering; 134 PetscInt nileaves, niranks; 135 const PetscInt *iroffset, *irmine, *degree; 136 PetscInt i, n, o, r; 137 PetscErrorCode ierr; 138 139 PetscFunctionBegin; 140 ierr = PetscSFGetGraph(imsf, NULL, &nileaves, NULL, NULL);CHKERRQ(ierr); 141 ierr = PetscSFGetRanks(imsf, &niranks, NULL, &iroffset, &irmine, NULL);CHKERRQ(ierr); 142 #if defined(PETSC_USE_DEBUG) 143 if (PetscUnlikely(nileaves != iroffset[niranks])) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"nileaves != iroffset[niranks])"); 144 #endif 145 ierr = PetscSFComputeDegreeBegin(sf, °ree);CHKERRQ(ierr); 146 ierr = PetscSFComputeDegreeEnd(sf, °ree);CHKERRQ(ierr); 147 ierr = PetscSFComputeMultiRootOriginalNumbering(sf, degree, &mRootsOrigNumbering);CHKERRQ(ierr); 148 ierr = PetscMalloc1(nileaves, irmine1);CHKERRQ(ierr); 149 for (r=0; r<niranks; r++) { 150 o = iroffset[r]; 151 n = iroffset[r+1] - o; 152 for (i=0; i<n; i++) (*irmine1)[o+i] = mRootsOrigNumbering[irmine[o+i]]; 153 } 154 ierr = PetscFree(mRootsOrigNumbering);CHKERRQ(ierr); 155 PetscFunctionReturn(0); 156 } 157 158 /* TODO documentation */ 159 PetscErrorCode DMPlexCheckConesConformOnInterfaces(DM dm) 160 { 161 PetscSF sf; 162 PetscInt nleaves, nranks, nroots; 163 const PetscInt *mine, *roffset, *rmine, *rremote; 164 const PetscSFNode *remote; 165 const PetscMPIInt *ranks; 166 PetscSF msf, imsf; 167 PetscInt nileaves, niranks; 168 const PetscMPIInt *iranks; 169 const PetscInt *iroffset, *irmine, *irremote; 170 PetscInt *rmine1, *rremote1; /* rmine and rremote copies simultaneously sorted by rank and rremote */ 171 PetscInt *mine_orig_numbering; 172 Vec *sntCoordinatesPerRank; 173 Vec *refCoordinatesPerRank; 174 Vec *recCoordinatesPerRank; 175 PetscInt r; 176 PetscMPIInt commsize, myrank; 177 PetscBool same; 178 MPI_Comm comm; 179 PetscErrorCode ierr; 180 181 PetscFunctionBegin; 182 ierr = PetscObjectGetComm((PetscObject)dm, &comm);CHKERRQ(ierr); 183 ierr = MPI_Comm_rank(comm, &myrank);CHKERRQ(ierr); 184 ierr = MPI_Comm_size(comm, &commsize);CHKERRQ(ierr); 185 if (commsize < 2) PetscFunctionReturn(0); 186 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr); 187 if (!sf) PetscFunctionReturn(0); 188 ierr = PetscSFGetGraph(sf, &nroots, &nleaves, &mine, &remote);CHKERRQ(ierr); 189 if (nroots < 0) PetscFunctionReturn(0); 190 ierr = PetscSFSetUp(sf);CHKERRQ(ierr); 191 ierr = PetscSFGetRanks(sf, &nranks, &ranks, &roffset, &rmine, &rremote);CHKERRQ(ierr); 192 193 /* Expand sent cones per rank */ 194 ierr = SortByRemote_Private(sf, &rmine1, &rremote1);CHKERRQ(ierr); 195 ierr = GetRecursiveConeCoordinatesPerRank_Private(dm, sf, rmine1, &sntCoordinatesPerRank);CHKERRQ(ierr); 196 197 /* Create inverse SF */ 198 ierr = PetscSFGetMultiSF(sf,&msf);CHKERRQ(ierr); 199 ierr = PetscSFCreateInverseSF(msf,&imsf);CHKERRQ(ierr); 200 ierr = PetscSFSetUp(imsf);CHKERRQ(ierr); 201 ierr = PetscSFGetGraph(imsf, NULL, &nileaves, NULL, NULL);CHKERRQ(ierr); 202 ierr = PetscSFGetRanks(imsf, &niranks, &iranks, &iroffset, &irmine, &irremote);CHKERRQ(ierr); 203 204 /* Compute original numbering of multi-roots (referenced points) */ 205 ierr = PetscSFComputeMultiRootOriginalNumberingByRank_Private(sf, imsf, &mine_orig_numbering);CHKERRQ(ierr); 206 207 /* Expand referred cones per rank */ 208 ierr = GetRecursiveConeCoordinatesPerRank_Private(dm, imsf, mine_orig_numbering, &refCoordinatesPerRank);CHKERRQ(ierr); 209 210 /* Send the coordinates */ 211 ierr = ExchangeVecByRank_Private((PetscObject)sf, nranks, ranks, sntCoordinatesPerRank, niranks, iranks, &recCoordinatesPerRank);CHKERRQ(ierr); 212 213 /* Compare recCoordinatesPerRank with refCoordinatesPerRank */ 214 for (r=0; r<niranks; r++) { 215 ierr = VecEqual(refCoordinatesPerRank[r], recCoordinatesPerRank[r], &same);CHKERRQ(ierr); 216 if (!same) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "interface cones do not conform for remote rank %d", iranks[r]); 217 } 218 219 /* destroy sent stuff */ 220 for (r=0; r<nranks; r++) { 221 ierr = VecDestroy(&sntCoordinatesPerRank[r]);CHKERRQ(ierr); 222 } 223 ierr = PetscFree(sntCoordinatesPerRank);CHKERRQ(ierr); 224 ierr = PetscFree2(rmine1, rremote1);CHKERRQ(ierr); 225 ierr = PetscSFDestroy(&imsf);CHKERRQ(ierr); 226 227 /* destroy referenced stuff */ 228 for (r=0; r<niranks; r++) { 229 ierr = VecDestroy(&refCoordinatesPerRank[r]);CHKERRQ(ierr); 230 } 231 ierr = PetscFree(refCoordinatesPerRank);CHKERRQ(ierr); 232 ierr = PetscFree(mine_orig_numbering);CHKERRQ(ierr); 233 234 /* destroy received stuff */ 235 for (r=0; r<niranks; r++) { 236 ierr = VecDestroy(&recCoordinatesPerRank[r]);CHKERRQ(ierr); 237 } 238 ierr = PetscFree(recCoordinatesPerRank);CHKERRQ(ierr); 239 PetscFunctionReturn(0); 240 } 241