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*/