Actual source code: ex26.c

      2: static char help[] = "Grad-Shafranov solver for one dimensional CHI equilibrium.\n\
      3: The command line options include:\n\
      4:   -n <n> number of grid points\n\
      5:   -psi_axis <axis> \n\
      6:   -r_min <min> \n\
      7:   -param <param> \n\n";

      9: /*T
     10:    Concepts: SNES^parallel CHI equilibrium
     11:    Concepts: DA^using distributed arrays;
     12:    Processors: n
     13: T*/

     15: /* ------------------------------------------------------------------------

     17:    Grad-Shafranov solver for one dimensional CHI equilibrium
     18:
     19:     A finite difference approximation with the usual 3-point stencil
     20:     is used to discretize the boundary value problem to obtain a nonlinear
     21:     system of equations.

     23:     Contributed by Xianzhu Tang, LANL

     25:     An interesting feature of this example is that as you refine the grid
     26:     (with a larger -n <n> you cannot force the residual norm as small. This
     27:     appears to be due to "NOISE" in the function, the FormFunctionLocal() cannot
     28:     be computed as accurately with a finer grid.
     29:   ------------------------------------------------------------------------- */

     31:  #include petscda.h
     32:  #include petscsnes.h

     34: /*
     35:    User-defined application context - contains data needed by the
     36:    application-provided call-back routines, FormJacobian() and
     37:    FormFunction().
     38: */
     39: typedef struct {
     40:   DA            da;               /* distributed array data structure */
     41:   Vec           psi,r;            /* solution, residual vectors */
     42:   Mat           A,J;              /* Jacobian matrix */
     43:   Vec           coordinates;      /* grid coordinates */
     44:   PassiveReal   psi_axis,psi_bdy;
     45:   PassiveReal   r_min;
     46:   PassiveReal   param;            /* test problem parameter */
     47: } AppCtx;

     49: #define GdGdPsi(r,u)      (((r) < 0.05) ? 0.0 : (user->param*((r)-0.05)*(1.0-(u)*(u))*(1.0-(u)*(u))))
     50: #define CurrentWire(r)    (((r) < .05) ? -3.E2 : 0.0)
     51: #define GdGdPsiPrime(r,u) (((r) < 0.05) ? 0.0 : -4.*(user->param*((r)-0.05)*(1.0-(u)*(u)))*u)

     53: /*
     54:    User-defined routines
     55: */

     62: int main(int argc,char **argv)
     63: {
     64:   SNES                   snes;                 /* nonlinear solver */
     65:   AppCtx                 user;                 /* user-defined work context */
     66:   PetscInt               its;                  /* iterations for convergence */
     67:   PetscTruth             fd_jacobian = PETSC_FALSE;
     68:   PetscTruth             adicmf_jacobian = PETSC_FALSE;
     69:   PetscInt               grids = 100, dof = 1, stencil_width = 1;
     70:   PetscErrorCode         ierr;
     71:   PetscReal              fnorm;
     72:   MatFDColoring          matfdcoloring = 0;
     73:   ISColoring             iscoloring;

     75:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     76:      Initialize program
     77:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

     79:   PetscInitialize(&argc,&argv,(char *)0,help);
     80:
     81:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     82:      Initialize problem parameters
     83:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

     85:   user.psi_axis=0.0;
     86:   PetscOptionsGetReal(PETSC_NULL,"-psi_axis",&user.psi_axis,PETSC_NULL);
     87:   user.psi_bdy=1.0;
     88:   PetscOptionsGetReal(PETSC_NULL,"-psi_bdy",&user.psi_bdy,PETSC_NULL);
     89:   user.r_min=0.0;
     90:   PetscOptionsGetReal(PETSC_NULL,"-r_min",&user.r_min,PETSC_NULL);
     91:   user.param=-1.E1;
     92:   PetscOptionsGetReal(PETSC_NULL,"-param",&user.param,PETSC_NULL);

     94:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     95:      Create nonlinear solver context
     96:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     97:   SNESCreate(PETSC_COMM_WORLD,&snes);
     98:
     99:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    100:      Create distributed array (DA) to manage parallel grid and vectors
    101:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    102:   PetscOptionsGetInt(PETSC_NULL,"-n",&grids,PETSC_NULL);
    103:   DACreate1d(PETSC_COMM_WORLD,DA_NONPERIODIC,grids,dof,stencil_width,PETSC_NULL,&user.da);
    104:
    105:   /*  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    106:      Extract global vectors from DA; then duplicate for remaining
    107:      vectors that are the same types
    108:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    109:   DACreateGlobalVector(user.da,&user.psi);
    110:   VecDuplicate(user.psi,&user.r);

    112:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    113:      Set local function evaluation routine
    114:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    116:   DASetLocalFunction(user.da,(DALocalFunction1)FormFunctionLocal);

    118:   SNESSetFunction(snes,user.r,SNESDAFormFunction,&user);

    120:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    121:      Create matrix data structure; set Jacobian evaluation routine

    123:      Set Jacobian matrix data structure and default Jacobian evaluation
    124:      routine. User can override with:
    125:      -snes_mf : matrix-free Newton-Krylov method with no preconditioning
    126:                 (unless user explicitly sets preconditioner)
    127:      -snes_mf_operator : form preconditioning matrix as set by the user,
    128:                          but use matrix-free approx for Jacobian-vector
    129:                          products within Newton-Krylov method

    131:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    132:   PetscOptionsGetTruth(PETSC_NULL,"-fd_jacobian",&fd_jacobian,0);
    133:   /*
    134:        Note that fd_jacobian DOES NOT compute the finite difference approximation to
    135:     the ENTIRE Jacobian. Rather it removes the global coupling from the Jacobian and
    136:     computes the finite difference approximation only for the "local" coupling.

    138:        Thus running with fd_jacobian and not -snes_mf_operator or -adicmf_jacobian
    139:     won't converge.
    140:   */
    141:   if (!fd_jacobian) {
    142:     MatCreate(PETSC_COMM_WORLD,&user.J);
    143:     MatSetSizes(user.J,PETSC_DECIDE,PETSC_DECIDE,grids,grids);
    144:     MatSetType(user.J,MATAIJ);
    145:     MatSetFromOptions(user.J);
    146:     MatSeqAIJSetPreallocation(user.J,5,PETSC_NULL);
    147:     MatMPIAIJSetPreallocation(user.J,5,PETSC_NULL,3,PETSC_NULL);
    148:     user.A    = user.J;
    149:   } else {
    150:     DAGetMatrix(user.da,MATAIJ,&user.J);
    151:     user.A    = user.J;
    152:   }

    154:   PetscOptionsGetTruth(PETSC_NULL,"-adicmf_jacobian",&adicmf_jacobian,0);
    155: #if defined(PETSC_HAVE_ADIC)
    156:   if (adicmf_jacobian) {
    157:     DASetLocalAdicMFFunction(user.da,admf_FormFunctionLocal);
    158:     MatRegisterDAAD();
    159:     MatCreateDAAD(user.da,&user.A);
    160:     MatDAADSetSNES(user.A,snes);
    161:     MatDAADSetCtx(user.A,&user);
    162:   }
    163: #endif

    165:   if (fd_jacobian) {
    166:     DAGetColoring(user.da,IS_COLORING_LOCAL,&iscoloring);
    167:     MatFDColoringCreate(user.J,iscoloring,&matfdcoloring);
    168:     ISColoringDestroy(iscoloring);
    169:     MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))SNESDAFormFunction,&user);
    170:     MatFDColoringSetFromOptions(matfdcoloring);
    171:     SNESSetJacobian(snes,user.A,user.J,SNESDefaultComputeJacobianColor,matfdcoloring);
    172:   } else {
    173:     SNESSetJacobian(snes,user.A,user.J,FormJacobian,&user);
    174:   }

    176:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    177:      Customize nonlinear solver; set runtime options
    178:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    179:   SNESSetFromOptions(snes);

    181:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    182:      Evaluate initial guess
    183:      Note: The user should initialize the vector, x, with the initial guess
    184:      for the nonlinear solver prior to calling SNESSolve().  In particular,
    185:      to employ an initial guess of zero, the user should explicitly set
    186:      this vector to zero by calling VecSet().
    187:   */
    188:   FormInitialGuess(&user,user.psi);

    190:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    191:      Solve nonlinear system
    192:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    193:   SNESSolve(snes,PETSC_NULL,user.psi);
    194:   SNESGetIterationNumber(snes,&its);

    196:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    197:      Explicitly check norm of the residual of the solution
    198:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    199:   SNESDAFormFunction(snes,user.psi,user.r,(void*)&user);
    200:   VecNorm(user.r,NORM_MAX,&fnorm);
    201:   PetscPrintf(PETSC_COMM_WORLD,"Number of Newton iterations = %D fnorm %g\n",its,fnorm);

    203:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    204:      Output the solution vector
    205:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    206:   {
    207:     PetscViewer view_out;
    208:     PetscViewerBinaryOpen(PETSC_COMM_WORLD,"psi.binary",PETSC_FILE_CREATE,&view_out);
    209:     VecView(user.psi,view_out);
    210:     PetscViewerDestroy(view_out);
    211:     PetscViewerASCIIOpen(PETSC_COMM_WORLD,"psi.out",&view_out);
    212:     VecView(user.psi,view_out);
    213:     PetscViewerDestroy(view_out);
    214:   }

    216:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    217:      Free work space.  All PETSc objects should be destroyed when they
    218:      are no longer needed.
    219:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    221:   if (user.A != user.J) {
    222:     MatDestroy(user.A);
    223:   }
    224:   MatDestroy(user.J);
    225:   if (matfdcoloring) {
    226:     MatFDColoringDestroy(matfdcoloring);
    227:   }
    228:   VecDestroy(user.psi);
    229:   VecDestroy(user.r);
    230:   SNESDestroy(snes);
    231:   DADestroy(user.da);
    232:   PetscFinalize();

    234:   return(0);
    235: }
    236: /* ------------------------------------------------------------------- */
    239: /*
    240:    FormInitialGuess - Forms initial approximation.

    242:    Input Parameters:
    243:    user - user-defined application context
    244:    X - vector

    246:    Output Parameter:
    247:    X - vector
    248:  */
    249: PetscErrorCode FormInitialGuess(AppCtx *user,Vec X)
    250: {
    252:   PetscInt       i,Mx,xs,xm;
    253:   PetscScalar    *x;

    256:   DAGetInfo(user->da,PETSC_IGNORE,&Mx,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
    257:                    PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);

    259:   /*
    260:      Get a pointer to vector data.
    261:        - For default PETSc vectors, VecGetArray() returns a pointer to
    262:          the data array.  Otherwise, the routine is implementation dependent.
    263:        - You MUST call VecRestoreArray() when you no longer need access to
    264:          the array.
    265:   */
    266:   DAVecGetArray(user->da,X,&x);

    268:   /*
    269:      Get local grid boundaries (for 2-dimensional DA):
    270:        xs, ys   - starting grid indices (no ghost points)
    271:        xm, ym   - widths of local grid (no ghost points)

    273:   */
    274:   DAGetCorners(user->da,&xs,PETSC_NULL,PETSC_NULL,&xm,PETSC_NULL,PETSC_NULL);

    276:   /*
    277:      Compute initial guess over the locally owned part of the grid
    278:   */
    279:   for (i=xs; i<xs+xm; i++) {
    280:     x[i] = user->psi_axis + i*(user->psi_bdy - user->psi_axis)/(PetscReal)(Mx-1);
    281:   }

    283:   /*
    284:      Restore vector
    285:   */
    286:   DAVecRestoreArray(user->da,X,&x);

    288:   /*
    289:      Check to see if we can import an initial guess from disk
    290:   */
    291:   {
    292:     char         filename[PETSC_MAX_PATH_LEN];
    293:     PetscTruth   flg;
    294:     PetscViewer  view_in;
    295:     PetscReal    fnorm;
    296:     Vec          Y;
    297:     PetscOptionsGetString(PETSC_NULL,"-initialGuess",filename,256,&flg);
    298:     if (flg) {
    299:       PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,PETSC_FILE_RDONLY,&view_in);
    300:       VecLoad(view_in,PETSC_NULL,&Y);
    301:       PetscViewerDestroy(view_in);
    302:       VecMax(Y,PETSC_NULL,&user->psi_bdy);
    303:       SNESDAFormFunction(PETSC_NULL,Y,user->r,(void*)user);
    304:       VecNorm(user->r,NORM_2,&fnorm);
    305:       PetscPrintf(PETSC_COMM_WORLD,"In initial guess: psi_bdy = %f, fnorm = %g.\n",user->psi_bdy,fnorm);
    306:       VecCopy(Y,X);
    307:       VecDestroy(Y);
    308:     }
    309:   }

    311:   return(0);
    312: }
    313: /* ------------------------------------------------------------------- */
    316: /*
    317:    FormFunctionLocal - Evaluates nonlinear function, F(x).
    318:
    319:    Process adiC(36): FormFunctionLocal
    320:
    321: */
    322: PetscErrorCode FormFunctionLocal(DALocalInfo *info,PetscScalar *x,PetscScalar *f,AppCtx *user)
    323: {
    325:   PetscInt       i,xint,xend;
    326:   PetscReal      hx,dhx,r;
    327:   PetscScalar    u,uxx,min = 100000.0,max = -10000.0;
    328:   PetscScalar    psi_0=0.0, psi_a=1.0;
    329:
    331:   for (i=info->xs; i<info->xs + info->xm; i++) {
    332:     if (x[i] > max) max = x[i];
    333:     if (x[i] < min) min = x[i];
    334:   }
    335:   PetscGlobalMax(&max,&psi_a,PETSC_COMM_WORLD);
    336:   PetscGlobalMin(&min,&psi_0,PETSC_COMM_WORLD);

    338:   hx     = 1.0/(PetscReal)(info->mx-1);  dhx    = 1.0/hx;
    339:
    340:   /*
    341:     Compute function over the locally owned part of the grid
    342:   */
    343:   if (info->xs == 0) {
    344:     xint = info->xs + 1; f[0] = (4.*x[1] - x[2] - 3.*x[0])*dhx; /* f[0] = x[0] - user->psi_axis; */
    345:   }
    346:   else {
    347:     xint = info->xs;
    348:   }
    349:   if ((info->xs+info->xm) == info->mx) {
    350:     xend = info->mx - 1; f[info->mx-1] = -(x[info->mx-1] - user->psi_bdy)*dhx;
    351:   }
    352:   else {
    353:     xend = info->xs + info->xm;
    354:   }
    355:
    356:   for (i=xint; i<xend; i++) {
    357:     r       = i*hx + user->r_min;   /* r coordinate */
    358:     u       = (x[i]-psi_0)/(psi_a - psi_0);
    359:     uxx     = ((r+0.5*hx)*(x[i+1]-x[i]) - (r-0.5*hx)*(x[i]-x[i-1]))*dhx; /* r*nabla^2\psi */
    360:     f[i] = uxx  + r*GdGdPsi(r,u)*hx  + r*CurrentWire(r)*hx ;
    361:   }

    363:   PetscLogFlops(11*info->ym*info->xm);
    364:   return(0);
    365: }
    366: /* ------------------------------------------------------------------- */
    369: /*
    370:    FormJacobian - Evaluates Jacobian matrix.

    372:    Input Parameters:
    373: .  snes - the SNES context
    374: .  x - input vector
    375: .  ptr - optional user-defined context, as set by SNESSetJacobian()

    377:    Output Parameters:
    378: .  A - Jacobian matrix
    379: .  B - optionally different preconditioning matrix
    380: .  flag - flag indicating matrix structure

    382: */
    383: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flag,void *ptr)
    384: {
    385:   AppCtx         *user = (AppCtx*)ptr;  /* user-defined application context */
    386:   Mat            jac = *B;                /* Jacobian matrix */
    387:   Vec            localX;
    389:   PetscInt       col[6],row,i,xs,xm,Mx,xint,xend,imin, imax;
    390:   PetscScalar    v[6],hx,dhx,*x;
    391:   PetscReal      r, u, psi_0=0.0, psi_a=1.0;
    392:   PetscTruth     assembled;

    395:   MatAssembled(*B,&assembled);
    396:   if (assembled) {
    397:     MatZeroEntries(*B);
    398:   }

    400:   DAGetLocalVector(user->da,&localX);
    401:   DAGetInfo(user->da,PETSC_IGNORE,&Mx,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
    402:                    PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);

    404:   hx     = 1.0/(PetscReal)(Mx-1);  dhx    = 1.0/hx;

    406:   imin = 0; imax = Mx-1;
    407:   VecMin(X,&imin,&psi_0);
    408:   VecMax(X,&imax,&psi_a);
    409:   PetscPrintf(PETSC_COMM_WORLD,"psi_0(%D)=%g, psi_a(%D)=%g.\n",imin,psi_0,imax,psi_a);

    411:   /*
    412:      Scatter ghost points to local vector, using the 2-step process
    413:         DAGlobalToLocalBegin(), DAGlobalToLocalEnd().
    414:      By placing code between these two statements, computations can be
    415:      done while messages are in transition.
    416:   */
    417:   DAGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX);
    418:   DAGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX);

    420:   /*
    421:      Get pointer to vector data
    422:   */
    423:   DAVecGetArray(user->da,localX,&x);

    425:   /*
    426:      Get local grid boundaries
    427:   */
    428:   DAGetCorners(user->da,&xs,PETSC_NULL,PETSC_NULL,&xm,PETSC_NULL,PETSC_NULL);

    430:   /*
    431:      Compute entries for the locally owned part of the Jacobian.
    432:       - Currently, all PETSc parallel matrix formats are partitioned by
    433:         contiguous chunks of rows across the processors.
    434:       - Each processor needs to insert only elements that it owns
    435:         locally (but any non-local elements will be sent to the
    436:         appropriate processor during matrix assembly).
    437:       - Here, we set all entries for a particular row at once.
    438:       - We can set matrix entries either using either
    439:         MatSetValuesLocal() or MatSetValues(), as discussed above.
    440:   */
    441:   if (xs == 0) {
    442:     xint = xs + 1; /* f[0] = 4.*x[1] - x[2] - 3.*x[0]; */
    443:     row  = 0;     /* first row */
    444:     v[0] = -3.0*dhx;                                              col[0]=row;
    445:     v[1] = 4.0*dhx;                                               col[1]=row+1;
    446:     v[2] = -1.0*dhx;                                              col[2]=row+2;
    447:     MatSetValues(jac,1,&row,3,col,v,ADD_VALUES);
    448:   }
    449:   else {
    450:     xint = xs;
    451:   }
    452:   if ((xs+xm) == Mx) {
    453:     xend = Mx - 1;   /* f[Mx-1] = x[Mx-1] - user->psi_bdy; */
    454:     row  = Mx - 1;  /* last row */
    455:     v[0] = -1.0*dhx;
    456:     MatSetValue(jac,row,row,v[0],ADD_VALUES);
    457:   }
    458:   else {
    459:     xend = xs + xm;
    460:   }

    462:   for (i=xint; i<xend; i++) {
    463:     r       = i*hx + user->r_min;   /* r coordinate */
    464:     u       = (x[i]-psi_0)/(psi_a - psi_0);
    465:     /* uxx     = ((r+0.5*hx)*(x[i+1]-x[i]) - (r-0.5*hx)*(x[i]-x[i-1]))*dhx*dhx; */ /* r*nabla^2\psi */
    466:     row  = i;
    467:     v[0] = (r-0.5*hx)*dhx;                                                              col[0] = i-1;
    468:     v[1] = -2.*r*dhx + hx*r*GdGdPsiPrime(r,u)/(psi_a - psi_0);                          col[1] = i;
    469:     v[2] = (r+0.5*hx)*dhx;                                                              col[2] = i+1;
    470:     v[3] = hx*r*GdGdPsiPrime(r,u)*(x[i] - psi_a)/((psi_a - psi_0)*(psi_a - psi_0));     col[3] = imin;
    471:     v[4] = hx*r*GdGdPsiPrime(r,u)*(psi_0 - x[i])/((psi_a - psi_0)*(psi_a - psi_0));     col[4] = imax;
    472:     MatSetValues(jac,1,&row,5,col,v,ADD_VALUES);
    473:   }
    474:
    475:   DAVecRestoreArray(user->da,localX,&x);
    476:   DARestoreLocalVector(user->da,&localX);

    478:   /*
    479:      Assemble matrix, using the 2-step process:
    480:        MatAssemblyBegin(), MatAssemblyEnd().
    481:   */
    482:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
    483:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);

    485:   /*
    486:      Normally since the matrix has already been assembled above; this
    487:      would do nothing. But in the matrix free mode -snes_mf_operator
    488:      this tells the "matrix-free" matrix that a new linear system solve
    489:      is about to be done.
    490:   */
    491:   MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
    492:   MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);

    494:   /*
    495:      Set flag to indicate that the Jacobian matrix retains an identical
    496:      nonzero structure throughout all nonlinear iterations (although the
    497:      values of the entries change). Thus, we can save some work in setting
    498:      up the preconditioner (e.g., no need to redo symbolic factorization for
    499:      ILU/ICC preconditioners).
    500:       - If the nonzero structure of the matrix is different during
    501:         successive linear solves, then the flag DIFFERENT_NONZERO_PATTERN
    502:         must be used instead.  If you are unsure whether the matrix
    503:         structure has changed or not, use the flag DIFFERENT_NONZERO_PATTERN.
    504:       - Caution:  If you specify SAME_NONZERO_PATTERN, PETSc
    505:         believes your assertion and does not check the structure
    506:         of the matrix.  If you erroneously claim that the structure
    507:         is the same when it actually is not, the new preconditioner
    508:         will not function correctly.  Thus, use this optimization
    509:         feature with caution!
    510:   */
    511:   *flag = SAME_NONZERO_PATTERN;

    513:   /*
    514:      Tell the matrix we will never add a new nonzero location to the
    515:      matrix. If we do, it will generate an error.
    516:   */

    518:   return(0);
    519: }