1 #define PETSCDM_DLL 2 #include <petsc/private/dmnetworkimpl.h> /*I "petscdmnetwork.h" I*/ 3 #include <petsc/private/vecimpl.h> 4 5 PetscErrorCode DMSetFromOptions_Network(PetscOptionItems *PetscOptionsObject,DM dm) 6 { 7 PetscErrorCode ierr; 8 9 PetscFunctionBegin; 10 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 11 ierr = PetscOptionsHead(PetscOptionsObject,"DMNetwork Options");CHKERRQ(ierr); 12 ierr = PetscOptionsTail();CHKERRQ(ierr); 13 PetscFunctionReturn(0); 14 } 15 16 /* External function declarations here */ 17 extern PetscErrorCode DMCreateMatrix_Network(DM, Mat*); 18 extern PetscErrorCode DMDestroy_Network(DM); 19 extern PetscErrorCode DMView_Network(DM, PetscViewer); 20 extern PetscErrorCode DMGlobalToLocalBegin_Network(DM, Vec, InsertMode, Vec); 21 extern PetscErrorCode DMGlobalToLocalEnd_Network(DM, Vec, InsertMode, Vec); 22 extern PetscErrorCode DMLocalToGlobalBegin_Network(DM, Vec, InsertMode, Vec); 23 extern PetscErrorCode DMLocalToGlobalEnd_Network(DM, Vec, InsertMode, Vec); 24 extern PetscErrorCode DMSetUp_Network(DM); 25 extern PetscErrorCode DMClone_Network(DM, DM*); 26 27 static PetscErrorCode VecArrayPrint_private(PetscViewer viewer,PetscInt n,const PetscScalar *xv) 28 { 29 PetscErrorCode ierr; 30 PetscInt i; 31 32 PetscFunctionBegin; 33 for (i=0; i<n; i++) { 34 #if defined(PETSC_USE_COMPLEX) 35 if (PetscImaginaryPart(xv[i]) > 0.0) { 36 ierr = PetscViewerASCIIPrintf(viewer," %g + %g i\n",(double)PetscRealPart(xv[i]),(double)PetscImaginaryPart(xv[i]));CHKERRQ(ierr); 37 } else if (PetscImaginaryPart(xv[i]) < 0.0) { 38 ierr = PetscViewerASCIIPrintf(viewer," %g - %g i\n",(double)PetscRealPart(xv[i]),-(double)PetscImaginaryPart(xv[i]));CHKERRQ(ierr); 39 } else { 40 ierr = PetscViewerASCIIPrintf(viewer," %g\n",(double)PetscRealPart(xv[i]));CHKERRQ(ierr); 41 } 42 #else 43 ierr = PetscViewerASCIIPrintf(viewer," %g\n",(double)xv[i]);CHKERRQ(ierr); 44 #endif 45 } 46 PetscFunctionReturn(0); 47 } 48 49 static PetscErrorCode VecView_Network_Seq(DM networkdm,Vec X,PetscViewer viewer) 50 { 51 PetscErrorCode ierr; 52 PetscInt e,v,Start,End,offset,nvar,id; 53 const PetscScalar *xv; 54 55 PetscFunctionBegin; 56 ierr = VecGetArrayRead(X,&xv);CHKERRQ(ierr); 57 58 /* iterate over edges */ 59 ierr = DMNetworkGetEdgeRange(networkdm,&Start,&End);CHKERRQ(ierr); 60 for (e=Start; e<End; e++) { 61 ierr = DMNetworkGetNumVariables(networkdm,e,&nvar);CHKERRQ(ierr); 62 if (!nvar) continue; 63 64 ierr = DMNetworkGetVariableOffset(networkdm,e,&offset);CHKERRQ(ierr); 65 ierr = DMNetworkGetGlobalEdgeIndex(networkdm,e,&id);CHKERRQ(ierr); 66 67 ierr = PetscViewerASCIIPrintf(viewer," Edge %D:\n",id);CHKERRQ(ierr); 68 ierr = VecArrayPrint_private(viewer,nvar,xv+offset);CHKERRQ(ierr); 69 } 70 71 /* iterate over vertices */ 72 ierr = DMNetworkGetVertexRange(networkdm,&Start,&End);CHKERRQ(ierr); 73 for (v=Start; v<End; v++) { 74 ierr = DMNetworkGetNumVariables(networkdm,v,&nvar);CHKERRQ(ierr); 75 if (!nvar) continue; 76 77 ierr = DMNetworkGetVariableOffset(networkdm,v,&offset);CHKERRQ(ierr); 78 ierr = DMNetworkGetGlobalVertexIndex(networkdm,v,&id);CHKERRQ(ierr); 79 80 ierr = PetscViewerASCIIPrintf(viewer," Vertex %D:\n",id);CHKERRQ(ierr); 81 ierr = VecArrayPrint_private(viewer,nvar,xv+offset);CHKERRQ(ierr); 82 } 83 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 84 ierr = VecRestoreArrayRead(X,&xv);CHKERRQ(ierr); 85 PetscFunctionReturn(0); 86 } 87 88 static PetscErrorCode VecView_Network_MPI(DM networkdm,Vec X,PetscViewer viewer) 89 { 90 PetscErrorCode ierr; 91 PetscInt i,e,v,eStart,eEnd,vStart,vEnd,offset,nvar,len_loc,len,k; 92 const PetscScalar *xv; 93 MPI_Comm comm; 94 PetscMPIInt size,rank,tag = ((PetscObject)viewer)->tag; 95 Vec localX; 96 PetscBool ghostvtex; 97 PetscScalar *values; 98 PetscInt j,ne,nv,id; 99 MPI_Status status; 100 101 PetscFunctionBegin; 102 ierr = PetscObjectGetComm((PetscObject)networkdm,&comm);CHKERRQ(ierr); 103 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 104 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 105 106 ierr = DMGetLocalVector(networkdm,&localX);CHKERRQ(ierr); 107 ierr = DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);CHKERRQ(ierr); 108 ierr = DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);CHKERRQ(ierr); 109 ierr = VecGetArrayRead(localX,&xv);CHKERRQ(ierr); 110 111 ierr = VecGetLocalSize(localX,&len_loc);CHKERRQ(ierr); 112 113 ierr = DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);CHKERRQ(ierr); 114 ierr = DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);CHKERRQ(ierr); 115 len_loc += 2*(1 + eEnd-eStart + vEnd-vStart); 116 117 /* values = [nedges, nvertices; id, nvar, xedge; ...; id, nvars, xvertex;...], to be sent to proc[0] */ 118 ierr = MPI_Allreduce(&len_loc,&len,1,MPIU_INT,MPI_MAX,comm);CHKERRQ(ierr); 119 ierr = PetscCalloc1(len,&values);CHKERRQ(ierr); 120 121 if (!rank) { 122 ierr = PetscViewerASCIIPrintf(viewer,"Process [%d]\n",rank);CHKERRQ(ierr); 123 } 124 125 /* iterate over edges */ 126 k = 2; 127 for (e=eStart; e<eEnd; e++) { 128 ierr = DMNetworkGetNumVariables(networkdm,e,&nvar);CHKERRQ(ierr); 129 if (!nvar) continue; 130 131 ierr = DMNetworkGetVariableOffset(networkdm,e,&offset);CHKERRQ(ierr); 132 ierr = DMNetworkGetGlobalEdgeIndex(networkdm,e,&id);CHKERRQ(ierr); 133 134 if (!rank) { /* print its own entries */ 135 ierr = PetscViewerASCIIPrintf(viewer," Edge %D:\n",id);CHKERRQ(ierr); 136 ierr = VecArrayPrint_private(viewer,nvar,xv+offset);CHKERRQ(ierr); 137 } else { 138 values[0] += 1; /* number of edges */ 139 values[k++] = id; 140 values[k++] = nvar; 141 for (i=offset; i< offset+nvar; i++) values[k++] = xv[i]; 142 } 143 } 144 145 /* iterate over vertices */ 146 for (v=vStart; v<vEnd; v++) { 147 ierr = DMNetworkIsGhostVertex(networkdm,v,&ghostvtex);CHKERRQ(ierr); 148 if (ghostvtex) continue; 149 ierr = DMNetworkGetNumVariables(networkdm,v,&nvar);CHKERRQ(ierr); 150 if (!nvar) continue; 151 152 ierr = DMNetworkGetVariableOffset(networkdm,v,&offset);CHKERRQ(ierr); 153 ierr = DMNetworkGetGlobalVertexIndex(networkdm,v,&id);CHKERRQ(ierr); 154 155 if (!rank) { 156 ierr = PetscViewerASCIIPrintf(viewer," Vertex %D:\n",id);CHKERRQ(ierr); 157 ierr = VecArrayPrint_private(viewer,nvar,xv+offset);CHKERRQ(ierr); 158 } else { 159 values[1] += 1; /* number of vertices */ 160 values[k++] = id; 161 values[k++] = nvar; 162 for (i=offset; i< offset+nvar; i++) values[k++] = xv[i]; 163 } 164 } 165 166 if (!rank) { 167 /* proc[0] receives and prints messages */ 168 for (j=1; j<size; j++) { 169 ierr = PetscViewerASCIIPrintf(viewer,"Process [%d]\n",j);CHKERRQ(ierr); 170 171 ierr = MPI_Recv(values,(PetscMPIInt)len,MPIU_SCALAR,j,tag,comm,&status);CHKERRQ(ierr); 172 173 ne = (PetscInt)PetscAbsScalar(values[0]); 174 nv = (PetscInt)PetscAbsScalar(values[1]); 175 176 /* print received edges */ 177 k = 2; 178 for (i=0; i<ne; i++) { 179 id = (PetscInt)PetscAbsScalar(values[k++]); 180 nvar = (PetscInt)PetscAbsScalar(values[k++]); 181 ierr = PetscViewerASCIIPrintf(viewer," Edge %D:\n",id);CHKERRQ(ierr); 182 ierr = VecArrayPrint_private(viewer,nvar,values+k);CHKERRQ(ierr); 183 k += nvar; 184 } 185 186 /* print received vertices */ 187 for (i=0; i<nv; i++) { 188 id = (PetscInt)PetscAbsScalar(values[k++]); 189 nvar = (PetscInt)PetscAbsScalar(values[k++]); 190 ierr = PetscViewerASCIIPrintf(viewer," Vertex %D:\n",id);CHKERRQ(ierr); 191 ierr = VecArrayPrint_private(viewer,nvar,values+k);CHKERRQ(ierr); 192 k += nvar; 193 } 194 } 195 } else { 196 /* sends values to proc[0] */ 197 ierr = MPI_Send((void*)values,k,MPIU_SCALAR,0,tag,comm);CHKERRQ(ierr); 198 } 199 200 ierr = PetscFree(values);CHKERRQ(ierr); 201 ierr = VecRestoreArrayRead(localX,&xv);CHKERRQ(ierr); 202 ierr = DMRestoreLocalVector(networkdm,&localX);CHKERRQ(ierr); 203 PetscFunctionReturn(0); 204 } 205 206 PETSC_EXTERN PetscErrorCode VecView_MPI(Vec,PetscViewer); 207 208 PetscErrorCode VecView_Network(Vec v,PetscViewer viewer) 209 { 210 DM dm; 211 PetscErrorCode ierr; 212 PetscBool isseq; 213 PetscBool iascii; 214 215 PetscFunctionBegin; 216 ierr = VecGetDM(v,&dm);CHKERRQ(ierr); 217 if (!dm) SETERRQ(PetscObjectComm((PetscObject)v),PETSC_ERR_ARG_WRONG,"Vector not generated from a DM"); 218 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 219 ierr = PetscObjectTypeCompare((PetscObject)v,VECSEQ,&isseq);CHKERRQ(ierr); 220 221 /* Use VecView_Network if the viewer is ASCII; use VecView_Seq/MPI for other viewer formats */ 222 if (iascii) { 223 if (isseq) { 224 ierr = VecView_Network_Seq(dm,v,viewer);CHKERRQ(ierr); 225 } else { 226 ierr = VecView_Network_MPI(dm,v,viewer);CHKERRQ(ierr); 227 } 228 } else { 229 if (isseq) { 230 ierr = VecView_Seq(v,viewer);CHKERRQ(ierr); 231 } else { 232 ierr = VecView_MPI(v,viewer);CHKERRQ(ierr); 233 } 234 } 235 PetscFunctionReturn(0); 236 } 237 238 static PetscErrorCode DMCreateGlobalVector_Network(DM dm,Vec *vec) 239 { 240 PetscErrorCode ierr; 241 DM_Network *network = (DM_Network*) dm->data; 242 243 PetscFunctionBegin; 244 ierr = DMCreateGlobalVector(network->plex,vec);CHKERRQ(ierr); 245 ierr = VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Network);CHKERRQ(ierr); 246 ierr = VecSetDM(*vec,dm);CHKERRQ(ierr); 247 PetscFunctionReturn(0); 248 } 249 250 static PetscErrorCode DMCreateLocalVector_Network(DM dm,Vec *vec) 251 { 252 PetscErrorCode ierr; 253 DM_Network *network = (DM_Network*) dm->data; 254 255 PetscFunctionBegin; 256 ierr = DMCreateLocalVector(network->plex,vec);CHKERRQ(ierr); 257 ierr = VecSetDM(*vec,dm);CHKERRQ(ierr); 258 PetscFunctionReturn(0); 259 } 260 261 PetscErrorCode DMInitialize_Network(DM dm) 262 { 263 PetscErrorCode ierr; 264 265 PetscFunctionBegin; 266 ierr = DMSetDimension(dm,1); CHKERRQ(ierr); 267 dm->ops->view = DMView_Network; 268 dm->ops->setfromoptions = DMSetFromOptions_Network; 269 dm->ops->clone = DMClone_Network; 270 dm->ops->setup = DMSetUp_Network; 271 dm->ops->createglobalvector = DMCreateGlobalVector_Network; 272 dm->ops->createlocalvector = DMCreateLocalVector_Network; 273 dm->ops->getlocaltoglobalmapping = NULL; 274 dm->ops->createfieldis = NULL; 275 dm->ops->createcoordinatedm = NULL; 276 dm->ops->getcoloring = 0; 277 dm->ops->creatematrix = DMCreateMatrix_Network; 278 dm->ops->createinterpolation = 0; 279 dm->ops->createinjection = 0; 280 dm->ops->refine = 0; 281 dm->ops->coarsen = 0; 282 dm->ops->refinehierarchy = 0; 283 dm->ops->coarsenhierarchy = 0; 284 dm->ops->globaltolocalbegin = DMGlobalToLocalBegin_Network; 285 dm->ops->globaltolocalend = DMGlobalToLocalEnd_Network; 286 dm->ops->localtoglobalbegin = DMLocalToGlobalBegin_Network; 287 dm->ops->localtoglobalend = DMLocalToGlobalEnd_Network; 288 dm->ops->destroy = DMDestroy_Network; 289 dm->ops->createsubdm = NULL; 290 dm->ops->locatepoints = NULL; 291 PetscFunctionReturn(0); 292 } 293 294 PetscErrorCode DMClone_Network(DM dm, DM *newdm) 295 { 296 DM_Network *network = (DM_Network *) dm->data; 297 PetscErrorCode ierr; 298 299 PetscFunctionBegin; 300 network->refct++; 301 (*newdm)->data = network; 302 ierr = PetscObjectChangeTypeName((PetscObject) *newdm, DMNETWORK);CHKERRQ(ierr); 303 ierr = DMInitialize_Network(*newdm);CHKERRQ(ierr); 304 PetscFunctionReturn(0); 305 } 306 307 /*MC 308 DMNETWORK = "network" - A DM object that encapsulates an unstructured network. The implementation is based on the DM object 309 DMPlex that manages unstructured grids. Distributed networks use a non-overlapping partitioning of 310 the edges. In the local representation, Vecs contain all unknowns in the interior and shared boundary. 311 This is specified by a PetscSection object. Ownership in the global representation is determined by 312 ownership of the underlying DMPlex points. This is specified by another PetscSection object. 313 314 Level: intermediate 315 316 .seealso: DMType, DMNetworkCreate(), DMCreate(), DMSetType() 317 M*/ 318 319 PETSC_EXTERN PetscErrorCode DMCreate_Network(DM dm) 320 { 321 DM_Network *network; 322 PetscErrorCode ierr; 323 324 PetscFunctionBegin; 325 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 326 ierr = PetscNewLog(dm,&network);CHKERRQ(ierr); 327 dm->data = network; 328 329 network->refct = 1; 330 network->NVertices = 0; 331 network->NEdges = 0; 332 network->nVertices = 0; 333 network->nEdges = 0; 334 network->nsubnet = 0; 335 336 337 ierr = DMInitialize_Network(dm);CHKERRQ(ierr); 338 PetscFunctionReturn(0); 339 } 340 341 /*@ 342 DMNetworkCreate - Creates a DMNetwork object, which encapsulates an unstructured network. 343 344 Collective 345 346 Input Parameter: 347 . comm - The communicator for the DMNetwork object 348 349 Output Parameter: 350 . network - The DMNetwork object 351 352 Level: beginner 353 354 @*/ 355 PetscErrorCode DMNetworkCreate(MPI_Comm comm, DM *network) 356 { 357 PetscErrorCode ierr; 358 359 PetscFunctionBegin; 360 PetscValidPointer(network,2); 361 ierr = DMCreate(comm, network);CHKERRQ(ierr); 362 ierr = DMSetType(*network, DMNETWORK);CHKERRQ(ierr); 363 PetscFunctionReturn(0); 364 } 365