Actual source code: ex10.c
petsc-3.4.2 2013-07-02

      2: /*
      3:    Include "petscsnes.h" so that we can use SNES solvers.  Note that this
      4:    file automatically includes:
      5:      petscsys.h       - base PETSc routines   petscvec.h - vectors
      6:      petscmat.h - matrices
      7:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
      8:      petscviewer.h - viewers               petscpc.h  - preconditioners
      9:      petscksp.h   - linear solvers
     10: */
     11: #include <petscsnes.h>
     12: #include <petscao.h>

     14: #if !defined(PETSC_USE_COMPLEX)

     16: static char help[] = "An Unstructured Grid Example.\n\
     17: This example demonstrates how to solve a nonlinear system in parallel\n\
     18: with SNES for an unstructured mesh. The mesh and partitioning information\n\
     19: is read in an application defined ordering,which is later transformed\n\
     20: into another convenient ordering (called the local ordering). The local\n\
     21: ordering, apart from being efficient on cpu cycles and memory, allows\n\
     22: the use of the SPMD model of parallel programming. After partitioning\n\
     23: is done, scatters are created between local (sequential)and global\n\
     24: (distributed) vectors. Finally, we set up the nonlinear solver context\n\
     25: in the usual way as a structured grid  (see\n\
     26: petsc/src/snes/examples/tutorials/ex5.c).\n\
     27: This example also illustrates the use of parallel matrix coloring.\n\
     28: The command line options include:\n\
     29:   -vert <Nv>, where Nv is the global number of nodes\n\
     30:   -elem <Ne>, where Ne is the global number of elements\n\
     31:   -nl_par <lambda>, where lambda is the multiplier for the non linear term (u*u) term\n\
     32:   -lin_par <alpha>, where alpha is the multiplier for the linear term (u)\n\
     33:   -fd_jacobian_coloring -mat_coloring_type lf\n";

     35: /*T
     36:    Concepts: SNES^unstructured grid
     37:    Concepts: AO^application to PETSc ordering or vice versa;
     38:    Concepts: VecScatter^using vector scatter operations;
     39:    Processors: n
     40: T*/

     42: /* ------------------------------------------------------------------------

     44:    PDE Solved : L(u) + lambda*u*u + alpha*u = 0 where L(u) is the Laplacian.

     46:    The Laplacian is approximated in the following way: each edge is given a weight
     47:    of one meaning that the diagonal term will have the weight equal to the degree
     48:    of a node. The off diagonal terms will get a weight of -1.

     50:    -----------------------------------------------------------------------*/

     53: #define MAX_ELEM      500  /* Maximum number of elements */
     54: #define MAX_VERT      100  /* Maximum number of vertices */
     55: #define MAX_VERT_ELEM   3  /* Vertices per element       */

     57: /*
     58:   Application-defined context for problem specific data
     59: */
     60: typedef struct {
     61:   PetscInt   Nvglobal,Nvlocal;              /* global and local number of vertices */
     62:   PetscInt   Neglobal,Nelocal;              /* global and local number of vertices */
     63:   PetscInt   AdjM[MAX_VERT][50];            /* adjacency list of a vertex */
     64:   PetscInt   itot[MAX_VERT];                /* total number of neighbors for a vertex */
     65:   PetscInt   icv[MAX_ELEM][MAX_VERT_ELEM];  /* vertices belonging to an element */
     66:   PetscInt   v2p[MAX_VERT];                 /* processor number for a vertex */
     67:   PetscInt   *locInd,*gloInd;               /* local and global orderings for a node */
     68:   Vec        localX,localF;                 /* local solution (u) and f(u) vectors */
     69:   PetscReal  non_lin_param;                 /* nonlinear parameter for the PDE */
     70:   PetscReal  lin_param;                     /* linear parameter for the PDE */
     71:   VecScatter scatter;                       /* scatter context for the local and
     72:                                                distributed vectors */
     73: } AppCtx;

     75: /*
     76:   User-defined routines
     77: */
     78: PetscErrorCode FormJacobian(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
     79: PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
     80: PetscErrorCode FormInitialGuess(AppCtx*,Vec);

     84: int main(int argc,char **argv)
     85: {
     86:   SNES                   snes;                 /* SNES context */
     87:   SNESType               type = SNESNEWTONLS;  /* default nonlinear solution method */
     88:   Vec                    x,r;                  /* solution, residual vectors */
     89:   Mat                    Jac;                  /* Jacobian matrix */
     90:   AppCtx                 user;                 /* user-defined application context */
     91:   AO                     ao;                   /* Application Ordering object */
     92:   IS                     isglobal,islocal;     /* global and local index sets */
     93:   PetscMPIInt            rank,size;            /* rank of a process, number of processors */
     94:   PetscInt               rstart;               /* starting index of PETSc ordering for a processor */
     95:   PetscInt               nfails;               /* number of unsuccessful Newton steps */
     96:   PetscInt               bs = 1;               /* block size for multicomponent systems */
     97:   PetscInt               nvertices;            /* number of local plus ghost nodes of a processor */
     98:   PetscInt               *pordering;           /* PETSc ordering */
     99:   PetscInt               *vertices;            /* list of all vertices (incl. ghost ones) on a processor */
    100:   PetscInt               *verticesmask;
    101:   PetscInt               *tmp;
    102:   PetscInt               i,j,jstart,inode,nb,nbrs,Nvneighborstotal = 0;
    103:   PetscErrorCode         ierr;
    104:   PetscInt               its,N;
    105:   PetscScalar            *xx;
    106:   char                   str[256],form[256],part_name[256];
    107:   FILE                   *fptr,*fptr1;
    108:   ISLocalToGlobalMapping isl2g;
    109:   int                    dtmp;
    110: #if defined(UNUSED_VARIABLES)
    111:   PetscDraw              draw;                 /* drawing context */
    112:   PetscScalar            *ff,*gg;
    113:   PetscReal              tiny = 1.0e-10,zero = 0.0,one = 1.0,big = 1.0e+10;
    114:   PetscInt               *tmp1,*tmp2;
    115: #endif
    116:   MatFDColoring          matfdcoloring = 0;
    117:   PetscBool              fd_jacobian_coloring = PETSC_FALSE;

    119:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    120:      Initialize program
    121:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    123:   PetscInitialize(&argc,&argv,"options.inf",help);
    124:   MPI_Comm_rank(MPI_COMM_WORLD,&rank);
    125:   MPI_Comm_size(MPI_COMM_WORLD,&size);

    127:   /* The current input file options.inf is for 2 proc run only */
    128:   if (size != 2) SETERRQ(PETSC_COMM_SELF,1,"This Example currently runs on 2 procs only.");

    130:   /*
    131:      Initialize problem parameters
    132:   */
    133:   user.Nvglobal = 16;      /*Global # of vertices  */
    134:   user.Neglobal = 18;      /*Global # of elements  */

    136:   PetscOptionsGetInt(NULL,"-vert",&user.Nvglobal,NULL);
    137:   PetscOptionsGetInt(NULL,"-elem",&user.Neglobal,NULL);

    139:   user.non_lin_param = 0.06;

    141:   PetscOptionsGetReal(NULL,"-nl_par",&user.non_lin_param,NULL);

    143:   user.lin_param = -1.0;

    145:   PetscOptionsGetReal(NULL,"-lin_par",&user.lin_param,NULL);

    147:   user.Nvlocal = 0;
    148:   user.Nelocal = 0;

    150:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    151:       Read the mesh and partitioning information
    152:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    154:   /*
    155:      Read the mesh and partitioning information from 'adj.in'.
    156:      The file format is as follows.
    157:      For each line the first entry is the processor rank where the
    158:      current node belongs. The second entry is the number of
    159:      neighbors of a node. The rest of the line is the adjacency
    160:      list of a node. Currently this file is set up to work on two
    161:      processors.

    163:      This is not a very good example of reading input. In the future,
    164:      we will put an example that shows the style that should be
    165:      used in a real application, where partitioning will be done
    166:      dynamically by calling partitioning routines (at present, we have
    167:      a  ready interface to ParMeTiS).
    168:    */
    169:   fptr = fopen("adj.in","r");
    170:   if (!fptr) SETERRQ(PETSC_COMM_SELF,0,"Could not open adj.in");

    172:   /*
    173:      Each processor writes to the file output.<rank> where rank is the
    174:      processor's rank.
    175:   */
    176:   sprintf(part_name,"output.%d",rank);
    177:   fptr1 = fopen(part_name,"w");
    178:   if (!fptr1) SETERRQ(PETSC_COMM_SELF,0,"Could no open output file");
    179:   PetscMalloc(user.Nvglobal*sizeof(PetscInt),&user.gloInd);
    180:   PetscFPrintf(PETSC_COMM_SELF,fptr1,"Rank is %D\n",rank);
    181:   for (inode = 0; inode < user.Nvglobal; inode++) {
    182:     if (!fgets(str,256,fptr)) SETERRQ(PETSC_COMM_SELF,1,"fgets read failed");
    183:     sscanf(str,"%d",&dtmp);user.v2p[inode] = dtmp;
    184:     if (user.v2p[inode] == rank) {
    185:       PetscFPrintf(PETSC_COMM_SELF,fptr1,"Node %D belongs to processor %D\n",inode,user.v2p[inode]);

    187:       user.gloInd[user.Nvlocal] = inode;
    188:       sscanf(str,"%*d %d",&dtmp);
    189:       nbrs = dtmp;
    190:       PetscFPrintf(PETSC_COMM_SELF,fptr1,"Number of neighbors for the vertex %D is %D\n",inode,nbrs);

    192:       user.itot[user.Nvlocal] = nbrs;
    193:       Nvneighborstotal       += nbrs;
    194:       for (i = 0; i < user.itot[user.Nvlocal]; i++) {
    195:         form[0]='\0';
    196:         for (j=0; j < i+2; j++) {
    197:           PetscStrcat(form,"%*d ");
    198:         }
    199:         PetscStrcat(form,"%d");

    201:         sscanf(str,form,&dtmp);
    202:         user.AdjM[user.Nvlocal][i] = dtmp;

    204:         PetscFPrintf(PETSC_COMM_SELF,fptr1,"%D ",user.AdjM[user.Nvlocal][i]);
    205:       }
    206:       PetscFPrintf(PETSC_COMM_SELF,fptr1,"\n");
    207:       user.Nvlocal++;
    208:     }
    209:   }
    210:   PetscFPrintf(PETSC_COMM_SELF,fptr1,"Total # of Local Vertices is %D \n",user.Nvlocal);

    212:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    213:      Create different orderings
    214:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    216:   /*
    217:     Create the local ordering list for vertices. First a list using the PETSc global
    218:     ordering is created. Then we use the AO object to get the PETSc-to-application and
    219:     application-to-PETSc mappings. Each vertex also gets a local index (stored in the
    220:     locInd array).
    221:   */
    222:   MPI_Scan(&user.Nvlocal,&rstart,1,MPIU_INT,MPI_SUM,PETSC_COMM_WORLD);
    223:   rstart -= user.Nvlocal;
    224:   PetscMalloc(user.Nvlocal*sizeof(PetscInt),&pordering);

    226:   for (i=0; i < user.Nvlocal; i++) pordering[i] = rstart + i;

    228:   /*
    229:     Create the AO object
    230:   */
    231:   AOCreateBasic(MPI_COMM_WORLD,user.Nvlocal,user.gloInd,pordering,&ao);
    232:   PetscFree(pordering);

    234:   /*
    235:     Keep the global indices for later use
    236:   */
    237:   PetscMalloc(user.Nvlocal*sizeof(PetscInt),&user.locInd);
    238:   PetscMalloc(Nvneighborstotal*sizeof(PetscInt),&tmp);

    240:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    241:     Demonstrate the use of AO functionality
    242:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    244:   PetscFPrintf(PETSC_COMM_SELF,fptr1,"Before AOApplicationToPetsc, local indices are : \n");
    245:   for (i=0; i < user.Nvlocal; i++) {
    246:     PetscFPrintf(PETSC_COMM_SELF,fptr1," %D ",user.gloInd[i]);

    248:     user.locInd[i] = user.gloInd[i];
    249:   }
    250:   PetscFPrintf(PETSC_COMM_SELF,fptr1,"\n");
    251:   jstart = 0;
    252:   for (i=0; i < user.Nvlocal; i++) {
    253:     PetscFPrintf(PETSC_COMM_SELF,fptr1,"Neghbors of local vertex %D are : ",user.gloInd[i]);
    254:     for (j=0; j < user.itot[i]; j++) {
    255:       PetscFPrintf(PETSC_COMM_SELF,fptr1,"%D ",user.AdjM[i][j]);

    257:       tmp[j + jstart] = user.AdjM[i][j];
    258:     }
    259:     jstart += user.itot[i];
    260:     PetscFPrintf(PETSC_COMM_SELF,fptr1,"\n");
    261:   }

    263:   /*
    264:     Now map the vlocal and neighbor lists to the PETSc ordering
    265:   */
    266:   AOApplicationToPetsc(ao,user.Nvlocal,user.locInd);
    267:   AOApplicationToPetsc(ao,Nvneighborstotal,tmp);
    268:   AODestroy(&ao);

    270:   PetscFPrintf(PETSC_COMM_SELF,fptr1,"After AOApplicationToPetsc, local indices are : \n");
    271:   for (i=0; i < user.Nvlocal; i++) {
    272:     PetscFPrintf(PETSC_COMM_SELF,fptr1," %D ",user.locInd[i]);
    273:   }
    274:   PetscFPrintf(PETSC_COMM_SELF,fptr1,"\n");

    276:   jstart = 0;
    277:   for (i=0; i < user.Nvlocal; i++) {
    278:     PetscFPrintf(PETSC_COMM_SELF,fptr1,"Neghbors of local vertex %D are : ",user.locInd[i]);
    279:     for (j=0; j < user.itot[i]; j++) {
    280:       user.AdjM[i][j] = tmp[j+jstart];

    282:       PetscFPrintf(PETSC_COMM_SELF,fptr1,"%D ",user.AdjM[i][j]);
    283:     }
    284:     jstart += user.itot[i];
    285:     PetscFPrintf(PETSC_COMM_SELF,fptr1,"\n");
    286:   }

    288:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    289:      Extract the ghost vertex information for each processor
    290:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    291:   /*
    292:    Next, we need to generate a list of vertices required for this processor
    293:    and a local numbering scheme for all vertices required on this processor.
    294:       vertices - integer array of all vertices needed on this processor in PETSc
    295:                  global numbering; this list consists of first the "locally owned"
    296:                  vertices followed by the ghost vertices.
    297:       verticesmask - integer array that for each global vertex lists its local
    298:                      vertex number (in vertices) + 1. If the global vertex is not
    299:                      represented on this processor, then the corresponding
    300:                      entry in verticesmask is zero

    302:       Note: vertices and verticesmask are both Nvglobal in length; this may
    303:     sound terribly non-scalable, but in fact is not so bad for a reasonable
    304:     number of processors. Importantly, it allows us to use NO SEARCHING
    305:     in setting up the data structures.
    306:   */
    307:   PetscMalloc(user.Nvglobal*sizeof(PetscInt),&vertices);
    308:   PetscMalloc(user.Nvglobal*sizeof(PetscInt),&verticesmask);
    309:   PetscMemzero(verticesmask,user.Nvglobal*sizeof(PetscInt));
    310:   nvertices = 0;

    312:   /*
    313:     First load "owned vertices" into list
    314:   */
    315:   for (i=0; i < user.Nvlocal; i++) {
    316:     vertices[nvertices++]        = user.locInd[i];
    317:     verticesmask[user.locInd[i]] = nvertices;
    318:   }

    320:   /*
    321:     Now load ghost vertices into list
    322:   */
    323:   for (i=0; i < user.Nvlocal; i++) {
    324:     for (j=0; j < user.itot[i]; j++) {
    325:       nb = user.AdjM[i][j];
    326:       if (!verticesmask[nb]) {
    327:         vertices[nvertices++] = nb;
    328:         verticesmask[nb]      = nvertices;
    329:       }
    330:     }
    331:   }

    333:   PetscFPrintf(PETSC_COMM_SELF,fptr1,"\n");
    334:   PetscFPrintf(PETSC_COMM_SELF,fptr1,"The array vertices is :\n");
    335:   for (i=0; i < nvertices; i++) {
    336:     PetscFPrintf(PETSC_COMM_SELF,fptr1,"%D ",vertices[i]);
    337:   }
    338:   PetscFPrintf(PETSC_COMM_SELF,fptr1,"\n");

    340:   /*
    341:      Map the vertices listed in the neighbors to the local numbering from
    342:     the global ordering that they contained initially.
    343:   */
    344:   PetscFPrintf(PETSC_COMM_SELF,fptr1,"\n");
    345:   PetscFPrintf(PETSC_COMM_SELF,fptr1,"After mapping neighbors in the local contiguous ordering\n");
    346:   for (i=0; i<user.Nvlocal; i++) {
    347:     PetscFPrintf(PETSC_COMM_SELF,fptr1,"Neghbors of local vertex %D are :\n",i);
    348:     for (j = 0; j < user.itot[i]; j++) {
    349:       nb              = user.AdjM[i][j];
    350:       user.AdjM[i][j] = verticesmask[nb] - 1;

    352:       PetscFPrintf(PETSC_COMM_SELF,fptr1,"%D ",user.AdjM[i][j]);
    353:     }
    354:     PetscFPrintf(PETSC_COMM_SELF,fptr1,"\n");
    355:   }

    357:   N = user.Nvglobal;

    359:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    360:      Create vector and matrix data structures
    361:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    363:   /*
    364:     Create vector data structures
    365:   */
    366:   VecCreate(MPI_COMM_WORLD,&x);
    367:   VecSetSizes(x,user.Nvlocal,N);
    368:   VecSetFromOptions(x);
    369:   VecDuplicate(x,&r);
    370:   VecCreateSeq(MPI_COMM_SELF,bs*nvertices,&user.localX);
    371:   VecDuplicate(user.localX,&user.localF);

    373:   /*
    374:     Create the scatter between the global representation and the
    375:     local representation
    376:   */
    377:   ISCreateStride(MPI_COMM_SELF,bs*nvertices,0,1,&islocal);
    378:   ISCreateBlock(MPI_COMM_SELF,bs,nvertices,vertices,PETSC_COPY_VALUES,&isglobal);
    379:   VecScatterCreate(x,isglobal,user.localX,islocal,&user.scatter);
    380:   ISDestroy(&isglobal);
    381:   ISDestroy(&islocal);

    383:   /*
    384:      Create matrix data structure; Just to keep the example simple, we have not done any
    385:      preallocation of memory for the matrix. In real application code with big matrices,
    386:      preallocation should always be done to expedite the matrix creation.
    387:   */
    388:   MatCreate(MPI_COMM_WORLD,&Jac);
    389:   MatSetSizes(Jac,PETSC_DECIDE,PETSC_DECIDE,N,N);
    390:   MatSetFromOptions(Jac);
    391:   MatSetUp(Jac);

    393:   /*
    394:     The following routine allows us to set the matrix values in local ordering
    395:   */
    396:   ISLocalToGlobalMappingCreate(MPI_COMM_SELF,bs*nvertices,vertices,PETSC_COPY_VALUES,&isl2g);
    397:   MatSetLocalToGlobalMapping(Jac,isl2g,isl2g);

    399:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    400:      Create nonlinear solver context
    401:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    403:   SNESCreate(MPI_COMM_WORLD,&snes);
    404:   SNESSetType(snes,type);

    406:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    407:     Set routines for function and Jacobian evaluation
    408:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    409:   SNESSetFunction(snes,r,FormFunction,(void*)&user);

    411:   PetscOptionsGetBool(NULL,"-fd_jacobian_coloring",&fd_jacobian_coloring,0);
    412:   if (!fd_jacobian_coloring) {
    413:     SNESSetJacobian(snes,Jac,Jac,FormJacobian,(void*)&user);
    414:   } else {  /* Use matfdcoloring */
    415:     ISColoring   iscoloring;
    416:     MatStructure flag;
    417:     /* Get the data structure of Jac */
    418:     FormJacobian(snes,x,&Jac,&Jac,&flag,&user);
    419:     /* Create coloring context */
    420:     MatGetColoring(Jac,MATCOLORINGSL,&iscoloring);
    421:     MatFDColoringCreate(Jac,iscoloring,&matfdcoloring);
    422:     MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))FormFunction,&user);
    423:     MatFDColoringSetFromOptions(matfdcoloring);
    424:     /* MatFDColoringView(matfdcoloring,PETSC_VIEWER_STDOUT_WORLD); */
    425:     SNESSetJacobian(snes,Jac,Jac,SNESComputeJacobianDefaultColor,matfdcoloring);
    426:     ISColoringDestroy(&iscoloring);
    427:   }

    429:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    430:      Customize nonlinear solver; set runtime options
    431:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    433:   SNESSetFromOptions(snes);

    435:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    436:      Evaluate initial guess; then solve nonlinear system
    437:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    439:   /*
    440:      Note: The user should initialize the vector, x, with the initial guess
    441:      for the nonlinear solver prior to calling SNESSolve().  In particular,
    442:      to employ an initial guess of zero, the user should explicitly set
    443:      this vector to zero by calling VecSet().
    444:   */
    445:   FormInitialGuess(&user,x);

    447:   /*
    448:     Print the initial guess
    449:   */
    450:   VecGetArray(x,&xx);
    451:   for (inode = 0; inode < user.Nvlocal; inode++) {
    452:     PetscFPrintf(PETSC_COMM_SELF,fptr1,"Initial Solution at node %D is %f \n",inode,xx[inode]);
    453:   }
    454:   VecRestoreArray(x,&xx);

    456:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    457:      Now solve the nonlinear system
    458:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    460:   SNESSolve(snes,NULL,x);
    461:   SNESGetIterationNumber(snes,&its);
    462:   SNESGetNonlinearStepFailures(snes,&nfails);

    464:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    465:     Print the output : solution vector and other information
    466:      Each processor writes to the file output.<rank> where rank is the
    467:      processor's rank.
    468:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    470:   VecGetArray(x,&xx);
    471:   for (inode = 0; inode < user.Nvlocal; inode++) {
    472:     PetscFPrintf(PETSC_COMM_SELF,fptr1,"Solution at node %D is %f \n",inode,xx[inode]);
    473:   }
    474:   VecRestoreArray(x,&xx);
    475:   fclose(fptr1);
    476:   PetscPrintf(MPI_COMM_WORLD,"number of SNES iterations = %D, ",its);
    477:   PetscPrintf(MPI_COMM_WORLD,"number of unsuccessful steps = %D\n",nfails);

    479:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    480:      Free work space.  All PETSc objects should be destroyed when they
    481:      are no longer needed.
    482:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    483:   PetscFree(user.gloInd);
    484:   PetscFree(user.locInd);
    485:   PetscFree(vertices);
    486:   PetscFree(verticesmask);
    487:   PetscFree(tmp);
    488:   VecScatterDestroy(&user.scatter);
    489:   ISLocalToGlobalMappingDestroy(&isl2g);
    490:   VecDestroy(&x);
    491:   VecDestroy(&r);
    492:   VecDestroy(&user.localX);
    493:   VecDestroy(&user.localF);
    494:   MatDestroy(&Jac);  SNESDestroy(&snes);
    495:   /*PetscDrawDestroy(draw);*/
    496:   if (fd_jacobian_coloring) {
    497:     MatFDColoringDestroy(&matfdcoloring);
    498:   }
    499:   PetscFinalize();

    501:   return 0;
    502: }
    505: /* --------------------  Form initial approximation ----------------- */

    507: /*
    508:    FormInitialGuess - Forms initial approximation.

    510:    Input Parameters:
    511:    user - user-defined application context
    512:    X - vector

    514:    Output Parameter:
    515:    X - vector
    516:  */
    517: PetscErrorCode FormInitialGuess(AppCtx *user,Vec X)
    518: {
    519:   PetscInt    i,Nvlocal,ierr;
    520:   PetscInt    *gloInd;
    521:   PetscScalar *x;
    522: #if defined(UNUSED_VARIABLES)
    523:   PetscReal temp1,temp,hx,hy,hxdhy,hydhx,sc;
    524:   PetscInt  Neglobal,Nvglobal,j,row;
    525:   PetscReal alpha,lambda;

    527:   Nvglobal = user->Nvglobal;
    528:   Neglobal = user->Neglobal;
    529:   lambda   = user->non_lin_param;
    530:   alpha    = user->lin_param;
    531: #endif

    533:   Nvlocal = user->Nvlocal;
    534:   gloInd  = user->gloInd;

    536:   /*
    537:      Get a pointer to vector data.
    538:        - For default PETSc vectors, VecGetArray() returns a pointer to
    539:          the data array.  Otherwise, the routine is implementation dependent.
    540:        - You MUST call VecRestoreArray() when you no longer need access to
    541:          the array.
    542:   */
    543:   VecGetArray(X,&x);

    545:   /*
    546:      Compute initial guess over the locally owned part of the grid
    547:   */
    548:   for (i=0; i < Nvlocal; i++) x[i] = (PetscReal)gloInd[i];

    550:   /*
    551:      Restore vector
    552:   */
    553:   VecRestoreArray(X,&x);
    554:   return 0;
    555: }
    558: /* --------------------  Evaluate Function F(x) --------------------- */
    559: /*
    560:    FormFunction - Evaluates nonlinear function, F(x).

    562:    Input Parameters:
    563: .  snes - the SNES context
    564: .  X - input vector
    565: .  ptr - optional user-defined context, as set by SNESSetFunction()

    567:    Output Parameter:
    568: .  F - function vector
    569:  */
    570: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ptr)
    571: {
    572:   AppCtx         *user = (AppCtx*)ptr;
    574:   PetscInt       i,j,Nvlocal;
    575:   PetscReal      alpha,lambda;
    576:   PetscScalar    *x,*f;
    577:   VecScatter     scatter;
    578:   Vec            localX = user->localX;
    579: #if defined(UNUSED_VARIABLES)
    580:   PetscScalar ut,ub,ul,ur,u,*g,sc,uyy,uxx;
    581:   PetscReal   hx,hy,hxdhy,hydhx;
    582:   PetscReal   two = 2.0,one = 1.0;
    583:   PetscInt    Nvglobal,Neglobal,row;
    584:   PetscInt    *gloInd;

    586:   Nvglobal = user->Nvglobal;
    587:   Neglobal = user->Neglobal;
    588:   gloInd   = user->gloInd;
    589: #endif

    591:   Nvlocal = user->Nvlocal;
    592:   lambda  = user->non_lin_param;
    593:   alpha   = user->lin_param;
    594:   scatter = user->scatter;

    596:   /*
    597:      PDE : L(u) + lambda*u*u +alpha*u = 0 where L(u) is the approximate Laplacian as
    598:      described in the beginning of this code

    600:      First scatter the distributed vector X into local vector localX (that includes
    601:      values for ghost nodes. If we wish,we can put some other work between
    602:      VecScatterBegin() and VecScatterEnd() to overlap the communication with
    603:      computation.
    604:  */
    605:   VecScatterBegin(scatter,X,localX,INSERT_VALUES,SCATTER_FORWARD);
    606:   VecScatterEnd(scatter,X,localX,INSERT_VALUES,SCATTER_FORWARD);

    608:   /*
    609:      Get pointers to vector data
    610:   */
    611:   VecGetArray(localX,&x);
    612:   VecGetArray(F,&f);

    614:   /*
    615:     Now compute the f(x). As mentioned earlier, the computed Laplacian is just an
    616:     approximate one chosen for illustrative purpose only. Another point to notice
    617:     is that this is a local (completly parallel) calculation. In practical application
    618:     codes, function calculation time is a dominat portion of the overall execution time.
    619:   */
    620:   for (i=0; i < Nvlocal; i++) {
    621:     f[i] = (user->itot[i] - alpha)*x[i] - lambda*x[i]*x[i];
    622:     for (j = 0; j < user->itot[i]; j++) f[i] -= x[user->AdjM[i][j]];
    623:   }

    625:   /*
    626:      Restore vectors
    627:   */
    628:   VecRestoreArray(localX,&x);
    629:   VecRestoreArray(F,&f);
    630:   /*VecView(F,PETSC_VIEWER_STDOUT_WORLD);*/

    632:   return 0;
    633: }

    637: /* --------------------  Evaluate Jacobian F'(x) -------------------- */
    638: /*
    639:    FormJacobian - Evaluates Jacobian matrix.

    641:    Input Parameters:
    642: .  snes - the SNES context
    643: .  X - input vector
    644: .  ptr - optional user-defined context, as set by SNESSetJacobian()

    646:    Output Parameters:
    647: .  A - Jacobian matrix
    648: .  B - optionally different preconditioning matrix
    649: .  flag - flag indicating matrix structure

    651: */
    652: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flag,void *ptr)
    653: {
    654:   AppCtx      *user = (AppCtx*)ptr;
    655:   Mat         jac   = *B;
    656:   PetscInt    i,j,Nvlocal,col[50],ierr;
    657:   PetscScalar alpha,lambda,value[50];
    658:   Vec         localX = user->localX;
    659:   VecScatter  scatter;
    660:   PetscScalar *x;
    661: #if defined(UNUSED_VARIABLES)
    662:   PetscScalar two = 2.0,one = 1.0;
    663:   PetscInt    row,Nvglobal,Neglobal;
    664:   PetscInt    *gloInd;

    666:   Nvglobal = user->Nvglobal;
    667:   Neglobal = user->Neglobal;
    668:   gloInd   = user->gloInd;
    669: #endif

    671:   /*printf("Entering into FormJacobian \n");*/
    672:   Nvlocal = user->Nvlocal;
    673:   lambda  = user->non_lin_param;
    674:   alpha   =  user->lin_param;
    675:   scatter = user->scatter;

    677:   /*
    678:      PDE : L(u) + lambda*u*u +alpha*u = 0 where L(u) is the approximate Laplacian as
    679:      described in the beginning of this code

    681:      First scatter the distributed vector X into local vector localX (that includes
    682:      values for ghost nodes. If we wish, we can put some other work between
    683:      VecScatterBegin() and VecScatterEnd() to overlap the communication with
    684:      computation.
    685:   */
    686:   VecScatterBegin(scatter,X,localX,INSERT_VALUES,SCATTER_FORWARD);
    687:   VecScatterEnd(scatter,X,localX,INSERT_VALUES,SCATTER_FORWARD);

    689:   /*
    690:      Get pointer to vector data
    691:   */
    692:   VecGetArray(localX,&x);

    694:   for (i=0; i < Nvlocal; i++) {
    695:     col[0]   = i;
    696:     value[0] = user->itot[i] - 2.0*lambda*x[i] - alpha;
    697:     for (j = 0; j < user->itot[i]; j++) {
    698:       col[j+1]   = user->AdjM[i][j];
    699:       value[j+1] = -1.0;
    700:     }

    702:     /*
    703:       Set the matrix values in the local ordering. Note that in order to use this
    704:       feature we must call the routine MatSetLocalToGlobalMapping() after the
    705:       matrix has been created.
    706:     */
    707:     MatSetValuesLocal(jac,1,&i,1+user->itot[i],col,value,INSERT_VALUES);
    708:   }

    710:   /*
    711:      Assemble matrix, using the 2-step process:
    712:        MatAssemblyBegin(), MatAssemblyEnd().
    713:      Between these two calls, the pointer to vector data has been restored to
    714:      demonstrate the use of overlapping communicationn with computation.
    715:   */
    716:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
    717:   VecRestoreArray(localX,&x);
    718:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);

    720:   /*
    721:      Set flag to indicate that the Jacobian matrix retains an identical
    722:      nonzero structure throughout all nonlinear iterations (although the
    723:      values of the entries change). Thus, we can save some work in setting
    724:      up the preconditioner (e.g., no need to redo symbolic factorization for
    725:      ILU/ICC preconditioners).
    726:       - If the nonzero structure of the matrix is different during
    727:         successive linear solves, then the flag DIFFERENT_NONZERO_PATTERN
    728:         must be used instead.  If you are unsure whether the matrix
    729:         structure has changed or not, use the flag DIFFERENT_NONZERO_PATTERN.
    730:       - Caution:  If you specify SAME_NONZERO_PATTERN, PETSc
    731:         believes your assertion and does not check the structure
    732:         of the matrix.  If you erroneously claim that the structure
    733:         is the same when it actually is not, the new preconditioner
    734:         will not function correctly.  Thus, use this optimization
    735:         feature with caution!
    736:   */
    737:   *flag = SAME_NONZERO_PATTERN;

    739:   /*
    740:      Tell the matrix we will never add a new nonzero location to the
    741:      matrix. If we do, it will generate an error.
    742:   */
    743:   MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
    744:   /* MatView(jac,PETSC_VIEWER_STDOUT_SELF); */
    745:   return 0;
    746: }
    747: #else

    749: int main(int argc,char **args)
    750: {
    751:   fprintf(stdout,"This example does not work for complex numbers.\n");
    752:   return 0;
    753: }
    754: #endif