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