Hi Florian, I was able to run your example with a small change (used PetscMalloc instead of PetscMalloc1). However, the solution on more than 1 processor is different presumably because of some work "TODO" for the ghost vertices. Answers to some of your questions in the code: I) You can change the non-zero pattern of the Jacobian matrix. But note that if the network topology then I think the circuit would need to be re-built.
II) If some of the Jacobian matrix entries are constant then you can use MatStoreValues http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatStoreVal ues.html#MatStoreValues to set the constant terms once before the Newton loop and retrieve these via MatRetrieveValues http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatRetrieve Values.html#MatRetrieveValues during the Jacobian evaluation. III) The Jacobian matrix for your application has a 2 X 2 structure. So you can set the entires in each block through one MatSetValues call. IV) If you want to reuse some term that calculated in the function evaluation (Pfree) during the Jacobian evaluation then you can either use your own user context or save the term in the component array. Thanks once again for contributing this example. Once you have fixed the ghost vertex part and get the same solution with different number of processors, I'll add it to the PETSc repository. Let us know if you have any more questions. Thanks, Shri On 2/14/14 11:25 AM, "Florian Meier" <[email protected]> wrote: >Ok, that was just a stupid typo. >I have now build the first attempt of an example. >I am sure that can be programmed in a much better way, >so I would be very glad about any comment. What is your >preferred review procedure? > >Greetings, >Florian > >static char help[] = "This example demonstrates the DMCircuit interface >for a radio network flow problem.\n\ > The available solver options are in the >alohaoptions file.\n\ > The file links.txt contains the network >description.\n\ > Run this program: mpiexec -n <n> ./PF\n > \ > mpiexec -n <n> ./PF -pfdata <filename>\n"; > >/** >This computes a greatly simplified model of a radio network. >Each node in the network generates packets conforming to a poisson >distribution (given rate R_gen). >These packets are routed along fixed paths to a sink. When two packets >are send at the same time, >they are dropped. No channel sensing or retransmission takes place. > >Each active link (i.e. each radio link that will eventually transmit >data) is modeled as a vertex. >Between those vertices there exist two kinds of relations: >A relation of type REL_INTERFERENCE(affected, source), implies that >whenever the source link tries >to send at the same time than the affected link, the packet will collide >and is not received properly. >A relation of type REL_INFLOW(affected, source) implies that each packet >that is send via the source >link and does not collide, adds up to the amount of packets send via the >affected link. > >The overall rate of accessing a radio link i is >R_access,i = R_gen,i + SUM_OVER_INFLOWS( R_success,j ) > >The probability that no other link interfering with i accesses the >channel at a given point in time is >P_free,i = PRODUCT_OVER_INTERFERERS( 1 - R_access,j ) > >Finally, the overall rate of successful transmissions over link i is >R_success,i = R_access,i * P_free,i > >The input file is structured as follows: >The first line contains L, the number of links, I, the number of >interferences and F, the number of inflow relations. >Then L lines follow with the packet generation rate for each link. >The subsequent I lines describe the affected and the source link of the >respective interference relation. >The final F lines represent the inflow relations. >*/ > >#include <petsc.h> >#include <petscdmcircuit.h> > >#include <fstream> >#include <iostream> > >using namespace std; > >enum { > REL_INTERFERENCE, > REL_INFLOW >}; > >enum { > VAR_ACCESS, > VAR_SUCCESS, > VAR_NVARS >}; > >struct _p_LINKDATA{ > PetscScalar packet_generation; >}; > >typedef struct _p_LINKDATA *LINKDATA; > >struct _p_RELATIONDATA{ > PetscInt source; > PetscInt affected; > PetscInt type; >}; > >typedef struct _p_RELATIONDATA *RELATIONDATA; > >typedef struct { > PetscInt nlinks,nrelations; /* # of nodes,relations */ > LINKDATA links; > RELATIONDATA relations; >} PFDATA; > > >PetscMPIInt rank; > >#undef __FUNCT__ >#define __FUNCT__ "FormFunction" >PetscErrorCode FormFunction(SNES snes,Vec X, Vec F,void *appctx) >{ > PetscErrorCode ierr; > DM circuitdm; > Vec localX,localF; > PetscInt e; > PetscInt v,vStart,vEnd,vaffected,vsource; > const PetscScalar *xarr; > PetscScalar *farr; > PetscInt offset,offsetsource,offsetlink,offsetrel; > DMCircuitComponentGenericDataType *arr; > > PetscFunctionBegin; > ierr = SNESGetDM(snes,&circuitdm);CHKERRQ(ierr); > ierr = DMGetLocalVector(circuitdm,&localX);CHKERRQ(ierr); > ierr = DMGetLocalVector(circuitdm,&localF);CHKERRQ(ierr); > ierr = VecSet(F,0.0);CHKERRQ(ierr); > > ierr = >DMGlobalToLocalBegin(circuitdm,X,INSERT_VALUES,localX);CHKERRQ(ierr); > ierr = >DMGlobalToLocalEnd(circuitdm,X,INSERT_VALUES,localX);CHKERRQ(ierr); > > ierr = >DMGlobalToLocalBegin(circuitdm,F,INSERT_VALUES,localF);CHKERRQ(ierr); > ierr = >DMGlobalToLocalEnd(circuitdm,F,INSERT_VALUES,localF);CHKERRQ(ierr); > > ierr = VecGetArrayRead(localX,&xarr);CHKERRQ(ierr); > ierr = VecGetArray(localF,&farr);CHKERRQ(ierr); > > ierr = DMCircuitGetVertexRange(circuitdm,&vStart,&vEnd);CHKERRQ(ierr); > ierr = DMCircuitGetComponentDataArray(circuitdm,&arr);CHKERRQ(ierr); > > for (v=vStart; v < vEnd; v++) { > PetscBool ghostvtex; > ierr = DMCircuitIsGhostVertex(circuitdm,v,&ghostvtex);CHKERRQ(ierr); > if (ghostvtex) { > continue; // TODO is that ok? > } > > PetscInt keyv; > ierr = >DMCircuitGetComponentTypeOffset(circuitdm,v,0,&keyv,&offsetlink);CHKERRQ(i >err); > > LINKDATA link = (LINKDATA)(arr+offsetlink); > > PetscScalar inflow = link->packet_generation; > PetscScalar Pfree = 1; > > PetscInt nconnedges; > const PetscInt *connedges; > ierr = >DMCircuitGetSupportingEdges(circuitdm,v,&nconnedges,&connedges);CHKERRQ(ie >rr); > for (PetscInt i = 0; i < nconnedges; i++) { > e = connedges[i]; > > const PetscInt *cone; > ierr = DMCircuitGetConnectedNodes(circuitdm,e,&cone);CHKERRQ(ierr); > vaffected = cone[0]; > vsource = cone[1]; > > if (vaffected == v) { > PetscInt keye; > ierr = >DMCircuitGetComponentTypeOffset(circuitdm,e,0,&keye,&offsetrel);CHKERRQ(ie >rr); > RELATIONDATA relation = (RELATIONDATA)(arr+offsetrel); > > ierr = >DMCircuitGetVariableOffset(circuitdm,vsource,&offsetsource);CHKERRQ(ierr); > > switch (relation->type) { > case REL_INTERFERENCE: > Pfree *= 1 - xarr[offsetsource+VAR_ACCESS]; > break; > case REL_INFLOW: > inflow += xarr[offsetsource+VAR_SUCCESS]; > break; > } > } > } > > ierr = DMCircuitGetVariableOffset(circuitdm,v,&offset);CHKERRQ(ierr); > farr[offset+VAR_ACCESS] = inflow - xarr[offset+VAR_ACCESS]; > farr[offset+VAR_SUCCESS] = xarr[offset+VAR_ACCESS]*Pfree - >xarr[offset+VAR_SUCCESS]; > } > > ierr = VecRestoreArrayRead(localX,&xarr);CHKERRQ(ierr); > ierr = VecRestoreArray(localF,&farr);CHKERRQ(ierr); > ierr = DMRestoreLocalVector(circuitdm,&localX);CHKERRQ(ierr); > > ierr = >DMLocalToGlobalBegin(circuitdm,localF,ADD_VALUES,F);CHKERRQ(ierr); > ierr = DMLocalToGlobalEnd(circuitdm,localF,ADD_VALUES,F);CHKERRQ(ierr); > ierr = DMRestoreLocalVector(circuitdm,&localF);CHKERRQ(ierr); > > PetscFunctionReturn(0); >} > > >#undef __FUNCT__ >#define __FUNCT__ "FormJacobian" >PetscErrorCode FormJacobian(SNES snes,Vec X, Mat *J,Mat >*Jpre,MatStructure *flg,void *appctx) >{ > PetscErrorCode ierr; > DM circuitdm; > Vec localX; > PetscInt e; > PetscInt v,vStart,vEnd,vaffected,vsource; > const PetscScalar *xarr; > PetscInt offsetrel; > DMCircuitComponentGenericDataType *arr; > PetscInt row[1],col[1]; > PetscScalar values[1]; > > PetscFunctionBegin; > *flg = SAME_NONZERO_PATTERN; // TODO ok for this problem? > ierr = MatZeroEntries(*J);CHKERRQ(ierr); > > ierr = SNESGetDM(snes,&circuitdm);CHKERRQ(ierr); > ierr = DMGetLocalVector(circuitdm,&localX);CHKERRQ(ierr); > > ierr = >DMGlobalToLocalBegin(circuitdm,X,INSERT_VALUES,localX);CHKERRQ(ierr); > ierr = >DMGlobalToLocalEnd(circuitdm,X,INSERT_VALUES,localX);CHKERRQ(ierr); > > ierr = VecGetArrayRead(localX,&xarr);CHKERRQ(ierr); > > ierr = DMCircuitGetVertexRange(circuitdm,&vStart,&vEnd);CHKERRQ(ierr); > ierr = DMCircuitGetComponentDataArray(circuitdm,&arr);CHKERRQ(ierr); > > for (v=vStart; v < vEnd; v++) { > PetscInt offset,goffset,offsetsource,goffsetsource; > PetscBool ghostvtex; > > ierr = DMCircuitIsGhostVertex(circuitdm,v,&ghostvtex);CHKERRQ(ierr); > if (ghostvtex) { > continue; // TODO is that ok? > } > > ierr = >DMCircuitGetVariableGlobalOffset(circuitdm,v,&goffset);CHKERRQ(ierr); > > // TODO some derivatives are constant, can this be handled in >SetInitialValues? > > // TODO can I combine these two MatSetValues? > row[0] = goffset+VAR_ACCESS; > col[0] = goffset+VAR_ACCESS; > values[0] = -1.0; > ierr = MatSetValues(*J,1,row,1,col,values,ADD_VALUES);CHKERRQ(ierr); > row[0] = goffset+VAR_SUCCESS; > col[0] = goffset+VAR_SUCCESS; > values[0] = -1.0; > ierr = MatSetValues(*J,1,row,1,col,values,ADD_VALUES);CHKERRQ(ierr); > > // Calculate Pfree > // TODO Pfree was already calculated in FromFunction, can I access >that? > PetscScalar Pfree = 1; > > PetscInt nconnedges; > const PetscInt *connedges; > ierr = >DMCircuitGetSupportingEdges(circuitdm,v,&nconnedges,&connedges);CHKERRQ(ie >rr); > for (PetscInt i = 0; i < nconnedges; i++) { > e = connedges[i]; > > const PetscInt *cone; > ierr = DMCircuitGetConnectedNodes(circuitdm,e,&cone);CHKERRQ(ierr); > vaffected = cone[0]; > vsource = cone[1]; > > if (vaffected == v) { > PetscInt keye; > ierr = >DMCircuitGetComponentTypeOffset(circuitdm,e,0,&keye,&offsetrel);CHKERRQ(ie >rr); > RELATIONDATA relation = (RELATIONDATA)(arr+offsetrel); > > ierr = >DMCircuitGetVariableOffset(circuitdm,vsource,&offsetsource);CHKERRQ(ierr); > > if (relation->type == REL_INTERFERENCE) { > Pfree *= 1 - xarr[offsetsource+VAR_ACCESS]; > } > } > } > > row[0] = goffset+VAR_SUCCESS; > col[0] = goffset+VAR_ACCESS; > values[0] = Pfree; > ierr = MatSetValues(*J,1,row,1,col,values,ADD_VALUES);CHKERRQ(ierr); > > // Set the derivatives > ierr = DMCircuitGetVariableOffset(circuitdm,v,&offset);CHKERRQ(ierr); > PetscScalar inflow = xarr[offset+VAR_ACCESS]; > > for (PetscInt i = 0; i < nconnedges; i++) { > e = connedges[i]; > > const PetscInt *cone; > ierr = DMCircuitGetConnectedNodes(circuitdm,e,&cone);CHKERRQ(ierr); > vaffected = cone[0]; > vsource = cone[1]; > > if (vaffected == v) { > PetscInt keye; > ierr = >DMCircuitGetComponentTypeOffset(circuitdm,e,0,&keye,&offsetrel);CHKERRQ(ie >rr); > RELATIONDATA relation = (RELATIONDATA)(arr+offsetrel); > > ierr = >DMCircuitGetVariableOffset(circuitdm,vsource,&offsetsource);CHKERRQ(ierr); > ierr = >DMCircuitGetVariableGlobalOffset(circuitdm,vsource,&goffsetsource);CHKERRQ >(ierr); > > switch (relation->type) { > case REL_INTERFERENCE: > // TODO Ok to set the entries one by one or is there a better >way? > row[0] = goffset+VAR_SUCCESS; > col[0] = goffsetsource+VAR_ACCESS; > values[0] = -inflow*(Pfree/(1-xarr[offsetsource+VAR_ACCESS])); > ierr = >MatSetValues(*J,1,row,1,col,values,ADD_VALUES);CHKERRQ(ierr); > break; > case REL_INFLOW: > row[0] = goffset+VAR_ACCESS; > col[0] = goffsetsource+VAR_SUCCESS; > values[0] = 1.0; > ierr = >MatSetValues(*J,1,row,1,col,values,ADD_VALUES);CHKERRQ(ierr); > break; > } > } > } > } > > ierr = VecRestoreArrayRead(localX,&xarr);CHKERRQ(ierr); > ierr = DMRestoreLocalVector(circuitdm,&localX);CHKERRQ(ierr); > > ierr = MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); > ierr = MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); > > PetscFunctionReturn(0); >} > > > >#undef __FUNCT__ >#define __FUNCT__ "SetInitialValues" >PetscErrorCode SetInitialValues(DM circuitdm,Vec X,void *appctx) >{ > PetscErrorCode ierr; > PetscInt v, vStart, vEnd; > Vec localX; > PetscScalar *xarr; > DMCircuitComponentGenericDataType *arr; > > PetscFunctionBegin; > ierr = DMCircuitGetVertexRange(circuitdm,&vStart,&vEnd);CHKERRQ(ierr); > > ierr = DMGetLocalVector(circuitdm,&localX);CHKERRQ(ierr); > > ierr = VecSet(X,0.0);CHKERRQ(ierr); > ierr = >DMGlobalToLocalBegin(circuitdm,X,INSERT_VALUES,localX);CHKERRQ(ierr); > ierr = >DMGlobalToLocalEnd(circuitdm,X,INSERT_VALUES,localX);CHKERRQ(ierr); > > ierr = DMCircuitGetComponentDataArray(circuitdm,&arr);CHKERRQ(ierr); > ierr = VecGetArray(localX,&xarr);CHKERRQ(ierr); > for (v = vStart; v < vEnd; v++) { > PetscBool ghostvtex; > ierr = DMCircuitIsGhostVertex(circuitdm,v,&ghostvtex);CHKERRQ(ierr); > if (ghostvtex) { > continue; // TODO is that ok? > } > > PetscInt offsetlink, offset; > PetscInt keyv; > ierr = DMCircuitGetVariableOffset(circuitdm,v,&offset);CHKERRQ(ierr); > ierr = >DMCircuitGetComponentTypeOffset(circuitdm,v,0,&keyv,&offsetlink);CHKERRQ(i >err); > > LINKDATA link = (LINKDATA)(arr+offsetlink); > > xarr[offset+VAR_ACCESS] = link->packet_generation; > xarr[offset+VAR_SUCCESS] = link->packet_generation; > } > > ierr = VecRestoreArray(localX,&xarr);CHKERRQ(ierr); > ierr = >DMLocalToGlobalBegin(circuitdm,localX,ADD_VALUES,X);CHKERRQ(ierr); > ierr = DMLocalToGlobalEnd(circuitdm,localX,ADD_VALUES,X);CHKERRQ(ierr); > ierr = DMRestoreLocalVector(circuitdm,&localX);CHKERRQ(ierr); > PetscFunctionReturn(0); >} > > >#undef __FUNCT__ >#define __FUNCT__ "main" >int main(int argc,char ** argv) >{ > string inputFile = "links.txt"; > PetscErrorCode ierr; > PFDATA pfdata; > PetscInt numEdges=0,numVertices=0; > int *edges = NULL; > PetscInt i; > DM circuitdm; > PetscInt componentkey[2]; > PetscLogStage stage1,stage2; > PetscInt size; > > PetscInitialize(&argc,&argv,"alohaoptions",help); > > ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr); > > /* Create an empty circuit object */ > ierr = DMCircuitCreate(PETSC_COMM_WORLD,&circuitdm);CHKERRQ(ierr); > > /* Register the components in the circuit */ > ierr = DMCircuitRegisterComponent(circuitdm,"linkstruct",sizeof(struct >_p_LINKDATA),&componentkey[0]);CHKERRQ(ierr); > ierr = >DMCircuitRegisterComponent(circuitdm,"relationstruct",sizeof(struct >_p_RELATIONDATA),&componentkey[1]);CHKERRQ(ierr); > > ierr = PetscLogStageRegister("Read Data",&stage1);CHKERRQ(ierr); > PetscLogStagePush(stage1); > > /* READ THE DATA */ > if (!rank) { > /* Only rank 0 reads the data */ > ifstream linksF(inputFile.c_str()); > > int ninterferences, ninflows; > linksF >> pfdata.nlinks >> ninterferences >> ninflows; > > numVertices = pfdata.nlinks; > numEdges = pfdata.nrelations = ninterferences + ninflows; > > ierr = PetscMalloc1(pfdata.nlinks,&pfdata.links);CHKERRQ(ierr); > ierr = >PetscMalloc1(pfdata.nrelations,&pfdata.relations);CHKERRQ(ierr); > > for(int i = 0; i < numVertices; i++) { > linksF >> pfdata.links[i].packet_generation; > } > > ierr = PetscMalloc1(2*numEdges,&edges);CHKERRQ(ierr); > > for(int i = 0; i < numEdges; i++) { > linksF >> pfdata.relations[i].affected >> >pfdata.relations[i].source; > > pfdata.relations[i].type = (i < ninterferences) ? REL_INTERFERENCE >: REL_INFLOW; > > edges[2*i] = pfdata.relations[i].affected; > edges[2*i+1] = pfdata.relations[i].source; > } > > linksF.close(); > } > > PetscLogStagePop(); > MPI_Barrier(PETSC_COMM_WORLD); > > ierr = PetscLogStageRegister("Create circuit",&stage2);CHKERRQ(ierr); > PetscLogStagePush(stage2); > /* Set number of nodes/edges */ > ierr = >DMCircuitSetSizes(circuitdm,numVertices,numEdges,PETSC_DETERMINE,PETSC_DET >ERMINE);CHKERRQ(ierr); > /* Add edge connectivity */ > ierr = DMCircuitSetEdgeList(circuitdm,edges);CHKERRQ(ierr); > /* Set up the circuit layout */ > ierr = DMCircuitLayoutSetUp(circuitdm);CHKERRQ(ierr); > > if (!rank) { > ierr = PetscFree(edges);CHKERRQ(ierr); > } > > /* Add circuit components */ > PetscInt eStart, eEnd, vStart, vEnd; > > ierr = DMCircuitGetVertexRange(circuitdm,&vStart,&vEnd);CHKERRQ(ierr); > for (i = vStart; i < vEnd; i++) { > ierr = >DMCircuitAddComponent(circuitdm,i,componentkey[0],&pfdata.links[i-vStart]) >;CHKERRQ(ierr); > > /* Add number of variables */ > ierr = DMCircuitAddNumVariables(circuitdm,i,VAR_NVARS);CHKERRQ(ierr); > } > > ierr = DMCircuitGetEdgeRange(circuitdm,&eStart,&eEnd);CHKERRQ(ierr); > for (i = eStart; i < eEnd; i++) { > ierr = >DMCircuitAddComponent(circuitdm,i,componentkey[1],&pfdata.relations[i-eSta >rt]);CHKERRQ(ierr); > } > > /* Set up DM for use */ > ierr = DMSetUp(circuitdm);CHKERRQ(ierr); > > if (!rank) { > ierr = PetscFree(pfdata.links);CHKERRQ(ierr); > ierr = PetscFree(pfdata.relations);CHKERRQ(ierr); > } > > > ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr); > if (size > 1) { > DM distcircuitdm; > /* Circuit partitioning and distribution of data */ > ierr = DMCircuitDistribute(circuitdm,&distcircuitdm);CHKERRQ(ierr); > ierr = DMDestroy(&circuitdm);CHKERRQ(ierr); > circuitdm = distcircuitdm; > } > > PetscLogStagePop(); > > Vec X,F; > ierr = DMCreateGlobalVector(circuitdm,&X);CHKERRQ(ierr); > ierr = VecDuplicate(X,&F);CHKERRQ(ierr); > > Mat J; > ierr = DMCreateMatrix(circuitdm,&J);CHKERRQ(ierr); > ierr = >MatSetOption(J,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);CHKERRQ(ierr); > > ierr = SetInitialValues(circuitdm,X,NULL);CHKERRQ(ierr); > > SNES snes; > /* HOOK UP SOLVER */ > ierr = SNESCreate(PETSC_COMM_WORLD,&snes);CHKERRQ(ierr); > ierr = SNESSetDM(snes,circuitdm);CHKERRQ(ierr); > ierr = SNESSetFunction(snes,F,FormFunction,NULL);CHKERRQ(ierr); > ierr = SNESSetJacobian(snes,J,J,FormJacobian,NULL);CHKERRQ(ierr); > ierr = SNESSetFromOptions(snes);CHKERRQ(ierr); > > ierr = SNESSolve(snes,NULL,X);CHKERRQ(ierr); > ierr = VecView(X,NULL);CHKERRQ(ierr); > > ierr = VecDestroy(&X);CHKERRQ(ierr); > ierr = VecDestroy(&F);CHKERRQ(ierr); > ierr = MatDestroy(&J);CHKERRQ(ierr); > > ierr = SNESDestroy(&snes);CHKERRQ(ierr); > ierr = DMDestroy(&circuitdm);CHKERRQ(ierr); > > PetscFinalize(); > return 0; >}
