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:      petscksp.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:   const 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);
434:   SNESGetIterationNumber(snes,&its);
435:   SNESGetNumberUnsuccessfulSteps(snes,&nfails);
436: 
437:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
438:     Print the output : solution vector and other information
439:      Each processor writes to the file output.<rank> where rank is the
440:      processor's rank.
441:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

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

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

464:   return 0;
465: }
468: /* --------------------  Form initial approximation ----------------- */

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

598:   return 0;
599: }

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

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

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

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

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

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

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

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

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

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

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