Lines Matching refs:user
75 PetscErrorCode ParabolicInitialize(AppCtx *user);
76 PetscErrorCode ParabolicDestroy(AppCtx *user);
100 AppCtx user; in main() local
109 user.mx = 8; in main()
111 …PetscCall(PetscOptionsInt("-mx", "Number of grid points in each direction", "", user.mx, &user.mx,… in main()
112 user.nt = 8; in main()
113 PetscCall(PetscOptionsInt("-nt", "Number of time steps", "", user.nt, &user.nt, NULL)); in main()
114 user.ndata = 64; in main()
115 …PetscCall(PetscOptionsInt("-ndata", "Numbers of data points per sample", "", user.ndata, &user.nda… in main()
116 user.ns = 8; in main()
117 PetscCall(PetscOptionsInt("-ns", "Number of samples", "", user.ns, &user.ns, NULL)); in main()
118 user.alpha = 1.0; in main()
119 …PetscCall(PetscOptionsReal("-alpha", "Regularization parameter", "", user.alpha, &user.alpha, NULL… in main()
120 user.beta = 0.01; in main()
121 …a", "Weight attributed to ||u||^2 in regularization functional", "", user.beta, &user.beta, NULL)); in main()
122 user.noise = 0.01; in main()
123 …PetscCall(PetscOptionsReal("-noise", "Amount of noise to add to data", "", user.noise, &user.noise… in main()
127 user.m = user.mx * user.mx * user.mx; /* number of constraints per time step */ in main()
128 user.n = user.m * (user.nt + 1); /* number of variables */ in main()
129 user.ht = (PetscReal)1 / user.nt; in main()
131 PetscCall(VecCreate(PETSC_COMM_WORLD, &user.u)); in main()
132 PetscCall(VecCreate(PETSC_COMM_WORLD, &user.y)); in main()
133 PetscCall(VecCreate(PETSC_COMM_WORLD, &user.c)); in main()
134 PetscCall(VecSetSizes(user.u, PETSC_DECIDE, user.n - user.m * user.nt)); in main()
135 PetscCall(VecSetSizes(user.y, PETSC_DECIDE, user.m * user.nt)); in main()
136 PetscCall(VecSetSizes(user.c, PETSC_DECIDE, user.m * user.nt)); in main()
137 PetscCall(VecSetFromOptions(user.u)); in main()
138 PetscCall(VecSetFromOptions(user.y)); in main()
139 PetscCall(VecSetFromOptions(user.c)); in main()
151 PetscCall(VecGetOwnershipRange(user.y, &lo, &hi)); in main()
152 PetscCall(VecGetOwnershipRange(user.u, &lo2, &hi2)); in main()
155 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi - lo, lo + lo2, 1, &user.s_is)); in main()
158 PetscCall(ISCreateStride(PETSC_COMM_SELF, hi2 - lo2, hi + lo2, 1, &user.d_is)); in main()
160 PetscCall(VecSetSizes(x, hi - lo + hi2 - lo2, user.n)); in main()
163 PetscCall(VecScatterCreate(x, user.s_is, user.y, is_allstate, &user.state_scatter)); in main()
164 PetscCall(VecScatterCreate(x, user.d_is, user.u, is_alldesign, &user.design_scatter)); in main()
173 PetscCall(ParabolicInitialize(&user)); in main()
175 PetscCall(Gather(x, user.y, user.state_scatter, user.u, user.design_scatter)); in main()
181 PetscCall(TaoSetObjective(tao, FormFunction, &user)); in main()
182 PetscCall(TaoSetGradient(tao, NULL, FormGradient, &user)); in main()
183 PetscCall(TaoSetConstraintsRoutine(tao, user.c, FormConstraints, &user)); in main()
185 …PetscCall(TaoSetJacobianStateRoutine(tao, user.Js, user.JsBlockPrec, user.JsInv, FormJacobianState… in main()
186 PetscCall(TaoSetJacobianDesignRoutine(tao, user.Jd, FormJacobianDesign, &user)); in main()
189 PetscCall(TaoSetStateDesignIS(tao, user.s_is, user.d_is)); in main()
194 user.ksp_its_initial = user.ksp_its; in main()
195 ksp_old = user.ksp_its; in main()
198 …PetscCall(PetscPrintf(PETSC_COMM_WORLD, "KSP Iterations = %" PetscInt_FMT "\n", user.ksp_its - ksp… in main()
205 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "\n", user.ksp_its_initial)); in main()
207 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "\n", user.ksp_its)); in main()
209 …PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "\n", (user.ksp_its - user.ksp_its_initia… in main()
214 PetscCall(ParabolicDestroy(&user)); in main()
229 AppCtx *user = (AppCtx *)ptr; in FormFunction() local
232 PetscCall(Scatter(X, user->y, user->state_scatter, user->u, user->design_scatter)); in FormFunction()
233 PetscCall(Scatter_i(user->y, user->yi, user->yi_scatter, user->nt)); in FormFunction()
234 for (j = 0; j < user->ns; j++) { in FormFunction()
235 i = user->sample_times[j]; in FormFunction()
236 PetscCall(MatMult(user->Qblock, user->yi[i], user->di[j])); in FormFunction()
238 PetscCall(Gather_i(user->dwork, user->di, user->di_scatter, user->ns)); in FormFunction()
239 PetscCall(VecAXPY(user->dwork, -1.0, user->d)); in FormFunction()
240 PetscCall(VecDot(user->dwork, user->dwork, &d1)); in FormFunction()
242 PetscCall(VecWAXPY(user->uwork, -1.0, user->ur, user->u)); in FormFunction()
243 PetscCall(MatMult(user->L, user->uwork, user->lwork)); in FormFunction()
244 PetscCall(VecDot(user->lwork, user->lwork, &d2)); in FormFunction()
246 *f = 0.5 * (d1 + user->alpha * d2); in FormFunction()
258 AppCtx *user = (AppCtx *)ptr; in FormGradient() local
261 PetscCall(Scatter(X, user->y, user->state_scatter, user->u, user->design_scatter)); in FormGradient()
262 PetscCall(Scatter_i(user->y, user->yi, user->yi_scatter, user->nt)); in FormGradient()
263 for (j = 0; j < user->ns; j++) { in FormGradient()
264 i = user->sample_times[j]; in FormGradient()
265 PetscCall(MatMult(user->Qblock, user->yi[i], user->di[j])); in FormGradient()
267 PetscCall(Gather_i(user->dwork, user->di, user->di_scatter, user->ns)); in FormGradient()
268 PetscCall(VecAXPY(user->dwork, -1.0, user->d)); in FormGradient()
269 PetscCall(Scatter_i(user->dwork, user->di, user->di_scatter, user->ns)); in FormGradient()
270 PetscCall(VecSet(user->ywork, 0.0)); in FormGradient()
271 PetscCall(Scatter_i(user->ywork, user->yiwork, user->yi_scatter, user->nt)); in FormGradient()
272 for (j = 0; j < user->ns; j++) { in FormGradient()
273 i = user->sample_times[j]; in FormGradient()
274 PetscCall(MatMult(user->QblockT, user->di[j], user->yiwork[i])); in FormGradient()
276 PetscCall(Gather_i(user->ywork, user->yiwork, user->yi_scatter, user->nt)); in FormGradient()
278 PetscCall(VecWAXPY(user->uwork, -1.0, user->ur, user->u)); in FormGradient()
279 PetscCall(MatMult(user->L, user->uwork, user->lwork)); in FormGradient()
280 PetscCall(MatMult(user->LT, user->lwork, user->uwork)); in FormGradient()
281 PetscCall(VecScale(user->uwork, user->alpha)); in FormGradient()
282 PetscCall(Gather(G, user->ywork, user->state_scatter, user->uwork, user->design_scatter)); in FormGradient()
290 AppCtx *user = (AppCtx *)ptr; in FormFunctionGradient() local
293 PetscCall(Scatter(X, user->y, user->state_scatter, user->u, user->design_scatter)); in FormFunctionGradient()
294 PetscCall(Scatter_i(user->y, user->yi, user->yi_scatter, user->nt)); in FormFunctionGradient()
295 for (j = 0; j < user->ns; j++) { in FormFunctionGradient()
296 i = user->sample_times[j]; in FormFunctionGradient()
297 PetscCall(MatMult(user->Qblock, user->yi[i], user->di[j])); in FormFunctionGradient()
299 PetscCall(Gather_i(user->dwork, user->di, user->di_scatter, user->ns)); in FormFunctionGradient()
300 PetscCall(VecAXPY(user->dwork, -1.0, user->d)); in FormFunctionGradient()
301 PetscCall(VecDot(user->dwork, user->dwork, &d1)); in FormFunctionGradient()
302 PetscCall(Scatter_i(user->dwork, user->di, user->di_scatter, user->ns)); in FormFunctionGradient()
303 PetscCall(VecSet(user->ywork, 0.0)); in FormFunctionGradient()
304 PetscCall(Scatter_i(user->ywork, user->yiwork, user->yi_scatter, user->nt)); in FormFunctionGradient()
305 for (j = 0; j < user->ns; j++) { in FormFunctionGradient()
306 i = user->sample_times[j]; in FormFunctionGradient()
307 PetscCall(MatMult(user->QblockT, user->di[j], user->yiwork[i])); in FormFunctionGradient()
309 PetscCall(Gather_i(user->ywork, user->yiwork, user->yi_scatter, user->nt)); in FormFunctionGradient()
311 PetscCall(VecWAXPY(user->uwork, -1.0, user->ur, user->u)); in FormFunctionGradient()
312 PetscCall(MatMult(user->L, user->uwork, user->lwork)); in FormFunctionGradient()
313 PetscCall(VecDot(user->lwork, user->lwork, &d2)); in FormFunctionGradient()
314 PetscCall(MatMult(user->LT, user->lwork, user->uwork)); in FormFunctionGradient()
315 PetscCall(VecScale(user->uwork, user->alpha)); in FormFunctionGradient()
316 *f = 0.5 * (d1 + user->alpha * d2); in FormFunctionGradient()
318 PetscCall(Gather(G, user->ywork, user->state_scatter, user->uwork, user->design_scatter)); in FormFunctionGradient()
328 AppCtx *user = (AppCtx *)ptr; in FormJacobianState() local
331 PetscCall(Scatter(X, user->y, user->state_scatter, user->u, user->design_scatter)); in FormJacobianState()
332 PetscCall(VecSet(user->uwork, 0)); in FormJacobianState()
333 PetscCall(VecAXPY(user->uwork, -1.0, user->u)); in FormJacobianState()
334 PetscCall(VecExp(user->uwork)); in FormJacobianState()
335 PetscCall(MatMult(user->Av, user->uwork, user->Av_u)); in FormJacobianState()
336 PetscCall(VecCopy(user->Av_u, user->Swork)); in FormJacobianState()
337 PetscCall(VecReciprocal(user->Swork)); in FormJacobianState()
338 PetscCall(MatCopy(user->Div, user->Divwork, SAME_NONZERO_PATTERN)); in FormJacobianState()
339 PetscCall(MatDiagonalScale(user->Divwork, NULL, user->Swork)); in FormJacobianState()
340 if (user->dsg_formed) { in FormJacobianState()
341 PetscCall(MatProductNumeric(user->DSG)); in FormJacobianState()
343 … PetscCall(MatMatMult(user->Divwork, user->Grad, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &user->DSG)); in FormJacobianState()
344 user->dsg_formed = PETSC_TRUE; in FormJacobianState()
348 PetscCall(MatScale(user->DSG, user->ht)); in FormJacobianState()
349 PetscCall(MatShift(user->DSG, 1.0)); in FormJacobianState()
357 AppCtx *user = (AppCtx *)ptr; in FormJacobianDesign() local
360 PetscCall(Scatter(X, user->y, user->state_scatter, user->u, user->design_scatter)); in FormJacobianDesign()
367 AppCtx *user; in StateMatMult() local
370 PetscCall(MatShellGetContext(J_shell, &user)); in StateMatMult()
371 PetscCall(Scatter_i(X, user->yi, user->yi_scatter, user->nt)); in StateMatMult()
372 PetscCall(MatMult(user->JsBlock, user->yi[0], user->yiwork[0])); in StateMatMult()
373 for (i = 1; i < user->nt; i++) { in StateMatMult()
374 PetscCall(MatMult(user->JsBlock, user->yi[i], user->yiwork[i])); in StateMatMult()
375 PetscCall(VecAXPY(user->yiwork[i], -1.0, user->yi[i - 1])); in StateMatMult()
377 PetscCall(Gather_i(Y, user->yiwork, user->yi_scatter, user->nt)); in StateMatMult()
384 AppCtx *user; in StateMatMultTranspose() local
387 PetscCall(MatShellGetContext(J_shell, &user)); in StateMatMultTranspose()
388 PetscCall(Scatter_i(X, user->yi, user->yi_scatter, user->nt)); in StateMatMultTranspose()
389 for (i = 0; i < user->nt - 1; i++) { in StateMatMultTranspose()
390 PetscCall(MatMult(user->JsBlock, user->yi[i], user->yiwork[i])); in StateMatMultTranspose()
391 PetscCall(VecAXPY(user->yiwork[i], -1.0, user->yi[i + 1])); in StateMatMultTranspose()
393 i = user->nt - 1; in StateMatMultTranspose()
394 PetscCall(MatMult(user->JsBlock, user->yi[i], user->yiwork[i])); in StateMatMultTranspose()
395 PetscCall(Gather_i(Y, user->yiwork, user->yi_scatter, user->nt)); in StateMatMultTranspose()
401 AppCtx *user; in StateMatBlockMult() local
404 PetscCall(MatShellGetContext(J_shell, &user)); in StateMatBlockMult()
405 PetscCall(MatMult(user->Grad, X, user->Swork)); in StateMatBlockMult()
406 PetscCall(VecPointwiseDivide(user->Swork, user->Swork, user->Av_u)); in StateMatBlockMult()
407 PetscCall(MatMult(user->Div, user->Swork, Y)); in StateMatBlockMult()
408 PetscCall(VecAYPX(Y, user->ht, X)); in StateMatBlockMult()
415 AppCtx *user; in DesignMatMult() local
418 PetscCall(MatShellGetContext(J_shell, &user)); in DesignMatMult()
421 PetscCall(VecSet(user->uwork, 0)); in DesignMatMult()
422 PetscCall(VecAXPY(user->uwork, -1.0, user->u)); in DesignMatMult()
423 PetscCall(VecExp(user->uwork)); in DesignMatMult()
426 PetscCall(MatMult(user->Av, user->uwork, user->Swork)); in DesignMatMult()
427 PetscCall(VecPointwiseMult(user->Swork, user->Swork, user->Swork)); in DesignMatMult()
428 PetscCall(VecReciprocal(user->Swork)); in DesignMatMult()
431 PetscCall(VecPointwiseMult(user->uwork, user->uwork, X)); in DesignMatMult()
432 PetscCall(MatMult(user->Av, user->uwork, user->Twork)); in DesignMatMult()
435 PetscCall(VecPointwiseMult(user->Swork, user->Twork, user->Swork)); in DesignMatMult()
437 PetscCall(Scatter_i(user->y, user->yi, user->yi_scatter, user->nt)); in DesignMatMult()
438 for (i = 0; i < user->nt; i++) { in DesignMatMult()
440 PetscCall(MatMult(user->Grad, user->yi[i], user->Twork)); in DesignMatMult()
443 PetscCall(VecPointwiseMult(user->Twork, user->Twork, user->Swork)); in DesignMatMult()
444 PetscCall(MatMult(user->Div, user->Twork, user->yiwork[i])); in DesignMatMult()
445 PetscCall(VecScale(user->yiwork[i], user->ht)); in DesignMatMult()
447 PetscCall(Gather_i(Y, user->yiwork, user->yi_scatter, user->nt)); in DesignMatMult()
454 AppCtx *user; in DesignMatMultTranspose() local
457 PetscCall(MatShellGetContext(J_shell, &user)); in DesignMatMultTranspose()
460 PetscCall(VecSet(user->uwork, 0)); in DesignMatMultTranspose()
461 PetscCall(VecAXPY(user->uwork, -1.0, user->u)); in DesignMatMultTranspose()
462 PetscCall(VecExp(user->uwork)); in DesignMatMultTranspose()
463 PetscCall(MatMult(user->Av, user->uwork, user->Rwork)); in DesignMatMultTranspose()
464 PetscCall(VecPointwiseMult(user->Rwork, user->Rwork, user->Rwork)); in DesignMatMultTranspose()
465 PetscCall(VecReciprocal(user->Rwork)); in DesignMatMultTranspose()
468 PetscCall(Scatter_i(user->y, user->yi, user->yi_scatter, user->nt)); in DesignMatMultTranspose()
469 PetscCall(Scatter_i(X, user->yiwork, user->yi_scatter, user->nt)); in DesignMatMultTranspose()
470 for (i = 0; i < user->nt; i++) { in DesignMatMultTranspose()
472 PetscCall(MatMult(user->Grad, user->yiwork[i], user->Swork)); in DesignMatMultTranspose()
475 PetscCall(MatMult(user->Grad, user->yi[i], user->Twork)); in DesignMatMultTranspose()
478 PetscCall(VecPointwiseMult(user->Twork, user->Swork, user->Twork)); in DesignMatMultTranspose()
481 PetscCall(VecPointwiseMult(user->Twork, user->Rwork, user->Twork)); in DesignMatMultTranspose()
484 PetscCall(MatMult(user->AvT, user->Twork, user->yiwork[i])); in DesignMatMultTranspose()
487 PetscCall(VecPointwiseMult(user->yiwork[i], user->uwork, user->yiwork[i])); in DesignMatMultTranspose()
488 PetscCall(VecAXPY(Y, user->ht, user->yiwork[i])); in DesignMatMultTranspose()
495 AppCtx *user; in StateMatBlockPrecMult() local
498 PetscCall(PCShellGetContext(PC_shell, &user)); in StateMatBlockPrecMult()
500 PetscCheck(user->dsg_formed, PETSC_COMM_SELF, PETSC_ERR_SUP, "DSG not formed"); in StateMatBlockPrecMult()
501 …PetscCall(MatSOR(user->DSG, X, 1.0, (MatSORType)(SOR_ZERO_INITIAL_GUESS | SOR_LOCAL_SYMMETRIC_SWEE… in StateMatBlockPrecMult()
507 AppCtx *user; in StateMatInvMult() local
511 PetscCall(MatShellGetContext(J_shell, &user)); in StateMatInvMult()
513 if (Y == user->ytrue) { in StateMatInvMult()
515 PetscCall(KSPSetTolerances(user->solver, 1e-8, PETSC_CURRENT, PETSC_CURRENT, PETSC_CURRENT)); in StateMatInvMult()
517 …PetscCall(KSPSetTolerances(user->solver, PETSC_CURRENT, PETSC_CURRENT, PETSC_CURRENT, PETSC_CURREN… in StateMatInvMult()
520 PetscCall(Scatter_i(X, user->yi, user->yi_scatter, user->nt)); in StateMatInvMult()
521 PetscCall(KSPSolve(user->solver, user->yi[0], user->yiwork[0])); in StateMatInvMult()
522 PetscCall(KSPGetIterationNumber(user->solver, &its)); in StateMatInvMult()
523 user->ksp_its = user->ksp_its + its; in StateMatInvMult()
525 for (i = 1; i < user->nt; i++) { in StateMatInvMult()
526 PetscCall(VecAXPY(user->yi[i], 1.0, user->yiwork[i - 1])); in StateMatInvMult()
527 PetscCall(KSPSolve(user->solver, user->yi[i], user->yiwork[i])); in StateMatInvMult()
528 PetscCall(KSPGetIterationNumber(user->solver, &its)); in StateMatInvMult()
529 user->ksp_its = user->ksp_its + its; in StateMatInvMult()
531 PetscCall(Gather_i(Y, user->yiwork, user->yi_scatter, user->nt)); in StateMatInvMult()
537 AppCtx *user; in StateMatInvTransposeMult() local
541 PetscCall(MatShellGetContext(J_shell, &user)); in StateMatInvTransposeMult()
543 PetscCall(Scatter_i(X, user->yi, user->yi_scatter, user->nt)); in StateMatInvTransposeMult()
545 i = user->nt - 1; in StateMatInvTransposeMult()
546 PetscCall(KSPSolve(user->solver, user->yi[i], user->yiwork[i])); in StateMatInvTransposeMult()
548 PetscCall(KSPGetIterationNumber(user->solver, &its)); in StateMatInvTransposeMult()
549 user->ksp_its = user->ksp_its + its; in StateMatInvTransposeMult()
551 for (i = user->nt - 2; i >= 0; i--) { in StateMatInvTransposeMult()
552 PetscCall(VecAXPY(user->yi[i], 1.0, user->yiwork[i + 1])); in StateMatInvTransposeMult()
553 PetscCall(KSPSolve(user->solver, user->yi[i], user->yiwork[i])); in StateMatInvTransposeMult()
555 PetscCall(KSPGetIterationNumber(user->solver, &its)); in StateMatInvTransposeMult()
556 user->ksp_its = user->ksp_its + its; in StateMatInvTransposeMult()
559 PetscCall(Gather_i(Y, user->yiwork, user->yi_scatter, user->nt)); in StateMatInvTransposeMult()
565 AppCtx *user; in StateMatDuplicate() local
568 PetscCall(MatShellGetContext(J_shell, &user)); in StateMatDuplicate()
570 …tCreateShell(PETSC_COMM_WORLD, PETSC_DETERMINE, PETSC_DETERMINE, user->m, user->m, user, new_shell… in StateMatDuplicate()
580 AppCtx *user; in StateMatGetDiagonal() local
583 PetscCall(MatShellGetContext(J_shell, &user)); in StateMatGetDiagonal()
584 PetscCall(VecCopy(user->js_diag, X)); in StateMatGetDiagonal()
597 AppCtx *user = (AppCtx *)ptr; in FormConstraints() local
600 PetscCall(Scatter(X, user->y, user->state_scatter, user->u, user->design_scatter)); in FormConstraints()
601 PetscCall(Scatter_i(user->y, user->yi, user->yi_scatter, user->nt)); in FormConstraints()
602 PetscCall(MatMult(user->JsBlock, user->yi[0], user->yiwork[0])); in FormConstraints()
603 for (i = 1; i < user->nt; i++) { in FormConstraints()
604 PetscCall(MatMult(user->JsBlock, user->yi[i], user->yiwork[i])); in FormConstraints()
605 PetscCall(VecAXPY(user->yiwork[i], -1.0, user->yi[i - 1])); in FormConstraints()
607 PetscCall(Gather_i(C, user->yiwork, user->yi_scatter, user->nt)); in FormConstraints()
608 PetscCall(VecAXPY(C, -1.0, user->q)); in FormConstraints()
656 PetscErrorCode ParabolicInitialize(AppCtx *user) in ParabolicInitialize() argument
681 PetscCall(PetscMalloc1(user->mx, &x)); in ParabolicInitialize()
682 PetscCall(PetscMalloc1(user->mx, &y)); in ParabolicInitialize()
683 PetscCall(PetscMalloc1(user->mx, &z)); in ParabolicInitialize()
684 user->jformed = PETSC_FALSE; in ParabolicInitialize()
685 user->dsg_formed = PETSC_FALSE; in ParabolicInitialize()
687 n = user->mx * user->mx * user->mx; in ParabolicInitialize()
688 m = 3 * user->mx * user->mx * (user->mx - 1); in ParabolicInitialize()
689 sqrt_beta = PetscSqrtScalar(user->beta); in ParabolicInitialize()
691 user->ksp_its = 0; in ParabolicInitialize()
692 user->ksp_its_initial = 0; in ParabolicInitialize()
694 stime = (PetscReal)user->nt / user->ns; in ParabolicInitialize()
695 PetscCall(PetscMalloc1(user->ns, &user->sample_times)); in ParabolicInitialize()
696 …for (i = 0; i < user->ns; i++) user->sample_times[i] = (PetscInt)(stime * ((PetscReal)i + 1.0) - 0… in ParabolicInitialize()
699 PetscCall(VecCreate(PETSC_COMM_WORLD, &user->q)); in ParabolicInitialize()
701 PetscCall(VecSetSizes(user->q, PETSC_DECIDE, n * user->nt)); in ParabolicInitialize()
703 PetscCall(VecSetFromOptions(user->q)); in ParabolicInitialize()
711 PetscCall(VecDuplicate(XX, &user->utrue)); in ParabolicInitialize()
715 h = 1.0 / user->mx; in ParabolicInitialize()
716 hinv = user->mx; in ParabolicInitialize()
721 i = linear_index % user->mx; in ParabolicInitialize()
722 j = ((linear_index - i) / user->mx) % user->mx; in ParabolicInitialize()
723 k = ((linear_index - i) / user->mx - j) / user->mx; in ParabolicInitialize()
759 PetscCall(VecCopy(XXwork, user->utrue)); in ParabolicInitialize()
760 PetscCall(VecAXPY(user->utrue, 1.0, YYwork)); in ParabolicInitialize()
761 PetscCall(VecAXPY(user->utrue, 1.0, ZZwork)); in ParabolicInitialize()
762 PetscCall(VecScale(user->utrue, -10.0)); in ParabolicInitialize()
763 PetscCall(VecExp(user->utrue)); in ParabolicInitialize()
764 PetscCall(VecShift(user->utrue, 0.5)); in ParabolicInitialize()
775 PetscCall(VecDuplicate(user->utrue, &user->ur)); in ParabolicInitialize()
776 PetscCall(VecSet(user->ur, 0.0)); in ParabolicInitialize()
779 PetscCall(MatCreate(PETSC_COMM_WORLD, &user->Grad)); in ParabolicInitialize()
780 PetscCall(MatSetSizes(user->Grad, PETSC_DECIDE, PETSC_DECIDE, m, n)); in ParabolicInitialize()
781 PetscCall(MatSetFromOptions(user->Grad)); in ParabolicInitialize()
782 PetscCall(MatMPIAIJSetPreallocation(user->Grad, 2, NULL, 2, NULL)); in ParabolicInitialize()
783 PetscCall(MatSeqAIJSetPreallocation(user->Grad, 2, NULL)); in ParabolicInitialize()
784 PetscCall(MatGetOwnershipRange(user->Grad, &istart, &iend)); in ParabolicInitialize()
788 iblock = i / (user->mx - 1); in ParabolicInitialize()
789 j = iblock * user->mx + (i % (user->mx - 1)); in ParabolicInitialize()
790 PetscCall(MatSetValues(user->Grad, 1, &i, 1, &j, &neg_hinv, INSERT_VALUES)); in ParabolicInitialize()
792 PetscCall(MatSetValues(user->Grad, 1, &i, 1, &j, &hinv, INSERT_VALUES)); in ParabolicInitialize()
795 iblock = (i - m / 3) / (user->mx * (user->mx - 1)); in ParabolicInitialize()
796 j = iblock * user->mx * user->mx + ((i - m / 3) % (user->mx * (user->mx - 1))); in ParabolicInitialize()
797 PetscCall(MatSetValues(user->Grad, 1, &i, 1, &j, &neg_hinv, INSERT_VALUES)); in ParabolicInitialize()
798 j = j + user->mx; in ParabolicInitialize()
799 PetscCall(MatSetValues(user->Grad, 1, &i, 1, &j, &hinv, INSERT_VALUES)); in ParabolicInitialize()
803 PetscCall(MatSetValues(user->Grad, 1, &i, 1, &j, &neg_hinv, INSERT_VALUES)); in ParabolicInitialize()
804 j = j + user->mx * user->mx; in ParabolicInitialize()
805 PetscCall(MatSetValues(user->Grad, 1, &i, 1, &j, &hinv, INSERT_VALUES)); in ParabolicInitialize()
809 PetscCall(MatAssemblyBegin(user->Grad, MAT_FINAL_ASSEMBLY)); in ParabolicInitialize()
810 PetscCall(MatAssemblyEnd(user->Grad, MAT_FINAL_ASSEMBLY)); in ParabolicInitialize()
813 PetscCall(MatCreate(PETSC_COMM_WORLD, &user->Av)); in ParabolicInitialize()
814 PetscCall(MatSetSizes(user->Av, PETSC_DECIDE, PETSC_DECIDE, m, n)); in ParabolicInitialize()
815 PetscCall(MatSetFromOptions(user->Av)); in ParabolicInitialize()
816 PetscCall(MatMPIAIJSetPreallocation(user->Av, 2, NULL, 2, NULL)); in ParabolicInitialize()
817 PetscCall(MatSeqAIJSetPreallocation(user->Av, 2, NULL)); in ParabolicInitialize()
818 PetscCall(MatGetOwnershipRange(user->Av, &istart, &iend)); in ParabolicInitialize()
822 iblock = i / (user->mx - 1); in ParabolicInitialize()
823 j = iblock * user->mx + (i % (user->mx - 1)); in ParabolicInitialize()
824 PetscCall(MatSetValues(user->Av, 1, &i, 1, &j, &half, INSERT_VALUES)); in ParabolicInitialize()
826 PetscCall(MatSetValues(user->Av, 1, &i, 1, &j, &half, INSERT_VALUES)); in ParabolicInitialize()
829 iblock = (i - m / 3) / (user->mx * (user->mx - 1)); in ParabolicInitialize()
830 j = iblock * user->mx * user->mx + ((i - m / 3) % (user->mx * (user->mx - 1))); in ParabolicInitialize()
831 PetscCall(MatSetValues(user->Av, 1, &i, 1, &j, &half, INSERT_VALUES)); in ParabolicInitialize()
832 j = j + user->mx; in ParabolicInitialize()
833 PetscCall(MatSetValues(user->Av, 1, &i, 1, &j, &half, INSERT_VALUES)); in ParabolicInitialize()
837 PetscCall(MatSetValues(user->Av, 1, &i, 1, &j, &half, INSERT_VALUES)); in ParabolicInitialize()
838 j = j + user->mx * user->mx; in ParabolicInitialize()
839 PetscCall(MatSetValues(user->Av, 1, &i, 1, &j, &half, INSERT_VALUES)); in ParabolicInitialize()
843 PetscCall(MatAssemblyBegin(user->Av, MAT_FINAL_ASSEMBLY)); in ParabolicInitialize()
844 PetscCall(MatAssemblyEnd(user->Av, MAT_FINAL_ASSEMBLY)); in ParabolicInitialize()
847 PetscCall(MatTranspose(user->Av, MAT_INITIAL_MATRIX, &user->AvT)); in ParabolicInitialize()
849 PetscCall(MatCreate(PETSC_COMM_WORLD, &user->L)); in ParabolicInitialize()
850 PetscCall(MatSetSizes(user->L, PETSC_DECIDE, PETSC_DECIDE, m + n, n)); in ParabolicInitialize()
851 PetscCall(MatSetFromOptions(user->L)); in ParabolicInitialize()
852 PetscCall(MatMPIAIJSetPreallocation(user->L, 2, NULL, 2, NULL)); in ParabolicInitialize()
853 PetscCall(MatSeqAIJSetPreallocation(user->L, 2, NULL)); in ParabolicInitialize()
854 PetscCall(MatGetOwnershipRange(user->L, &istart, &iend)); in ParabolicInitialize()
858 iblock = i / (user->mx - 1); in ParabolicInitialize()
859 j = iblock * user->mx + (i % (user->mx - 1)); in ParabolicInitialize()
860 PetscCall(MatSetValues(user->L, 1, &i, 1, &j, &neg_hinv, INSERT_VALUES)); in ParabolicInitialize()
862 PetscCall(MatSetValues(user->L, 1, &i, 1, &j, &hinv, INSERT_VALUES)); in ParabolicInitialize()
865 iblock = (i - m / 3) / (user->mx * (user->mx - 1)); in ParabolicInitialize()
866 j = iblock * user->mx * user->mx + ((i - m / 3) % (user->mx * (user->mx - 1))); in ParabolicInitialize()
867 PetscCall(MatSetValues(user->L, 1, &i, 1, &j, &neg_hinv, INSERT_VALUES)); in ParabolicInitialize()
868 j = j + user->mx; in ParabolicInitialize()
869 PetscCall(MatSetValues(user->L, 1, &i, 1, &j, &hinv, INSERT_VALUES)); in ParabolicInitialize()
873 PetscCall(MatSetValues(user->L, 1, &i, 1, &j, &neg_hinv, INSERT_VALUES)); in ParabolicInitialize()
874 j = j + user->mx * user->mx; in ParabolicInitialize()
875 PetscCall(MatSetValues(user->L, 1, &i, 1, &j, &hinv, INSERT_VALUES)); in ParabolicInitialize()
879 PetscCall(MatSetValues(user->L, 1, &i, 1, &j, &sqrt_beta, INSERT_VALUES)); in ParabolicInitialize()
883 PetscCall(MatAssemblyBegin(user->L, MAT_FINAL_ASSEMBLY)); in ParabolicInitialize()
884 PetscCall(MatAssemblyEnd(user->L, MAT_FINAL_ASSEMBLY)); in ParabolicInitialize()
886 PetscCall(MatScale(user->L, PetscPowScalar(h, 1.5))); in ParabolicInitialize()
889 PetscCall(MatTranspose(user->Grad, MAT_INITIAL_MATRIX, &user->Div)); in ParabolicInitialize()
892 PetscCall(VecCreate(PETSC_COMM_WORLD, &user->S)); in ParabolicInitialize()
893 PetscCall(VecSetSizes(user->S, PETSC_DECIDE, user->mx * user->mx * (user->mx - 1) * 3)); in ParabolicInitialize()
894 PetscCall(VecSetFromOptions(user->S)); in ParabolicInitialize()
896 PetscCall(VecCreate(PETSC_COMM_WORLD, &user->lwork)); in ParabolicInitialize()
897 PetscCall(VecSetSizes(user->lwork, PETSC_DECIDE, m + user->mx * user->mx * user->mx)); in ParabolicInitialize()
898 PetscCall(VecSetFromOptions(user->lwork)); in ParabolicInitialize()
900 PetscCall(MatDuplicate(user->Div, MAT_SHARE_NONZERO_PATTERN, &user->Divwork)); in ParabolicInitialize()
901 PetscCall(MatDuplicate(user->Av, MAT_SHARE_NONZERO_PATTERN, &user->Avwork)); in ParabolicInitialize()
903 PetscCall(VecCreate(PETSC_COMM_WORLD, &user->d)); in ParabolicInitialize()
904 PetscCall(VecSetSizes(user->d, PETSC_DECIDE, user->ndata * user->nt)); in ParabolicInitialize()
905 PetscCall(VecSetFromOptions(user->d)); in ParabolicInitialize()
907 PetscCall(VecDuplicate(user->S, &user->Swork)); in ParabolicInitialize()
908 PetscCall(VecDuplicate(user->S, &user->Av_u)); in ParabolicInitialize()
909 PetscCall(VecDuplicate(user->S, &user->Twork)); in ParabolicInitialize()
910 PetscCall(VecDuplicate(user->S, &user->Rwork)); in ParabolicInitialize()
911 PetscCall(VecDuplicate(user->y, &user->ywork)); in ParabolicInitialize()
912 PetscCall(VecDuplicate(user->u, &user->uwork)); in ParabolicInitialize()
913 PetscCall(VecDuplicate(user->u, &user->js_diag)); in ParabolicInitialize()
914 PetscCall(VecDuplicate(user->c, &user->cwork)); in ParabolicInitialize()
915 PetscCall(VecDuplicate(user->d, &user->dwork)); in ParabolicInitialize()
918 …ETSC_COMM_WORLD, PETSC_DETERMINE, PETSC_DETERMINE, user->m * user->nt, user->m * user->nt, user, &… in ParabolicInitialize()
919 PetscCall(MatShellSetOperation(user->Js, MATOP_MULT, (PetscErrorCodeFn *)StateMatMult)); in ParabolicInitialize()
920 PetscCall(MatShellSetOperation(user->Js, MATOP_DUPLICATE, (PetscErrorCodeFn *)StateMatDuplicate)); in ParabolicInitialize()
921 …PetscCall(MatShellSetOperation(user->Js, MATOP_MULT_TRANSPOSE, (PetscErrorCodeFn *)StateMatMultTra… in ParabolicInitialize()
922 …PetscCall(MatShellSetOperation(user->Js, MATOP_GET_DIAGONAL, (PetscErrorCodeFn *)StateMatGetDiagon… in ParabolicInitialize()
925 …eateShell(PETSC_COMM_WORLD, PETSC_DETERMINE, PETSC_DETERMINE, user->m, user->m, user, &user->JsBlo… in ParabolicInitialize()
926 PetscCall(MatShellSetOperation(user->JsBlock, MATOP_MULT, (PetscErrorCodeFn *)StateMatBlockMult)); in ParabolicInitialize()
928 …PetscCall(MatShellSetOperation(user->JsBlock, MATOP_MULT_TRANSPOSE, (PetscErrorCodeFn *)StateMatBl… in ParabolicInitialize()
933 …eateShell(PETSC_COMM_WORLD, PETSC_DETERMINE, PETSC_DETERMINE, user->m, user->m, user, &user->JsBlo… in ParabolicInitialize()
934 …PetscCall(MatShellSetOperation(user->JsBlockPrec, MATOP_MULT, (PetscErrorCodeFn *)StateMatBlockPre… in ParabolicInitialize()
936 …PetscCall(MatShellSetOperation(user->JsBlockPrec, MATOP_MULT_TRANSPOSE, (PetscErrorCodeFn *)StateM… in ParabolicInitialize()
937 PetscCall(MatSetOption(user->JsBlockPrec, MAT_SYMMETRIC, PETSC_TRUE)); in ParabolicInitialize()
938 PetscCall(MatSetOption(user->JsBlockPrec, MAT_SYMMETRY_ETERNAL, PETSC_TRUE)); in ParabolicInitialize()
941 …ell(PETSC_COMM_WORLD, PETSC_DETERMINE, PETSC_DETERMINE, user->m * user->nt, user->m, user, &user->… in ParabolicInitialize()
942 PetscCall(MatShellSetOperation(user->Jd, MATOP_MULT, (PetscErrorCodeFn *)DesignMatMult)); in ParabolicInitialize()
943 …PetscCall(MatShellSetOperation(user->Jd, MATOP_MULT_TRANSPOSE, (PetscErrorCodeFn *)DesignMatMultTr… in ParabolicInitialize()
946 …ETSC_COMM_WORLD, PETSC_DETERMINE, PETSC_DETERMINE, user->m * user->nt, user->m * user->nt, user, &… in ParabolicInitialize()
947 PetscCall(MatShellSetOperation(user->JsInv, MATOP_MULT, (PetscErrorCodeFn *)StateMatInvMult)); in ParabolicInitialize()
948 …PetscCall(MatShellSetOperation(user->JsInv, MATOP_MULT_TRANSPOSE, (PetscErrorCodeFn *)StateMatInvT… in ParabolicInitialize()
951 PetscCall(KSPCreate(PETSC_COMM_WORLD, &user->solver)); in ParabolicInitialize()
952 PetscCall(KSPSetType(user->solver, KSPCG)); in ParabolicInitialize()
953 PetscCall(KSPSetOperators(user->solver, user->JsBlock, user->JsBlockPrec)); in ParabolicInitialize()
954 PetscCall(KSPSetInitialGuessNonzero(user->solver, PETSC_FALSE)); in ParabolicInitialize()
955 PetscCall(KSPSetTolerances(user->solver, 1e-4, 1e-20, 1e3, 500)); in ParabolicInitialize()
956 PetscCall(KSPSetFromOptions(user->solver)); in ParabolicInitialize()
957 PetscCall(KSPGetPC(user->solver, &user->prec)); in ParabolicInitialize()
958 PetscCall(PCSetType(user->prec, PCSHELL)); in ParabolicInitialize()
960 PetscCall(PCShellSetApply(user->prec, StateMatBlockPrecMult)); in ParabolicInitialize()
961 PetscCall(PCShellSetApplyTranspose(user->prec, StateMatBlockPrecMult)); in ParabolicInitialize()
962 PetscCall(PCShellSetContext(user->prec, user)); in ParabolicInitialize()
965 PetscCall(PetscMalloc1(user->nt * user->m, &user->yi_scatter)); in ParabolicInitialize()
967 PetscCall(VecSetSizes(yi, PETSC_DECIDE, user->mx * user->mx * user->mx)); in ParabolicInitialize()
969 PetscCall(VecDuplicateVecs(yi, user->nt, &user->yi)); in ParabolicInitialize()
970 PetscCall(VecDuplicateVecs(yi, user->nt, &user->yiwork)); in ParabolicInitialize()
972 PetscCall(VecGetOwnershipRange(user->y, &lo2, &hi2)); in ParabolicInitialize()
974 for (i = 0; i < user->nt; i++) { in ParabolicInitialize()
975 PetscCall(VecGetOwnershipRange(user->yi[i], &lo, &hi)); in ParabolicInitialize()
978 PetscCall(VecScatterCreate(user->y, is_from_y, user->yi[i], is_to_yi, &user->yi_scatter[i])); in ParabolicInitialize()
986 PetscCall(PetscMalloc1(user->ns * user->ndata, &user->di_scatter)); in ParabolicInitialize()
988 PetscCall(VecSetSizes(di, PETSC_DECIDE, user->ndata)); in ParabolicInitialize()
990 PetscCall(VecDuplicateVecs(di, user->ns, &user->di)); in ParabolicInitialize()
991 PetscCall(VecGetOwnershipRange(user->d, &lo2, &hi2)); in ParabolicInitialize()
993 for (i = 0; i < user->ns; i++) { in ParabolicInitialize()
994 PetscCall(VecGetOwnershipRange(user->di[i], &lo, &hi)); in ParabolicInitialize()
997 PetscCall(VecScatterCreate(user->d, is_from_d, user->di[i], is_to_di, &user->di_scatter[i])); in ParabolicInitialize()
1005 PetscCall(VecCopy(bc, user->yiwork[0])); in ParabolicInitialize()
1006 for (i = 1; i < user->nt; i++) PetscCall(VecSet(user->yiwork[i], 0.0)); in ParabolicInitialize()
1007 PetscCall(Gather_i(user->q, user->yiwork, user->yi_scatter, user->nt)); in ParabolicInitialize()
1011 PetscCall(VecCreate(PETSC_COMM_WORLD, &user->ytrue)); in ParabolicInitialize()
1012 PetscCall(VecSetSizes(user->ytrue, PETSC_DECIDE, n * user->nt)); in ParabolicInitialize()
1013 PetscCall(VecSetFromOptions(user->ytrue)); in ParabolicInitialize()
1016 PetscCall(VecSet(user->uwork, 0)); in ParabolicInitialize()
1017 PetscCall(VecAXPY(user->uwork, -1.0, user->utrue)); /* Note: user->utrue */ in ParabolicInitialize()
1018 PetscCall(VecExp(user->uwork)); in ParabolicInitialize()
1019 PetscCall(MatMult(user->Av, user->uwork, user->Av_u)); in ParabolicInitialize()
1022 PetscCall(MatProductCreate(user->Div, user->Grad, NULL, &user->DSG)); in ParabolicInitialize()
1023 PetscCall(MatProductSetType(user->DSG, MATPRODUCT_AB)); in ParabolicInitialize()
1024 PetscCall(MatProductSetAlgorithm(user->DSG, "default")); in ParabolicInitialize()
1025 PetscCall(MatProductSetFill(user->DSG, PETSC_DETERMINE)); in ParabolicInitialize()
1026 PetscCall(MatProductSetFromOptions(user->DSG)); in ParabolicInitialize()
1027 PetscCall(MatProductSymbolic(user->DSG)); in ParabolicInitialize()
1029 user->dsg_formed = PETSC_TRUE; in ParabolicInitialize()
1032 PetscCall(MatCopy(user->Div, user->Divwork, SAME_NONZERO_PATTERN)); in ParabolicInitialize()
1033 PetscCall(MatDiagonalScale(user->Divwork, NULL, user->Av_u)); in ParabolicInitialize()
1034 if (user->dsg_formed) { in ParabolicInitialize()
1035 PetscCall(MatProductNumeric(user->DSG)); in ParabolicInitialize()
1037 PetscCall(MatMatMult(user->Div, user->Grad, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &user->DSG)); in ParabolicInitialize()
1038 user->dsg_formed = PETSC_TRUE; in ParabolicInitialize()
1041 PetscCall(MatScale(user->DSG, user->ht)); in ParabolicInitialize()
1042 PetscCall(MatShift(user->DSG, 1.0)); in ParabolicInitialize()
1045 PetscCall(StateMatInvMult(user->Js, user->q, user->ytrue)); in ParabolicInitialize()
1050 PetscCall(VecSet(user->uwork, 0)); in ParabolicInitialize()
1051 PetscCall(VecAXPY(user->uwork, -1.0, user->u)); /* Note: user->u */ in ParabolicInitialize()
1052 PetscCall(VecExp(user->uwork)); in ParabolicInitialize()
1053 PetscCall(MatMult(user->Av, user->uwork, user->Av_u)); in ParabolicInitialize()
1056 PetscCall(MatCopy(user->Div, user->Divwork, SAME_NONZERO_PATTERN)); in ParabolicInitialize()
1057 PetscCall(MatDiagonalScale(user->Divwork, NULL, user->Av_u)); in ParabolicInitialize()
1058 if (user->dsg_formed) { in ParabolicInitialize()
1059 PetscCall(MatProductNumeric(user->DSG)); in ParabolicInitialize()
1061 PetscCall(MatMatMult(user->Div, user->Grad, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &user->DSG)); in ParabolicInitialize()
1063 user->dsg_formed = PETSC_TRUE; in ParabolicInitialize()
1066 PetscCall(MatScale(user->DSG, user->ht)); in ParabolicInitialize()
1067 PetscCall(MatShift(user->DSG, 1.0)); in ParabolicInitialize()
1070 PetscCall(StateMatInvMult(user->Js, user->q, user->y)); in ParabolicInitialize()
1073 PetscCall(MatCreate(PETSC_COMM_WORLD, &user->Qblock)); in ParabolicInitialize()
1074 PetscCall(MatSetSizes(user->Qblock, PETSC_DECIDE, PETSC_DECIDE, user->ndata, n)); in ParabolicInitialize()
1075 PetscCall(MatSetFromOptions(user->Qblock)); in ParabolicInitialize()
1076 PetscCall(MatMPIAIJSetPreallocation(user->Qblock, 8, NULL, 8, NULL)); in ParabolicInitialize()
1077 PetscCall(MatSeqAIJSetPreallocation(user->Qblock, 8, NULL)); in ParabolicInitialize()
1079 for (i = 0; i < user->mx; i++) { in ParabolicInitialize()
1085 PetscCall(MatGetOwnershipRange(user->Qblock, &istart, &iend)); in ParabolicInitialize()
1086 nx = user->mx; in ParabolicInitialize()
1087 ny = user->mx; in ParabolicInitialize()
1088 nz = user->mx; in ParabolicInitialize()
1132 PetscCall(MatSetValues(user->Qblock, 1, &i, 1, &j, &v, INSERT_VALUES)); in ParabolicInitialize()
1136 PetscCall(MatSetValues(user->Qblock, 1, &i, 1, &j, &v, INSERT_VALUES)); in ParabolicInitialize()
1140 PetscCall(MatSetValues(user->Qblock, 1, &i, 1, &j, &v, INSERT_VALUES)); in ParabolicInitialize()
1144 PetscCall(MatSetValues(user->Qblock, 1, &i, 1, &j, &v, INSERT_VALUES)); in ParabolicInitialize()
1148 PetscCall(MatSetValues(user->Qblock, 1, &i, 1, &j, &v, INSERT_VALUES)); in ParabolicInitialize()
1152 PetscCall(MatSetValues(user->Qblock, 1, &i, 1, &j, &v, INSERT_VALUES)); in ParabolicInitialize()
1156 PetscCall(MatSetValues(user->Qblock, 1, &i, 1, &j, &v, INSERT_VALUES)); in ParabolicInitialize()
1160 PetscCall(MatSetValues(user->Qblock, 1, &i, 1, &j, &v, INSERT_VALUES)); in ParabolicInitialize()
1162 PetscCall(MatAssemblyBegin(user->Qblock, MAT_FINAL_ASSEMBLY)); in ParabolicInitialize()
1163 PetscCall(MatAssemblyEnd(user->Qblock, MAT_FINAL_ASSEMBLY)); in ParabolicInitialize()
1165 PetscCall(MatTranspose(user->Qblock, MAT_INITIAL_MATRIX, &user->QblockT)); in ParabolicInitialize()
1166 PetscCall(MatTranspose(user->L, MAT_INITIAL_MATRIX, &user->LT)); in ParabolicInitialize()
1169 PetscCall(VecSet(user->ywork, 1.0)); in ParabolicInitialize()
1170 PetscCall(VecAYPX(user->ywork, user->noise, user->ytrue)); in ParabolicInitialize()
1171 PetscCall(Scatter_i(user->ywork, user->yiwork, user->yi_scatter, user->nt)); in ParabolicInitialize()
1172 for (j = 0; j < user->ns; j++) { in ParabolicInitialize()
1173 i = user->sample_times[j]; in ParabolicInitialize()
1174 PetscCall(MatMult(user->Qblock, user->yiwork[i], user->di[j])); in ParabolicInitialize()
1176 PetscCall(Gather_i(user->d, user->di, user->di_scatter, user->ns)); in ParabolicInitialize()
1179 PetscCall(KSPSetFromOptions(user->solver)); in ParabolicInitialize()
1186 PetscErrorCode ParabolicDestroy(AppCtx *user) in ParabolicDestroy() argument
1191 PetscCall(MatDestroy(&user->Qblock)); in ParabolicDestroy()
1192 PetscCall(MatDestroy(&user->QblockT)); in ParabolicDestroy()
1193 PetscCall(MatDestroy(&user->Div)); in ParabolicDestroy()
1194 PetscCall(MatDestroy(&user->Divwork)); in ParabolicDestroy()
1195 PetscCall(MatDestroy(&user->Grad)); in ParabolicDestroy()
1196 PetscCall(MatDestroy(&user->Av)); in ParabolicDestroy()
1197 PetscCall(MatDestroy(&user->Avwork)); in ParabolicDestroy()
1198 PetscCall(MatDestroy(&user->AvT)); in ParabolicDestroy()
1199 PetscCall(MatDestroy(&user->DSG)); in ParabolicDestroy()
1200 PetscCall(MatDestroy(&user->L)); in ParabolicDestroy()
1201 PetscCall(MatDestroy(&user->LT)); in ParabolicDestroy()
1202 PetscCall(KSPDestroy(&user->solver)); in ParabolicDestroy()
1203 PetscCall(MatDestroy(&user->Js)); in ParabolicDestroy()
1204 PetscCall(MatDestroy(&user->Jd)); in ParabolicDestroy()
1205 PetscCall(MatDestroy(&user->JsInv)); in ParabolicDestroy()
1206 PetscCall(MatDestroy(&user->JsBlock)); in ParabolicDestroy()
1207 PetscCall(MatDestroy(&user->JsBlockPrec)); in ParabolicDestroy()
1208 PetscCall(VecDestroy(&user->u)); in ParabolicDestroy()
1209 PetscCall(VecDestroy(&user->uwork)); in ParabolicDestroy()
1210 PetscCall(VecDestroy(&user->utrue)); in ParabolicDestroy()
1211 PetscCall(VecDestroy(&user->y)); in ParabolicDestroy()
1212 PetscCall(VecDestroy(&user->ywork)); in ParabolicDestroy()
1213 PetscCall(VecDestroy(&user->ytrue)); in ParabolicDestroy()
1214 PetscCall(VecDestroyVecs(user->nt, &user->yi)); in ParabolicDestroy()
1215 PetscCall(VecDestroyVecs(user->nt, &user->yiwork)); in ParabolicDestroy()
1216 PetscCall(VecDestroyVecs(user->ns, &user->di)); in ParabolicDestroy()
1217 PetscCall(PetscFree(user->yi)); in ParabolicDestroy()
1218 PetscCall(PetscFree(user->yiwork)); in ParabolicDestroy()
1219 PetscCall(PetscFree(user->di)); in ParabolicDestroy()
1220 PetscCall(VecDestroy(&user->c)); in ParabolicDestroy()
1221 PetscCall(VecDestroy(&user->cwork)); in ParabolicDestroy()
1222 PetscCall(VecDestroy(&user->ur)); in ParabolicDestroy()
1223 PetscCall(VecDestroy(&user->q)); in ParabolicDestroy()
1224 PetscCall(VecDestroy(&user->d)); in ParabolicDestroy()
1225 PetscCall(VecDestroy(&user->dwork)); in ParabolicDestroy()
1226 PetscCall(VecDestroy(&user->lwork)); in ParabolicDestroy()
1227 PetscCall(VecDestroy(&user->S)); in ParabolicDestroy()
1228 PetscCall(VecDestroy(&user->Swork)); in ParabolicDestroy()
1229 PetscCall(VecDestroy(&user->Av_u)); in ParabolicDestroy()
1230 PetscCall(VecDestroy(&user->Twork)); in ParabolicDestroy()
1231 PetscCall(VecDestroy(&user->Rwork)); in ParabolicDestroy()
1232 PetscCall(VecDestroy(&user->js_diag)); in ParabolicDestroy()
1233 PetscCall(ISDestroy(&user->s_is)); in ParabolicDestroy()
1234 PetscCall(ISDestroy(&user->d_is)); in ParabolicDestroy()
1235 PetscCall(VecScatterDestroy(&user->state_scatter)); in ParabolicDestroy()
1236 PetscCall(VecScatterDestroy(&user->design_scatter)); in ParabolicDestroy()
1237 for (i = 0; i < user->nt; i++) PetscCall(VecScatterDestroy(&user->yi_scatter[i])); in ParabolicDestroy()
1238 for (i = 0; i < user->ns; i++) PetscCall(VecScatterDestroy(&user->di_scatter[i])); in ParabolicDestroy()
1239 PetscCall(PetscFree(user->yi_scatter)); in ParabolicDestroy()
1240 PetscCall(PetscFree(user->di_scatter)); in ParabolicDestroy()
1241 PetscCall(PetscFree(user->sample_times)); in ParabolicDestroy()
1249 AppCtx *user = (AppCtx *)ptr; in ParabolicMonitor() local
1253 PetscCall(Scatter(X, user->ywork, user->state_scatter, user->uwork, user->design_scatter)); in ParabolicMonitor()
1254 PetscCall(VecAXPY(user->ywork, -1.0, user->ytrue)); in ParabolicMonitor()
1255 PetscCall(VecAXPY(user->uwork, -1.0, user->utrue)); in ParabolicMonitor()
1256 PetscCall(VecNorm(user->uwork, NORM_2, &unorm)); in ParabolicMonitor()
1257 PetscCall(VecNorm(user->ywork, NORM_2, &ynorm)); in ParabolicMonitor()