Actual source code: elliptic.c

  1: #include <petsc/private/taoimpl.h>

  3: typedef struct {
  4:   PetscInt n; /* Number of total variables */
  5:   PetscInt m; /* Number of constraints */
  6:   PetscInt nstate;
  7:   PetscInt ndesign;
  8:   PetscInt mx; /* grid points in each direction */
  9:   PetscInt ns; /* Number of data samples (1<=ns<=8)
 10:                   Currently only ns=1 is supported */
 11:   PetscInt ndata; /* Number of data points per sample */
 12:   IS       s_is;
 13:   IS       d_is;

 15:   VecScatter state_scatter;
 16:   VecScatter design_scatter;
 17:   VecScatter *yi_scatter, *di_scatter;
 18:   Vec        suby,subq,subd;
 19:   Mat        Js,Jd,JsPrec,JsInv,JsBlock;

 21:   PetscReal alpha; /* Regularization parameter */
 22:   PetscReal beta; /* Weight attributed to ||u||^2 in regularization functional */
 23:   PetscReal noise; /* Amount of noise to add to data */
 24:   PetscReal *ones;
 25:   Mat       Q;
 26:   Mat       MQ;
 27:   Mat       L;

 29:   Mat Grad;
 30:   Mat Av,Avwork;
 31:   Mat Div, Divwork;
 32:   Mat DSG;
 33:   Mat Diag,Ones;

 35:   Vec q;
 36:   Vec ur; /* reference */

 38:   Vec d;
 39:   Vec dwork;

 41:   Vec x; /* super vec of y,u */

 43:   Vec y; /* state variables */
 44:   Vec ywork;

 46:   Vec ytrue;

 48:   Vec u; /* design variables */
 49:   Vec uwork;

 51:   Vec utrue;

 53:   Vec js_diag;

 55:   Vec c; /* constraint vector */
 56:   Vec cwork;

 58:   Vec lwork;
 59:   Vec S;
 60:   Vec Swork,Twork,Sdiag,Ywork;
 61:   Vec Av_u;

 63:   KSP solver;
 64:   PC  prec;

 66:   PetscReal tola,tolb,tolc,told;
 67:   PetscInt  ksp_its;
 68:   PetscInt  ksp_its_initial;
 69:   PetscLogStage stages[10];
 70:   PetscBool use_ptap;
 71:   PetscBool use_lrc;
 72: } AppCtx;

 74: PetscErrorCode FormFunction(Tao, Vec, PetscReal*, void*);
 75: PetscErrorCode FormGradient(Tao, Vec, Vec, void*);
 76: PetscErrorCode FormFunctionGradient(Tao, Vec, PetscReal*, Vec, void*);
 77: PetscErrorCode FormJacobianState(Tao, Vec, Mat, Mat, Mat, void*);
 78: PetscErrorCode FormJacobianDesign(Tao, Vec, Mat,void*);
 79: PetscErrorCode FormConstraints(Tao, Vec, Vec, void*);
 80: PetscErrorCode FormHessian(Tao, Vec, Mat, Mat, void*);
 81: PetscErrorCode Gather(Vec, Vec, VecScatter, Vec, VecScatter);
 82: PetscErrorCode Scatter(Vec, Vec, VecScatter, Vec, VecScatter);
 83: PetscErrorCode EllipticInitialize(AppCtx*);
 84: PetscErrorCode EllipticDestroy(AppCtx*);
 85: PetscErrorCode EllipticMonitor(Tao, void*);

 87: PetscErrorCode StateBlockMatMult(Mat,Vec,Vec);
 88: PetscErrorCode StateMatMult(Mat,Vec,Vec);

 90: PetscErrorCode StateInvMatMult(Mat,Vec,Vec);
 91: PetscErrorCode DesignMatMult(Mat,Vec,Vec);
 92: PetscErrorCode DesignMatMultTranspose(Mat,Vec,Vec);

 94: PetscErrorCode QMatMult(Mat,Vec,Vec);
 95: PetscErrorCode QMatMultTranspose(Mat,Vec,Vec);

 97: static  char help[]="";

 99: int main(int argc, char **argv)
100: {
101:   PetscErrorCode     ierr;
102:   Vec                x0;
103:   Tao                tao;
104:   AppCtx             user;
105:   PetscInt           ntests = 1;
106:   PetscInt           i;

108:   PetscInitialize(&argc, &argv, (char*)0,help);
109:   user.mx = 8;
110:   PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"elliptic example",NULL);
111:   PetscOptionsInt("-mx","Number of grid points in each direction","",user.mx,&user.mx,NULL);
112:   user.ns = 6;
113:   PetscOptionsInt("-ns","Number of data samples (1<=ns<=8)","",user.ns,&user.ns,NULL);
114:   user.ndata = 64;
115:   PetscOptionsInt("-ndata","Numbers of data points per sample","",user.ndata,&user.ndata,NULL);
116:   user.alpha = 0.1;
117:   PetscOptionsReal("-alpha","Regularization parameter","",user.alpha,&user.alpha,NULL);
118:   user.beta = 0.00001;
119:   PetscOptionsReal("-beta","Weight attributed to ||u||^2 in regularization functional","",user.beta,&user.beta,NULL);
120:   user.noise = 0.01;
121:   PetscOptionsReal("-noise","Amount of noise to add to data","",user.noise,&user.noise,NULL);

123:   user.use_ptap = PETSC_FALSE;
124:   PetscOptionsBool("-use_ptap","Use ptap matrix for DSG","",user.use_ptap,&user.use_ptap,NULL);
125:   user.use_lrc = PETSC_FALSE;
126:   PetscOptionsBool("-use_lrc","Use lrc matrix for Js","",user.use_lrc,&user.use_lrc,NULL);
127:   PetscOptionsInt("-ntests","Number of times to repeat TaoSolve","",ntests,&ntests,NULL);
128:   PetscOptionsEnd();

130:   user.m = user.ns*user.mx*user.mx*user.mx; /* number of constraints */
131:   user.nstate =  user.m;
132:   user.ndesign = user.mx*user.mx*user.mx;
133:   user.n = user.nstate + user.ndesign; /* number of variables */

135:   /* Create TAO solver and set desired solution method */
136:   TaoCreate(PETSC_COMM_WORLD,&tao);
137:   TaoSetType(tao,TAOLCL);

139:   /* Set up initial vectors and matrices */
140:   EllipticInitialize(&user);

142:   Gather(user.x,user.y,user.state_scatter,user.u,user.design_scatter);
143:   VecDuplicate(user.x,&x0);
144:   VecCopy(user.x,x0);

146:   /* Set solution vector with an initial guess */
147:   TaoSetSolution(tao,user.x);
148:   TaoSetObjective(tao, FormFunction, &user);
149:   TaoSetGradient(tao, NULL, FormGradient, &user);
150:   TaoSetConstraintsRoutine(tao, user.c, FormConstraints, &user);

152:   TaoSetJacobianStateRoutine(tao, user.Js, NULL, user.JsInv, FormJacobianState, &user);
153:   TaoSetJacobianDesignRoutine(tao, user.Jd, FormJacobianDesign, &user);

155:   TaoSetStateDesignIS(tao,user.s_is,user.d_is);
156:   TaoSetFromOptions(tao);

158:   /* SOLVE THE APPLICATION */
159:   PetscLogStageRegister("Trials",&user.stages[1]);
160:   PetscLogStagePush(user.stages[1]);
161:   for (i=0; i<ntests; i++) {
162:     TaoSolve(tao);
163:     PetscPrintf(PETSC_COMM_WORLD,"KSP Iterations = %D\n",user.ksp_its);
164:     VecCopy(x0,user.x);
165:   }
166:   PetscLogStagePop();
167:   PetscBarrier((PetscObject)user.x);
168:   PetscPrintf(PETSC_COMM_WORLD,"KSP iterations within initialization: ");
169:   PetscPrintf(PETSC_COMM_WORLD,"%D\n",user.ksp_its_initial);

171:   TaoDestroy(&tao);
172:   VecDestroy(&x0);
173:   EllipticDestroy(&user);
174:   PetscFinalize();
175:   return 0;
176: }
177: /* ------------------------------------------------------------------- */
178: /*
179:    dwork = Qy - d
180:    lwork = L*(u-ur)
181:    f = 1/2 * (dwork.dwork + alpha*lwork.lwork)
182: */
183: PetscErrorCode FormFunction(Tao tao,Vec X,PetscReal *f,void *ptr)
184: {
185:   PetscReal      d1=0,d2=0;
186:   AppCtx         *user = (AppCtx*)ptr;

188:   Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);
189:   MatMult(user->MQ,user->y,user->dwork);
190:   VecAXPY(user->dwork,-1.0,user->d);
191:   VecDot(user->dwork,user->dwork,&d1);
192:   VecWAXPY(user->uwork,-1.0,user->ur,user->u);
193:   MatMult(user->L,user->uwork,user->lwork);
194:   VecDot(user->lwork,user->lwork,&d2);
195:   *f = 0.5 * (d1 + user->alpha*d2);
196:   return 0;
197: }

199: /* ------------------------------------------------------------------- */
200: /*
201:     state: g_s = Q' *(Qy - d)
202:     design: g_d = alpha*L'*L*(u-ur)
203: */
204: PetscErrorCode FormGradient(Tao tao,Vec X,Vec G,void *ptr)
205: {
206:   AppCtx         *user = (AppCtx*)ptr;

208:   Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);
209:   MatMult(user->MQ,user->y,user->dwork);
210:   VecAXPY(user->dwork,-1.0,user->d);
211:   MatMultTranspose(user->MQ,user->dwork,user->ywork);
212:   VecWAXPY(user->uwork,-1.0,user->ur,user->u);
213:   MatMult(user->L,user->uwork,user->lwork);
214:   MatMultTranspose(user->L,user->lwork,user->uwork);
215:   VecScale(user->uwork, user->alpha);
216:   Gather(G,user->ywork,user->state_scatter,user->uwork,user->design_scatter);
217:   return 0;
218: }

220: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *f, Vec G, void *ptr)
221: {
222:   PetscReal      d1,d2;
223:   AppCtx         *user = (AppCtx*)ptr;

225:   Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);
226:   MatMult(user->MQ,user->y,user->dwork);
227:   VecAXPY(user->dwork,-1.0,user->d);
228:   VecDot(user->dwork,user->dwork,&d1);
229:   MatMultTranspose(user->MQ,user->dwork,user->ywork);

231:   VecWAXPY(user->uwork,-1.0,user->ur,user->u);
232:   MatMult(user->L,user->uwork,user->lwork);
233:   VecDot(user->lwork,user->lwork,&d2);
234:   MatMultTranspose(user->L,user->lwork,user->uwork);
235:   VecScale(user->uwork, user->alpha);
236:   *f = 0.5 * (d1 + user->alpha*d2);
237:   Gather(G,user->ywork,user->state_scatter,user->uwork,user->design_scatter);
238:   return 0;
239: }

241: /* ------------------------------------------------------------------- */
242: /* A
243: MatShell object
244: */
245: PetscErrorCode FormJacobianState(Tao tao, Vec X, Mat J, Mat JPre, Mat JInv, void *ptr)
246: {
247:   AppCtx         *user = (AppCtx*)ptr;

249:   Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);
250:   /* DSG = Div * (1/Av_u) * Grad */
251:   VecSet(user->uwork,0);
252:   VecAXPY(user->uwork,-1.0,user->u);
253:   VecExp(user->uwork);
254:   MatMult(user->Av,user->uwork,user->Av_u);
255:   VecCopy(user->Av_u,user->Swork);
256:   VecReciprocal(user->Swork);
257:   if (user->use_ptap) {
258:     MatDiagonalSet(user->Diag,user->Swork,INSERT_VALUES);
259:     MatPtAP(user->Diag,user->Grad,MAT_REUSE_MATRIX,1.0,&user->DSG);
260:   } else {
261:     MatCopy(user->Div,user->Divwork,SAME_NONZERO_PATTERN);
262:     MatDiagonalScale(user->Divwork,NULL,user->Swork);
263:     MatProductNumeric(user->DSG);
264:   }
265:   return 0;
266: }
267: /* ------------------------------------------------------------------- */
268: /* B */
269: PetscErrorCode FormJacobianDesign(Tao tao, Vec X, Mat J, void *ptr)
270: {
271:   AppCtx         *user = (AppCtx*)ptr;

273:   Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);
274:   return 0;
275: }

277: PetscErrorCode StateBlockMatMult(Mat J_shell, Vec X, Vec Y)
278: {
279:   PetscReal      sum;
280:   AppCtx         *user;

282:   MatShellGetContext(J_shell,&user);
283:   MatMult(user->DSG,X,Y);
284:   VecSum(X,&sum);
285:   sum /= user->ndesign;
286:   VecShift(Y,sum);
287:   return 0;
288: }

290: PetscErrorCode StateMatMult(Mat J_shell, Vec X, Vec Y)
291: {
292:   PetscInt       i;
293:   AppCtx         *user;

295:   MatShellGetContext(J_shell,&user);
296:   if (user->ns == 1) {
297:     MatMult(user->JsBlock,X,Y);
298:   } else {
299:     for (i=0;i<user->ns;i++) {
300:       Scatter(X,user->subq,user->yi_scatter[i],0,0);
301:       Scatter(Y,user->suby,user->yi_scatter[i],0,0);
302:       MatMult(user->JsBlock,user->subq,user->suby);
303:       Gather(Y,user->suby,user->yi_scatter[i],0,0);
304:     }
305:   }
306:   return 0;
307: }

309: PetscErrorCode StateInvMatMult(Mat J_shell, Vec X, Vec Y)
310: {
311:   PetscInt       its,i;
312:   AppCtx         *user;

314:   MatShellGetContext(J_shell,&user);
315:   KSPSetOperators(user->solver,user->JsBlock,user->DSG);
316:   if (Y == user->ytrue) {
317:     /* First solve is done using true solution to set up problem */
318:     KSPSetTolerances(user->solver,1e-8,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
319:   } else {
320:     KSPSetTolerances(user->solver,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
321:   }
322:   if (user->ns == 1) {
323:     KSPSolve(user->solver,X,Y);
324:     KSPGetIterationNumber(user->solver,&its);
325:     user->ksp_its+=its;
326:   } else {
327:     for (i=0;i<user->ns;i++) {
328:       Scatter(X,user->subq,user->yi_scatter[i],0,0);
329:       Scatter(Y,user->suby,user->yi_scatter[i],0,0);
330:       KSPSolve(user->solver,user->subq,user->suby);
331:       KSPGetIterationNumber(user->solver,&its);
332:       user->ksp_its+=its;
333:       Gather(Y,user->suby,user->yi_scatter[i],0,0);
334:     }
335:   }
336:   return 0;
337: }
338: PetscErrorCode QMatMult(Mat J_shell, Vec X, Vec Y)
339: {
340:   AppCtx         *user;
341:   PetscInt       i;

343:   MatShellGetContext(J_shell,&user);
344:   if (user->ns == 1) {
345:     MatMult(user->Q,X,Y);
346:   } else {
347:     for (i=0;i<user->ns;i++) {
348:       Scatter(X,user->subq,user->yi_scatter[i],0,0);
349:       Scatter(Y,user->subd,user->di_scatter[i],0,0);
350:       MatMult(user->Q,user->subq,user->subd);
351:       Gather(Y,user->subd,user->di_scatter[i],0,0);
352:     }
353:   }
354:   return 0;
355: }

357: PetscErrorCode QMatMultTranspose(Mat J_shell, Vec X, Vec Y)
358: {
359:   AppCtx         *user;
360:   PetscInt       i;

362:   MatShellGetContext(J_shell,&user);
363:   if (user->ns == 1) {
364:     MatMultTranspose(user->Q,X,Y);
365:   } else {
366:     for (i=0;i<user->ns;i++) {
367:       Scatter(X,user->subd,user->di_scatter[i],0,0);
368:       Scatter(Y,user->suby,user->yi_scatter[i],0,0);
369:       MatMultTranspose(user->Q,user->subd,user->suby);
370:       Gather(Y,user->suby,user->yi_scatter[i],0,0);
371:     }
372:   }
373:   return 0;
374: }

376: PetscErrorCode DesignMatMult(Mat J_shell, Vec X, Vec Y)
377: {
378:   PetscInt       i;
379:   AppCtx         *user;

381:   MatShellGetContext(J_shell,&user);

383:   /* sdiag(1./v) */
384:   VecSet(user->uwork,0);
385:   VecAXPY(user->uwork,-1.0,user->u);
386:   VecExp(user->uwork);

388:   /* sdiag(1./((Av*(1./v)).^2)) */
389:   MatMult(user->Av,user->uwork,user->Swork);
390:   VecPointwiseMult(user->Swork,user->Swork,user->Swork);
391:   VecReciprocal(user->Swork);

393:   /* (Av * (sdiag(1./v) * b)) */
394:   VecPointwiseMult(user->uwork,user->uwork,X);
395:   MatMult(user->Av,user->uwork,user->Twork);

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

400:   if (user->ns == 1) {
401:     /* (sdiag(Grad*y(:,i)) */
402:     MatMult(user->Grad,user->y,user->Twork);

404:     /* Div * (sdiag(Grad*y(:,i)) * (sdiag(1./((Av*(1./v)).^2)) * (Av * (sdiag(1./v) * b)))) */
405:     VecPointwiseMult(user->Swork,user->Twork,user->Swork);
406:     MatMultTranspose(user->Grad,user->Swork,Y);
407:   } else {
408:     for (i=0;i<user->ns;i++) {
409:       Scatter(user->y,user->suby,user->yi_scatter[i],0,0);
410:       Scatter(Y,user->subq,user->yi_scatter[i],0,0);

412:       MatMult(user->Grad,user->suby,user->Twork);
413:       VecPointwiseMult(user->Twork,user->Twork,user->Swork);
414:       MatMultTranspose(user->Grad,user->Twork,user->subq);
415:       Gather(user->y,user->suby,user->yi_scatter[i],0,0);
416:       Gather(Y,user->subq,user->yi_scatter[i],0,0);
417:     }
418:   }
419:   return 0;
420: }

422: PetscErrorCode DesignMatMultTranspose(Mat J_shell, Vec X, Vec Y)
423: {
424:   PetscInt       i;
425:   AppCtx         *user;

427:   MatShellGetContext(J_shell,&user);
428:   VecZeroEntries(Y);

430:   /* Sdiag = 1./((Av*(1./v)).^2) */
431:   VecSet(user->uwork,0);
432:   VecAXPY(user->uwork,-1.0,user->u);
433:   VecExp(user->uwork);
434:   MatMult(user->Av,user->uwork,user->Swork);
435:   VecPointwiseMult(user->Sdiag,user->Swork,user->Swork);
436:   VecReciprocal(user->Sdiag);

438:   for (i=0;i<user->ns;i++) {
439:     Scatter(X,user->subq,user->yi_scatter[i],0,0);
440:     Scatter(user->y,user->suby,user->yi_scatter[i],0,0);

442:     /* Swork = (Div' * b(:,i)) */
443:     MatMult(user->Grad,user->subq,user->Swork);

445:     /* Twork = Grad*y(:,i) */
446:     MatMult(user->Grad,user->suby,user->Twork);

448:     /* Twork = sdiag(Twork) * Swork */
449:     VecPointwiseMult(user->Twork,user->Swork,user->Twork);

451:     /* Swork = pointwisemult(Sdiag,Twork) */
452:     VecPointwiseMult(user->Swork,user->Twork,user->Sdiag);

454:     /* Ywork = Av' * Swork */
455:     MatMultTranspose(user->Av,user->Swork,user->Ywork);

457:     /* Ywork = pointwisemult(uwork,Ywork) */
458:     VecPointwiseMult(user->Ywork,user->uwork,user->Ywork);
459:     VecAXPY(Y,1.0,user->Ywork);
460:     Gather(user->y,user->suby,user->yi_scatter[i],0,0);
461:   }
462:   return 0;
463: }

465: PetscErrorCode FormConstraints(Tao tao, Vec X, Vec C, void *ptr)
466: {
467:    /* C=Ay - q      A = Div * Sigma * Grad + hx*hx*hx*ones(n,n) */
468:    PetscReal      sum;
469:    PetscInt       i;
470:    AppCtx         *user = (AppCtx*)ptr;

472:    Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);
473:    if (user->ns == 1) {
474:      MatMult(user->Grad,user->y,user->Swork);
475:      VecPointwiseDivide(user->Swork,user->Swork,user->Av_u);
476:      MatMultTranspose(user->Grad,user->Swork,C);
477:      VecSum(user->y,&sum);
478:      sum /= user->ndesign;
479:      VecShift(C,sum);
480:    } else {
481:      for (i=0;i<user->ns;i++) {
482:       Scatter(user->y,user->suby,user->yi_scatter[i],0,0);
483:       Scatter(C,user->subq,user->yi_scatter[i],0,0);
484:       MatMult(user->Grad,user->suby,user->Swork);
485:       VecPointwiseDivide(user->Swork,user->Swork,user->Av_u);
486:       MatMultTranspose(user->Grad,user->Swork,user->subq);

488:       VecSum(user->suby,&sum);
489:       sum /= user->ndesign;
490:       VecShift(user->subq,sum);

492:       Gather(user->y,user->suby,user->yi_scatter[i],0,0);
493:       Gather(C,user->subq,user->yi_scatter[i],0,0);
494:      }
495:    }
496:    VecAXPY(C,-1.0,user->q);
497:    return 0;
498: }

500: PetscErrorCode Scatter(Vec x, Vec sub1, VecScatter scat1, Vec sub2, VecScatter scat2)
501: {
502:   VecScatterBegin(scat1,x,sub1,INSERT_VALUES,SCATTER_FORWARD);
503:   VecScatterEnd(scat1,x,sub1,INSERT_VALUES,SCATTER_FORWARD);
504:   if (sub2) {
505:     VecScatterBegin(scat2,x,sub2,INSERT_VALUES,SCATTER_FORWARD);
506:     VecScatterEnd(scat2,x,sub2,INSERT_VALUES,SCATTER_FORWARD);
507:   }
508:   return 0;
509: }

511: PetscErrorCode Gather(Vec x, Vec sub1, VecScatter scat1, Vec sub2, VecScatter scat2)
512: {
513:   VecScatterBegin(scat1,sub1,x,INSERT_VALUES,SCATTER_REVERSE);
514:   VecScatterEnd(scat1,sub1,x,INSERT_VALUES,SCATTER_REVERSE);
515:   if (sub2) {
516:     VecScatterBegin(scat2,sub2,x,INSERT_VALUES,SCATTER_REVERSE);
517:     VecScatterEnd(scat2,sub2,x,INSERT_VALUES,SCATTER_REVERSE);
518:   }
519:   return 0;
520: }

522: PetscErrorCode EllipticInitialize(AppCtx *user)
523: {
524:   PetscInt       m,n,i,j,k,l,linear_index,is,js,ks,ls,istart,iend,iblock;
525:   Vec            XX,YY,ZZ,XXwork,YYwork,ZZwork,UTwork;
526:   PetscReal      *x,*y,*z;
527:   PetscReal      h,meanut;
528:   PetscScalar    hinv,neg_hinv,half = 0.5,sqrt_beta;
529:   PetscInt       im,indx1,indx2,indy1,indy2,indz1,indz2,nx,ny,nz;
530:   IS             is_alldesign,is_allstate;
531:   IS             is_from_d;
532:   IS             is_from_y;
533:   PetscInt       lo,hi,hi2,lo2,ysubnlocal,dsubnlocal;
534:   const PetscInt *ranges, *subranges;
535:   PetscMPIInt    size;
536:   PetscReal      xri,yri,zri,xim,yim,zim,dx1,dx2,dy1,dy2,dz1,dz2,Dx,Dy,Dz;
537:   PetscScalar    v,vx,vy,vz;
538:   PetscInt       offset,subindex,subvec,nrank,kk;

540:   PetscScalar xr[64] = {0.4970,     0.8498,     0.7814,     0.6268,     0.7782,     0.6402,     0.3617,     0.3160,
541:                         0.3610,     0.5298,     0.6987,     0.3331,     0.7962,     0.5596,     0.3866,     0.6774,
542:                         0.5407,     0.4518,     0.6702,     0.6061,     0.7580,     0.8997,     0.5198,     0.8326,
543:                         0.2138,     0.9198,     0.3000,     0.2833,     0.8288,     0.7076,     0.1820,     0.0728,
544:                         0.8447,     0.2367,     0.3239,     0.6413,     0.3114,     0.4731,     0.1192,     0.9273,
545:                         0.5724,     0.4331,     0.5136,     0.3547,     0.4413,     0.2602,     0.5698,     0.7278,
546:                         0.5261,     0.6230,     0.2454,     0.3948,     0.7479,     0.6582,     0.4660,     0.5594,
547:                         0.7574,     0.1143,     0.5900,     0.1065,     0.4260,     0.3294,     0.8276,     0.0756};

549:   PetscScalar yr[64] = {0.7345,     0.9120,     0.9288,     0.7528,     0.4463,     0.4985,     0.2497,     0.6256,
550:                         0.3425,     0.9026,     0.6983,     0.4230,     0.7140,     0.2970,     0.4474,     0.8792,
551:                         0.6604,     0.2485,     0.7968,     0.6127,     0.1796,     0.2437,     0.5938,     0.6137,
552:                         0.3867,     0.5658,     0.4575,     0.1009,     0.0863,     0.3361,     0.0738,     0.3985,
553:                         0.6602,     0.1437,     0.0934,     0.5983,     0.5950,     0.0763,     0.0768,     0.2288,
554:                         0.5761,     0.1129,     0.3841,     0.6150,     0.6904,     0.6686,     0.1361,     0.4601,
555:                         0.4491,     0.3716,     0.1969,     0.6537,     0.6743,     0.6991,     0.4811,     0.5480,
556:                         0.1684,     0.4569,     0.6889,     0.8437,     0.3015,     0.2854,     0.8199,     0.2658};

558:   PetscScalar zr[64] = {0.7668,     0.8573,     0.2654,     0.2719,     0.1060,     0.1311,     0.6232,     0.2295,
559:                         0.8009,     0.2147,     0.2119,     0.9325,     0.4473,     0.3600,     0.3374,     0.3819,
560:                         0.4066,     0.5801,     0.1673,     0.0959,     0.4638,     0.8236,     0.8800,     0.2939,
561:                         0.2028,     0.8262,     0.2706,     0.6276,     0.9085,     0.6443,     0.8241,     0.0712,
562:                         0.1824,     0.7789,     0.4389,     0.8415,     0.7055,     0.6639,     0.3653,     0.2078,
563:                         0.1987,     0.2297,     0.4321,     0.8115,     0.4915,     0.7764,     0.4657,     0.4627,
564:                         0.4569,     0.4232,     0.8514,     0.0674,     0.3227,     0.1055,     0.6690,     0.6313,
565:                         0.9226,     0.5461,     0.4126,     0.2364,     0.6096,     0.7042,     0.3914,     0.0711};

567:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
568:   PetscLogStageRegister("Elliptic Setup",&user->stages[0]);
569:   PetscLogStagePush(user->stages[0]);

571:   /* Create u,y,c,x */
572:   VecCreate(PETSC_COMM_WORLD,&user->u);
573:   VecCreate(PETSC_COMM_WORLD,&user->y);
574:   VecCreate(PETSC_COMM_WORLD,&user->c);
575:   VecSetSizes(user->u,PETSC_DECIDE,user->ndesign);
576:   VecSetFromOptions(user->u);
577:   VecGetLocalSize(user->u,&ysubnlocal);
578:   VecSetSizes(user->y,ysubnlocal*user->ns,user->nstate);
579:   VecSetSizes(user->c,ysubnlocal*user->ns,user->m);
580:   VecSetFromOptions(user->y);
581:   VecSetFromOptions(user->c);

583:   /*
584:      *******************************
585:      Create scatters for x <-> y,u
586:      *******************************

588:      If the state vector y and design vector u are partitioned as
589:      [y_1; y_2; ...; y_np] and [u_1; u_2; ...; u_np] (with np = # of processors),
590:      then the solution vector x is organized as
591:      [y_1; u_1; y_2; u_2; ...; y_np; u_np].
592:      The index sets user->s_is and user->d_is correspond to the indices of the
593:      state and design variables owned by the current processor.
594:   */
595:   VecCreate(PETSC_COMM_WORLD,&user->x);

597:   VecGetOwnershipRange(user->y,&lo,&hi);
598:   VecGetOwnershipRange(user->u,&lo2,&hi2);

600:   ISCreateStride(PETSC_COMM_SELF,hi-lo,lo,1,&is_allstate);
601:   ISCreateStride(PETSC_COMM_SELF,hi-lo,lo+lo2,1,&user->s_is);
602:   ISCreateStride(PETSC_COMM_SELF,hi2-lo2,lo2,1,&is_alldesign);
603:   ISCreateStride(PETSC_COMM_SELF,hi2-lo2,hi+lo2,1,&user->d_is);

605:   VecSetSizes(user->x,hi-lo+hi2-lo2,user->n);
606:   VecSetFromOptions(user->x);

608:   VecScatterCreate(user->x,user->s_is,user->y,is_allstate,&user->state_scatter);
609:   VecScatterCreate(user->x,user->d_is,user->u,is_alldesign,&user->design_scatter);
610:   ISDestroy(&is_alldesign);
611:   ISDestroy(&is_allstate);
612:   /*
613:      *******************************
614:      Create scatter from y to y_1,y_2,...,y_ns
615:      *******************************
616:   */
617:   PetscMalloc1(user->ns,&user->yi_scatter);
618:   VecDuplicate(user->u,&user->suby);
619:   VecDuplicate(user->u,&user->subq);

621:   VecGetOwnershipRange(user->y,&lo2,&hi2);
622:   istart = 0;
623:   for (i=0; i<user->ns; i++) {
624:     VecGetOwnershipRange(user->suby,&lo,&hi);
625:     ISCreateStride(PETSC_COMM_SELF,hi-lo,lo2+istart,1,&is_from_y);
626:     VecScatterCreate(user->y,is_from_y,user->suby,NULL,&user->yi_scatter[i]);
627:     istart = istart + hi-lo;
628:     ISDestroy(&is_from_y);
629:   }
630:   /*
631:      *******************************
632:      Create scatter from d to d_1,d_2,...,d_ns
633:      *******************************
634:   */
635:   VecCreate(PETSC_COMM_WORLD,&user->subd);
636:   VecSetSizes(user->subd,PETSC_DECIDE,user->ndata);
637:   VecSetFromOptions(user->subd);
638:   VecCreate(PETSC_COMM_WORLD,&user->d);
639:   VecGetLocalSize(user->subd,&dsubnlocal);
640:   VecSetSizes(user->d,dsubnlocal*user->ns,user->ndata*user->ns);
641:   VecSetFromOptions(user->d);
642:   PetscMalloc1(user->ns,&user->di_scatter);

644:   VecGetOwnershipRange(user->d,&lo2,&hi2);
645:   istart = 0;
646:   for (i=0; i<user->ns; i++) {
647:     VecGetOwnershipRange(user->subd,&lo,&hi);
648:     ISCreateStride(PETSC_COMM_SELF,hi-lo,lo2+istart,1,&is_from_d);
649:     VecScatterCreate(user->d,is_from_d,user->subd,NULL,&user->di_scatter[i]);
650:     istart = istart + hi-lo;
651:     ISDestroy(&is_from_d);
652:   }

654:   PetscMalloc1(user->mx,&x);
655:   PetscMalloc1(user->mx,&y);
656:   PetscMalloc1(user->mx,&z);

658:   user->ksp_its = 0;
659:   user->ksp_its_initial = 0;

661:   n = user->mx * user->mx * user->mx;
662:   m = 3 * user->mx * user->mx * (user->mx-1);
663:   sqrt_beta = PetscSqrtScalar(user->beta);

665:   VecCreate(PETSC_COMM_WORLD,&XX);
666:   VecCreate(PETSC_COMM_WORLD,&user->q);
667:   VecSetSizes(XX,ysubnlocal,n);
668:   VecSetSizes(user->q,ysubnlocal*user->ns,user->m);
669:   VecSetFromOptions(XX);
670:   VecSetFromOptions(user->q);

672:   VecDuplicate(XX,&YY);
673:   VecDuplicate(XX,&ZZ);
674:   VecDuplicate(XX,&XXwork);
675:   VecDuplicate(XX,&YYwork);
676:   VecDuplicate(XX,&ZZwork);
677:   VecDuplicate(XX,&UTwork);
678:   VecDuplicate(XX,&user->utrue);

680:   /* map for striding q */
681:   VecGetOwnershipRanges(user->q,&ranges);
682:   VecGetOwnershipRanges(user->u,&subranges);

684:   VecGetOwnershipRange(user->q,&lo2,&hi2);
685:   VecGetOwnershipRange(user->u,&lo,&hi);
686:   /* Generate 3D grid, and collect ns (1<=ns<=8) right-hand-side vectors into user->q */
687:   h = 1.0/user->mx;
688:   hinv = user->mx;
689:   neg_hinv = -hinv;

691:   VecGetOwnershipRange(XX,&istart,&iend);
692:   for (linear_index=istart; linear_index<iend; linear_index++) {
693:     i = linear_index % user->mx;
694:     j = ((linear_index-i)/user->mx) % user->mx;
695:     k = ((linear_index-i)/user->mx-j) / user->mx;
696:     vx = h*(i+0.5);
697:     vy = h*(j+0.5);
698:     vz = h*(k+0.5);
699:     VecSetValues(XX,1,&linear_index,&vx,INSERT_VALUES);
700:     VecSetValues(YY,1,&linear_index,&vy,INSERT_VALUES);
701:     VecSetValues(ZZ,1,&linear_index,&vz,INSERT_VALUES);
702:     for (is=0; is<2; is++) {
703:       for (js=0; js<2; js++) {
704:         for (ks=0; ks<2; ks++) {
705:           ls = is*4 + js*2 + ks;
706:           if (ls<user->ns) {
707:             l =ls*n + linear_index;
708:             /* remap */
709:             subindex = l%n;
710:             subvec = l/n;
711:             nrank=0;
712:             while (subindex >= subranges[nrank+1]) nrank++;
713:             offset = subindex - subranges[nrank];
714:             istart=0;
715:             for (kk=0;kk<nrank;kk++) istart+=user->ns*(subranges[kk+1]-subranges[kk]);
716:             istart += (subranges[nrank+1]-subranges[nrank])*subvec;
717:             l = istart+offset;
718:             v = 100*PetscSinScalar(2*PETSC_PI*(vx+0.25*is))*PetscSinScalar(2*PETSC_PI*(vy+0.25*js))*PetscSinScalar(2*PETSC_PI*(vz+0.25*ks));
719:             VecSetValues(user->q,1,&l,&v,INSERT_VALUES);
720:           }
721:         }
722:       }
723:     }
724:   }

726:   VecAssemblyBegin(XX);
727:   VecAssemblyEnd(XX);
728:   VecAssemblyBegin(YY);
729:   VecAssemblyEnd(YY);
730:   VecAssemblyBegin(ZZ);
731:   VecAssemblyEnd(ZZ);
732:   VecAssemblyBegin(user->q);
733:   VecAssemblyEnd(user->q);

735:   /* Compute true parameter function
736:      ut = exp(-((x-0.25)^2+(y-0.25)^2+(z-0.25)^2)/0.05) - exp((x-0.75)^2-(y-0.75)^2-(z-0.75))^2/0.05) */
737:   VecCopy(XX,XXwork);
738:   VecCopy(YY,YYwork);
739:   VecCopy(ZZ,ZZwork);

741:   VecShift(XXwork,-0.25);
742:   VecShift(YYwork,-0.25);
743:   VecShift(ZZwork,-0.25);

745:   VecPointwiseMult(XXwork,XXwork,XXwork);
746:   VecPointwiseMult(YYwork,YYwork,YYwork);
747:   VecPointwiseMult(ZZwork,ZZwork,ZZwork);

749:   VecCopy(XXwork,UTwork);
750:   VecAXPY(UTwork,1.0,YYwork);
751:   VecAXPY(UTwork,1.0,ZZwork);
752:   VecScale(UTwork,-20.0);
753:   VecExp(UTwork);
754:   VecCopy(UTwork,user->utrue);

756:   VecCopy(XX,XXwork);
757:   VecCopy(YY,YYwork);
758:   VecCopy(ZZ,ZZwork);

760:   VecShift(XXwork,-0.75);
761:   VecShift(YYwork,-0.75);
762:   VecShift(ZZwork,-0.75);

764:   VecPointwiseMult(XXwork,XXwork,XXwork);
765:   VecPointwiseMult(YYwork,YYwork,YYwork);
766:   VecPointwiseMult(ZZwork,ZZwork,ZZwork);

768:   VecCopy(XXwork,UTwork);
769:   VecAXPY(UTwork,1.0,YYwork);
770:   VecAXPY(UTwork,1.0,ZZwork);
771:   VecScale(UTwork,-20.0);
772:   VecExp(UTwork);

774:   VecAXPY(user->utrue,-1.0,UTwork);

776:   VecDestroy(&XX);
777:   VecDestroy(&YY);
778:   VecDestroy(&ZZ);
779:   VecDestroy(&XXwork);
780:   VecDestroy(&YYwork);
781:   VecDestroy(&ZZwork);
782:   VecDestroy(&UTwork);

784:   /* Initial guess and reference model */
785:   VecDuplicate(user->utrue,&user->ur);
786:   VecSum(user->utrue,&meanut);
787:   meanut = meanut / n;
788:   VecSet(user->ur,meanut);
789:   VecCopy(user->ur,user->u);

791:   /* Generate Grad matrix */
792:   MatCreate(PETSC_COMM_WORLD,&user->Grad);
793:   MatSetSizes(user->Grad,PETSC_DECIDE,ysubnlocal,m,n);
794:   MatSetFromOptions(user->Grad);
795:   MatMPIAIJSetPreallocation(user->Grad,2,NULL,2,NULL);
796:   MatSeqAIJSetPreallocation(user->Grad,2,NULL);
797:   MatGetOwnershipRange(user->Grad,&istart,&iend);

799:   for (i=istart; i<iend; i++) {
800:     if (i<m/3) {
801:       iblock = i / (user->mx-1);
802:       j = iblock*user->mx + (i % (user->mx-1));
803:       MatSetValues(user->Grad,1,&i,1,&j,&neg_hinv,INSERT_VALUES);
804:       j = j+1;
805:       MatSetValues(user->Grad,1,&i,1,&j,&hinv,INSERT_VALUES);
806:     }
807:     if (i>=m/3 && i<2*m/3) {
808:       iblock = (i-m/3) / (user->mx*(user->mx-1));
809:       j = iblock*user->mx*user->mx + ((i-m/3) % (user->mx*(user->mx-1)));
810:       MatSetValues(user->Grad,1,&i,1,&j,&neg_hinv,INSERT_VALUES);
811:       j = j + user->mx;
812:       MatSetValues(user->Grad,1,&i,1,&j,&hinv,INSERT_VALUES);
813:     }
814:     if (i>=2*m/3) {
815:       j = i-2*m/3;
816:       MatSetValues(user->Grad,1,&i,1,&j,&neg_hinv,INSERT_VALUES);
817:       j = j + user->mx*user->mx;
818:       MatSetValues(user->Grad,1,&i,1,&j,&hinv,INSERT_VALUES);
819:     }
820:   }

822:   MatAssemblyBegin(user->Grad,MAT_FINAL_ASSEMBLY);
823:   MatAssemblyEnd(user->Grad,MAT_FINAL_ASSEMBLY);

825:   /* Generate arithmetic averaging matrix Av */
826:   MatCreate(PETSC_COMM_WORLD,&user->Av);
827:   MatSetSizes(user->Av,PETSC_DECIDE,ysubnlocal,m,n);
828:   MatSetFromOptions(user->Av);
829:   MatMPIAIJSetPreallocation(user->Av,2,NULL,2,NULL);
830:   MatSeqAIJSetPreallocation(user->Av,2,NULL);
831:   MatGetOwnershipRange(user->Av,&istart,&iend);

833:   for (i=istart; i<iend; i++) {
834:     if (i<m/3) {
835:       iblock = i / (user->mx-1);
836:       j = iblock*user->mx + (i % (user->mx-1));
837:       MatSetValues(user->Av,1,&i,1,&j,&half,INSERT_VALUES);
838:       j = j+1;
839:       MatSetValues(user->Av,1,&i,1,&j,&half,INSERT_VALUES);
840:     }
841:     if (i>=m/3 && i<2*m/3) {
842:       iblock = (i-m/3) / (user->mx*(user->mx-1));
843:       j = iblock*user->mx*user->mx + ((i-m/3) % (user->mx*(user->mx-1)));
844:       MatSetValues(user->Av,1,&i,1,&j,&half,INSERT_VALUES);
845:       j = j + user->mx;
846:       MatSetValues(user->Av,1,&i,1,&j,&half,INSERT_VALUES);
847:     }
848:     if (i>=2*m/3) {
849:       j = i-2*m/3;
850:       MatSetValues(user->Av,1,&i,1,&j,&half,INSERT_VALUES);
851:       j = j + user->mx*user->mx;
852:       MatSetValues(user->Av,1,&i,1,&j,&half,INSERT_VALUES);
853:     }
854:   }

856:   MatAssemblyBegin(user->Av,MAT_FINAL_ASSEMBLY);
857:   MatAssemblyEnd(user->Av,MAT_FINAL_ASSEMBLY);

859:   MatCreate(PETSC_COMM_WORLD,&user->L);
860:   MatSetSizes(user->L,PETSC_DECIDE,ysubnlocal,m+n,n);
861:   MatSetFromOptions(user->L);
862:   MatMPIAIJSetPreallocation(user->L,2,NULL,2,NULL);
863:   MatSeqAIJSetPreallocation(user->L,2,NULL);
864:   MatGetOwnershipRange(user->L,&istart,&iend);

866:   for (i=istart; i<iend; i++) {
867:     if (i<m/3) {
868:       iblock = i / (user->mx-1);
869:       j = iblock*user->mx + (i % (user->mx-1));
870:       MatSetValues(user->L,1,&i,1,&j,&neg_hinv,INSERT_VALUES);
871:       j = j+1;
872:       MatSetValues(user->L,1,&i,1,&j,&hinv,INSERT_VALUES);
873:     }
874:     if (i>=m/3 && i<2*m/3) {
875:       iblock = (i-m/3) / (user->mx*(user->mx-1));
876:       j = iblock*user->mx*user->mx + ((i-m/3) % (user->mx*(user->mx-1)));
877:       MatSetValues(user->L,1,&i,1,&j,&neg_hinv,INSERT_VALUES);
878:       j = j + user->mx;
879:       MatSetValues(user->L,1,&i,1,&j,&hinv,INSERT_VALUES);
880:     }
881:     if (i>=2*m/3 && i<m) {
882:       j = i-2*m/3;
883:       MatSetValues(user->L,1,&i,1,&j,&neg_hinv,INSERT_VALUES);
884:       j = j + user->mx*user->mx;
885:       MatSetValues(user->L,1,&i,1,&j,&hinv,INSERT_VALUES);
886:     }
887:     if (i>=m) {
888:       j = i - m;
889:       MatSetValues(user->L,1,&i,1,&j,&sqrt_beta,INSERT_VALUES);
890:     }
891:   }
892:   MatAssemblyBegin(user->L,MAT_FINAL_ASSEMBLY);
893:   MatAssemblyEnd(user->L,MAT_FINAL_ASSEMBLY);
894:   MatScale(user->L,PetscPowScalar(h,1.5));

896:   /* Generate Div matrix */
897:   if (!user->use_ptap) {
898:     /* Generate Div matrix */
899:     MatCreate(PETSC_COMM_WORLD,&user->Div);
900:     MatSetSizes(user->Div,ysubnlocal,PETSC_DECIDE,n,m);
901:     MatSetFromOptions(user->Div);
902:     MatMPIAIJSetPreallocation(user->Div,4,NULL,4,NULL);
903:     MatSeqAIJSetPreallocation(user->Div,6,NULL);
904:     MatGetOwnershipRange(user->Grad,&istart,&iend);

906:     for (i=istart; i<iend; i++) {
907:       if (i<m/3) {
908:         iblock = i / (user->mx-1);
909:         j = iblock*user->mx + (i % (user->mx-1));
910:         MatSetValues(user->Div,1,&j,1,&i,&neg_hinv,INSERT_VALUES);
911:         j = j+1;
912:         MatSetValues(user->Div,1,&j,1,&i,&hinv,INSERT_VALUES);
913:       }
914:       if (i>=m/3 && i<2*m/3) {
915:         iblock = (i-m/3) / (user->mx*(user->mx-1));
916:         j = iblock*user->mx*user->mx + ((i-m/3) % (user->mx*(user->mx-1)));
917:         MatSetValues(user->Div,1,&j,1,&i,&neg_hinv,INSERT_VALUES);
918:         j = j + user->mx;
919:         MatSetValues(user->Div,1,&j,1,&i,&hinv,INSERT_VALUES);
920:       }
921:       if (i>=2*m/3) {
922:         j = i-2*m/3;
923:         MatSetValues(user->Div,1,&j,1,&i,&neg_hinv,INSERT_VALUES);
924:         j = j + user->mx*user->mx;
925:         MatSetValues(user->Div,1,&j,1,&i,&hinv,INSERT_VALUES);
926:       }
927:     }

929:     MatAssemblyBegin(user->Div,MAT_FINAL_ASSEMBLY);
930:     MatAssemblyEnd(user->Div,MAT_FINAL_ASSEMBLY);
931:     MatDuplicate(user->Div,MAT_SHARE_NONZERO_PATTERN,&user->Divwork);
932:   } else {
933:     MatCreate(PETSC_COMM_WORLD,&user->Diag);
934:     MatSetSizes(user->Diag,PETSC_DECIDE,PETSC_DECIDE,m,m);
935:     MatSetFromOptions(user->Diag);
936:     MatMPIAIJSetPreallocation(user->Diag,1,NULL,0,NULL);
937:     MatSeqAIJSetPreallocation(user->Diag,1,NULL);
938:   }

940:   /* Build work vectors and matrices */
941:   VecCreate(PETSC_COMM_WORLD,&user->S);
942:   VecSetSizes(user->S, PETSC_DECIDE, m);
943:   VecSetFromOptions(user->S);

945:   VecCreate(PETSC_COMM_WORLD,&user->lwork);
946:   VecSetSizes(user->lwork,PETSC_DECIDE,m+user->mx*user->mx*user->mx);
947:   VecSetFromOptions(user->lwork);

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

951:   VecDuplicate(user->S,&user->Swork);
952:   VecDuplicate(user->S,&user->Sdiag);
953:   VecDuplicate(user->S,&user->Av_u);
954:   VecDuplicate(user->S,&user->Twork);
955:   VecDuplicate(user->y,&user->ywork);
956:   VecDuplicate(user->u,&user->Ywork);
957:   VecDuplicate(user->u,&user->uwork);
958:   VecDuplicate(user->u,&user->js_diag);
959:   VecDuplicate(user->c,&user->cwork);
960:   VecDuplicate(user->d,&user->dwork);

962:   /* Create a matrix-free shell user->Jd for computing B*x */
963:   MatCreateShell(PETSC_COMM_WORLD,ysubnlocal*user->ns,ysubnlocal,user->nstate,user->ndesign,user,&user->Jd);
964:   MatShellSetOperation(user->Jd,MATOP_MULT,(void(*)(void))DesignMatMult);
965:   MatShellSetOperation(user->Jd,MATOP_MULT_TRANSPOSE,(void(*)(void))DesignMatMultTranspose);

967:   /* Compute true state function ytrue given utrue */
968:   VecDuplicate(user->y,&user->ytrue);

970:   /* First compute Av_u = Av*exp(-u) */
971:   VecSet(user->uwork, 0);
972:   VecAXPY(user->uwork,-1.0,user->utrue); /* Note: user->utrue */
973:   VecExp(user->uwork);
974:   MatMult(user->Av,user->uwork,user->Av_u);

976:   /* Next form DSG = Div*S*Grad */
977:   VecCopy(user->Av_u,user->Swork);
978:   VecReciprocal(user->Swork);
979:   if (user->use_ptap) {
980:     MatDiagonalSet(user->Diag,user->Swork,INSERT_VALUES);
981:     MatPtAP(user->Diag,user->Grad,MAT_INITIAL_MATRIX,1.0,&user->DSG);
982:   } else {
983:     MatCopy(user->Div,user->Divwork,SAME_NONZERO_PATTERN);
984:     MatDiagonalScale(user->Divwork,NULL,user->Swork);

986:     MatMatMult(user->Divwork,user->Grad,MAT_INITIAL_MATRIX,1.0,&user->DSG);
987:   }

989:   MatSetOption(user->DSG,MAT_SYMMETRIC,PETSC_TRUE);
990:   MatSetOption(user->DSG,MAT_SYMMETRY_ETERNAL,PETSC_TRUE);

992:   if (user->use_lrc == PETSC_TRUE) {
993:     v=PetscSqrtReal(1.0 /user->ndesign);
994:     PetscMalloc1(user->ndesign,&user->ones);

996:     for (i=0;i<user->ndesign;i++) {
997:       user->ones[i]=v;
998:     }
999:     MatCreateDense(PETSC_COMM_WORLD,ysubnlocal,PETSC_DECIDE,user->ndesign,1,user->ones,&user->Ones);
1000:     MatAssemblyBegin(user->Ones, MAT_FINAL_ASSEMBLY);
1001:     MatAssemblyEnd(user->Ones, MAT_FINAL_ASSEMBLY);
1002:     MatCreateLRC(user->DSG,user->Ones,NULL,user->Ones,&user->JsBlock);
1003:     MatSetUp(user->JsBlock);
1004:   } else {
1005:     /* Create matrix-free shell user->Js for computing (A + h^3*e*e^T)*x */
1006:     MatCreateShell(PETSC_COMM_WORLD,ysubnlocal,ysubnlocal,user->ndesign,user->ndesign,user,&user->JsBlock);
1007:     MatShellSetOperation(user->JsBlock,MATOP_MULT,(void(*)(void))StateBlockMatMult);
1008:     MatShellSetOperation(user->JsBlock,MATOP_MULT_TRANSPOSE,(void(*)(void))StateBlockMatMult);
1009:   }
1010:   MatSetOption(user->JsBlock,MAT_SYMMETRIC,PETSC_TRUE);
1011:   MatSetOption(user->JsBlock,MAT_SYMMETRY_ETERNAL,PETSC_TRUE);
1012:   MatCreateShell(PETSC_COMM_WORLD,ysubnlocal*user->ns,ysubnlocal*user->ns,user->nstate,user->nstate,user,&user->Js);
1013:   MatShellSetOperation(user->Js,MATOP_MULT,(void(*)(void))StateMatMult);
1014:   MatShellSetOperation(user->Js,MATOP_MULT_TRANSPOSE,(void(*)(void))StateMatMult);
1015:   MatSetOption(user->Js,MAT_SYMMETRIC,PETSC_TRUE);
1016:   MatSetOption(user->Js,MAT_SYMMETRY_ETERNAL,PETSC_TRUE);

1018:   MatCreateShell(PETSC_COMM_WORLD,ysubnlocal*user->ns,ysubnlocal*user->ns,user->nstate,user->nstate,user,&user->JsInv);
1019:   MatShellSetOperation(user->JsInv,MATOP_MULT,(void(*)(void))StateInvMatMult);
1020:   MatShellSetOperation(user->JsInv,MATOP_MULT_TRANSPOSE,(void(*)(void))StateInvMatMult);
1021:   MatSetOption(user->JsInv,MAT_SYMMETRIC,PETSC_TRUE);
1022:   MatSetOption(user->JsInv,MAT_SYMMETRY_ETERNAL,PETSC_TRUE);

1024:   MatSetOption(user->DSG,MAT_SYMMETRIC,PETSC_TRUE);
1025:   MatSetOption(user->DSG,MAT_SYMMETRY_ETERNAL,PETSC_TRUE);
1026:   /* Now solve for ytrue */
1027:   KSPCreate(PETSC_COMM_WORLD,&user->solver);
1028:   KSPSetFromOptions(user->solver);

1030:   KSPSetOperators(user->solver,user->JsBlock,user->DSG);

1032:   MatMult(user->JsInv,user->q,user->ytrue);
1033:   /* First compute Av_u = Av*exp(-u) */
1034:   VecSet(user->uwork,0);
1035:   VecAXPY(user->uwork,-1.0,user->u); /* Note: user->u */
1036:   VecExp(user->uwork);
1037:   MatMult(user->Av,user->uwork,user->Av_u);

1039:   /* Next update DSG = Div*S*Grad  with user->u */
1040:   VecCopy(user->Av_u,user->Swork);
1041:   VecReciprocal(user->Swork);
1042:   if (user->use_ptap) {
1043:     MatDiagonalSet(user->Diag,user->Swork,INSERT_VALUES);
1044:     MatPtAP(user->Diag,user->Grad,MAT_REUSE_MATRIX,1.0,&user->DSG);
1045:   } else {
1046:     MatCopy(user->Div,user->Divwork,SAME_NONZERO_PATTERN);
1047:     MatDiagonalScale(user->Divwork,NULL,user->Av_u);
1048:     MatProductNumeric(user->DSG);
1049:   }

1051:   /* Now solve for y */

1053:   MatMult(user->JsInv,user->q,user->y);

1055:   user->ksp_its_initial = user->ksp_its;
1056:   user->ksp_its = 0;
1057:   /* Construct projection matrix Q (blocks) */
1058:   MatCreate(PETSC_COMM_WORLD,&user->Q);
1059:   MatSetSizes(user->Q,dsubnlocal,ysubnlocal,user->ndata,user->ndesign);
1060:   MatSetFromOptions(user->Q);
1061:   MatMPIAIJSetPreallocation(user->Q,8,NULL,8,NULL);
1062:   MatSeqAIJSetPreallocation(user->Q,8,NULL);

1064:   for (i=0; i<user->mx; i++) {
1065:     x[i] = h*(i+0.5);
1066:     y[i] = h*(i+0.5);
1067:     z[i] = h*(i+0.5);
1068:   }
1069:   MatGetOwnershipRange(user->Q,&istart,&iend);

1071:   nx = user->mx; ny = user->mx; nz = user->mx;
1072:   for (i=istart; i<iend; i++) {

1074:     xri = xr[i];
1075:     im = 0;
1076:     xim = x[im];
1077:     while (xri>xim && im<nx) {
1078:       im = im+1;
1079:       xim = x[im];
1080:     }
1081:     indx1 = im-1;
1082:     indx2 = im;
1083:     dx1 = xri - x[indx1];
1084:     dx2 = x[indx2] - xri;

1086:     yri = yr[i];
1087:     im = 0;
1088:     yim = y[im];
1089:     while (yri>yim && im<ny) {
1090:       im = im+1;
1091:       yim = y[im];
1092:     }
1093:     indy1 = im-1;
1094:     indy2 = im;
1095:     dy1 = yri - y[indy1];
1096:     dy2 = y[indy2] - yri;

1098:     zri = zr[i];
1099:     im = 0;
1100:     zim = z[im];
1101:     while (zri>zim && im<nz) {
1102:       im = im+1;
1103:       zim = z[im];
1104:     }
1105:     indz1 = im-1;
1106:     indz2 = im;
1107:     dz1 = zri - z[indz1];
1108:     dz2 = z[indz2] - zri;

1110:     Dx = x[indx2] - x[indx1];
1111:     Dy = y[indy2] - y[indy1];
1112:     Dz = z[indz2] - z[indz1];

1114:     j = indx1 + indy1*nx + indz1*nx*ny;
1115:     v = (1-dx1/Dx)*(1-dy1/Dy)*(1-dz1/Dz);
1116:     MatSetValues(user->Q,1,&i,1,&j,&v,INSERT_VALUES);

1118:     j = indx1 + indy1*nx + indz2*nx*ny;
1119:     v = (1-dx1/Dx)*(1-dy1/Dy)*(1-dz2/Dz);
1120:     MatSetValues(user->Q,1,&i,1,&j,&v,INSERT_VALUES);

1122:     j = indx1 + indy2*nx + indz1*nx*ny;
1123:     v = (1-dx1/Dx)*(1-dy2/Dy)*(1-dz1/Dz);
1124:     MatSetValues(user->Q,1,&i,1,&j,&v,INSERT_VALUES);

1126:     j = indx1 + indy2*nx + indz2*nx*ny;
1127:     v = (1-dx1/Dx)*(1-dy2/Dy)*(1-dz2/Dz);
1128:     MatSetValues(user->Q,1,&i,1,&j,&v,INSERT_VALUES);

1130:     j = indx2 + indy1*nx + indz1*nx*ny;
1131:     v = (1-dx2/Dx)*(1-dy1/Dy)*(1-dz1/Dz);
1132:     MatSetValues(user->Q,1,&i,1,&j,&v,INSERT_VALUES);

1134:     j = indx2 + indy1*nx + indz2*nx*ny;
1135:     v = (1-dx2/Dx)*(1-dy1/Dy)*(1-dz2/Dz);
1136:     MatSetValues(user->Q,1,&i,1,&j,&v,INSERT_VALUES);

1138:     j = indx2 + indy2*nx + indz1*nx*ny;
1139:     v = (1-dx2/Dx)*(1-dy2/Dy)*(1-dz1/Dz);
1140:     MatSetValues(user->Q,1,&i,1,&j,&v,INSERT_VALUES);

1142:     j = indx2 + indy2*nx + indz2*nx*ny;
1143:     v = (1-dx2/Dx)*(1-dy2/Dy)*(1-dz2/Dz);
1144:     MatSetValues(user->Q,1,&i,1,&j,&v,INSERT_VALUES);
1145:   }

1147:   MatAssemblyBegin(user->Q,MAT_FINAL_ASSEMBLY);
1148:   MatAssemblyEnd(user->Q,MAT_FINAL_ASSEMBLY);
1149:   /* Create MQ (composed of blocks of Q */
1150:   MatCreateShell(PETSC_COMM_WORLD,dsubnlocal*user->ns,PETSC_DECIDE,user->ndata*user->ns,user->nstate,user,&user->MQ);
1151:   MatShellSetOperation(user->MQ,MATOP_MULT,(void(*)(void))QMatMult);
1152:   MatShellSetOperation(user->MQ,MATOP_MULT_TRANSPOSE,(void(*)(void))QMatMultTranspose);

1154:   /* Add noise to the measurement data */
1155:   VecSet(user->ywork,1.0);
1156:   VecAYPX(user->ywork,user->noise,user->ytrue);
1157:   MatMult(user->MQ,user->ywork,user->d);

1159:   /* Now that initial conditions have been set, let the user pass tolerance options to the KSP solver */
1160:   PetscFree(x);
1161:   PetscFree(y);
1162:   PetscFree(z);
1163:   PetscLogStagePop();
1164:   return 0;
1165: }

1167: PetscErrorCode EllipticDestroy(AppCtx *user)
1168: {
1169:   PetscInt       i;

1171:   MatDestroy(&user->DSG);
1172:   KSPDestroy(&user->solver);
1173:   MatDestroy(&user->Q);
1174:   MatDestroy(&user->MQ);
1175:   if (!user->use_ptap) {
1176:     MatDestroy(&user->Div);
1177:     MatDestroy(&user->Divwork);
1178:   } else {
1179:     MatDestroy(&user->Diag);
1180:   }
1181:   if (user->use_lrc) {
1182:     MatDestroy(&user->Ones);
1183:   }

1185:   MatDestroy(&user->Grad);
1186:   MatDestroy(&user->Av);
1187:   MatDestroy(&user->Avwork);
1188:   MatDestroy(&user->L);
1189:   MatDestroy(&user->Js);
1190:   MatDestroy(&user->Jd);
1191:   MatDestroy(&user->JsBlock);
1192:   MatDestroy(&user->JsInv);

1194:   VecDestroy(&user->x);
1195:   VecDestroy(&user->u);
1196:   VecDestroy(&user->uwork);
1197:   VecDestroy(&user->utrue);
1198:   VecDestroy(&user->y);
1199:   VecDestroy(&user->ywork);
1200:   VecDestroy(&user->ytrue);
1201:   VecDestroy(&user->c);
1202:   VecDestroy(&user->cwork);
1203:   VecDestroy(&user->ur);
1204:   VecDestroy(&user->q);
1205:   VecDestroy(&user->d);
1206:   VecDestroy(&user->dwork);
1207:   VecDestroy(&user->lwork);
1208:   VecDestroy(&user->S);
1209:   VecDestroy(&user->Swork);
1210:   VecDestroy(&user->Sdiag);
1211:   VecDestroy(&user->Ywork);
1212:   VecDestroy(&user->Twork);
1213:   VecDestroy(&user->Av_u);
1214:   VecDestroy(&user->js_diag);
1215:   ISDestroy(&user->s_is);
1216:   ISDestroy(&user->d_is);
1217:   VecDestroy(&user->suby);
1218:   VecDestroy(&user->subd);
1219:   VecDestroy(&user->subq);
1220:   VecScatterDestroy(&user->state_scatter);
1221:   VecScatterDestroy(&user->design_scatter);
1222:   for (i=0;i<user->ns;i++) {
1223:     VecScatterDestroy(&user->yi_scatter[i]);
1224:     VecScatterDestroy(&user->di_scatter[i]);
1225:   }
1226:   PetscFree(user->yi_scatter);
1227:   PetscFree(user->di_scatter);
1228:   if (user->use_lrc) {
1229:     PetscFree(user->ones);
1230:     MatDestroy(&user->Ones);
1231:   }
1232:   return 0;
1233: }

1235: PetscErrorCode EllipticMonitor(Tao tao, void *ptr)
1236: {
1237:   Vec            X;
1238:   PetscReal      unorm,ynorm;
1239:   AppCtx         *user = (AppCtx*)ptr;

1241:   TaoGetSolution(tao,&X);
1242:   Scatter(X,user->ywork,user->state_scatter,user->uwork,user->design_scatter);
1243:   VecAXPY(user->ywork,-1.0,user->ytrue);
1244:   VecAXPY(user->uwork,-1.0,user->utrue);
1245:   VecNorm(user->uwork,NORM_2,&unorm);
1246:   VecNorm(user->ywork,NORM_2,&ynorm);
1247:   PetscPrintf(MPI_COMM_WORLD, "||u-ut||=%g ||y-yt||=%g\n",(double)unorm,(double)ynorm);
1248:   return 0;
1249: }

1251: /*TEST

1253:    build:
1254:       requires: !complex

1256:    test:
1257:       args: -tao_cmonitor -ns 1 -tao_type lcl -tao_gatol 1.e-3 -tao_max_it 11
1258:       requires: !single

1260:    test:
1261:       suffix: 2
1262:       args: -tao_cmonitor -tao_type lcl -tao_max_it 11 -use_ptap -use_lrc -ns 1 -tao_gatol 1.e-3
1263:       requires: !single

1265: TEST*/