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