Actual source code: parabolic.c

  1: #include "private/taosolver_impl.h"
  2: #include "src/pde_constrained/impls/lcl/lcl.h"
  3: #include "petsctime.h"

  5: /*T
  6:    Concepts: TAO - Solving a system of nonlinear equations, nonlinear ;east squares
  7:    Routines: TaoInitialize(); TaoFinalize(); 
  8:    Routines: TaoCreate();
  9:    Routines: TaoSetType(); 
 10:    Routines: TaoSetInitialVector();
 11:    Routines: TaoSetObjectiveRoutine();
 12:    Routines: TaoSetGradientRoutine();
 13:    Routines: TaoSetConstraintsRoutine();
 14:    Routines: TaoSetJacobianStateRoutine();
 15:    Routines: TaoSetJacobianDesignRoutine();
 16:    Routines: TaoSetStateDesignIS();
 17:    Routines: TaoSetFromOptions();
 18:    Routines: TaoSetHistory(); TaoGetHistory();
 19:    Routines: TaoSolve();
 20:    Routines: TaoGetTerminationReason(); TaoDestroy(); 
 21:    Processors: n
 22: T*/

 24: typedef struct {
 25:   PetscInt n; /*  Number of variables */
 26:   PetscInt m; /*  Number of constraints per time step */
 27:   PetscInt mx; /*  grid points in each direction */
 28:   PetscInt nt; /*  Number of time steps; as of now, must be divisible by 8 */
 29:   PetscInt ndata; /*  Number of data points per sample */
 30:   PetscInt ns; /*  Number of samples */
 31:   PetscInt *sample_times; /*  Times of samples */
 32:   IS s_is;
 33:   IS d_is;
 34:   VecScatter state_scatter;
 35:   VecScatter design_scatter;
 36:   VecScatter *yi_scatter;
 37:   VecScatter *di_scatter;

 39:   Mat Js,Jd,JsBlockPrec,JsInv,JsBlock;
 40:   PetscBool jformed,dsg_formed;

 42:   PetscReal alpha; /*  Regularization parameter */
 43:   PetscReal beta; /*  Weight attributed to ||u||^2 in regularization functional */
 44:   PetscReal noise; /*  Amount of noise to add to data */
 45:   PetscReal ht; /*  Time step */

 47:   Mat Qblock,QblockT;
 48:   Mat L,LT;
 49:   Mat Div,Divwork;
 50:   Mat Grad;
 51:   Mat Av,Avwork,AvT;
 52:   Mat DSG;
 53:   Vec q;
 54:   Vec ur; /*  reference */

 56:   Vec d;
 57:   Vec dwork;
 58:   Vec *di;

 60:   Vec y; /*  state variables */
 61:   Vec ywork;

 63:   Vec ytrue;
 64:   Vec *yi,*yiwork;

 66:   Vec u; /*  design variables */
 67:   Vec uwork;

 69:   Vec utrue;
 70:  
 71:   Vec js_diag;
 72:   
 73:   Vec c; /*  constraint vector */
 74:   Vec cwork;
 75:   
 76:   Vec lwork;
 77:   Vec S;
 78:   Vec Rwork,Swork,Twork;
 79:   Vec Av_u;

 81:   KSP solver;
 82:   PC prec;

 84:   TAO_LCL *lcl;
 85:   PetscReal tau[4];
 86:   PetscInt ksp_its;
 87:   PetscInt ksp_its_initial;

 89: } AppCtx;


 92: PetscErrorCode FormFunction(TaoSolver, Vec, PetscReal*, void*);
 93: PetscErrorCode FormGradient(TaoSolver, Vec, Vec, void*);
 94: PetscErrorCode FormFunctionGradient(TaoSolver, Vec, PetscReal*, Vec, void*);
 95: PetscErrorCode FormJacobianState(TaoSolver, Vec, Mat*, Mat*, Mat*, MatStructure*,void*);
 96: PetscErrorCode FormJacobianDesign(TaoSolver, Vec, Mat*, void*);
 97: PetscErrorCode FormConstraints(TaoSolver, Vec, Vec, void*);
 98: PetscErrorCode FormHessian(TaoSolver, Vec, Mat*, Mat*, MatStructure*, void*);
 99: PetscErrorCode Gather(Vec x, Vec state, VecScatter s_scat, Vec design, VecScatter d_scat);
100: PetscErrorCode Scatter(Vec x, Vec state, VecScatter s_scat, Vec design, VecScatter d_scat);
101: PetscErrorCode ParabolicInitialize(AppCtx *user);
102: PetscErrorCode ParabolicDestroy(AppCtx *user);
103: PetscErrorCode ParabolicMonitor(TaoSolver, void*);

105: PetscErrorCode StateMatMult(Mat,Vec,Vec);
106: PetscErrorCode StateMatBlockMult(Mat,Vec,Vec);
107: PetscErrorCode StateMatMultTranspose(Mat,Vec,Vec);
108: PetscErrorCode StateMatGetDiagonal(Mat,Vec); 
109: PetscErrorCode StateMatDuplicate(Mat,MatDuplicateOption,Mat*);
110: PetscErrorCode StateMatInvMult(Mat,Vec,Vec);
111: PetscErrorCode StateMatInvTransposeMult(Mat,Vec,Vec);
112: PetscErrorCode StateMatBlockPrecMult(PC,Vec,Vec);

114: PetscErrorCode DesignMatMult(Mat,Vec,Vec);
115: PetscErrorCode DesignMatMultTranspose(Mat,Vec,Vec);

117: PetscErrorCode Gather_i(Vec,Vec*,VecScatter*,PetscInt);
118: PetscErrorCode Scatter_i(Vec,Vec*,VecScatter*,PetscInt);

120: static  char help[]="";

124: int main(int argc, char **argv)
125: {
127:   Vec x,x0;
128:   
129:   TaoSolver tao;
130:   TaoSolverTerminationReason reason;
131:   AppCtx user;
132:   IS is_allstate,is_alldesign;
133:   PetscInt lo,hi,hi2,lo2,ksp_old;
134:   PetscBool flag,showtime;
135:   PetscInt ntests = 1;
136:   PetscLogDouble v1,v2;
137:   PetscInt i;
138:   int stages[1];
139:   

141:   PetscInitialize(&argc, &argv, (char*)0,help);
142:   TaoInitialize(&argc, &argv, (char*)0,help);

144:   showtime = PETSC_FALSE;
145:   PetscOptionsBool("-showtime","Display time elapsed","",showtime,&showtime,&flag); 
146:   user.mx = 8;
147:   PetscOptionsInt("-mx","Number of grid points in each direction","",user.mx,&user.mx,&flag); 
148:   user.nt = 8;
149:   PetscOptionsInt("-nt","Number of time steps","",user.nt,&user.nt,&flag); 
150:   user.ndata = 64;
151:   PetscOptionsInt("-ndata","Numbers of data points per sample","",user.ndata,&user.ndata,&flag); 
152:   user.ns = 8;
153:   PetscOptionsInt("-ns","Number of samples","",user.ns,&user.ns,&flag); 
154:   user.alpha = 1.0;
155:   PetscOptionsReal("-alpha","Regularization parameter","",user.alpha,&user.alpha,&flag); 
156:   user.beta = 0.01;
157:   PetscOptionsReal("-beta","Weight attributed to ||u||^2 in regularization functional","",user.beta,&user.beta,&flag); 
158:   user.noise = 0.01;
159:   PetscOptionsReal("-noise","Amount of noise to add to data","",user.noise,&user.noise,&flag); 

161:   user.tau[0] = user.tau[1] = user.tau[2] = user.tau[3] = 1.0e-4;
162:   PetscOptionsReal("-tola","Tolerance for first forward solve","",user.tau[0],&user.tau[0],&flag); 
163:   PetscOptionsReal("-tolb","Tolerance for first adjoint solve","",user.tau[1],&user.tau[1],&flag); 
164:   PetscOptionsReal("-tolc","Tolerance for second forward solve","",user.tau[2],&user.tau[2],&flag); 
165:   PetscOptionsReal("-told","Tolerance for second adjoint solve","",user.tau[3],&user.tau[3],&flag); 

167:   user.m = user.mx*user.mx*user.mx; /*  number of constraints per time step */
168:   user.n = user.m*(user.nt+1); /*  number of variables */
169:   user.ht = (PetscReal)1/user.nt;

171:   VecCreate(PETSC_COMM_WORLD,&user.u); 
172:   VecCreate(PETSC_COMM_WORLD,&user.y); 
173:   VecCreate(PETSC_COMM_WORLD,&user.c); 
174:   VecSetSizes(user.u,PETSC_DECIDE,user.n-user.m*user.nt); 
175:   VecSetSizes(user.y,PETSC_DECIDE,user.m*user.nt); 
176:   VecSetSizes(user.c,PETSC_DECIDE,user.m*user.nt); 
177:   VecSetFromOptions(user.u); 
178:   VecSetFromOptions(user.y); 
179:   VecSetFromOptions(user.c); 

181:   /* Create scatters for reduced spaces.
182:      If the state vector y and design vector u are partitioned as 
183:      [y_1; y_2; ...; y_np] and [u_1; u_2; ...; u_np] (with np = # of processors),
184:      then the solution vector x is organized as
185:      [y_1; u_1; y_2; u_2; ...; y_np; u_np]. 
186:      The index sets user.s_is and user.d_is correspond to the indices of the
187:      state and design variables owned by the current processor.
188:   */
189:   VecCreate(PETSC_COMM_WORLD,&x); 

191:   VecGetOwnershipRange(user.y,&lo,&hi); 
192:   VecGetOwnershipRange(user.u,&lo2,&hi2);  

194:   ISCreateStride(PETSC_COMM_SELF,hi-lo,lo,1,&is_allstate); 
195:   ISCreateStride(PETSC_COMM_SELF,hi-lo,lo+lo2,1,&user.s_is); 

197:   ISCreateStride(PETSC_COMM_SELF,hi2-lo2,lo2,1,&is_alldesign); 
198:   ISCreateStride(PETSC_COMM_SELF,hi2-lo2,hi+lo2,1,&user.d_is); 

200:   VecSetSizes(x,hi-lo+hi2-lo2,user.n); 
201:   VecSetFromOptions(x); 

203:   VecScatterCreate(x,user.s_is,user.y,is_allstate,&user.state_scatter); 
204:   VecScatterCreate(x,user.d_is,user.u,is_alldesign,&user.design_scatter); 
205:   ISDestroy(&is_alldesign); 
206:   ISDestroy(&is_allstate); 



210:   /* Create TAO solver and set desired solution method */
211:   TaoCreate(PETSC_COMM_WORLD,&tao); 
212:   TaoSetType(tao,"tao_lcl"); 
213:   user.lcl = (TAO_LCL*)(tao->data); 
214:   /* Set up initial vectors and matrices */
215:   ParabolicInitialize(&user); 

217:   Gather(x,user.y,user.state_scatter,user.u,user.design_scatter); 
218:   VecDuplicate(x,&x0); 
219:   VecCopy(x,x0); 

221:   /* Set solution vector with an initial guess */
222:   TaoSetInitialVector(tao,x); 
223:   TaoSetObjectiveRoutine(tao, FormFunction, (void *)&user); 
224:   TaoSetGradientRoutine(tao, FormGradient, (void *)&user); 
225:   TaoSetConstraintsRoutine(tao, user.c, FormConstraints, (void *)&user); 

227:   TaoSetJacobianStateRoutine(tao, user.Js, user.JsBlockPrec, user.JsInv, FormJacobianState, (void *)&user);  

229:   TaoSetJacobianDesignRoutine(tao, user.Jd, FormJacobianDesign, (void *)&user); 

231:   TaoSetFromOptions(tao); 
232:   TaoSetStateDesignIS(tao,user.s_is,user.d_is); 

234:  /* SOLVE THE APPLICATION */
235:   PetscOptionsInt("-ntests","Number of times to repeat TaoSolve","",ntests,&ntests,&flag); 
236:   PetscLogStageRegister("Trials",&stages[0]); 
237:   PetscLogStagePush(stages[0]); 
238:   user.ksp_its_initial = user.ksp_its;
239:   ksp_old = user.ksp_its;
240:   for (i=0; i<ntests; i++){
241:     PetscGetTime(&v1); 
242:     TaoSolve(tao);  
243:     PetscGetTime(&v2); 
244:     if (showtime) {
245:       PetscPrintf(PETSC_COMM_WORLD,"Elapsed time = %G\n",v2-v1); 
246:     }
247:     PetscPrintf(PETSC_COMM_WORLD,"KSP Iterations = %D\n",user.ksp_its-ksp_old); 
248:     VecCopy(x0,x); 
249:     TaoSetInitialVector(tao,x); 
250:   }
251:   PetscLogStagePop(); 
252:   PetscBarrier((PetscObject)x); 
253:   PetscPrintf(PETSC_COMM_WORLD,"KSP iterations within initialization: ");
254:   PetscPrintf(PETSC_COMM_WORLD,"%D\n",user.ksp_its_initial);
255:   PetscPrintf(PETSC_COMM_WORLD,"Total KSP iterations over %D trial(s): ",ntests);
256:   PetscPrintf(PETSC_COMM_WORLD,"%D\n",user.ksp_its);
257:   PetscPrintf(PETSC_COMM_WORLD,"KSP iterations per trial: ");
258:   PetscPrintf(PETSC_COMM_WORLD,"%D\n",(user.ksp_its-user.ksp_its_initial)/ntests);

260:   TaoGetTerminationReason(tao,&reason); 

262:   if (reason < 0)
263:   {
264:     PetscPrintf(MPI_COMM_WORLD, "TAO failed to converge.\n");
265:   }
266:   else
267:   {
268:     PetscPrintf(MPI_COMM_WORLD, "Optimization terminated with status %D.\n", reason);
269:   }


272:   /* Free TAO data structures */
273:   TaoDestroy(&tao); 

275:   /* Free PETSc data structures */
276:   VecDestroy(&x); 
277:   VecDestroy(&x0); 
278:   ParabolicDestroy(&user); 

280:   /* Finalize TAO, PETSc */
281:   TaoFinalize();
282:   PetscFinalize();

284:   return 0;
285: }
286: /* ------------------------------------------------------------------- */
289: /* 
290:    dwork = Qy - d  
291:    lwork = L*(u-ur)
292:    f = 1/2 * (dwork.dork + alpha*lwork.lwork)
293: */
294: PetscErrorCode FormFunction(TaoSolver tao,Vec X,PetscReal *f,void *ptr)
295: {
297:   PetscReal d1=0,d2=0;
298:   PetscInt i,j;
299:   AppCtx *user = (AppCtx*)ptr;
301:   Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter); 
302:   Scatter_i(user->y,user->yi,user->yi_scatter,user->nt); 
303:   for (j=0; j<user->ns; j++){
304:     i = user->sample_times[j];
305:     MatMult(user->Qblock,user->yi[i],user->di[j]); 
306:   }
307:   Gather_i(user->dwork,user->di,user->di_scatter,user->ns); 
308:   VecAXPY(user->dwork,-1.0,user->d); 
309:   VecDot(user->dwork,user->dwork,&d1); 

311:   VecWAXPY(user->uwork,-1.0,user->ur,user->u); 
312:   MatMult(user->L,user->uwork,user->lwork); 
313:   VecDot(user->lwork,user->lwork,&d2); 

315:   *f = 0.5 * (d1 + user->alpha*d2);
316:   return(0);
317: }

319: /* ------------------------------------------------------------------- */
322: /*  
323:     state: g_s = Q' *(Qy - d)
324:     design: g_d = alpha*L'*L*(u-ur)
325: */
326: PetscErrorCode FormGradient(TaoSolver tao,Vec X,Vec G,void *ptr)
327: {
329:   PetscInt i,j;
330:   AppCtx *user = (AppCtx*)ptr;
332:   CHKMEMQ;
333:   Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter); 
334:   Scatter_i(user->y,user->yi,user->yi_scatter,user->nt); 
335:   for (j=0; j<user->ns; j++){
336:     i = user->sample_times[j];
337:     MatMult(user->Qblock,user->yi[i],user->di[j]); 
338:   }
339:   Gather_i(user->dwork,user->di,user->di_scatter,user->ns); 
340:   VecAXPY(user->dwork,-1.0,user->d); 
341:   Scatter_i(user->dwork,user->di,user->di_scatter,user->ns); 
342:   VecSet(user->ywork,0.0); 
343:   Scatter_i(user->ywork,user->yiwork,user->yi_scatter,user->nt); 
344:   for (j=0; j<user->ns; j++){
345:     i = user->sample_times[j];
346:     MatMult(user->QblockT,user->di[j],user->yiwork[i]); 
347:   }
348:   Gather_i(user->ywork,user->yiwork,user->yi_scatter,user->nt); 
349:   
350:   VecWAXPY(user->uwork,-1.0,user->ur,user->u); 
351:   MatMult(user->L,user->uwork,user->lwork); 
352:   MatMult(user->LT,user->lwork,user->uwork); 
353:   VecScale(user->uwork, user->alpha); 

355:                       
356:   Gather(G,user->ywork,user->state_scatter,user->uwork,user->design_scatter); 
357:   CHKMEMQ;
358:   return(0);
359: }

363: PetscErrorCode FormFunctionGradient(TaoSolver tao, Vec X, PetscReal *f, Vec G, void *ptr)
364: {
366:   PetscReal d1,d2;
367:   PetscInt i,j;
368:   AppCtx *user = (AppCtx*)ptr;
370:   Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter); 

372:   Scatter_i(user->y,user->yi,user->yi_scatter,user->nt); 
373:   for (j=0; j<user->ns; j++){
374:     i = user->sample_times[j];
375:     MatMult(user->Qblock,user->yi[i],user->di[j]); 
376:   }
377:   Gather_i(user->dwork,user->di,user->di_scatter,user->ns); 
378:   VecAXPY(user->dwork,-1.0,user->d); 
379:   VecDot(user->dwork,user->dwork,&d1); 
380:   Scatter_i(user->dwork,user->di,user->di_scatter,user->ns); 
381:   VecSet(user->ywork,0.0); 
382:   Scatter_i(user->ywork,user->yiwork,user->yi_scatter,user->nt); 
383:   for (j=0; j<user->ns; j++){
384:     i = user->sample_times[j];
385:     MatMult(user->QblockT,user->di[j],user->yiwork[i]); 
386:   }
387:   Gather_i(user->ywork,user->yiwork,user->yi_scatter,user->nt); 

389:   VecWAXPY(user->uwork,-1.0,user->ur,user->u); 
390:   MatMult(user->L,user->uwork,user->lwork); 
391:   VecDot(user->lwork,user->lwork,&d2); 
392:   MatMult(user->LT,user->lwork,user->uwork); 
393:   VecScale(user->uwork, user->alpha); 
394:   *f = 0.5 * (d1 + user->alpha*d2); 
395:   
396:   Gather(G,user->ywork,user->state_scatter,user->uwork,user->design_scatter); 
397:   return(0);

399: }


402: /* ------------------------------------------------------------------- */
405: /* A 
406: MatShell object
407: */
408: PetscErrorCode FormJacobianState(TaoSolver tao, Vec X, Mat *J, Mat* JPre, Mat* JInv, MatStructure* flag, void *ptr)
409: {
411:   AppCtx *user = (AppCtx*)ptr;
413:   Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter); 
414:   VecSet(user->uwork,0); 
415:   VecAXPY(user->uwork,-1.0,user->u); 
416:   VecExp(user->uwork); 
417:   MatMult(user->Av,user->uwork,user->Av_u); 
418:   VecCopy(user->Av_u,user->Swork);  
419:   VecReciprocal(user->Swork); 
420:   MatCopy(user->Div,user->Divwork,SAME_NONZERO_PATTERN); 
421:   MatDiagonalScale(user->Divwork,PETSC_NULL,user->Swork); 
422:   if (user->dsg_formed) {
423:     MatMatMultNumeric(user->Divwork,user->Grad,user->DSG); 
424:   } else {
425:     MatMatMultSymbolic(user->Divwork,user->Grad,PETSC_DECIDE,&user->DSG); 
426:     MatMatMultNumeric(user->Divwork,user->Grad,user->DSG); 
427:     user->dsg_formed = PETSC_TRUE;
428:   }
429:   
430:   /* B = speye(nx^3) + ht*DSG; */
431:   MatScale(user->DSG,user->ht); 
432:   MatShift(user->DSG,1.0); 
433:     
434:   *JInv = user->JsInv;

436:   return(0);
437: }
438: /* ------------------------------------------------------------------- */
441: /* B */
442: PetscErrorCode FormJacobianDesign(TaoSolver tao, Vec X, Mat *J, void *ptr)
443: {
445:   AppCtx *user = (AppCtx*)ptr;

448:   Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter); 

450:   *J = user->Jd;

452:   return(0);
453: }



459: PetscErrorCode StateMatMult(Mat J_shell, Vec X, Vec Y) 
460: {
462:   PetscInt i;
463:   void *ptr;
464:   AppCtx *user;
466:   MatShellGetContext(J_shell,&ptr); 
467:   user = (AppCtx*)ptr;

469:   Scatter_i(X,user->yi,user->yi_scatter,user->nt); 
470:   
471:   MatMult(user->JsBlock,user->yi[0],user->yiwork[0]); 

473:   for (i=1; i<user->nt; i++){
474:     MatMult(user->JsBlock,user->yi[i],user->yiwork[i]); 
475:     VecAXPY(user->yiwork[i],-1.0,user->yi[i-1]); 
476:   }

478:   Gather_i(Y,user->yiwork,user->yi_scatter,user->nt); 

480:   return(0);
481: }

485: PetscErrorCode StateMatMultTranspose(Mat J_shell, Vec X, Vec Y) 
486: {
488:   PetscInt i;
489:   void *ptr;
490:   AppCtx *user;
492:   MatShellGetContext(J_shell,&ptr); 
493:   user = (AppCtx*)ptr;

495:   Scatter_i(X,user->yi,user->yi_scatter,user->nt); 

497:   for (i=0; i<user->nt-1; i++){
498:     MatMult(user->JsBlock,user->yi[i],user->yiwork[i]); 
499:     VecAXPY(user->yiwork[i],-1.0,user->yi[i+1]); 
500:   }

502:   i = user->nt-1;
503:   MatMult(user->JsBlock,user->yi[i],user->yiwork[i]); 

505:   Gather_i(Y,user->yiwork,user->yi_scatter,user->nt); 

507:   return(0);
508: }

512: PetscErrorCode StateMatBlockMult(Mat J_shell, Vec X, Vec Y) 
513: {
515:   void *ptr;
516:   AppCtx *user;
518:   MatShellGetContext(J_shell,&ptr); 
519:   user = (AppCtx*)ptr;
520:    
521:   MatMult(user->Grad,X,user->Swork);  
522:   VecPointwiseDivide(user->Swork,user->Swork,user->Av_u); 
523:   MatMult(user->Div,user->Swork,Y);  
524:   VecAYPX(Y,user->ht,X); 

526:   return(0);
527: }

531: PetscErrorCode DesignMatMult(Mat J_shell, Vec X, Vec Y) 
532: {
534:   void *ptr;
535:   PetscInt i;
536:   AppCtx *user;
538:   MatShellGetContext(J_shell,&ptr); 
539:   user = (AppCtx*)ptr;
540:  
541:   /* sdiag(1./v) */ 
542:   VecSet(user->uwork,0); 
543:   VecAXPY(user->uwork,-1.0,user->u); 
544:   VecExp(user->uwork);   

546:   /* sdiag(1./((Av*(1./v)).^2)) */
547:   MatMult(user->Av,user->uwork,user->Swork); 
548:   VecPointwiseMult(user->Swork,user->Swork,user->Swork); 
549:   VecReciprocal(user->Swork);  

551:   /* (Av * (sdiag(1./v) * b)) */ 
552:   VecPointwiseMult(user->uwork,user->uwork,X); 
553:   MatMult(user->Av,user->uwork,user->Twork); 

555:   /* (sdiag(1./((Av*(1./v)).^2)) * (Av * (sdiag(1./v) * b))) */
556:   VecPointwiseMult(user->Swork,user->Twork,user->Swork);  

558:   Scatter_i(user->y,user->yi,user->yi_scatter,user->nt); 
559:   for (i=0; i<user->nt; i++){
560:     /* (sdiag(Grad*y(:,i)) */
561:     MatMult(user->Grad,user->yi[i],user->Twork); 
562:   
563:     /* ht * Div * (sdiag(Grad*y(:,i)) * (sdiag(1./((Av*(1./v)).^2)) * (Av * (sdiag(1./v) * b)))) */
564:     VecPointwiseMult(user->Twork,user->Twork,user->Swork);  
565:     MatMult(user->Div,user->Twork,user->yiwork[i]); 
566:     VecScale(user->yiwork[i],user->ht); 
567:   }
568:   Gather_i(Y,user->yiwork,user->yi_scatter,user->nt); 

570:   return(0);
571: }

575: PetscErrorCode DesignMatMultTranspose(Mat J_shell, Vec X, Vec Y) 
576: {
578:   void *ptr;
579:   PetscInt i;
580:   AppCtx *user;
582:   MatShellGetContext(J_shell,&ptr); 
583:   user = (AppCtx*)ptr;

585:   /* sdiag(1./((Av*(1./v)).^2)) */
586:   VecSet(user->uwork,0); 
587:   VecAXPY(user->uwork,-1.0,user->u); 
588:   VecExp(user->uwork); 
589:   MatMult(user->Av,user->uwork,user->Rwork); 
590:   VecPointwiseMult(user->Rwork,user->Rwork,user->Rwork); 
591:   VecReciprocal(user->Rwork); 

593:   VecSet(Y,0.0); 
594:   Scatter_i(user->y,user->yi,user->yi_scatter,user->nt); 
595:   Scatter_i(X,user->yiwork,user->yi_scatter,user->nt); 
596:   for (i=0; i<user->nt; i++){
597:     /* (Div' * b(:,i)) */
598:     MatMult(user->Grad,user->yiwork[i],user->Swork); 

600:     /* sdiag(Grad*y(:,i)) */
601:     MatMult(user->Grad,user->yi[i],user->Twork); 

603:     /* (sdiag(Grad*y(:,i)) * (Div' * b(:,i))) */
604:     VecPointwiseMult(user->Twork,user->Swork,user->Twork); 

606:     /* (sdiag(1./((Av*(1./v)).^2)) * (sdiag(Grad*y(:,i)) * (Div' * b(:,i)))) */
607:     VecPointwiseMult(user->Twork,user->Rwork,user->Twork); 

609:     /* (Av' * (sdiag(1./((Av*(1./v)).^2)) * (sdiag(Grad*y(:,i)) * (Div' * b(:,i))))) */
610:     MatMult(user->AvT,user->Twork,user->yiwork[i]); 
611:   
612:     /* sdiag(1./v) * (Av' * (sdiag(1./((Av*(1./v)).^2)) * (sdiag(Grad*y(:,i)) * (Div' * b(:,i))))) */
613:     VecPointwiseMult(user->yiwork[i],user->uwork,user->yiwork[i]); 
614:     VecAXPY(Y,user->ht,user->yiwork[i]); 
615:   }

617:   return(0);
618: }

622: PetscErrorCode StateMatBlockPrecMult(PC PC_shell, Vec X, Vec Y) 
623: {
625:   void *ptr;
626:   AppCtx *user;
628:   PCShellGetContext(PC_shell,&ptr); 
629:   user = (AppCtx*)ptr;

631:   if (user->dsg_formed) {
632:     MatSOR(user->DSG,X,1.0,(MatSORType)(SOR_ZERO_INITIAL_GUESS | SOR_LOCAL_SYMMETRIC_SWEEP),0.0,1,1,Y); 
633:   }
634:   else {
635:     printf("DSG not formed"); abort();
636:   }

638:   return(0);
639: }

643: PetscErrorCode StateMatInvMult(Mat J_shell, Vec X, Vec Y)
644: {
646:   void *ptr;
647:   AppCtx *user;
648:   PetscReal tau;
649:   PetscInt its,i;
651:   MatShellGetContext(J_shell,&ptr); 
652:   user = (AppCtx*)ptr;

654:   if (Y == user->ytrue) {
655:     KSPSetTolerances(user->solver,1e-8,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
656:   } else if (user->lcl) {
657:     tau = user->tau[user->lcl->solve_type];
658:     KSPSetTolerances(user->solver,tau,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
659:   }

661:   Scatter_i(X,user->yi,user->yi_scatter,user->nt); 
662:   
663:   KSPSolve(user->solver,user->yi[0],user->yiwork[0]); 

665:   KSPGetIterationNumber(user->solver,&its); 
666:   user->ksp_its = user->ksp_its + its;


669:   for (i=1; i<user->nt; i++){
670:     VecAXPY(user->yi[i],1.0,user->yiwork[i-1]); 
671:     KSPSolve(user->solver,user->yi[i],user->yiwork[i]); 

673:     KSPGetIterationNumber(user->solver,&its); 
674:     user->ksp_its = user->ksp_its + its;

676:   }

678:   Gather_i(Y,user->yiwork,user->yi_scatter,user->nt); 


681:   return(0);
682: }

686: PetscErrorCode StateMatInvTransposeMult(Mat J_shell, Vec X, Vec Y)
687: {
689:   void *ptr;
690:   PetscReal tau;
691:   AppCtx *user;
692:   PetscInt its,i;

695:   MatShellGetContext(J_shell,&ptr); 
696:   user = (AppCtx*)ptr;
697:   if (user->lcl) {
698:     tau = user->tau[user->lcl->solve_type];
699:     KSPSetTolerances(user->solver,tau,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
700:   }


703:   Scatter_i(X,user->yi,user->yi_scatter,user->nt); 
704:   
705:   i = user->nt - 1;
706:   KSPSolve(user->solver,user->yi[i],user->yiwork[i]); 
707:  
708:   KSPGetIterationNumber(user->solver,&its); 
709:   user->ksp_its = user->ksp_its + its;

711:   for (i=user->nt-2; i>=0; i--){
712:     VecAXPY(user->yi[i],1.0,user->yiwork[i+1]); 
713:     KSPSolve(user->solver,user->yi[i],user->yiwork[i]); 

715:     KSPGetIterationNumber(user->solver,&its); 
716:     user->ksp_its = user->ksp_its + its;

718:   }

720:   Gather_i(Y,user->yiwork,user->yi_scatter,user->nt); 
721:   return(0);
722: }

726: PetscErrorCode StateMatDuplicate(Mat J_shell, MatDuplicateOption opt, Mat *new_shell)
727: {
729:   void *ptr;
730:   AppCtx *user;
732:   MatShellGetContext(J_shell,&ptr); 
733:   user = (AppCtx*)ptr;

735:   MatCreateShell(PETSC_COMM_WORLD,PETSC_DETERMINE,PETSC_DETERMINE,user->m,user->m,user,new_shell); 
736:   MatShellSetOperation(*new_shell,MATOP_MULT,(void(*)(void))StateMatMult); 
737:   MatShellSetOperation(*new_shell,MATOP_DUPLICATE,(void(*)(void))StateMatDuplicate); 
738:   MatShellSetOperation(*new_shell,MATOP_MULT_TRANSPOSE,(void(*)(void))StateMatMultTranspose); 
739:   MatShellSetOperation(*new_shell,MATOP_GET_DIAGONAL,(void(*)(void))StateMatGetDiagonal); 
740:   
741:   return(0);
742: }

746: PetscErrorCode StateMatGetDiagonal(Mat J_shell, Vec X)
747: {
749:   void *ptr;
750:   AppCtx *user;
752:   MatShellGetContext(J_shell,&ptr); 
753:   user = (AppCtx*)ptr;
754:    VecCopy(user->js_diag,X); 
755:   return(0);
756:   
757: }

761: PetscErrorCode FormConstraints(TaoSolver tao, Vec X, Vec C, void *ptr)
762: {
763:   /* con = Ay - q, A = [B  0  0 ... 0; 
764:                        -I  B  0 ... 0; 
765:                         0 -I  B ... 0;
766:                              ...     ;
767:                         0    ... -I B] 
768:      B = ht * Div * Sigma * Grad + eye */
770:   PetscInt i;
771:   AppCtx *user = (AppCtx*)ptr;
773:    
774:   Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter); 
775:   Scatter_i(user->y,user->yi,user->yi_scatter,user->nt); 
776:    
777:   MatMult(user->JsBlock,user->yi[0],user->yiwork[0]); 

779:   for (i=1; i<user->nt; i++){
780:     MatMult(user->JsBlock,user->yi[i],user->yiwork[i]); 
781:     VecAXPY(user->yiwork[i],-1.0,user->yi[i-1]);                     
782:   }

784:   Gather_i(C,user->yiwork,user->yi_scatter,user->nt); 
785:   VecAXPY(C,-1.0,user->q); 

787:   return(0);
788: }


793: PetscErrorCode Scatter(Vec x, Vec state, VecScatter s_scat, Vec design, VecScatter d_scat)
794: {
797:   VecScatterBegin(s_scat,x,state,INSERT_VALUES,SCATTER_FORWARD); 
798:   VecScatterEnd(s_scat,x,state,INSERT_VALUES,SCATTER_FORWARD); 

800:   VecScatterBegin(d_scat,x,design,INSERT_VALUES,SCATTER_FORWARD); 
801:   VecScatterEnd(d_scat,x,design,INSERT_VALUES,SCATTER_FORWARD); 
802:   return(0);
803: }

807: PetscErrorCode Scatter_i(Vec y, Vec *yi, VecScatter *scat, PetscInt nt)
808: {
810:   PetscInt i;
812:   for (i=0; i<nt; i++){
813:     VecScatterBegin(scat[i],y,yi[i],INSERT_VALUES,SCATTER_FORWARD); 
814:     VecScatterEnd(scat[i],y,yi[i],INSERT_VALUES,SCATTER_FORWARD); 
815:   }
816:   return(0);
817: }


822: PetscErrorCode Gather(Vec x, Vec state, VecScatter s_scat, Vec design, VecScatter d_scat)
823: {
826:   VecScatterBegin(s_scat,state,x,INSERT_VALUES,SCATTER_REVERSE); 
827:   VecScatterEnd(s_scat,state,x,INSERT_VALUES,SCATTER_REVERSE); 
828:   VecScatterBegin(d_scat,design,x,INSERT_VALUES,SCATTER_REVERSE); 
829:   VecScatterEnd(d_scat,design,x,INSERT_VALUES,SCATTER_REVERSE); 
830:   return(0);
831: }

835: PetscErrorCode Gather_i(Vec y, Vec *yi, VecScatter *scat, PetscInt nt)
836: {
838:   PetscInt i;
840:   for (i=0; i<nt; i++){
841:     VecScatterBegin(scat[i],yi[i],y,INSERT_VALUES,SCATTER_REVERSE); 
842:     VecScatterEnd(scat[i],yi[i],y,INSERT_VALUES,SCATTER_REVERSE); 
843:   }
844:   return(0);
845: }
846:   
847:     
850: PetscErrorCode ParabolicInitialize(AppCtx *user)
851: {
853:   PetscInt m,n,i,j,k,linear_index,istart,iend,iblock,lo,hi,lo2,hi2;
854:   Vec XX,YY,ZZ,XXwork,YYwork,ZZwork,UTwork,yi,di,bc;
855:   PetscReal *x, *y, *z;
856:   PetscReal h,stime;
857:   PetscScalar hinv,neg_hinv,half = 0.5,sqrt_beta;
858:   PetscInt im,indx1,indx2,indy1,indy2,indz1,indz2,nx,ny,nz;
859:   PetscReal xri,yri,zri,xim,yim,zim,dx1,dx2,dy1,dy2,dz1,dz2,Dx,Dy,Dz;
860:   PetscScalar v,vx,vy,vz;
861:   IS is_from_y,is_to_yi,is_from_d,is_to_di;
862:   /* Data locations */
863:   PetscScalar xr[64] = {0.4970,     0.8498,     0.7814,     0.6268,     0.7782,     0.6402,     0.3617,     0.3160,     
864:                         0.3610,     0.5298,     0.6987,     0.3331,     0.7962,     0.5596,     0.3866,     0.6774,     
865:                         0.5407,     0.4518,     0.6702,     0.6061,     0.7580,     0.8997,     0.5198,     0.8326,     
866:                         0.2138,     0.9198,     0.3000,     0.2833,     0.8288,     0.7076,     0.1820,     0.0728,     
867:                         0.8447,     0.2367,     0.3239,     0.6413,     0.3114,     0.4731,     0.1192,     0.9273,     
868:                         0.5724,     0.4331,     0.5136,     0.3547,     0.4413,     0.2602,     0.5698,     0.7278,     
869:                         0.5261,     0.6230,     0.2454,     0.3948,     0.7479,     0.6582,     0.4660,     0.5594,     
870:                         0.7574,     0.1143,     0.5900,     0.1065,     0.4260,     0.3294,     0.8276,     0.0756};
871:   
872:   PetscScalar yr[64] = {0.7345,     0.9120,     0.9288,     0.7528,     0.4463,     0.4985,     0.2497,     0.6256,     
873:                         0.3425,     0.9026,     0.6983,     0.4230,     0.7140,     0.2970,     0.4474,     0.8792,     
874:                         0.6604,     0.2485,     0.7968,     0.6127,     0.1796,     0.2437,     0.5938,     0.6137,     
875:                         0.3867,     0.5658,     0.4575,     0.1009,     0.0863,     0.3361,     0.0738,     0.3985,     
876:                         0.6602,     0.1437,     0.0934,     0.5983,     0.5950,     0.0763,     0.0768,     0.2288,     
877:                         0.5761,     0.1129,     0.3841,     0.6150,     0.6904,     0.6686,     0.1361,     0.4601,     
878:                         0.4491,     0.3716,     0.1969,     0.6537,     0.6743,     0.6991,     0.4811,     0.5480,     
879:                         0.1684,     0.4569,     0.6889,     0.8437,     0.3015,     0.2854,     0.8199,     0.2658};
880:   
881:   PetscScalar zr[64] = {0.7668,     0.8573,     0.2654,     0.2719,     0.1060,     0.1311,     0.6232,     0.2295,     
882:                         0.8009,     0.2147,     0.2119,     0.9325,     0.4473,     0.3600,     0.3374,     0.3819,     
883:                         0.4066,     0.5801,     0.1673,     0.0959,     0.4638,     0.8236,     0.8800,     0.2939,     
884:                         0.2028,     0.8262,     0.2706,     0.6276,     0.9085,     0.6443,     0.8241,     0.0712,     
885:                         0.1824,     0.7789,     0.4389,     0.8415,     0.7055,     0.6639,     0.3653,     0.2078,     
886:                         0.1987,     0.2297,     0.4321,     0.8115,     0.4915,     0.7764,     0.4657,     0.4627,     
887:                         0.4569,     0.4232,     0.8514,     0.0674,     0.3227,     0.1055,     0.6690,     0.6313,     
888:                         0.9226,     0.5461,     0.4126,     0.2364,     0.6096,     0.7042,     0.3914,     0.0711};

891:   PetscMalloc(user->mx*sizeof(PetscReal),&x); 
892:   PetscMalloc(user->mx*sizeof(PetscReal),&y); 
893:   PetscMalloc(user->mx*sizeof(PetscReal),&z); 
894:   user->jformed = PETSC_FALSE;
895:   user->dsg_formed = PETSC_FALSE;

897:   n = user->mx * user->mx * user->mx;
898:   m = 3 * user->mx * user->mx * (user->mx-1);
899:   sqrt_beta = PetscSqrtScalar(user->beta);

901:   user->ksp_its = 0;
902:   user->ksp_its_initial = 0;

904:   stime = (PetscReal)user->nt/user->ns;
905:   PetscMalloc(user->ns*sizeof(PetscInt),&user->sample_times); 
906:   for (i=0; i<user->ns; i++){
907:     /* round((i+1)*stime)-1 */
908:     user->sample_times[i] = (PetscInt)(stime*((PetscReal)i+1.0)-0.5);
909:   }

911:   VecCreate(PETSC_COMM_WORLD,&XX); 
912:   VecCreate(PETSC_COMM_WORLD,&user->q); 
913:   VecSetSizes(XX,PETSC_DECIDE,n); 
914:   VecSetSizes(user->q,PETSC_DECIDE,n*user->nt); 
915:   VecSetFromOptions(XX); 
916:   VecSetFromOptions(user->q); 

918:   VecDuplicate(XX,&YY); 
919:   VecDuplicate(XX,&ZZ); 
920:   VecDuplicate(XX,&XXwork); 
921:   VecDuplicate(XX,&YYwork); 
922:   VecDuplicate(XX,&ZZwork); 
923:   VecDuplicate(XX,&UTwork); 
924:   VecDuplicate(XX,&user->utrue); 
925:   VecDuplicate(XX,&bc); 

927:   /* Generate 3D grid, and collect ns (1<=ns<=8) right-hand-side vectors into user->q */
928:   h = 1.0/user->mx; 
929:   hinv = user->mx;
930:   neg_hinv = -hinv;
931:  
932:   VecGetOwnershipRange(XX,&istart,&iend); 
933:   for (linear_index=istart; linear_index<iend; linear_index++){
934:     i = linear_index % user->mx;
935:     j = ((linear_index-i)/user->mx) % user->mx;
936:     k = ((linear_index-i)/user->mx-j) / user->mx;
937:     vx = h*(i+0.5);
938:     vy = h*(j+0.5);
939:     vz = h*(k+0.5);        
940:     VecSetValues(XX,1,&linear_index,&vx,INSERT_VALUES); 
941:     VecSetValues(YY,1,&linear_index,&vy,INSERT_VALUES); 
942:     VecSetValues(ZZ,1,&linear_index,&vz,INSERT_VALUES); 
943:     if ((vx<0.6) && (vx>0.4) && (vy<0.6) && (vy>0.4) && (vy<0.6) && (vz<0.6) && (vz>0.4)){
944:       v = 1000.0;
945:       VecSetValues(bc,1,&linear_index,&v,INSERT_VALUES); 
946:     }
947:   }

949:   VecAssemblyBegin(XX); 
950:   VecAssemblyEnd(XX); 
951:   VecAssemblyBegin(YY); 
952:   VecAssemblyEnd(YY); 
953:   VecAssemblyBegin(ZZ); 
954:   VecAssemblyEnd(ZZ); 
955:   VecAssemblyBegin(bc); 
956:   VecAssemblyEnd(bc);   

958:   /* Compute true parameter function
959:      ut = 0.5 + exp(-10*((x-0.5)^2+(y-0.5)^2+(z-0.5)^2)) */
960:   VecCopy(XX,XXwork); 
961:   VecCopy(YY,YYwork); 
962:   VecCopy(ZZ,ZZwork); 

964:   VecShift(XXwork,-0.5); 
965:   VecShift(YYwork,-0.5); 
966:   VecShift(ZZwork,-0.5); 

968:   VecPointwiseMult(XXwork,XXwork,XXwork); 
969:   VecPointwiseMult(YYwork,YYwork,YYwork); 
970:   VecPointwiseMult(ZZwork,ZZwork,ZZwork); 

972:   VecCopy(XXwork,user->utrue); 
973:   VecAXPY(user->utrue,1.0,YYwork); 
974:   VecAXPY(user->utrue,1.0,ZZwork); 
975:   VecScale(user->utrue,-10.0); 
976:   VecExp(user->utrue); 
977:   VecShift(user->utrue,0.5); 

979:   VecDestroy(&XX); 
980:   VecDestroy(&YY); 
981:   VecDestroy(&ZZ); 
982:   VecDestroy(&XXwork); 
983:   VecDestroy(&YYwork); 
984:   VecDestroy(&ZZwork); 
985:   VecDestroy(&UTwork); 
986:  
987:   /* Initial guess and reference model */
988:   VecDuplicate(user->utrue,&user->ur); 
989:   VecSet(user->ur,0.0); 

991:   /* Generate Grad matrix */
992:   MatCreate(PETSC_COMM_WORLD,&user->Grad); 
993:   MatSetSizes(user->Grad,PETSC_DECIDE,PETSC_DECIDE,m,n); 
994:   MatSetFromOptions(user->Grad); 
995:   MatMPIAIJSetPreallocation(user->Grad,2,PETSC_NULL,2,PETSC_NULL); 
996:   MatSeqAIJSetPreallocation(user->Grad,2,PETSC_NULL); 
997:   MatGetOwnershipRange(user->Grad,&istart,&iend); 

999:   for (i=istart; i<iend; i++){
1000:     if (i<m/3){
1001:       iblock = i / (user->mx-1);
1002:       j = iblock*user->mx + (i % (user->mx-1));
1003:       MatSetValues(user->Grad,1,&i,1,&j,&neg_hinv,INSERT_VALUES); 
1004:       j = j+1;
1005:       MatSetValues(user->Grad,1,&i,1,&j,&hinv,INSERT_VALUES); 
1006:     }
1007:     if (i>=m/3 && i<2*m/3){
1008:       iblock = (i-m/3) / (user->mx*(user->mx-1));
1009:       j = iblock*user->mx*user->mx + ((i-m/3) % (user->mx*(user->mx-1)));
1010:       MatSetValues(user->Grad,1,&i,1,&j,&neg_hinv,INSERT_VALUES); 
1011:       j = j + user->mx;
1012:       MatSetValues(user->Grad,1,&i,1,&j,&hinv,INSERT_VALUES); 
1013:     }
1014:     if (i>=2*m/3){
1015:       j = i-2*m/3;
1016:       MatSetValues(user->Grad,1,&i,1,&j,&neg_hinv,INSERT_VALUES); 
1017:       j = j + user->mx*user->mx;
1018:       MatSetValues(user->Grad,1,&i,1,&j,&hinv,INSERT_VALUES); 
1019:     }
1020:   }

1022:   MatAssemblyBegin(user->Grad,MAT_FINAL_ASSEMBLY); 
1023:   MatAssemblyEnd(user->Grad,MAT_FINAL_ASSEMBLY); 


1026:   /* Generate arithmetic averaging matrix Av */
1027:   MatCreate(PETSC_COMM_WORLD,&user->Av); 
1028:   MatSetSizes(user->Av,PETSC_DECIDE,PETSC_DECIDE,m,n); 
1029:   MatSetFromOptions(user->Av); 
1030:   MatMPIAIJSetPreallocation(user->Av,2,PETSC_NULL,2,PETSC_NULL); 
1031:   MatSeqAIJSetPreallocation(user->Av,2,PETSC_NULL); 
1032:   MatGetOwnershipRange(user->Av,&istart,&iend); 

1034:   for (i=istart; i<iend; i++){
1035:     if (i<m/3){
1036:       iblock = i / (user->mx-1);
1037:       j = iblock*user->mx + (i % (user->mx-1));
1038:       MatSetValues(user->Av,1,&i,1,&j,&half,INSERT_VALUES); 
1039:       j = j+1;
1040:       MatSetValues(user->Av,1,&i,1,&j,&half,INSERT_VALUES); 
1041:     }
1042:     if (i>=m/3 && i<2*m/3){
1043:       iblock = (i-m/3) / (user->mx*(user->mx-1));
1044:       j = iblock*user->mx*user->mx + ((i-m/3) % (user->mx*(user->mx-1)));
1045:       MatSetValues(user->Av,1,&i,1,&j,&half,INSERT_VALUES); 
1046:       j = j + user->mx;
1047:       MatSetValues(user->Av,1,&i,1,&j,&half,INSERT_VALUES); 
1048:     }
1049:     if (i>=2*m/3){
1050:       j = i-2*m/3;
1051:       MatSetValues(user->Av,1,&i,1,&j,&half,INSERT_VALUES); 
1052:       j = j + user->mx*user->mx;
1053:       MatSetValues(user->Av,1,&i,1,&j,&half,INSERT_VALUES); 
1054:     }
1055:   }

1057:   MatAssemblyBegin(user->Av,MAT_FINAL_ASSEMBLY); 
1058:   MatAssemblyEnd(user->Av,MAT_FINAL_ASSEMBLY); 

1060:   /* Generate transpose of averaging matrix Av */
1061:   MatTranspose(user->Av,MAT_INITIAL_MATRIX,&user->AvT); 



1065:   MatCreate(PETSC_COMM_WORLD,&user->L); 
1066:   MatSetSizes(user->L,PETSC_DECIDE,PETSC_DECIDE,m+n,n); 
1067:   MatSetFromOptions(user->L); 
1068:   MatMPIAIJSetPreallocation(user->L,2,PETSC_NULL,2,PETSC_NULL); 
1069:   MatSeqAIJSetPreallocation(user->L,2,PETSC_NULL); 
1070:   MatGetOwnershipRange(user->L,&istart,&iend);

1072:   for (i=istart; i<iend; i++){
1073:     if (i<m/3){
1074:       iblock = i / (user->mx-1);
1075:       j = iblock*user->mx + (i % (user->mx-1));
1076:       MatSetValues(user->L,1,&i,1,&j,&neg_hinv,INSERT_VALUES); 
1077:       j = j+1;
1078:       MatSetValues(user->L,1,&i,1,&j,&hinv,INSERT_VALUES); 
1079:     }
1080:     if (i>=m/3 && i<2*m/3){
1081:       iblock = (i-m/3) / (user->mx*(user->mx-1));
1082:       j = iblock*user->mx*user->mx + ((i-m/3) % (user->mx*(user->mx-1)));
1083:       MatSetValues(user->L,1,&i,1,&j,&neg_hinv,INSERT_VALUES); 
1084:       j = j + user->mx;
1085:       MatSetValues(user->L,1,&i,1,&j,&hinv,INSERT_VALUES); 
1086:     }
1087:     if (i>=2*m/3 && i<m){
1088:       j = i-2*m/3;
1089:       MatSetValues(user->L,1,&i,1,&j,&neg_hinv,INSERT_VALUES); 
1090:       j = j + user->mx*user->mx;
1091:       MatSetValues(user->L,1,&i,1,&j,&hinv,INSERT_VALUES); 
1092:     }
1093:     if (i>=m){
1094:       j = i - m;
1095:       MatSetValues(user->L,1,&i,1,&j,&sqrt_beta,INSERT_VALUES); 
1096:     }
1097:   }

1099:   MatAssemblyBegin(user->L,MAT_FINAL_ASSEMBLY); 
1100:   MatAssemblyEnd(user->L,MAT_FINAL_ASSEMBLY); 

1102:   MatScale(user->L,PetscPowScalar(h,1.5)); 

1104:   /* Generate Div matrix */
1105:   MatTranspose(user->Grad,MAT_INITIAL_MATRIX,&user->Div);


1108:   /* Build work vectors and matrices */
1109:   VecCreate(PETSC_COMM_WORLD,&user->S); 
1110:   VecSetSizes(user->S, PETSC_DECIDE, user->mx*user->mx*(user->mx-1)*3); 
1111:   VecSetFromOptions(user->S); 

1113:   VecCreate(PETSC_COMM_WORLD,&user->lwork); 
1114:   VecSetSizes(user->lwork,PETSC_DECIDE,m+user->mx*user->mx*user->mx);  
1115:   VecSetFromOptions(user->lwork); 

1117:   MatDuplicate(user->Div,MAT_SHARE_NONZERO_PATTERN,&user->Divwork); 

1119:   MatDuplicate(user->Av,MAT_SHARE_NONZERO_PATTERN,&user->Avwork); 

1121:   VecCreate(PETSC_COMM_WORLD,&user->d); 
1122:   VecSetSizes(user->d,PETSC_DECIDE,user->ndata*user->nt); 
1123:   VecSetFromOptions(user->d); 

1125:   VecDuplicate(user->S,&user->Swork); 
1126:   VecDuplicate(user->S,&user->Av_u); 
1127:   VecDuplicate(user->S,&user->Twork); 
1128:   VecDuplicate(user->S,&user->Rwork); 
1129:   VecDuplicate(user->y,&user->ywork); 
1130:   VecDuplicate(user->u,&user->uwork); 
1131:   VecDuplicate(user->u,&user->js_diag); 
1132:   VecDuplicate(user->c,&user->cwork); 
1133:   VecDuplicate(user->d,&user->dwork); 

1135:   /* Create matrix-free shell user->Js for computing A*x */
1136:   MatCreateShell(PETSC_COMM_WORLD,PETSC_DETERMINE,PETSC_DETERMINE,user->m*user->nt,user->m*user->nt,user,&user->Js); 
1137:   MatShellSetOperation(user->Js,MATOP_MULT,(void(*)(void))StateMatMult); 
1138:   MatShellSetOperation(user->Js,MATOP_DUPLICATE,(void(*)(void))StateMatDuplicate); 
1139:   MatShellSetOperation(user->Js,MATOP_MULT_TRANSPOSE,(void(*)(void))StateMatMultTranspose); 
1140:   MatShellSetOperation(user->Js,MATOP_GET_DIAGONAL,(void(*)(void))StateMatGetDiagonal); 

1142:   /* Diagonal blocks of user->Js */
1143:   MatCreateShell(PETSC_COMM_WORLD,PETSC_DETERMINE,PETSC_DETERMINE,user->m,user->m,user,&user->JsBlock); 
1144:   MatShellSetOperation(user->JsBlock,MATOP_MULT,(void(*)(void))StateMatBlockMult); 
1145:   /* Blocks are symmetric */
1146:   MatShellSetOperation(user->JsBlock,MATOP_MULT_TRANSPOSE,(void(*)(void))StateMatBlockMult); 

1148:   /* Create a matrix-free shell user->JsBlockPrec for computing (U+D)\D*(L+D)\x, where JsBlock = L+D+U,
1149:      D is diagonal, L is strictly lower triangular, and U is strictly upper triangular.
1150:      This is an SSOR preconditioner for user->JsBlock. */
1151:   MatCreateShell(PETSC_COMM_WORLD,PETSC_DETERMINE,PETSC_DETERMINE,user->m,user->m,user,&user->JsBlockPrec); 
1152:   MatShellSetOperation(user->JsBlockPrec,MATOP_MULT,(void(*)(void))StateMatBlockPrecMult);  
1153:   /* JsBlockPrec is symmetric */
1154:   MatShellSetOperation(user->JsBlockPrec,MATOP_MULT_TRANSPOSE,(void(*)(void))StateMatBlockPrecMult); 
1155:   MatSetOption(user->JsBlockPrec,MAT_SYMMETRY_ETERNAL,PETSC_TRUE); 
1156:   
1157:   /* Create a matrix-free shell user->Jd for computing B*x */
1158:   MatCreateShell(PETSC_COMM_WORLD,PETSC_DETERMINE,PETSC_DETERMINE,user->m*user->nt,user->m,user,&user->Jd); 
1159:   MatShellSetOperation(user->Jd,MATOP_MULT,(void(*)(void))DesignMatMult); 
1160:   MatShellSetOperation(user->Jd,MATOP_MULT_TRANSPOSE,(void(*)(void))DesignMatMultTranspose); 

1162:   /* User-defined routines for computing user->Js\x and user->Js^T\x*/
1163:   MatCreateShell(PETSC_COMM_WORLD,PETSC_DETERMINE,PETSC_DETERMINE,user->m*user->nt,user->m*user->nt,user,&user->JsInv); 
1164:   MatShellSetOperation(user->JsInv,MATOP_MULT,(void(*)(void))StateMatInvMult); 
1165:   MatShellSetOperation(user->JsInv,MATOP_MULT_TRANSPOSE,(void(*)(void))StateMatInvTransposeMult); 

1167:   /* Solver options and tolerances */
1168:   KSPCreate(PETSC_COMM_WORLD,&user->solver); 
1169:   KSPSetType(user->solver,KSPCG); 
1170:   KSPSetOperators(user->solver,user->JsBlock,user->JsBlockPrec,DIFFERENT_NONZERO_PATTERN); 
1171:   KSPSetInitialGuessNonzero(user->solver,PETSC_FALSE); 
1172:   KSPSetTolerances(user->solver,1e-4,1e-20,1e3,500); 
1173:   KSPGetPC(user->solver,&user->prec); 
1174:   PCSetType(user->prec,PCSHELL); 

1176:   PCShellSetApply(user->prec,StateMatBlockPrecMult); 
1177:   PCShellSetApplyTranspose(user->prec,StateMatBlockPrecMult); 
1178:   PCShellSetContext(user->prec,user); 


1181:   /* Create scatter from y to y_1,y_2,...,y_nt */
1182:   PetscMalloc(user->nt*user->m*sizeof(PetscInt),&user->yi_scatter);
1183:   VecCreate(PETSC_COMM_WORLD,&yi); 
1184:   VecSetSizes(yi,PETSC_DECIDE,user->mx*user->mx*user->mx); 
1185:   VecSetFromOptions(yi); 
1186:   VecDuplicateVecs(yi,user->nt,&user->yi); 
1187:   VecDuplicateVecs(yi,user->nt,&user->yiwork); 

1189:   VecGetOwnershipRange(user->y,&lo2,&hi2); 
1190:   istart = 0;
1191:   for (i=0; i<user->nt; i++){
1192:     VecGetOwnershipRange(user->yi[i],&lo,&hi); 
1193:     ISCreateStride(PETSC_COMM_SELF,hi-lo,lo,1,&is_to_yi); 
1194:     ISCreateStride(PETSC_COMM_SELF,hi-lo,lo2+istart,1,&is_from_y); 
1195:     VecScatterCreate(user->y,is_from_y,user->yi[i],is_to_yi,&user->yi_scatter[i]); 
1196:     istart = istart + hi-lo;
1197:     ISDestroy(&is_to_yi); 
1198:     ISDestroy(&is_from_y); 
1199:   }
1200:   VecDestroy(&yi); 

1202:   /* Create scatter from d to d_1,d_2,...,d_ns */
1203:   PetscMalloc(user->ns*user->ndata*sizeof(PetscInt),&user->di_scatter);
1204:   VecCreate(PETSC_COMM_WORLD,&di); 
1205:   VecSetSizes(di,PETSC_DECIDE,user->ndata); 
1206:   VecSetFromOptions(di); 
1207:   VecDuplicateVecs(di,user->ns,&user->di); 
1208:   VecGetOwnershipRange(user->d,&lo2,&hi2); 
1209:   istart = 0;
1210:   for (i=0; i<user->ns; i++){
1211:     VecGetOwnershipRange(user->di[i],&lo,&hi); 
1212:     ISCreateStride(PETSC_COMM_SELF,hi-lo,lo,1,&is_to_di); 
1213:     ISCreateStride(PETSC_COMM_SELF,hi-lo,lo2+istart,1,&is_from_d); 
1214:     VecScatterCreate(user->d,is_from_d,user->di[i],is_to_di,&user->di_scatter[i]); 
1215:     istart = istart + hi-lo;
1216:     ISDestroy(&is_to_di); 
1217:     ISDestroy(&is_from_d); 
1218:   }
1219:   VecDestroy(&di); 

1221:   /* Assemble RHS of forward problem */
1222:   VecCopy(bc,user->yiwork[0]);
1223:   for (i=1; i<user->nt; i++){
1224:     VecSet(user->yiwork[i],0.0); 
1225:   }
1226:   Gather_i(user->q,user->yiwork,user->yi_scatter,user->nt); 
1227:   VecDestroy(&bc); 

1229:   /* Compute true state function yt given ut */
1230:   VecCreate(PETSC_COMM_WORLD,&user->ytrue); 
1231:   VecSetSizes(user->ytrue,PETSC_DECIDE,n*user->nt); 
1232:   VecSetFromOptions(user->ytrue); 

1234:   /* First compute Av_u = Av*exp(-u) */
1235:   VecSet(user->uwork,0);
1236:   VecAXPY(user->uwork,-1.0,user->utrue); /*  Note: user->utrue */
1237:   VecExp(user->uwork); 
1238:   MatMult(user->Av,user->uwork,user->Av_u); 

1240:   MatMatMultSymbolic(user->Div,user->Grad,PETSC_DEFAULT,&user->DSG); 
1241:   user->dsg_formed = PETSC_TRUE;

1243:   /* Next form DSG = Div*S*Grad */
1244:   MatCopy(user->Div,user->Divwork,SAME_NONZERO_PATTERN); 
1245:   MatDiagonalScale(user->Divwork,PETSC_NULL,user->Av_u); 
1246:   if (user->dsg_formed) {
1247:     MatMatMultNumeric(user->Divwork,user->Grad,user->DSG); 
1248:   } else {
1249:     MatMatMultSymbolic(user->Div,user->Grad,PETSC_DEFAULT,&user->DSG); 
1250:     MatMatMultNumeric(user->Div,user->Grad,user->DSG); 
1251:     user->dsg_formed = PETSC_TRUE;
1252:   }
1253:   /* B = speye(nx^3) + ht*DSG; */
1254:   MatScale(user->DSG,user->ht); 
1255:   MatShift(user->DSG,1.0); 

1257:   /* Now solve for ytrue */
1258:   StateMatInvMult(user->Js,user->q,user->ytrue); 

1260:   /* Initial guess y0 for state given u0 */

1262:   /* First compute Av_u = Av*exp(-u) */
1263:   VecSet(user->uwork,0);
1264:   VecAXPY(user->uwork,-1.0,user->u); /*  Note: user->u */
1265:   VecExp(user->uwork); 
1266:   MatMult(user->Av,user->uwork,user->Av_u); 

1268:   /* Next form DSG = Div*S*Grad */
1269:   MatCopy(user->Div,user->Divwork,SAME_NONZERO_PATTERN); 
1270:   MatDiagonalScale(user->Divwork,PETSC_NULL,user->Av_u); 
1271:   if (user->dsg_formed) {
1272:     MatMatMultNumeric(user->Divwork,user->Grad,user->DSG); 
1273:   } else {
1274:     MatMatMultSymbolic(user->Div,user->Grad,PETSC_DEFAULT,&user->DSG); 
1275:     MatMatMultNumeric(user->Div,user->Grad,user->DSG); 
1276:     user->dsg_formed = PETSC_TRUE;
1277:   }
1278:   /* B = speye(nx^3) + ht*DSG; */
1279:   MatScale(user->DSG,user->ht); 
1280:   MatShift(user->DSG,1.0); 

1282:   /* Now solve for y */
1283:   user->lcl->solve_type = LCL_FORWARD1;
1284:   StateMatInvMult(user->Js,user->q,user->y); 
1285:   
1286:   /* Construct projection matrix Q, a block diagonal matrix consisting of nt copies of Qblock along the diagonal */
1287:   MatCreate(PETSC_COMM_WORLD,&user->Qblock); 
1288:   MatSetSizes(user->Qblock,PETSC_DECIDE,PETSC_DECIDE,user->ndata,n); 
1289:   MatSetFromOptions(user->Qblock); 
1290:   MatMPIAIJSetPreallocation(user->Qblock,8,PETSC_NULL,8,PETSC_NULL); 
1291:   MatSeqAIJSetPreallocation(user->Qblock,8,PETSC_NULL); 

1293:  
1294:   for (i=0; i<user->mx; i++){
1295:     x[i] = h*(i+0.5);
1296:     y[i] = h*(i+0.5);
1297:     z[i] = h*(i+0.5);
1298:   }
1299:   
1300:   MatGetOwnershipRange(user->Qblock,&istart,&iend);

1302:   nx = user->mx; ny = user->mx; nz = user->mx;
1303:   for (i=istart; i<iend; i++){
1304:     
1305:     xri = xr[i];
1306:     im = 0;
1307:     xim = x[im];
1308:     while (xri>xim && im<nx){
1309:       im = im+1;
1310:       xim = x[im];
1311:     }
1312:     indx1 = im-1;
1313:     indx2 = im;
1314:     dx1 = xri - x[indx1];
1315:     dx2 = x[indx2] - xri;

1317:     yri = yr[i];
1318:     im = 0;
1319:     yim = y[im];
1320:     while (yri>yim && im<ny){
1321:       im = im+1;
1322:       yim = y[im];
1323:     }
1324:     indy1 = im-1;
1325:     indy2 = im;
1326:     dy1 = yri - y[indy1];
1327:     dy2 = y[indy2] - yri;
1328:     
1329:     zri = zr[i];
1330:     im = 0;
1331:     zim = z[im];
1332:     while (zri>zim && im<nz){
1333:       im = im+1;
1334:       zim = z[im];
1335:     }
1336:     indz1 = im-1;
1337:     indz2 = im;
1338:     dz1 = zri - z[indz1];
1339:     dz2 = z[indz2] - zri;

1341:     Dx = x[indx2] - x[indx1];
1342:     Dy = y[indy2] - y[indy1];
1343:     Dz = z[indz2] - z[indz1];

1345:     j = indx1 + indy1*nx + indz1*nx*ny;
1346:     v = (1-dx1/Dx)*(1-dy1/Dy)*(1-dz1/Dz);
1347:     MatSetValues(user->Qblock,1,&i,1,&j,&v,INSERT_VALUES); 

1349:     j = indx1 + indy1*nx + indz2*nx*ny;
1350:     v = (1-dx1/Dx)*(1-dy1/Dy)*(1-dz2/Dz);
1351:     MatSetValues(user->Qblock,1,&i,1,&j,&v,INSERT_VALUES); 

1353:     j = indx1 + indy2*nx + indz1*nx*ny;
1354:     v = (1-dx1/Dx)*(1-dy2/Dy)*(1-dz1/Dz);
1355:     MatSetValues(user->Qblock,1,&i,1,&j,&v,INSERT_VALUES); 

1357:     j = indx1 + indy2*nx + indz2*nx*ny;
1358:     v = (1-dx1/Dx)*(1-dy2/Dy)*(1-dz2/Dz);
1359:     MatSetValues(user->Qblock,1,&i,1,&j,&v,INSERT_VALUES); 

1361:     j = indx2 + indy1*nx + indz1*nx*ny;
1362:     v = (1-dx2/Dx)*(1-dy1/Dy)*(1-dz1/Dz);
1363:     MatSetValues(user->Qblock,1,&i,1,&j,&v,INSERT_VALUES); 

1365:     j = indx2 + indy1*nx + indz2*nx*ny;
1366:     v = (1-dx2/Dx)*(1-dy1/Dy)*(1-dz2/Dz);
1367:     MatSetValues(user->Qblock,1,&i,1,&j,&v,INSERT_VALUES); 

1369:     j = indx2 + indy2*nx + indz1*nx*ny;
1370:     v = (1-dx2/Dx)*(1-dy2/Dy)*(1-dz1/Dz);
1371:     MatSetValues(user->Qblock,1,&i,1,&j,&v,INSERT_VALUES); 

1373:     j = indx2 + indy2*nx + indz2*nx*ny;
1374:     v = (1-dx2/Dx)*(1-dy2/Dy)*(1-dz2/Dz);
1375:     MatSetValues(user->Qblock,1,&i,1,&j,&v,INSERT_VALUES); 

1377:   }

1379:   MatAssemblyBegin(user->Qblock,MAT_FINAL_ASSEMBLY); 
1380:   MatAssemblyEnd(user->Qblock,MAT_FINAL_ASSEMBLY); 


1383:   MatTranspose(user->Qblock,MAT_INITIAL_MATRIX,&user->QblockT); 
1384:  

1386:   MatTranspose(user->L,MAT_INITIAL_MATRIX,&user->LT); 

1388:   /* Add noise to the measurement data */
1389:   VecSet(user->ywork,1.0); 
1390:   VecAYPX(user->ywork,user->noise,user->ytrue); 
1391:   Scatter_i(user->ywork,user->yiwork,user->yi_scatter,user->nt); 
1392:   for (j=0; j<user->ns; j++){
1393:     i = user->sample_times[j];
1394:     MatMult(user->Qblock,user->yiwork[i],user->di[j]);
1395:   }
1396:   Gather_i(user->d,user->di,user->di_scatter,user->ns); 

1398:   /* Now that initial conditions have been set, let the user pass tolerance options to the KSP solver */
1399:   KSPSetFromOptions(user->solver); 
1400:   PetscFree(x); 
1401:   PetscFree(y); 
1402:   PetscFree(z); 


1405:   return(0);
1406: }

1410: PetscErrorCode ParabolicDestroy(AppCtx *user)
1411: {
1413:   PetscInt i;
1415:   MatDestroy(&user->Qblock); 
1416:   MatDestroy(&user->QblockT); 
1417:   MatDestroy(&user->Div); 
1418:   MatDestroy(&user->Divwork); 
1419:   MatDestroy(&user->Grad); 
1420:   MatDestroy(&user->Av); 
1421:   MatDestroy(&user->Avwork); 
1422:   MatDestroy(&user->AvT); 
1423:   MatDestroy(&user->DSG); 
1424:   MatDestroy(&user->L); 
1425:   MatDestroy(&user->LT); 
1426:   KSPDestroy(&user->solver); 
1427:   MatDestroy(&user->Js); 
1428:   MatDestroy(&user->Jd); 
1429:   MatDestroy(&user->JsInv); 
1430:   MatDestroy(&user->JsBlock); 
1431:   MatDestroy(&user->JsBlockPrec); 
1432:   VecDestroy(&user->u); 
1433:   VecDestroy(&user->uwork); 
1434:   VecDestroy(&user->utrue); 
1435:   VecDestroy(&user->y); 
1436:   VecDestroy(&user->ywork); 
1437:   VecDestroy(&user->ytrue); 
1438:   VecDestroyVecs(user->nt,&user->yi); 
1439:   VecDestroyVecs(user->nt,&user->yiwork); 
1440:   VecDestroyVecs(user->ns,&user->di); 
1441:   PetscFree(user->yi); 
1442:   PetscFree(user->yiwork); 
1443:   PetscFree(user->di); 
1444:   VecDestroy(&user->c); 
1445:   VecDestroy(&user->cwork); 
1446:   VecDestroy(&user->ur); 
1447:   VecDestroy(&user->q); 
1448:   VecDestroy(&user->d); 
1449:   VecDestroy(&user->dwork); 
1450:   VecDestroy(&user->lwork); 
1451:   VecDestroy(&user->S); 
1452:   VecDestroy(&user->Swork); 
1453:   VecDestroy(&user->Av_u); 
1454:   VecDestroy(&user->Twork); 
1455:   VecDestroy(&user->Rwork); 
1456:   VecDestroy(&user->js_diag); 
1457:   ISDestroy(&user->s_is); 
1458:   ISDestroy(&user->d_is); 
1459:   VecScatterDestroy(&user->state_scatter); 
1460:   VecScatterDestroy(&user->design_scatter); 
1461:   for (i=0; i<user->nt; i++){
1462:     VecScatterDestroy(&user->yi_scatter[i]); 
1463:   }
1464:   for (i=0; i<user->ns; i++){
1465:     VecScatterDestroy(&user->di_scatter[i]); 
1466:   }
1467:   PetscFree(user->yi_scatter); 
1468:   PetscFree(user->di_scatter); 
1469:   PetscFree(user->sample_times); 
1470:   return(0);
1471: }

1475: PetscErrorCode ParabolicMonitor(TaoSolver tao, void *ptr)
1476: {
1478:   Vec X;
1479:   PetscReal unorm,ynorm;
1480:   AppCtx *user = (AppCtx*)ptr;
1482:   TaoGetSolutionVector(tao,&X); 
1483:   Scatter(X,user->ywork,user->state_scatter,user->uwork,user->design_scatter); 
1484:   VecAXPY(user->ywork,-1.0,user->ytrue); 
1485:   VecAXPY(user->uwork,-1.0,user->utrue); 
1486:   VecNorm(user->uwork,NORM_2,&unorm); 
1487:   VecNorm(user->ywork,NORM_2,&ynorm); 
1488:   PetscPrintf(MPI_COMM_WORLD, "||u-ut||=%G ||y-yt||=%G\n",unorm,ynorm); 
1489:   return(0);
1490: }