static char help[] = "Newton's method for a two-variable system, sequential with variable switching\n\n";

/*T
   Concepts: SNES^basic example
T*/

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

/*
   User-defined application context - contains data needed by the
   application-provided call-back routines, FormFunction() and
   FormUpdate().
*/

typedef struct {
  Vec prevreg, region;
} AppCtx;

/*
   User-defined routines
*/
extern PetscErrorCode FormJacobian(SNES,Vec,Mat,Mat,void*);
extern PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
extern PetscErrorCode FormUpdate(SNES,PetscInt);
extern PetscErrorCode CalcInitialRegion(Vec,Vec);

#undef __FUNCT__
#define __FUNCT__ "main"
int main(int argc,char **argv)
{
  SNES           snes;         /* nonlinear solver context */
  KSP            ksp;          /* linear solver context */
  PC             pc;           /* preconditioner context */
  Vec            x,r;          /* solution, residual vectors */
  Mat            J;            /* Jacobian matrix */
  AppCtx         user;         /* user context with region, prevreg vectors */
  PetscErrorCode ierr;
  PetscInt       its;
  PetscMPIInt    size;
  PetscScalar    pfive = 50.;

  PetscInitialize(&argc,&argv,(char*)0,help);
  ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
  if (size > 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Example is only for sequential runs");

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     Create nonlinear solver context
     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  ierr = SNESCreate(PETSC_COMM_WORLD,&snes);CHKERRQ(ierr);

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     Create matrix and vector data structures; set corresponding routines
     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  /*
     Create vectors for solution and nonlinear function
  */
  ierr = VecCreate(PETSC_COMM_WORLD,&x);CHKERRQ(ierr);
  ierr = VecSetSizes(x,PETSC_DECIDE,2);CHKERRQ(ierr);
  ierr = VecSetFromOptions(x);CHKERRQ(ierr);
  ierr = VecDuplicate(x,&r);CHKERRQ(ierr);
  ierr = VecDuplicate(x,&(user.region));CHKERRQ(ierr);
  ierr = VecDuplicate(x,&(user.prevreg));CHKERRQ(ierr);

  /*
     Create Jacobian matrix data structure
  */
  ierr = MatCreate(PETSC_COMM_WORLD,&J);CHKERRQ(ierr);
  ierr = MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,2,2);CHKERRQ(ierr);
  ierr = MatSetFromOptions(J);CHKERRQ(ierr);
  ierr = MatSetUp(J);CHKERRQ(ierr);

  /*
   Set function evaluation routine and vector.
   */
  ierr = SNESSetFunction(snes,r,FormFunction,NULL);CHKERRQ(ierr);

  /*
   Set Jacobian matrix data structure and Jacobian evaluation routine
   */
  //  ierr = SNESSetJacobian(snes,J,J,FormJacobian,NULL);CHKERRQ(ierr);
   ierr = SNESSetJacobian(snes,J,J,SNESComputeJacobianDefault,NULL);CHKERRQ(ierr);

  /*
   Set update.
   */
  ierr = SNESSetUpdate(snes,FormUpdate);CHKERRQ(ierr);
  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     Customize nonlinear solver; set runtime options
   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  /*
     Set linear solver defaults for this problem. By extracting the
     KSP and PC contexts from the SNES context, we can then
     directly call any KSP and PC routines to set various options.
  */
  ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr);
  ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
  ierr = PCSetType(pc,PCNONE);CHKERRQ(ierr);
  ierr = KSPSetTolerances(ksp,1.e-4,PETSC_DEFAULT,PETSC_DEFAULT,20);CHKERRQ(ierr);

  /*
     Set SNES/KSP/KSP/PC runtime options, e.g.,
         -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
     These options will override those specified above as long as
     SNESSetFromOptions() is called _after_ any other customization
     routines.
  */
  ierr = SNESSetFromOptions(snes);CHKERRQ(ierr);

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     Evaluate initial guess; then solve nonlinear system
   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  ierr = VecSet(x,pfive);CHKERRQ(ierr);
  ierr = CalcInitialRegion(x,user.region);CHKERRQ(ierr);
  ierr = VecView(user.region,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
  ierr = SNESSetApplicationContext(snes,&user);CHKERRQ(ierr);

  /*
     Note: The user should initialize the vector, x, with the initial guess
     for the nonlinear solver prior to calling SNESSolve().  In particular,
     to employ an initial guess of zero, the user should explicitly set
     this vector to zero by calling VecSet().
  */

  ierr = SNESSolve(snes,NULL,x);CHKERRQ(ierr);
  ierr = SNESGetIterationNumber(snes,&its);CHKERRQ(ierr);
  ierr = VecView(x,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
  ierr = VecView(r,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);

  ierr = PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %D\n",its);CHKERRQ(ierr);

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     Free work space.  All PETSc objects should be destroyed when they
     are no longer needed.
   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

  ierr = VecDestroy(&x);CHKERRQ(ierr); ierr = VecDestroy(&r);CHKERRQ(ierr);
  ierr = MatDestroy(&J);CHKERRQ(ierr); ierr = SNESDestroy(&snes);CHKERRQ(ierr);
  ierr = PetscFinalize();
  return 0;
}
/* ------------------------------------------------------------------- */
#undef __FUNCT__
#define __FUNCT__ "FormFunction"
/*
   FormFunction - Evaluates nonlinear function, F(x).

   Input Parameters:
.  snes - the SNES context
.  x    - input vector
.  ctx  - optional user-defined context

   Output Parameter:
.  f - function vector
 */
PetscErrorCode FormFunction(SNES snes,Vec x,Vec f,void *ctx)
{
  PetscErrorCode    ierr;
  const PetscScalar *xx, *rr;
  PetscScalar       *ff;
  AppCtx            *user;

  /*
   Get pointers to vector data.
      - For default PETSc vectors, VecGetArray() returns a pointer to
        the data array.  Otherwise, the routine is implementation dependent.
      - You MUST call VecRestoreArray() when you no longer need access to
        the array.
   */
  ierr = PetscPrintf(PETSC_COMM_WORLD,"In Function\n");CHKERRQ(ierr);
  ierr = VecView(x,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
  ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr);
  ierr = SNESGetApplicationContext(snes,&user);CHKERRQ(ierr);
  ierr = VecView(user->region,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
  ierr = VecGetArrayRead(user->region,&rr);CHKERRQ(ierr);
  ierr = VecGetArray(f,&ff);CHKERRQ(ierr);

  /* Compute function */
  if (rr[0] == 1.0)
    ff[0] = xx[0]*xx[0]*xx[0] - 50.0;
  else
    ff[0] = (xx[0]-3.0)*(xx[0]-3.0)*(xx[0]-3.0) - 50.0;
  ff[1] = xx[1]*xx[1]*xx[1] - 50.0;

  /* Restore vectors */
  ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr);
  ierr = VecRestoreArrayRead(user->region,&rr);CHKERRQ(ierr);
  ierr = VecRestoreArray(f,&ff);CHKERRQ(ierr);
  return 0;
}

/* ------------------------------------------------------------------- */
#undef __FUNCT__
#define __FUNCT__ "FormUpdate"
/*
   FormUpdate - Use for switching variables in equation 1.

   Input Parameters:
.  snes - the SNES context
.  step    - step

   Output Parameter:
 */
PetscErrorCode FormUpdate(SNES snes,PetscInt step)
{
  PetscErrorCode    ierr;
  AppCtx            *user;
  Vec               x;
  PetscScalar       *xx,*rr,*pr;
  PetscInt          i;

  ierr = PetscPrintf(PETSC_COMM_WORLD,"=========================\n",step);CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_WORLD,"In Update, step = %D\n",step);CHKERRQ(ierr);
  ierr = SNESGetSolution(snes,&x);CHKERRQ(ierr);
  ierr = VecView(x,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
  ierr = SNESGetApplicationContext(snes,&user);CHKERRQ(ierr);
  ierr = VecView(user->region,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);

  ierr = VecGetArray(x,&xx);CHKERRQ(ierr);
  ierr = VecGetArray(user->region,&rr);CHKERRQ(ierr);
  ierr = VecGetArray(user->prevreg,&pr);CHKERRQ(ierr);

  for (i = 0; i < 1; i++)
  {
  /* Replace prevregion with region */
    pr[i] = rr[i];
    if (pr[i] == 1.0)
    {
      if (xx[i] > 5.0)
  /* Change from region 1 to 2 and update x */
      {
        rr[i] = 2.0;
        //        xx[i] += 5.0;
  /* Don't allow the solution vector to move too far across the region boundary */
        xx[i] = 8.0 + 1e-6;
      }
    }
    else
    {
      if (xx[i] <= 8.0)
  /* Change from region 2 to 1 and update x */
      {
        rr[i] = 1.0;
        //        xx[i] -= 5.0;
  /* Don't allow the solution vector to move too far across the region boundary */
        xx[i] = 5.0 - 1e-6;
      }
    }
  }

  /* Restore vectors */
  ierr = VecRestoreArray(x,&xx);CHKERRQ(ierr);
  ierr = VecRestoreArray(user->region,&rr);CHKERRQ(ierr);
  ierr = VecRestoreArray(user->prevreg,&pr);CHKERRQ(ierr);

  ierr = PetscPrintf(PETSC_COMM_WORLD,"Leaving Update, step = %D\n",step);CHKERRQ(ierr);
  ierr = VecView(x,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
  ierr = VecView(user->region,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_WORLD,"=========================\n",step);CHKERRQ(ierr);

  return 0;
}

/* ------------------------------------------------------------------- */
#undef __FUNCT__
#define __FUNCT__ "CalcInitialRegion"
/*
   CalcInitialRegion - Calculates which region the components of the X vector are in initially.

   Input Parameters:
.  x    - input vector

   Output Parameter:
.  reg - function vector
 */
PetscErrorCode CalcInitialRegion(Vec x,Vec reg)
{
  PetscErrorCode    ierr;
  PetscScalar       *xx;
  PetscScalar       *rr;
  PetscInt          i;

  /*
   Get pointers to vector data.
      - For default PETSc vectors, VecGetArray() returns a pointer to
        the data array.  Otherwise, the routine is implementation dependent.
      - You MUST call VecRestoreArray() when you no longer need access to
        the array.
   */
  ierr = VecGetArray(x,&xx);CHKERRQ(ierr);
  ierr = VecGetArray(reg,&rr);CHKERRQ(ierr);

  /* Compute function */
  for (i = 0; i < 1; i++)
    if (xx[i] <= 5.0)
      rr[i] = 1.0;
    else
    {
      rr[i] = 2.0;
      xx[i] += 3.0;
    }

  /* Restore vectors */
  ierr = VecRestoreArray(x,&xx);CHKERRQ(ierr);
  ierr = VecRestoreArray(reg,&rr);CHKERRQ(ierr);
  return 0;
}