Actual source code: ex14.c

  2: /* Program usage:  mpiexec -n <procs> ex14 [-help] [all PETSc options] */

  4: static char help[] = "Solves a nonlinear system in parallel with a user-defined Newton method.\n\
  5: Uses KSP to solve the linearized Newton sytems.  This solver\n\
  6: is a very simplistic inexact Newton method.  The intent of this code is to\n\
  7: demonstrate the repeated solution of linear sytems with the same nonzero pattern.\n\
  8: \n\
  9: This is NOT the recommended approach for solving nonlinear problems with PETSc!\n\
 10: We urge users to employ the SNES component for solving nonlinear problems whenever\n\
 11: possible, as it offers many advantages over coding nonlinear solvers independently.\n\
 12: \n\
 13: We solve the  Bratu (SFI - solid fuel ignition) problem in a 2D rectangular\n\
 14: domain, using distributed arrays (DMDAs) to partition the parallel grid.\n\
 15: The command line options include:\n\
 16:   -par <parameter>, where <parameter> indicates the problem's nonlinearity\n\
 17:      problem SFI:  <parameter> = Bratu parameter (0 <= par <= 6.81)\n\
 18:   -mx <xg>, where <xg> = number of grid points in the x-direction\n\
 19:   -my <yg>, where <yg> = number of grid points in the y-direction\n\
 20:   -Nx <npx>, where <npx> = number of processors in the x-direction\n\
 21:   -Ny <npy>, where <npy> = number of processors in the y-direction\n\n";

 23: /*T
 24:    Concepts: KSP^writing a user-defined nonlinear solver (parallel Bratu example);
 25:    Concepts: DMDA^using distributed arrays;
 26:    Processors: n
 27: T*/

 29: /* ------------------------------------------------------------------------

 31:     Solid Fuel Ignition (SFI) problem.  This problem is modeled by
 32:     the partial differential equation
 33:   
 34:             -Laplacian u - lambda*exp(u) = 0,  0 < x,y < 1,
 35:   
 36:     with boundary conditions
 37:    
 38:              u = 0  for  x = 0, x = 1, y = 0, y = 1.
 39:   
 40:     A finite difference approximation with the usual 5-point stencil
 41:     is used to discretize the boundary value problem to obtain a nonlinear 
 42:     system of equations.

 44:     The SNES version of this problem is:  snes/examples/tutorials/ex5.c
 45:     We urge users to employ the SNES component for solving nonlinear
 46:     problems whenever possible, as it offers many advantages over coding 
 47:     nonlinear solvers independently.

 49:   ------------------------------------------------------------------------- */

 51: /* 
 52:    Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
 53:    Include "petscksp.h" so that we can use KSP solvers.  Note that this
 54:    file automatically includes:
 55:      petscsys.h       - base PETSc routines   petscvec.h - vectors
 56:      petscmat.h - matrices
 57:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
 58:      petscviewer.h - viewers               petscpc.h  - preconditioners
 59: */
 60: #include <petscdmda.h>
 61: #include <petscksp.h>

 63: /* 
 64:    User-defined application context - contains data needed by the 
 65:    application-provided call-back routines, ComputeJacobian() and
 66:    ComputeFunction().
 67: */
 68: typedef struct {
 69:    PetscReal   param;          /* test problem parameter */
 70:    PetscInt    mx,my;          /* discretization in x,y directions */
 71:    Vec         localX,localF; /* ghosted local vector */
 72:    DM          da;             /* distributed array data structure */
 73:    PetscInt    rank;           /* processor rank */
 74: } AppCtx;

 76: /* 
 77:    User-defined routines
 78: */

 84: int main(int argc,char **argv)
 85: {
 86:   /* -------------- Data to define application problem ---------------- */
 87:   MPI_Comm       comm;                /* communicator */
 88:   KSP            ksp;                /* linear solver */
 89:   Vec            X,Y,F;             /* solution, update, residual vectors */
 90:   Mat            J;                   /* Jacobian matrix */
 91:   AppCtx         user;                /* user-defined work context */
 92:   PetscInt       Nx,Ny;              /* number of preocessors in x- and y- directions */
 93:   PetscMPIInt    size;                /* number of processors */
 94:   PetscReal      bratu_lambda_max = 6.81,bratu_lambda_min = 0.;
 95:   PetscInt       m,N;

 98:   /* --------------- Data to define nonlinear solver -------------- */
 99:   PetscReal      rtol = 1.e-8;        /* relative convergence tolerance */
100:   PetscReal      xtol = 1.e-8;        /* step convergence tolerance */
101:   PetscReal      ttol;                /* convergence tolerance */
102:   PetscReal      fnorm,ynorm,xnorm; /* various vector norms */
103:   PetscInt       max_nonlin_its = 10; /* maximum number of iterations for nonlinear solver */
104:   PetscInt       max_functions = 50;  /* maximum number of function evaluations */
105:   PetscInt       lin_its;             /* number of linear solver iterations for each step */
106:   PetscInt       i;                   /* nonlinear solve iteration number */
107:   MatStructure   mat_flag;        /* flag indicating structure of preconditioner matrix */
108:   PetscBool      no_output = PETSC_FALSE;           /* flag indicating whether to surpress output */

110:   PetscInitialize(&argc,&argv,(char *)0,help);
111:   comm = PETSC_COMM_WORLD;
112:   MPI_Comm_rank(comm,&user.rank);
113:   PetscOptionsGetBool(PETSC_NULL,"-no_output",&no_output,PETSC_NULL);

115:   /*
116:      Initialize problem parameters
117:   */
118:   user.mx = 4; user.my = 4; user.param = 6.0;
119:   PetscOptionsGetInt(PETSC_NULL,"-mx",&user.mx,PETSC_NULL);
120:   PetscOptionsGetInt(PETSC_NULL,"-my",&user.my,PETSC_NULL);
121:   PetscOptionsGetReal(PETSC_NULL,"-par",&user.param,PETSC_NULL);
122:   if (user.param >= bratu_lambda_max || user.param <= bratu_lambda_min) {
123:     SETERRQ(PETSC_COMM_WORLD,1,"Lambda is out of range");
124:   }
125:   N = user.mx*user.my;

127:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
128:      Create linear solver context
129:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

131:   KSPCreate(comm,&ksp);

133:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
134:      Create vector data structures
135:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

137:   /*
138:      Create distributed array (DMDA) to manage parallel grid and vectors
139:   */
140:   MPI_Comm_size(comm,&size);
141:   Nx = PETSC_DECIDE; Ny = PETSC_DECIDE;
142:   PetscOptionsGetInt(PETSC_NULL,"-Nx",&Nx,PETSC_NULL);
143:   PetscOptionsGetInt(PETSC_NULL,"-Ny",&Ny,PETSC_NULL);
144:   if (Nx*Ny != size && (Nx != PETSC_DECIDE || Ny != PETSC_DECIDE))
145:     SETERRQ(PETSC_COMM_WORLD,1,"Incompatible number of processors:  Nx * Ny != size");
146:   DMDACreate2d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,user.mx,
147:                     user.my,Nx,Ny,1,1,PETSC_NULL,PETSC_NULL,&user.da);

149:   /*
150:      Extract global and local vectors from DMDA; then duplicate for remaining
151:      vectors that are the same types
152:   */
153:   DMCreateGlobalVector(user.da,&X);
154:   DMCreateLocalVector(user.da,&user.localX);
155:   VecDuplicate(X,&F);
156:   VecDuplicate(X,&Y);
157:   VecDuplicate(user.localX,&user.localF);


160:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
161:      Create matrix data structure for Jacobian
162:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
163:   /*
164:      Note:  For the parallel case, vectors and matrices MUST be partitioned
165:      accordingly.  When using distributed arrays (DMDAs) to create vectors,
166:      the DMDAs determine the problem partitioning.  We must explicitly
167:      specify the local matrix dimensions upon its creation for compatibility
168:      with the vector distribution.  Thus, the generic MatCreate() routine
169:      is NOT sufficient when working with distributed arrays.

171:      Note: Here we only approximately preallocate storage space for the
172:      Jacobian.  See the users manual for a discussion of better techniques
173:      for preallocating matrix memory.
174:   */
175:   if (size == 1) {
176:     MatCreateSeqAIJ(comm,N,N,5,PETSC_NULL,&J);
177:   } else {
178:     VecGetLocalSize(X,&m);
179:     MatCreateMPIAIJ(comm,m,m,N,N,5,PETSC_NULL,3,PETSC_NULL,&J);
180:   }

182:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
183:      Customize linear solver; set runtime options
184:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

186:   /*
187:      Set runtime options (e.g.,-ksp_monitor -ksp_rtol <rtol> -ksp_type <type>)
188:   */
189:   KSPSetFromOptions(ksp);

191:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
192:      Evaluate initial guess
193:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

195:   FormInitialGuess(&user,X);
196:   ComputeFunction(&user,X,F);   /* Compute F(X)    */
197:   VecNorm(F,NORM_2,&fnorm);     /* fnorm = || F || */
198:   ttol = fnorm*rtol;
199:   if (!no_output) PetscPrintf(comm,"Initial function norm = %G\n",fnorm);

201:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
202:      Solve nonlinear system with a user-defined method
203:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

205:   /* 
206:       This solver is a very simplistic inexact Newton method, with no
207:       no damping strategies or bells and whistles. The intent of this code
208:       is  merely to demonstrate the repeated solution with KSP of linear
209:       sytems with the same nonzero structure.

211:       This is NOT the recommended approach for solving nonlinear problems
212:       with PETSc!  We urge users to employ the SNES component for solving
213:       nonlinear problems whenever possible with application codes, as it
214:       offers many advantages over coding nonlinear solvers independently.
215:    */

217:   for (i=0; i<max_nonlin_its; i++) {

219:     /* 
220:         Compute the Jacobian matrix.  See the comments in this routine for
221:         important information about setting the flag mat_flag.
222:      */
223:     ComputeJacobian(&user,X,J,&mat_flag);

225:     /* 
226:         Solve J Y = F, where J is the Jacobian matrix.
227:           - First, set the KSP linear operators.  Here the matrix that
228:             defines the linear system also serves as the preconditioning
229:             matrix.
230:           - Then solve the Newton system.
231:      */
232:     KSPSetOperators(ksp,J,J,mat_flag);
233:     KSPSolve(ksp,F,Y);
234:     KSPGetIterationNumber(ksp,&lin_its);

236:     /* 
237:        Compute updated iterate
238:      */
239:     VecNorm(Y,NORM_2,&ynorm);       /* ynorm = || Y || */
240:     VecAYPX(Y,-1.0,X);              /* Y <- X - Y      */
241:     VecCopy(Y,X);                   /* X <- Y          */
242:     VecNorm(X,NORM_2,&xnorm);       /* xnorm = || X || */
243:     if (!no_output) {
244:       PetscPrintf(comm,"   linear solve iterations = %D, xnorm=%G, ynorm=%G\n",lin_its,xnorm,ynorm);
245:     }

247:     /* 
248:        Evaluate new nonlinear function
249:      */
250:     ComputeFunction(&user,X,F);     /* Compute F(X)    */
251:     VecNorm(F,NORM_2,&fnorm);       /* fnorm = || F || */
252:     if (!no_output) {
253:       PetscPrintf(comm,"Iteration %D, function norm = %G\n",i+1,fnorm);
254:     }

256:     /*
257:        Test for convergence
258:      */
259:     if (fnorm <= ttol) {
260:       if (!no_output) {
261:          PetscPrintf(comm,"Converged due to function norm %G < %G (relative tolerance)\n",fnorm,ttol);
262:       }
263:       break;
264:     }
265:     if (ynorm < xtol*(xnorm)) {
266:       if (!no_output) {
267:          PetscPrintf(comm,"Converged due to small update length: %G < %G * %G\n",ynorm,xtol,xnorm);
268:       }
269:       break;
270:     }
271:     if (i > max_functions) {
272:       if (!no_output) {
273:         PetscPrintf(comm,"Exceeded maximum number of function evaluations: %D > %D\n",i,max_functions);
274:       }
275:       break;
276:     }
277:   }
278:   PetscPrintf(comm,"Number of Newton iterations = %D\n",i+1);

280:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
281:      Free work space.  All PETSc objects should be destroyed when they
282:      are no longer needed.
283:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

285:   MatDestroy(&J);           VecDestroy(&Y);
286:   VecDestroy(&user.localX); VecDestroy(&X);
287:   VecDestroy(&user.localF); VecDestroy(&F);
288:   KSPDestroy(&ksp);  DMDestroy(&user.da);
289:   PetscFinalize();

291:   return 0;
292: }
293: /* ------------------------------------------------------------------- */
296: /* 
297:    FormInitialGuess - Forms initial approximation.

299:    Input Parameters:
300:    user - user-defined application context
301:    X - vector

303:    Output Parameter:
304:    X - vector
305:  */
306: PetscErrorCode FormInitialGuess(AppCtx *user,Vec X)
307: {
308:   PetscInt     i,j,row,mx,my,ierr,xs,ys,xm,ym,gxm,gym,gxs,gys;
309:   PetscReal    one = 1.0,lambda,temp1,temp,hx,hy;
310:   PetscScalar  *x;
311:   Vec          localX = user->localX;

313:   mx = user->mx;            my = user->my;            lambda = user->param;
314:   hx = one/(PetscReal)(mx-1);  hy = one/(PetscReal)(my-1);
315:   temp1 = lambda/(lambda + one);

317:   /*
318:      Get a pointer to vector data.
319:        - For default PETSc vectors, VecGetArray() returns a pointer to
320:          the data array.  Otherwise, the routine is implementation dependent.
321:        - You MUST call VecRestoreArray() when you no longer need access to
322:          the array.
323:   */
324:   VecGetArray(localX,&x);

326:   /*
327:      Get local grid boundaries (for 2-dimensional DMDA):
328:        xs, ys   - starting grid indices (no ghost points)
329:        xm, ym   - widths of local grid (no ghost points)
330:        gxs, gys - starting grid indices (including ghost points)
331:        gxm, gym - widths of local grid (including ghost points)
332:   */
333:   DMDAGetCorners(user->da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
334:   DMDAGetGhostCorners(user->da,&gxs,&gys,PETSC_NULL,&gxm,&gym,PETSC_NULL);

336:   /*
337:      Compute initial guess over the locally owned part of the grid
338:   */
339:   for (j=ys; j<ys+ym; j++) {
340:     temp = (PetscReal)(PetscMin(j,my-j-1))*hy;
341:     for (i=xs; i<xs+xm; i++) {
342:       row = i - gxs + (j - gys)*gxm;
343:       if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
344:         x[row] = 0.0;
345:         continue;
346:       }
347:       x[row] = temp1*PetscSqrtScalar(PetscMin((PetscReal)(PetscMin(i,mx-i-1))*hx,temp));
348:     }
349:   }

351:   /*
352:      Restore vector
353:   */
354:   VecRestoreArray(localX,&x);

356:   /*
357:      Insert values into global vector
358:   */
359:   DMLocalToGlobalBegin(user->da,localX,INSERT_VALUES,X);
360:   DMLocalToGlobalEnd(user->da,localX,INSERT_VALUES,X);
361:   return 0;
362: }
363: /* ------------------------------------------------------------------- */
366: /* 
367:    ComputeFunction - Evaluates nonlinear function, F(x).

369:    Input Parameters:
370: .  X - input vector
371: .  user - user-defined application context

373:    Output Parameter:
374: .  F - function vector
375:  */
376: PetscErrorCode ComputeFunction(AppCtx *user,Vec X,Vec F)
377: {
379:   PetscInt       i,j,row,mx,my,xs,ys,xm,ym,gxs,gys,gxm,gym;
380:   PetscReal      two = 2.0,one = 1.0,lambda,hx,hy,hxdhy,hydhx,sc;
381:   PetscScalar    u,uxx,uyy,*x,*f;
382:   Vec            localX = user->localX,localF = user->localF;

384:   mx = user->mx;            my = user->my;            lambda = user->param;
385:   hx = one/(PetscReal)(mx-1);  hy = one/(PetscReal)(my-1);
386:   sc = hx*hy*lambda;        hxdhy = hx/hy;            hydhx = hy/hx;

388:   /*
389:      Scatter ghost points to local vector, using the 2-step process
390:         DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
391:      By placing code between these two statements, computations can be
392:      done while messages are in transition.
393:   */
394:   DMGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX);
395:   DMGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX);

397:   /*
398:      Get pointers to vector data
399:   */
400:   VecGetArray(localX,&x);
401:   VecGetArray(localF,&f);

403:   /*
404:      Get local grid boundaries
405:   */
406:   DMDAGetCorners(user->da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
407:   DMDAGetGhostCorners(user->da,&gxs,&gys,PETSC_NULL,&gxm,&gym,PETSC_NULL);

409:   /*
410:      Compute function over the locally owned part of the grid
411:   */
412:   for (j=ys; j<ys+ym; j++) {
413:     row = (j - gys)*gxm + xs - gxs - 1;
414:     for (i=xs; i<xs+xm; i++) {
415:       row++;
416:       if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
417:         f[row] = x[row];
418:         continue;
419:       }
420:       u = x[row];
421:       uxx = (two*u - x[row-1] - x[row+1])*hydhx;
422:       uyy = (two*u - x[row-gxm] - x[row+gxm])*hxdhy;
423:       f[row] = uxx + uyy - sc*PetscExpScalar(u);
424:     }
425:   }

427:   /*
428:      Restore vectors
429:   */
430:   VecRestoreArray(localX,&x);
431:   VecRestoreArray(localF,&f);

433:   /*
434:      Insert values into global vector
435:   */
436:   DMLocalToGlobalBegin(user->da,localF,INSERT_VALUES,F);
437:   DMLocalToGlobalEnd(user->da,localF,INSERT_VALUES,F);
438:   PetscLogFlops(11.0*ym*xm);
439:   return 0;
440: }
441: /* ------------------------------------------------------------------- */
444: /*
445:    ComputeJacobian - Evaluates Jacobian matrix.

447:    Input Parameters:
448: .  x - input vector
449: .  user - user-defined application context

451:    Output Parameters:
452: .  jac - Jacobian matrix
453: .  flag - flag indicating matrix structure

455:    Notes:
456:    Due to grid point reordering with DMDAs, we must always work
457:    with the local grid points, and then transform them to the new
458:    global numbering with the "ltog" mapping (via DMDAGetGlobalIndices()).
459:    We cannot work directly with the global numbers for the original
460:    uniprocessor grid!
461: */
462: PetscErrorCode ComputeJacobian(AppCtx *user,Vec X,Mat jac,MatStructure *flag)
463: {
465:   Vec            localX = user->localX;   /* local vector */
466:   PetscInt       *ltog;                   /* local-to-global mapping */
467:   PetscInt       i,j,row,mx,my,col[5];
468:   PetscInt       nloc,xs,ys,xm,ym,gxs,gys,gxm,gym,grow;
469:   PetscScalar    two = 2.0,one = 1.0,lambda,v[5],hx,hy,hxdhy,hydhx,sc,*x;

471:   mx = user->mx;            my = user->my;            lambda = user->param;
472:   hx = one/(PetscReal)(mx-1);  hy = one/(PetscReal)(my-1);
473:   sc = hx*hy;               hxdhy = hx/hy;            hydhx = hy/hx;

475:   /*
476:      Scatter ghost points to local vector, using the 2-step process
477:         DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
478:      By placing code between these two statements, computations can be
479:      done while messages are in transition.
480:   */
481:   DMGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX);
482:   DMGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX);

484:   /*
485:      Get pointer to vector data
486:   */
487:   VecGetArray(localX,&x);

489:   /*
490:      Get local grid boundaries
491:   */
492:   DMDAGetCorners(user->da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
493:   DMDAGetGhostCorners(user->da,&gxs,&gys,PETSC_NULL,&gxm,&gym,PETSC_NULL);

495:   /*
496:      Get the global node numbers for all local nodes, including ghost points
497:   */
498:   DMDAGetGlobalIndices(user->da,&nloc,&ltog);

500:   /* 
501:      Compute entries for the locally owned part of the Jacobian.
502:       - Currently, all PETSc parallel matrix formats are partitioned by
503:         contiguous chunks of rows across the processors. The "grow"
504:         parameter computed below specifies the global row number 
505:         corresponding to each local grid point.
506:       - Each processor needs to insert only elements that it owns
507:         locally (but any non-local elements will be sent to the
508:         appropriate processor during matrix assembly). 
509:       - Always specify global row and columns of matrix entries.
510:       - Here, we set all entries for a particular row at once.
511:   */
512:   for (j=ys; j<ys+ym; j++) {
513:     row = (j - gys)*gxm + xs - gxs - 1;
514:     for (i=xs; i<xs+xm; i++) {
515:       row++;
516:       grow = ltog[row];
517:       /* boundary points */
518:       if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
519:         MatSetValues(jac,1,&grow,1,&grow,&one,INSERT_VALUES);
520:         continue;
521:       }
522:       /* interior grid points */
523:       v[0] = -hxdhy; col[0] = ltog[row - gxm];
524:       v[1] = -hydhx; col[1] = ltog[row - 1];
525:       v[2] = two*(hydhx + hxdhy) - sc*lambda*PetscExpScalar(x[row]); col[2] = grow;
526:       v[3] = -hydhx; col[3] = ltog[row + 1];
527:       v[4] = -hxdhy; col[4] = ltog[row + gxm];
528:       MatSetValues(jac,1,&grow,5,col,v,INSERT_VALUES);
529:     }
530:   }

532:   /* 
533:      Assemble matrix, using the 2-step process:
534:        MatAssemblyBegin(), MatAssemblyEnd().
535:      By placing code between these two statements, computations can be
536:      done while messages are in transition.
537:   */
538:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
539:   VecRestoreArray(localX,&x);
540:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);

542:   /*
543:      Set flag to indicate that the Jacobian matrix retains an identical
544:      nonzero structure throughout all nonlinear iterations (although the
545:      values of the entries change). Thus, we can save some work in setting
546:      up the preconditioner (e.g., no need to redo symbolic factorization for
547:      ILU/ICC preconditioners).
548:       - If the nonzero structure of the matrix is different during
549:         successive linear solves, then the flag DIFFERENT_NONZERO_PATTERN
550:         must be used instead.  If you are unsure whether the matrix
551:         structure has changed or not, use the flag DIFFERENT_NONZERO_PATTERN.
552:       - Caution:  If you specify SAME_NONZERO_PATTERN, PETSc
553:         believes your assertion and does not check the structure
554:         of the matrix.  If you erroneously claim that the structure
555:         is the same when it actually is not, the new preconditioner
556:         will not function correctly.  Thus, use this optimization
557:         feature with caution!
558:   */
559:   *flag = SAME_NONZERO_PATTERN;
560:   return 0;
561: }