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

      1: static const char help[] = "Newton's method to solve a two-variable system, sequentially.\n"
      2:                            "The same problem is solved twice - i) fully assembled system + ii) block system\n\n";

      4: /*T
      5: Concepts: SNES^basic uniprocessor example, block objects
      6: Processors: 1
      7: T*/

      9: /*
     10: Include "petscsnes.h" so that we can use SNES solvers.  Note that this
     11: file automatically includes:
     12: petscsys.h       - base PETSc routines   petscvec.h - vectors
     13: petscsys.h    - system routines       petscmat.h - matrices
     14: petscis.h     - index sets            petscksp.h - Krylov subspace methods
     15: petscviewer.h - viewers               petscpc.h  - preconditioners
     16: petscksp.h   - linear solvers
     17: */
     18: #include <petscsnes.h>

     20: /*
     21: This example is block version of the test found at
     22:   ${PETSC_DIR}/src/snes/examples/tutorials/ex1.c
     23: In this test we replace the Jacobian systems
     24:   [J]{x} = {F}
     25: where

     27: [J] = (j_00, j_01),  {x} = (x_0, x_1)^T,   {F} = (f_0, f_1)^T
     28:       (j_10, j_11)
     29: where [J] \in \mathbb^{2 \times 2}, {x},{F} \in \mathbb^{2 \times 1},

     31: with a block system in which each block is of length 1.
     32: i.e. The block system is thus

     34: [J] = ([j00], [j01]),  {x} = ({x0}, {x1})^T, {F} = ({f0}, {f1})^T
     35:       ([j10], [j11])
     36: where
     37: [j00], [j01], [j10], [j11] \in \mathbb^{1 \times 1}
     38: {x0}, {x1}, {f0}, {f1} \in \mathbb^{1 \times 1}

     40: In practice we would not bother defing blocks of size one, and would instead assemble the
     41: full system. This is just a simple test to illustrate how to manipulate the blocks and
     42: to confirm the implementation is correct.
     43: */

     45: /*
     46: User-defined routines
     47: */
     48: static PetscErrorCode FormJacobian1(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
     49: static PetscErrorCode FormFunction1(SNES,Vec,Vec,void*);
     50: static PetscErrorCode FormJacobian2(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
     51: static PetscErrorCode FormFunction2(SNES,Vec,Vec,void*);
     52: static PetscErrorCode FormJacobian1_block(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
     53: static PetscErrorCode FormFunction1_block(SNES,Vec,Vec,void*);
     54: static PetscErrorCode FormJacobian2_block(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
     55: static PetscErrorCode FormFunction2_block(SNES,Vec,Vec,void*);

     60: static PetscErrorCode assembled_system(void)
     61: {
     62:   SNES           snes;         /* nonlinear solver context */
     63:   KSP            ksp;         /* linear solver context */
     64:   PC             pc;           /* preconditioner context */
     65:   Vec            x,r;         /* solution, residual vectors */
     66:   Mat            J;            /* Jacobian matrix */
     68:   PetscInt       its;
     69:   PetscScalar    pfive = .5,*xx;
     70:   PetscBool      flg;

     73:   PetscPrintf(PETSC_COMM_WORLD, "\n\n========================= Assembled system =========================\n\n");

     75:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     76:   Create nonlinear solver context
     77:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

     79:   SNESCreate(PETSC_COMM_WORLD,&snes);

     81:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     82:   Create matrix and vector data structures; set corresponding routines
     83:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

     85:   /*
     86:   Create vectors for solution and nonlinear function
     87:   */
     88:   VecCreateSeq(PETSC_COMM_SELF,2,&x);
     89:   VecDuplicate(x,&r);

     91:   /*
     92:   Create Jacobian matrix data structure
     93:   */
     94:   MatCreate(PETSC_COMM_SELF,&J);
     95:   MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,2,2);
     96:   MatSetFromOptions(J);
     97:   MatSetUp(J);

     99:   PetscOptionsHasName(NULL,"-hard",&flg);
    100:   if (!flg) {
    101:     /*
    102:     Set function evaluation routine and vector.
    103:     */
    104:     SNESSetFunction(snes,r,FormFunction1,NULL);

    106:     /*
    107:     Set Jacobian matrix data structure and Jacobian evaluation routine
    108:     */
    109:     SNESSetJacobian(snes,J,J,FormJacobian1,NULL);
    110:   } else {
    111:     SNESSetFunction(snes,r,FormFunction2,NULL);
    112:     SNESSetJacobian(snes,J,J,FormJacobian2,NULL);
    113:   }

    115:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    116:   Customize nonlinear solver; set runtime options
    117:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    119:   /*
    120:   Set linear solver defaults for this problem. By extracting the
    121:   KSP, KSP, and PC contexts from the SNES context, we can then
    122:   directly call any KSP, KSP, and PC routines to set various options.
    123:   */
    124:   SNESGetKSP(snes,&ksp);
    125:   KSPGetPC(ksp,&pc);
    126:   PCSetType(pc,PCNONE);
    127:   KSPSetTolerances(ksp,1.e-4,PETSC_DEFAULT,PETSC_DEFAULT,20);

    129:   /*
    130:   Set SNES/KSP/KSP/PC runtime options, e.g.,
    131:   -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
    132:   These options will override those specified above as long as
    133:   SNESSetFromOptions() is called _after_ any other customization
    134:   routines.
    135:   */
    136:   SNESSetFromOptions(snes);

    138:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    139:   Evaluate initial guess; then solve nonlinear system
    140:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    141:   if (!flg) {
    142:     VecSet(x,pfive);
    143:   } else {
    144:     VecGetArray(x,&xx);
    145:     xx[0] = 2.0; xx[1] = 3.0;
    146:     VecRestoreArray(x,&xx);
    147:   }
    148:   /*
    149:   Note: The user should initialize the vector, x, with the initial guess
    150:   for the nonlinear solver prior to calling SNESSolve().  In particular,
    151:   to employ an initial guess of zero, the user should explicitly set
    152:   this vector to zero by calling VecSet().
    153:   */

    155:   SNESSolve(snes,NULL,x);
    156:   SNESGetIterationNumber(snes,&its);
    157:   if (flg) {
    158:     Vec f;
    159:     VecView(x,PETSC_VIEWER_STDOUT_WORLD);
    160:     SNESGetFunction(snes,&f,0,0);
    161:     VecView(r,PETSC_VIEWER_STDOUT_WORLD);
    162:   }

    164:   PetscPrintf(PETSC_COMM_SELF,"number of SNES iterations = %D\n\n",its);

    166:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    167:   Free work space.  All PETSc objects should be destroyed when they
    168:   are no longer needed.
    169:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    171:   VecDestroy(&x); VecDestroy(&r);
    172:   MatDestroy(&J); SNESDestroy(&snes);
    173:   return(0);
    174: }

    176: /* ------------------------------------------------------------------- */
    179: /*
    180: FormFunction1 - Evaluates nonlinear function, F(x).

    182: Input Parameters:
    183: .  snes - the SNES context
    184: .  x - input vector
    185: .  dummy - optional user-defined context (not used here)

    187: Output Parameter:
    188: .  f - function vector
    189: */
    190: static PetscErrorCode FormFunction1(SNES snes,Vec x,Vec f,void *dummy)
    191: {
    193:   PetscScalar    *xx,*ff;

    196:   /*
    197:   Get pointers to vector data.
    198:   - For default PETSc vectors, VecGetArray() returns a pointer to
    199:   the data array.  Otherwise, the routine is implementation dependent.
    200:   - You MUST call VecRestoreArray() when you no longer need access to
    201:   the array.
    202:   */
    203:   VecGetArray(x,&xx);
    204:   VecGetArray(f,&ff);

    206:   /*
    207:   Compute function
    208:   */
    209:   ff[0] = xx[0]*xx[0] + xx[0]*xx[1] - 3.0;
    210:   ff[1] = xx[0]*xx[1] + xx[1]*xx[1] - 6.0;

    213:   /*
    214:   Restore vectors
    215:   */
    216:   VecRestoreArray(x,&xx);
    217:   VecRestoreArray(f,&ff);
    218:   return(0);
    219: }
    220: /* ------------------------------------------------------------------- */
    223: /*
    224: FormJacobian1 - Evaluates Jacobian matrix.

    226: Input Parameters:
    227: .  snes - the SNES context
    228: .  x - input vector
    229: .  dummy - optional user-defined context (not used here)

    231: Output Parameters:
    232: .  jac - Jacobian matrix
    233: .  B - optionally different preconditioning matrix
    234: .  flag - flag indicating matrix structure
    235: */
    236: static PetscErrorCode FormJacobian1(SNES snes,Vec x,Mat *jac,Mat *B,MatStructure *flag,void *dummy)
    237: {
    238:   PetscScalar    *xx,A[4];
    240:   PetscInt       idx[2] = {0,1};

    243:   /*
    244:   Get pointer to vector data
    245:   */
    246:   VecGetArray(x,&xx);

    248:   /*
    249:   Compute Jacobian entries and insert into matrix.
    250:   - Since this is such a small problem, we set all entries for
    251:   the matrix at once.
    252:   */
    253:   A[0]  = 2.0*xx[0] + xx[1]; A[1] = xx[0];
    254:   A[2]  = xx[1]; A[3] = xx[0] + 2.0*xx[1];
    255:   MatSetValues(*jac,2,idx,2,idx,A,INSERT_VALUES);
    256:   *flag = SAME_NONZERO_PATTERN;

    258:   /*
    259:   Restore vector
    260:   */
    261:   VecRestoreArray(x,&xx);

    263:   /*
    264:   Assemble matrix
    265:   */
    266:   MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
    267:   MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);
    268:   return(0);
    269: }

    272: /* ------------------------------------------------------------------- */
    275: static PetscErrorCode FormFunction2(SNES snes,Vec x,Vec f,void *dummy)
    276: {
    278:   PetscScalar    *xx,*ff;

    281:   /*
    282:   Get pointers to vector data.
    283:   - For default PETSc vectors, VecGetArray() returns a pointer to
    284:   the data array.  Otherwise, the routine is implementation dependent.
    285:   - You MUST call VecRestoreArray() when you no longer need access to
    286:   the array.
    287:   */
    288:   VecGetArray(x,&xx);
    289:   VecGetArray(f,&ff);

    291:   /*
    292:   Compute function
    293:   */
    294:   ff[0] = PetscSinScalar(3.0*xx[0]) + xx[0];
    295:   ff[1] = xx[1];

    297:   /*
    298:   Restore vectors
    299:   */
    300:   VecRestoreArray(x,&xx);
    301:   VecRestoreArray(f,&ff);
    302:   return(0);
    303: }

    305: /* ------------------------------------------------------------------- */
    308: static PetscErrorCode FormJacobian2(SNES snes,Vec x,Mat *jac,Mat *B,MatStructure *flag,void *dummy)
    309: {
    310:   PetscScalar    *xx,A[4];
    312:   PetscInt       idx[2] = {0,1};

    315:   /*
    316:   Get pointer to vector data
    317:   */
    318:   VecGetArray(x,&xx);

    320:   /*
    321:   Compute Jacobian entries and insert into matrix.
    322:   - Since this is such a small problem, we set all entries for
    323:   the matrix at once.
    324:   */
    325:   A[0]  = 3.0*PetscCosScalar(3.0*xx[0]) + 1.0; A[1] = 0.0;
    326:   A[2]  = 0.0;                     A[3] = 1.0;
    327:   MatSetValues(*jac,2,idx,2,idx,A,INSERT_VALUES);
    328:   *flag = SAME_NONZERO_PATTERN;

    330:   /*
    331:   Restore vector
    332:   */
    333:   VecRestoreArray(x,&xx);

    335:   /*
    336:   Assemble matrix
    337:   */
    338:   MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
    339:   MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);
    340:   return(0);
    341: }

    345: static int block_system(void)
    346: {
    347:   SNES           snes;         /* nonlinear solver context */
    348:   KSP            ksp;         /* linear solver context */
    349:   PC             pc;           /* preconditioner context */
    350:   Vec            x,r;         /* solution, residual vectors */
    351:   Mat            J;            /* Jacobian matrix */
    353:   PetscInt       its;
    354:   PetscScalar    pfive = .5;
    355:   PetscBool      flg;

    357:   Mat j11, j12, j21, j22;
    358:   Vec x1, x2, r1, r2;
    359:   Vec bv;
    360:   Vec bx[2];
    361:   Mat bA[2][2];

    364:   PetscPrintf(PETSC_COMM_WORLD, "\n\n========================= Block system =========================\n\n");

    366:   SNESCreate(PETSC_COMM_WORLD,&snes);

    368:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    369:   Create matrix and vector data structures; set corresponding routines
    370:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    372:   /*
    373:   Create sub vectors for solution and nonlinear function
    374:   */
    375:   VecCreateSeq(PETSC_COMM_SELF,1,&x1);
    376:   VecDuplicate(x1,&r1);

    378:   VecCreateSeq(PETSC_COMM_SELF,1,&x2);
    379:   VecDuplicate(x2,&r2);

    381:   /*
    382:   Create the block vectors
    383:   */
    384:   bx[0] = x1;
    385:   bx[1] = x2;
    386:   VecCreateNest(PETSC_COMM_WORLD,2,NULL,bx,&x);
    387:   VecAssemblyBegin(x);
    388:   VecAssemblyEnd(x);
    389:   VecDestroy(&x1);
    390:   VecDestroy(&x2);

    392:   bx[0] = r1;
    393:   bx[1] = r2;
    394:   VecCreateNest(PETSC_COMM_WORLD,2,NULL,bx,&r);
    395:   VecDestroy(&r1);
    396:   VecDestroy(&r2);
    397:   VecAssemblyBegin(r);
    398:   VecAssemblyEnd(r);

    400:   /*
    401:   Create sub Jacobian matrix data structure
    402:   */
    403:   MatCreate(PETSC_COMM_WORLD, &j11);
    404:   MatSetSizes(j11, 1, 1, 1, 1);
    405:   MatSetType(j11, MATSEQAIJ);
    406:   MatSetUp(j11);

    408:   MatCreate(PETSC_COMM_WORLD, &j12);
    409:   MatSetSizes(j12, 1, 1, 1, 1);
    410:   MatSetType(j12, MATSEQAIJ);
    411:   MatSetUp(j12);

    413:   MatCreate(PETSC_COMM_WORLD, &j21);
    414:   MatSetSizes(j21, 1, 1, 1, 1);
    415:   MatSetType(j21, MATSEQAIJ);
    416:   MatSetUp(j21);

    418:   MatCreate(PETSC_COMM_WORLD, &j22);
    419:   MatSetSizes(j22, PETSC_DECIDE, PETSC_DECIDE, 1, 1);
    420:   MatSetType(j22, MATSEQAIJ);
    421:   MatSetUp(j22);
    422:   /*
    423:   Create block Jacobian matrix data structure
    424:   */
    425:   bA[0][0] = j11;
    426:   bA[0][1] = j12;
    427:   bA[1][0] = j21;
    428:   bA[1][1] = j22;

    430:   MatCreateNest(PETSC_COMM_WORLD,2,NULL,2,NULL,&bA[0][0],&J);
    431:   MatSetUp(J);
    432:   MatNestSetVecType(J,VECNEST);
    433:   MatDestroy(&j11);
    434:   MatDestroy(&j12);
    435:   MatDestroy(&j21);
    436:   MatDestroy(&j22);

    438:   MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
    439:   MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);

    441:   PetscOptionsHasName(NULL,"-hard",&flg);
    442:   if (!flg) {
    443:     /*
    444:     Set function evaluation routine and vector.
    445:     */
    446:     SNESSetFunction(snes,r,FormFunction1_block,NULL);

    448:     /*
    449:     Set Jacobian matrix data structure and Jacobian evaluation routine
    450:     */
    451:     SNESSetJacobian(snes,J,J,FormJacobian1_block,NULL);
    452:   } else {
    453:     SNESSetFunction(snes,r,FormFunction2_block,NULL);
    454:     SNESSetJacobian(snes,J,J,FormJacobian2_block,NULL);
    455:   }

    457:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    458:   Customize nonlinear solver; set runtime options
    459:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    461:   /*
    462:   Set linear solver defaults for this problem. By extracting the
    463:   KSP, KSP, and PC contexts from the SNES context, we can then
    464:   directly call any KSP, KSP, and PC routines to set various options.
    465:   */
    466:   SNESGetKSP(snes,&ksp);
    467:   KSPGetPC(ksp,&pc);
    468:   PCSetType(pc,PCNONE);
    469:   KSPSetTolerances(ksp,1.e-4,PETSC_DEFAULT,PETSC_DEFAULT,20);

    471:   /*
    472:   Set SNES/KSP/KSP/PC runtime options, e.g.,
    473:   -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
    474:   These options will override those specified above as long as
    475:   SNESSetFromOptions() is called _after_ any other customization
    476:   routines.
    477:   */
    478:   SNESSetFromOptions(snes);

    480:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    481:   Evaluate initial guess; then solve nonlinear system
    482:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    483:   if (!flg) {
    484:     VecSet(x,pfive);
    485:   } else {
    486:     Vec *vecs;
    487:     VecNestGetSubVecs(x, NULL, &vecs);
    488:     bv   = vecs[0];
    489: /*    VecBlockGetSubVec(x, 0, &bv); */
    490:     VecSetValue(bv, 0, 2.0, INSERT_VALUES);  /* xx[0] = 2.0; */
    491:     VecAssemblyBegin(bv);
    492:     VecAssemblyEnd(bv);

    494: /*    VecBlockGetSubVec(x, 1, &bv); */
    495:     bv   = vecs[1];
    496:     VecSetValue(bv, 0, 3.0, INSERT_VALUES);  /* xx[1] = 3.0; */
    497:     VecAssemblyBegin(bv);
    498:     VecAssemblyEnd(bv);
    499:   }
    500:   /*
    501:   Note: The user should initialize the vector, x, with the initial guess
    502:   for the nonlinear solver prior to calling SNESSolve().  In particular,
    503:   to employ an initial guess of zero, the user should explicitly set
    504:   this vector to zero by calling VecSet().
    505:   */
    506:   SNESSolve(snes,NULL,x);
    507:   SNESGetIterationNumber(snes,&its);
    508:   if (flg) {
    509:     Vec f;
    510:     VecView(x,PETSC_VIEWER_STDOUT_WORLD);
    511:     SNESGetFunction(snes,&f,0,0);
    512:     VecView(r,PETSC_VIEWER_STDOUT_WORLD);
    513:   }

    515:   PetscPrintf(PETSC_COMM_SELF,"number of SNES iterations = %D\n\n",its);

    517:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    518:   Free work space.  All PETSc objects should be destroyed when they
    519:   are no longer needed.
    520:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    521:   VecDestroy(&x); VecDestroy(&r);
    522:   MatDestroy(&J); SNESDestroy(&snes);
    523:   return(0);
    524: }

    526: /* ------------------------------------------------------------------- */
    529: static PetscErrorCode FormFunction1_block(SNES snes,Vec x,Vec f,void *dummy)
    530: {
    531:   Vec            *xx, *ff, x1,x2, f1,f2;
    532:   PetscScalar    ff_0, ff_1;
    533:   PetscScalar    xx_0, xx_1;
    534:   PetscInt       index,nb;

    538:   /* get blocks for function */
    539:   VecNestGetSubVecs(f, &nb, &ff);
    540:   f1   = ff[0];  f2 = ff[1];

    542:   /* get blocks for solution */
    543:   VecNestGetSubVecs(x, &nb, &xx);
    544:   x1   = xx[0];  x2 = xx[1];

    546:   /* get solution values */
    547:   index = 0;
    548:   VecGetValues(x1,1, &index, &xx_0);
    549:   VecGetValues(x2,1, &index, &xx_1);

    551:   /* Compute function */
    552:   ff_0 = xx_0*xx_0 + xx_0*xx_1 - 3.0;
    553:   ff_1 = xx_0*xx_1 + xx_1*xx_1 - 6.0;

    555:   /* set function values */
    556:   VecSetValue(f1, index, ff_0, INSERT_VALUES);

    558:   VecSetValue(f2, index, ff_1, INSERT_VALUES);

    560:   VecAssemblyBegin(f);
    561:   VecAssemblyEnd(f);
    562:   return(0);
    563: }

    565: /* ------------------------------------------------------------------- */
    568: static PetscErrorCode FormJacobian1_block(SNES snes,Vec x,Mat *jac,Mat *B,MatStructure *flag,void *dummy)
    569: {
    570:   Vec            *xx, x1,x2;
    571:   PetscScalar    xx_0, xx_1;
    572:   PetscInt       index,nb;
    573:   PetscScalar    A_00, A_01, A_10, A_11;
    574:   Mat            j11, j12, j21, j22;
    575:   Mat            **mats;

    579:   /* get blocks for solution */
    580:   VecNestGetSubVecs(x, &nb, &xx);
    581:   x1   = xx[0];  x2 = xx[1];

    583:   /* get solution values */
    584:   index = 0;
    585:   VecGetValues(x1,1, &index, &xx_0);
    586:   VecGetValues(x2,1, &index, &xx_1);

    588:   /* get block matrices */
    589:   MatNestGetSubMats(*jac,NULL,NULL,&mats);
    590:   j11  = mats[0][0];
    591:   j12  = mats[0][1];
    592:   j21  = mats[1][0];
    593:   j22  = mats[1][1];

    595:   /* compute jacobian entries */
    596:   A_00 = 2.0*xx_0 + xx_1;
    597:   A_01 = xx_0;
    598:   A_10 = xx_1;
    599:   A_11 = xx_0 + 2.0*xx_1;

    601:   /* set jacobian values */
    602:   MatSetValue(j11, 0,0, A_00, INSERT_VALUES);
    603:   MatSetValue(j12, 0,0, A_01, INSERT_VALUES);
    604:   MatSetValue(j21, 0,0, A_10, INSERT_VALUES);
    605:   MatSetValue(j22, 0,0, A_11, INSERT_VALUES);

    607:   *flag = SAME_NONZERO_PATTERN;

    609:   /* Assemble sub matrix */
    610:   MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
    611:   MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);
    612:   return(0);
    613: }

    615: /* ------------------------------------------------------------------- */
    618: static PetscErrorCode FormFunction2_block(SNES snes,Vec x,Vec f,void *dummy)
    619: {
    621:   PetscScalar    *xx,*ff;

    624:   /*
    625:   Get pointers to vector data.
    626:   - For default PETSc vectors, VecGetArray() returns a pointer to
    627:   the data array.  Otherwise, the routine is implementation dependent.
    628:   - You MUST call VecRestoreArray() when you no longer need access to
    629:   the array.
    630:   */
    631:   VecGetArray(x,&xx);
    632:   VecGetArray(f,&ff);

    634:   /*
    635:   Compute function
    636:   */
    637:   ff[0] = PetscSinScalar(3.0*xx[0]) + xx[0];
    638:   ff[1] = xx[1];

    640:   /*
    641:   Restore vectors
    642:   */
    643:   VecRestoreArray(x,&xx);
    644:   VecRestoreArray(f,&ff);
    645:   return(0);
    646: }

    648: /* ------------------------------------------------------------------- */
    651: static PetscErrorCode FormJacobian2_block(SNES snes,Vec x,Mat *jac,Mat *B,MatStructure *flag,void *dummy)
    652: {
    653:   PetscScalar    *xx,A[4];
    655:   PetscInt       idx[2] = {0,1};

    658:   /*
    659:   Get pointer to vector data
    660:   */
    661:   VecGetArray(x,&xx);

    663:   /*
    664:   Compute Jacobian entries and insert into matrix.
    665:   - Since this is such a small problem, we set all entries for
    666:   the matrix at once.
    667:   */
    668:   A[0]  = 3.0*PetscCosScalar(3.0*xx[0]) + 1.0; A[1] = 0.0;
    669:   A[2]  = 0.0;                     A[3] = 1.0;
    670:   MatSetValues(*jac,2,idx,2,idx,A,INSERT_VALUES);
    671:   *flag = SAME_NONZERO_PATTERN;

    673:   /*
    674:   Restore vector
    675:   */
    676:   VecRestoreArray(x,&xx);

    678:   /*
    679:   Assemble matrix
    680:   */
    681:   MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
    682:   MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);
    683:   return(0);
    684: }

    688: int main(int argc,char **argv)
    689: {
    690:   PetscMPIInt    size;

    693:   PetscInitialize(&argc,&argv,(char*)0,help);

    695:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
    696:   if (size != 1) SETERRQ(PETSC_COMM_WORLD, 1,"This is a uniprocessor example only!");

    698:   assembled_system();

    700:   block_system();

    702:   PetscFinalize();
    703:   return 0;
    704: }