C++ PetscOptionsGetReal函数代码示例

PetscErrorCode GetParams(AppCtx *user)
{
  PetscErrorCode ierr;
  PetscBool      flg;

  PetscFunctionBeginUser;
  /* Set default parameters */
  user->xmin       = 0.0; user->xmax = 128.0;
  user->Mv         = 1.0;
  user->L          = 1.0;
  user->kaeta      = 1.0;
  user->kav        = 0.5;
  user->Evf        = 9.09;
  user->A          = 9.09;
  user->B          = 9.09;
  user->cv0        = 1.13e-4;
  user->Sv         = 500.0;
  user->dt         = 1.0e-5;
  user->T          = 1.0e-2;
  user->graphics   = PETSC_TRUE;
  user->periodic   = PETSC_FALSE;
  user->lumpedmass = PETSC_FALSE;

  ierr = PetscOptionsGetReal(NULL,"-xmin",&user->xmin,&flg);CHKERRQ(ierr);
  ierr = PetscOptionsGetReal(NULL,"-xmax",&user->xmax,&flg);CHKERRQ(ierr);
  ierr = PetscOptionsGetReal(NULL,"-T",&user->T,&flg);CHKERRQ(ierr);
  ierr = PetscOptionsGetReal(NULL,"-dt",&user->dt,&flg);CHKERRQ(ierr);
  ierr = PetscOptionsBool("-graphics","Contour plot solutions at each timestep\n","None",user->graphics,&user->graphics,&flg);CHKERRQ(ierr);
  ierr = PetscOptionsBool("-periodic","Use periodic boundary conditions\n","None",user->periodic,&user->periodic,&flg);CHKERRQ(ierr);
  ierr = PetscOptionsBool("-lumpedmass","Use lumped mass matrix\n","None",user->lumpedmass,&user->lumpedmass,&flg);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

PetscErrorCode Parameter_settings(AppCtx *user)
{
  PetscErrorCode ierr;
  PetscBool      flg;

  PetscFunctionBeginUser;

  /* Set default parameters */
  user->ws     = 1.0;
  user->H      = 5.0;
  user->D      = 0.0;
  user->Pmax = user->Pmax_s  = 2.1;
  user->PM_min = 1.0;
  user->lambda = 0.1;
  user->q      = 1.0;
  user->mux    = PetscAsinScalar(user->PM_min/user->Pmax);
  user->sigmax = 0.1;
  user->sigmay = 0.1;
  user->rho    = 0.0;
  user->xmin   = -PETSC_PI;
  user->xmax   = PETSC_PI;
  user->bx     = DM_BOUNDARY_PERIODIC;
  user->by     = DM_BOUNDARY_GHOSTED;
  user->tf = user->tcl = -1;
  user->ymin   = -2.0;
  user->ymax   = 2.0;
  user->st_width = 1;

  ierr = PetscOptionsGetScalar(NULL,"-ws",&user->ws,&flg);CHKERRQ(ierr);
  ierr = PetscOptionsGetScalar(NULL,"-Inertia",&user->H,&flg);CHKERRQ(ierr);
  ierr = PetscOptionsGetScalar(NULL,"-D",&user->D,&flg);CHKERRQ(ierr);
  ierr = PetscOptionsGetScalar(NULL,"-Pmax",&user->Pmax,&flg);CHKERRQ(ierr);
  ierr = PetscOptionsGetScalar(NULL,"-PM_min",&user->PM_min,&flg);CHKERRQ(ierr);
  ierr = PetscOptionsGetScalar(NULL,"-lambda",&user->lambda,&flg);CHKERRQ(ierr);
  ierr = PetscOptionsGetScalar(NULL,"-q",&user->q,&flg);CHKERRQ(ierr);
  ierr = PetscOptionsGetScalar(NULL,"-mux",&user->mux,&flg);CHKERRQ(ierr);
  ierr = PetscOptionsGetScalar(NULL,"-muy",&user->muy,&flg);CHKERRQ(ierr);
  if (flg == 0) {
    user->muy = user->ws;
  }
  ierr = PetscOptionsGetScalar(NULL,"-xmin",&user->xmin,&flg);CHKERRQ(ierr);
  ierr = PetscOptionsGetScalar(NULL,"-xmax",&user->xmax,&flg);CHKERRQ(ierr);
  ierr = PetscOptionsGetScalar(NULL,"-ymin",&user->ymin,&flg);CHKERRQ(ierr);
  ierr = PetscOptionsGetScalar(NULL,"-ymax",&user->ymax,&flg);CHKERRQ(ierr);
  ierr = PetscOptionsGetInt(NULL,"-stencil_width",&user->st_width,&flg);CHKERRQ(ierr);
  ierr = PetscOptionsGetEnum("","-bx",BoundaryTypes,(PetscEnum*)&user->bx,&flg);CHKERRQ(ierr);
  ierr = PetscOptionsGetEnum("","-by",BoundaryTypes,(PetscEnum*)&user->by,&flg);CHKERRQ(ierr);
  ierr = PetscOptionsGetReal(NULL,"-tf",&user->tf,&flg);CHKERRQ(ierr);
  ierr = PetscOptionsGetReal(NULL,"-tcl",&user->tcl,&flg);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

PetscErrorCode FluidFieldCreate(MPI_Comm comm, FluidField *fluid)
{
  FluidField f;
  PetscErrorCode ierr;
  
  PetscFunctionBegin;
  ierr = PetscNew(struct _FluidField, &f); CHKERRQ(ierr);

  // Default values

  f->comm = comm;
  f->mu = 1e-3; // pN s / um^2
  ierr = PetscOptionsGetReal(0,"-fluid_viscosity",&f->mu,0); CHKERRQ(ierr);

  int nmax = 3;
  f->lens.x = 7;
  f->lens.y = 7;
  f->lens.z = 0;
  ierr = PetscOptionsGetRealArray(0,"-fluid_lens", &f->lens.x, &nmax, 0); CHKERRQ(ierr);
  f->is3D = (nmax == 3);

  PetscReal dx = 1;
  ierr = PetscOptionsGetReal(0,"-fluid_dx",&dx,0); CHKERRQ(ierr);
  f->dh.x = dx;
  f->dh.y = dx;
  f->dh.z = dx;
  
  nmax = 3;
  f->dims.x = f->lens.x / f->dh.x;
  f->dims.y = f->lens.y / f->dh.y;
  f->dims.z = f->lens.z / f->dh.z;
  ierr = PetscOptionsGetIntArray(0,"-fluid_dims", &f->dims.x, &nmax, 0); CHKERRQ(ierr);

  if( !f->is3D ) {
    f->lens.z = 0;
    f->dims.z = 0;
    f->dh.z = 0;
  }

  // Create BC index set
  ierr = ArrayCreate("dirichletBC",sizeof(MatStencil),&f->dirichletBC); CHKERRQ(ierr);

  ierr = FluidFieldRegisterEvents(); CHKERRQ(ierr);

  *fluid = f;

  PetscFunctionReturn(0);
}

// ---- Setter ----
PetscErrorCode BPSSetFromOptions(BPS self) {

  PetscBool find;
  PetscReal zmax = 20.0;
  PetscInt num_zs = 3;
  char type[10] = "line";
  PetscErrorCode ierr;
  
  //  ierr = PetscOptionsGetInt(NULL, "-bps_num_zs", &num_zs, &find); CHKERRQ(ierr);
  //  ierr = PetscOptionsGetReal(NULL, "-bps_zmax", &zmax, &find); CHKERRQ(ierr);
  //  ierr = PetscOptionsGetString(NULL, "-bps_type", type, 10, &find); CHKERRQ(ierr);
  ierr = PetscOptionsGetInt(NULL, NULL, "-bps_num_zs", &num_zs, &find); CHKERRQ(ierr);
  ierr = PetscOptionsGetReal(NULL, NULL, "-bps_zmax", &zmax, &find); CHKERRQ(ierr);
  ierr = PetscOptionsGetString(NULL, NULL, "-bps_type", type, 10, &find); CHKERRQ(ierr);

  if(strcmp(type, "line") == 0 ) {
    ierr = BPSSetLine(self, zmax, num_zs); CHKERRQ(ierr); 
  } else if(strcmp(type, "exp") == 0) {
    ierr = BPSSetExp(self, zmax, num_zs, 5.0); CHKERRQ(ierr);
  } else
    SETERRQ(self->comm, 1, "-bps_type must be line or exp");

  return 0;

}

PetscErrorCode MatLUFactorNumeric_Essl(Mat F,Mat A,const MatFactorInfo *info)
{
  Mat_SeqAIJ     *aa  =(Mat_SeqAIJ*)(A)->data;
  Mat_Essl       *essl=(Mat_Essl*)(F)->spptr;
  PetscErrorCode ierr;
  int            nessl,i,one = 1;

  PetscFunctionBegin;
  ierr = PetscBLASIntCast(A->rmap->n,&nessl);CHKERRQ(ierr);
  /* copy matrix data into silly ESSL data structure (1-based Frotran style) */
  for (i=0; i<A->rmap->n+1; i++) essl->ia[i] = aa->i[i] + 1;
  for (i=0; i<aa->nz; i++) essl->ja[i] = aa->j[i] + 1;

  ierr = PetscMemcpy(essl->a,aa->a,(aa->nz)*sizeof(PetscScalar));CHKERRQ(ierr);

  /* set Essl options */
  essl->iparm[0] = 1;
  essl->iparm[1] = 5;
  essl->iparm[2] = 1;
  essl->iparm[3] = 0;
  essl->rparm[0] = 1.e-12;
  essl->rparm[1] = 1.0;

  ierr = PetscOptionsGetReal(((PetscObject)A)->prefix,"-matessl_lu_threshold",&essl->rparm[1],NULL);CHKERRQ(ierr);

  dgsf(&one,&nessl,&essl->nz,essl->a,essl->ia,essl->ja,&essl->lna,essl->iparm,essl->rparm,essl->oparm,essl->aux,&essl->naux);

  F->ops->solve     = MatSolve_Essl;
  (F)->assembled    = PETSC_TRUE;
  (F)->preallocated = PETSC_TRUE;
  PetscFunctionReturn(0);
}

/*@C
   PetscBagRegisterReal - add a real value to the bag

   Logically Collective on PetscBag

   Input Parameter:
+  bag - the bag of values
.  addr - location of double in struct
.  mdefault - the initial value
.  name - name of the variable
-  help - longer string with more information about the value

   Level: beginner

.seealso: PetscBag, PetscBagSetName(), PetscBagView(), PetscBagLoad(), PetscBagGetData()
           PetscBagRegisterInt(), PetscBagRegisterBool(), PetscBagRegisterScalar()
           PetscBagSetFromOptions(), PetscBagCreate(), PetscBagGetName(), PetscBagRegisterEnum()

@*/
PetscErrorCode PetscBagRegisterReal(PetscBag bag,void *addr,PetscReal mdefault, const char *name, const char *help)
{
  PetscErrorCode ierr;
  PetscBagItem   item;
  char           nname[PETSC_BAG_NAME_LENGTH+1];
  PetscBool      printhelp;

  PetscFunctionBegin;
  nname[0] = '-';
  nname[1] = 0;
  ierr     = PetscStrncat(nname,name,PETSC_BAG_NAME_LENGTH-1);CHKERRQ(ierr);
  ierr     = PetscOptionsHasName(NULL,"-help",&printhelp);CHKERRQ(ierr);
  if (printhelp) {
    ierr = (*PetscHelpPrintf)(bag->bagcomm,"  -%s%s <%g>: %s \n",bag->bagprefix ? bag->bagprefix : "",name,(double)mdefault,help);CHKERRQ(ierr);
  }
  ierr = PetscOptionsGetReal(bag->bagprefix,nname,&mdefault,NULL);CHKERRQ(ierr);

  ierr         = PetscNew(&item);CHKERRQ(ierr);
  item->dtype  = PETSC_REAL;
  item->offset = ((char*)addr) - ((char*)bag);
  if (item->offset > bag->bagsize) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Registered item %s %s is not in bag memory space",name,help);
  item->next        = 0;
  item->msize       = 1;
  *(PetscReal*)addr = mdefault;
  ierr              = PetscBagRegister_Private(bag,item,name,help);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

int main(int argc, char **argv){

	PetscInitialize(&argc, &argv, PETSC_NULL, PETSC_NULL);

	int Mx, My, ix, iy, LowerPML;
	double hx, hy, EpsUniformr, EpsUniformi, x, y, R;

	// todo: read M, i, R, LowerPML, EpsUniform
 	PetscOptionsGetInt(PETSC_NULL,"-LowerPML", &LowerPML,PETSC_NULL);
 	PetscOptionsGetInt(PETSC_NULL,"-Mx", &Mx,PETSC_NULL);
 	PetscOptionsGetInt(PETSC_NULL,"-My", &My,PETSC_NULL);

 	PetscOptionsGetReal(PETSC_NULL,"-hx", &hx,PETSC_NULL);
 	PetscOptionsGetReal(PETSC_NULL,"-hy", &hy,PETSC_NULL);
 	PetscOptionsGetReal(PETSC_NULL,"-epsuniformr", &EpsUniformr,PETSC_NULL);
 	PetscOptionsGetReal(PETSC_NULL,"-epsuniformi", &EpsUniformi,PETSC_NULL);
 	PetscOptionsGetReal(PETSC_NULL,"-R", &R,PETSC_NULL);

	FILE *ptfEpsr, *ptfEpsi, *ptfProfile;
	ptfEpsr = fopen("geometry.epsr","w"); 
	ptfEpsi = fopen("geometry.epsi","w"); 
	ptfProfile = fopen("geometry.prof","w");

	// column major, so iy on the inside
	for (ix = 0; ix < Mx; ix++) for (iy = 0; iy < My; iy++) {

		x = ComputeGridpoint(ix, Mx, hx, LowerPML); 
		y = ComputeGridpoint(iy, My, hy, LowerPML);
		
		if (x*x + y*y < R*R){
			fprintf(ptfEpsr,"%g\n", EpsUniformr-1);
			fprintf(ptfEpsi,"%g\n", EpsUniformi);
			fprintf(ptfProfile,"%g\n", 1.0);
		}else{
			fprintf(ptfEpsr,"%g\n", 0.0);
			fprintf(ptfEpsi,"%g\n", 0.0);
			fprintf(ptfProfile,"%g\n", 0.0);
		}
	}


	fclose(ptfEpsr);
	fclose(ptfEpsi);
	fclose(ptfProfile);
	PetscFinalize();
	return 0;
}

/*
  AppCtxInitialize - Sets initial values for the application context parameters

  Input:
    ptr - void user-defined application context

  Output:
    ptr - user-defined application context with the default or user-provided
             parameters
*/
static int AppCtxInitialize(void *ptr) {

  AppCtx *user = (AppCtx*)ptr;
  PetscTruth flg;
  int info;

  /* Specify default parameters */
  user->mx = user->my = 11;
  user->b = -0.5;
  user->t = 0.5;
  user->l = -0.5;
  user->r = 0.5;
  user->fx=0.5;
  user->fy=0.5;
  user->bheight=0.0;

  /* Check for command line arguments that override defaults */
  info = PetscOptionsGetInt(TAO_NULL, "-mx", &user->mx, &flg); CHKERRQ(info);
  info = PetscOptionsGetInt(TAO_NULL, "-my", &user->my, &flg); CHKERRQ(info);
  info = PetscOptionsGetReal(TAO_NULL, "-bottom", &user->b, &flg); CHKERRQ(info);
  info = PetscOptionsGetReal(TAO_NULL, "-top", &user->t, &flg); CHKERRQ(info);
  info = PetscOptionsGetReal(TAO_NULL, "-left", &user->l, &flg); CHKERRQ(info);
  info = PetscOptionsGetReal(TAO_NULL, "-right", &user->r, &flg); CHKERRQ(info);
  info = PetscOptionsGetReal(PETSC_NULL,"-bmx",&user->fx,&flg); CHKERRQ(info);
  info = PetscOptionsGetReal(PETSC_NULL,"-bmy",&user->fy,&flg); CHKERRQ(info);
  info = PetscOptionsGetReal(PETSC_NULL,"-bheight",&user->bheight,&flg); CHKERRQ(info);

  user->hx = (user->r - user->l) / (user->mx - 1);
  user->hy = (user->t - user->b) / (user->my - 1);
  user->area = 0.5 * user->hx * user->hy;
  info = PetscLogFlops(8); CHKERRQ(info);

  return 0;

} /* AppCtxInitialize */

PetscErrorCode PotSetFromOptions(Pot self) {

  MPI_Comm comm; PetscObjectGetComm((PetscObject)self, &comm);
  
  char type[10];
  PetscErrorCode ierr;
  PetscBool find;

  ierr = PetscOptionsGetString(NULL, NULL, "-pot_type", type, 10, &find); CHKERRQ(ierr);

  if(!find) 
    SETERRQ(comm, 1, "options -pot_type is not found");
  
  if(strcmp(type, "harm") == 0) {
    PetscReal a;
    ierr = PetscOptionsGetReal(NULL, NULL, "-pot_a", &a, &find); CHKERRQ(ierr);
    if(!find)
      SETERRQ(comm, 1, "-pot_a is not found");
    PotSetHarm(self, a);

  } else if(strcmp(type, "slater") == 0) {
    PetscReal v0, z;
    PetscInt n;
    ierr = PetscOptionsGetReal(NULL, NULL, "-pot_v0", &v0, &find); CHKERRQ(ierr);
    if(!find)
      SETERRQ(comm, 1, "-pot_v0 is not found");
    ierr = PetscOptionsGetReal(NULL, NULL, "-pot_z", &z, &find); CHKERRQ(ierr);
    if(!find)
      SETERRQ(comm, 1, "-pot_z is not found");
    ierr = PetscOptionsGetInt(NULL, NULL, "-pot_n", &n, &find); CHKERRQ(ierr);
    if(!find)
      SETERRQ(comm, 1, "-pot_n is not found");
    PotSetSlater(self, v0, n, z);

  } else if(strcmp(type, "morse") == 0) {
    PetscReal D0, a, Re;
    ierr = PetscOptionsGetReal(NULL, NULL, "-pot_D0", &D0, &find); CHKERRQ(ierr);
    if(!find)
      SETERRQ(comm, 1, "-pot_D0 is not found");
    ierr = PetscOptionsGetReal(NULL, NULL, "-pot_a", &a, &find); CHKERRQ(ierr);
    if(!find)
      SETERRQ(comm, 1, "-pot_z is not found");
    ierr = PetscOptionsGetReal(NULL, NULL, "-pot_Re", &Re, &find); CHKERRQ(ierr);
    if(!find)
      SETERRQ(comm, 1, "-pot_Re is not found");

    PotSetMorse(self, D0, a, Re);

  } else if(strcmp(type, "coulomb") == 0) {
    PetscReal z = -1.0;
    ierr = PetscOptionsGetReal(NULL, NULL, "-pot_z", &z, &find); CHKERRQ(ierr);
    ierr = PotSetCoulombNE(self, 0, 0.0, z); CHKERRQ(ierr);
  } else {
    char msg[100]; sprintf(msg, "unsupported pot_type: %s", type);
    SETERRQ(comm, 1, msg);
  }
  return 0;
}

int main(int argc,char **args)
{
  Mat            C;
  Vec            u,x,b,e;
  PetscInt       i,n = 10,midx[3];
  PetscErrorCode ierr;
  PetscScalar    v[3];
  PetscReal      omega = 1.0,norm;

  PetscInitialize(&argc,&args,(char*)0,help);
  ierr = PetscOptionsGetReal(NULL,"-omega",&omega,NULL);CHKERRQ(ierr);
  ierr = PetscOptionsGetInt(NULL,"-n",&n,NULL);CHKERRQ(ierr);

  ierr = MatCreate(PETSC_COMM_SELF,&C);CHKERRQ(ierr);
  ierr = MatSetSizes(C,n,n,n,n);CHKERRQ(ierr);
  ierr = MatSetType(C,MATSEQDENSE);CHKERRQ(ierr);
  ierr = MatSetUp(C);CHKERRQ(ierr);
  ierr = VecCreateSeq(PETSC_COMM_SELF,n,&b);CHKERRQ(ierr);
  ierr = VecCreateSeq(PETSC_COMM_SELF,n,&x);CHKERRQ(ierr);
  ierr = VecCreateSeq(PETSC_COMM_SELF,n,&u);CHKERRQ(ierr);
  ierr = VecCreateSeq(PETSC_COMM_SELF,n,&e);CHKERRQ(ierr);
  ierr = VecSet(u,1.0);CHKERRQ(ierr);
  ierr = VecSet(x,0.0);CHKERRQ(ierr);

  v[0] = -1.; v[1] = 2.; v[2] = -1.;
  for (i=1; i<n-1; i++) {
    midx[0] = i-1; midx[1] = i; midx[2] = i+1;
    ierr    = MatSetValues(C,1,&i,3,midx,v,INSERT_VALUES);CHKERRQ(ierr);
  }
  i    = 0; midx[0] = 0; midx[1] = 1;
  v[0] = 2.0; v[1] = -1.;
  ierr = MatSetValues(C,1,&i,2,midx,v,INSERT_VALUES);CHKERRQ(ierr);
  i    = n-1; midx[0] = n-2; midx[1] = n-1;
  v[0] = -1.0; v[1] = 2.;
  ierr = MatSetValues(C,1,&i,2,midx,v,INSERT_VALUES);CHKERRQ(ierr);

  ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);

  ierr = MatMult(C,u,b);CHKERRQ(ierr);

  for (i=0; i<n; i++) {
    ierr = MatSOR(C,b,omega,SOR_FORWARD_SWEEP,0.0,1,1,x);CHKERRQ(ierr);
    ierr = VecWAXPY(e,-1.0,x,u);CHKERRQ(ierr);
    ierr = VecNorm(e,NORM_2,&norm);CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_SELF,"2-norm of error %g\n",(double)norm);CHKERRQ(ierr);
  }
  ierr = MatDestroy(&C);CHKERRQ(ierr);
  ierr = VecDestroy(&x);CHKERRQ(ierr);
  ierr = VecDestroy(&b);CHKERRQ(ierr);
  ierr = VecDestroy(&u);CHKERRQ(ierr);
  ierr = VecDestroy(&e);CHKERRQ(ierr);
  ierr = PetscFinalize();
  return 0;
}

/*
   MSA_InitialPoint - Calculates the initial guess in one of three ways. 

   Input Parameters:
.  user - user-defined application context
.  X - vector for initial guess

   Output Parameters:
.  X - newly computed initial guess
*/
static int MSA_InitialPoint(AppCtx * user, Vec X)
{
  int      info;
  PetscInt   start2=-1,i,j;
  PetscReal   start1=0;
  PetscTruth flg1,flg2;

  info = PetscOptionsGetReal(PETSC_NULL,"-start",&start1,&flg1); CHKERRQ(info);
  info = PetscOptionsGetInt(PETSC_NULL,"-random",&start2,&flg2); CHKERRQ(info);

  if (flg1){ /* The zero vector is reasonable */
 
    info = VecSet(X, start1); CHKERRQ(info);

  } else if (flg2 && start2>0){ /* Try a random start between -0.5 and 0.5 */

    PetscRandom rctx;  PetscScalar np5=-0.5;

    info = PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
    CHKERRQ(info);
    for (i=0; i<start2; i++){
      info = VecSetRandom(X, rctx); CHKERRQ(info);
    }
    info = PetscRandomDestroy(rctx); CHKERRQ(info);
    info = VecShift(X, np5); CHKERRQ(info);

  } else { /* Take an average of the boundary conditions */

    PetscInt xs,xm,ys,ym;
    PetscInt mx=user->mx,my=user->my;
    PetscScalar **x;
    
    /* Get local mesh boundaries */
    info = DAGetCorners(user->da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL); CHKERRQ(info);
    
    /* Get pointers to vector data */
    info = DAVecGetArray(user->da,X,(void**)&x);

    /* Perform local computations */    
    for (j=ys; j<ys+ym; j++){
      for (i=xs; i< xs+xm; i++){
	x[j][i] = ( ((j+1)*user->bottom[i-xs+1]+(my-j+1)*user->top[i-xs+1])/(my+2)+
		   ((i+1)*user->left[j-ys+1]+(mx-i+1)*user->right[j-ys+1])/(mx+2))/2.0; 
      }
    }
    
    /* Restore vectors */
    info = DAVecRestoreArray(user->da,X,(void**)&x);  CHKERRQ(info);

    info = PetscLogFlops(9*xm*ym); CHKERRQ(info);
    
  }
  return 0;
}

int main(int argc,char **argv)
{
  PetscErrorCode     ierr;                  /* used to check for functions returning nonzeros */
  PetscReal          zero=0.0;
  Vec                x;                     /* solution vector */
  Mat                H;
  Tao                tao;                   /* Tao solver context */
  PetscBool          flg;
  PetscMPIInt        size,rank;                  /* number of processes running */
  AppCtx             user;                  /* user-defined application context */

  /* Initialize TAO and PETSc */
  PetscInitialize(&argc,&argv,(char*)0,help);
  ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
  ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
  if (size >1) SETERRQ(PETSC_COMM_SELF,1,"Incorrect number of processors");

  /* Initialize problem parameters */
  user.n = 2; user.alpha = 99.0;
  /* Check for command line arguments to override defaults */
  ierr = PetscOptionsGetInt(NULL,NULL,"-n",&user.n,&flg);CHKERRQ(ierr);
  ierr = PetscOptionsGetReal(NULL,NULL,"-alpha",&user.alpha,&flg);CHKERRQ(ierr);

  /* Allocate vectors for the solution and gradient */
  ierr = VecCreateSeq(PETSC_COMM_SELF,user.n,&x);CHKERRQ(ierr);
  ierr = MatCreateSeqBAIJ(PETSC_COMM_SELF,2,user.n,user.n,1,NULL,&H);CHKERRQ(ierr);

  /* The TAO code begins here */

  /* Create TAO solver with desired solution method */
  ierr = TaoCreate(PETSC_COMM_SELF,&tao);CHKERRQ(ierr);
  ierr = TaoSetType(tao,TAOLMVM);CHKERRQ(ierr);

  /* Set solution vec and an initial guess */
  ierr = VecSet(x, zero);CHKERRQ(ierr);
  ierr = TaoSetInitialVector(tao,x);CHKERRQ(ierr);

  /* Set routines for function, gradient, hessian evaluation */
  ierr = TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,&user);CHKERRQ(ierr);
  ierr = TaoSetHessianRoutine(tao,H,H,FormHessian,&user);CHKERRQ(ierr);

  /* Check for TAO command line options */
  ierr = TaoSetFromOptions(tao);CHKERRQ(ierr);

  /* SOLVE THE APPLICATION */
  ierr = TaoSolve(tao);CHKERRQ(ierr);

  ierr = TaoDestroy(&tao);CHKERRQ(ierr);
  ierr = VecDestroy(&x);CHKERRQ(ierr);
  ierr = MatDestroy(&H);CHKERRQ(ierr);

  PetscFinalize();
  return 0;
}

PetscErrorCode  PetscDrawUtilitySetCmap(const char colormap[],int mapsize,unsigned char R[],unsigned char G[],unsigned char B[])
{
  int             i,j;
  const char      *cmap_name_list[sizeof(PetscDrawCmapTable)/sizeof(PetscDrawCmapTable[0])];
  PetscInt        id = 0, count = (PetscInt)(sizeof(cmap_name_list)/sizeof(char*));
  PetscBool       reverse = PETSC_FALSE, brighten = PETSC_FALSE;
  PetscReal       beta = 0;
  PetscErrorCode  ierr;

  PetscFunctionBegin;
  for (i=0; i<count; i++) cmap_name_list[i] = PetscDrawCmapTable[i].name;
  if (colormap && colormap[0]) {
    PetscBool match = PETSC_FALSE;
    for (id=0; !match && id<count; id++) {ierr = PetscStrcasecmp(colormap,cmap_name_list[id],&match);CHKERRQ(ierr);}
    if (!match) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Colormap '%s' not found",colormap);
  }
  ierr = PetscOptionsGetEList(NULL,NULL,"-draw_cmap",cmap_name_list,count,&id,NULL);CHKERRQ(ierr);
  ierr = PetscOptionsGetBool(NULL,NULL,"-draw_cmap_reverse",&reverse,NULL);CHKERRQ(ierr);
  ierr = PetscOptionsGetReal(NULL,NULL,"-draw_cmap_brighten",&beta,&brighten);CHKERRQ(ierr);
  if (brighten && (beta <= (PetscReal)-1 || beta >= (PetscReal)+1)) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"brighten parameter %g must be in the range (-1,1)",(double)beta);

  if (PetscDrawCmapTable[id].cmap) {
    ierr = PetscDrawCmapTable[id].cmap(mapsize,R,G,B);CHKERRQ(ierr);
  } else {
    const unsigned char (*rgb)[3] = PetscDrawCmapTable[id].data;
    if (mapsize != 256-PETSC_DRAW_BASIC_COLORS) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Colormap '%s' with size %d not supported",cmap_name_list[id],mapsize);
    for (i=0; i<mapsize; i++) {R[i] = rgb[i][0]; G[i] = rgb[i][1]; B[i] = rgb[i][2];}
  }

  if (reverse) {
    i = 0; j = mapsize-1;
    while(i < j) {
#define SWAP(a,i,j) do { unsigned char t = a[i]; a[i] = a[j]; a[j] = t; } while (0)
      SWAP(R,i,j);
      SWAP(G,i,j);
      SWAP(B,i,j);
#undef SWAP
      i++; j--;
    }
  }

  if (brighten) {
    PetscReal gamma = (beta > 0.0) ? (1 - beta) : (1 / (1 + beta));
    for (i=0; i<mapsize; i++) {
      PetscReal r = PetscPowReal((PetscReal)R[i]/255,gamma);
      PetscReal g = PetscPowReal((PetscReal)G[i]/255,gamma);
      PetscReal b = PetscPowReal((PetscReal)B[i]/255,gamma);
      R[i] = (unsigned char)(255*PetscMin(r,(PetscReal)1.0));
      G[i] = (unsigned char)(255*PetscMin(g,(PetscReal)1.0));
      B[i] = (unsigned char)(255*PetscMin(b,(PetscReal)1.0));
    }
  }
  PetscFunctionReturn(0);
}

PetscErrorCode GetParams(AppCtx *user)
{
  PetscErrorCode ierr;
  PetscBool      flg;

  PetscFunctionBeginUser;
  /* Set default parameters */
  user->xmin  = 0.0; user->xmax = 1.0;
  user->ymin  = 0.0; user->ymax = 1.0;
  user->Dv    = 1.0; user->Di=4.0;
  user->Evf   = 0.8; user->Eif = 1.2;
  user->A     = 1.0;
  user->kBT   = 0.11;
  user->kav   = 1.0; user->kai = 1.0; user->kaeta = 1.0;
  user->Rsurf = 10.0; user->Rbulk = 1.0;
  user->L     = 10.0; user->P_casc = 0.05;
  user->T     = 1.0e-2;    user->dt = 1.0e-4;
  user->VG    = 100.0;

  ierr = PetscOptionsGetReal(NULL,"-xmin",&user->xmin,&flg);CHKERRQ(ierr);
  ierr = PetscOptionsGetReal(NULL,"-xmax",&user->xmax,&flg);CHKERRQ(ierr);
  ierr = PetscOptionsGetReal(NULL,"-T",&user->T,&flg);CHKERRQ(ierr);
  ierr = PetscOptionsGetReal(NULL,"-dt",&user->dt,&flg);CHKERRQ(ierr);
  ierr = PetscOptionsGetReal(NULL,"-P_casc",&user->P_casc,&flg);CHKERRQ(ierr);
  ierr = PetscOptionsGetReal(NULL,"-VG",&user->VG,&flg);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

PetscErrorCode GetParams(AppCtx* user)
{
  PetscErrorCode ierr;
  PetscBool      flg;
  
  PetscFunctionBegin;

  /* Set default parameters */
  user->tsmonitor = PETSC_FALSE;
  user->xmin = 0.0; user->xmax = 1.0;
  user->ymin = 0.0; user->ymax = 1.0;
  user->T = 0.0002;    user->dt = 0.0001;
  user->gamma = 3.2E-4; user->theta_c = 0;

  ierr = PetscOptionsGetBool(PETSC_NULL,"-ts_monitor",&user->tsmonitor,PETSC_NULL);CHKERRQ(ierr);
  ierr = PetscOptionsGetReal(PETSC_NULL,"-xmin",&user->xmin,&flg);CHKERRQ(ierr);
  ierr = PetscOptionsGetReal(PETSC_NULL,"-xmax",&user->xmax,&flg);CHKERRQ(ierr);
  ierr = PetscOptionsGetReal(PETSC_NULL,"-ymin",&user->ymin,&flg);CHKERRQ(ierr);
  ierr = PetscOptionsGetReal(PETSC_NULL,"-ymax",&user->ymax,&flg);CHKERRQ(ierr);
  ierr = PetscOptionsGetReal(PETSC_NULL,"-T",&user->T,&flg);CHKERRQ(ierr);
  ierr = PetscOptionsGetReal(PETSC_NULL,"-dt",&user->dt,&flg);CHKERRQ(ierr);
  ierr = PetscOptionsGetScalar(PETSC_NULL,"-gamma",&user->gamma,&flg);CHKERRQ(ierr);
  ierr = PetscOptionsGetScalar(PETSC_NULL,"-theta_c",&user->theta_c,&flg);CHKERRQ(ierr);

  PetscFunctionReturn(0);
}