Actual source code: elliptic.c
1: #include <petsc/private/taoimpl.h>
3: /*T
4: Concepts: TAO^Solving a system of nonlinear equations, nonlinear least squares
5: Routines: TaoCreate();
6: Routines: TaoSetType();
7: Routines: TaoSetInitialVector();
8: Routines: TaoSetObjectiveRoutine();
9: Routines: TaoSetGradientRoutine();
10: Routines: TaoSetConstraintsRoutine();
11: Routines: TaoSetJacobianStateRoutine();
12: Routines: TaoSetJacobianDesignRoutine();
13: Routines: TaoSetStateDesignIS();
14: Routines: TaoSetFromOptions();
15: Routines: TaoSolve();
16: Routines: TaoDestroy();
17: Processors: n
18: T*/
20: typedef struct {
21: PetscInt n; /* Number of total variables */
22: PetscInt m; /* Number of constraints */
23: PetscInt nstate;
24: PetscInt ndesign;
25: PetscInt mx; /* grid points in each direction */
26: PetscInt ns; /* Number of data samples (1<=ns<=8)
27: Currently only ns=1 is supported */
28: PetscInt ndata; /* Number of data points per sample */
29: IS s_is;
30: IS d_is;
32: VecScatter state_scatter;
33: VecScatter design_scatter;
34: VecScatter *yi_scatter, *di_scatter;
35: Vec suby,subq,subd;
36: Mat Js,Jd,JsPrec,JsInv,JsBlock;
38: PetscReal alpha; /* Regularization parameter */
39: PetscReal beta; /* Weight attributed to ||u||^2 in regularization functional */
40: PetscReal noise; /* Amount of noise to add to data */
41: PetscReal *ones;
42: Mat Q;
43: Mat MQ;
44: Mat L;
46: Mat Grad;
47: Mat Av,Avwork;
48: Mat Div, Divwork;
49: Mat DSG;
50: Mat Diag,Ones;
52: Vec q;
53: Vec ur; /* reference */
55: Vec d;
56: Vec dwork;
58: Vec x; /* super vec of y,u */
60: Vec y; /* state variables */
61: Vec ywork;
63: Vec ytrue;
65: Vec u; /* design variables */
66: Vec uwork;
68: Vec utrue;
70: Vec js_diag;
72: Vec c; /* constraint vector */
73: Vec cwork;
75: Vec lwork;
76: Vec S;
77: Vec Swork,Twork,Sdiag,Ywork;
78: Vec Av_u;
80: KSP solver;
81: PC prec;
83: PetscReal tola,tolb,tolc,told;
84: PetscInt ksp_its;
85: PetscInt ksp_its_initial;
86: PetscLogStage stages[10];
87: PetscBool use_ptap;
88: PetscBool use_lrc;
89: } AppCtx;
91: PetscErrorCode FormFunction(Tao, Vec, PetscReal*, void*);
92: PetscErrorCode FormGradient(Tao, Vec, Vec, void*);
93: PetscErrorCode FormFunctionGradient(Tao, Vec, PetscReal*, Vec, void*);
94: PetscErrorCode FormJacobianState(Tao, Vec, Mat, Mat, Mat, void*);
95: PetscErrorCode FormJacobianDesign(Tao, Vec, Mat,void*);
96: PetscErrorCode FormConstraints(Tao, Vec, Vec, void*);
97: PetscErrorCode FormHessian(Tao, Vec, Mat, Mat, void*);
98: PetscErrorCode Gather(Vec, Vec, VecScatter, Vec, VecScatter);
99: PetscErrorCode Scatter(Vec, Vec, VecScatter, Vec, VecScatter);
100: PetscErrorCode EllipticInitialize(AppCtx*);
101: PetscErrorCode EllipticDestroy(AppCtx*);
102: PetscErrorCode EllipticMonitor(Tao, void*);
104: PetscErrorCode StateBlockMatMult(Mat,Vec,Vec);
105: PetscErrorCode StateMatMult(Mat,Vec,Vec);
107: PetscErrorCode StateInvMatMult(Mat,Vec,Vec);
108: PetscErrorCode DesignMatMult(Mat,Vec,Vec);
109: PetscErrorCode DesignMatMultTranspose(Mat,Vec,Vec);
111: PetscErrorCode QMatMult(Mat,Vec,Vec);
112: PetscErrorCode QMatMultTranspose(Mat,Vec,Vec);
114: static char help[]="";
116: int main(int argc, char **argv)
117: {
118: PetscErrorCode ierr;
119: Vec x0;
120: Tao tao;
121: AppCtx user;
122: PetscInt ntests = 1;
123: PetscInt i;
125: PetscInitialize(&argc, &argv, (char*)0,help);if (ierr) return ierr;
126: user.mx = 8;
127: PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"elliptic example",NULL);
128: PetscOptionsInt("-mx","Number of grid points in each direction","",user.mx,&user.mx,NULL);
129: user.ns = 6;
130: PetscOptionsInt("-ns","Number of data samples (1<=ns<=8)","",user.ns,&user.ns,NULL);
131: user.ndata = 64;
132: PetscOptionsInt("-ndata","Numbers of data points per sample","",user.ndata,&user.ndata,NULL);
133: user.alpha = 0.1;
134: PetscOptionsReal("-alpha","Regularization parameter","",user.alpha,&user.alpha,NULL);
135: user.beta = 0.00001;
136: PetscOptionsReal("-beta","Weight attributed to ||u||^2 in regularization functional","",user.beta,&user.beta,NULL);
137: user.noise = 0.01;
138: PetscOptionsReal("-noise","Amount of noise to add to data","",user.noise,&user.noise,NULL);
140: user.use_ptap = PETSC_FALSE;
141: PetscOptionsBool("-use_ptap","Use ptap matrix for DSG","",user.use_ptap,&user.use_ptap,NULL);
142: user.use_lrc = PETSC_FALSE;
143: PetscOptionsBool("-use_lrc","Use lrc matrix for Js","",user.use_lrc,&user.use_lrc,NULL);
144: PetscOptionsInt("-ntests","Number of times to repeat TaoSolve","",ntests,&ntests,NULL);
145: PetscOptionsEnd();
147: user.m = user.ns*user.mx*user.mx*user.mx; /* number of constraints */
148: user.nstate = user.m;
149: user.ndesign = user.mx*user.mx*user.mx;
150: user.n = user.nstate + user.ndesign; /* number of variables */
152: /* Create TAO solver and set desired solution method */
153: TaoCreate(PETSC_COMM_WORLD,&tao);
154: TaoSetType(tao,TAOLCL);
156: /* Set up initial vectors and matrices */
157: EllipticInitialize(&user);
159: Gather(user.x,user.y,user.state_scatter,user.u,user.design_scatter);
160: VecDuplicate(user.x,&x0);
161: VecCopy(user.x,x0);
163: /* Set solution vector with an initial guess */
164: TaoSetInitialVector(tao,user.x);
165: TaoSetObjectiveRoutine(tao, FormFunction, &user);
166: TaoSetGradientRoutine(tao, FormGradient, &user);
167: TaoSetConstraintsRoutine(tao, user.c, FormConstraints, &user);
169: TaoSetJacobianStateRoutine(tao, user.Js, NULL, user.JsInv, FormJacobianState, &user);
170: TaoSetJacobianDesignRoutine(tao, user.Jd, FormJacobianDesign, &user);
172: TaoSetStateDesignIS(tao,user.s_is,user.d_is);
173: TaoSetFromOptions(tao);
175: /* SOLVE THE APPLICATION */
176: PetscLogStageRegister("Trials",&user.stages[1]);
177: PetscLogStagePush(user.stages[1]);
178: for (i=0; i<ntests; i++) {
179: TaoSolve(tao);
180: PetscPrintf(PETSC_COMM_WORLD,"KSP Iterations = %D\n",user.ksp_its);
181: VecCopy(x0,user.x);
182: }
183: PetscLogStagePop();
184: PetscBarrier((PetscObject)user.x);
185: PetscPrintf(PETSC_COMM_WORLD,"KSP iterations within initialization: ");
186: PetscPrintf(PETSC_COMM_WORLD,"%D\n",user.ksp_its_initial);
188: TaoDestroy(&tao);
189: VecDestroy(&x0);
190: EllipticDestroy(&user);
191: PetscFinalize();
192: return ierr;
193: }
194: /* ------------------------------------------------------------------- */
195: /*
196: dwork = Qy - d
197: lwork = L*(u-ur)
198: f = 1/2 * (dwork.dwork + alpha*lwork.lwork)
199: */
200: PetscErrorCode FormFunction(Tao tao,Vec X,PetscReal *f,void *ptr)
201: {
203: PetscReal d1=0,d2=0;
204: AppCtx *user = (AppCtx*)ptr;
207: Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);
208: MatMult(user->MQ,user->y,user->dwork);
209: VecAXPY(user->dwork,-1.0,user->d);
210: VecDot(user->dwork,user->dwork,&d1);
211: VecWAXPY(user->uwork,-1.0,user->ur,user->u);
212: MatMult(user->L,user->uwork,user->lwork);
213: VecDot(user->lwork,user->lwork,&d2);
214: *f = 0.5 * (d1 + user->alpha*d2);
215: return(0);
216: }
218: /* ------------------------------------------------------------------- */
219: /*
220: state: g_s = Q' *(Qy - d)
221: design: g_d = alpha*L'*L*(u-ur)
222: */
223: PetscErrorCode FormGradient(Tao tao,Vec X,Vec G,void *ptr)
224: {
226: AppCtx *user = (AppCtx*)ptr;
229: Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);
230: MatMult(user->MQ,user->y,user->dwork);
231: VecAXPY(user->dwork,-1.0,user->d);
232: MatMultTranspose(user->MQ,user->dwork,user->ywork);
233: VecWAXPY(user->uwork,-1.0,user->ur,user->u);
234: MatMult(user->L,user->uwork,user->lwork);
235: MatMultTranspose(user->L,user->lwork,user->uwork);
236: VecScale(user->uwork, user->alpha);
237: Gather(G,user->ywork,user->state_scatter,user->uwork,user->design_scatter);
238: return(0);
239: }
241: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *f, Vec G, void *ptr)
242: {
244: PetscReal d1,d2;
245: AppCtx *user = (AppCtx*)ptr;
248: Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);
249: MatMult(user->MQ,user->y,user->dwork);
250: VecAXPY(user->dwork,-1.0,user->d);
251: VecDot(user->dwork,user->dwork,&d1);
252: MatMultTranspose(user->MQ,user->dwork,user->ywork);
254: VecWAXPY(user->uwork,-1.0,user->ur,user->u);
255: MatMult(user->L,user->uwork,user->lwork);
256: VecDot(user->lwork,user->lwork,&d2);
257: MatMultTranspose(user->L,user->lwork,user->uwork);
258: VecScale(user->uwork, user->alpha);
259: *f = 0.5 * (d1 + user->alpha*d2);
260: Gather(G,user->ywork,user->state_scatter,user->uwork,user->design_scatter);
261: return(0);
262: }
264: /* ------------------------------------------------------------------- */
265: /* A
266: MatShell object
267: */
268: PetscErrorCode FormJacobianState(Tao tao, Vec X, Mat J, Mat JPre, Mat JInv, void *ptr)
269: {
271: AppCtx *user = (AppCtx*)ptr;
274: Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);
275: /* DSG = Div * (1/Av_u) * Grad */
276: VecSet(user->uwork,0);
277: VecAXPY(user->uwork,-1.0,user->u);
278: VecExp(user->uwork);
279: MatMult(user->Av,user->uwork,user->Av_u);
280: VecCopy(user->Av_u,user->Swork);
281: VecReciprocal(user->Swork);
282: if (user->use_ptap) {
283: MatDiagonalSet(user->Diag,user->Swork,INSERT_VALUES);
284: MatPtAP(user->Diag,user->Grad,MAT_REUSE_MATRIX,1.0,&user->DSG);
285: } else {
286: MatCopy(user->Div,user->Divwork,SAME_NONZERO_PATTERN);
287: MatDiagonalScale(user->Divwork,NULL,user->Swork);
288: MatProductNumeric(user->DSG);
289: }
290: return(0);
291: }
292: /* ------------------------------------------------------------------- */
293: /* B */
294: PetscErrorCode FormJacobianDesign(Tao tao, Vec X, Mat J, void *ptr)
295: {
297: AppCtx *user = (AppCtx*)ptr;
300: Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);
301: return(0);
302: }
304: PetscErrorCode StateBlockMatMult(Mat J_shell, Vec X, Vec Y)
305: {
307: PetscReal sum;
308: AppCtx *user;
311: MatShellGetContext(J_shell,&user);
312: MatMult(user->DSG,X,Y);
313: VecSum(X,&sum);
314: sum /= user->ndesign;
315: VecShift(Y,sum);
316: return(0);
317: }
319: PetscErrorCode StateMatMult(Mat J_shell, Vec X, Vec Y)
320: {
322: PetscInt i;
323: AppCtx *user;
326: MatShellGetContext(J_shell,&user);
327: if (user->ns == 1) {
328: MatMult(user->JsBlock,X,Y);
329: } else {
330: for (i=0;i<user->ns;i++) {
331: Scatter(X,user->subq,user->yi_scatter[i],0,0);
332: Scatter(Y,user->suby,user->yi_scatter[i],0,0);
333: MatMult(user->JsBlock,user->subq,user->suby);
334: Gather(Y,user->suby,user->yi_scatter[i],0,0);
335: }
336: }
337: return(0);
338: }
340: PetscErrorCode StateInvMatMult(Mat J_shell, Vec X, Vec Y)
341: {
343: PetscInt its,i;
344: AppCtx *user;
347: MatShellGetContext(J_shell,&user);
348: KSPSetOperators(user->solver,user->JsBlock,user->DSG);
349: if (Y == user->ytrue) {
350: /* First solve is done using true solution to set up problem */
351: KSPSetTolerances(user->solver,1e-8,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
352: } else {
353: KSPSetTolerances(user->solver,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
354: }
355: if (user->ns == 1) {
356: KSPSolve(user->solver,X,Y);
357: KSPGetIterationNumber(user->solver,&its);
358: user->ksp_its+=its;
359: } else {
360: for (i=0;i<user->ns;i++) {
361: Scatter(X,user->subq,user->yi_scatter[i],0,0);
362: Scatter(Y,user->suby,user->yi_scatter[i],0,0);
363: KSPSolve(user->solver,user->subq,user->suby);
364: KSPGetIterationNumber(user->solver,&its);
365: user->ksp_its+=its;
366: Gather(Y,user->suby,user->yi_scatter[i],0,0);
367: }
368: }
369: return(0);
370: }
371: PetscErrorCode QMatMult(Mat J_shell, Vec X, Vec Y)
372: {
374: AppCtx *user;
375: PetscInt i;
378: MatShellGetContext(J_shell,&user);
379: if (user->ns == 1) {
380: MatMult(user->Q,X,Y);
381: } else {
382: for (i=0;i<user->ns;i++) {
383: Scatter(X,user->subq,user->yi_scatter[i],0,0);
384: Scatter(Y,user->subd,user->di_scatter[i],0,0);
385: MatMult(user->Q,user->subq,user->subd);
386: Gather(Y,user->subd,user->di_scatter[i],0,0);
387: }
388: }
389: return(0);
390: }
392: PetscErrorCode QMatMultTranspose(Mat J_shell, Vec X, Vec Y)
393: {
395: AppCtx *user;
396: PetscInt i;
399: MatShellGetContext(J_shell,&user);
400: if (user->ns == 1) {
401: MatMultTranspose(user->Q,X,Y);
402: } else {
403: for (i=0;i<user->ns;i++) {
404: Scatter(X,user->subd,user->di_scatter[i],0,0);
405: Scatter(Y,user->suby,user->yi_scatter[i],0,0);
406: MatMultTranspose(user->Q,user->subd,user->suby);
407: Gather(Y,user->suby,user->yi_scatter[i],0,0);
408: }
409: }
410: return(0);
411: }
413: PetscErrorCode DesignMatMult(Mat J_shell, Vec X, Vec Y)
414: {
416: PetscInt i;
417: AppCtx *user;
420: MatShellGetContext(J_shell,&user);
422: /* sdiag(1./v) */
423: VecSet(user->uwork,0);
424: VecAXPY(user->uwork,-1.0,user->u);
425: VecExp(user->uwork);
427: /* sdiag(1./((Av*(1./v)).^2)) */
428: MatMult(user->Av,user->uwork,user->Swork);
429: VecPointwiseMult(user->Swork,user->Swork,user->Swork);
430: VecReciprocal(user->Swork);
432: /* (Av * (sdiag(1./v) * b)) */
433: VecPointwiseMult(user->uwork,user->uwork,X);
434: MatMult(user->Av,user->uwork,user->Twork);
436: /* (sdiag(1./((Av*(1./v)).^2)) * (Av * (sdiag(1./v) * b))) */
437: VecPointwiseMult(user->Swork,user->Twork,user->Swork);
439: if (user->ns == 1) {
440: /* (sdiag(Grad*y(:,i)) */
441: MatMult(user->Grad,user->y,user->Twork);
443: /* Div * (sdiag(Grad*y(:,i)) * (sdiag(1./((Av*(1./v)).^2)) * (Av * (sdiag(1./v) * b)))) */
444: VecPointwiseMult(user->Swork,user->Twork,user->Swork);
445: MatMultTranspose(user->Grad,user->Swork,Y);
446: } else {
447: for (i=0;i<user->ns;i++) {
448: Scatter(user->y,user->suby,user->yi_scatter[i],0,0);
449: Scatter(Y,user->subq,user->yi_scatter[i],0,0);
451: MatMult(user->Grad,user->suby,user->Twork);
452: VecPointwiseMult(user->Twork,user->Twork,user->Swork);
453: MatMultTranspose(user->Grad,user->Twork,user->subq);
454: Gather(user->y,user->suby,user->yi_scatter[i],0,0);
455: Gather(Y,user->subq,user->yi_scatter[i],0,0);
456: }
457: }
458: return(0);
459: }
461: PetscErrorCode DesignMatMultTranspose(Mat J_shell, Vec X, Vec Y)
462: {
464: PetscInt i;
465: AppCtx *user;
468: MatShellGetContext(J_shell,&user);
469: VecZeroEntries(Y);
471: /* Sdiag = 1./((Av*(1./v)).^2) */
472: VecSet(user->uwork,0);
473: VecAXPY(user->uwork,-1.0,user->u);
474: VecExp(user->uwork);
475: MatMult(user->Av,user->uwork,user->Swork);
476: VecPointwiseMult(user->Sdiag,user->Swork,user->Swork);
477: VecReciprocal(user->Sdiag);
479: for (i=0;i<user->ns;i++) {
480: Scatter(X,user->subq,user->yi_scatter[i],0,0);
481: Scatter(user->y,user->suby,user->yi_scatter[i],0,0);
483: /* Swork = (Div' * b(:,i)) */
484: MatMult(user->Grad,user->subq,user->Swork);
486: /* Twork = Grad*y(:,i) */
487: MatMult(user->Grad,user->suby,user->Twork);
489: /* Twork = sdiag(Twork) * Swork */
490: VecPointwiseMult(user->Twork,user->Swork,user->Twork);
492: /* Swork = pointwisemult(Sdiag,Twork) */
493: VecPointwiseMult(user->Swork,user->Twork,user->Sdiag);
495: /* Ywork = Av' * Swork */
496: MatMultTranspose(user->Av,user->Swork,user->Ywork);
498: /* Ywork = pointwisemult(uwork,Ywork) */
499: VecPointwiseMult(user->Ywork,user->uwork,user->Ywork);
500: VecAXPY(Y,1.0,user->Ywork);
501: Gather(user->y,user->suby,user->yi_scatter[i],0,0);
502: }
503: return(0);
504: }
506: PetscErrorCode FormConstraints(Tao tao, Vec X, Vec C, void *ptr)
507: {
508: /* C=Ay - q A = Div * Sigma * Grad + hx*hx*hx*ones(n,n) */
510: PetscReal sum;
511: PetscInt i;
512: AppCtx *user = (AppCtx*)ptr;
515: Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);
516: if (user->ns == 1) {
517: MatMult(user->Grad,user->y,user->Swork);
518: VecPointwiseDivide(user->Swork,user->Swork,user->Av_u);
519: MatMultTranspose(user->Grad,user->Swork,C);
520: VecSum(user->y,&sum);
521: sum /= user->ndesign;
522: VecShift(C,sum);
523: } else {
524: for (i=0;i<user->ns;i++) {
525: Scatter(user->y,user->suby,user->yi_scatter[i],0,0);
526: Scatter(C,user->subq,user->yi_scatter[i],0,0);
527: MatMult(user->Grad,user->suby,user->Swork);
528: VecPointwiseDivide(user->Swork,user->Swork,user->Av_u);
529: MatMultTranspose(user->Grad,user->Swork,user->subq);
531: VecSum(user->suby,&sum);
532: sum /= user->ndesign;
533: VecShift(user->subq,sum);
535: Gather(user->y,user->suby,user->yi_scatter[i],0,0);
536: Gather(C,user->subq,user->yi_scatter[i],0,0);
537: }
538: }
539: VecAXPY(C,-1.0,user->q);
540: return(0);
541: }
543: PetscErrorCode Scatter(Vec x, Vec sub1, VecScatter scat1, Vec sub2, VecScatter scat2)
544: {
548: VecScatterBegin(scat1,x,sub1,INSERT_VALUES,SCATTER_FORWARD);
549: VecScatterEnd(scat1,x,sub1,INSERT_VALUES,SCATTER_FORWARD);
550: if (sub2) {
551: VecScatterBegin(scat2,x,sub2,INSERT_VALUES,SCATTER_FORWARD);
552: VecScatterEnd(scat2,x,sub2,INSERT_VALUES,SCATTER_FORWARD);
553: }
554: return(0);
555: }
557: PetscErrorCode Gather(Vec x, Vec sub1, VecScatter scat1, Vec sub2, VecScatter scat2)
558: {
562: VecScatterBegin(scat1,sub1,x,INSERT_VALUES,SCATTER_REVERSE);
563: VecScatterEnd(scat1,sub1,x,INSERT_VALUES,SCATTER_REVERSE);
564: if (sub2) {
565: VecScatterBegin(scat2,sub2,x,INSERT_VALUES,SCATTER_REVERSE);
566: VecScatterEnd(scat2,sub2,x,INSERT_VALUES,SCATTER_REVERSE);
567: }
568: return(0);
569: }
571: PetscErrorCode EllipticInitialize(AppCtx *user)
572: {
574: PetscInt m,n,i,j,k,l,linear_index,is,js,ks,ls,istart,iend,iblock;
575: Vec XX,YY,ZZ,XXwork,YYwork,ZZwork,UTwork;
576: PetscReal *x,*y,*z;
577: PetscReal h,meanut;
578: PetscScalar hinv,neg_hinv,half = 0.5,sqrt_beta;
579: PetscInt im,indx1,indx2,indy1,indy2,indz1,indz2,nx,ny,nz;
580: IS is_alldesign,is_allstate;
581: IS is_from_d;
582: IS is_from_y;
583: PetscInt lo,hi,hi2,lo2,ysubnlocal,dsubnlocal;
584: const PetscInt *ranges, *subranges;
585: PetscMPIInt size;
586: PetscReal xri,yri,zri,xim,yim,zim,dx1,dx2,dy1,dy2,dz1,dz2,Dx,Dy,Dz;
587: PetscScalar v,vx,vy,vz;
588: PetscInt offset,subindex,subvec,nrank,kk;
590: PetscScalar xr[64] = {0.4970, 0.8498, 0.7814, 0.6268, 0.7782, 0.6402, 0.3617, 0.3160,
591: 0.3610, 0.5298, 0.6987, 0.3331, 0.7962, 0.5596, 0.3866, 0.6774,
592: 0.5407, 0.4518, 0.6702, 0.6061, 0.7580, 0.8997, 0.5198, 0.8326,
593: 0.2138, 0.9198, 0.3000, 0.2833, 0.8288, 0.7076, 0.1820, 0.0728,
594: 0.8447, 0.2367, 0.3239, 0.6413, 0.3114, 0.4731, 0.1192, 0.9273,
595: 0.5724, 0.4331, 0.5136, 0.3547, 0.4413, 0.2602, 0.5698, 0.7278,
596: 0.5261, 0.6230, 0.2454, 0.3948, 0.7479, 0.6582, 0.4660, 0.5594,
597: 0.7574, 0.1143, 0.5900, 0.1065, 0.4260, 0.3294, 0.8276, 0.0756};
599: PetscScalar yr[64] = {0.7345, 0.9120, 0.9288, 0.7528, 0.4463, 0.4985, 0.2497, 0.6256,
600: 0.3425, 0.9026, 0.6983, 0.4230, 0.7140, 0.2970, 0.4474, 0.8792,
601: 0.6604, 0.2485, 0.7968, 0.6127, 0.1796, 0.2437, 0.5938, 0.6137,
602: 0.3867, 0.5658, 0.4575, 0.1009, 0.0863, 0.3361, 0.0738, 0.3985,
603: 0.6602, 0.1437, 0.0934, 0.5983, 0.5950, 0.0763, 0.0768, 0.2288,
604: 0.5761, 0.1129, 0.3841, 0.6150, 0.6904, 0.6686, 0.1361, 0.4601,
605: 0.4491, 0.3716, 0.1969, 0.6537, 0.6743, 0.6991, 0.4811, 0.5480,
606: 0.1684, 0.4569, 0.6889, 0.8437, 0.3015, 0.2854, 0.8199, 0.2658};
608: PetscScalar zr[64] = {0.7668, 0.8573, 0.2654, 0.2719, 0.1060, 0.1311, 0.6232, 0.2295,
609: 0.8009, 0.2147, 0.2119, 0.9325, 0.4473, 0.3600, 0.3374, 0.3819,
610: 0.4066, 0.5801, 0.1673, 0.0959, 0.4638, 0.8236, 0.8800, 0.2939,
611: 0.2028, 0.8262, 0.2706, 0.6276, 0.9085, 0.6443, 0.8241, 0.0712,
612: 0.1824, 0.7789, 0.4389, 0.8415, 0.7055, 0.6639, 0.3653, 0.2078,
613: 0.1987, 0.2297, 0.4321, 0.8115, 0.4915, 0.7764, 0.4657, 0.4627,
614: 0.4569, 0.4232, 0.8514, 0.0674, 0.3227, 0.1055, 0.6690, 0.6313,
615: 0.9226, 0.5461, 0.4126, 0.2364, 0.6096, 0.7042, 0.3914, 0.0711};
618: MPI_Comm_size(PETSC_COMM_WORLD,&size);
619: PetscLogStageRegister("Elliptic Setup",&user->stages[0]);
620: PetscLogStagePush(user->stages[0]);
622: /* Create u,y,c,x */
623: VecCreate(PETSC_COMM_WORLD,&user->u);
624: VecCreate(PETSC_COMM_WORLD,&user->y);
625: VecCreate(PETSC_COMM_WORLD,&user->c);
626: VecSetSizes(user->u,PETSC_DECIDE,user->ndesign);
627: VecSetFromOptions(user->u);
628: VecGetLocalSize(user->u,&ysubnlocal);
629: VecSetSizes(user->y,ysubnlocal*user->ns,user->nstate);
630: VecSetSizes(user->c,ysubnlocal*user->ns,user->m);
631: VecSetFromOptions(user->y);
632: VecSetFromOptions(user->c);
634: /*
635: *******************************
636: Create scatters for x <-> y,u
637: *******************************
639: If the state vector y and design vector u are partitioned as
640: [y_1; y_2; ...; y_np] and [u_1; u_2; ...; u_np] (with np = # of processors),
641: then the solution vector x is organized as
642: [y_1; u_1; y_2; u_2; ...; y_np; u_np].
643: The index sets user->s_is and user->d_is correspond to the indices of the
644: state and design variables owned by the current processor.
645: */
646: VecCreate(PETSC_COMM_WORLD,&user->x);
648: VecGetOwnershipRange(user->y,&lo,&hi);
649: VecGetOwnershipRange(user->u,&lo2,&hi2);
651: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo,1,&is_allstate);
652: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo+lo2,1,&user->s_is);
653: ISCreateStride(PETSC_COMM_SELF,hi2-lo2,lo2,1,&is_alldesign);
654: ISCreateStride(PETSC_COMM_SELF,hi2-lo2,hi+lo2,1,&user->d_is);
656: VecSetSizes(user->x,hi-lo+hi2-lo2,user->n);
657: VecSetFromOptions(user->x);
659: VecScatterCreate(user->x,user->s_is,user->y,is_allstate,&user->state_scatter);
660: VecScatterCreate(user->x,user->d_is,user->u,is_alldesign,&user->design_scatter);
661: ISDestroy(&is_alldesign);
662: ISDestroy(&is_allstate);
663: /*
664: *******************************
665: Create scatter from y to y_1,y_2,...,y_ns
666: *******************************
667: */
668: PetscMalloc1(user->ns,&user->yi_scatter);
669: VecDuplicate(user->u,&user->suby);
670: VecDuplicate(user->u,&user->subq);
672: VecGetOwnershipRange(user->y,&lo2,&hi2);
673: istart = 0;
674: for (i=0; i<user->ns; i++) {
675: VecGetOwnershipRange(user->suby,&lo,&hi);
676: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo2+istart,1,&is_from_y);
677: VecScatterCreate(user->y,is_from_y,user->suby,NULL,&user->yi_scatter[i]);
678: istart = istart + hi-lo;
679: ISDestroy(&is_from_y);
680: }
681: /*
682: *******************************
683: Create scatter from d to d_1,d_2,...,d_ns
684: *******************************
685: */
686: VecCreate(PETSC_COMM_WORLD,&user->subd);
687: VecSetSizes(user->subd,PETSC_DECIDE,user->ndata);
688: VecSetFromOptions(user->subd);
689: VecCreate(PETSC_COMM_WORLD,&user->d);
690: VecGetLocalSize(user->subd,&dsubnlocal);
691: VecSetSizes(user->d,dsubnlocal*user->ns,user->ndata*user->ns);
692: VecSetFromOptions(user->d);
693: PetscMalloc1(user->ns,&user->di_scatter);
695: VecGetOwnershipRange(user->d,&lo2,&hi2);
696: istart = 0;
697: for (i=0; i<user->ns; i++) {
698: VecGetOwnershipRange(user->subd,&lo,&hi);
699: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo2+istart,1,&is_from_d);
700: VecScatterCreate(user->d,is_from_d,user->subd,NULL,&user->di_scatter[i]);
701: istart = istart + hi-lo;
702: ISDestroy(&is_from_d);
703: }
705: PetscMalloc1(user->mx,&x);
706: PetscMalloc1(user->mx,&y);
707: PetscMalloc1(user->mx,&z);
709: user->ksp_its = 0;
710: user->ksp_its_initial = 0;
712: n = user->mx * user->mx * user->mx;
713: m = 3 * user->mx * user->mx * (user->mx-1);
714: sqrt_beta = PetscSqrtScalar(user->beta);
716: VecCreate(PETSC_COMM_WORLD,&XX);
717: VecCreate(PETSC_COMM_WORLD,&user->q);
718: VecSetSizes(XX,ysubnlocal,n);
719: VecSetSizes(user->q,ysubnlocal*user->ns,user->m);
720: VecSetFromOptions(XX);
721: VecSetFromOptions(user->q);
723: VecDuplicate(XX,&YY);
724: VecDuplicate(XX,&ZZ);
725: VecDuplicate(XX,&XXwork);
726: VecDuplicate(XX,&YYwork);
727: VecDuplicate(XX,&ZZwork);
728: VecDuplicate(XX,&UTwork);
729: VecDuplicate(XX,&user->utrue);
731: /* map for striding q */
732: VecGetOwnershipRanges(user->q,&ranges);
733: VecGetOwnershipRanges(user->u,&subranges);
735: VecGetOwnershipRange(user->q,&lo2,&hi2);
736: VecGetOwnershipRange(user->u,&lo,&hi);
737: /* Generate 3D grid, and collect ns (1<=ns<=8) right-hand-side vectors into user->q */
738: h = 1.0/user->mx;
739: hinv = user->mx;
740: neg_hinv = -hinv;
742: VecGetOwnershipRange(XX,&istart,&iend);
743: for (linear_index=istart; linear_index<iend; linear_index++) {
744: i = linear_index % user->mx;
745: j = ((linear_index-i)/user->mx) % user->mx;
746: k = ((linear_index-i)/user->mx-j) / user->mx;
747: vx = h*(i+0.5);
748: vy = h*(j+0.5);
749: vz = h*(k+0.5);
750: VecSetValues(XX,1,&linear_index,&vx,INSERT_VALUES);
751: VecSetValues(YY,1,&linear_index,&vy,INSERT_VALUES);
752: VecSetValues(ZZ,1,&linear_index,&vz,INSERT_VALUES);
753: for (is=0; is<2; is++) {
754: for (js=0; js<2; js++) {
755: for (ks=0; ks<2; ks++) {
756: ls = is*4 + js*2 + ks;
757: if (ls<user->ns) {
758: l =ls*n + linear_index;
759: /* remap */
760: subindex = l%n;
761: subvec = l/n;
762: nrank=0;
763: while (subindex >= subranges[nrank+1]) nrank++;
764: offset = subindex - subranges[nrank];
765: istart=0;
766: for (kk=0;kk<nrank;kk++) istart+=user->ns*(subranges[kk+1]-subranges[kk]);
767: istart += (subranges[nrank+1]-subranges[nrank])*subvec;
768: l = istart+offset;
769: 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));
770: VecSetValues(user->q,1,&l,&v,INSERT_VALUES);
771: }
772: }
773: }
774: }
775: }
777: VecAssemblyBegin(XX);
778: VecAssemblyEnd(XX);
779: VecAssemblyBegin(YY);
780: VecAssemblyEnd(YY);
781: VecAssemblyBegin(ZZ);
782: VecAssemblyEnd(ZZ);
783: VecAssemblyBegin(user->q);
784: VecAssemblyEnd(user->q);
786: /* Compute true parameter function
787: 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) */
788: VecCopy(XX,XXwork);
789: VecCopy(YY,YYwork);
790: VecCopy(ZZ,ZZwork);
792: VecShift(XXwork,-0.25);
793: VecShift(YYwork,-0.25);
794: VecShift(ZZwork,-0.25);
796: VecPointwiseMult(XXwork,XXwork,XXwork);
797: VecPointwiseMult(YYwork,YYwork,YYwork);
798: VecPointwiseMult(ZZwork,ZZwork,ZZwork);
800: VecCopy(XXwork,UTwork);
801: VecAXPY(UTwork,1.0,YYwork);
802: VecAXPY(UTwork,1.0,ZZwork);
803: VecScale(UTwork,-20.0);
804: VecExp(UTwork);
805: VecCopy(UTwork,user->utrue);
807: VecCopy(XX,XXwork);
808: VecCopy(YY,YYwork);
809: VecCopy(ZZ,ZZwork);
811: VecShift(XXwork,-0.75);
812: VecShift(YYwork,-0.75);
813: VecShift(ZZwork,-0.75);
815: VecPointwiseMult(XXwork,XXwork,XXwork);
816: VecPointwiseMult(YYwork,YYwork,YYwork);
817: VecPointwiseMult(ZZwork,ZZwork,ZZwork);
819: VecCopy(XXwork,UTwork);
820: VecAXPY(UTwork,1.0,YYwork);
821: VecAXPY(UTwork,1.0,ZZwork);
822: VecScale(UTwork,-20.0);
823: VecExp(UTwork);
825: VecAXPY(user->utrue,-1.0,UTwork);
827: VecDestroy(&XX);
828: VecDestroy(&YY);
829: VecDestroy(&ZZ);
830: VecDestroy(&XXwork);
831: VecDestroy(&YYwork);
832: VecDestroy(&ZZwork);
833: VecDestroy(&UTwork);
835: /* Initial guess and reference model */
836: VecDuplicate(user->utrue,&user->ur);
837: VecSum(user->utrue,&meanut);
838: meanut = meanut / n;
839: VecSet(user->ur,meanut);
840: VecCopy(user->ur,user->u);
842: /* Generate Grad matrix */
843: MatCreate(PETSC_COMM_WORLD,&user->Grad);
844: MatSetSizes(user->Grad,PETSC_DECIDE,ysubnlocal,m,n);
845: MatSetFromOptions(user->Grad);
846: MatMPIAIJSetPreallocation(user->Grad,2,NULL,2,NULL);
847: MatSeqAIJSetPreallocation(user->Grad,2,NULL);
848: MatGetOwnershipRange(user->Grad,&istart,&iend);
850: for (i=istart; i<iend; i++) {
851: if (i<m/3) {
852: iblock = i / (user->mx-1);
853: j = iblock*user->mx + (i % (user->mx-1));
854: MatSetValues(user->Grad,1,&i,1,&j,&neg_hinv,INSERT_VALUES);
855: j = j+1;
856: MatSetValues(user->Grad,1,&i,1,&j,&hinv,INSERT_VALUES);
857: }
858: if (i>=m/3 && i<2*m/3) {
859: iblock = (i-m/3) / (user->mx*(user->mx-1));
860: j = iblock*user->mx*user->mx + ((i-m/3) % (user->mx*(user->mx-1)));
861: MatSetValues(user->Grad,1,&i,1,&j,&neg_hinv,INSERT_VALUES);
862: j = j + user->mx;
863: MatSetValues(user->Grad,1,&i,1,&j,&hinv,INSERT_VALUES);
864: }
865: if (i>=2*m/3) {
866: j = i-2*m/3;
867: MatSetValues(user->Grad,1,&i,1,&j,&neg_hinv,INSERT_VALUES);
868: j = j + user->mx*user->mx;
869: MatSetValues(user->Grad,1,&i,1,&j,&hinv,INSERT_VALUES);
870: }
871: }
873: MatAssemblyBegin(user->Grad,MAT_FINAL_ASSEMBLY);
874: MatAssemblyEnd(user->Grad,MAT_FINAL_ASSEMBLY);
876: /* Generate arithmetic averaging matrix Av */
877: MatCreate(PETSC_COMM_WORLD,&user->Av);
878: MatSetSizes(user->Av,PETSC_DECIDE,ysubnlocal,m,n);
879: MatSetFromOptions(user->Av);
880: MatMPIAIJSetPreallocation(user->Av,2,NULL,2,NULL);
881: MatSeqAIJSetPreallocation(user->Av,2,NULL);
882: MatGetOwnershipRange(user->Av,&istart,&iend);
884: for (i=istart; i<iend; i++) {
885: if (i<m/3) {
886: iblock = i / (user->mx-1);
887: j = iblock*user->mx + (i % (user->mx-1));
888: MatSetValues(user->Av,1,&i,1,&j,&half,INSERT_VALUES);
889: j = j+1;
890: MatSetValues(user->Av,1,&i,1,&j,&half,INSERT_VALUES);
891: }
892: if (i>=m/3 && i<2*m/3) {
893: iblock = (i-m/3) / (user->mx*(user->mx-1));
894: j = iblock*user->mx*user->mx + ((i-m/3) % (user->mx*(user->mx-1)));
895: MatSetValues(user->Av,1,&i,1,&j,&half,INSERT_VALUES);
896: j = j + user->mx;
897: MatSetValues(user->Av,1,&i,1,&j,&half,INSERT_VALUES);
898: }
899: if (i>=2*m/3) {
900: j = i-2*m/3;
901: MatSetValues(user->Av,1,&i,1,&j,&half,INSERT_VALUES);
902: j = j + user->mx*user->mx;
903: MatSetValues(user->Av,1,&i,1,&j,&half,INSERT_VALUES);
904: }
905: }
907: MatAssemblyBegin(user->Av,MAT_FINAL_ASSEMBLY);
908: MatAssemblyEnd(user->Av,MAT_FINAL_ASSEMBLY);
910: MatCreate(PETSC_COMM_WORLD,&user->L);
911: MatSetSizes(user->L,PETSC_DECIDE,ysubnlocal,m+n,n);
912: MatSetFromOptions(user->L);
913: MatMPIAIJSetPreallocation(user->L,2,NULL,2,NULL);
914: MatSeqAIJSetPreallocation(user->L,2,NULL);
915: MatGetOwnershipRange(user->L,&istart,&iend);
917: for (i=istart; i<iend; i++) {
918: if (i<m/3) {
919: iblock = i / (user->mx-1);
920: j = iblock*user->mx + (i % (user->mx-1));
921: MatSetValues(user->L,1,&i,1,&j,&neg_hinv,INSERT_VALUES);
922: j = j+1;
923: MatSetValues(user->L,1,&i,1,&j,&hinv,INSERT_VALUES);
924: }
925: if (i>=m/3 && i<2*m/3) {
926: iblock = (i-m/3) / (user->mx*(user->mx-1));
927: j = iblock*user->mx*user->mx + ((i-m/3) % (user->mx*(user->mx-1)));
928: MatSetValues(user->L,1,&i,1,&j,&neg_hinv,INSERT_VALUES);
929: j = j + user->mx;
930: MatSetValues(user->L,1,&i,1,&j,&hinv,INSERT_VALUES);
931: }
932: if (i>=2*m/3 && i<m) {
933: j = i-2*m/3;
934: MatSetValues(user->L,1,&i,1,&j,&neg_hinv,INSERT_VALUES);
935: j = j + user->mx*user->mx;
936: MatSetValues(user->L,1,&i,1,&j,&hinv,INSERT_VALUES);
937: }
938: if (i>=m) {
939: j = i - m;
940: MatSetValues(user->L,1,&i,1,&j,&sqrt_beta,INSERT_VALUES);
941: }
942: }
943: MatAssemblyBegin(user->L,MAT_FINAL_ASSEMBLY);
944: MatAssemblyEnd(user->L,MAT_FINAL_ASSEMBLY);
945: MatScale(user->L,PetscPowScalar(h,1.5));
947: /* Generate Div matrix */
948: if (!user->use_ptap) {
949: /* Generate Div matrix */
950: MatCreate(PETSC_COMM_WORLD,&user->Div);
951: MatSetSizes(user->Div,ysubnlocal,PETSC_DECIDE,n,m);
952: MatSetFromOptions(user->Div);
953: MatMPIAIJSetPreallocation(user->Div,4,NULL,4,NULL);
954: MatSeqAIJSetPreallocation(user->Div,6,NULL);
955: MatGetOwnershipRange(user->Grad,&istart,&iend);
957: for (i=istart; i<iend; i++) {
958: if (i<m/3) {
959: iblock = i / (user->mx-1);
960: j = iblock*user->mx + (i % (user->mx-1));
961: MatSetValues(user->Div,1,&j,1,&i,&neg_hinv,INSERT_VALUES);
962: j = j+1;
963: MatSetValues(user->Div,1,&j,1,&i,&hinv,INSERT_VALUES);
964: }
965: if (i>=m/3 && i<2*m/3) {
966: iblock = (i-m/3) / (user->mx*(user->mx-1));
967: j = iblock*user->mx*user->mx + ((i-m/3) % (user->mx*(user->mx-1)));
968: MatSetValues(user->Div,1,&j,1,&i,&neg_hinv,INSERT_VALUES);
969: j = j + user->mx;
970: MatSetValues(user->Div,1,&j,1,&i,&hinv,INSERT_VALUES);
971: }
972: if (i>=2*m/3) {
973: j = i-2*m/3;
974: MatSetValues(user->Div,1,&j,1,&i,&neg_hinv,INSERT_VALUES);
975: j = j + user->mx*user->mx;
976: MatSetValues(user->Div,1,&j,1,&i,&hinv,INSERT_VALUES);
977: }
978: }
980: MatAssemblyBegin(user->Div,MAT_FINAL_ASSEMBLY);
981: MatAssemblyEnd(user->Div,MAT_FINAL_ASSEMBLY);
982: MatDuplicate(user->Div,MAT_SHARE_NONZERO_PATTERN,&user->Divwork);
983: } else {
984: MatCreate(PETSC_COMM_WORLD,&user->Diag);
985: MatSetSizes(user->Diag,PETSC_DECIDE,PETSC_DECIDE,m,m);
986: MatSetFromOptions(user->Diag);
987: MatMPIAIJSetPreallocation(user->Diag,1,NULL,0,NULL);
988: MatSeqAIJSetPreallocation(user->Diag,1,NULL);
989: }
991: /* Build work vectors and matrices */
992: VecCreate(PETSC_COMM_WORLD,&user->S);
993: VecSetSizes(user->S, PETSC_DECIDE, m);
994: VecSetFromOptions(user->S);
996: VecCreate(PETSC_COMM_WORLD,&user->lwork);
997: VecSetSizes(user->lwork,PETSC_DECIDE,m+user->mx*user->mx*user->mx);
998: VecSetFromOptions(user->lwork);
1000: MatDuplicate(user->Av,MAT_SHARE_NONZERO_PATTERN,&user->Avwork);
1002: VecDuplicate(user->S,&user->Swork);
1003: VecDuplicate(user->S,&user->Sdiag);
1004: VecDuplicate(user->S,&user->Av_u);
1005: VecDuplicate(user->S,&user->Twork);
1006: VecDuplicate(user->y,&user->ywork);
1007: VecDuplicate(user->u,&user->Ywork);
1008: VecDuplicate(user->u,&user->uwork);
1009: VecDuplicate(user->u,&user->js_diag);
1010: VecDuplicate(user->c,&user->cwork);
1011: VecDuplicate(user->d,&user->dwork);
1013: /* Create a matrix-free shell user->Jd for computing B*x */
1014: MatCreateShell(PETSC_COMM_WORLD,ysubnlocal*user->ns,ysubnlocal,user->nstate,user->ndesign,user,&user->Jd);
1015: MatShellSetOperation(user->Jd,MATOP_MULT,(void(*)(void))DesignMatMult);
1016: MatShellSetOperation(user->Jd,MATOP_MULT_TRANSPOSE,(void(*)(void))DesignMatMultTranspose);
1018: /* Compute true state function ytrue given utrue */
1019: VecDuplicate(user->y,&user->ytrue);
1021: /* First compute Av_u = Av*exp(-u) */
1022: VecSet(user->uwork, 0);
1023: VecAXPY(user->uwork,-1.0,user->utrue); /* Note: user->utrue */
1024: VecExp(user->uwork);
1025: MatMult(user->Av,user->uwork,user->Av_u);
1027: /* Next form DSG = Div*S*Grad */
1028: VecCopy(user->Av_u,user->Swork);
1029: VecReciprocal(user->Swork);
1030: if (user->use_ptap) {
1031: MatDiagonalSet(user->Diag,user->Swork,INSERT_VALUES);
1032: MatPtAP(user->Diag,user->Grad,MAT_INITIAL_MATRIX,1.0,&user->DSG);
1033: } else {
1034: MatCopy(user->Div,user->Divwork,SAME_NONZERO_PATTERN);
1035: MatDiagonalScale(user->Divwork,NULL,user->Swork);
1037: MatMatMult(user->Divwork,user->Grad,MAT_INITIAL_MATRIX,1.0,&user->DSG);
1038: }
1040: MatSetOption(user->DSG,MAT_SYMMETRIC,PETSC_TRUE);
1041: MatSetOption(user->DSG,MAT_SYMMETRY_ETERNAL,PETSC_TRUE);
1043: if (user->use_lrc == PETSC_TRUE) {
1044: v=PetscSqrtReal(1.0 /user->ndesign);
1045: PetscMalloc1(user->ndesign,&user->ones);
1047: for (i=0;i<user->ndesign;i++) {
1048: user->ones[i]=v;
1049: }
1050: MatCreateDense(PETSC_COMM_WORLD,ysubnlocal,PETSC_DECIDE,user->ndesign,1,user->ones,&user->Ones);
1051: MatAssemblyBegin(user->Ones, MAT_FINAL_ASSEMBLY);
1052: MatAssemblyEnd(user->Ones, MAT_FINAL_ASSEMBLY);
1053: MatCreateLRC(user->DSG,user->Ones,NULL,user->Ones,&user->JsBlock);
1054: MatSetUp(user->JsBlock);
1055: } else {
1056: /* Create matrix-free shell user->Js for computing (A + h^3*e*e^T)*x */
1057: MatCreateShell(PETSC_COMM_WORLD,ysubnlocal,ysubnlocal,user->ndesign,user->ndesign,user,&user->JsBlock);
1058: MatShellSetOperation(user->JsBlock,MATOP_MULT,(void(*)(void))StateBlockMatMult);
1059: MatShellSetOperation(user->JsBlock,MATOP_MULT_TRANSPOSE,(void(*)(void))StateBlockMatMult);
1060: }
1061: MatSetOption(user->JsBlock,MAT_SYMMETRIC,PETSC_TRUE);
1062: MatSetOption(user->JsBlock,MAT_SYMMETRY_ETERNAL,PETSC_TRUE);
1063: MatCreateShell(PETSC_COMM_WORLD,ysubnlocal*user->ns,ysubnlocal*user->ns,user->nstate,user->nstate,user,&user->Js);
1064: MatShellSetOperation(user->Js,MATOP_MULT,(void(*)(void))StateMatMult);
1065: MatShellSetOperation(user->Js,MATOP_MULT_TRANSPOSE,(void(*)(void))StateMatMult);
1066: MatSetOption(user->Js,MAT_SYMMETRIC,PETSC_TRUE);
1067: MatSetOption(user->Js,MAT_SYMMETRY_ETERNAL,PETSC_TRUE);
1069: MatCreateShell(PETSC_COMM_WORLD,ysubnlocal*user->ns,ysubnlocal*user->ns,user->nstate,user->nstate,user,&user->JsInv);
1070: MatShellSetOperation(user->JsInv,MATOP_MULT,(void(*)(void))StateInvMatMult);
1071: MatShellSetOperation(user->JsInv,MATOP_MULT_TRANSPOSE,(void(*)(void))StateInvMatMult);
1072: MatSetOption(user->JsInv,MAT_SYMMETRIC,PETSC_TRUE);
1073: MatSetOption(user->JsInv,MAT_SYMMETRY_ETERNAL,PETSC_TRUE);
1075: MatSetOption(user->DSG,MAT_SYMMETRIC,PETSC_TRUE);
1076: MatSetOption(user->DSG,MAT_SYMMETRY_ETERNAL,PETSC_TRUE);
1077: /* Now solve for ytrue */
1078: KSPCreate(PETSC_COMM_WORLD,&user->solver);
1079: KSPSetFromOptions(user->solver);
1081: KSPSetOperators(user->solver,user->JsBlock,user->DSG);
1083: MatMult(user->JsInv,user->q,user->ytrue);
1084: /* First compute Av_u = Av*exp(-u) */
1085: VecSet(user->uwork,0);
1086: VecAXPY(user->uwork,-1.0,user->u); /* Note: user->u */
1087: VecExp(user->uwork);
1088: MatMult(user->Av,user->uwork,user->Av_u);
1090: /* Next update DSG = Div*S*Grad with user->u */
1091: VecCopy(user->Av_u,user->Swork);
1092: VecReciprocal(user->Swork);
1093: if (user->use_ptap) {
1094: MatDiagonalSet(user->Diag,user->Swork,INSERT_VALUES);
1095: MatPtAP(user->Diag,user->Grad,MAT_REUSE_MATRIX,1.0,&user->DSG);
1096: } else {
1097: MatCopy(user->Div,user->Divwork,SAME_NONZERO_PATTERN);
1098: MatDiagonalScale(user->Divwork,NULL,user->Av_u);
1099: MatProductNumeric(user->DSG);
1100: }
1102: /* Now solve for y */
1104: MatMult(user->JsInv,user->q,user->y);
1106: user->ksp_its_initial = user->ksp_its;
1107: user->ksp_its = 0;
1108: /* Construct projection matrix Q (blocks) */
1109: MatCreate(PETSC_COMM_WORLD,&user->Q);
1110: MatSetSizes(user->Q,dsubnlocal,ysubnlocal,user->ndata,user->ndesign);
1111: MatSetFromOptions(user->Q);
1112: MatMPIAIJSetPreallocation(user->Q,8,NULL,8,NULL);
1113: MatSeqAIJSetPreallocation(user->Q,8,NULL);
1115: for (i=0; i<user->mx; i++) {
1116: x[i] = h*(i+0.5);
1117: y[i] = h*(i+0.5);
1118: z[i] = h*(i+0.5);
1119: }
1120: MatGetOwnershipRange(user->Q,&istart,&iend);
1122: nx = user->mx; ny = user->mx; nz = user->mx;
1123: for (i=istart; i<iend; i++) {
1125: xri = xr[i];
1126: im = 0;
1127: xim = x[im];
1128: while (xri>xim && im<nx) {
1129: im = im+1;
1130: xim = x[im];
1131: }
1132: indx1 = im-1;
1133: indx2 = im;
1134: dx1 = xri - x[indx1];
1135: dx2 = x[indx2] - xri;
1137: yri = yr[i];
1138: im = 0;
1139: yim = y[im];
1140: while (yri>yim && im<ny) {
1141: im = im+1;
1142: yim = y[im];
1143: }
1144: indy1 = im-1;
1145: indy2 = im;
1146: dy1 = yri - y[indy1];
1147: dy2 = y[indy2] - yri;
1149: zri = zr[i];
1150: im = 0;
1151: zim = z[im];
1152: while (zri>zim && im<nz) {
1153: im = im+1;
1154: zim = z[im];
1155: }
1156: indz1 = im-1;
1157: indz2 = im;
1158: dz1 = zri - z[indz1];
1159: dz2 = z[indz2] - zri;
1161: Dx = x[indx2] - x[indx1];
1162: Dy = y[indy2] - y[indy1];
1163: Dz = z[indz2] - z[indz1];
1165: j = indx1 + indy1*nx + indz1*nx*ny;
1166: v = (1-dx1/Dx)*(1-dy1/Dy)*(1-dz1/Dz);
1167: MatSetValues(user->Q,1,&i,1,&j,&v,INSERT_VALUES);
1169: j = indx1 + indy1*nx + indz2*nx*ny;
1170: v = (1-dx1/Dx)*(1-dy1/Dy)*(1-dz2/Dz);
1171: MatSetValues(user->Q,1,&i,1,&j,&v,INSERT_VALUES);
1173: j = indx1 + indy2*nx + indz1*nx*ny;
1174: v = (1-dx1/Dx)*(1-dy2/Dy)*(1-dz1/Dz);
1175: MatSetValues(user->Q,1,&i,1,&j,&v,INSERT_VALUES);
1177: j = indx1 + indy2*nx + indz2*nx*ny;
1178: v = (1-dx1/Dx)*(1-dy2/Dy)*(1-dz2/Dz);
1179: MatSetValues(user->Q,1,&i,1,&j,&v,INSERT_VALUES);
1181: j = indx2 + indy1*nx + indz1*nx*ny;
1182: v = (1-dx2/Dx)*(1-dy1/Dy)*(1-dz1/Dz);
1183: MatSetValues(user->Q,1,&i,1,&j,&v,INSERT_VALUES);
1185: j = indx2 + indy1*nx + indz2*nx*ny;
1186: v = (1-dx2/Dx)*(1-dy1/Dy)*(1-dz2/Dz);
1187: MatSetValues(user->Q,1,&i,1,&j,&v,INSERT_VALUES);
1189: j = indx2 + indy2*nx + indz1*nx*ny;
1190: v = (1-dx2/Dx)*(1-dy2/Dy)*(1-dz1/Dz);
1191: MatSetValues(user->Q,1,&i,1,&j,&v,INSERT_VALUES);
1193: j = indx2 + indy2*nx + indz2*nx*ny;
1194: v = (1-dx2/Dx)*(1-dy2/Dy)*(1-dz2/Dz);
1195: MatSetValues(user->Q,1,&i,1,&j,&v,INSERT_VALUES);
1196: }
1198: MatAssemblyBegin(user->Q,MAT_FINAL_ASSEMBLY);
1199: MatAssemblyEnd(user->Q,MAT_FINAL_ASSEMBLY);
1200: /* Create MQ (composed of blocks of Q */
1201: MatCreateShell(PETSC_COMM_WORLD,dsubnlocal*user->ns,PETSC_DECIDE,user->ndata*user->ns,user->nstate,user,&user->MQ);
1202: MatShellSetOperation(user->MQ,MATOP_MULT,(void(*)(void))QMatMult);
1203: MatShellSetOperation(user->MQ,MATOP_MULT_TRANSPOSE,(void(*)(void))QMatMultTranspose);
1205: /* Add noise to the measurement data */
1206: VecSet(user->ywork,1.0);
1207: VecAYPX(user->ywork,user->noise,user->ytrue);
1208: MatMult(user->MQ,user->ywork,user->d);
1210: /* Now that initial conditions have been set, let the user pass tolerance options to the KSP solver */
1211: PetscFree(x);
1212: PetscFree(y);
1213: PetscFree(z);
1214: PetscLogStagePop();
1215: return(0);
1216: }
1218: PetscErrorCode EllipticDestroy(AppCtx *user)
1219: {
1221: PetscInt i;
1224: MatDestroy(&user->DSG);
1225: KSPDestroy(&user->solver);
1226: MatDestroy(&user->Q);
1227: MatDestroy(&user->MQ);
1228: if (!user->use_ptap) {
1229: MatDestroy(&user->Div);
1230: MatDestroy(&user->Divwork);
1231: } else {
1232: MatDestroy(&user->Diag);
1233: }
1234: if (user->use_lrc) {
1235: MatDestroy(&user->Ones);
1236: }
1238: MatDestroy(&user->Grad);
1239: MatDestroy(&user->Av);
1240: MatDestroy(&user->Avwork);
1241: MatDestroy(&user->L);
1242: MatDestroy(&user->Js);
1243: MatDestroy(&user->Jd);
1244: MatDestroy(&user->JsBlock);
1245: MatDestroy(&user->JsInv);
1247: VecDestroy(&user->x);
1248: VecDestroy(&user->u);
1249: VecDestroy(&user->uwork);
1250: VecDestroy(&user->utrue);
1251: VecDestroy(&user->y);
1252: VecDestroy(&user->ywork);
1253: VecDestroy(&user->ytrue);
1254: VecDestroy(&user->c);
1255: VecDestroy(&user->cwork);
1256: VecDestroy(&user->ur);
1257: VecDestroy(&user->q);
1258: VecDestroy(&user->d);
1259: VecDestroy(&user->dwork);
1260: VecDestroy(&user->lwork);
1261: VecDestroy(&user->S);
1262: VecDestroy(&user->Swork);
1263: VecDestroy(&user->Sdiag);
1264: VecDestroy(&user->Ywork);
1265: VecDestroy(&user->Twork);
1266: VecDestroy(&user->Av_u);
1267: VecDestroy(&user->js_diag);
1268: ISDestroy(&user->s_is);
1269: ISDestroy(&user->d_is);
1270: VecDestroy(&user->suby);
1271: VecDestroy(&user->subd);
1272: VecDestroy(&user->subq);
1273: VecScatterDestroy(&user->state_scatter);
1274: VecScatterDestroy(&user->design_scatter);
1275: for (i=0;i<user->ns;i++) {
1276: VecScatterDestroy(&user->yi_scatter[i]);
1277: VecScatterDestroy(&user->di_scatter[i]);
1278: }
1279: PetscFree(user->yi_scatter);
1280: PetscFree(user->di_scatter);
1281: if (user->use_lrc) {
1282: PetscFree(user->ones);
1283: MatDestroy(&user->Ones);
1284: }
1285: return(0);
1286: }
1288: PetscErrorCode EllipticMonitor(Tao tao, void *ptr)
1289: {
1291: Vec X;
1292: PetscReal unorm,ynorm;
1293: AppCtx *user = (AppCtx*)ptr;
1296: TaoGetSolutionVector(tao,&X);
1297: Scatter(X,user->ywork,user->state_scatter,user->uwork,user->design_scatter);
1298: VecAXPY(user->ywork,-1.0,user->ytrue);
1299: VecAXPY(user->uwork,-1.0,user->utrue);
1300: VecNorm(user->uwork,NORM_2,&unorm);
1301: VecNorm(user->ywork,NORM_2,&ynorm);
1302: PetscPrintf(MPI_COMM_WORLD, "||u-ut||=%g ||y-yt||=%g\n",(double)unorm,(double)ynorm);
1303: return(0);
1304: }
1306: /*TEST
1308: build:
1309: requires: !complex
1311: test:
1312: args: -tao_cmonitor -ns 1 -tao_type lcl -tao_gatol 1.e-3 -tao_max_it 11
1313: requires: !single
1315: test:
1316: suffix: 2
1317: args: -tao_cmonitor -tao_type lcl -tao_max_it 11 -use_ptap -use_lrc -ns 1 -tao_gatol 1.e-3
1318: requires: !single
1320: TEST*/