Lines Matching refs:user
95 PetscErrorCode ParabolicInitialize(AppCtx *user);
96 PetscErrorCode ParabolicDestroy(AppCtx *user);
121 AppCtx user; in main() local
131 user.mx = 8; in main()
133 …ierr = PetscOptionsInt("-mx","Number of grid points in each direction","",user.mx,&user.mx,NULL);C… in main()
134 user.nt = 8; in main()
135 ierr = PetscOptionsInt("-nt","Number of time steps","",user.nt,&user.nt,NULL);CHKERRQ(ierr); in main()
136 user.ndata = 64; in main()
137 …ierr = PetscOptionsInt("-ndata","Numbers of data points per sample","",user.ndata,&user.ndata,NULL… in main()
138 user.ns = 8; in main()
139 ierr = PetscOptionsInt("-ns","Number of samples","",user.ns,&user.ns,NULL);CHKERRQ(ierr); in main()
140 user.alpha = 1.0; in main()
141 …ierr = PetscOptionsReal("-alpha","Regularization parameter","",user.alpha,&user.alpha,NULL);CHKERR… in main()
142 user.beta = 0.01; in main()
143 …beta","Weight attributed to ||u||^2 in regularization functional","",user.beta,&user.beta,NULL);CH… in main()
144 user.noise = 0.01; in main()
145 …ierr = PetscOptionsReal("-noise","Amount of noise to add to data","",user.noise,&user.noise,NULL);… in main()
149 user.m = user.mx*user.mx*user.mx; /* number of constraints per time step */ in main()
150 user.n = user.m*(user.nt+1); /* number of variables */ in main()
151 user.ht = (PetscReal)1/user.nt; in main()
153 ierr = VecCreate(PETSC_COMM_WORLD,&user.u);CHKERRQ(ierr); in main()
154 ierr = VecCreate(PETSC_COMM_WORLD,&user.y);CHKERRQ(ierr); in main()
155 ierr = VecCreate(PETSC_COMM_WORLD,&user.c);CHKERRQ(ierr); in main()
156 ierr = VecSetSizes(user.u,PETSC_DECIDE,user.n-user.m*user.nt);CHKERRQ(ierr); in main()
157 ierr = VecSetSizes(user.y,PETSC_DECIDE,user.m*user.nt);CHKERRQ(ierr); in main()
158 ierr = VecSetSizes(user.c,PETSC_DECIDE,user.m*user.nt);CHKERRQ(ierr); in main()
159 ierr = VecSetFromOptions(user.u);CHKERRQ(ierr); in main()
160 ierr = VecSetFromOptions(user.y);CHKERRQ(ierr); in main()
161 ierr = VecSetFromOptions(user.c);CHKERRQ(ierr); in main()
173 ierr = VecGetOwnershipRange(user.y,&lo,&hi);CHKERRQ(ierr); in main()
174 ierr = VecGetOwnershipRange(user.u,&lo2,&hi2);CHKERRQ(ierr); in main()
177 ierr = ISCreateStride(PETSC_COMM_SELF,hi-lo,lo+lo2,1,&user.s_is);CHKERRQ(ierr); in main()
180 ierr = ISCreateStride(PETSC_COMM_SELF,hi2-lo2,hi+lo2,1,&user.d_is);CHKERRQ(ierr); in main()
182 ierr = VecSetSizes(x,hi-lo+hi2-lo2,user.n);CHKERRQ(ierr); in main()
185 ierr = VecScatterCreate(x,user.s_is,user.y,is_allstate,&user.state_scatter);CHKERRQ(ierr); in main()
186 ierr = VecScatterCreate(x,user.d_is,user.u,is_alldesign,&user.design_scatter);CHKERRQ(ierr); in main()
195 ierr = ParabolicInitialize(&user);CHKERRQ(ierr); in main()
197 ierr = Gather(x,user.y,user.state_scatter,user.u,user.design_scatter);CHKERRQ(ierr); in main()
203 ierr = TaoSetObjectiveRoutine(tao, FormFunction, (void *)&user);CHKERRQ(ierr); in main()
204 ierr = TaoSetGradientRoutine(tao, FormGradient, (void *)&user);CHKERRQ(ierr); in main()
205 ierr = TaoSetConstraintsRoutine(tao, user.c, FormConstraints, (void *)&user);CHKERRQ(ierr); in main()
207 …ierr = TaoSetJacobianStateRoutine(tao, user.Js, user.JsBlockPrec, user.JsInv, FormJacobianState, (… in main()
208 ierr = TaoSetJacobianDesignRoutine(tao, user.Jd, FormJacobianDesign, (void *)&user);CHKERRQ(ierr); in main()
211 ierr = TaoSetStateDesignIS(tao,user.s_is,user.d_is);CHKERRQ(ierr); in main()
216 user.ksp_its_initial = user.ksp_its; in main()
217 ksp_old = user.ksp_its; in main()
220 ierr = PetscPrintf(PETSC_COMM_WORLD,"KSP Iterations = %D\n",user.ksp_its-ksp_old);CHKERRQ(ierr); in main()
227 ierr = PetscPrintf(PETSC_COMM_WORLD,"%D\n",user.ksp_its_initial);CHKERRQ(ierr); in main()
229 ierr = PetscPrintf(PETSC_COMM_WORLD,"%D\n",user.ksp_its);CHKERRQ(ierr); in main()
231 …ierr = PetscPrintf(PETSC_COMM_WORLD,"%D\n",(user.ksp_its-user.ksp_its_initial)/ntests);CHKERRQ(ier… in main()
236 ierr = ParabolicDestroy(&user);CHKERRQ(ierr); in main()
252 AppCtx *user = (AppCtx*)ptr; in FormFunction() local
255 ierr = Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);CHKERRQ(ierr); in FormFunction()
256 ierr = Scatter_i(user->y,user->yi,user->yi_scatter,user->nt);CHKERRQ(ierr); in FormFunction()
257 for (j=0; j<user->ns; j++){ in FormFunction()
258 i = user->sample_times[j]; in FormFunction()
259 ierr = MatMult(user->Qblock,user->yi[i],user->di[j]);CHKERRQ(ierr); in FormFunction()
261 ierr = Gather_i(user->dwork,user->di,user->di_scatter,user->ns);CHKERRQ(ierr); in FormFunction()
262 ierr = VecAXPY(user->dwork,-1.0,user->d);CHKERRQ(ierr); in FormFunction()
263 ierr = VecDot(user->dwork,user->dwork,&d1);CHKERRQ(ierr); in FormFunction()
265 ierr = VecWAXPY(user->uwork,-1.0,user->ur,user->u);CHKERRQ(ierr); in FormFunction()
266 ierr = MatMult(user->L,user->uwork,user->lwork);CHKERRQ(ierr); in FormFunction()
267 ierr = VecDot(user->lwork,user->lwork,&d2);CHKERRQ(ierr); in FormFunction()
269 *f = 0.5 * (d1 + user->alpha*d2); in FormFunction()
282 AppCtx *user = (AppCtx*)ptr; in FormGradient() local
285 ierr = Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);CHKERRQ(ierr); in FormGradient()
286 ierr = Scatter_i(user->y,user->yi,user->yi_scatter,user->nt);CHKERRQ(ierr); in FormGradient()
287 for (j=0; j<user->ns; j++){ in FormGradient()
288 i = user->sample_times[j]; in FormGradient()
289 ierr = MatMult(user->Qblock,user->yi[i],user->di[j]);CHKERRQ(ierr); in FormGradient()
291 ierr = Gather_i(user->dwork,user->di,user->di_scatter,user->ns);CHKERRQ(ierr); in FormGradient()
292 ierr = VecAXPY(user->dwork,-1.0,user->d);CHKERRQ(ierr); in FormGradient()
293 ierr = Scatter_i(user->dwork,user->di,user->di_scatter,user->ns);CHKERRQ(ierr); in FormGradient()
294 ierr = VecSet(user->ywork,0.0);CHKERRQ(ierr); in FormGradient()
295 ierr = Scatter_i(user->ywork,user->yiwork,user->yi_scatter,user->nt);CHKERRQ(ierr); in FormGradient()
296 for (j=0; j<user->ns; j++){ in FormGradient()
297 i = user->sample_times[j]; in FormGradient()
298 ierr = MatMult(user->QblockT,user->di[j],user->yiwork[i]);CHKERRQ(ierr); in FormGradient()
300 ierr = Gather_i(user->ywork,user->yiwork,user->yi_scatter,user->nt);CHKERRQ(ierr); in FormGradient()
302 ierr = VecWAXPY(user->uwork,-1.0,user->ur,user->u);CHKERRQ(ierr); in FormGradient()
303 ierr = MatMult(user->L,user->uwork,user->lwork);CHKERRQ(ierr); in FormGradient()
304 ierr = MatMult(user->LT,user->lwork,user->uwork);CHKERRQ(ierr); in FormGradient()
305 ierr = VecScale(user->uwork, user->alpha);CHKERRQ(ierr); in FormGradient()
306 ierr = Gather(G,user->ywork,user->state_scatter,user->uwork,user->design_scatter);CHKERRQ(ierr); in FormGradient()
315 AppCtx *user = (AppCtx*)ptr; in FormFunctionGradient() local
318 ierr = Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);CHKERRQ(ierr); in FormFunctionGradient()
319 ierr = Scatter_i(user->y,user->yi,user->yi_scatter,user->nt);CHKERRQ(ierr); in FormFunctionGradient()
320 for (j=0; j<user->ns; j++){ in FormFunctionGradient()
321 i = user->sample_times[j]; in FormFunctionGradient()
322 ierr = MatMult(user->Qblock,user->yi[i],user->di[j]);CHKERRQ(ierr); in FormFunctionGradient()
324 ierr = Gather_i(user->dwork,user->di,user->di_scatter,user->ns);CHKERRQ(ierr); in FormFunctionGradient()
325 ierr = VecAXPY(user->dwork,-1.0,user->d);CHKERRQ(ierr); in FormFunctionGradient()
326 ierr = VecDot(user->dwork,user->dwork,&d1);CHKERRQ(ierr); in FormFunctionGradient()
327 ierr = Scatter_i(user->dwork,user->di,user->di_scatter,user->ns);CHKERRQ(ierr); in FormFunctionGradient()
328 ierr = VecSet(user->ywork,0.0);CHKERRQ(ierr); in FormFunctionGradient()
329 ierr = Scatter_i(user->ywork,user->yiwork,user->yi_scatter,user->nt);CHKERRQ(ierr); in FormFunctionGradient()
330 for (j=0; j<user->ns; j++){ in FormFunctionGradient()
331 i = user->sample_times[j]; in FormFunctionGradient()
332 ierr = MatMult(user->QblockT,user->di[j],user->yiwork[i]);CHKERRQ(ierr); in FormFunctionGradient()
334 ierr = Gather_i(user->ywork,user->yiwork,user->yi_scatter,user->nt);CHKERRQ(ierr); in FormFunctionGradient()
336 ierr = VecWAXPY(user->uwork,-1.0,user->ur,user->u);CHKERRQ(ierr); in FormFunctionGradient()
337 ierr = MatMult(user->L,user->uwork,user->lwork);CHKERRQ(ierr); in FormFunctionGradient()
338 ierr = VecDot(user->lwork,user->lwork,&d2);CHKERRQ(ierr); in FormFunctionGradient()
339 ierr = MatMult(user->LT,user->lwork,user->uwork);CHKERRQ(ierr); in FormFunctionGradient()
340 ierr = VecScale(user->uwork, user->alpha);CHKERRQ(ierr); in FormFunctionGradient()
341 *f = 0.5 * (d1 + user->alpha*d2); in FormFunctionGradient()
343 ierr = Gather(G,user->ywork,user->state_scatter,user->uwork,user->design_scatter);CHKERRQ(ierr); in FormFunctionGradient()
354 AppCtx *user = (AppCtx*)ptr; in FormJacobianState() local
357 ierr = Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);CHKERRQ(ierr); in FormJacobianState()
358 ierr = VecSet(user->uwork,0);CHKERRQ(ierr); in FormJacobianState()
359 ierr = VecAXPY(user->uwork,-1.0,user->u);CHKERRQ(ierr); in FormJacobianState()
360 ierr = VecExp(user->uwork);CHKERRQ(ierr); in FormJacobianState()
361 ierr = MatMult(user->Av,user->uwork,user->Av_u);CHKERRQ(ierr); in FormJacobianState()
362 ierr = VecCopy(user->Av_u,user->Swork);CHKERRQ(ierr); in FormJacobianState()
363 ierr = VecReciprocal(user->Swork);CHKERRQ(ierr); in FormJacobianState()
364 ierr = MatCopy(user->Div,user->Divwork,SAME_NONZERO_PATTERN);CHKERRQ(ierr); in FormJacobianState()
365 ierr = MatDiagonalScale(user->Divwork,NULL,user->Swork);CHKERRQ(ierr); in FormJacobianState()
366 if (user->dsg_formed) { in FormJacobianState()
367 ierr = MatProductNumeric(user->DSG);CHKERRQ(ierr); in FormJacobianState()
369 …ierr = MatMatMult(user->Divwork,user->Grad,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&user->DSG);CHKERRQ(ie… in FormJacobianState()
370 user->dsg_formed = PETSC_TRUE; in FormJacobianState()
374 ierr = MatScale(user->DSG,user->ht);CHKERRQ(ierr); in FormJacobianState()
375 ierr = MatShift(user->DSG,1.0);CHKERRQ(ierr); in FormJacobianState()
384 AppCtx *user = (AppCtx*)ptr; in FormJacobianDesign() local
387 ierr = Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);CHKERRQ(ierr); in FormJacobianDesign()
395 AppCtx *user; in StateMatMult() local
398 ierr = MatShellGetContext(J_shell,(void**)&user);CHKERRQ(ierr); in StateMatMult()
399 ierr = Scatter_i(X,user->yi,user->yi_scatter,user->nt);CHKERRQ(ierr); in StateMatMult()
400 ierr = MatMult(user->JsBlock,user->yi[0],user->yiwork[0]);CHKERRQ(ierr); in StateMatMult()
401 for (i=1; i<user->nt; i++){ in StateMatMult()
402 ierr = MatMult(user->JsBlock,user->yi[i],user->yiwork[i]);CHKERRQ(ierr); in StateMatMult()
403 ierr = VecAXPY(user->yiwork[i],-1.0,user->yi[i-1]);CHKERRQ(ierr); in StateMatMult()
405 ierr = Gather_i(Y,user->yiwork,user->yi_scatter,user->nt);CHKERRQ(ierr); in StateMatMult()
413 AppCtx *user; in StateMatMultTranspose() local
416 ierr = MatShellGetContext(J_shell,(void**)&user);CHKERRQ(ierr); in StateMatMultTranspose()
417 ierr = Scatter_i(X,user->yi,user->yi_scatter,user->nt);CHKERRQ(ierr); in StateMatMultTranspose()
418 for (i=0; i<user->nt-1; i++){ in StateMatMultTranspose()
419 ierr = MatMult(user->JsBlock,user->yi[i],user->yiwork[i]);CHKERRQ(ierr); in StateMatMultTranspose()
420 ierr = VecAXPY(user->yiwork[i],-1.0,user->yi[i+1]);CHKERRQ(ierr); in StateMatMultTranspose()
422 i = user->nt-1; in StateMatMultTranspose()
423 ierr = MatMult(user->JsBlock,user->yi[i],user->yiwork[i]);CHKERRQ(ierr); in StateMatMultTranspose()
424 ierr = Gather_i(Y,user->yiwork,user->yi_scatter,user->nt);CHKERRQ(ierr); in StateMatMultTranspose()
431 AppCtx *user; in StateMatBlockMult() local
434 ierr = MatShellGetContext(J_shell,(void**)&user);CHKERRQ(ierr); in StateMatBlockMult()
435 ierr = MatMult(user->Grad,X,user->Swork);CHKERRQ(ierr); in StateMatBlockMult()
436 ierr = VecPointwiseDivide(user->Swork,user->Swork,user->Av_u);CHKERRQ(ierr); in StateMatBlockMult()
437 ierr = MatMult(user->Div,user->Swork,Y);CHKERRQ(ierr); in StateMatBlockMult()
438 ierr = VecAYPX(Y,user->ht,X);CHKERRQ(ierr); in StateMatBlockMult()
446 AppCtx *user; in DesignMatMult() local
449 ierr = MatShellGetContext(J_shell,(void**)&user);CHKERRQ(ierr); in DesignMatMult()
452 ierr = VecSet(user->uwork,0);CHKERRQ(ierr); in DesignMatMult()
453 ierr = VecAXPY(user->uwork,-1.0,user->u);CHKERRQ(ierr); in DesignMatMult()
454 ierr = VecExp(user->uwork);CHKERRQ(ierr); in DesignMatMult()
457 ierr = MatMult(user->Av,user->uwork,user->Swork);CHKERRQ(ierr); in DesignMatMult()
458 ierr = VecPointwiseMult(user->Swork,user->Swork,user->Swork);CHKERRQ(ierr); in DesignMatMult()
459 ierr = VecReciprocal(user->Swork);CHKERRQ(ierr); in DesignMatMult()
462 ierr = VecPointwiseMult(user->uwork,user->uwork,X);CHKERRQ(ierr); in DesignMatMult()
463 ierr = MatMult(user->Av,user->uwork,user->Twork);CHKERRQ(ierr); in DesignMatMult()
466 ierr = VecPointwiseMult(user->Swork,user->Twork,user->Swork);CHKERRQ(ierr); in DesignMatMult()
468 ierr = Scatter_i(user->y,user->yi,user->yi_scatter,user->nt);CHKERRQ(ierr); in DesignMatMult()
469 for (i=0; i<user->nt; i++){ in DesignMatMult()
471 ierr = MatMult(user->Grad,user->yi[i],user->Twork);CHKERRQ(ierr); in DesignMatMult()
474 ierr = VecPointwiseMult(user->Twork,user->Twork,user->Swork);CHKERRQ(ierr); in DesignMatMult()
475 ierr = MatMult(user->Div,user->Twork,user->yiwork[i]);CHKERRQ(ierr); in DesignMatMult()
476 ierr = VecScale(user->yiwork[i],user->ht);CHKERRQ(ierr); in DesignMatMult()
478 ierr = Gather_i(Y,user->yiwork,user->yi_scatter,user->nt);CHKERRQ(ierr); in DesignMatMult()
487 AppCtx *user; in DesignMatMultTranspose() local
490 ierr = MatShellGetContext(J_shell,(void**)&user);CHKERRQ(ierr); in DesignMatMultTranspose()
493 ierr = VecSet(user->uwork,0);CHKERRQ(ierr); in DesignMatMultTranspose()
494 ierr = VecAXPY(user->uwork,-1.0,user->u);CHKERRQ(ierr); in DesignMatMultTranspose()
495 ierr = VecExp(user->uwork);CHKERRQ(ierr); in DesignMatMultTranspose()
496 ierr = MatMult(user->Av,user->uwork,user->Rwork);CHKERRQ(ierr); in DesignMatMultTranspose()
497 ierr = VecPointwiseMult(user->Rwork,user->Rwork,user->Rwork);CHKERRQ(ierr); in DesignMatMultTranspose()
498 ierr = VecReciprocal(user->Rwork);CHKERRQ(ierr); in DesignMatMultTranspose()
501 ierr = Scatter_i(user->y,user->yi,user->yi_scatter,user->nt);CHKERRQ(ierr); in DesignMatMultTranspose()
502 ierr = Scatter_i(X,user->yiwork,user->yi_scatter,user->nt);CHKERRQ(ierr); in DesignMatMultTranspose()
503 for (i=0; i<user->nt; i++){ in DesignMatMultTranspose()
505 ierr = MatMult(user->Grad,user->yiwork[i],user->Swork);CHKERRQ(ierr); in DesignMatMultTranspose()
508 ierr = MatMult(user->Grad,user->yi[i],user->Twork);CHKERRQ(ierr); in DesignMatMultTranspose()
511 ierr = VecPointwiseMult(user->Twork,user->Swork,user->Twork);CHKERRQ(ierr); in DesignMatMultTranspose()
514 ierr = VecPointwiseMult(user->Twork,user->Rwork,user->Twork);CHKERRQ(ierr); in DesignMatMultTranspose()
517 ierr = MatMult(user->AvT,user->Twork,user->yiwork[i]);CHKERRQ(ierr); in DesignMatMultTranspose()
520 ierr = VecPointwiseMult(user->yiwork[i],user->uwork,user->yiwork[i]);CHKERRQ(ierr); in DesignMatMultTranspose()
521 ierr = VecAXPY(Y,user->ht,user->yiwork[i]);CHKERRQ(ierr); in DesignMatMultTranspose()
529 AppCtx *user; in StateMatBlockPrecMult() local
532 ierr = PCShellGetContext(PC_shell,(void**)&user);CHKERRQ(ierr); in StateMatBlockPrecMult()
534 if (user->dsg_formed) { in StateMatBlockPrecMult()
535 …ierr = MatSOR(user->DSG,X,1.0,(MatSORType)(SOR_ZERO_INITIAL_GUESS | SOR_LOCAL_SYMMETRIC_SWEEP),0.0… in StateMatBlockPrecMult()
543 AppCtx *user; in StateMatInvMult() local
547 ierr = MatShellGetContext(J_shell,(void**)&user);CHKERRQ(ierr); in StateMatInvMult()
549 if (Y == user->ytrue) { in StateMatInvMult()
551 …ierr = KSPSetTolerances(user->solver,1e-8,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr); in StateMatInvMult()
553 …ierr = KSPSetTolerances(user->solver,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);CHKE… in StateMatInvMult()
556 ierr = Scatter_i(X,user->yi,user->yi_scatter,user->nt);CHKERRQ(ierr); in StateMatInvMult()
557 ierr = KSPSolve(user->solver,user->yi[0],user->yiwork[0]);CHKERRQ(ierr); in StateMatInvMult()
558 ierr = KSPGetIterationNumber(user->solver,&its);CHKERRQ(ierr); in StateMatInvMult()
559 user->ksp_its = user->ksp_its + its; in StateMatInvMult()
561 for (i=1; i<user->nt; i++){ in StateMatInvMult()
562 ierr = VecAXPY(user->yi[i],1.0,user->yiwork[i-1]);CHKERRQ(ierr); in StateMatInvMult()
563 ierr = KSPSolve(user->solver,user->yi[i],user->yiwork[i]);CHKERRQ(ierr); in StateMatInvMult()
564 ierr = KSPGetIterationNumber(user->solver,&its);CHKERRQ(ierr); in StateMatInvMult()
565 user->ksp_its = user->ksp_its + its; in StateMatInvMult()
567 ierr = Gather_i(Y,user->yiwork,user->yi_scatter,user->nt);CHKERRQ(ierr); in StateMatInvMult()
574 AppCtx *user; in StateMatInvTransposeMult() local
578 ierr = MatShellGetContext(J_shell,(void**)&user);CHKERRQ(ierr); in StateMatInvTransposeMult()
580 ierr = Scatter_i(X,user->yi,user->yi_scatter,user->nt);CHKERRQ(ierr); in StateMatInvTransposeMult()
582 i = user->nt - 1; in StateMatInvTransposeMult()
583 ierr = KSPSolve(user->solver,user->yi[i],user->yiwork[i]);CHKERRQ(ierr); in StateMatInvTransposeMult()
585 ierr = KSPGetIterationNumber(user->solver,&its);CHKERRQ(ierr); in StateMatInvTransposeMult()
586 user->ksp_its = user->ksp_its + its; in StateMatInvTransposeMult()
588 for (i=user->nt-2; i>=0; i--){ in StateMatInvTransposeMult()
589 ierr = VecAXPY(user->yi[i],1.0,user->yiwork[i+1]);CHKERRQ(ierr); in StateMatInvTransposeMult()
590 ierr = KSPSolve(user->solver,user->yi[i],user->yiwork[i]);CHKERRQ(ierr); in StateMatInvTransposeMult()
592 ierr = KSPGetIterationNumber(user->solver,&its);CHKERRQ(ierr); in StateMatInvTransposeMult()
593 user->ksp_its = user->ksp_its + its; in StateMatInvTransposeMult()
597 ierr = Gather_i(Y,user->yiwork,user->yi_scatter,user->nt);CHKERRQ(ierr); in StateMatInvTransposeMult()
604 AppCtx *user; in StateMatDuplicate() local
607 ierr = MatShellGetContext(J_shell,(void**)&user);CHKERRQ(ierr); in StateMatDuplicate()
609 …ierr = MatCreateShell(PETSC_COMM_WORLD,PETSC_DETERMINE,PETSC_DETERMINE,user->m,user->m,user,new_sh… in StateMatDuplicate()
620 AppCtx *user; in StateMatGetDiagonal() local
623 ierr = MatShellGetContext(J_shell,(void**)&user);CHKERRQ(ierr); in StateMatGetDiagonal()
624 ierr = VecCopy(user->js_diag,X);CHKERRQ(ierr); in StateMatGetDiagonal()
639 AppCtx *user = (AppCtx*)ptr; in FormConstraints() local
642 ierr = Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter);CHKERRQ(ierr); in FormConstraints()
643 ierr = Scatter_i(user->y,user->yi,user->yi_scatter,user->nt);CHKERRQ(ierr); in FormConstraints()
644 ierr = MatMult(user->JsBlock,user->yi[0],user->yiwork[0]);CHKERRQ(ierr); in FormConstraints()
645 for (i=1; i<user->nt; i++){ in FormConstraints()
646 ierr = MatMult(user->JsBlock,user->yi[i],user->yiwork[i]);CHKERRQ(ierr); in FormConstraints()
647 ierr = VecAXPY(user->yiwork[i],-1.0,user->yi[i-1]);CHKERRQ(ierr); in FormConstraints()
649 ierr = Gather_i(C,user->yiwork,user->yi_scatter,user->nt);CHKERRQ(ierr); in FormConstraints()
650 ierr = VecAXPY(C,-1.0,user->q);CHKERRQ(ierr); in FormConstraints()
706 PetscErrorCode ParabolicInitialize(AppCtx *user) in ParabolicInitialize() argument
747 ierr = PetscMalloc1(user->mx,&x);CHKERRQ(ierr); in ParabolicInitialize()
748 ierr = PetscMalloc1(user->mx,&y);CHKERRQ(ierr); in ParabolicInitialize()
749 ierr = PetscMalloc1(user->mx,&z);CHKERRQ(ierr); in ParabolicInitialize()
750 user->jformed = PETSC_FALSE; in ParabolicInitialize()
751 user->dsg_formed = PETSC_FALSE; in ParabolicInitialize()
753 n = user->mx * user->mx * user->mx; in ParabolicInitialize()
754 m = 3 * user->mx * user->mx * (user->mx-1); in ParabolicInitialize()
755 sqrt_beta = PetscSqrtScalar(user->beta); in ParabolicInitialize()
757 user->ksp_its = 0; in ParabolicInitialize()
758 user->ksp_its_initial = 0; in ParabolicInitialize()
760 stime = (PetscReal)user->nt/user->ns; in ParabolicInitialize()
761 ierr = PetscMalloc1(user->ns,&user->sample_times);CHKERRQ(ierr); in ParabolicInitialize()
762 for (i=0; i<user->ns; i++){ in ParabolicInitialize()
763 user->sample_times[i] = (PetscInt)(stime*((PetscReal)i+1.0)-0.5); in ParabolicInitialize()
767 ierr = VecCreate(PETSC_COMM_WORLD,&user->q);CHKERRQ(ierr); in ParabolicInitialize()
769 ierr = VecSetSizes(user->q,PETSC_DECIDE,n*user->nt);CHKERRQ(ierr); in ParabolicInitialize()
771 ierr = VecSetFromOptions(user->q);CHKERRQ(ierr); in ParabolicInitialize()
779 ierr = VecDuplicate(XX,&user->utrue);CHKERRQ(ierr); in ParabolicInitialize()
783 h = 1.0/user->mx; in ParabolicInitialize()
784 hinv = user->mx; in ParabolicInitialize()
789 i = linear_index % user->mx; in ParabolicInitialize()
790 j = ((linear_index-i)/user->mx) % user->mx; in ParabolicInitialize()
791 k = ((linear_index-i)/user->mx-j) / user->mx; in ParabolicInitialize()
827 ierr = VecCopy(XXwork,user->utrue);CHKERRQ(ierr); in ParabolicInitialize()
828 ierr = VecAXPY(user->utrue,1.0,YYwork);CHKERRQ(ierr); in ParabolicInitialize()
829 ierr = VecAXPY(user->utrue,1.0,ZZwork);CHKERRQ(ierr); in ParabolicInitialize()
830 ierr = VecScale(user->utrue,-10.0);CHKERRQ(ierr); in ParabolicInitialize()
831 ierr = VecExp(user->utrue);CHKERRQ(ierr); in ParabolicInitialize()
832 ierr = VecShift(user->utrue,0.5);CHKERRQ(ierr); in ParabolicInitialize()
843 ierr = VecDuplicate(user->utrue,&user->ur);CHKERRQ(ierr); in ParabolicInitialize()
844 ierr = VecSet(user->ur,0.0);CHKERRQ(ierr); in ParabolicInitialize()
847 ierr = MatCreate(PETSC_COMM_WORLD,&user->Grad);CHKERRQ(ierr); in ParabolicInitialize()
848 ierr = MatSetSizes(user->Grad,PETSC_DECIDE,PETSC_DECIDE,m,n);CHKERRQ(ierr); in ParabolicInitialize()
849 ierr = MatSetFromOptions(user->Grad);CHKERRQ(ierr); in ParabolicInitialize()
850 ierr = MatMPIAIJSetPreallocation(user->Grad,2,NULL,2,NULL);CHKERRQ(ierr); in ParabolicInitialize()
851 ierr = MatSeqAIJSetPreallocation(user->Grad,2,NULL);CHKERRQ(ierr); in ParabolicInitialize()
852 ierr = MatGetOwnershipRange(user->Grad,&istart,&iend);CHKERRQ(ierr); in ParabolicInitialize()
856 iblock = i / (user->mx-1); in ParabolicInitialize()
857 j = iblock*user->mx + (i % (user->mx-1)); in ParabolicInitialize()
858 ierr = MatSetValues(user->Grad,1,&i,1,&j,&neg_hinv,INSERT_VALUES);CHKERRQ(ierr); in ParabolicInitialize()
860 ierr = MatSetValues(user->Grad,1,&i,1,&j,&hinv,INSERT_VALUES);CHKERRQ(ierr); in ParabolicInitialize()
863 iblock = (i-m/3) / (user->mx*(user->mx-1)); in ParabolicInitialize()
864 j = iblock*user->mx*user->mx + ((i-m/3) % (user->mx*(user->mx-1))); in ParabolicInitialize()
865 ierr = MatSetValues(user->Grad,1,&i,1,&j,&neg_hinv,INSERT_VALUES);CHKERRQ(ierr); in ParabolicInitialize()
866 j = j + user->mx; in ParabolicInitialize()
867 ierr = MatSetValues(user->Grad,1,&i,1,&j,&hinv,INSERT_VALUES);CHKERRQ(ierr); in ParabolicInitialize()
871 ierr = MatSetValues(user->Grad,1,&i,1,&j,&neg_hinv,INSERT_VALUES);CHKERRQ(ierr); in ParabolicInitialize()
872 j = j + user->mx*user->mx; in ParabolicInitialize()
873 ierr = MatSetValues(user->Grad,1,&i,1,&j,&hinv,INSERT_VALUES);CHKERRQ(ierr); in ParabolicInitialize()
877 ierr = MatAssemblyBegin(user->Grad,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); in ParabolicInitialize()
878 ierr = MatAssemblyEnd(user->Grad,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); in ParabolicInitialize()
882 ierr = MatCreate(PETSC_COMM_WORLD,&user->Av);CHKERRQ(ierr); in ParabolicInitialize()
883 ierr = MatSetSizes(user->Av,PETSC_DECIDE,PETSC_DECIDE,m,n);CHKERRQ(ierr); in ParabolicInitialize()
884 ierr = MatSetFromOptions(user->Av);CHKERRQ(ierr); in ParabolicInitialize()
885 ierr = MatMPIAIJSetPreallocation(user->Av,2,NULL,2,NULL);CHKERRQ(ierr); in ParabolicInitialize()
886 ierr = MatSeqAIJSetPreallocation(user->Av,2,NULL);CHKERRQ(ierr); in ParabolicInitialize()
887 ierr = MatGetOwnershipRange(user->Av,&istart,&iend);CHKERRQ(ierr); in ParabolicInitialize()
891 iblock = i / (user->mx-1); in ParabolicInitialize()
892 j = iblock*user->mx + (i % (user->mx-1)); in ParabolicInitialize()
893 ierr = MatSetValues(user->Av,1,&i,1,&j,&half,INSERT_VALUES);CHKERRQ(ierr); in ParabolicInitialize()
895 ierr = MatSetValues(user->Av,1,&i,1,&j,&half,INSERT_VALUES);CHKERRQ(ierr); in ParabolicInitialize()
898 iblock = (i-m/3) / (user->mx*(user->mx-1)); in ParabolicInitialize()
899 j = iblock*user->mx*user->mx + ((i-m/3) % (user->mx*(user->mx-1))); in ParabolicInitialize()
900 ierr = MatSetValues(user->Av,1,&i,1,&j,&half,INSERT_VALUES);CHKERRQ(ierr); in ParabolicInitialize()
901 j = j + user->mx; in ParabolicInitialize()
902 ierr = MatSetValues(user->Av,1,&i,1,&j,&half,INSERT_VALUES);CHKERRQ(ierr); in ParabolicInitialize()
906 ierr = MatSetValues(user->Av,1,&i,1,&j,&half,INSERT_VALUES);CHKERRQ(ierr); in ParabolicInitialize()
907 j = j + user->mx*user->mx; in ParabolicInitialize()
908 ierr = MatSetValues(user->Av,1,&i,1,&j,&half,INSERT_VALUES);CHKERRQ(ierr); in ParabolicInitialize()
912 ierr = MatAssemblyBegin(user->Av,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); in ParabolicInitialize()
913 ierr = MatAssemblyEnd(user->Av,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); in ParabolicInitialize()
916 ierr = MatTranspose(user->Av,MAT_INITIAL_MATRIX,&user->AvT);CHKERRQ(ierr); in ParabolicInitialize()
918 ierr = MatCreate(PETSC_COMM_WORLD,&user->L);CHKERRQ(ierr); in ParabolicInitialize()
919 ierr = MatSetSizes(user->L,PETSC_DECIDE,PETSC_DECIDE,m+n,n);CHKERRQ(ierr); in ParabolicInitialize()
920 ierr = MatSetFromOptions(user->L);CHKERRQ(ierr); in ParabolicInitialize()
921 ierr = MatMPIAIJSetPreallocation(user->L,2,NULL,2,NULL);CHKERRQ(ierr); in ParabolicInitialize()
922 ierr = MatSeqAIJSetPreallocation(user->L,2,NULL);CHKERRQ(ierr); in ParabolicInitialize()
923 ierr = MatGetOwnershipRange(user->L,&istart,&iend);CHKERRQ(ierr); in ParabolicInitialize()
927 iblock = i / (user->mx-1); in ParabolicInitialize()
928 j = iblock*user->mx + (i % (user->mx-1)); in ParabolicInitialize()
929 ierr = MatSetValues(user->L,1,&i,1,&j,&neg_hinv,INSERT_VALUES);CHKERRQ(ierr); in ParabolicInitialize()
931 ierr = MatSetValues(user->L,1,&i,1,&j,&hinv,INSERT_VALUES);CHKERRQ(ierr); in ParabolicInitialize()
934 iblock = (i-m/3) / (user->mx*(user->mx-1)); in ParabolicInitialize()
935 j = iblock*user->mx*user->mx + ((i-m/3) % (user->mx*(user->mx-1))); in ParabolicInitialize()
936 ierr = MatSetValues(user->L,1,&i,1,&j,&neg_hinv,INSERT_VALUES);CHKERRQ(ierr); in ParabolicInitialize()
937 j = j + user->mx; in ParabolicInitialize()
938 ierr = MatSetValues(user->L,1,&i,1,&j,&hinv,INSERT_VALUES);CHKERRQ(ierr); in ParabolicInitialize()
942 ierr = MatSetValues(user->L,1,&i,1,&j,&neg_hinv,INSERT_VALUES);CHKERRQ(ierr); in ParabolicInitialize()
943 j = j + user->mx*user->mx; in ParabolicInitialize()
944 ierr = MatSetValues(user->L,1,&i,1,&j,&hinv,INSERT_VALUES);CHKERRQ(ierr); in ParabolicInitialize()
948 ierr = MatSetValues(user->L,1,&i,1,&j,&sqrt_beta,INSERT_VALUES);CHKERRQ(ierr); in ParabolicInitialize()
952 ierr = MatAssemblyBegin(user->L,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); in ParabolicInitialize()
953 ierr = MatAssemblyEnd(user->L,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); in ParabolicInitialize()
955 ierr = MatScale(user->L,PetscPowScalar(h,1.5));CHKERRQ(ierr); in ParabolicInitialize()
958 ierr = MatTranspose(user->Grad,MAT_INITIAL_MATRIX,&user->Div);CHKERRQ(ierr); in ParabolicInitialize()
961 ierr = VecCreate(PETSC_COMM_WORLD,&user->S);CHKERRQ(ierr); in ParabolicInitialize()
962 ierr = VecSetSizes(user->S, PETSC_DECIDE, user->mx*user->mx*(user->mx-1)*3);CHKERRQ(ierr); in ParabolicInitialize()
963 ierr = VecSetFromOptions(user->S);CHKERRQ(ierr); in ParabolicInitialize()
965 ierr = VecCreate(PETSC_COMM_WORLD,&user->lwork);CHKERRQ(ierr); in ParabolicInitialize()
966 ierr = VecSetSizes(user->lwork,PETSC_DECIDE,m+user->mx*user->mx*user->mx);CHKERRQ(ierr); in ParabolicInitialize()
967 ierr = VecSetFromOptions(user->lwork);CHKERRQ(ierr); in ParabolicInitialize()
969 ierr = MatDuplicate(user->Div,MAT_SHARE_NONZERO_PATTERN,&user->Divwork);CHKERRQ(ierr); in ParabolicInitialize()
970 ierr = MatDuplicate(user->Av,MAT_SHARE_NONZERO_PATTERN,&user->Avwork);CHKERRQ(ierr); in ParabolicInitialize()
972 ierr = VecCreate(PETSC_COMM_WORLD,&user->d);CHKERRQ(ierr); in ParabolicInitialize()
973 ierr = VecSetSizes(user->d,PETSC_DECIDE,user->ndata*user->nt);CHKERRQ(ierr); in ParabolicInitialize()
974 ierr = VecSetFromOptions(user->d);CHKERRQ(ierr); in ParabolicInitialize()
976 ierr = VecDuplicate(user->S,&user->Swork);CHKERRQ(ierr); in ParabolicInitialize()
977 ierr = VecDuplicate(user->S,&user->Av_u);CHKERRQ(ierr); in ParabolicInitialize()
978 ierr = VecDuplicate(user->S,&user->Twork);CHKERRQ(ierr); in ParabolicInitialize()
979 ierr = VecDuplicate(user->S,&user->Rwork);CHKERRQ(ierr); in ParabolicInitialize()
980 ierr = VecDuplicate(user->y,&user->ywork);CHKERRQ(ierr); in ParabolicInitialize()
981 ierr = VecDuplicate(user->u,&user->uwork);CHKERRQ(ierr); in ParabolicInitialize()
982 ierr = VecDuplicate(user->u,&user->js_diag);CHKERRQ(ierr); in ParabolicInitialize()
983 ierr = VecDuplicate(user->c,&user->cwork);CHKERRQ(ierr); in ParabolicInitialize()
984 ierr = VecDuplicate(user->d,&user->dwork);CHKERRQ(ierr); in ParabolicInitialize()
987 …hell(PETSC_COMM_WORLD,PETSC_DETERMINE,PETSC_DETERMINE,user->m*user->nt,user->m*user->nt,user,&user… in ParabolicInitialize()
988 ierr = MatShellSetOperation(user->Js,MATOP_MULT,(void(*)(void))StateMatMult);CHKERRQ(ierr); in ParabolicInitialize()
989 …ierr = MatShellSetOperation(user->Js,MATOP_DUPLICATE,(void(*)(void))StateMatDuplicate);CHKERRQ(ier… in ParabolicInitialize()
990 …ierr = MatShellSetOperation(user->Js,MATOP_MULT_TRANSPOSE,(void(*)(void))StateMatMultTranspose);CH… in ParabolicInitialize()
991 …ierr = MatShellSetOperation(user->Js,MATOP_GET_DIAGONAL,(void(*)(void))StateMatGetDiagonal);CHKERR… in ParabolicInitialize()
994 …ierr = MatCreateShell(PETSC_COMM_WORLD,PETSC_DETERMINE,PETSC_DETERMINE,user->m,user->m,user,&user-… in ParabolicInitialize()
995 …ierr = MatShellSetOperation(user->JsBlock,MATOP_MULT,(void(*)(void))StateMatBlockMult);CHKERRQ(ier… in ParabolicInitialize()
997 …ierr = MatShellSetOperation(user->JsBlock,MATOP_MULT_TRANSPOSE,(void(*)(void))StateMatBlockMult);C… in ParabolicInitialize()
1002 …ierr = MatCreateShell(PETSC_COMM_WORLD,PETSC_DETERMINE,PETSC_DETERMINE,user->m,user->m,user,&user-… in ParabolicInitialize()
1003 …ierr = MatShellSetOperation(user->JsBlockPrec,MATOP_MULT,(void(*)(void))StateMatBlockPrecMult);CHK… in ParabolicInitialize()
1005 …ierr = MatShellSetOperation(user->JsBlockPrec,MATOP_MULT_TRANSPOSE,(void(*)(void))StateMatBlockPre… in ParabolicInitialize()
1006 ierr = MatSetOption(user->JsBlockPrec,MAT_SYMMETRY_ETERNAL,PETSC_TRUE);CHKERRQ(ierr); in ParabolicInitialize()
1009 …eateShell(PETSC_COMM_WORLD,PETSC_DETERMINE,PETSC_DETERMINE,user->m*user->nt,user->m,user,&user->Jd… in ParabolicInitialize()
1010 ierr = MatShellSetOperation(user->Jd,MATOP_MULT,(void(*)(void))DesignMatMult);CHKERRQ(ierr); in ParabolicInitialize()
1011 …ierr = MatShellSetOperation(user->Jd,MATOP_MULT_TRANSPOSE,(void(*)(void))DesignMatMultTranspose);C… in ParabolicInitialize()
1014 …hell(PETSC_COMM_WORLD,PETSC_DETERMINE,PETSC_DETERMINE,user->m*user->nt,user->m*user->nt,user,&user… in ParabolicInitialize()
1015 ierr = MatShellSetOperation(user->JsInv,MATOP_MULT,(void(*)(void))StateMatInvMult);CHKERRQ(ierr); in ParabolicInitialize()
1016 …ierr = MatShellSetOperation(user->JsInv,MATOP_MULT_TRANSPOSE,(void(*)(void))StateMatInvTransposeMu… in ParabolicInitialize()
1019 ierr = KSPCreate(PETSC_COMM_WORLD,&user->solver);CHKERRQ(ierr); in ParabolicInitialize()
1020 ierr = KSPSetType(user->solver,KSPCG);CHKERRQ(ierr); in ParabolicInitialize()
1021 ierr = KSPSetOperators(user->solver,user->JsBlock,user->JsBlockPrec);CHKERRQ(ierr); in ParabolicInitialize()
1022 ierr = KSPSetInitialGuessNonzero(user->solver,PETSC_FALSE);CHKERRQ(ierr); in ParabolicInitialize()
1023 ierr = KSPSetTolerances(user->solver,1e-4,1e-20,1e3,500);CHKERRQ(ierr); in ParabolicInitialize()
1024 ierr = KSPSetFromOptions(user->solver);CHKERRQ(ierr); in ParabolicInitialize()
1025 ierr = KSPGetPC(user->solver,&user->prec);CHKERRQ(ierr); in ParabolicInitialize()
1026 ierr = PCSetType(user->prec,PCSHELL);CHKERRQ(ierr); in ParabolicInitialize()
1028 ierr = PCShellSetApply(user->prec,StateMatBlockPrecMult);CHKERRQ(ierr); in ParabolicInitialize()
1029 ierr = PCShellSetApplyTranspose(user->prec,StateMatBlockPrecMult);CHKERRQ(ierr); in ParabolicInitialize()
1030 ierr = PCShellSetContext(user->prec,user);CHKERRQ(ierr); in ParabolicInitialize()
1033 ierr = PetscMalloc1(user->nt*user->m,&user->yi_scatter);CHKERRQ(ierr); in ParabolicInitialize()
1035 ierr = VecSetSizes(yi,PETSC_DECIDE,user->mx*user->mx*user->mx);CHKERRQ(ierr); in ParabolicInitialize()
1037 ierr = VecDuplicateVecs(yi,user->nt,&user->yi);CHKERRQ(ierr); in ParabolicInitialize()
1038 ierr = VecDuplicateVecs(yi,user->nt,&user->yiwork);CHKERRQ(ierr); in ParabolicInitialize()
1040 ierr = VecGetOwnershipRange(user->y,&lo2,&hi2);CHKERRQ(ierr); in ParabolicInitialize()
1042 for (i=0; i<user->nt; i++){ in ParabolicInitialize()
1043 ierr = VecGetOwnershipRange(user->yi[i],&lo,&hi);CHKERRQ(ierr); in ParabolicInitialize()
1046 …ierr = VecScatterCreate(user->y,is_from_y,user->yi[i],is_to_yi,&user->yi_scatter[i]);CHKERRQ(ierr); in ParabolicInitialize()
1054 ierr = PetscMalloc1(user->ns*user->ndata,&user->di_scatter);CHKERRQ(ierr); in ParabolicInitialize()
1056 ierr = VecSetSizes(di,PETSC_DECIDE,user->ndata);CHKERRQ(ierr); in ParabolicInitialize()
1058 ierr = VecDuplicateVecs(di,user->ns,&user->di);CHKERRQ(ierr); in ParabolicInitialize()
1059 ierr = VecGetOwnershipRange(user->d,&lo2,&hi2);CHKERRQ(ierr); in ParabolicInitialize()
1061 for (i=0; i<user->ns; i++){ in ParabolicInitialize()
1062 ierr = VecGetOwnershipRange(user->di[i],&lo,&hi);CHKERRQ(ierr); in ParabolicInitialize()
1065 …ierr = VecScatterCreate(user->d,is_from_d,user->di[i],is_to_di,&user->di_scatter[i]);CHKERRQ(ierr); in ParabolicInitialize()
1073 ierr = VecCopy(bc,user->yiwork[0]);CHKERRQ(ierr); in ParabolicInitialize()
1074 for (i=1; i<user->nt; i++){ in ParabolicInitialize()
1075 ierr = VecSet(user->yiwork[i],0.0);CHKERRQ(ierr); in ParabolicInitialize()
1077 ierr = Gather_i(user->q,user->yiwork,user->yi_scatter,user->nt);CHKERRQ(ierr); in ParabolicInitialize()
1081 ierr = VecCreate(PETSC_COMM_WORLD,&user->ytrue);CHKERRQ(ierr); in ParabolicInitialize()
1082 ierr = VecSetSizes(user->ytrue,PETSC_DECIDE,n*user->nt);CHKERRQ(ierr); in ParabolicInitialize()
1083 ierr = VecSetFromOptions(user->ytrue);CHKERRQ(ierr); in ParabolicInitialize()
1086 ierr = VecSet(user->uwork,0);CHKERRQ(ierr); in ParabolicInitialize()
1087 ierr = VecAXPY(user->uwork,-1.0,user->utrue);CHKERRQ(ierr); /* Note: user->utrue */ in ParabolicInitialize()
1088 ierr = VecExp(user->uwork);CHKERRQ(ierr); in ParabolicInitialize()
1089 ierr = MatMult(user->Av,user->uwork,user->Av_u);CHKERRQ(ierr); in ParabolicInitialize()
1092 ierr = MatProductCreate(user->Div,user->Grad,NULL,&user->DSG);CHKERRQ(ierr); in ParabolicInitialize()
1093 ierr = MatProductSetType(user->DSG,MATPRODUCT_AB);CHKERRQ(ierr); in ParabolicInitialize()
1094 ierr = MatProductSetAlgorithm(user->DSG,"default");CHKERRQ(ierr); in ParabolicInitialize()
1095 ierr = MatProductSetFill(user->DSG,PETSC_DEFAULT);CHKERRQ(ierr); in ParabolicInitialize()
1096 ierr = MatProductSetFromOptions(user->DSG);CHKERRQ(ierr); in ParabolicInitialize()
1097 ierr = MatProductSymbolic(user->DSG);CHKERRQ(ierr); in ParabolicInitialize()
1099 user->dsg_formed = PETSC_TRUE; in ParabolicInitialize()
1102 ierr = MatCopy(user->Div,user->Divwork,SAME_NONZERO_PATTERN);CHKERRQ(ierr); in ParabolicInitialize()
1103 ierr = MatDiagonalScale(user->Divwork,NULL,user->Av_u);CHKERRQ(ierr); in ParabolicInitialize()
1104 if (user->dsg_formed) { in ParabolicInitialize()
1105 ierr = MatProductNumeric(user->DSG);CHKERRQ(ierr); in ParabolicInitialize()
1107 … ierr = MatMatMult(user->Div,user->Grad,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&user->DSG);CHKERRQ(ierr); in ParabolicInitialize()
1108 user->dsg_formed = PETSC_TRUE; in ParabolicInitialize()
1111 ierr = MatScale(user->DSG,user->ht);CHKERRQ(ierr); in ParabolicInitialize()
1112 ierr = MatShift(user->DSG,1.0);CHKERRQ(ierr); in ParabolicInitialize()
1115 ierr = StateMatInvMult(user->Js,user->q,user->ytrue);CHKERRQ(ierr); in ParabolicInitialize()
1120 ierr = VecSet(user->uwork,0);CHKERRQ(ierr); in ParabolicInitialize()
1121 ierr = VecAXPY(user->uwork,-1.0,user->u);CHKERRQ(ierr); /* Note: user->u */ in ParabolicInitialize()
1122 ierr = VecExp(user->uwork);CHKERRQ(ierr); in ParabolicInitialize()
1123 ierr = MatMult(user->Av,user->uwork,user->Av_u);CHKERRQ(ierr); in ParabolicInitialize()
1126 ierr = MatCopy(user->Div,user->Divwork,SAME_NONZERO_PATTERN);CHKERRQ(ierr); in ParabolicInitialize()
1127 ierr = MatDiagonalScale(user->Divwork,NULL,user->Av_u);CHKERRQ(ierr); in ParabolicInitialize()
1128 if (user->dsg_formed) { in ParabolicInitialize()
1129 ierr = MatProductNumeric(user->DSG);CHKERRQ(ierr); in ParabolicInitialize()
1131 … ierr = MatMatMult(user->Div,user->Grad,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&user->DSG);CHKERRQ(ierr); in ParabolicInitialize()
1133 user->dsg_formed = PETSC_TRUE; in ParabolicInitialize()
1136 ierr = MatScale(user->DSG,user->ht);CHKERRQ(ierr); in ParabolicInitialize()
1137 ierr = MatShift(user->DSG,1.0);CHKERRQ(ierr); in ParabolicInitialize()
1140 ierr = StateMatInvMult(user->Js,user->q,user->y);CHKERRQ(ierr); in ParabolicInitialize()
1143 ierr = MatCreate(PETSC_COMM_WORLD,&user->Qblock);CHKERRQ(ierr); in ParabolicInitialize()
1144 ierr = MatSetSizes(user->Qblock,PETSC_DECIDE,PETSC_DECIDE,user->ndata,n);CHKERRQ(ierr); in ParabolicInitialize()
1145 ierr = MatSetFromOptions(user->Qblock);CHKERRQ(ierr); in ParabolicInitialize()
1146 ierr = MatMPIAIJSetPreallocation(user->Qblock,8,NULL,8,NULL);CHKERRQ(ierr); in ParabolicInitialize()
1147 ierr = MatSeqAIJSetPreallocation(user->Qblock,8,NULL);CHKERRQ(ierr); in ParabolicInitialize()
1149 for (i=0; i<user->mx; i++){ in ParabolicInitialize()
1155 ierr = MatGetOwnershipRange(user->Qblock,&istart,&iend);CHKERRQ(ierr); in ParabolicInitialize()
1156 nx = user->mx; ny = user->mx; nz = user->mx; in ParabolicInitialize()
1200 ierr = MatSetValues(user->Qblock,1,&i,1,&j,&v,INSERT_VALUES);CHKERRQ(ierr); in ParabolicInitialize()
1204 ierr = MatSetValues(user->Qblock,1,&i,1,&j,&v,INSERT_VALUES);CHKERRQ(ierr); in ParabolicInitialize()
1208 ierr = MatSetValues(user->Qblock,1,&i,1,&j,&v,INSERT_VALUES);CHKERRQ(ierr); in ParabolicInitialize()
1212 ierr = MatSetValues(user->Qblock,1,&i,1,&j,&v,INSERT_VALUES);CHKERRQ(ierr); in ParabolicInitialize()
1216 ierr = MatSetValues(user->Qblock,1,&i,1,&j,&v,INSERT_VALUES);CHKERRQ(ierr); in ParabolicInitialize()
1220 ierr = MatSetValues(user->Qblock,1,&i,1,&j,&v,INSERT_VALUES);CHKERRQ(ierr); in ParabolicInitialize()
1224 ierr = MatSetValues(user->Qblock,1,&i,1,&j,&v,INSERT_VALUES);CHKERRQ(ierr); in ParabolicInitialize()
1228 ierr = MatSetValues(user->Qblock,1,&i,1,&j,&v,INSERT_VALUES);CHKERRQ(ierr); in ParabolicInitialize()
1230 ierr = MatAssemblyBegin(user->Qblock,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); in ParabolicInitialize()
1231 ierr = MatAssemblyEnd(user->Qblock,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); in ParabolicInitialize()
1233 ierr = MatTranspose(user->Qblock,MAT_INITIAL_MATRIX,&user->QblockT);CHKERRQ(ierr); in ParabolicInitialize()
1234 ierr = MatTranspose(user->L,MAT_INITIAL_MATRIX,&user->LT);CHKERRQ(ierr); in ParabolicInitialize()
1237 ierr = VecSet(user->ywork,1.0);CHKERRQ(ierr); in ParabolicInitialize()
1238 ierr = VecAYPX(user->ywork,user->noise,user->ytrue);CHKERRQ(ierr); in ParabolicInitialize()
1239 ierr = Scatter_i(user->ywork,user->yiwork,user->yi_scatter,user->nt);CHKERRQ(ierr); in ParabolicInitialize()
1240 for (j=0; j<user->ns; j++){ in ParabolicInitialize()
1241 i = user->sample_times[j]; in ParabolicInitialize()
1242 ierr = MatMult(user->Qblock,user->yiwork[i],user->di[j]);CHKERRQ(ierr); in ParabolicInitialize()
1244 ierr = Gather_i(user->d,user->di,user->di_scatter,user->ns);CHKERRQ(ierr); in ParabolicInitialize()
1247 ierr = KSPSetFromOptions(user->solver);CHKERRQ(ierr); in ParabolicInitialize()
1254 PetscErrorCode ParabolicDestroy(AppCtx *user) in ParabolicDestroy() argument
1260 ierr = MatDestroy(&user->Qblock);CHKERRQ(ierr); in ParabolicDestroy()
1261 ierr = MatDestroy(&user->QblockT);CHKERRQ(ierr); in ParabolicDestroy()
1262 ierr = MatDestroy(&user->Div);CHKERRQ(ierr); in ParabolicDestroy()
1263 ierr = MatDestroy(&user->Divwork);CHKERRQ(ierr); in ParabolicDestroy()
1264 ierr = MatDestroy(&user->Grad);CHKERRQ(ierr); in ParabolicDestroy()
1265 ierr = MatDestroy(&user->Av);CHKERRQ(ierr); in ParabolicDestroy()
1266 ierr = MatDestroy(&user->Avwork);CHKERRQ(ierr); in ParabolicDestroy()
1267 ierr = MatDestroy(&user->AvT);CHKERRQ(ierr); in ParabolicDestroy()
1268 ierr = MatDestroy(&user->DSG);CHKERRQ(ierr); in ParabolicDestroy()
1269 ierr = MatDestroy(&user->L);CHKERRQ(ierr); in ParabolicDestroy()
1270 ierr = MatDestroy(&user->LT);CHKERRQ(ierr); in ParabolicDestroy()
1271 ierr = KSPDestroy(&user->solver);CHKERRQ(ierr); in ParabolicDestroy()
1272 ierr = MatDestroy(&user->Js);CHKERRQ(ierr); in ParabolicDestroy()
1273 ierr = MatDestroy(&user->Jd);CHKERRQ(ierr); in ParabolicDestroy()
1274 ierr = MatDestroy(&user->JsInv);CHKERRQ(ierr); in ParabolicDestroy()
1275 ierr = MatDestroy(&user->JsBlock);CHKERRQ(ierr); in ParabolicDestroy()
1276 ierr = MatDestroy(&user->JsBlockPrec);CHKERRQ(ierr); in ParabolicDestroy()
1277 ierr = VecDestroy(&user->u);CHKERRQ(ierr); in ParabolicDestroy()
1278 ierr = VecDestroy(&user->uwork);CHKERRQ(ierr); in ParabolicDestroy()
1279 ierr = VecDestroy(&user->utrue);CHKERRQ(ierr); in ParabolicDestroy()
1280 ierr = VecDestroy(&user->y);CHKERRQ(ierr); in ParabolicDestroy()
1281 ierr = VecDestroy(&user->ywork);CHKERRQ(ierr); in ParabolicDestroy()
1282 ierr = VecDestroy(&user->ytrue);CHKERRQ(ierr); in ParabolicDestroy()
1283 ierr = VecDestroyVecs(user->nt,&user->yi);CHKERRQ(ierr); in ParabolicDestroy()
1284 ierr = VecDestroyVecs(user->nt,&user->yiwork);CHKERRQ(ierr); in ParabolicDestroy()
1285 ierr = VecDestroyVecs(user->ns,&user->di);CHKERRQ(ierr); in ParabolicDestroy()
1286 ierr = PetscFree(user->yi);CHKERRQ(ierr); in ParabolicDestroy()
1287 ierr = PetscFree(user->yiwork);CHKERRQ(ierr); in ParabolicDestroy()
1288 ierr = PetscFree(user->di);CHKERRQ(ierr); in ParabolicDestroy()
1289 ierr = VecDestroy(&user->c);CHKERRQ(ierr); in ParabolicDestroy()
1290 ierr = VecDestroy(&user->cwork);CHKERRQ(ierr); in ParabolicDestroy()
1291 ierr = VecDestroy(&user->ur);CHKERRQ(ierr); in ParabolicDestroy()
1292 ierr = VecDestroy(&user->q);CHKERRQ(ierr); in ParabolicDestroy()
1293 ierr = VecDestroy(&user->d);CHKERRQ(ierr); in ParabolicDestroy()
1294 ierr = VecDestroy(&user->dwork);CHKERRQ(ierr); in ParabolicDestroy()
1295 ierr = VecDestroy(&user->lwork);CHKERRQ(ierr); in ParabolicDestroy()
1296 ierr = VecDestroy(&user->S);CHKERRQ(ierr); in ParabolicDestroy()
1297 ierr = VecDestroy(&user->Swork);CHKERRQ(ierr); in ParabolicDestroy()
1298 ierr = VecDestroy(&user->Av_u);CHKERRQ(ierr); in ParabolicDestroy()
1299 ierr = VecDestroy(&user->Twork);CHKERRQ(ierr); in ParabolicDestroy()
1300 ierr = VecDestroy(&user->Rwork);CHKERRQ(ierr); in ParabolicDestroy()
1301 ierr = VecDestroy(&user->js_diag);CHKERRQ(ierr); in ParabolicDestroy()
1302 ierr = ISDestroy(&user->s_is);CHKERRQ(ierr); in ParabolicDestroy()
1303 ierr = ISDestroy(&user->d_is);CHKERRQ(ierr); in ParabolicDestroy()
1304 ierr = VecScatterDestroy(&user->state_scatter);CHKERRQ(ierr); in ParabolicDestroy()
1305 ierr = VecScatterDestroy(&user->design_scatter);CHKERRQ(ierr); in ParabolicDestroy()
1306 for (i=0; i<user->nt; i++){ in ParabolicDestroy()
1307 ierr = VecScatterDestroy(&user->yi_scatter[i]);CHKERRQ(ierr); in ParabolicDestroy()
1309 for (i=0; i<user->ns; i++){ in ParabolicDestroy()
1310 ierr = VecScatterDestroy(&user->di_scatter[i]);CHKERRQ(ierr); in ParabolicDestroy()
1312 ierr = PetscFree(user->yi_scatter);CHKERRQ(ierr); in ParabolicDestroy()
1313 ierr = PetscFree(user->di_scatter);CHKERRQ(ierr); in ParabolicDestroy()
1314 ierr = PetscFree(user->sample_times);CHKERRQ(ierr); in ParabolicDestroy()
1323 AppCtx *user = (AppCtx*)ptr; in ParabolicMonitor() local
1327 ierr = Scatter(X,user->ywork,user->state_scatter,user->uwork,user->design_scatter);CHKERRQ(ierr); in ParabolicMonitor()
1328 ierr = VecAXPY(user->ywork,-1.0,user->ytrue);CHKERRQ(ierr); in ParabolicMonitor()
1329 ierr = VecAXPY(user->uwork,-1.0,user->utrue);CHKERRQ(ierr); in ParabolicMonitor()
1330 ierr = VecNorm(user->uwork,NORM_2,&unorm);CHKERRQ(ierr); in ParabolicMonitor()
1331 ierr = VecNorm(user->ywork,NORM_2,&ynorm);CHKERRQ(ierr); in ParabolicMonitor()