Thank you very much for your review. On 02/17/2014 02:50 AM, Abhyankar, Shrirang G. wrote: > Hi Florian, > I was able to run your example with a small change (used PetscMalloc > instead of PetscMalloc1).
That is strange. When I use PetscMalloc, it says "Memory corruption!". > However, the solution on more than 1 processor > is different presumably because of some work "TODO" for the ghost vertices. I do not really understand these ghost vertices. I thought it would be ok to just skip those when I iterate over the vertices, but not when I iterate over the edges, but apparently that is not the case. > 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. I do not know if I really change the non-zero pattern. What I mean is the following: Some entries in the Jacobian MIGHT be zero not because of the structure of the problem is changing, but just because the partial derivative for a special x is (near) zero by chance. > 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. That looks great! Do I have to set the non-constant elements in MatStoreValues, too, so that the structure does not change? > 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. It has? I would say that is not the case since all partial derivatives might be non-zero when there is a relation (look into the PDF I send you). But maybe I just do not know what you mean by 2 X 2 structure or what is mean to NOT have a 2 X 2 structure. > 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. That is a good idea. I will try to save it into the component array. Can I assume that FormFunction always gets called before FormJacobian with the same x? > 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 again, Florian > 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; >> } >
