Actual source code: elliptic.c
petsc-3.6.1 2015-08-06
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: TaoSetHistory(); TaoGetHistory();
16: Routines: TaoSolve();
17: Routines: TaoGetConvergedReason(); TaoDestroy();
18: Processors: n
19: T*/
21: typedef struct {
22: PetscInt n; /* Number of total variables */
23: PetscInt m; /* Number of constraints */
24: PetscInt nstate;
25: PetscInt ndesign;
26: PetscInt mx; /* grid points in each direction */
27: PetscInt ns; /* Number of data samples (1<=ns<=8)
28: Currently only ns=1 is supported */
29: PetscInt ndata; /* Number of data points per sample */
30: IS s_is;
31: IS d_is;
33: VecScatter state_scatter;
34: VecScatter design_scatter;
35: VecScatter *yi_scatter, *di_scatter;
36: Vec suby,subq,subd;
37: Mat Js,Jd,JsPrec,JsInv,JsBlock;
39: PetscReal alpha; /* Regularization parameter */
40: PetscReal beta; /* Weight attributed to ||u||^2 in regularization functional */
41: PetscReal noise; /* Amount of noise to add to data */
42: PetscReal *ones;
43: Mat Q;
44: Mat MQ;
45: Mat L;
47: Mat Grad;
48: Mat Av,Avwork;
49: Mat Div, Divwork;
50: Mat DSG;
51: Mat Diag,Ones;
54: Vec q;
55: Vec ur; /* reference */
57: Vec d;
58: Vec dwork;
60: Vec x; /* super vec of y,u */
62: Vec y; /* state variables */
63: Vec ywork;
65: Vec ytrue;
67: Vec u; /* design variables */
68: Vec uwork;
70: Vec utrue;
72: Vec js_diag;
74: Vec c; /* constraint vector */
75: Vec cwork;
77: Vec lwork;
78: Vec S;
79: Vec Swork,Twork,Sdiag,Ywork;
80: Vec Av_u;
82: KSP solver;
83: PC prec;
85: PetscReal tola,tolb,tolc,told;
86: PetscInt ksp_its;
87: PetscInt ksp_its_initial;
88: int stages[10];
89: PetscBool use_ptap;
90: PetscBool use_lrc;
91: } AppCtx;
93: PetscErrorCode FormFunction(Tao, Vec, PetscReal*, void*);
94: PetscErrorCode FormGradient(Tao, Vec, Vec, void*);
95: PetscErrorCode FormFunctionGradient(Tao, Vec, PetscReal*, Vec, void*);
96: PetscErrorCode FormJacobianState(Tao, Vec, Mat, Mat, Mat, void*);
97: PetscErrorCode FormJacobianDesign(Tao, Vec, Mat,void*);
98: PetscErrorCode FormConstraints(Tao, Vec, Vec, void*);
99: PetscErrorCode FormHessian(Tao, Vec, Mat, Mat, void*);
100: PetscErrorCode Gather(Vec, Vec, VecScatter, Vec, VecScatter);
101: PetscErrorCode Scatter(Vec, Vec, VecScatter, Vec, VecScatter);
102: PetscErrorCode EllipticInitialize(AppCtx*);
103: PetscErrorCode EllipticDestroy(AppCtx*);
104: PetscErrorCode EllipticMonitor(Tao, void*);
106: PetscErrorCode StateBlockMatMult(Mat,Vec,Vec);
107: PetscErrorCode StateMatMult(Mat,Vec,Vec);
109: PetscErrorCode StateInvMatMult(Mat,Vec,Vec);
110: PetscErrorCode DesignMatMult(Mat,Vec,Vec);
111: PetscErrorCode DesignMatMultTranspose(Mat,Vec,Vec);
113: PetscErrorCode QMatMult(Mat,Vec,Vec);
114: PetscErrorCode QMatMultTranspose(Mat,Vec,Vec);
116: static char help[]="";
120: int main(int argc, char **argv)
121: {
122: PetscErrorCode ierr;
123: Vec x0;
124: Tao tao;
125: TaoConvergedReason reason;
126: AppCtx user;
127: PetscInt ntests = 1;
128: PetscInt i;
130: PetscInitialize(&argc, &argv, (char*)0,help);
131: user.mx = 8;
132: PetscOptionsBegin(PETSC_COMM_WORLD,NULL,NULL,NULL);
133: PetscOptionsInt("-mx","Number of grid points in each direction","",user.mx,&user.mx,NULL);
134: user.ns = 6;
135: PetscOptionsInt("-ns","Number of data samples (1<=ns<=8)","",user.ns,&user.ns,NULL);
136: user.ndata = 64;
137: PetscOptionsInt("-ndata","Numbers of data points per sample","",user.ndata,&user.ndata,NULL);
138: user.alpha = 0.1;
139: PetscOptionsReal("-alpha","Regularization parameter","",user.alpha,&user.alpha,NULL);
140: user.beta = 0.00001;
141: PetscOptionsReal("-beta","Weight attributed to ||u||^2 in regularization functional","",user.beta,&user.beta,NULL);
142: user.noise = 0.01;
143: PetscOptionsReal("-noise","Amount of noise to add to data","",user.noise,&user.noise,NULL);
145: user.use_ptap = PETSC_FALSE;
146: PetscOptionsBool("-use_ptap","Use ptap matrix for DSG","",user.use_ptap,&user.use_ptap,NULL);
147: user.use_lrc = PETSC_FALSE;
148: PetscOptionsBool("-use_lrc","Use lrc matrix for Js","",user.use_lrc,&user.use_lrc,NULL);
149: user.m = user.ns*user.mx*user.mx*user.mx; /* number of constraints */
150: user.nstate = user.m;
151: user.ndesign = user.mx*user.mx*user.mx;
152: user.n = user.nstate + user.ndesign; /* number of variables */
154: /* Create TAO solver and set desired solution method */
155: TaoCreate(PETSC_COMM_WORLD,&tao);
156: TaoSetType(tao,TAOLCL);
158: /* Set up initial vectors and matrices */
159: EllipticInitialize(&user);
161: Gather(user.x,user.y,user.state_scatter,user.u,user.design_scatter);
162: VecDuplicate(user.x,&x0);
163: VecCopy(user.x,x0);
165: /* Set solution vector with an initial guess */
166: TaoSetInitialVector(tao,user.x);
167: TaoSetObjectiveRoutine(tao, FormFunction, (void *)&user);
168: TaoSetGradientRoutine(tao, FormGradient, (void *)&user);
169: TaoSetConstraintsRoutine(tao, user.c, FormConstraints, (void *)&user);
171: TaoSetJacobianStateRoutine(tao, user.Js, NULL, user.JsInv, FormJacobianState, (void *)&user);
172: TaoSetJacobianDesignRoutine(tao, user.Jd, FormJacobianDesign, (void *)&user);
174: TaoSetStateDesignIS(tao,user.s_is,user.d_is);
175: TaoSetFromOptions(tao);
177: /* SOLVE THE APPLICATION */
178: PetscOptionsInt("-ntests","Number of times to repeat TaoSolve","",ntests,&ntests,NULL);
179: PetscOptionsEnd();
180: PetscLogStageRegister("Trials",&user.stages[1]);
181: PetscLogStagePush(user.stages[1]);
182: for (i=0; i<ntests; i++){
183: TaoSolve(tao);
184: PetscPrintf(PETSC_COMM_WORLD,"KSP Iterations = %D\n",user.ksp_its);
185: VecCopy(x0,user.x);
186: }
187: PetscLogStagePop();
188: PetscBarrier((PetscObject)user.x);
189: PetscPrintf(PETSC_COMM_WORLD,"KSP iterations within initialization: ");
190: PetscPrintf(PETSC_COMM_WORLD,"%D\n",user.ksp_its_initial);
192: TaoGetConvergedReason(tao,&reason);
193: if (reason < 0) {
194: PetscPrintf(MPI_COMM_WORLD, "TAO failed to converge.\n");
195: } else {
196: PetscPrintf(MPI_COMM_WORLD, "Optimization terminated with status %D.\n", reason);
197: }
199: TaoDestroy(&tao);
200: VecDestroy(&x0);
201: EllipticDestroy(&user);
202: PetscFinalize();
203: return 0;
204: }
205: /* ------------------------------------------------------------------- */
208: /*
209: dwork = Qy - d
210: lwork = L*(u-ur)
211: f = 1/2 * (dwork.dwork + alpha*lwork.lwork)
212: */
213: PetscErrorCode FormFunction(Tao tao,Vec X,PetscReal *f,void *ptr)
214: {
216: PetscReal d1=0,d2=0;
217: AppCtx *user = (AppCtx*)ptr;
220: Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);
221: MatMult(user->MQ,user->y,user->dwork);
222: VecAXPY(user->dwork,-1.0,user->d);
223: VecDot(user->dwork,user->dwork,&d1);
224: VecWAXPY(user->uwork,-1.0,user->ur,user->u);
225: MatMult(user->L,user->uwork,user->lwork);
226: VecDot(user->lwork,user->lwork,&d2);
227: *f = 0.5 * (d1 + user->alpha*d2);
228: return(0);
229: }
231: /* ------------------------------------------------------------------- */
234: /*
235: state: g_s = Q' *(Qy - d)
236: design: g_d = alpha*L'*L*(u-ur)
237: */
238: PetscErrorCode FormGradient(Tao tao,Vec X,Vec G,void *ptr)
239: {
241: AppCtx *user = (AppCtx*)ptr;
244: Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);
245: MatMult(user->MQ,user->y,user->dwork);
246: VecAXPY(user->dwork,-1.0,user->d);
247: MatMultTranspose(user->MQ,user->dwork,user->ywork);
248: VecWAXPY(user->uwork,-1.0,user->ur,user->u);
249: MatMult(user->L,user->uwork,user->lwork);
250: MatMultTranspose(user->L,user->lwork,user->uwork);
251: VecScale(user->uwork, user->alpha);
252: Gather(G,user->ywork,user->state_scatter,user->uwork,user->design_scatter);
253: return(0);
254: }
258: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *f, Vec G, void *ptr)
259: {
261: PetscReal d1,d2;
262: AppCtx *user = (AppCtx*)ptr;
265: Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);
266: MatMult(user->MQ,user->y,user->dwork);
267: VecAXPY(user->dwork,-1.0,user->d);
268: VecDot(user->dwork,user->dwork,&d1);
269: MatMultTranspose(user->MQ,user->dwork,user->ywork);
271: VecWAXPY(user->uwork,-1.0,user->ur,user->u);
272: MatMult(user->L,user->uwork,user->lwork);
273: VecDot(user->lwork,user->lwork,&d2);
274: MatMultTranspose(user->L,user->lwork,user->uwork);
275: VecScale(user->uwork, user->alpha);
276: *f = 0.5 * (d1 + user->alpha*d2);
277: Gather(G,user->ywork,user->state_scatter,user->uwork,user->design_scatter);
278: return(0);
279: }
281: /* ------------------------------------------------------------------- */
284: /* A
285: MatShell object
286: */
287: PetscErrorCode FormJacobianState(Tao tao, Vec X, Mat J, Mat JPre, Mat JInv, void *ptr)
288: {
290: AppCtx *user = (AppCtx*)ptr;
293: Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);
294: /* DSG = Div * (1/Av_u) * Grad */
295: VecSet(user->uwork,0);
296: VecAXPY(user->uwork,-1.0,user->u);
297: VecExp(user->uwork);
298: MatMult(user->Av,user->uwork,user->Av_u);
299: VecCopy(user->Av_u,user->Swork);
300: VecReciprocal(user->Swork);
301: if (user->use_ptap) {
302: MatDiagonalSet(user->Diag,user->Swork,INSERT_VALUES);
303: MatPtAP(user->Diag,user->Grad,MAT_REUSE_MATRIX,1.0,&user->DSG);
304: } else {
305: MatCopy(user->Div,user->Divwork,SAME_NONZERO_PATTERN);
306: MatDiagonalScale(user->Divwork,NULL,user->Swork);
307: MatMatMultNumeric(user->Divwork,user->Grad,user->DSG);
308: }
309: return(0);
310: }
311: /* ------------------------------------------------------------------- */
314: /* B */
315: PetscErrorCode FormJacobianDesign(Tao tao, Vec X, Mat J, void *ptr)
316: {
318: AppCtx *user = (AppCtx*)ptr;
321: Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);
322: return(0);
323: }
327: PetscErrorCode StateBlockMatMult(Mat J_shell, Vec X, Vec Y)
328: {
330: PetscReal sum;
331: AppCtx *user;
334: MatShellGetContext(J_shell,(void**)&user);
335: MatMult(user->DSG,X,Y);
336: VecSum(X,&sum);
337: sum /= user->ndesign;
338: VecShift(Y,sum);
339: return(0);
340: }
344: PetscErrorCode StateMatMult(Mat J_shell, Vec X, Vec Y)
345: {
347: PetscInt i;
348: AppCtx *user;
351: MatShellGetContext(J_shell,(void**)&user);
352: if (user->ns == 1) {
353: MatMult(user->JsBlock,X,Y);
354: } else {
355: for (i=0;i<user->ns;i++) {
356: Scatter(X,user->subq,user->yi_scatter[i],0,0);
357: Scatter(Y,user->suby,user->yi_scatter[i],0,0);
358: MatMult(user->JsBlock,user->subq,user->suby);
359: Gather(Y,user->suby,user->yi_scatter[i],0,0);
360: }
361: }
362: return(0);
363: }
367: PetscErrorCode StateInvMatMult(Mat J_shell, Vec X, Vec Y)
368: {
370: PetscInt its,i;
371: AppCtx *user;
374: MatShellGetContext(J_shell,(void**)&user);
375: KSPSetOperators(user->solver,user->JsBlock,user->DSG);
376: if (Y == user->ytrue) {
377: /* First solve is done using true solution to set up problem */
378: KSPSetTolerances(user->solver,1e-8,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
379: } else {
380: KSPSetTolerances(user->solver,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
381: }
382: if (user->ns == 1) {
383: KSPSolve(user->solver,X,Y);
384: KSPGetIterationNumber(user->solver,&its);
385: user->ksp_its+=its;
386: } else {
387: for (i=0;i<user->ns;i++) {
388: Scatter(X,user->subq,user->yi_scatter[i],0,0);
389: Scatter(Y,user->suby,user->yi_scatter[i],0,0);
390: KSPSolve(user->solver,user->subq,user->suby);
391: KSPGetIterationNumber(user->solver,&its);
392: user->ksp_its+=its;
393: Gather(Y,user->suby,user->yi_scatter[i],0,0);
394: }
395: }
396: return(0);
397: }
400: PetscErrorCode QMatMult(Mat J_shell, Vec X, Vec Y)
401: {
403: AppCtx *user;
404: PetscInt i;
407: MatShellGetContext(J_shell,(void**)&user);
408: if (user->ns == 1) {
409: MatMult(user->Q,X,Y);
410: } else {
411: for (i=0;i<user->ns;i++) {
412: Scatter(X,user->subq,user->yi_scatter[i],0,0);
413: Scatter(Y,user->subd,user->di_scatter[i],0,0);
414: MatMult(user->Q,user->subq,user->subd);
415: Gather(Y,user->subd,user->di_scatter[i],0,0);
416: }
417: }
418: return(0);
419: }
423: PetscErrorCode QMatMultTranspose(Mat J_shell, Vec X, Vec Y)
424: {
426: AppCtx *user;
427: PetscInt i;
430: MatShellGetContext(J_shell,(void**)&user);
431: if (user->ns == 1) {
432: MatMultTranspose(user->Q,X,Y);
433: } else {
434: for (i=0;i<user->ns;i++) {
435: Scatter(X,user->subd,user->di_scatter[i],0,0);
436: Scatter(Y,user->suby,user->yi_scatter[i],0,0);
437: MatMultTranspose(user->Q,user->subd,user->suby);
438: Gather(Y,user->suby,user->yi_scatter[i],0,0);
439: }
440: }
441: return(0);
442: }
446: PetscErrorCode DesignMatMult(Mat J_shell, Vec X, Vec Y)
447: {
449: PetscInt i;
450: AppCtx *user;
453: MatShellGetContext(J_shell,(void**)&user);
455: /* sdiag(1./v) */
456: VecSet(user->uwork,0);
457: VecAXPY(user->uwork,-1.0,user->u);
458: VecExp(user->uwork);
460: /* sdiag(1./((Av*(1./v)).^2)) */
461: MatMult(user->Av,user->uwork,user->Swork);
462: VecPointwiseMult(user->Swork,user->Swork,user->Swork);
463: VecReciprocal(user->Swork);
465: /* (Av * (sdiag(1./v) * b)) */
466: VecPointwiseMult(user->uwork,user->uwork,X);
467: MatMult(user->Av,user->uwork,user->Twork);
469: /* (sdiag(1./((Av*(1./v)).^2)) * (Av * (sdiag(1./v) * b))) */
470: VecPointwiseMult(user->Swork,user->Twork,user->Swork);
472: if (user->ns == 1) {
473: /* (sdiag(Grad*y(:,i)) */
474: MatMult(user->Grad,user->y,user->Twork);
476: /* Div * (sdiag(Grad*y(:,i)) * (sdiag(1./((Av*(1./v)).^2)) * (Av * (sdiag(1./v) * b)))) */
477: VecPointwiseMult(user->Swork,user->Twork,user->Swork);
478: MatMultTranspose(user->Grad,user->Swork,Y);
479: } else {
480: for (i=0;i<user->ns;i++) {
481: Scatter(user->y,user->suby,user->yi_scatter[i],0,0);
482: Scatter(Y,user->subq,user->yi_scatter[i],0,0);
484: MatMult(user->Grad,user->suby,user->Twork);
485: VecPointwiseMult(user->Twork,user->Twork,user->Swork);
486: MatMultTranspose(user->Grad,user->Twork,user->subq);
487: Gather(user->y,user->suby,user->yi_scatter[i],0,0);
488: Gather(Y,user->subq,user->yi_scatter[i],0,0);
489: }
490: }
491: return(0);
492: }
496: PetscErrorCode DesignMatMultTranspose(Mat J_shell, Vec X, Vec Y)
497: {
499: PetscInt i;
500: AppCtx *user;
503: MatShellGetContext(J_shell,(void**)&user);
504: VecZeroEntries(Y);
506: /* Sdiag = 1./((Av*(1./v)).^2) */
507: VecSet(user->uwork,0);
508: VecAXPY(user->uwork,-1.0,user->u);
509: VecExp(user->uwork);
510: MatMult(user->Av,user->uwork,user->Swork);
511: VecPointwiseMult(user->Sdiag,user->Swork,user->Swork);
512: VecReciprocal(user->Sdiag);
514: for (i=0;i<user->ns;i++) {
515: Scatter(X,user->subq,user->yi_scatter[i],0,0);
516: Scatter(user->y,user->suby,user->yi_scatter[i],0,0);
518: /* Swork = (Div' * b(:,i)) */
519: MatMult(user->Grad,user->subq,user->Swork);
521: /* Twork = Grad*y(:,i) */
522: MatMult(user->Grad,user->suby,user->Twork);
524: /* Twork = sdiag(Twork) * Swork */
525: VecPointwiseMult(user->Twork,user->Swork,user->Twork);
528: /* Swork = pointwisemult(Sdiag,Twork) */
529: VecPointwiseMult(user->Swork,user->Twork,user->Sdiag);
531: /* Ywork = Av' * Swork */
532: MatMultTranspose(user->Av,user->Swork,user->Ywork);
534: /* Ywork = pointwisemult(uwork,Ywork) */
535: VecPointwiseMult(user->Ywork,user->uwork,user->Ywork);
536: VecAXPY(Y,1.0,user->Ywork);
537: Gather(user->y,user->suby,user->yi_scatter[i],0,0);
538: }
539: return(0);
540: }
544: PetscErrorCode FormConstraints(Tao tao, Vec X, Vec C, void *ptr)
545: {
546: /* C=Ay - q A = Div * Sigma * Grad + hx*hx*hx*ones(n,n) */
548: PetscReal sum;
549: PetscInt i;
550: AppCtx *user = (AppCtx*)ptr;
553: Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);
554: if (user->ns == 1) {
555: MatMult(user->Grad,user->y,user->Swork);
556: VecPointwiseDivide(user->Swork,user->Swork,user->Av_u);
557: MatMultTranspose(user->Grad,user->Swork,C);
558: VecSum(user->y,&sum);
559: sum /= user->ndesign;
560: VecShift(C,sum);
561: } else {
562: for (i=0;i<user->ns;i++) {
563: Scatter(user->y,user->suby,user->yi_scatter[i],0,0);
564: Scatter(C,user->subq,user->yi_scatter[i],0,0);
565: MatMult(user->Grad,user->suby,user->Swork);
566: VecPointwiseDivide(user->Swork,user->Swork,user->Av_u);
567: MatMultTranspose(user->Grad,user->Swork,user->subq);
569: VecSum(user->suby,&sum);
570: sum /= user->ndesign;
571: VecShift(user->subq,sum);
573: Gather(user->y,user->suby,user->yi_scatter[i],0,0);
574: Gather(C,user->subq,user->yi_scatter[i],0,0);
575: }
576: }
577: VecAXPY(C,-1.0,user->q);
578: return(0);
579: }
583: PetscErrorCode Scatter(Vec x, Vec sub1, VecScatter scat1, Vec sub2, VecScatter scat2)
584: {
588: VecScatterBegin(scat1,x,sub1,INSERT_VALUES,SCATTER_FORWARD);
589: VecScatterEnd(scat1,x,sub1,INSERT_VALUES,SCATTER_FORWARD);
590: if (sub2) {
591: VecScatterBegin(scat2,x,sub2,INSERT_VALUES,SCATTER_FORWARD);
592: VecScatterEnd(scat2,x,sub2,INSERT_VALUES,SCATTER_FORWARD);
593: }
594: return(0);
595: }
599: PetscErrorCode Gather(Vec x, Vec sub1, VecScatter scat1, Vec sub2, VecScatter scat2)
600: {
604: VecScatterBegin(scat1,sub1,x,INSERT_VALUES,SCATTER_REVERSE);
605: VecScatterEnd(scat1,sub1,x,INSERT_VALUES,SCATTER_REVERSE);
606: if (sub2) {
607: VecScatterBegin(scat2,sub2,x,INSERT_VALUES,SCATTER_REVERSE);
608: VecScatterEnd(scat2,sub2,x,INSERT_VALUES,SCATTER_REVERSE);
609: }
610: return(0);
611: }
615: PetscErrorCode EllipticInitialize(AppCtx *user)
616: {
618: PetscInt m,n,i,j,k,l,linear_index,is,js,ks,ls,istart,iend,iblock;
619: Vec XX,YY,ZZ,XXwork,YYwork,ZZwork,UTwork;
620: PetscReal *x,*y,*z;
621: PetscReal h,meanut;
622: PetscScalar hinv,neg_hinv,half = 0.5,sqrt_beta;
623: PetscInt im,indx1,indx2,indy1,indy2,indz1,indz2,nx,ny,nz;
624: IS is_alldesign,is_allstate;
625: IS is_from_d;
626: IS is_from_y;
627: PetscInt lo,hi,hi2,lo2,ysubnlocal,dsubnlocal;
628: const PetscInt *ranges, *subranges;
629: PetscMPIInt size;
630: PetscReal xri,yri,zri,xim,yim,zim,dx1,dx2,dy1,dy2,dz1,dz2,Dx,Dy,Dz;
631: PetscScalar v,vx,vy,vz;
632: PetscInt offset,subindex,subvec,nrank,kk;
634: PetscScalar xr[64] = {0.4970, 0.8498, 0.7814, 0.6268, 0.7782, 0.6402, 0.3617, 0.3160,
635: 0.3610, 0.5298, 0.6987, 0.3331, 0.7962, 0.5596, 0.3866, 0.6774,
636: 0.5407, 0.4518, 0.6702, 0.6061, 0.7580, 0.8997, 0.5198, 0.8326,
637: 0.2138, 0.9198, 0.3000, 0.2833, 0.8288, 0.7076, 0.1820, 0.0728,
638: 0.8447, 0.2367, 0.3239, 0.6413, 0.3114, 0.4731, 0.1192, 0.9273,
639: 0.5724, 0.4331, 0.5136, 0.3547, 0.4413, 0.2602, 0.5698, 0.7278,
640: 0.5261, 0.6230, 0.2454, 0.3948, 0.7479, 0.6582, 0.4660, 0.5594,
641: 0.7574, 0.1143, 0.5900, 0.1065, 0.4260, 0.3294, 0.8276, 0.0756};
643: PetscScalar yr[64] = {0.7345, 0.9120, 0.9288, 0.7528, 0.4463, 0.4985, 0.2497, 0.6256,
644: 0.3425, 0.9026, 0.6983, 0.4230, 0.7140, 0.2970, 0.4474, 0.8792,
645: 0.6604, 0.2485, 0.7968, 0.6127, 0.1796, 0.2437, 0.5938, 0.6137,
646: 0.3867, 0.5658, 0.4575, 0.1009, 0.0863, 0.3361, 0.0738, 0.3985,
647: 0.6602, 0.1437, 0.0934, 0.5983, 0.5950, 0.0763, 0.0768, 0.2288,
648: 0.5761, 0.1129, 0.3841, 0.6150, 0.6904, 0.6686, 0.1361, 0.4601,
649: 0.4491, 0.3716, 0.1969, 0.6537, 0.6743, 0.6991, 0.4811, 0.5480,
650: 0.1684, 0.4569, 0.6889, 0.8437, 0.3015, 0.2854, 0.8199, 0.2658};
652: PetscScalar zr[64] = {0.7668, 0.8573, 0.2654, 0.2719, 0.1060, 0.1311, 0.6232, 0.2295,
653: 0.8009, 0.2147, 0.2119, 0.9325, 0.4473, 0.3600, 0.3374, 0.3819,
654: 0.4066, 0.5801, 0.1673, 0.0959, 0.4638, 0.8236, 0.8800, 0.2939,
655: 0.2028, 0.8262, 0.2706, 0.6276, 0.9085, 0.6443, 0.8241, 0.0712,
656: 0.1824, 0.7789, 0.4389, 0.8415, 0.7055, 0.6639, 0.3653, 0.2078,
657: 0.1987, 0.2297, 0.4321, 0.8115, 0.4915, 0.7764, 0.4657, 0.4627,
658: 0.4569, 0.4232, 0.8514, 0.0674, 0.3227, 0.1055, 0.6690, 0.6313,
659: 0.9226, 0.5461, 0.4126, 0.2364, 0.6096, 0.7042, 0.3914, 0.0711};
662: MPI_Comm_size(PETSC_COMM_WORLD,&size);
663: PetscLogStageRegister("Elliptic Setup",&user->stages[0]);
664: PetscLogStagePush(user->stages[0]);
666: /* Create u,y,c,x */
667: VecCreate(PETSC_COMM_WORLD,&user->u);
668: VecCreate(PETSC_COMM_WORLD,&user->y);
669: VecCreate(PETSC_COMM_WORLD,&user->c);
670: VecSetSizes(user->u,PETSC_DECIDE,user->ndesign);
671: VecSetFromOptions(user->u);
672: VecGetLocalSize(user->u,&ysubnlocal);
673: VecSetSizes(user->y,ysubnlocal*user->ns,user->nstate);
674: VecSetSizes(user->c,ysubnlocal*user->ns,user->m);
675: VecSetFromOptions(user->y);
676: VecSetFromOptions(user->c);
678: /*
679: *******************************
680: Create scatters for x <-> y,u
681: *******************************
683: If the state vector y and design vector u are partitioned as
684: [y_1; y_2; ...; y_np] and [u_1; u_2; ...; u_np] (with np = # of processors),
685: then the solution vector x is organized as
686: [y_1; u_1; y_2; u_2; ...; y_np; u_np].
687: The index sets user->s_is and user->d_is correspond to the indices of the
688: state and design variables owned by the current processor.
689: */
690: VecCreate(PETSC_COMM_WORLD,&user->x);
692: VecGetOwnershipRange(user->y,&lo,&hi);
693: VecGetOwnershipRange(user->u,&lo2,&hi2);
695: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo,1,&is_allstate);
696: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo+lo2,1,&user->s_is);
697: ISCreateStride(PETSC_COMM_SELF,hi2-lo2,lo2,1,&is_alldesign);
698: ISCreateStride(PETSC_COMM_SELF,hi2-lo2,hi+lo2,1,&user->d_is);
700: VecSetSizes(user->x,hi-lo+hi2-lo2,user->n);
701: VecSetFromOptions(user->x);
703: VecScatterCreate(user->x,user->s_is,user->y,is_allstate,&user->state_scatter);
704: VecScatterCreate(user->x,user->d_is,user->u,is_alldesign,&user->design_scatter);
705: ISDestroy(&is_alldesign);
706: ISDestroy(&is_allstate);
707: /*
708: *******************************
709: Create scatter from y to y_1,y_2,...,y_ns
710: *******************************
711: */
712: PetscMalloc1(user->ns,&user->yi_scatter);
713: VecDuplicate(user->u,&user->suby);
714: VecDuplicate(user->u,&user->subq);
716: VecGetOwnershipRange(user->y,&lo2,&hi2);
717: istart = 0;
718: for (i=0; i<user->ns; i++){
719: VecGetOwnershipRange(user->suby,&lo,&hi);
720: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo2+istart,1,&is_from_y);
721: VecScatterCreate(user->y,is_from_y,user->suby,NULL,&user->yi_scatter[i]);
722: istart = istart + hi-lo;
723: ISDestroy(&is_from_y);
724: }
725: /*
726: *******************************
727: Create scatter from d to d_1,d_2,...,d_ns
728: *******************************
729: */
730: VecCreate(PETSC_COMM_WORLD,&user->subd);
731: VecSetSizes(user->subd,PETSC_DECIDE,user->ndata);
732: VecSetFromOptions(user->subd);
733: VecCreate(PETSC_COMM_WORLD,&user->d);
734: VecGetLocalSize(user->subd,&dsubnlocal);
735: VecSetSizes(user->d,dsubnlocal*user->ns,user->ndata*user->ns);
736: VecSetFromOptions(user->d);
737: PetscMalloc1(user->ns,&user->di_scatter);
739: VecGetOwnershipRange(user->d,&lo2,&hi2);
740: istart = 0;
741: for (i=0; i<user->ns; i++){
742: VecGetOwnershipRange(user->subd,&lo,&hi);
743: ISCreateStride(PETSC_COMM_SELF,hi-lo,lo2+istart,1,&is_from_d);
744: VecScatterCreate(user->d,is_from_d,user->subd,NULL,&user->di_scatter[i]);
745: istart = istart + hi-lo;
746: ISDestroy(&is_from_d);
747: }
749: PetscMalloc1(user->mx,&x);
750: PetscMalloc1(user->mx,&y);
751: PetscMalloc1(user->mx,&z);
753: user->ksp_its = 0;
754: user->ksp_its_initial = 0;
756: n = user->mx * user->mx * user->mx;
757: m = 3 * user->mx * user->mx * (user->mx-1);
758: sqrt_beta = PetscSqrtScalar(user->beta);
760: VecCreate(PETSC_COMM_WORLD,&XX);
761: VecCreate(PETSC_COMM_WORLD,&user->q);
762: VecSetSizes(XX,ysubnlocal,n);
763: VecSetSizes(user->q,ysubnlocal*user->ns,user->m);
764: VecSetFromOptions(XX);
765: VecSetFromOptions(user->q);
767: VecDuplicate(XX,&YY);
768: VecDuplicate(XX,&ZZ);
769: VecDuplicate(XX,&XXwork);
770: VecDuplicate(XX,&YYwork);
771: VecDuplicate(XX,&ZZwork);
772: VecDuplicate(XX,&UTwork);
773: VecDuplicate(XX,&user->utrue);
775: /* map for striding q */
776: VecGetOwnershipRanges(user->q,&ranges);
777: VecGetOwnershipRanges(user->u,&subranges);
779: VecGetOwnershipRange(user->q,&lo2,&hi2);
780: VecGetOwnershipRange(user->u,&lo,&hi);
781: /* Generate 3D grid, and collect ns (1<=ns<=8) right-hand-side vectors into user->q */
782: h = 1.0/user->mx;
783: hinv = user->mx;
784: neg_hinv = -hinv;
786: VecGetOwnershipRange(XX,&istart,&iend);
787: for (linear_index=istart; linear_index<iend; linear_index++){
788: i = linear_index % user->mx;
789: j = ((linear_index-i)/user->mx) % user->mx;
790: k = ((linear_index-i)/user->mx-j) / user->mx;
791: vx = h*(i+0.5);
792: vy = h*(j+0.5);
793: vz = h*(k+0.5);
794: VecSetValues(XX,1,&linear_index,&vx,INSERT_VALUES);
795: VecSetValues(YY,1,&linear_index,&vy,INSERT_VALUES);
796: VecSetValues(ZZ,1,&linear_index,&vz,INSERT_VALUES);
797: for (is=0; is<2; is++){
798: for (js=0; js<2; js++){
799: for (ks=0; ks<2; ks++){
800: ls = is*4 + js*2 + ks;
801: if (ls<user->ns){
802: l =ls*n + linear_index;
803: /* remap */
804: subindex = l%n;
805: subvec = l/n;
806: nrank=0;
807: while (subindex >= subranges[nrank+1]) nrank++;
808: offset = subindex - subranges[nrank];
809: istart=0;
810: for (kk=0;kk<nrank;kk++) istart+=user->ns*(subranges[kk+1]-subranges[kk]);
811: istart += (subranges[nrank+1]-subranges[nrank])*subvec;
812: l = istart+offset;
813: v = 100*PetscSinScalar(2*PETSC_PI*(vx+0.25*is))*sin(2*PETSC_PI*(vy+0.25*js))*sin(2*PETSC_PI*(vz+0.25*ks));
814: VecSetValues(user->q,1,&l,&v,INSERT_VALUES);
815: }
816: }
817: }
818: }
819: }
821: VecAssemblyBegin(XX);
822: VecAssemblyEnd(XX);
823: VecAssemblyBegin(YY);
824: VecAssemblyEnd(YY);
825: VecAssemblyBegin(ZZ);
826: VecAssemblyEnd(ZZ);
827: VecAssemblyBegin(user->q);
828: VecAssemblyEnd(user->q);
830: /* Compute true parameter function
831: 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) */
832: VecCopy(XX,XXwork);
833: VecCopy(YY,YYwork);
834: VecCopy(ZZ,ZZwork);
836: VecShift(XXwork,-0.25);
837: VecShift(YYwork,-0.25);
838: VecShift(ZZwork,-0.25);
840: VecPointwiseMult(XXwork,XXwork,XXwork);
841: VecPointwiseMult(YYwork,YYwork,YYwork);
842: VecPointwiseMult(ZZwork,ZZwork,ZZwork);
844: VecCopy(XXwork,UTwork);
845: VecAXPY(UTwork,1.0,YYwork);
846: VecAXPY(UTwork,1.0,ZZwork);
847: VecScale(UTwork,-20.0);
848: VecExp(UTwork);
849: VecCopy(UTwork,user->utrue);
851: VecCopy(XX,XXwork);
852: VecCopy(YY,YYwork);
853: VecCopy(ZZ,ZZwork);
855: VecShift(XXwork,-0.75);
856: VecShift(YYwork,-0.75);
857: VecShift(ZZwork,-0.75);
859: VecPointwiseMult(XXwork,XXwork,XXwork);
860: VecPointwiseMult(YYwork,YYwork,YYwork);
861: VecPointwiseMult(ZZwork,ZZwork,ZZwork);
863: VecCopy(XXwork,UTwork);
864: VecAXPY(UTwork,1.0,YYwork);
865: VecAXPY(UTwork,1.0,ZZwork);
866: VecScale(UTwork,-20.0);
867: VecExp(UTwork);
869: VecAXPY(user->utrue,-1.0,UTwork);
871: VecDestroy(&XX);
872: VecDestroy(&YY);
873: VecDestroy(&ZZ);
874: VecDestroy(&XXwork);
875: VecDestroy(&YYwork);
876: VecDestroy(&ZZwork);
877: VecDestroy(&UTwork);
879: /* Initial guess and reference model */
880: VecDuplicate(user->utrue,&user->ur);
881: VecSum(user->utrue,&meanut);
882: meanut = meanut / n;
883: VecSet(user->ur,meanut);
884: VecCopy(user->ur,user->u);
886: /* Generate Grad matrix */
887: MatCreate(PETSC_COMM_WORLD,&user->Grad);
888: MatSetSizes(user->Grad,PETSC_DECIDE,ysubnlocal,m,n);
889: MatSetFromOptions(user->Grad);
890: MatMPIAIJSetPreallocation(user->Grad,2,NULL,2,NULL);
891: MatSeqAIJSetPreallocation(user->Grad,2,NULL);
892: MatGetOwnershipRange(user->Grad,&istart,&iend);
894: for (i=istart; i<iend; i++){
895: if (i<m/3){
896: iblock = i / (user->mx-1);
897: j = iblock*user->mx + (i % (user->mx-1));
898: MatSetValues(user->Grad,1,&i,1,&j,&neg_hinv,INSERT_VALUES);
899: j = j+1;
900: MatSetValues(user->Grad,1,&i,1,&j,&hinv,INSERT_VALUES);
901: }
902: if (i>=m/3 && i<2*m/3){
903: iblock = (i-m/3) / (user->mx*(user->mx-1));
904: j = iblock*user->mx*user->mx + ((i-m/3) % (user->mx*(user->mx-1)));
905: MatSetValues(user->Grad,1,&i,1,&j,&neg_hinv,INSERT_VALUES);
906: j = j + user->mx;
907: MatSetValues(user->Grad,1,&i,1,&j,&hinv,INSERT_VALUES);
908: }
909: if (i>=2*m/3){
910: j = i-2*m/3;
911: MatSetValues(user->Grad,1,&i,1,&j,&neg_hinv,INSERT_VALUES);
912: j = j + user->mx*user->mx;
913: MatSetValues(user->Grad,1,&i,1,&j,&hinv,INSERT_VALUES);
914: }
915: }
917: MatAssemblyBegin(user->Grad,MAT_FINAL_ASSEMBLY);
918: MatAssemblyEnd(user->Grad,MAT_FINAL_ASSEMBLY);
920: /* Generate arithmetic averaging matrix Av */
921: MatCreate(PETSC_COMM_WORLD,&user->Av);
922: MatSetSizes(user->Av,PETSC_DECIDE,ysubnlocal,m,n);
923: MatSetFromOptions(user->Av);
924: MatMPIAIJSetPreallocation(user->Av,2,NULL,2,NULL);
925: MatSeqAIJSetPreallocation(user->Av,2,NULL);
926: MatGetOwnershipRange(user->Av,&istart,&iend);
928: for (i=istart; i<iend; i++){
929: if (i<m/3){
930: iblock = i / (user->mx-1);
931: j = iblock*user->mx + (i % (user->mx-1));
932: MatSetValues(user->Av,1,&i,1,&j,&half,INSERT_VALUES);
933: j = j+1;
934: MatSetValues(user->Av,1,&i,1,&j,&half,INSERT_VALUES);
935: }
936: if (i>=m/3 && i<2*m/3){
937: iblock = (i-m/3) / (user->mx*(user->mx-1));
938: j = iblock*user->mx*user->mx + ((i-m/3) % (user->mx*(user->mx-1)));
939: MatSetValues(user->Av,1,&i,1,&j,&half,INSERT_VALUES);
940: j = j + user->mx;
941: MatSetValues(user->Av,1,&i,1,&j,&half,INSERT_VALUES);
942: }
943: if (i>=2*m/3){
944: j = i-2*m/3;
945: MatSetValues(user->Av,1,&i,1,&j,&half,INSERT_VALUES);
946: j = j + user->mx*user->mx;
947: MatSetValues(user->Av,1,&i,1,&j,&half,INSERT_VALUES);
948: }
949: }
951: MatAssemblyBegin(user->Av,MAT_FINAL_ASSEMBLY);
952: MatAssemblyEnd(user->Av,MAT_FINAL_ASSEMBLY);
954: MatCreate(PETSC_COMM_WORLD,&user->L);
955: MatSetSizes(user->L,PETSC_DECIDE,ysubnlocal,m+n,n);
956: MatSetFromOptions(user->L);
957: MatMPIAIJSetPreallocation(user->L,2,NULL,2,NULL);
958: MatSeqAIJSetPreallocation(user->L,2,NULL);
959: MatGetOwnershipRange(user->L,&istart,&iend);
961: for (i=istart; i<iend; i++){
962: if (i<m/3){
963: iblock = i / (user->mx-1);
964: j = iblock*user->mx + (i % (user->mx-1));
965: MatSetValues(user->L,1,&i,1,&j,&neg_hinv,INSERT_VALUES);
966: j = j+1;
967: MatSetValues(user->L,1,&i,1,&j,&hinv,INSERT_VALUES);
968: }
969: if (i>=m/3 && i<2*m/3){
970: iblock = (i-m/3) / (user->mx*(user->mx-1));
971: j = iblock*user->mx*user->mx + ((i-m/3) % (user->mx*(user->mx-1)));
972: MatSetValues(user->L,1,&i,1,&j,&neg_hinv,INSERT_VALUES);
973: j = j + user->mx;
974: MatSetValues(user->L,1,&i,1,&j,&hinv,INSERT_VALUES);
975: }
976: if (i>=2*m/3 && i<m){
977: j = i-2*m/3;
978: MatSetValues(user->L,1,&i,1,&j,&neg_hinv,INSERT_VALUES);
979: j = j + user->mx*user->mx;
980: MatSetValues(user->L,1,&i,1,&j,&hinv,INSERT_VALUES);
981: }
982: if (i>=m){
983: j = i - m;
984: MatSetValues(user->L,1,&i,1,&j,&sqrt_beta,INSERT_VALUES);
985: }
986: }
987: MatAssemblyBegin(user->L,MAT_FINAL_ASSEMBLY);
988: MatAssemblyEnd(user->L,MAT_FINAL_ASSEMBLY);
989: MatScale(user->L,PetscPowScalar(h,1.5));
991: /* Generate Div matrix */
992: if (!user->use_ptap) {
993: /* Generate Div matrix */
994: MatCreate(PETSC_COMM_WORLD,&user->Div);
995: MatSetSizes(user->Div,ysubnlocal,PETSC_DECIDE,n,m);
996: MatSetFromOptions(user->Div);
997: MatMPIAIJSetPreallocation(user->Div,4,NULL,4,NULL);
998: MatSeqAIJSetPreallocation(user->Div,6,NULL);
999: MatGetOwnershipRange(user->Grad,&istart,&iend);
1001: for (i=istart; i<iend; i++){
1002: if (i<m/3){
1003: iblock = i / (user->mx-1);
1004: j = iblock*user->mx + (i % (user->mx-1));
1005: MatSetValues(user->Div,1,&j,1,&i,&neg_hinv,INSERT_VALUES);
1006: j = j+1;
1007: MatSetValues(user->Div,1,&j,1,&i,&hinv,INSERT_VALUES);
1008: }
1009: if (i>=m/3 && i<2*m/3){
1010: iblock = (i-m/3) / (user->mx*(user->mx-1));
1011: j = iblock*user->mx*user->mx + ((i-m/3) % (user->mx*(user->mx-1)));
1012: MatSetValues(user->Div,1,&j,1,&i,&neg_hinv,INSERT_VALUES);
1013: j = j + user->mx;
1014: MatSetValues(user->Div,1,&j,1,&i,&hinv,INSERT_VALUES);
1015: }
1016: if (i>=2*m/3){
1017: j = i-2*m/3;
1018: MatSetValues(user->Div,1,&j,1,&i,&neg_hinv,INSERT_VALUES);
1019: j = j + user->mx*user->mx;
1020: MatSetValues(user->Div,1,&j,1,&i,&hinv,INSERT_VALUES);
1021: }
1022: }
1024: MatAssemblyBegin(user->Div,MAT_FINAL_ASSEMBLY);
1025: MatAssemblyEnd(user->Div,MAT_FINAL_ASSEMBLY);
1026: MatDuplicate(user->Div,MAT_SHARE_NONZERO_PATTERN,&user->Divwork);
1027: } else {
1028: MatCreate(PETSC_COMM_WORLD,&user->Diag);
1029: MatSetSizes(user->Diag,PETSC_DECIDE,PETSC_DECIDE,m,m);
1030: MatSetFromOptions(user->Diag);
1031: MatMPIAIJSetPreallocation(user->Diag,1,NULL,0,NULL);
1032: MatSeqAIJSetPreallocation(user->Diag,1,NULL);
1033: }
1035: /* Build work vectors and matrices */
1036: VecCreate(PETSC_COMM_WORLD,&user->S);
1037: VecSetSizes(user->S, PETSC_DECIDE, m);
1038: VecSetFromOptions(user->S);
1040: VecCreate(PETSC_COMM_WORLD,&user->lwork);
1041: VecSetSizes(user->lwork,PETSC_DECIDE,m+user->mx*user->mx*user->mx);
1042: VecSetFromOptions(user->lwork);
1044: MatDuplicate(user->Av,MAT_SHARE_NONZERO_PATTERN,&user->Avwork);
1046: VecDuplicate(user->S,&user->Swork);
1047: VecDuplicate(user->S,&user->Sdiag);
1048: VecDuplicate(user->S,&user->Av_u);
1049: VecDuplicate(user->S,&user->Twork);
1050: VecDuplicate(user->y,&user->ywork);
1051: VecDuplicate(user->u,&user->Ywork);
1052: VecDuplicate(user->u,&user->uwork);
1053: VecDuplicate(user->u,&user->js_diag);
1054: VecDuplicate(user->c,&user->cwork);
1055: VecDuplicate(user->d,&user->dwork);
1057: /* Create a matrix-free shell user->Jd for computing B*x */
1058: MatCreateShell(PETSC_COMM_WORLD,ysubnlocal*user->ns,ysubnlocal,user->nstate,user->ndesign,user,&user->Jd);
1059: MatShellSetOperation(user->Jd,MATOP_MULT,(void(*)(void))DesignMatMult);
1060: MatShellSetOperation(user->Jd,MATOP_MULT_TRANSPOSE,(void(*)(void))DesignMatMultTranspose);
1062: /* Compute true state function ytrue given utrue */
1063: VecDuplicate(user->y,&user->ytrue);
1065: /* First compute Av_u = Av*exp(-u) */
1066: VecSet(user->uwork, 0);
1067: VecAXPY(user->uwork,-1.0,user->utrue); /* Note: user->utrue */
1068: VecExp(user->uwork);
1069: MatMult(user->Av,user->uwork,user->Av_u);
1071: /* Next form DSG = Div*S*Grad */
1072: VecCopy(user->Av_u,user->Swork);
1073: VecReciprocal(user->Swork);
1074: if (user->use_ptap) {
1075: MatDiagonalSet(user->Diag,user->Swork,INSERT_VALUES);
1076: MatPtAP(user->Diag,user->Grad,MAT_INITIAL_MATRIX,1.0,&user->DSG);
1077: } else {
1078: MatCopy(user->Div,user->Divwork,SAME_NONZERO_PATTERN);
1079: MatDiagonalScale(user->Divwork,NULL,user->Swork);
1080: MatMatMultSymbolic(user->Divwork,user->Grad,1.0,&user->DSG);
1081: MatMatMultNumeric(user->Divwork,user->Grad,user->DSG);
1082: }
1084: MatSetOption(user->DSG,MAT_SYMMETRIC,PETSC_TRUE);
1085: MatSetOption(user->DSG,MAT_SYMMETRY_ETERNAL,PETSC_TRUE);
1087: if (user->use_lrc == PETSC_TRUE) {
1088: v=PetscSqrtReal(1.0 /user->ndesign);
1089: PetscMalloc1(user->ndesign,&user->ones);
1091: for (i=0;i<user->ndesign;i++) {
1092: user->ones[i]=v;
1093: }
1094: MatCreateDense(PETSC_COMM_WORLD,ysubnlocal,PETSC_DECIDE,user->ndesign,1,user->ones,&user->Ones);
1095: MatAssemblyBegin(user->Ones, MAT_FINAL_ASSEMBLY);
1096: MatAssemblyEnd(user->Ones, MAT_FINAL_ASSEMBLY);
1097: MatCreateLRC(user->DSG,user->Ones,user->Ones,&user->JsBlock);
1098: MatSetUp(user->JsBlock);
1099: } else {
1100: /* Create matrix-free shell user->Js for computing (A + h^3*e*e^T)*x */
1101: MatCreateShell(PETSC_COMM_WORLD,ysubnlocal,ysubnlocal,user->ndesign,user->ndesign,user,&user->JsBlock);
1102: MatShellSetOperation(user->JsBlock,MATOP_MULT,(void(*)(void))StateBlockMatMult);
1103: MatShellSetOperation(user->JsBlock,MATOP_MULT_TRANSPOSE,(void(*)(void))StateBlockMatMult);
1104: }
1105: MatSetOption(user->JsBlock,MAT_SYMMETRIC,PETSC_TRUE);
1106: MatSetOption(user->JsBlock,MAT_SYMMETRY_ETERNAL,PETSC_TRUE);
1107: MatCreateShell(PETSC_COMM_WORLD,ysubnlocal*user->ns,ysubnlocal*user->ns,user->nstate,user->nstate,user,&user->Js);
1108: MatShellSetOperation(user->Js,MATOP_MULT,(void(*)(void))StateMatMult);
1109: MatShellSetOperation(user->Js,MATOP_MULT_TRANSPOSE,(void(*)(void))StateMatMult);
1110: MatSetOption(user->Js,MAT_SYMMETRIC,PETSC_TRUE);
1111: MatSetOption(user->Js,MAT_SYMMETRY_ETERNAL,PETSC_TRUE);
1113: MatCreateShell(PETSC_COMM_WORLD,ysubnlocal*user->ns,ysubnlocal*user->ns,user->nstate,user->nstate,user,&user->JsInv);
1114: MatShellSetOperation(user->JsInv,MATOP_MULT,(void(*)(void))StateInvMatMult);
1115: MatShellSetOperation(user->JsInv,MATOP_MULT_TRANSPOSE,(void(*)(void))StateInvMatMult);
1116: MatSetOption(user->JsInv,MAT_SYMMETRIC,PETSC_TRUE);
1117: MatSetOption(user->JsInv,MAT_SYMMETRY_ETERNAL,PETSC_TRUE);
1119: MatSetOption(user->DSG,MAT_SYMMETRIC,PETSC_TRUE);
1120: MatSetOption(user->DSG,MAT_SYMMETRY_ETERNAL,PETSC_TRUE);
1121: /* Now solve for ytrue */
1122: KSPCreate(PETSC_COMM_WORLD,&user->solver);
1123: KSPSetFromOptions(user->solver);
1125: KSPSetOperators(user->solver,user->JsBlock,user->DSG);
1127: MatMult(user->JsInv,user->q,user->ytrue);
1128: /* First compute Av_u = Av*exp(-u) */
1129: VecSet(user->uwork,0);
1130: VecAXPY(user->uwork,-1.0,user->u); /* Note: user->u */
1131: VecExp(user->uwork);
1132: MatMult(user->Av,user->uwork,user->Av_u);
1134: /* Next update DSG = Div*S*Grad with user->u */
1135: VecCopy(user->Av_u,user->Swork);
1136: VecReciprocal(user->Swork);
1137: if (user->use_ptap) {
1138: MatDiagonalSet(user->Diag,user->Swork,INSERT_VALUES);
1139: MatPtAP(user->Diag,user->Grad,MAT_REUSE_MATRIX,1.0,&user->DSG);
1140: } else {
1141: MatCopy(user->Div,user->Divwork,SAME_NONZERO_PATTERN);
1142: MatDiagonalScale(user->Divwork,NULL,user->Av_u);
1143: MatMatMultNumeric(user->Divwork,user->Grad,user->DSG);
1144: }
1146: /* Now solve for y */
1148: MatMult(user->JsInv,user->q,user->y);
1150: user->ksp_its_initial = user->ksp_its;
1151: user->ksp_its = 0;
1152: /* Construct projection matrix Q (blocks) */
1153: MatCreate(PETSC_COMM_WORLD,&user->Q);
1154: MatSetSizes(user->Q,dsubnlocal,ysubnlocal,user->ndata,user->ndesign);
1155: MatSetFromOptions(user->Q);
1156: MatMPIAIJSetPreallocation(user->Q,8,NULL,8,NULL);
1157: MatSeqAIJSetPreallocation(user->Q,8,NULL);
1159: for (i=0; i<user->mx; i++){
1160: x[i] = h*(i+0.5);
1161: y[i] = h*(i+0.5);
1162: z[i] = h*(i+0.5);
1163: }
1164: MatGetOwnershipRange(user->Q,&istart,&iend);
1166: nx = user->mx; ny = user->mx; nz = user->mx;
1167: for (i=istart; i<iend; i++){
1169: xri = xr[i];
1170: im = 0;
1171: xim = x[im];
1172: while (xri>xim && im<nx){
1173: im = im+1;
1174: xim = x[im];
1175: }
1176: indx1 = im-1;
1177: indx2 = im;
1178: dx1 = xri - x[indx1];
1179: dx2 = x[indx2] - xri;
1181: yri = yr[i];
1182: im = 0;
1183: yim = y[im];
1184: while (yri>yim && im<ny){
1185: im = im+1;
1186: yim = y[im];
1187: }
1188: indy1 = im-1;
1189: indy2 = im;
1190: dy1 = yri - y[indy1];
1191: dy2 = y[indy2] - yri;
1193: zri = zr[i];
1194: im = 0;
1195: zim = z[im];
1196: while (zri>zim && im<nz){
1197: im = im+1;
1198: zim = z[im];
1199: }
1200: indz1 = im-1;
1201: indz2 = im;
1202: dz1 = zri - z[indz1];
1203: dz2 = z[indz2] - zri;
1205: Dx = x[indx2] - x[indx1];
1206: Dy = y[indy2] - y[indy1];
1207: Dz = z[indz2] - z[indz1];
1209: j = indx1 + indy1*nx + indz1*nx*ny;
1210: v = (1-dx1/Dx)*(1-dy1/Dy)*(1-dz1/Dz);
1211: MatSetValues(user->Q,1,&i,1,&j,&v,INSERT_VALUES);
1213: j = indx1 + indy1*nx + indz2*nx*ny;
1214: v = (1-dx1/Dx)*(1-dy1/Dy)*(1-dz2/Dz);
1215: MatSetValues(user->Q,1,&i,1,&j,&v,INSERT_VALUES);
1217: j = indx1 + indy2*nx + indz1*nx*ny;
1218: v = (1-dx1/Dx)*(1-dy2/Dy)*(1-dz1/Dz);
1219: MatSetValues(user->Q,1,&i,1,&j,&v,INSERT_VALUES);
1221: j = indx1 + indy2*nx + indz2*nx*ny;
1222: v = (1-dx1/Dx)*(1-dy2/Dy)*(1-dz2/Dz);
1223: MatSetValues(user->Q,1,&i,1,&j,&v,INSERT_VALUES);
1225: j = indx2 + indy1*nx + indz1*nx*ny;
1226: v = (1-dx2/Dx)*(1-dy1/Dy)*(1-dz1/Dz);
1227: MatSetValues(user->Q,1,&i,1,&j,&v,INSERT_VALUES);
1229: j = indx2 + indy1*nx + indz2*nx*ny;
1230: v = (1-dx2/Dx)*(1-dy1/Dy)*(1-dz2/Dz);
1231: MatSetValues(user->Q,1,&i,1,&j,&v,INSERT_VALUES);
1233: j = indx2 + indy2*nx + indz1*nx*ny;
1234: v = (1-dx2/Dx)*(1-dy2/Dy)*(1-dz1/Dz);
1235: MatSetValues(user->Q,1,&i,1,&j,&v,INSERT_VALUES);
1237: j = indx2 + indy2*nx + indz2*nx*ny;
1238: v = (1-dx2/Dx)*(1-dy2/Dy)*(1-dz2/Dz);
1239: MatSetValues(user->Q,1,&i,1,&j,&v,INSERT_VALUES);
1240: }
1242: MatAssemblyBegin(user->Q,MAT_FINAL_ASSEMBLY);
1243: MatAssemblyEnd(user->Q,MAT_FINAL_ASSEMBLY);
1244: /* Create MQ (composed of blocks of Q */
1245: MatCreateShell(PETSC_COMM_WORLD,dsubnlocal*user->ns,PETSC_DECIDE,user->ndata*user->ns,user->nstate,user,&user->MQ);
1246: MatShellSetOperation(user->MQ,MATOP_MULT,(void(*)(void))QMatMult);
1247: MatShellSetOperation(user->MQ,MATOP_MULT_TRANSPOSE,(void(*)(void))QMatMultTranspose);
1249: /* Add noise to the measurement data */
1250: VecSet(user->ywork,1.0);
1251: VecAYPX(user->ywork,user->noise,user->ytrue);
1252: MatMult(user->MQ,user->ywork,user->d);
1254: /* Now that initial conditions have been set, let the user pass tolerance options to the KSP solver */
1255: PetscFree(x);
1256: PetscFree(y);
1257: PetscFree(z);
1258: PetscLogStagePop();
1259: return(0);
1260: }
1264: PetscErrorCode EllipticDestroy(AppCtx *user)
1265: {
1267: PetscInt i;
1270: MatDestroy(&user->DSG);
1271: KSPDestroy(&user->solver);
1272: MatDestroy(&user->Q);
1273: MatDestroy(&user->MQ);
1274: if (!user->use_ptap) {
1275: MatDestroy(&user->Div);
1276: MatDestroy(&user->Divwork);
1277: } else {
1278: MatDestroy(&user->Diag);
1279: }
1280: if (user->use_lrc) {
1281: MatDestroy(&user->Ones);
1282: }
1284: MatDestroy(&user->Grad);
1285: MatDestroy(&user->Av);
1286: MatDestroy(&user->Avwork);
1287: MatDestroy(&user->L);
1288: MatDestroy(&user->Js);
1289: MatDestroy(&user->Jd);
1290: MatDestroy(&user->JsBlock);
1291: MatDestroy(&user->JsInv);
1293: VecDestroy(&user->x);
1294: VecDestroy(&user->u);
1295: VecDestroy(&user->uwork);
1296: VecDestroy(&user->utrue);
1297: VecDestroy(&user->y);
1298: VecDestroy(&user->ywork);
1299: VecDestroy(&user->ytrue);
1300: VecDestroy(&user->c);
1301: VecDestroy(&user->cwork);
1302: VecDestroy(&user->ur);
1303: VecDestroy(&user->q);
1304: VecDestroy(&user->d);
1305: VecDestroy(&user->dwork);
1306: VecDestroy(&user->lwork);
1307: VecDestroy(&user->S);
1308: VecDestroy(&user->Swork);
1309: VecDestroy(&user->Sdiag);
1310: VecDestroy(&user->Ywork);
1311: VecDestroy(&user->Twork);
1312: VecDestroy(&user->Av_u);
1313: VecDestroy(&user->js_diag);
1314: ISDestroy(&user->s_is);
1315: ISDestroy(&user->d_is);
1316: VecDestroy(&user->suby);
1317: VecDestroy(&user->subd);
1318: VecDestroy(&user->subq);
1319: VecScatterDestroy(&user->state_scatter);
1320: VecScatterDestroy(&user->design_scatter);
1321: for (i=0;i<user->ns;i++) {
1322: VecScatterDestroy(&user->yi_scatter[i]);
1323: VecScatterDestroy(&user->di_scatter[i]);
1324: }
1325: PetscFree(user->yi_scatter);
1326: PetscFree(user->di_scatter);
1327: if (user->use_lrc) {
1328: PetscFree(user->ones);
1329: MatDestroy(&user->Ones);
1330: }
1331: return(0);
1332: }
1336: PetscErrorCode EllipticMonitor(Tao tao, void *ptr)
1337: {
1339: Vec X;
1340: PetscReal unorm,ynorm;
1341: AppCtx *user = (AppCtx*)ptr;
1344: TaoGetSolutionVector(tao,&X);
1345: Scatter(X,user->ywork,user->state_scatter,user->uwork,user->design_scatter);
1346: VecAXPY(user->ywork,-1.0,user->ytrue);
1347: VecAXPY(user->uwork,-1.0,user->utrue);
1348: VecNorm(user->uwork,NORM_2,&unorm);
1349: VecNorm(user->ywork,NORM_2,&ynorm);
1350: PetscPrintf(MPI_COMM_WORLD, "||u-ut||=%g ||y-yt||=%g\n",(double)unorm,(double)ynorm);
1351: return(0);
1352: }