2: static char help[] = "u`` + u^{2} = f. \n\
  3: Demonstrates the use of matrix-free Newton-Krylov methods in conjunction\n\
  4: with a user-provided preconditioner.  Simplified from ex6.c\n\
  5: Input arguments are:\n\
  6:    -snes_mf :      Use matrix-free Newton methods\n\
  7:    -user_precond : Employ a user-defined preconditioner.  Used only with\n\
  8:                    matrix-free methods in this example.\n\n";

 10: /*
 11:   Example: 
 12:   ./ex44 -snes_mf -snes_monitor -snes_view
 13:   ./ex44 -snes_mf -user_precond -snes_monitor -snes_view
 14:  */

 16: /*T
 17:    Concepts: SNES^different matrices for the Jacobian and preconditioner;
 18:    Concepts: SNES^matrix-free methods
 19:    Concepts: SNES^user-provided preconditioner;
 20:    Concepts: matrix-free methods
 21:    Concepts: user-provided preconditioner;
 22:    Processors: 1
 23: T*/

 25: #include <petscsnes.h>

 27: /* User-defined routines */
 28: PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
 29: PetscErrorCode MatrixFreePreconditioner(PC,Vec,Vec);

 31: int main(int argc,char **argv)
 32: {
 33:   SNES           snes;                /* SNES context */
 34:   KSP            ksp;                 /* KSP context */
 35:   PC             pc;                  /* PC context */
 36:   Vec            x,r,F;               /* vectors */
 38:   PetscInt       it,n = 5,i;
 39:   PetscMPIInt    size;
 40:   PetscReal      h,xp = 0.0;
 41:   PetscScalar    v,pfive = .5;
 42:   PetscBool      flg;

 44:   PetscInitialize(&argc,&argv,(char *)0,help);
 45:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 46:   if (size != 1) SETERRQ(PETSC_COMM_SELF,1,"This is a uniprocessor example only!");
 47:   PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);
 48:   h = 1.0/(n-1);

 50:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 51:      Create nonlinear solver context
 52:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 53:   SNESCreate(PETSC_COMM_WORLD,&snes);

 55:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 56:      Create vector data structures; set function evaluation routine
 57:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 58:   VecCreate(PETSC_COMM_SELF,&x);
 59:   VecSetSizes(x,PETSC_DECIDE,n);
 60:   VecSetFromOptions(x);
 61:   VecDuplicate(x,&r);
 62:   VecDuplicate(x,&F);

 64:   SNESSetFunction(snes,r,FormFunction,(void*)F);

 66:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 67:      Customize nonlinear solver; set runtime options
 68:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 69:   /* Set preconditioner for matrix-free method */
 70:   flg  = PETSC_FALSE;
 71:   PetscOptionsGetBool(PETSC_NULL,"-snes_mf",&flg,PETSC_NULL);
 72:   if (flg) {
 73:     SNESGetKSP(snes,&ksp);
 74:     KSPGetPC(ksp,&pc);
 75:     PetscOptionsHasName(PETSC_NULL,"-user_precond",&flg);
 76:     if (flg) { /* user-defined precond */
 77:       PCSetType(pc,PCSHELL);
 78:       PCShellSetApply(pc,MatrixFreePreconditioner);
 79:     } else {PCSetType(pc,PCNONE);}
 80:   }

 82:   SNESSetFromOptions(snes);

 84:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 85:      Initialize application:
 86:      Store right-hand-side of PDE and exact solution
 87:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 88:   //#if defined(TMP)
 89:   xp = 0.0;
 90:   for (i=0; i<n; i++) {
 91:     v = 6.0*xp + pow(xp+1.e-12,6.0); /* +1.e-12 is to prevent 0^6 */
 92:     VecSetValues(F,1,&i,&v,INSERT_VALUES);
 93:     xp += h;
 94:   }
 95:   //#endif
 96:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 97:      Evaluate initial guess; then solve nonlinear system
 98:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 99:   VecSet(x,pfive);
100:   SNESSolve(snes,PETSC_NULL,x);
101:   SNESGetIterationNumber(snes,&it);
102:   PetscPrintf(PETSC_COMM_SELF,"SNES iterations = %D\n\n",it);

104:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
105:      Free work space.  
106:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
107:   VecDestroy(&x);
108:   VecDestroy(&r);
109:   VecDestroy(&F);
110:   SNESDestroy(&snes);
111:   PetscFinalize();
112:   return 0;
113: }

115: /* ------------------------------------------------------------------- */
116: /* 
117:    FormInitialGuess - Forms initial approximation.

119:    Input Parameters:
120:    user - user-defined application context
121:    X - vector

123:    Output Parameter:
124:    X - vector
125:  */
126: PetscErrorCode FormFunction(SNES snes,Vec x,Vec f,void *dummy)
127: {
128:   PetscScalar    *xx,*ff,*FF,d;
130:   PetscInt       i,n;

132:   VecGetArray(x,&xx);
133:   VecGetArray(f,&ff);
134:   VecGetArray((Vec)dummy,&FF);
135:   VecGetSize(x,&n);
136:   d = (PetscReal)(n - 1); d = d*d;
137:   ff[0]   = xx[0];
138:   for (i=1; i<n-1; i++) {
139:     ff[i] = d*(xx[i-1] - 2.0*xx[i] + xx[i+1]) + xx[i]*xx[i] - FF[i];
140:   }
141:   ff[n-1] = xx[n-1] - 1.0;
142:   VecRestoreArray(x,&xx);
143:   VecRestoreArray(f,&ff);
144:   VecRestoreArray((Vec)dummy,&FF);
145:   return 0;
146: }

148: /* ------------------------------------------------------------------- */
149: /*
150:    MatrixFreePreconditioner - This routine demonstrates the use of a
151:    user-provided preconditioner.  This code implements just the null
152:    preconditioner, which of course is not recommended for general use.

154:    Input Parameters:
155: +  pc - preconditioner
156: -  x - input vector

158:    Output Parameter:
159: .  y - preconditioned vector
160: */
161: PetscErrorCode MatrixFreePreconditioner(PC pc,Vec x,Vec y)
162: {
164:   VecCopy(x,y);
165:   return 0;
166: }