[petsc-users] Modeling a Network with DMPlex?
Abhyankar, Shrirang G.
abhyshr at mcs.anl.gov
Sun Feb 16 19:50:45 CST 2014
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" <florian.meier at koalo.de> 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;
>}
More information about the petsc-users
mailing list