Actual source code: ex1.c

      2: static char help[] = "Newton's method to solve a two-variable system, sequentially.\n\n";

      4: /*T
      5:    Concepts: SNES^basic uniprocessor example
      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:      petsc.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:    User-defined routines
     22: */

     30: int main(int argc,char **argv)
     31: {
     32:   SNES           snes;         /* nonlinear solver context */
     33:   KSP            ksp;         /* linear solver context */
     34:   PC             pc;           /* preconditioner context */
     35:   Vec            x,r;         /* solution, residual vectors */
     36:   Mat            J;            /* Jacobian matrix */
     38:   PetscInt       its;
     39:   PetscMPIInt    size;
     40:   PetscScalar    pfive = .5,*xx;
     41:   PetscTruth     flg;

     43:   PetscInitialize(&argc,&argv,(char *)0,help);
     44:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
     45:   if (size != 1) SETERRQ(1,"This is a uniprocessor example only!");

     47:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     48:      Create nonlinear solver context
     49:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

     51:   SNESCreate(PETSC_COMM_WORLD,&snes);

     53:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     54:      Create matrix and vector data structures; set corresponding routines
     55:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

     57:   /*
     58:      Create vectors for solution and nonlinear function
     59:   */
     60:   VecCreateSeq(PETSC_COMM_SELF,2,&x);
     61:   VecDuplicate(x,&r);

     63:   /*
     64:      Create Jacobian matrix data structure
     65:   */
     66:   MatCreate(PETSC_COMM_SELF,&J);
     67:   MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,2,2);
     68:   MatSetFromOptions(J);

     70:   PetscOptionsHasName(PETSC_NULL,"-hard",&flg);
     71:   if (!flg) {
     72:     /*
     73:      Set function evaluation routine and vector.
     74:     */
     75:     SNESSetFunction(snes,r,FormFunction1,PETSC_NULL);

     77:     /*
     78:      Set Jacobian matrix data structure and Jacobian evaluation routine
     79:     */
     80:     SNESSetJacobian(snes,J,J,FormJacobian1,PETSC_NULL);
     81:   } else {
     82:     SNESSetFunction(snes,r,FormFunction2,PETSC_NULL);
     83:     SNESSetJacobian(snes,J,J,FormJacobian2,PETSC_NULL);
     84:   }

     86:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     87:      Customize nonlinear solver; set runtime options
     88:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

     90:   /*
     91:      Set linear solver defaults for this problem. By extracting the
     92:      KSP, KSP, and PC contexts from the SNES context, we can then
     93:      directly call any KSP, KSP, and PC routines to set various options.
     94:   */
     95:   SNESGetKSP(snes,&ksp);
     96:   KSPGetPC(ksp,&pc);
     97:   PCSetType(pc,PCNONE);
     98:   KSPSetTolerances(ksp,1.e-4,PETSC_DEFAULT,PETSC_DEFAULT,20);

    100:   /*
    101:      Set SNES/KSP/KSP/PC runtime options, e.g.,
    102:          -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
    103:      These options will override those specified above as long as
    104:      SNESSetFromOptions() is called _after_ any other customization
    105:      routines.
    106:   */
    107:   SNESSetFromOptions(snes);

    109:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    110:      Evaluate initial guess; then solve nonlinear system
    111:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    112:   if (!flg) {
    113:     VecSet(x,pfive);
    114:   } else {
    115:     VecGetArray(x,&xx);
    116:     xx[0] = 2.0; xx[1] = 3.0;
    117:     VecRestoreArray(x,&xx);
    118:   }
    119:   /*
    120:      Note: The user should initialize the vector, x, with the initial guess
    121:      for the nonlinear solver prior to calling SNESSolve().  In particular,
    122:      to employ an initial guess of zero, the user should explicitly set
    123:      this vector to zero by calling VecSet().
    124:   */

    126:   SNESSolve(snes,PETSC_NULL,x);
    127:   SNESGetIterationNumber(snes,&its);
    128:   if (flg) {
    129:     Vec f;
    130:     VecView(x,PETSC_VIEWER_STDOUT_WORLD);
    131:     SNESGetFunction(snes,&f,0,0);
    132:     VecView(r,PETSC_VIEWER_STDOUT_WORLD);
    133:   }

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

    137:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    138:      Free work space.  All PETSc objects should be destroyed when they
    139:      are no longer needed.
    140:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    142:   VecDestroy(x); VecDestroy(r);
    143:   MatDestroy(J); SNESDestroy(snes);

    145:   PetscFinalize();
    146:   return 0;
    147: }
    148: /* ------------------------------------------------------------------- */
    151: /*
    152:    FormFunction1 - Evaluates nonlinear function, F(x).

    154:    Input Parameters:
    155: .  snes - the SNES context
    156: .  x - input vector
    157: .  dummy - optional user-defined context (not used here)

    159:    Output Parameter:
    160: .  f - function vector
    161:  */
    162: PetscErrorCode FormFunction1(SNES snes,Vec x,Vec f,void *dummy)
    163: {
    165:   PetscScalar    *xx,*ff;

    167:   /*
    168:      Get pointers to vector data.
    169:        - For default PETSc vectors, VecGetArray() returns a pointer to
    170:          the data array.  Otherwise, the routine is implementation dependent.
    171:        - You MUST call VecRestoreArray() when you no longer need access to
    172:          the array.
    173:   */
    174:   VecGetArray(x,&xx);
    175:   VecGetArray(f,&ff);

    177:   /*
    178:      Compute function
    179:   */
    180:   ff[0] = xx[0]*xx[0] + xx[0]*xx[1] - 3.0;
    181:   ff[1] = xx[0]*xx[1] + xx[1]*xx[1] - 6.0;

    183:   /*
    184:      Restore vectors
    185:   */
    186:   VecRestoreArray(x,&xx);
    187:   VecRestoreArray(f,&ff);

    189:   return 0;
    190: }
    191: /* ------------------------------------------------------------------- */
    194: /*
    195:    FormJacobian1 - Evaluates Jacobian matrix.

    197:    Input Parameters:
    198: .  snes - the SNES context
    199: .  x - input vector
    200: .  dummy - optional user-defined context (not used here)

    202:    Output Parameters:
    203: .  jac - Jacobian matrix
    204: .  B - optionally different preconditioning matrix
    205: .  flag - flag indicating matrix structure
    206: */
    207: PetscErrorCode FormJacobian1(SNES snes,Vec x,Mat *jac,Mat *B,MatStructure *flag,void *dummy)
    208: {
    209:   PetscScalar    *xx,A[4];
    211:   PetscInt       idx[2] = {0,1};

    213:   /*
    214:      Get pointer to vector data
    215:   */
    216:   VecGetArray(x,&xx);

    218:   /*
    219:      Compute Jacobian entries and insert into matrix.
    220:       - Since this is such a small problem, we set all entries for
    221:         the matrix at once.
    222:   */
    223:   A[0] = 2.0*xx[0] + xx[1]; A[1] = xx[0];
    224:   A[2] = xx[1]; A[3] = xx[0] + 2.0*xx[1];
    225:   MatSetValues(*jac,2,idx,2,idx,A,INSERT_VALUES);
    226:   *flag = SAME_NONZERO_PATTERN;

    228:   /*
    229:      Restore vector
    230:   */
    231:   VecRestoreArray(x,&xx);

    233:   /*
    234:      Assemble matrix
    235:   */
    236:   MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
    237:   MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);

    239:   return 0;
    240: }

    243: /* ------------------------------------------------------------------- */
    246: PetscErrorCode FormFunction2(SNES snes,Vec x,Vec f,void *dummy)
    247: {
    249:   PetscScalar    *xx,*ff;

    251:   /*
    252:      Get pointers to vector data.
    253:        - For default PETSc vectors, VecGetArray() returns a pointer to
    254:          the data array.  Otherwise, the routine is implementation dependent.
    255:        - You MUST call VecRestoreArray() when you no longer need access to
    256:          the array.
    257:   */
    258:   VecGetArray(x,&xx);
    259:   VecGetArray(f,&ff);

    261:   /*
    262:      Compute function
    263:   */
    264:   ff[0] = PetscSinScalar(3.0*xx[0]) + xx[0];
    265:   ff[1] = xx[1];

    267:   /*
    268:      Restore vectors
    269:   */
    270:   VecRestoreArray(x,&xx);
    271:   VecRestoreArray(f,&ff);

    273:   return 0;
    274: }
    275: /* ------------------------------------------------------------------- */
    278: PetscErrorCode FormJacobian2(SNES snes,Vec x,Mat *jac,Mat *B,MatStructure *flag,void *dummy)
    279: {
    280:   PetscScalar    *xx,A[4];
    282:   PetscInt       idx[2] = {0,1};

    284:   /*
    285:      Get pointer to vector data
    286:   */
    287:   VecGetArray(x,&xx);

    289:   /*
    290:      Compute Jacobian entries and insert into matrix.
    291:       - Since this is such a small problem, we set all entries for
    292:         the matrix at once.
    293:   */
    294:   A[0] = 3.0*PetscCosScalar(3.0*xx[0]) + 1.0; A[1] = 0.0;
    295:   A[2] = 0.0;                     A[3] = 1.0;
    296:   MatSetValues(*jac,2,idx,2,idx,A,INSERT_VALUES);
    297:   *flag = SAME_NONZERO_PATTERN;

    299:   /*
    300:      Restore vector
    301:   */
    302:   VecRestoreArray(x,&xx);

    304:   /*
    305:      Assemble matrix
    306:   */
    307:   MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
    308:   MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);

    310:   return 0;
    311: }