Commits

sarich committed 9ea0e1c

tao: update pde examples for separated state/design variables

Comments (1)

  1. Matt Knepley

    I am interested in doing a version of these which uses Plex. I think it would be much shorter and clearer what is going on. Are the equations and a sketch of the method written out anywhere?

Files changed (3)

src/tao/pde_constrained/examples/tutorials/elliptic.c

   /* Set up initial vectors and matrices */
   ierr = EllipticInitialize(&user); CHKERRQ(ierr);
 
-
-  //  ierr = Gather(user.x,user.y,user.state_scatter,user.u,user.design_scatter); CHKERRQ(ierr);
   ierr = VecDuplicate(user.y,&y0); CHKERRQ(ierr);
   ierr = VecDuplicate(user.y,&state); CHKERRQ(ierr);
   ierr = VecCopy(user.y,y0); CHKERRQ(ierr);
   ierr = VecScale(user->uwork, user->alpha); CHKERRQ(ierr);
   *f = 0.5 * (d1 + user->alpha*d2);
 
-  ierr = Gather(G,user->ywork,user->state_scatter,user->uwork,user->design_scatter); CHKERRQ(ierr);
   PetscFunctionReturn(0);
 
 }

src/tao/pde_constrained/examples/tutorials/hyperbolic.c

   PetscInt mx; /*  grid points in each direction */
   PetscInt nt; /*  Number of time steps */
   PetscInt ndata; /*  Number of data points per sample */
-  IS       s_is;
-  IS       d_is;
-  VecScatter state_scatter;
-  VecScatter design_scatter;
   VecScatter *uxi_scatter,*uyi_scatter,*ux_scatter,*uy_scatter,*ui_scatter;
   VecScatter *yi_scatter;
 
 } AppCtx;
 
 
-PetscErrorCode FormFunction(Tao, Vec, PetscReal*, void*);
-PetscErrorCode FormGradient(Tao, Vec, Vec, void*);
-PetscErrorCode FormFunctionGradient(Tao, Vec, PetscReal*, Vec, void*);
-PetscErrorCode FormJacobianState(Tao, Vec, Mat, Mat, Mat, void*);
-PetscErrorCode FormJacobianDesign(Tao, Vec, Mat,void*);
-PetscErrorCode FormConstraints(Tao, Vec, Vec, void*);
-PetscErrorCode FormHessian(Tao, Vec, Mat, Mat, void*);
-PetscErrorCode Gather(Vec x, Vec state, VecScatter s_scat, Vec design, VecScatter d_scat);
-PetscErrorCode Scatter(Vec x, Vec state, VecScatter s_scat, Vec design, VecScatter d_scat);
+PetscErrorCode FormFunction(Tao, Vec, Vec, PetscReal*, void*);
+PetscErrorCode FormGradient(Tao, Vec, Vec, Vec, Vec, void*);
+PetscErrorCode FormFunctionGradient(Tao, Vec, Vec, PetscReal*, Vec, Vec, void*);
+PetscErrorCode FormJacobian(Tao, Vec, Vec, Mat, Mat, Mat, Mat, void*);
+PetscErrorCode FormConstraints(Tao, Vec, Vec, Vec, void*);
+
 PetscErrorCode HyperbolicInitialize(AppCtx *user);
 PetscErrorCode HyperbolicDestroy(AppCtx *user);
 PetscErrorCode HyperbolicMonitor(Tao, void*);
 PetscErrorCode DesignMatMult(Mat,Vec,Vec);
 PetscErrorCode DesignMatMultTranspose(Mat,Vec,Vec);
 
-PetscErrorCode Scatter_yi(Vec,Vec*,VecScatter*,PetscInt); /*  y to y1,y2,...,y_nt */
+/*  y to y1,y2,...,y_nt */
+PetscErrorCode Scatter_yi(Vec,Vec*,VecScatter*,PetscInt);
 PetscErrorCode Gather_yi(Vec,Vec*,VecScatter*,PetscInt);
-PetscErrorCode Scatter_uxi_uyi(Vec,Vec*,VecScatter*,Vec*,VecScatter*,PetscInt); /*  u to ux_1,uy_1,ux_2,uy_2,...,u */
+
+/*  u to ux_1,uy_1,ux_2,uy_2,...,u */
+PetscErrorCode Scatter_uxi_uyi(Vec,Vec*,VecScatter*,Vec*,VecScatter*,PetscInt);
 PetscErrorCode Gather_uxi_uyi(Vec,Vec*,VecScatter*,Vec*,VecScatter*,PetscInt);
+PetscErrorCode Gather(Vec x, Vec subx, VecScatter scatter);
+PetscErrorCode Scatter(Vec x, Vec subx, VecScatter scatter);
+
 
 static  char help[]="";
 
 int main(int argc, char **argv)
 {
   PetscErrorCode     ierr;
-  Vec                x,x0;
+  Vec                state,y0;
+  Vec                design,u0;
   Tao                tao;
   TaoConvergedReason reason;
   AppCtx             user;
-  IS                 is_allstate,is_alldesign;
-  PetscInt           lo,hi,hi2,lo2,ksp_old;
+  PetscInt           ksp_old;
   PetscBool          flag;
   PetscInt           ntests = 1;
   PetscInt           i;
   user.ht = user.T/user.nt; /*  Time step */
   user.gamma = user.T*user.ht / (user.mx*user.mx);
 
-  ierr = VecCreate(PETSC_COMM_WORLD,&user.u);CHKERRQ(ierr);
   ierr = VecCreate(PETSC_COMM_WORLD,&user.y);CHKERRQ(ierr);
+  ierr = VecCreate(PETSC_COMM_WORLD,&user.u);CHKERRQ(ierr);
   ierr = VecCreate(PETSC_COMM_WORLD,&user.c);CHKERRQ(ierr);
-  ierr = VecSetSizes(user.u,PETSC_DECIDE,user.n-user.m);CHKERRQ(ierr);
   ierr = VecSetSizes(user.y,PETSC_DECIDE,user.m);CHKERRQ(ierr);
+  ierr = VecSetSizes(user.u,PETSC_DECIDE,user.n-user.m);CHKERRQ(ierr);
   ierr = VecSetSizes(user.c,PETSC_DECIDE,user.m);CHKERRQ(ierr);
-  ierr = VecSetFromOptions(user.u);CHKERRQ(ierr);
   ierr = VecSetFromOptions(user.y);CHKERRQ(ierr);
+  ierr = VecSetFromOptions(user.u);CHKERRQ(ierr);
   ierr = VecSetFromOptions(user.c);CHKERRQ(ierr);
 
-  /* Create scatters for reduced spaces.
-     If the state vector y and design vector u are partitioned as
-     [y_1; y_2; ...; y_np] and [u_1; u_2; ...; u_np] (with np = # of processors),
-     then the solution vector x is organized as
-     [y_1; u_1; y_2; u_2; ...; y_np; u_np].
-     The index sets user.s_is and user.d_is correspond to the indices of the
-     state and design variables owned by the current processor.
-  */
-  ierr = VecCreate(PETSC_COMM_WORLD,&x);CHKERRQ(ierr);
-
-  ierr = VecGetOwnershipRange(user.y,&lo,&hi);CHKERRQ(ierr);
-  ierr = VecGetOwnershipRange(user.u,&lo2,&hi2);CHKERRQ(ierr);
-
-  ierr = ISCreateStride(PETSC_COMM_SELF,hi-lo,lo,1,&is_allstate);CHKERRQ(ierr);
-  ierr = ISCreateStride(PETSC_COMM_SELF,hi-lo,lo+lo2,1,&user.s_is);CHKERRQ(ierr);
-
-  ierr = ISCreateStride(PETSC_COMM_SELF,hi2-lo2,lo2,1,&is_alldesign);CHKERRQ(ierr);
-  ierr = ISCreateStride(PETSC_COMM_SELF,hi2-lo2,hi+lo2,1,&user.d_is);CHKERRQ(ierr);
-
-  ierr = VecSetSizes(x,hi-lo+hi2-lo2,user.n);CHKERRQ(ierr);
-  ierr = VecSetFromOptions(x);CHKERRQ(ierr);
-
-  ierr = VecScatterCreate(x,user.s_is,user.y,is_allstate,&user.state_scatter);CHKERRQ(ierr);
-  ierr = VecScatterCreate(x,user.d_is,user.u,is_alldesign,&user.design_scatter);CHKERRQ(ierr);
-  ierr = ISDestroy(&is_alldesign);CHKERRQ(ierr);
-  ierr = ISDestroy(&is_allstate);CHKERRQ(ierr);
 
   /* Create TAO solver and set desired solution method */
   ierr = TaoCreate(PETSC_COMM_WORLD,&tao);CHKERRQ(ierr);
   /* Set up initial vectors and matrices */
   ierr = HyperbolicInitialize(&user);CHKERRQ(ierr);
 
-  ierr = Gather(x,user.y,user.state_scatter,user.u,user.design_scatter);CHKERRQ(ierr);
-  ierr = VecDuplicate(x,&x0);CHKERRQ(ierr);
-  ierr = VecCopy(x,x0);CHKERRQ(ierr);
+  ierr = VecDuplicate(user.y,&state);CHKERRQ(ierr);
+  ierr = VecDuplicate(user.y,&y0);CHKERRQ(ierr);
+  ierr = VecCopy(user.y,state);CHKERRQ(ierr);
+  ierr = VecCopy(user.y,y0);CHKERRQ(ierr);
 
-  /* Set solution vector with an initial guess */
-  ierr = TaoSetInitialVector(tao,x);CHKERRQ(ierr);
-  ierr = TaoSetObjectiveRoutine(tao, FormFunction, (void *)&user);CHKERRQ(ierr);
-  ierr = TaoSetGradientRoutine(tao, FormGradient, (void *)&user);CHKERRQ(ierr);
-  ierr = TaoSetConstraintsRoutine(tao, user.c, FormConstraints, (void *)&user);CHKERRQ(ierr);
+  ierr = VecDuplicate(user.u,&design);CHKERRQ(ierr);
+  ierr = VecDuplicate(user.u,&u0);CHKERRQ(ierr);
+  ierr = VecCopy(user.u,design);CHKERRQ(ierr);
+  ierr = VecCopy(user.u,u0);CHKERRQ(ierr);
 
-  ierr = TaoSetJacobianStateRoutine(tao, user.Js, user.Js, user.JsInv, FormJacobianState, (void *)&user);CHKERRQ(ierr);
+  /* Set solution vector with an initial guess */
+  ierr = TaoSetOptControlInitialVector(tao,state,design);CHKERRQ(ierr);
+  ierr = TaoSetOptControlObjectiveRoutine(tao, FormFunction, (void *)&user);CHKERRQ(ierr);
+  ierr = TaoSetOptControlGradientRoutine(tao, FormGradient, (void *)&user);CHKERRQ(ierr);
+  ierr = TaoSetOptControlConstraintsRoutine(tao, user.c, FormConstraints, (void *)&user);CHKERRQ(ierr);
 
-  ierr = TaoSetJacobianDesignRoutine(tao, user.Jd, FormJacobianDesign, (void *)&user);CHKERRQ(ierr);
+  ierr = TaoSetOptControlJacobianRoutine(tao, user.Js, user.Js, user.JsInv, user.Jd, FormJacobian, (void *)&user);CHKERRQ(ierr);
 
   ierr = TaoSetFromOptions(tao);CHKERRQ(ierr);
-  ierr = TaoSetStateDesignIS(tao,user.s_is,user.d_is);CHKERRQ(ierr);
 
   /* SOLVE THE APPLICATION */
   ierr = PetscOptionsInt("-ntests","Number of times to repeat TaoSolve","",ntests,&ntests,&flag);CHKERRQ(ierr);
   for (i=0; i<ntests; i++){
     ierr = TaoSolve(tao);CHKERRQ(ierr);
     ierr = PetscPrintf(PETSC_COMM_WORLD,"KSP Iterations = %D\n",user.ksp_its-ksp_old);CHKERRQ(ierr);
-    ierr = VecCopy(x0,x);CHKERRQ(ierr);
-    ierr = TaoSetInitialVector(tao,x);CHKERRQ(ierr);
+    ierr = VecCopy(y0,state);CHKERRQ(ierr);
+    ierr = VecCopy(u0,design);CHKERRQ(ierr);
+    ierr = TaoSetOptControlInitialVector(tao,state,design);CHKERRQ(ierr);
   }
   ierr = PetscLogStagePop();CHKERRQ(ierr);
-  ierr = PetscBarrier((PetscObject)x);CHKERRQ(ierr);
+  ierr = PetscBarrier((PetscObject)state);CHKERRQ(ierr);
   ierr = PetscPrintf(PETSC_COMM_WORLD,"KSP iterations within initialization: ");CHKERRQ(ierr);
   ierr = PetscPrintf(PETSC_COMM_WORLD,"%D\n",user.ksp_its_initial);CHKERRQ(ierr);
   ierr = PetscPrintf(PETSC_COMM_WORLD,"Total KSP iterations over %D trial(s): ",ntests);CHKERRQ(ierr);
     ierr = PetscPrintf(MPI_COMM_WORLD, "Optimization terminated with status %D.\n", reason);CHKERRQ(ierr);
   }
 
+  ierr = VecDestroy(&design);CHKERRQ(ierr);
+  ierr = VecDestroy(&u0);CHKERRQ(ierr);
+  ierr = VecDestroy(&state);CHKERRQ(ierr);
+  ierr = VecDestroy(&y0);CHKERRQ(ierr);
   ierr = TaoDestroy(&tao);CHKERRQ(ierr);
-  ierr = VecDestroy(&x);CHKERRQ(ierr);
-  ierr = VecDestroy(&x0);CHKERRQ(ierr);
   ierr = HyperbolicDestroy(&user);CHKERRQ(ierr);
   PetscFinalize();
   return 0;
    lwork = L*(u-ur).^2
    f = 1/2 * (dwork.dork + alpha*y.lwork)
 */
-PetscErrorCode FormFunction(Tao tao,Vec X,PetscReal *f,void *ptr)
+PetscErrorCode FormFunction(Tao tao,Vec Y,Vec U,PetscReal *f,void *ptr)
 {
   PetscErrorCode ierr;
   PetscReal      d1=0,d2=0;
   AppCtx         *user = (AppCtx*)ptr;
 
   PetscFunctionBegin;
-  ierr = Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);CHKERRQ(ierr);
+  ierr = VecCopy(Y,user->y);CHKERRQ(ierr);
+  ierr = VecCopy(U,user->u);CHKERRQ(ierr);
+
   ierr = MatMult(user->Q,user->y,user->dwork);CHKERRQ(ierr);
   ierr = VecAXPY(user->dwork,-1.0,user->d);CHKERRQ(ierr);
   ierr = VecDot(user->dwork,user->dwork,&d1);CHKERRQ(ierr);
     state: g_s = Q' *(Qy - d) + 0.5*alpha*L*(u-ur).^2
     design: g_d = alpha*(L'y).*(u-ur)
 */
-PetscErrorCode FormGradient(Tao tao,Vec X,Vec G,void *ptr)
+PetscErrorCode FormGradient(Tao tao,Vec Y,Vec U,Vec GY,Vec GU,void *ptr)
 {
-  PetscErrorCode ierr;
   AppCtx         *user = (AppCtx*)ptr;
+  PetscErrorCode ierr;
 
   PetscFunctionBegin;
-  ierr = Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);CHKERRQ(ierr);
+
+  ierr = VecCopy(Y,user->y);CHKERRQ(ierr);
+  ierr = VecCopy(U,user->u);CHKERRQ(ierr);
+
   ierr = MatMult(user->Q,user->y,user->dwork);CHKERRQ(ierr);
   ierr = VecAXPY(user->dwork,-1.0,user->d);CHKERRQ(ierr);
 
-  ierr = MatMult(user->QT,user->dwork,user->ywork);CHKERRQ(ierr);
+  ierr = MatMult(user->QT,user->dwork,GY);CHKERRQ(ierr);
 
-  ierr = MatMult(user->LT,user->y,user->uwork);CHKERRQ(ierr);
+  ierr = MatMult(user->LT,user->y,GU);CHKERRQ(ierr);
   ierr = VecWAXPY(user->vwork,-1.0,user->ur,user->u);CHKERRQ(ierr);
-  ierr = VecPointwiseMult(user->uwork,user->vwork,user->uwork);CHKERRQ(ierr);
-  ierr = VecScale(user->uwork,user->alpha);CHKERRQ(ierr);
+  ierr = VecPointwiseMult(GU,user->vwork,GU);CHKERRQ(ierr);
+  ierr = VecScale(GU,user->alpha);CHKERRQ(ierr);
 
   ierr = VecPointwiseMult(user->vwork,user->vwork,user->vwork);CHKERRQ(ierr);
   ierr = MatMult(user->L,user->vwork,user->lwork);CHKERRQ(ierr);
-  ierr = VecAXPY(user->ywork,0.5*user->alpha,user->lwork);CHKERRQ(ierr);
+  ierr = VecAXPY(GY,0.5*user->alpha,user->lwork);CHKERRQ(ierr);
 
-  ierr = Gather(G,user->ywork,user->state_scatter,user->uwork,user->design_scatter);CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }
 
 #undef __FUNCT__
 #define __FUNCT__ "FormFunctionGradient"
-PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *f, Vec G, void *ptr)
+PetscErrorCode FormFunctionGradient(Tao tao, Vec Y, Vec U, PetscReal *f, Vec GY, Vec GU, void *ptr)
 {
   PetscErrorCode ierr;
   PetscReal      d1,d2;
   AppCtx         *user = (AppCtx*)ptr;
 
   PetscFunctionBegin;
-  ierr = Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);CHKERRQ(ierr);
+  ierr = VecCopy(Y,user->y);CHKERRQ(ierr);
+  ierr = VecCopy(U,user->u);CHKERRQ(ierr);
   ierr = MatMult(user->Q,user->y,user->dwork);CHKERRQ(ierr);
   ierr = VecAXPY(user->dwork,-1.0,user->d);CHKERRQ(ierr);
 
-  ierr = MatMult(user->QT,user->dwork,user->ywork);CHKERRQ(ierr);
+  ierr = MatMult(user->QT,user->dwork,GY);CHKERRQ(ierr);
 
   ierr = VecDot(user->dwork,user->dwork,&d1);CHKERRQ(ierr);
 
-  ierr = MatMult(user->LT,user->y,user->uwork);CHKERRQ(ierr);
+  ierr = MatMult(user->LT,user->y,GU);CHKERRQ(ierr);
   ierr = VecWAXPY(user->vwork,-1.0,user->ur,user->u);CHKERRQ(ierr);
-  ierr = VecPointwiseMult(user->uwork,user->vwork,user->uwork);CHKERRQ(ierr);
+  ierr = VecPointwiseMult(GU,user->vwork,GU);CHKERRQ(ierr);
   ierr = VecScale(user->uwork,user->alpha);CHKERRQ(ierr);
 
   ierr = VecPointwiseMult(user->vwork,user->vwork,user->vwork);CHKERRQ(ierr);
   ierr = MatMult(user->L,user->vwork,user->lwork);CHKERRQ(ierr);
-  ierr = VecAXPY(user->ywork,0.5*user->alpha,user->lwork);CHKERRQ(ierr);
+  ierr = VecAXPY(GY,0.5*user->alpha,user->lwork);CHKERRQ(ierr);
 
   ierr = VecDot(user->y,user->lwork,&d2);CHKERRQ(ierr);
 
   *f = 0.5 * (d1 + user->alpha*d2);
-  ierr = Gather(G,user->ywork,user->state_scatter,user->uwork,user->design_scatter);CHKERRQ(ierr);
+
   PetscFunctionReturn(0);
 }
 
 /* ------------------------------------------------------------------- */
 #undef __FUNCT__
-#define __FUNCT__ "FormJacobianState"
-/* A
-MatShell object
-*/
-PetscErrorCode FormJacobianState(Tao tao, Vec X, Mat J, Mat JPre, Mat JInv, void *ptr)
+#define __FUNCT__ "FormJacobian"
+PetscErrorCode FormJacobian(Tao tao, Vec Y, Vec U, Mat Js, Mat JsPre, Mat JsInv, Mat Jd, void *ptr)
 {
   PetscErrorCode ierr;
   PetscInt       i;
   AppCtx         *user = (AppCtx*)ptr;
 
   PetscFunctionBegin;
-  ierr = Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);CHKERRQ(ierr);
+  ierr = VecCopy(Y,user->y);CHKERRQ(ierr);
+  ierr = VecCopy(U,user->u);CHKERRQ(ierr);
+
   ierr = Scatter_yi(user->u,user->ui,user->ui_scatter,user->nt);CHKERRQ(ierr);
   ierr = Scatter_uxi_uyi(user->u,user->uxi,user->uxi_scatter,user->uyi,user->uyi_scatter,user->nt);CHKERRQ(ierr);
   for (i=0; i<user->nt; i++){
   PetscFunctionReturn(0);
 }
 
-/* ------------------------------------------------------------------- */
-#undef __FUNCT__
-#define __FUNCT__ "FormJacobianDesign"
-/* B */
-PetscErrorCode FormJacobianDesign(Tao tao, Vec X, Mat J, void *ptr)
-{
-  PetscErrorCode ierr;
-  AppCtx         *user = (AppCtx*)ptr;
-
-  PetscFunctionBegin;
-  ierr = Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);CHKERRQ(ierr);
-  PetscFunctionReturn(0);
-}
 
 #undef __FUNCT__
 #define __FUNCT__ "StateMatMult"
   i = user->block_index;
   ierr = VecPointwiseMult(user->uxiwork[i],X,user->uxi[i]);CHKERRQ(ierr);
   ierr = VecPointwiseMult(user->uyiwork[i],X,user->uyi[i]);CHKERRQ(ierr);
-  ierr = Gather(user->uiwork[i],user->uxiwork[i],user->ux_scatter[i],user->uyiwork[i],user->uy_scatter[i]);CHKERRQ(ierr);
+  ierr = Gather(user->uiwork[i],user->uxiwork[i],user->ux_scatter[i]);CHKERRQ(ierr);
+  ierr = Gather(user->uiwork[i],user->uyiwork[i],user->uy_scatter[i]);CHKERRQ(ierr);
   ierr = MatMult(user->Div,user->uiwork[i],Y);CHKERRQ(ierr);
   ierr = VecAYPX(Y,user->ht,X);CHKERRQ(ierr);
   PetscFunctionReturn(0);
   ierr = MatShellGetContext(J_shell,(void**)&user);CHKERRQ(ierr);
   i = user->block_index;
   ierr = MatMult(user->Grad,X,user->uiwork[i]);CHKERRQ(ierr);
-  ierr = Scatter(user->uiwork[i],user->uxiwork[i],user->ux_scatter[i],user->uyiwork[i],user->uy_scatter[i]);CHKERRQ(ierr);
+  ierr = Scatter(user->uiwork[i],user->uxiwork[i],user->ux_scatter[i]);CHKERRQ(ierr);
+  ierr = Scatter(user->uiwork[i],user->uyiwork[i],user->uy_scatter[i]);CHKERRQ(ierr);
   ierr = VecPointwiseMult(user->uxiwork[i],user->uxi[i],user->uxiwork[i]);CHKERRQ(ierr);
   ierr = VecPointwiseMult(user->uyiwork[i],user->uyi[i],user->uyiwork[i]);CHKERRQ(ierr);
   ierr = VecWAXPY(Y,1.0,user->uxiwork[i],user->uyiwork[i]);CHKERRQ(ierr);
   for (i=0; i<user->nt; i++){
     ierr = VecPointwiseMult(user->uxiwork[i],user->yi[i],user->uxiwork[i]);CHKERRQ(ierr);
     ierr = VecPointwiseMult(user->uyiwork[i],user->yi[i],user->uyiwork[i]);CHKERRQ(ierr);
-    ierr = Gather(user->uiwork[i],user->uxiwork[i],user->ux_scatter[i],user->uyiwork[i],user->uy_scatter[i]);CHKERRQ(ierr);
+    ierr = Gather(user->uiwork[i],user->uxiwork[i],user->ux_scatter[i]);CHKERRQ(ierr);
+    ierr = Gather(user->uiwork[i],user->uyiwork[i],user->uy_scatter[i]);CHKERRQ(ierr);
     ierr = MatMult(user->Div,user->uiwork[i],user->ziwork[i]);CHKERRQ(ierr);
     ierr = VecScale(user->ziwork[i],user->ht);CHKERRQ(ierr);
   }
   ierr = Scatter_yi(X,user->yiwork,user->yi_scatter,user->nt);CHKERRQ(ierr);
   for (i=0; i<user->nt; i++){
     ierr = MatMult(user->Grad,user->yiwork[i],user->uiwork[i]);CHKERRQ(ierr);
-    ierr = Scatter(user->uiwork[i],user->uxiwork[i],user->ux_scatter[i],user->uyiwork[i],user->uy_scatter[i]);CHKERRQ(ierr);
+    ierr = Scatter(user->uiwork[i],user->uxiwork[i],user->ux_scatter[i]);CHKERRQ(ierr);
+    ierr = Scatter(user->uiwork[i],user->uyiwork[i],user->uy_scatter[i]);CHKERRQ(ierr);
     ierr = VecPointwiseMult(user->uxiwork[i],user->yi[i],user->uxiwork[i]);CHKERRQ(ierr);
     ierr = VecPointwiseMult(user->uyiwork[i],user->yi[i],user->uyiwork[i]);CHKERRQ(ierr);
-    ierr = Gather(user->uiwork[i],user->uxiwork[i],user->ux_scatter[i],user->uyiwork[i],user->uy_scatter[i]);CHKERRQ(ierr);
+    ierr = Gather(user->uiwork[i],user->uxiwork[i],user->ux_scatter[i]);CHKERRQ(ierr);
+    ierr = Gather(user->uiwork[i],user->uyiwork[i],user->uy_scatter[i]);CHKERRQ(ierr);
     ierr = VecScale(user->uiwork[i],user->ht);CHKERRQ(ierr);
   }
   ierr = Gather_yi(Y,user->uiwork,user->ui_scatter,user->nt);CHKERRQ(ierr);
 
 #undef __FUNCT__
 #define __FUNCT__ "FormConstraints"
-PetscErrorCode FormConstraints(Tao tao, Vec X, Vec C, void *ptr)
+PetscErrorCode FormConstraints(Tao tao, Vec Y, Vec U, Vec CE, void *ptr)
 {
-  /* con = Ay - q, A = [C(u1)  0     0     ...   0;
-                         -M  C(u2)   0     ...   0;
-                          0   -M   C(u3)   ...   0;
+  /* con = Ay - q, A = [CE(u1)  0     0     ...   0;
+                         -M  CE(u2)   0     ...   0;
+                          0   -M   CE(u3)   ...   0;
                                       ...         ;
-                          0    ...      -M C(u_nt)]
-     C(u) = eye + ht*Div*[diag(u1); diag(u2)]       */
+                          0    ...      -M CE(u_nt)]
+     CE(u) = eye + ht*Div*[diag(u1); diag(u2)]       */
   PetscErrorCode ierr;
   PetscInt       i;
   AppCtx         *user = (AppCtx*)ptr;
 
   PetscFunctionBegin;
-  ierr = Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);CHKERRQ(ierr);
+  ierr = VecCopy(Y,user->y);CHKERRQ(ierr);
+  ierr = VecCopy(U,user->u);CHKERRQ(ierr);
   ierr = Scatter_yi(user->y,user->yi,user->yi_scatter,user->nt);CHKERRQ(ierr);
   ierr = Scatter_uxi_uyi(user->u,user->uxi,user->uxi_scatter,user->uyi,user->uyi_scatter,user->nt);
 
     ierr = VecAXPY(user->yiwork[i],-1.0,user->ziwork[i-1]);CHKERRQ(ierr);
   }
 
-  ierr = Gather_yi(C,user->yiwork,user->yi_scatter,user->nt);CHKERRQ(ierr);
-  ierr = VecAXPY(C,-1.0,user->q);CHKERRQ(ierr);
-
-  PetscFunctionReturn(0);
-}
-
-
-#undef __FUNCT__
-#define __FUNCT__ "Scatter"
-PetscErrorCode Scatter(Vec x, Vec state, VecScatter s_scat, Vec design, VecScatter d_scat)
-{
-  PetscErrorCode ierr;
+  ierr = Gather_yi(CE,user->yiwork,user->yi_scatter,user->nt);CHKERRQ(ierr);
+  ierr = VecAXPY(CE,-1.0,user->q);CHKERRQ(ierr);
 
-  PetscFunctionBegin;
-  ierr = VecScatterBegin(s_scat,x,state,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
-  ierr = VecScatterEnd(s_scat,x,state,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
-  ierr = VecScatterBegin(d_scat,x,design,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
-  ierr = VecScatterEnd(d_scat,x,design,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }
 
   PetscFunctionReturn(0);
 }
 
-#undef __FUNCT__
-#define __FUNCT__ "Gather"
-PetscErrorCode Gather(Vec x, Vec state, VecScatter s_scat, Vec design, VecScatter d_scat)
-{
-  PetscErrorCode ierr;
-
-  PetscFunctionBegin;
-  ierr = VecScatterBegin(s_scat,state,x,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
-  ierr = VecScatterEnd(s_scat,state,x,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
-  ierr = VecScatterBegin(d_scat,design,x,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
-  ierr = VecScatterEnd(d_scat,design,x,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
-  PetscFunctionReturn(0);
-}
 
 #undef __FUNCT__
 #define __FUNCT__ "Gather_uxi_uyi"
 }
 
 #undef __FUNCT__
+#define __FUNCT__ "Scatter"
+PetscErrorCode Scatter(Vec v, Vec subv, VecScatter scatter)
+{
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+
+  ierr = VecScatterBegin(scatter,v,subv,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
+  ierr = VecScatterEnd(scatter,v,subv,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
+
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "Gather"
+PetscErrorCode Gather(Vec v, Vec subv, VecScatter scatter)
+{
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+
+  ierr = VecScatterBegin(scatter,subv,v,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
+  ierr = VecScatterEnd(scatter,subv,v,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
+
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
 #define __FUNCT__ "HyperbolicInitialize"
 PetscErrorCode HyperbolicInitialize(AppCtx *user)
 {
   ierr = VecScale(bc,1.0/(h*h*sum));CHKERRQ(ierr);
 
   /* Create scatter from y to y_1,y_2,...,y_nt */
-  /*  TODO: Reorder for better parallelism. (This will require reordering Q and L as well.) */
+  /*  Should reorder for better parallelism. (This will require reordering Q and L as well.) */
   ierr = PetscMalloc1(user->nt*user->mx*user->mx,&user->yi_scatter);
   ierr = VecCreate(PETSC_COMM_WORLD,&yi);CHKERRQ(ierr);
   ierr = VecSetSizes(yi,PETSC_DECIDE,user->mx*user->mx);CHKERRQ(ierr);
   }
 
   /* Create scatter from u to ux_1,uy_1,ux_2,uy_2,...,ux_nt,uy_nt */
-  /*  TODO: reorder for better parallelism */
+  /*  Should reorder for better parallelism */
   ierr = PetscMalloc1(user->nt*user->mx*user->mx,&user->uxi_scatter);CHKERRQ(ierr);
   ierr = PetscMalloc1(user->nt*user->mx*user->mx,&user->uyi_scatter);CHKERRQ(ierr);
   ierr = PetscMalloc1(user->nt*user->mx*user->mx,&user->ux_scatter);CHKERRQ(ierr);
   ierr = VecDestroy(&user->dwork);CHKERRQ(ierr);
   ierr = VecDestroy(&user->lwork);CHKERRQ(ierr);
   ierr = VecDestroy(&user->js_diag);CHKERRQ(ierr);
-  ierr = ISDestroy(&user->s_is);CHKERRQ(ierr);
-  ierr = ISDestroy(&user->d_is);CHKERRQ(ierr);
-  ierr = VecScatterDestroy(&user->state_scatter);CHKERRQ(ierr);
-  ierr = VecScatterDestroy(&user->design_scatter);CHKERRQ(ierr);
   for (i=0; i<user->nt; i++){
     ierr = VecScatterDestroy(&user->uxi_scatter[i]);CHKERRQ(ierr);
     ierr = VecScatterDestroy(&user->uyi_scatter[i]);CHKERRQ(ierr);
 PetscErrorCode HyperbolicMonitor(Tao tao, void *ptr)
 {
   PetscErrorCode ierr;
-  Vec            X;
+  Vec            Y,U;
   PetscReal      unorm,ynorm;
   AppCtx         *user = (AppCtx*)ptr;
 
   PetscFunctionBegin;
-  ierr = TaoGetSolutionVector(tao,&X);CHKERRQ(ierr);
-  ierr = Scatter(X,user->ywork,user->state_scatter,user->uwork,user->design_scatter);CHKERRQ(ierr);
-  ierr = VecAXPY(user->ywork,-1.0,user->ytrue);CHKERRQ(ierr);
-  ierr = VecAXPY(user->uwork,-1.0,user->utrue);CHKERRQ(ierr);
-  ierr = VecNorm(user->uwork,NORM_2,&unorm);CHKERRQ(ierr);
-  ierr = VecNorm(user->ywork,NORM_2,&ynorm);CHKERRQ(ierr);
+
+  ierr = TaoGetOptControlSolutionVector(tao,&Y,&U);CHKERRQ(ierr);
+
+  ierr = VecAXPY(Y,-1.0,user->ytrue);CHKERRQ(ierr);
+  ierr = VecNorm(Y,NORM_2,&ynorm);CHKERRQ(ierr);
+
+  ierr = VecAXPY(U,-1.0,user->utrue);CHKERRQ(ierr);
+  ierr = VecNorm(U,NORM_2,&unorm);CHKERRQ(ierr);
   ierr = PetscPrintf(MPI_COMM_WORLD, "||u-ut||=%g ||y-yt||=%g\n",(double)unorm,(double)ynorm);CHKERRQ(ierr);
+
   PetscFunctionReturn(0);
 }

src/tao/pde_constrained/examples/tutorials/parabolic.c

   PetscInt ndata; /*  Number of data points per sample */
   PetscInt ns; /*  Number of samples */
   PetscInt *sample_times; /*  Times of samples */
-  IS       s_is;
-  IS       d_is;
 
   VecScatter state_scatter;
   VecScatter design_scatter;
 } AppCtx;
 
 
-PetscErrorCode FormFunction(Tao, Vec, PetscReal*, void*);
-PetscErrorCode FormGradient(Tao, Vec, Vec, void*);
-PetscErrorCode FormFunctionGradient(Tao, Vec, PetscReal*, Vec, void*);
-PetscErrorCode FormJacobianState(Tao, Vec, Mat, Mat, Mat, void*);
-PetscErrorCode FormJacobianDesign(Tao, Vec, Mat, void*);
-PetscErrorCode FormConstraints(Tao, Vec, Vec, void*);
-PetscErrorCode FormHessian(Tao, Vec, Mat, Mat, void*);
-PetscErrorCode Gather(Vec x, Vec state, VecScatter s_scat, Vec design, VecScatter d_scat);
-PetscErrorCode Scatter(Vec x, Vec state, VecScatter s_scat, Vec design, VecScatter d_scat);
+PetscErrorCode FormFunction(Tao, Vec, Vec, PetscReal*, void*);
+PetscErrorCode FormGradient(Tao, Vec, Vec, Vec, Vec, void*);
+PetscErrorCode FormFunctionGradient(Tao, Vec, Vec, PetscReal*, Vec, Vec, void*);
+PetscErrorCode FormJacobian(Tao, Vec, Vec, Mat, Mat, Mat, Mat, void*);
+PetscErrorCode FormConstraints(Tao, Vec, Vec, Vec, void*);
+
 PetscErrorCode ParabolicInitialize(AppCtx *user);
 PetscErrorCode ParabolicDestroy(AppCtx *user);
 PetscErrorCode ParabolicMonitor(Tao, void*);
 PetscErrorCode DesignMatMult(Mat,Vec,Vec);
 PetscErrorCode DesignMatMultTranspose(Mat,Vec,Vec);
 
-PetscErrorCode Gather_i(Vec,Vec*,VecScatter*,PetscInt);
-PetscErrorCode Scatter_i(Vec,Vec*,VecScatter*,PetscInt);
+PetscErrorCode Gather(Vec,Vec*,VecScatter*,PetscInt);
+PetscErrorCode Scatter(Vec,Vec*,VecScatter*,PetscInt);
 
 static  char help[]="";
 
 int main(int argc, char **argv)
 {
   PetscErrorCode     ierr;
-  Vec                x,x0;
+  Vec                state, design, y0, u0;
   Tao                tao;
   TaoConvergedReason reason;
   AppCtx             user;
-  IS                 is_allstate,is_alldesign;
-  PetscInt           lo,hi,hi2,lo2,ksp_old;
+  PetscInt           ksp_old;
   PetscBool          flag;
   PetscInt           ntests = 1;
   PetscInt           i;
-  int                stages[1];
+  PetscInt           stages[1];
 
   PetscInitialize(&argc, &argv, (char*)0,help);
 
   user.n = user.m*(user.nt+1); /*  number of variables */
   user.ht = (PetscReal)1/user.nt;
 
-  ierr = VecCreate(PETSC_COMM_WORLD,&user.u);CHKERRQ(ierr);
   ierr = VecCreate(PETSC_COMM_WORLD,&user.y);CHKERRQ(ierr);
+  ierr = VecCreate(PETSC_COMM_WORLD,&user.u);CHKERRQ(ierr);
   ierr = VecCreate(PETSC_COMM_WORLD,&user.c);CHKERRQ(ierr);
-  ierr = VecSetSizes(user.u,PETSC_DECIDE,user.n-user.m*user.nt);CHKERRQ(ierr);
   ierr = VecSetSizes(user.y,PETSC_DECIDE,user.m*user.nt);CHKERRQ(ierr);
+  ierr = VecSetSizes(user.u,PETSC_DECIDE,user.n-user.m*user.nt);CHKERRQ(ierr);
   ierr = VecSetSizes(user.c,PETSC_DECIDE,user.m*user.nt);CHKERRQ(ierr);
-  ierr = VecSetFromOptions(user.u);CHKERRQ(ierr);
   ierr = VecSetFromOptions(user.y);CHKERRQ(ierr);
+  ierr = VecSetFromOptions(user.u);CHKERRQ(ierr);
   ierr = VecSetFromOptions(user.c);CHKERRQ(ierr);
 
-  /* Create scatters for reduced spaces.
-     If the state vector y and design vector u are partitioned as
-     [y_1; y_2; ...; y_np] and [u_1; u_2; ...; u_np] (with np = # of processors),
-     then the solution vector x is organized as
-     [y_1; u_1; y_2; u_2; ...; y_np; u_np].
-     The index sets user.s_is and user.d_is correspond to the indices of the
-     state and design variables owned by the current processor.
-  */
-  ierr = VecCreate(PETSC_COMM_WORLD,&x);CHKERRQ(ierr);
-
-  ierr = VecGetOwnershipRange(user.y,&lo,&hi);CHKERRQ(ierr);
-  ierr = VecGetOwnershipRange(user.u,&lo2,&hi2);CHKERRQ(ierr);
-
-  ierr = ISCreateStride(PETSC_COMM_SELF,hi-lo,lo,1,&is_allstate);CHKERRQ(ierr);
-  ierr = ISCreateStride(PETSC_COMM_SELF,hi-lo,lo+lo2,1,&user.s_is);CHKERRQ(ierr);
-
-  ierr = ISCreateStride(PETSC_COMM_SELF,hi2-lo2,lo2,1,&is_alldesign);CHKERRQ(ierr);
-  ierr = ISCreateStride(PETSC_COMM_SELF,hi2-lo2,hi+lo2,1,&user.d_is);CHKERRQ(ierr);
-
-  ierr = VecSetSizes(x,hi-lo+hi2-lo2,user.n);CHKERRQ(ierr);
-  ierr = VecSetFromOptions(x);CHKERRQ(ierr);
-
-  ierr = VecScatterCreate(x,user.s_is,user.y,is_allstate,&user.state_scatter);CHKERRQ(ierr);
-  ierr = VecScatterCreate(x,user.d_is,user.u,is_alldesign,&user.design_scatter);CHKERRQ(ierr);
-  ierr = ISDestroy(&is_alldesign);CHKERRQ(ierr);
-  ierr = ISDestroy(&is_allstate);CHKERRQ(ierr);
 
   /* Create TAO solver and set desired solution method */
   ierr = TaoCreate(PETSC_COMM_WORLD,&tao);CHKERRQ(ierr);
   /* Set up initial vectors and matrices */
   ierr = ParabolicInitialize(&user);CHKERRQ(ierr);
 
-  ierr = Gather(x,user.y,user.state_scatter,user.u,user.design_scatter);CHKERRQ(ierr);
-  ierr = VecDuplicate(x,&x0);CHKERRQ(ierr);
-  ierr = VecCopy(x,x0);CHKERRQ(ierr);
+  ierr = VecDuplicate(user.y,&state);CHKERRQ(ierr);
+  ierr = VecDuplicate(user.y,&y0);CHKERRQ(ierr);
+  ierr = VecCopy(user.y,state);CHKERRQ(ierr);
+  ierr = VecCopy(user.y,y0);CHKERRQ(ierr);
+  ierr = VecDuplicate(user.u,&design);CHKERRQ(ierr);
+  ierr = VecDuplicate(user.u,&u0);CHKERRQ(ierr);
+  ierr = VecCopy(user.u,design);CHKERRQ(ierr);
+  ierr = VecCopy(user.u,u0);CHKERRQ(ierr);
 
   /* Set solution vector with an initial guess */
-  ierr = TaoSetInitialVector(tao,x);CHKERRQ(ierr);
-  ierr = TaoSetObjectiveRoutine(tao, FormFunction, (void *)&user);CHKERRQ(ierr);
-  ierr = TaoSetGradientRoutine(tao, FormGradient, (void *)&user);CHKERRQ(ierr);
-  ierr = TaoSetConstraintsRoutine(tao, user.c, FormConstraints, (void *)&user);CHKERRQ(ierr);
+  ierr = TaoSetOptControlInitialVector(tao,state,design);CHKERRQ(ierr);
+  ierr = TaoSetOptControlObjectiveRoutine(tao, FormFunction, (void *)&user);CHKERRQ(ierr);
+  ierr = TaoSetOptControlGradientRoutine(tao, FormGradient, (void *)&user);CHKERRQ(ierr);
+  ierr = TaoSetOptControlConstraintsRoutine(tao, user.c, FormConstraints, (void *)&user);CHKERRQ(ierr);
 
-  ierr = TaoSetJacobianStateRoutine(tao, user.Js, user.JsBlockPrec, user.JsInv, FormJacobianState, (void *)&user);CHKERRQ(ierr);
-  ierr = TaoSetJacobianDesignRoutine(tao, user.Jd, FormJacobianDesign, (void *)&user);CHKERRQ(ierr);
+  ierr = TaoSetOptControlJacobianRoutine(tao, user.Js, user.JsBlockPrec, user.JsInv, user.Jd, FormJacobian, (void *)&user);CHKERRQ(ierr);
 
   ierr = TaoSetFromOptions(tao);CHKERRQ(ierr);
-  ierr = TaoSetStateDesignIS(tao,user.s_is,user.d_is);CHKERRQ(ierr);
 
  /* SOLVE THE APPLICATION */
   ierr = PetscOptionsInt("-ntests","Number of times to repeat TaoSolve","",ntests,&ntests,&flag);CHKERRQ(ierr);
   for (i=0; i<ntests; i++){
     ierr = TaoSolve(tao);CHKERRQ(ierr);
     ierr = PetscPrintf(PETSC_COMM_WORLD,"KSP Iterations = %D\n",user.ksp_its-ksp_old);CHKERRQ(ierr);
-    ierr = VecCopy(x0,x);CHKERRQ(ierr);
-    ierr = TaoSetInitialVector(tao,x);CHKERRQ(ierr);
+    ierr = VecCopy(y0,state);CHKERRQ(ierr);
+    ierr = VecCopy(u0,design);CHKERRQ(ierr);
   }
   ierr = PetscLogStagePop();CHKERRQ(ierr);
-  ierr = PetscBarrier((PetscObject)x);CHKERRQ(ierr);
+  ierr = PetscBarrier((PetscObject)state);CHKERRQ(ierr);
   ierr = PetscPrintf(PETSC_COMM_WORLD,"KSP iterations within initialization: ");CHKERRQ(ierr);
   ierr = PetscPrintf(PETSC_COMM_WORLD,"%D\n",user.ksp_its_initial);CHKERRQ(ierr);
   ierr = PetscPrintf(PETSC_COMM_WORLD,"Total KSP iterations over %D trial(s): ",ntests);CHKERRQ(ierr);
     ierr = PetscPrintf(MPI_COMM_WORLD, "Optimization terminated with status %D.\n", reason);CHKERRQ(ierr);
   }
   ierr = TaoDestroy(&tao);CHKERRQ(ierr);
-  ierr = VecDestroy(&x);CHKERRQ(ierr);
-  ierr = VecDestroy(&x0);CHKERRQ(ierr);
+  ierr = VecDestroy(&state);CHKERRQ(ierr);
+  ierr = VecDestroy(&design);CHKERRQ(ierr);
+  ierr = VecDestroy(&y0);CHKERRQ(ierr);
+  ierr = VecDestroy(&u0);CHKERRQ(ierr);
+
   ierr = ParabolicDestroy(&user);CHKERRQ(ierr);
 
   PetscFinalize();
    lwork = L*(u-ur)
    f = 1/2 * (dwork.dork + alpha*lwork.lwork)
 */
-PetscErrorCode FormFunction(Tao tao,Vec X,PetscReal *f,void *ptr)
+PetscErrorCode FormFunction(Tao tao,Vec Y,Vec U,PetscReal *f,void *ptr)
 {
   PetscErrorCode ierr;
   PetscReal      d1=0,d2=0;
   AppCtx         *user = (AppCtx*)ptr;
 
   PetscFunctionBegin;
-  ierr = Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);CHKERRQ(ierr);
-  ierr = Scatter_i(user->y,user->yi,user->yi_scatter,user->nt);CHKERRQ(ierr);
+
+  ierr = VecCopy(Y,user->y);CHKERRQ(ierr);
+  ierr = VecCopy(U,user->u);CHKERRQ(ierr);
+  ierr = Scatter(user->y,user->yi,user->yi_scatter,user->nt);CHKERRQ(ierr);
   for (j=0; j<user->ns; j++){
     i = user->sample_times[j];
     ierr = MatMult(user->Qblock,user->yi[i],user->di[j]);CHKERRQ(ierr);
   }
-  ierr = Gather_i(user->dwork,user->di,user->di_scatter,user->ns);CHKERRQ(ierr);
+  ierr = Gather(user->dwork,user->di,user->di_scatter,user->ns);CHKERRQ(ierr);
   ierr = VecAXPY(user->dwork,-1.0,user->d);CHKERRQ(ierr);
   ierr = VecDot(user->dwork,user->dwork,&d1);CHKERRQ(ierr);
 
     state: g_s = Q' *(Qy - d)
     design: g_d = alpha*L'*L*(u-ur)
 */
-PetscErrorCode FormGradient(Tao tao,Vec X,Vec G,void *ptr)
+PetscErrorCode FormGradient(Tao tao,Vec Y,Vec U,Vec GY,Vec GU,void *ptr)
 {
   PetscErrorCode ierr;
   PetscInt       i,j;
   AppCtx         *user = (AppCtx*)ptr;
 
   PetscFunctionBegin;
-  ierr = Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);CHKERRQ(ierr);
-  ierr = Scatter_i(user->y,user->yi,user->yi_scatter,user->nt);CHKERRQ(ierr);
+  ierr = VecCopy(Y,user->y);CHKERRQ(ierr);
+  ierr = VecCopy(U,user->u);CHKERRQ(ierr);
+  ierr = Scatter(user->y,user->yi,user->yi_scatter,user->nt);CHKERRQ(ierr);
   for (j=0; j<user->ns; j++){
     i = user->sample_times[j];
     ierr = MatMult(user->Qblock,user->yi[i],user->di[j]);CHKERRQ(ierr);
   }
-  ierr = Gather_i(user->dwork,user->di,user->di_scatter,user->ns);CHKERRQ(ierr);
+  ierr = Gather(user->dwork,user->di,user->di_scatter,user->ns);CHKERRQ(ierr);
   ierr = VecAXPY(user->dwork,-1.0,user->d);CHKERRQ(ierr);
-  ierr = Scatter_i(user->dwork,user->di,user->di_scatter,user->ns);CHKERRQ(ierr);
-  ierr = VecSet(user->ywork,0.0);CHKERRQ(ierr);
-  ierr = Scatter_i(user->ywork,user->yiwork,user->yi_scatter,user->nt);CHKERRQ(ierr);
+  ierr = Scatter(user->dwork,user->di,user->di_scatter,user->ns);CHKERRQ(ierr);
+  ierr = VecSet(GY,0.0);CHKERRQ(ierr);
+  ierr = Scatter(GY,user->yiwork,user->yi_scatter,user->nt);CHKERRQ(ierr);
   for (j=0; j<user->ns; j++){
     i = user->sample_times[j];
     ierr = MatMult(user->QblockT,user->di[j],user->yiwork[i]);CHKERRQ(ierr);
   }
-  ierr = Gather_i(user->ywork,user->yiwork,user->yi_scatter,user->nt);CHKERRQ(ierr);
+  ierr = Gather(GY,user->yiwork,user->yi_scatter,user->nt);CHKERRQ(ierr);
+
+  ierr = VecWAXPY(GU,-1.0,user->ur,user->u);CHKERRQ(ierr);
+  ierr = MatMult(user->L,GU,user->lwork);CHKERRQ(ierr);
+  ierr = MatMult(user->LT,user->lwork,GU);CHKERRQ(ierr);
+  ierr = VecScale(GU, user->alpha);CHKERRQ(ierr);
 
-  ierr = VecWAXPY(user->uwork,-1.0,user->ur,user->u);CHKERRQ(ierr);
-  ierr = MatMult(user->L,user->uwork,user->lwork);CHKERRQ(ierr);
-  ierr = MatMult(user->LT,user->lwork,user->uwork);CHKERRQ(ierr);
-  ierr = VecScale(user->uwork, user->alpha);CHKERRQ(ierr);
-  ierr = Gather(G,user->ywork,user->state_scatter,user->uwork,user->design_scatter);CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }
-
+/*
 #undef __FUNCT__
 #define __FUNCT__ "FormFunctionGradient"
 PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *f, Vec G, void *ptr)
 
   PetscFunctionBegin;
   ierr = Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);CHKERRQ(ierr);
-  ierr = Scatter_i(user->y,user->yi,user->yi_scatter,user->nt);CHKERRQ(ierr);
+  ierr = Scatter(user->y,user->yi,user->yi_scatter,user->nt);CHKERRQ(ierr);
   for (j=0; j<user->ns; j++){
     i = user->sample_times[j];
     ierr = MatMult(user->Qblock,user->yi[i],user->di[j]);CHKERRQ(ierr);
   }
-  ierr = Gather_i(user->dwork,user->di,user->di_scatter,user->ns);CHKERRQ(ierr);
+  ierr = Gather(user->dwork,user->di,user->di_scatter,user->ns);CHKERRQ(ierr);
   ierr = VecAXPY(user->dwork,-1.0,user->d);CHKERRQ(ierr);
   ierr = VecDot(user->dwork,user->dwork,&d1);CHKERRQ(ierr);
-  ierr = Scatter_i(user->dwork,user->di,user->di_scatter,user->ns);CHKERRQ(ierr);
+  ierr = Scatter(user->dwork,user->di,user->di_scatter,user->ns);CHKERRQ(ierr);
   ierr = VecSet(user->ywork,0.0);CHKERRQ(ierr);
-  ierr = Scatter_i(user->ywork,user->yiwork,user->yi_scatter,user->nt);CHKERRQ(ierr);
+  ierr = Scatter(user->ywork,user->yiwork,user->yi_scatter,user->nt);CHKERRQ(ierr);
   for (j=0; j<user->ns; j++){
     i = user->sample_times[j];
     ierr = MatMult(user->QblockT,user->di[j],user->yiwork[i]);CHKERRQ(ierr);
   }
-  ierr = Gather_i(user->ywork,user->yiwork,user->yi_scatter,user->nt);CHKERRQ(ierr);
+  ierr = Gather(user->ywork,user->yiwork,user->yi_scatter,user->nt);CHKERRQ(ierr);
 
   ierr = VecWAXPY(user->uwork,-1.0,user->ur,user->u);CHKERRQ(ierr);
   ierr = MatMult(user->L,user->uwork,user->lwork);CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }
 
+ */
+
 /* ------------------------------------------------------------------- */
 #undef __FUNCT__
-#define __FUNCT__ "FormJacobianState"
+#define __FUNCT__ "FormJacobian"
 /* A
 MatShell object
 */
-PetscErrorCode FormJacobianState(Tao tao, Vec X, Mat J, Mat JPre, Mat JInv, void *ptr)
+PetscErrorCode FormJacobian(Tao tao,Vec Y,Vec U,Mat J,Mat JPre,Mat JInv,Mat Jd,void *ptr)
 {
   PetscErrorCode ierr;
   AppCtx         *user = (AppCtx*)ptr;
 
   PetscFunctionBegin;
-  ierr = Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);CHKERRQ(ierr);
+
+  ierr = VecCopy(Y,user->y);CHKERRQ(ierr);
+  ierr = VecCopy(U,user->u);CHKERRQ(ierr);
   ierr = VecSet(user->uwork,0);CHKERRQ(ierr);
   ierr = VecAXPY(user->uwork,-1.0,user->u);CHKERRQ(ierr);
   ierr = VecExp(user->uwork);CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }
 
-/* ------------------------------------------------------------------- */
-#undef __FUNCT__
-#define __FUNCT__ "FormJacobianDesign"
-/* B */
-PetscErrorCode FormJacobianDesign(Tao tao, Vec X, Mat J, void *ptr)
-{
-  PetscErrorCode ierr;
-  AppCtx         *user = (AppCtx*)ptr;
-
-  PetscFunctionBegin;
-  ierr = Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);CHKERRQ(ierr);
-  PetscFunctionReturn(0);
-}
-
 #undef __FUNCT__
 #define __FUNCT__ "StateMatMult"
 PetscErrorCode StateMatMult(Mat J_shell, Vec X, Vec Y)
 
   PetscFunctionBegin;
   ierr = MatShellGetContext(J_shell,(void**)&user);CHKERRQ(ierr);
-  ierr = Scatter_i(X,user->yi,user->yi_scatter,user->nt);CHKERRQ(ierr);
+  ierr = Scatter(X,user->yi,user->yi_scatter,user->nt);CHKERRQ(ierr);
   ierr = MatMult(user->JsBlock,user->yi[0],user->yiwork[0]);CHKERRQ(ierr);
   for (i=1; i<user->nt; i++){
     ierr = MatMult(user->JsBlock,user->yi[i],user->yiwork[i]);CHKERRQ(ierr);
     ierr = VecAXPY(user->yiwork[i],-1.0,user->yi[i-1]);CHKERRQ(ierr);
   }
-  ierr = Gather_i(Y,user->yiwork,user->yi_scatter,user->nt);CHKERRQ(ierr);
+  ierr = Gather(Y,user->yiwork,user->yi_scatter,user->nt);CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }
 
 
   PetscFunctionBegin;
   ierr = MatShellGetContext(J_shell,(void**)&user);CHKERRQ(ierr);
-  ierr = Scatter_i(X,user->yi,user->yi_scatter,user->nt);CHKERRQ(ierr);
+  ierr = Scatter(X,user->yi,user->yi_scatter,user->nt);CHKERRQ(ierr);
   for (i=0; i<user->nt-1; i++){
     ierr = MatMult(user->JsBlock,user->yi[i],user->yiwork[i]);CHKERRQ(ierr);
     ierr = VecAXPY(user->yiwork[i],-1.0,user->yi[i+1]);CHKERRQ(ierr);
   }
   i = user->nt-1;
   ierr = MatMult(user->JsBlock,user->yi[i],user->yiwork[i]);CHKERRQ(ierr);
-  ierr = Gather_i(Y,user->yiwork,user->yi_scatter,user->nt);CHKERRQ(ierr);
+  ierr = Gather(Y,user->yiwork,user->yi_scatter,user->nt);CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }
 
   /* (sdiag(1./((Av*(1./v)).^2)) * (Av * (sdiag(1./v) * b))) */
   ierr = VecPointwiseMult(user->Swork,user->Twork,user->Swork);CHKERRQ(ierr);
 
-  ierr = Scatter_i(user->y,user->yi,user->yi_scatter,user->nt);CHKERRQ(ierr);
+  ierr = Scatter(user->y,user->yi,user->yi_scatter,user->nt);CHKERRQ(ierr);
   for (i=0; i<user->nt; i++){
     /* (sdiag(Grad*y(:,i)) */
     ierr = MatMult(user->Grad,user->yi[i],user->Twork);CHKERRQ(ierr);
     ierr = MatMult(user->Div,user->Twork,user->yiwork[i]);CHKERRQ(ierr);
     ierr = VecScale(user->yiwork[i],user->ht);CHKERRQ(ierr);
   }
-  ierr = Gather_i(Y,user->yiwork,user->yi_scatter,user->nt);CHKERRQ(ierr);
+  ierr = Gather(Y,user->yiwork,user->yi_scatter,user->nt);CHKERRQ(ierr);
 
   PetscFunctionReturn(0);
 }
   ierr = VecReciprocal(user->Rwork);CHKERRQ(ierr);
 
   ierr = VecSet(Y,0.0);CHKERRQ(ierr);
-  ierr = Scatter_i(user->y,user->yi,user->yi_scatter,user->nt);CHKERRQ(ierr);
-  ierr = Scatter_i(X,user->yiwork,user->yi_scatter,user->nt);CHKERRQ(ierr);
+  ierr = Scatter(user->y,user->yi,user->yi_scatter,user->nt);CHKERRQ(ierr);
+  ierr = Scatter(X,user->yiwork,user->yi_scatter,user->nt);CHKERRQ(ierr);
   for (i=0; i<user->nt; i++){
     /* (Div' * b(:,i)) */
     ierr = MatMult(user->Grad,user->yiwork[i],user->Swork);CHKERRQ(ierr);
     ierr = KSPSetTolerances(user->solver,tau,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr);
   }
 
-  ierr = Scatter_i(X,user->yi,user->yi_scatter,user->nt);CHKERRQ(ierr);
+  ierr = Scatter(X,user->yi,user->yi_scatter,user->nt);CHKERRQ(ierr);
   ierr = KSPSolve(user->solver,user->yi[0],user->yiwork[0]);CHKERRQ(ierr);
   ierr = KSPGetIterationNumber(user->solver,&its);CHKERRQ(ierr);
   user->ksp_its = user->ksp_its + its;
     ierr = KSPGetIterationNumber(user->solver,&its);CHKERRQ(ierr);
     user->ksp_its = user->ksp_its + its;
   }
-  ierr = Gather_i(Y,user->yiwork,user->yi_scatter,user->nt);CHKERRQ(ierr);
+  ierr = Gather(Y,user->yiwork,user->yi_scatter,user->nt);CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }
 
     tau = user->lcl->tau[user->lcl->solve_type];
     ierr = KSPSetTolerances(user->solver,tau,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr);
   }
-  ierr = Scatter_i(X,user->yi,user->yi_scatter,user->nt);CHKERRQ(ierr);
+  ierr = Scatter(X,user->yi,user->yi_scatter,user->nt);CHKERRQ(ierr);
 
   i = user->nt - 1;
   ierr = KSPSolve(user->solver,user->yi[i],user->yiwork[i]);CHKERRQ(ierr);
 
   }
 
-  ierr = Gather_i(Y,user->yiwork,user->yi_scatter,user->nt);CHKERRQ(ierr);
+  ierr = Gather(Y,user->yiwork,user->yi_scatter,user->nt);CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }
 
 
 #undef __FUNCT__
 #define __FUNCT__ "FormConstraints"
-PetscErrorCode FormConstraints(Tao tao, Vec X, Vec C, void *ptr)
+PetscErrorCode FormConstraints(Tao tao, Vec Y, Vec U, Vec CE, void *ptr)
 {
   /* con = Ay - q, A = [B  0  0 ... 0;
                        -I  B  0 ... 0;
   AppCtx         *user = (AppCtx*)ptr;
 
   PetscFunctionBegin;
-  ierr = Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);CHKERRQ(ierr);
-  ierr = Scatter_i(user->y,user->yi,user->yi_scatter,user->nt);CHKERRQ(ierr);
+
+  ierr = VecCopy(Y,user->y);CHKERRQ(ierr);
+  ierr = VecCopy(U,user->u);CHKERRQ(ierr);
+  ierr = Scatter(user->y,user->yi,user->yi_scatter,user->nt);CHKERRQ(ierr);
   ierr = MatMult(user->JsBlock,user->yi[0],user->yiwork[0]);CHKERRQ(ierr);
   for (i=1; i<user->nt; i++){
     ierr = MatMult(user->JsBlock,user->yi[i],user->yiwork[i]);CHKERRQ(ierr);
     ierr = VecAXPY(user->yiwork[i],-1.0,user->yi[i-1]);CHKERRQ(ierr);
   }
-  ierr = Gather_i(C,user->yiwork,user->yi_scatter,user->nt);CHKERRQ(ierr);
-  ierr = VecAXPY(C,-1.0,user->q);CHKERRQ(ierr);
+  ierr = Gather(CE,user->yiwork,user->yi_scatter,user->nt);CHKERRQ(ierr);
+  ierr = VecAXPY(CE,-1.0,user->q);CHKERRQ(ierr);
+
   PetscFunctionReturn(0);
 }
 
 
 #undef __FUNCT__
 #define __FUNCT__ "Scatter"
-PetscErrorCode Scatter(Vec x, Vec state, VecScatter s_scat, Vec design, VecScatter d_scat)
-{
-  PetscErrorCode ierr;
-
-  PetscFunctionBegin;
-  ierr = VecScatterBegin(s_scat,x,state,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
-  ierr = VecScatterEnd(s_scat,x,state,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
-  ierr = VecScatterBegin(d_scat,x,design,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
-  ierr = VecScatterEnd(d_scat,x,design,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
-  PetscFunctionReturn(0);
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "Scatter_i"
-PetscErrorCode Scatter_i(Vec y, Vec *yi, VecScatter *scat, PetscInt nt)
+PetscErrorCode Scatter(Vec y, Vec *yi, VecScatter *scat, PetscInt nt)
 {
   PetscErrorCode ierr;
   PetscInt       i;
   PetscFunctionReturn(0);
 }
 
-
 #undef __FUNCT__
 #define __FUNCT__ "Gather"
-PetscErrorCode Gather(Vec x, Vec state, VecScatter s_scat, Vec design, VecScatter d_scat)
-{
-  PetscErrorCode ierr;
-
-  PetscFunctionBegin;
-  ierr = VecScatterBegin(s_scat,state,x,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
-  ierr = VecScatterEnd(s_scat,state,x,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
-  ierr = VecScatterBegin(d_scat,design,x,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
-  ierr = VecScatterEnd(d_scat,design,x,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
-  PetscFunctionReturn(0);
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "Gather_i"
-PetscErrorCode Gather_i(Vec y, Vec *yi, VecScatter *scat, PetscInt nt)
+PetscErrorCode Gather(Vec y, Vec *yi, VecScatter *scat, PetscInt nt)
 {
   PetscErrorCode ierr;
   PetscInt       i;
   for (i=1; i<user->nt; i++){
     ierr = VecSet(user->yiwork[i],0.0);CHKERRQ(ierr);
   }
-  ierr = Gather_i(user->q,user->yiwork,user->yi_scatter,user->nt);CHKERRQ(ierr);
+  ierr = Gather(user->q,user->yiwork,user->yi_scatter,user->nt);CHKERRQ(ierr);
   ierr = VecDestroy(&bc);CHKERRQ(ierr);
 
   /* Compute true state function yt given ut */
   /* Add noise to the measurement data */
   ierr = VecSet(user->ywork,1.0);CHKERRQ(ierr);
   ierr = VecAYPX(user->ywork,user->noise,user->ytrue);CHKERRQ(ierr);
-  ierr = Scatter_i(user->ywork,user->yiwork,user->yi_scatter,user->nt);CHKERRQ(ierr);
+  ierr = Scatter(user->ywork,user->yiwork,user->yi_scatter,user->nt);CHKERRQ(ierr);
   for (j=0; j<user->ns; j++){
     i = user->sample_times[j];
     ierr = MatMult(user->Qblock,user->yiwork[i],user->di[j]);
   }
-  ierr = Gather_i(user->d,user->di,user->di_scatter,user->ns);CHKERRQ(ierr);
+  ierr = Gather(user->d,user->di,user->di_scatter,user->ns);CHKERRQ(ierr);
 
   /* Now that initial conditions have been set, let the user pass tolerance options to the KSP solver */
   ierr = KSPSetFromOptions(user->solver);CHKERRQ(ierr);
   ierr = VecDestroy(&user->Twork);CHKERRQ(ierr);
   ierr = VecDestroy(&user->Rwork);CHKERRQ(ierr);
   ierr = VecDestroy(&user->js_diag);CHKERRQ(ierr);
-  ierr = ISDestroy(&user->s_is);CHKERRQ(ierr);
-  ierr = ISDestroy(&user->d_is);CHKERRQ(ierr);
-  ierr = VecScatterDestroy(&user->state_scatter);CHKERRQ(ierr);
-  ierr = VecScatterDestroy(&user->design_scatter);CHKERRQ(ierr);
+
   for (i=0; i<user->nt; i++){
     ierr = VecScatterDestroy(&user->yi_scatter[i]);CHKERRQ(ierr);
   }
 PetscErrorCode ParabolicMonitor(Tao tao, void *ptr)
 {
   PetscErrorCode ierr;
-  Vec            X;
+  Vec            Y,U;
   PetscReal      unorm,ynorm;
   AppCtx         *user = (AppCtx*)ptr;
 
   PetscFunctionBegin;
-  ierr = TaoGetSolutionVector(tao,&X);CHKERRQ(ierr);
-  ierr = Scatter(X,user->ywork,user->state_scatter,user->uwork,user->design_scatter);CHKERRQ(ierr);
-  ierr = VecAXPY(user->ywork,-1.0,user->ytrue);CHKERRQ(ierr);
-  ierr = VecAXPY(user->uwork,-1.0,user->utrue);CHKERRQ(ierr);
+
+  ierr = TaoGetOptControlSolutionVector(tao,&Y,&U);CHKERRQ(ierr);
+  ierr = VecWAXPY(user->ywork,-1.0,user->ytrue,Y);CHKERRQ(ierr);
+  ierr = VecWAXPY(user->uwork,-1.0,user->utrue,U);CHKERRQ(ierr);
   ierr = VecNorm(user->uwork,NORM_2,&unorm);CHKERRQ(ierr);
   ierr = VecNorm(user->ywork,NORM_2,&ynorm);CHKERRQ(ierr);
   ierr = PetscPrintf(MPI_COMM_WORLD, "||u-ut||=%g ||y-yt||=%g\n",(double)unorm,(double)ynorm);CHKERRQ(ierr);
+
   PetscFunctionReturn(0);
 }