Actual source code: ex10.c

  1: /*$Id: ex10.c,v 1.29 2001/09/11 16:34:10 bsmith Exp $*/

  3: /* 
  4:   Program usage:  mpirun -np <procs> usg [-help] [all PETSc options] 
  5: */

  7: #if !defined(PETSC_USE_COMPLEX)

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

 26: /*T
 27:    Concepts: SNES^unstructured grid
 28:    Concepts: AO^application to PETSc ordering or vice versa;
 29:    Concepts: VecScatter^using vector scatter operations;
 30:    Processors: n
 31: T*/

 33: /* ------------------------------------------------------------------------

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

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

 41:    -----------------------------------------------------------------------*/

 43: /*
 44:    Include petscao.h so that we can use AO (Application Ordering) object's services.
 45:    Include "petscsnes.h" so that we can use SNES solvers.  Note that this
 46:    file automatically includes:
 47:      petsc.h       - base PETSc routines   petscvec.h - vectors
 48:      petscsys.h    - system routines       petscmat.h - matrices
 49:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
 50:      petscviewer.h - viewers               petscpc.h  - preconditioners
 51:      petscsles.h   - linear solvers
 52: */
 53: #include "petscao.h"
 54: #include "petscsnes.h"


 57: #define MAX_ELEM      500  /* Maximum number of elements */
 58: #define MAX_VERT      100  /* Maximum number of vertices */
 59: #define MAX_VERT_ELEM   3  /* Vertices per element       */

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

 79: /*
 80:   User-defined routines
 81: */
 82: int  FormJacobian(SNES,Vec,Mat*,Mat*,MatStructure*,void*),
 83:      FormFunction(SNES,Vec,Vec,void*),
 84:      FormInitialGuess(AppCtx*,Vec);

 88: int main(int argc,char **argv)
 89: {
 90:   SNES     snes;                 /* SNES context */
 91:   SNESType type = SNESLS;        /* default nonlinear solution method */
 92:   Vec      x,r;                  /* solution, residual vectors */
 93:   Mat      Jac;                  /* Jacobian matrix */
 94:   AppCtx   user;                 /* user-defined application context */
 95:   AO       ao;                   /* Application Ordering object */
 96:   IS       isglobal,islocal;     /* global and local index sets */
 97:   int           rank,size;            /* rank of a process, number of processors */
 98:   int      rstart;               /* starting index of PETSc ordering for a processor */
 99:   int      nfails;               /* number of unsuccessful Newton steps */
100:   int      bs = 1;               /* block size for multicomponent systems */
101:   int      nvertices;            /* number of local plus ghost nodes of a processor */
102:   int            *pordering;           /* PETSc ordering */
103:   int      *vertices;            /* list of all vertices (incl. ghost ones) 
104:                                     on a processor */
105:   int      *verticesmask,*svertices;
106:   int      *tmp;
107:   int      i,j,jstart,inode,nb,nbrs,Nvneighborstotal = 0;
108:   int      ierr,its,N;
109:   PetscScalar   *xx;
110:   char     str[256],form[256],part_name[256];
111:   FILE     *fptr,*fptr1;
112:   ISLocalToGlobalMapping isl2g;
113: #if defined (UNUSED_VARIABLES)
114:   PetscDraw    draw;                 /* drawing context */
115:   PetscScalar  *ff,*gg;
116:   PetscReal  tiny = 1.0e-10,zero = 0.0,one = 1.0,big = 1.0e+10;
117:   int     *tmp1,*tmp2;
118: #endif
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(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  */
135:   PetscOptionsGetInt(PETSC_NULL,"-vert",&user.Nvglobal,PETSC_NULL);
136:   PetscOptionsGetInt(PETSC_NULL,"-elem",&user.Neglobal,PETSC_NULL);
137:   user.non_lin_param = 0.06;
138:   PetscOptionsGetReal(PETSC_NULL,"-nl_par",&user.non_lin_param,PETSC_NULL);
139:   user.lin_param = -1.0;
140:   PetscOptionsGetReal(PETSC_NULL,"-lin_par",&user.lin_param,PETSC_NULL);
141:   user.Nvlocal = 0;
142:   user.Nelocal = 0;

144:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
145:       Read the mesh and partitioning information  
146:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
147: 
148:   /*
149:      Read the mesh and partitioning information from 'adj.in'.
150:      The file format is as follows.
151:      For each line the first entry is the processor rank where the 
152:      current node belongs. The second entry is the number of 
153:      neighbors of a node. The rest of the line is the adjacency
154:      list of a node. Currently this file is set up to work on two 
155:      processors.

157:      This is not a very good example of reading input. In the future, 
158:      we will put an example that shows the style that should be 
159:      used in a real application, where partitioning will be done 
160:      dynamically by calling partitioning routines (at present, we have 
161:      a  ready interface to ParMeTiS). 
162:    */
163:   fptr = fopen("adj.in","r");
164:   if (!fptr) {
165:       SETERRQ(0,"Could not open adj.in")
166:   }
167: 
168:   /*
169:      Each processor writes to the file output.<rank> where rank is the
170:      processor's rank.
171:   */
172:   sprintf(part_name,"output.%d",rank);
173:   fptr1 = fopen(part_name,"w");
174:   if (!fptr1) {
175:       SETERRQ(0,"Could no open output file");
176:   }
177:   PetscMalloc(user.Nvglobal*sizeof(int),&user.gloInd);
178:   fprintf(fptr1,"Rank is %d\n",rank);
179:   for (inode = 0; inode < user.Nvglobal; inode++) {
180:     fgets(str,256,fptr);
181:     sscanf(str,"%d",&user.v2p[inode]);
182:     if (user.v2p[inode] == rank) {
183:        fprintf(fptr1,"Node %d belongs to processor %d\n",inode,user.v2p[inode]);
184:        user.gloInd[user.Nvlocal] = inode;
185:        sscanf(str,"%*d %d",&nbrs);
186:        fprintf(fptr1,"Number of neighbors for the vertex %d is %d\n",inode,nbrs);
187:        user.itot[user.Nvlocal] = nbrs;
188:        Nvneighborstotal += nbrs;
189:        for (i = 0; i < user.itot[user.Nvlocal]; i++){
190:          form[0]='\0';
191:          for (j=0; j < i+2; j++){
192:            PetscStrcat(form,"%*d ");
193:          }
194:            PetscStrcat(form,"%d");
195:            sscanf(str,form,&user.AdjM[user.Nvlocal][i]);
196:            fprintf(fptr1,"%d ",user.AdjM[user.Nvlocal][i]);
197:         }
198:         fprintf(fptr1,"\n");
199:         user.Nvlocal++;
200:      }
201:    }
202:   fprintf(fptr1,"Total # of Local Vertices is %d \n",user.Nvlocal);

204:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
205:      Create different orderings
206:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

208:   /*
209:     Create the local ordering list for vertices. First a list using the PETSc global 
210:     ordering is created. Then we use the AO object to get the PETSc-to-application and 
211:     application-to-PETSc mappings. Each vertex also gets a local index (stored in the 
212:     locInd array). 
213:   */
214:   MPI_Scan(&user.Nvlocal,&rstart,1,MPI_INT,MPI_SUM,MPI_COMM_WORLD);
215:   rstart -= user.Nvlocal;
216:   PetscMalloc(user.Nvlocal*sizeof(int),&pordering);

218:   for (i=0; i < user.Nvlocal; i++) {
219:     pordering[i] = rstart + i;
220:   }

222:   /* 
223:     Create the AO object 
224:   */
225:   AOCreateBasic(MPI_COMM_WORLD,user.Nvlocal,user.gloInd,pordering,&ao);
226:   PetscFree(pordering);
227: 
228:   /* 
229:     Keep the global indices for later use 
230:   */
231:   PetscMalloc(user.Nvlocal*sizeof(int),&user.locInd);
232:   PetscMalloc(Nvneighborstotal*sizeof(int),&tmp);
233: 
234:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
235:     Demonstrate the use of AO functionality 
236:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

238:   fprintf(fptr1,"Before AOApplicationToPetsc, local indices are : \n");
239:   for (i=0; i < user.Nvlocal; i++) {
240:    fprintf(fptr1," %d ",user.gloInd[i]);
241:    user.locInd[i] = user.gloInd[i];
242:   }
243:   fprintf(fptr1,"\n");
244:   jstart = 0;
245:   for (i=0; i < user.Nvlocal; i++) {
246:    fprintf(fptr1,"Neghbors of local vertex %d are : ",user.gloInd[i]);
247:    for (j=0; j < user.itot[i]; j++) {
248:     fprintf(fptr1,"%d ",user.AdjM[i][j]);
249:     tmp[j + jstart] = user.AdjM[i][j];
250:    }
251:    jstart += user.itot[i];
252:    fprintf(fptr1,"\n");
253:   }

255:   /* 
256:     Now map the vlocal and neighbor lists to the PETSc ordering 
257:   */
258:   AOApplicationToPetsc(ao,user.Nvlocal,user.locInd);
259:   AOApplicationToPetsc(ao,Nvneighborstotal,tmp);
260:   AODestroy(ao);

262:   fprintf(fptr1,"After AOApplicationToPetsc, local indices are : \n");
263:   for (i=0; i < user.Nvlocal; i++) {
264:    fprintf(fptr1," %d ",user.locInd[i]);
265:   }
266:   fprintf(fptr1,"\n");

268:   jstart = 0;
269:   for (i=0; i < user.Nvlocal; i++) {
270:    fprintf(fptr1,"Neghbors of local vertex %d are : ",user.locInd[i]);
271:    for (j=0; j < user.itot[i]; j++) {
272:     user.AdjM[i][j] = tmp[j+jstart];
273:     fprintf(fptr1,"%d ",user.AdjM[i][j]);
274:    }
275:    jstart += user.itot[i];
276:    fprintf(fptr1,"\n");
277:   }

279:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
280:      Extract the ghost vertex information for each processor
281:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
282:   /*
283:    Next, we need to generate a list of vertices required for this processor
284:    and a local numbering scheme for all vertices required on this processor.
285:       vertices - integer array of all vertices needed on this processor in PETSc 
286:                  global numbering; this list consists of first the "locally owned" 
287:                  vertices followed by the ghost vertices.
288:       verticesmask - integer array that for each global vertex lists its local 
289:                      vertex number (in vertices) + 1. If the global vertex is not
290:                      represented on this processor, then the corresponding
291:                      entry in verticesmask is zero
292:  
293:       Note: vertices and verticesmask are both Nvglobal in length; this may
294:     sound terribly non-scalable, but in fact is not so bad for a reasonable
295:     number of processors. Importantly, it allows us to use NO SEARCHING
296:     in setting up the data structures.
297:   */
298:   PetscMalloc(user.Nvglobal*sizeof(int),&vertices);
299:   PetscMalloc(user.Nvglobal*sizeof(int),&verticesmask);
300:   PetscMemzero(verticesmask,user.Nvglobal*sizeof(int));
301:   nvertices = 0;
302: 
303:   /* 
304:     First load "owned vertices" into list 
305:   */
306:   for (i=0; i < user.Nvlocal; i++) {
307:     vertices[nvertices++]   = user.locInd[i];
308:     verticesmask[user.locInd[i]] = nvertices;
309:   }
310: 
311:   /* 
312:     Now load ghost vertices into list 
313:   */
314:   for (i=0; i < user.Nvlocal; i++) {
315:     for (j=0; j < user.itot[i]; j++) {
316:       nb = user.AdjM[i][j];
317:       if (!verticesmask[nb]) {
318:         vertices[nvertices++] = nb;
319:         verticesmask[nb]      = nvertices;
320:       }
321:     }
322:   }

324:   fprintf(fptr1,"\n");
325:   fprintf(fptr1,"The array vertices is :\n");
326:   for (i=0; i < nvertices; i++) {
327:    fprintf(fptr1,"%d ",vertices[i]);
328:    }
329:   fprintf(fptr1,"\n");
330: 
331:   /*
332:      Map the vertices listed in the neighbors to the local numbering from
333:     the global ordering that they contained initially.
334:   */
335:   fprintf(fptr1,"\n");
336:   fprintf(fptr1,"After mapping neighbors in the local contiguous ordering\n");
337:   for (i=0; i<user.Nvlocal; i++) {
338:     fprintf(fptr1,"Neghbors of local vertex %d are :\n",i);
339:     for (j = 0; j < user.itot[i]; j++) {
340:       nb = user.AdjM[i][j];
341:       user.AdjM[i][j] = verticesmask[nb] - 1;
342:       fprintf(fptr1,"%d ",user.AdjM[i][j]);
343:     }
344:    fprintf(fptr1,"\n");
345:   }

347:   N = user.Nvglobal;
348: 
349:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
350:      Create vector and matrix data structures
351:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

353:   /* 
354:     Create vector data structures 
355:   */
356:   VecCreate(MPI_COMM_WORLD,&x);
357:   VecSetSizes(x,user.Nvlocal,N);
358:   VecSetFromOptions(x);
359:   VecDuplicate(x,&r);
360:   VecCreateSeq(MPI_COMM_SELF,bs*nvertices,&user.localX);
361:   VecDuplicate(user.localX,&user.localF);

363:   /*
364:     Create the scatter between the global representation and the 
365:     local representation
366:   */
367:   ISCreateStride(MPI_COMM_SELF,bs*nvertices,0,1,&islocal);
368:   PetscMalloc(nvertices*sizeof(int),&svertices);
369:   for (i=0; i<nvertices; i++) svertices[i] = bs*vertices[i];
370:   ISCreateBlock(MPI_COMM_SELF,bs,nvertices,svertices,&isglobal);
371:   PetscFree(svertices);
372:   VecScatterCreate(x,isglobal,user.localX,islocal,&user.scatter);

374:   /* 
375:      Create matrix data structure; Just to keep the example simple, we have not done any 
376:      preallocation of memory for the matrix. In real application code with big matrices,
377:      preallocation should always be done to expedite the matrix creation. 
378:   */
379:   MatCreate(MPI_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,N,N,&Jac);
380:   MatSetFromOptions(Jac);

382:   /* 
383:     The following routine allows us to set the matrix values in local ordering 
384:   */
385:   ISLocalToGlobalMappingCreate(MPI_COMM_SELF,bs*nvertices,vertices,&isl2g);
386:   MatSetLocalToGlobalMapping(Jac,isl2g);

388:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
389:      Create nonlinear solver context
390:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

392:   SNESCreate(MPI_COMM_WORLD,&snes);
393:   SNESSetType(snes,type);

395:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
396:     Set routines for function and Jacobian evaluation
397:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

399:    FormInitialGuess(&user,x);
400:    SNESSetFunction(snes,r,FormFunction,(void *)&user);
401:    SNESSetJacobian(snes,Jac,Jac,FormJacobian,(void *)&user);

403:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
404:      Customize nonlinear solver; set runtime options
405:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

407:   SNESSetFromOptions(snes);

409:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
410:      Evaluate initial guess; then solve nonlinear system
411:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

413:   /*
414:      Note: The user should initialize the vector, x, with the initial guess
415:      for the nonlinear solver prior to calling SNESSolve().  In particular,
416:      to employ an initial guess of zero, the user should explicitly set
417:      this vector to zero by calling VecSet().
418:   */
419:    FormInitialGuess(&user,x);

421:    /* 
422:      Print the initial guess 
423:    */
424:    VecGetArray(x,&xx);
425:    for (inode = 0; inode < user.Nvlocal; inode++)
426:     fprintf(fptr1,"Initial Solution at node %d is %f \n",inode,xx[inode]);
427:    VecRestoreArray(x,&xx);

429:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
430:      Now solve the nonlinear system
431:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

433:   SNESSolve(snes,x,&its);
434:   SNESGetNumberUnsuccessfulSteps(snes,&nfails);
435: 
436:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
437:     Print the output : solution vector and other information
438:      Each processor writes to the file output.<rank> where rank is the
439:      processor's rank.
440:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

442:   VecGetArray(x,&xx);
443:   for (inode = 0; inode < user.Nvlocal; inode++)
444:    fprintf(fptr1,"Solution at node %d is %f \n",inode,xx[inode]);
445:   VecRestoreArray(x,&xx);
446:   fclose(fptr1);
447:   PetscPrintf(MPI_COMM_WORLD,"number of Newton iterations = %d, ",its);
448:   PetscPrintf(MPI_COMM_WORLD,"number of unsuccessful steps = %d\n",nfails);

450:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
451:      Free work space.  All PETSc objects should be destroyed when they
452:      are no longer needed.
453:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

455:   VecDestroy(x);
456:   VecDestroy(r);
457:   VecDestroy(user.localX);
458:   VecDestroy(user.localF);
459:   MatDestroy(Jac);  SNESDestroy(snes);
460:   /*PetscDrawDestroy(draw);*/
461:   PetscFinalize();

463:   return 0;
464: }
467: /* --------------------  Form initial approximation ----------------- */

469: /* 
470:    FormInitialGuess - Forms initial approximation.

472:    Input Parameters:
473:    user - user-defined application context
474:    X - vector

476:    Output Parameter:
477:    X - vector
478:  */
479: int FormInitialGuess(AppCtx *user,Vec X)
480: {
481:   int     i,Nvlocal,ierr;
482:   int     *gloInd;
483:   PetscScalar  *x;
484: #if defined (UNUSED_VARIABLES)
485:   PetscReal  temp1,temp,hx,hy,hxdhy,hydhx,sc;
486:   int     Neglobal,Nvglobal,j,row;
487:   PetscReal  alpha,lambda;

489:   Nvglobal = user->Nvglobal;
490:   Neglobal = user->Neglobal;
491:   lambda   = user->non_lin_param;
492:   alpha    =  user->lin_param;
493: #endif

495:   Nvlocal  = user->Nvlocal;
496:   gloInd   = user->gloInd;

498:   /*
499:      Get a pointer to vector data.
500:        - For default PETSc vectors, VecGetArray() returns a pointer to
501:          the data array.  Otherwise, the routine is implementation dependent.
502:        - You MUST call VecRestoreArray() when you no longer need access to
503:          the array.
504:   */
505:   VecGetArray(X,&x);

507:   /*
508:      Compute initial guess over the locally owned part of the grid
509:   */
510:   for (i=0; i < Nvlocal; i++) {
511:     x[i] = (PetscReal)gloInd[i];
512:   }

514:   /*
515:      Restore vector
516:   */
517:   VecRestoreArray(X,&x);
518:   return 0;
519: }
522: /* --------------------  Evaluate Function F(x) --------------------- */
523: /* 
524:    FormFunction - Evaluates nonlinear function, F(x).

526:    Input Parameters:
527: .  snes - the SNES context
528: .  X - input vector
529: .  ptr - optional user-defined context, as set by SNESSetFunction()

531:    Output Parameter:
532: .  F - function vector
533:  */
534: int FormFunction(SNES snes,Vec X,Vec F,void *ptr)
535: {
536:   AppCtx       *user = (AppCtx*)ptr;
537:   int          ierr,i,j,Nvlocal;
538:   PetscReal       alpha,lambda;
539:   PetscScalar   *x,*f;
540:   VecScatter   scatter;
541:   Vec          localX = user->localX;
542: #if defined (UNUSED_VARIABLES)
543:   PetscScalar  ut,ub,ul,ur,u,*g,sc,uyy,uxx;
544:   PetscReal       hx,hy,hxdhy,hydhx;
545:   PetscReal       two = 2.0,one = 1.0;
546:   int          Nvglobal,Neglobal,row;
547:   int          *gloInd;

549:   Nvglobal = user->Nvglobal;
550:   Neglobal = user->Neglobal;
551:   gloInd   = user->gloInd;
552: #endif

554:   Nvlocal  = user->Nvlocal;
555:   lambda   = user->non_lin_param;
556:   alpha    =  user->lin_param;
557:   scatter  = user->scatter;

559:   /* 
560:      PDE : L(u) + lambda*u*u +alpha*u = 0 where L(u) is the approximate Laplacian as
561:      described in the beginning of this code 
562:                                                                                    
563:      First scatter the distributed vector X into local vector localX (that includes
564:      values for ghost nodes. If we wish,we can put some other work between 
565:      VecScatterBegin() and VecScatterEnd() to overlap the communication with
566:      computation.
567:  */
568:   VecScatterBegin(X,localX,INSERT_VALUES,SCATTER_FORWARD,scatter);
569:   VecScatterEnd(X,localX,INSERT_VALUES,SCATTER_FORWARD,scatter);

571:   /*
572:      Get pointers to vector data
573:   */
574:   VecGetArray(localX,&x);
575:   VecGetArray(F,&f);

577:   /* 
578:     Now compute the f(x). As mentioned earlier, the computed Laplacian is just an 
579:     approximate one chosen for illustrative purpose only. Another point to notice 
580:     is that this is a local (completly parallel) calculation. In practical application 
581:     codes, function calculation time is a dominat portion of the overall execution time.
582:   */
583:   for (i=0; i < Nvlocal; i++) {
584:     f[i] = (user->itot[i] - alpha)*x[i] - lambda*x[i]*x[i];
585:     for (j = 0; j < user->itot[i]; j++) {
586:       f[i] -= x[user->AdjM[i][j]];
587:     }
588:   }

590:   /*
591:      Restore vectors
592:   */
593:   VecRestoreArray(localX,&x);
594:   VecRestoreArray(F,&f);
595:   /*VecView(F,PETSC_VIEWER_STDOUT_WORLD);*/

597:   return 0;
598: }

602: /* --------------------  Evaluate Jacobian F'(x) -------------------- */
603: /*
604:    FormJacobian - Evaluates Jacobian matrix.

606:    Input Parameters:
607: .  snes - the SNES context
608: .  X - input vector
609: .  ptr - optional user-defined context, as set by SNESSetJacobian()

611:    Output Parameters:
612: .  A - Jacobian matrix
613: .  B - optionally different preconditioning matrix
614: .  flag - flag indicating matrix structure

616: */
617: int FormJacobian(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flag,void *ptr)
618: {
619:   AppCtx *user = (AppCtx*)ptr;
620:   Mat     jac = *B;
621:   int     i,j,Nvlocal,col[50],ierr;
622:   PetscScalar  alpha,lambda,value[50];
623:   Vec     localX = user->localX;
624:   VecScatter scatter;
625:   PetscScalar  *x;
626: #if defined (UNUSED_VARIABLES)
627:   PetscScalar  two = 2.0,one = 1.0;
628:   int     row,Nvglobal,Neglobal;
629:   int     *gloInd;

631:   Nvglobal = user->Nvglobal;
632:   Neglobal = user->Neglobal;
633:   gloInd   = user->gloInd;
634: #endif
635: 
636:   /*printf("Entering into FormJacobian \n");*/
637:   Nvlocal  = user->Nvlocal;
638:   lambda   = user->non_lin_param;
639:   alpha    =  user->lin_param;
640:   scatter  = user->scatter;

642:   /* 
643:      PDE : L(u) + lambda*u*u +alpha*u = 0 where L(u) is the approximate Laplacian as
644:      described in the beginning of this code 
645:                                                                                    
646:      First scatter the distributed vector X into local vector localX (that includes
647:      values for ghost nodes. If we wish, we can put some other work between 
648:      VecScatterBegin() and VecScatterEnd() to overlap the communication with
649:      computation.
650:   */
651:   VecScatterBegin(X,localX,INSERT_VALUES,SCATTER_FORWARD,scatter);
652:   VecScatterEnd(X,localX,INSERT_VALUES,SCATTER_FORWARD,scatter);
653: 
654:   /*
655:      Get pointer to vector data
656:   */
657:   VecGetArray(localX,&x);

659:   for (i=0; i < Nvlocal; i++) {
660:     col[0] = i;
661:     value[0] = user->itot[i] - 2.0*lambda*x[i] - alpha;
662:     for (j = 0; j < user->itot[i]; j++) {
663:       col[j+1] = user->AdjM[i][j];
664:       value[j+1] = -1.0;
665:     }

667:   /* 
668:     Set the matrix values in the local ordering. Note that in order to use this
669:     feature we must call the routine MatSetLocalToGlobalMapping() after the 
670:     matrix has been created. 
671:   */
672:     MatSetValuesLocal(jac,1,&i,1+user->itot[i],col,value,INSERT_VALUES);
673:   }

675:   /* 
676:      Assemble matrix, using the 2-step process:
677:        MatAssemblyBegin(), MatAssemblyEnd(). 
678:      Between these two calls, the pointer to vector data has been restored to
679:      demonstrate the use of overlapping communicationn with computation.
680:   */
681:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
682:   VecRestoreArray(localX,&x);
683:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);

685:   /*
686:      Set flag to indicate that the Jacobian matrix retains an identical
687:      nonzero structure throughout all nonlinear iterations (although the
688:      values of the entries change). Thus, we can save some work in setting
689:      up the preconditioner (e.g., no need to redo symbolic factorization for
690:      ILU/ICC preconditioners).
691:       - If the nonzero structure of the matrix is different during
692:         successive linear solves, then the flag DIFFERENT_NONZERO_PATTERN
693:         must be used instead.  If you are unsure whether the matrix
694:         structure has changed or not, use the flag DIFFERENT_NONZERO_PATTERN.
695:       - Caution:  If you specify SAME_NONZERO_PATTERN, PETSc
696:         believes your assertion and does not check the structure
697:         of the matrix.  If you erroneously claim that the structure
698:         is the same when it actually is not, the new preconditioner
699:         will not function correctly.  Thus, use this optimization
700:         feature with caution!
701:   */
702:   *flag = SAME_NONZERO_PATTERN;

704:   /*
705:      Tell the matrix we will never add a new nonzero location to the
706:      matrix. If we do, it will generate an error.
707:   */
708:   MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR);
709:   /* MatView(jac,PETSC_VIEWER_STDOUT_SELF); */
710:   return 0;
711: }
712: #else
713: #include <stdio.h>
714: int main(int argc,char **args)
715: {
716:   fprintf(stdout,"This example does not work for complex numbers.\n");
717:   return 0;
718: }
719: #endif