1 #include <petsc/private/taoimpl.h> 2 3 /*T 4 Concepts: TAO^Solving a system of nonlinear equations, nonlinear least squares 5 Routines: TaoCreate(); 6 Routines: TaoSetType(); 7 Routines: TaoSetSolution(); 8 Routines: TaoSetObjective(); 9 Routines: TaoSetGradient(); 10 Routines: TaoSetConstraintsRoutine(); 11 Routines: TaoSetJacobianStateRoutine(); 12 Routines: TaoSetJacobianDesignRoutine(); 13 Routines: TaoSetStateDesignIS(); 14 Routines: TaoSetFromOptions(); 15 Routines: TaoSolve(); 16 Routines: TaoDestroy(); 17 Processors: n 18 T*/ 19 20 typedef struct { 21 PetscInt n; /* Number of variables */ 22 PetscInt m; /* Number of constraints per time step */ 23 PetscInt mx; /* grid points in each direction */ 24 PetscInt nt; /* Number of time steps; as of now, must be divisible by 8 */ 25 PetscInt ndata; /* Number of data points per sample */ 26 PetscInt ns; /* Number of samples */ 27 PetscInt *sample_times; /* Times of samples */ 28 IS s_is; 29 IS d_is; 30 31 VecScatter state_scatter; 32 VecScatter design_scatter; 33 VecScatter *yi_scatter; 34 VecScatter *di_scatter; 35 36 Mat Js,Jd,JsBlockPrec,JsInv,JsBlock; 37 PetscBool jformed,dsg_formed; 38 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 ht; /* Time step */ 43 44 Mat Qblock,QblockT; 45 Mat L,LT; 46 Mat Div,Divwork; 47 Mat Grad; 48 Mat Av,Avwork,AvT; 49 Mat DSG; 50 Vec q; 51 Vec ur; /* reference */ 52 53 Vec d; 54 Vec dwork; 55 Vec *di; 56 57 Vec y; /* state variables */ 58 Vec ywork; 59 60 Vec ytrue; 61 Vec *yi,*yiwork; 62 63 Vec u; /* design variables */ 64 Vec uwork; 65 66 Vec utrue; 67 Vec js_diag; 68 Vec c; /* constraint vector */ 69 Vec cwork; 70 71 Vec lwork; 72 Vec S; 73 Vec Rwork,Swork,Twork; 74 Vec Av_u; 75 76 KSP solver; 77 PC prec; 78 79 PetscInt ksp_its; 80 PetscInt ksp_its_initial; 81 } AppCtx; 82 83 PetscErrorCode FormFunction(Tao, Vec, PetscReal*, void*); 84 PetscErrorCode FormGradient(Tao, Vec, Vec, void*); 85 PetscErrorCode FormFunctionGradient(Tao, Vec, PetscReal*, Vec, void*); 86 PetscErrorCode FormJacobianState(Tao, Vec, Mat, Mat, Mat, void*); 87 PetscErrorCode FormJacobianDesign(Tao, Vec, Mat, void*); 88 PetscErrorCode FormConstraints(Tao, Vec, Vec, void*); 89 PetscErrorCode FormHessian(Tao, Vec, Mat, Mat, void*); 90 PetscErrorCode Gather(Vec x, Vec state, VecScatter s_scat, Vec design, VecScatter d_scat); 91 PetscErrorCode Scatter(Vec x, Vec state, VecScatter s_scat, Vec design, VecScatter d_scat); 92 PetscErrorCode ParabolicInitialize(AppCtx *user); 93 PetscErrorCode ParabolicDestroy(AppCtx *user); 94 PetscErrorCode ParabolicMonitor(Tao, void*); 95 96 PetscErrorCode StateMatMult(Mat,Vec,Vec); 97 PetscErrorCode StateMatBlockMult(Mat,Vec,Vec); 98 PetscErrorCode StateMatMultTranspose(Mat,Vec,Vec); 99 PetscErrorCode StateMatGetDiagonal(Mat,Vec); 100 PetscErrorCode StateMatDuplicate(Mat,MatDuplicateOption,Mat*); 101 PetscErrorCode StateMatInvMult(Mat,Vec,Vec); 102 PetscErrorCode StateMatInvTransposeMult(Mat,Vec,Vec); 103 PetscErrorCode StateMatBlockPrecMult(PC,Vec,Vec); 104 105 PetscErrorCode DesignMatMult(Mat,Vec,Vec); 106 PetscErrorCode DesignMatMultTranspose(Mat,Vec,Vec); 107 108 PetscErrorCode Gather_i(Vec,Vec*,VecScatter*,PetscInt); 109 PetscErrorCode Scatter_i(Vec,Vec*,VecScatter*,PetscInt); 110 111 static char help[]=""; 112 113 int main(int argc, char **argv) 114 { 115 PetscErrorCode ierr; 116 Vec x,x0; 117 Tao tao; 118 AppCtx user; 119 IS is_allstate,is_alldesign; 120 PetscInt lo,hi,hi2,lo2,ksp_old; 121 PetscInt ntests = 1; 122 PetscInt i; 123 #if defined(PETSC_USE_LOG) 124 PetscLogStage stages[1]; 125 #endif 126 127 PetscCall(PetscInitialize(&argc, &argv, (char*)0,help)); 128 user.mx = 8; 129 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"parabolic example",NULL);PetscCall(ierr); 130 PetscCall(PetscOptionsInt("-mx","Number of grid points in each direction","",user.mx,&user.mx,NULL)); 131 user.nt = 8; 132 PetscCall(PetscOptionsInt("-nt","Number of time steps","",user.nt,&user.nt,NULL)); 133 user.ndata = 64; 134 PetscCall(PetscOptionsInt("-ndata","Numbers of data points per sample","",user.ndata,&user.ndata,NULL)); 135 user.ns = 8; 136 PetscCall(PetscOptionsInt("-ns","Number of samples","",user.ns,&user.ns,NULL)); 137 user.alpha = 1.0; 138 PetscCall(PetscOptionsReal("-alpha","Regularization parameter","",user.alpha,&user.alpha,NULL)); 139 user.beta = 0.01; 140 PetscCall(PetscOptionsReal("-beta","Weight attributed to ||u||^2 in regularization functional","",user.beta,&user.beta,NULL)); 141 user.noise = 0.01; 142 PetscCall(PetscOptionsReal("-noise","Amount of noise to add to data","",user.noise,&user.noise,NULL)); 143 PetscCall(PetscOptionsInt("-ntests","Number of times to repeat TaoSolve","",ntests,&ntests,NULL)); 144 ierr = PetscOptionsEnd();PetscCall(ierr); 145 146 user.m = user.mx*user.mx*user.mx; /* number of constraints per time step */ 147 user.n = user.m*(user.nt+1); /* number of variables */ 148 user.ht = (PetscReal)1/user.nt; 149 150 PetscCall(VecCreate(PETSC_COMM_WORLD,&user.u)); 151 PetscCall(VecCreate(PETSC_COMM_WORLD,&user.y)); 152 PetscCall(VecCreate(PETSC_COMM_WORLD,&user.c)); 153 PetscCall(VecSetSizes(user.u,PETSC_DECIDE,user.n-user.m*user.nt)); 154 PetscCall(VecSetSizes(user.y,PETSC_DECIDE,user.m*user.nt)); 155 PetscCall(VecSetSizes(user.c,PETSC_DECIDE,user.m*user.nt)); 156 PetscCall(VecSetFromOptions(user.u)); 157 PetscCall(VecSetFromOptions(user.y)); 158 PetscCall(VecSetFromOptions(user.c)); 159 160 /* Create scatters for reduced spaces. 161 If the state vector y and design vector u are partitioned as 162 [y_1; y_2; ...; y_np] and [u_1; u_2; ...; u_np] (with np = # of processors), 163 then the solution vector x is organized as 164 [y_1; u_1; y_2; u_2; ...; y_np; u_np]. 165 The index sets user.s_is and user.d_is correspond to the indices of the 166 state and design variables owned by the current processor. 167 */ 168 PetscCall(VecCreate(PETSC_COMM_WORLD,&x)); 169 170 PetscCall(VecGetOwnershipRange(user.y,&lo,&hi)); 171 PetscCall(VecGetOwnershipRange(user.u,&lo2,&hi2)); 172 173 PetscCall(ISCreateStride(PETSC_COMM_SELF,hi-lo,lo,1,&is_allstate)); 174 PetscCall(ISCreateStride(PETSC_COMM_SELF,hi-lo,lo+lo2,1,&user.s_is)); 175 176 PetscCall(ISCreateStride(PETSC_COMM_SELF,hi2-lo2,lo2,1,&is_alldesign)); 177 PetscCall(ISCreateStride(PETSC_COMM_SELF,hi2-lo2,hi+lo2,1,&user.d_is)); 178 179 PetscCall(VecSetSizes(x,hi-lo+hi2-lo2,user.n)); 180 PetscCall(VecSetFromOptions(x)); 181 182 PetscCall(VecScatterCreate(x,user.s_is,user.y,is_allstate,&user.state_scatter)); 183 PetscCall(VecScatterCreate(x,user.d_is,user.u,is_alldesign,&user.design_scatter)); 184 PetscCall(ISDestroy(&is_alldesign)); 185 PetscCall(ISDestroy(&is_allstate)); 186 187 /* Create TAO solver and set desired solution method */ 188 PetscCall(TaoCreate(PETSC_COMM_WORLD,&tao)); 189 PetscCall(TaoSetType(tao,TAOLCL)); 190 191 /* Set up initial vectors and matrices */ 192 PetscCall(ParabolicInitialize(&user)); 193 194 PetscCall(Gather(x,user.y,user.state_scatter,user.u,user.design_scatter)); 195 PetscCall(VecDuplicate(x,&x0)); 196 PetscCall(VecCopy(x,x0)); 197 198 /* Set solution vector with an initial guess */ 199 PetscCall(TaoSetSolution(tao,x)); 200 PetscCall(TaoSetObjective(tao, FormFunction, &user)); 201 PetscCall(TaoSetGradient(tao, NULL, FormGradient, &user)); 202 PetscCall(TaoSetConstraintsRoutine(tao, user.c, FormConstraints, &user)); 203 204 PetscCall(TaoSetJacobianStateRoutine(tao, user.Js, user.JsBlockPrec, user.JsInv, FormJacobianState, &user)); 205 PetscCall(TaoSetJacobianDesignRoutine(tao, user.Jd, FormJacobianDesign, &user)); 206 207 PetscCall(TaoSetFromOptions(tao)); 208 PetscCall(TaoSetStateDesignIS(tao,user.s_is,user.d_is)); 209 210 /* SOLVE THE APPLICATION */ 211 PetscCall(PetscLogStageRegister("Trials",&stages[0])); 212 PetscCall(PetscLogStagePush(stages[0])); 213 user.ksp_its_initial = user.ksp_its; 214 ksp_old = user.ksp_its; 215 for (i=0; i<ntests; i++) { 216 PetscCall(TaoSolve(tao)); 217 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"KSP Iterations = %D\n",user.ksp_its-ksp_old)); 218 PetscCall(VecCopy(x0,x)); 219 PetscCall(TaoSetSolution(tao,x)); 220 } 221 PetscCall(PetscLogStagePop()); 222 PetscCall(PetscBarrier((PetscObject)x)); 223 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"KSP iterations within initialization: ")); 224 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"%D\n",user.ksp_its_initial)); 225 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Total KSP iterations over %D trial(s): ",ntests)); 226 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"%D\n",user.ksp_its)); 227 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"KSP iterations per trial: ")); 228 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"%D\n",(user.ksp_its-user.ksp_its_initial)/ntests)); 229 230 PetscCall(TaoDestroy(&tao)); 231 PetscCall(VecDestroy(&x)); 232 PetscCall(VecDestroy(&x0)); 233 PetscCall(ParabolicDestroy(&user)); 234 235 PetscCall(PetscFinalize()); 236 return 0; 237 } 238 /* ------------------------------------------------------------------- */ 239 /* 240 dwork = Qy - d 241 lwork = L*(u-ur) 242 f = 1/2 * (dwork.dork + alpha*lwork.lwork) 243 */ 244 PetscErrorCode FormFunction(Tao tao,Vec X,PetscReal *f,void *ptr) 245 { 246 PetscReal d1=0,d2=0; 247 PetscInt i,j; 248 AppCtx *user = (AppCtx*)ptr; 249 250 PetscFunctionBegin; 251 PetscCall(Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter)); 252 PetscCall(Scatter_i(user->y,user->yi,user->yi_scatter,user->nt)); 253 for (j=0; j<user->ns; j++) { 254 i = user->sample_times[j]; 255 PetscCall(MatMult(user->Qblock,user->yi[i],user->di[j])); 256 } 257 PetscCall(Gather_i(user->dwork,user->di,user->di_scatter,user->ns)); 258 PetscCall(VecAXPY(user->dwork,-1.0,user->d)); 259 PetscCall(VecDot(user->dwork,user->dwork,&d1)); 260 261 PetscCall(VecWAXPY(user->uwork,-1.0,user->ur,user->u)); 262 PetscCall(MatMult(user->L,user->uwork,user->lwork)); 263 PetscCall(VecDot(user->lwork,user->lwork,&d2)); 264 265 *f = 0.5 * (d1 + user->alpha*d2); 266 PetscFunctionReturn(0); 267 } 268 269 /* ------------------------------------------------------------------- */ 270 /* 271 state: g_s = Q' *(Qy - d) 272 design: g_d = alpha*L'*L*(u-ur) 273 */ 274 PetscErrorCode FormGradient(Tao tao,Vec X,Vec G,void *ptr) 275 { 276 PetscInt i,j; 277 AppCtx *user = (AppCtx*)ptr; 278 279 PetscFunctionBegin; 280 PetscCall(Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter)); 281 PetscCall(Scatter_i(user->y,user->yi,user->yi_scatter,user->nt)); 282 for (j=0; j<user->ns; j++) { 283 i = user->sample_times[j]; 284 PetscCall(MatMult(user->Qblock,user->yi[i],user->di[j])); 285 } 286 PetscCall(Gather_i(user->dwork,user->di,user->di_scatter,user->ns)); 287 PetscCall(VecAXPY(user->dwork,-1.0,user->d)); 288 PetscCall(Scatter_i(user->dwork,user->di,user->di_scatter,user->ns)); 289 PetscCall(VecSet(user->ywork,0.0)); 290 PetscCall(Scatter_i(user->ywork,user->yiwork,user->yi_scatter,user->nt)); 291 for (j=0; j<user->ns; j++) { 292 i = user->sample_times[j]; 293 PetscCall(MatMult(user->QblockT,user->di[j],user->yiwork[i])); 294 } 295 PetscCall(Gather_i(user->ywork,user->yiwork,user->yi_scatter,user->nt)); 296 297 PetscCall(VecWAXPY(user->uwork,-1.0,user->ur,user->u)); 298 PetscCall(MatMult(user->L,user->uwork,user->lwork)); 299 PetscCall(MatMult(user->LT,user->lwork,user->uwork)); 300 PetscCall(VecScale(user->uwork, user->alpha)); 301 PetscCall(Gather(G,user->ywork,user->state_scatter,user->uwork,user->design_scatter)); 302 PetscFunctionReturn(0); 303 } 304 305 PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *f, Vec G, void *ptr) 306 { 307 PetscReal d1,d2; 308 PetscInt i,j; 309 AppCtx *user = (AppCtx*)ptr; 310 311 PetscFunctionBegin; 312 PetscCall(Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter)); 313 PetscCall(Scatter_i(user->y,user->yi,user->yi_scatter,user->nt)); 314 for (j=0; j<user->ns; j++) { 315 i = user->sample_times[j]; 316 PetscCall(MatMult(user->Qblock,user->yi[i],user->di[j])); 317 } 318 PetscCall(Gather_i(user->dwork,user->di,user->di_scatter,user->ns)); 319 PetscCall(VecAXPY(user->dwork,-1.0,user->d)); 320 PetscCall(VecDot(user->dwork,user->dwork,&d1)); 321 PetscCall(Scatter_i(user->dwork,user->di,user->di_scatter,user->ns)); 322 PetscCall(VecSet(user->ywork,0.0)); 323 PetscCall(Scatter_i(user->ywork,user->yiwork,user->yi_scatter,user->nt)); 324 for (j=0; j<user->ns; j++) { 325 i = user->sample_times[j]; 326 PetscCall(MatMult(user->QblockT,user->di[j],user->yiwork[i])); 327 } 328 PetscCall(Gather_i(user->ywork,user->yiwork,user->yi_scatter,user->nt)); 329 330 PetscCall(VecWAXPY(user->uwork,-1.0,user->ur,user->u)); 331 PetscCall(MatMult(user->L,user->uwork,user->lwork)); 332 PetscCall(VecDot(user->lwork,user->lwork,&d2)); 333 PetscCall(MatMult(user->LT,user->lwork,user->uwork)); 334 PetscCall(VecScale(user->uwork, user->alpha)); 335 *f = 0.5 * (d1 + user->alpha*d2); 336 337 PetscCall(Gather(G,user->ywork,user->state_scatter,user->uwork,user->design_scatter)); 338 PetscFunctionReturn(0); 339 } 340 341 /* ------------------------------------------------------------------- */ 342 /* A 343 MatShell object 344 */ 345 PetscErrorCode FormJacobianState(Tao tao, Vec X, Mat J, Mat JPre, Mat JInv, void *ptr) 346 { 347 AppCtx *user = (AppCtx*)ptr; 348 349 PetscFunctionBegin; 350 PetscCall(Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter)); 351 PetscCall(VecSet(user->uwork,0)); 352 PetscCall(VecAXPY(user->uwork,-1.0,user->u)); 353 PetscCall(VecExp(user->uwork)); 354 PetscCall(MatMult(user->Av,user->uwork,user->Av_u)); 355 PetscCall(VecCopy(user->Av_u,user->Swork)); 356 PetscCall(VecReciprocal(user->Swork)); 357 PetscCall(MatCopy(user->Div,user->Divwork,SAME_NONZERO_PATTERN)); 358 PetscCall(MatDiagonalScale(user->Divwork,NULL,user->Swork)); 359 if (user->dsg_formed) { 360 PetscCall(MatProductNumeric(user->DSG)); 361 } else { 362 PetscCall(MatMatMult(user->Divwork,user->Grad,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&user->DSG)); 363 user->dsg_formed = PETSC_TRUE; 364 } 365 366 /* B = speye(nx^3) + ht*DSG; */ 367 PetscCall(MatScale(user->DSG,user->ht)); 368 PetscCall(MatShift(user->DSG,1.0)); 369 PetscFunctionReturn(0); 370 } 371 372 /* ------------------------------------------------------------------- */ 373 /* B */ 374 PetscErrorCode FormJacobianDesign(Tao tao, Vec X, Mat J, void *ptr) 375 { 376 AppCtx *user = (AppCtx*)ptr; 377 378 PetscFunctionBegin; 379 PetscCall(Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter)); 380 PetscFunctionReturn(0); 381 } 382 383 PetscErrorCode StateMatMult(Mat J_shell, Vec X, Vec Y) 384 { 385 PetscInt i; 386 AppCtx *user; 387 388 PetscFunctionBegin; 389 PetscCall(MatShellGetContext(J_shell,&user)); 390 PetscCall(Scatter_i(X,user->yi,user->yi_scatter,user->nt)); 391 PetscCall(MatMult(user->JsBlock,user->yi[0],user->yiwork[0])); 392 for (i=1; i<user->nt; i++) { 393 PetscCall(MatMult(user->JsBlock,user->yi[i],user->yiwork[i])); 394 PetscCall(VecAXPY(user->yiwork[i],-1.0,user->yi[i-1])); 395 } 396 PetscCall(Gather_i(Y,user->yiwork,user->yi_scatter,user->nt)); 397 PetscFunctionReturn(0); 398 } 399 400 PetscErrorCode StateMatMultTranspose(Mat J_shell, Vec X, Vec Y) 401 { 402 PetscInt i; 403 AppCtx *user; 404 405 PetscFunctionBegin; 406 PetscCall(MatShellGetContext(J_shell,&user)); 407 PetscCall(Scatter_i(X,user->yi,user->yi_scatter,user->nt)); 408 for (i=0; i<user->nt-1; i++) { 409 PetscCall(MatMult(user->JsBlock,user->yi[i],user->yiwork[i])); 410 PetscCall(VecAXPY(user->yiwork[i],-1.0,user->yi[i+1])); 411 } 412 i = user->nt-1; 413 PetscCall(MatMult(user->JsBlock,user->yi[i],user->yiwork[i])); 414 PetscCall(Gather_i(Y,user->yiwork,user->yi_scatter,user->nt)); 415 PetscFunctionReturn(0); 416 } 417 418 PetscErrorCode StateMatBlockMult(Mat J_shell, Vec X, Vec Y) 419 { 420 AppCtx *user; 421 422 PetscFunctionBegin; 423 PetscCall(MatShellGetContext(J_shell,&user)); 424 PetscCall(MatMult(user->Grad,X,user->Swork)); 425 PetscCall(VecPointwiseDivide(user->Swork,user->Swork,user->Av_u)); 426 PetscCall(MatMult(user->Div,user->Swork,Y)); 427 PetscCall(VecAYPX(Y,user->ht,X)); 428 PetscFunctionReturn(0); 429 } 430 431 PetscErrorCode DesignMatMult(Mat J_shell, Vec X, Vec Y) 432 { 433 PetscInt i; 434 AppCtx *user; 435 436 PetscFunctionBegin; 437 PetscCall(MatShellGetContext(J_shell,&user)); 438 439 /* sdiag(1./v) */ 440 PetscCall(VecSet(user->uwork,0)); 441 PetscCall(VecAXPY(user->uwork,-1.0,user->u)); 442 PetscCall(VecExp(user->uwork)); 443 444 /* sdiag(1./((Av*(1./v)).^2)) */ 445 PetscCall(MatMult(user->Av,user->uwork,user->Swork)); 446 PetscCall(VecPointwiseMult(user->Swork,user->Swork,user->Swork)); 447 PetscCall(VecReciprocal(user->Swork)); 448 449 /* (Av * (sdiag(1./v) * b)) */ 450 PetscCall(VecPointwiseMult(user->uwork,user->uwork,X)); 451 PetscCall(MatMult(user->Av,user->uwork,user->Twork)); 452 453 /* (sdiag(1./((Av*(1./v)).^2)) * (Av * (sdiag(1./v) * b))) */ 454 PetscCall(VecPointwiseMult(user->Swork,user->Twork,user->Swork)); 455 456 PetscCall(Scatter_i(user->y,user->yi,user->yi_scatter,user->nt)); 457 for (i=0; i<user->nt; i++) { 458 /* (sdiag(Grad*y(:,i)) */ 459 PetscCall(MatMult(user->Grad,user->yi[i],user->Twork)); 460 461 /* ht * Div * (sdiag(Grad*y(:,i)) * (sdiag(1./((Av*(1./v)).^2)) * (Av * (sdiag(1./v) * b)))) */ 462 PetscCall(VecPointwiseMult(user->Twork,user->Twork,user->Swork)); 463 PetscCall(MatMult(user->Div,user->Twork,user->yiwork[i])); 464 PetscCall(VecScale(user->yiwork[i],user->ht)); 465 } 466 PetscCall(Gather_i(Y,user->yiwork,user->yi_scatter,user->nt)); 467 468 PetscFunctionReturn(0); 469 } 470 471 PetscErrorCode DesignMatMultTranspose(Mat J_shell, Vec X, Vec Y) 472 { 473 PetscInt i; 474 AppCtx *user; 475 476 PetscFunctionBegin; 477 PetscCall(MatShellGetContext(J_shell,&user)); 478 479 /* sdiag(1./((Av*(1./v)).^2)) */ 480 PetscCall(VecSet(user->uwork,0)); 481 PetscCall(VecAXPY(user->uwork,-1.0,user->u)); 482 PetscCall(VecExp(user->uwork)); 483 PetscCall(MatMult(user->Av,user->uwork,user->Rwork)); 484 PetscCall(VecPointwiseMult(user->Rwork,user->Rwork,user->Rwork)); 485 PetscCall(VecReciprocal(user->Rwork)); 486 487 PetscCall(VecSet(Y,0.0)); 488 PetscCall(Scatter_i(user->y,user->yi,user->yi_scatter,user->nt)); 489 PetscCall(Scatter_i(X,user->yiwork,user->yi_scatter,user->nt)); 490 for (i=0; i<user->nt; i++) { 491 /* (Div' * b(:,i)) */ 492 PetscCall(MatMult(user->Grad,user->yiwork[i],user->Swork)); 493 494 /* sdiag(Grad*y(:,i)) */ 495 PetscCall(MatMult(user->Grad,user->yi[i],user->Twork)); 496 497 /* (sdiag(Grad*y(:,i)) * (Div' * b(:,i))) */ 498 PetscCall(VecPointwiseMult(user->Twork,user->Swork,user->Twork)); 499 500 /* (sdiag(1./((Av*(1./v)).^2)) * (sdiag(Grad*y(:,i)) * (Div' * b(:,i)))) */ 501 PetscCall(VecPointwiseMult(user->Twork,user->Rwork,user->Twork)); 502 503 /* (Av' * (sdiag(1./((Av*(1./v)).^2)) * (sdiag(Grad*y(:,i)) * (Div' * b(:,i))))) */ 504 PetscCall(MatMult(user->AvT,user->Twork,user->yiwork[i])); 505 506 /* sdiag(1./v) * (Av' * (sdiag(1./((Av*(1./v)).^2)) * (sdiag(Grad*y(:,i)) * (Div' * b(:,i))))) */ 507 PetscCall(VecPointwiseMult(user->yiwork[i],user->uwork,user->yiwork[i])); 508 PetscCall(VecAXPY(Y,user->ht,user->yiwork[i])); 509 } 510 PetscFunctionReturn(0); 511 } 512 513 PetscErrorCode StateMatBlockPrecMult(PC PC_shell, Vec X, Vec Y) 514 { 515 AppCtx *user; 516 517 PetscFunctionBegin; 518 PetscCall(PCShellGetContext(PC_shell,&user)); 519 520 if (user->dsg_formed) { 521 PetscCall(MatSOR(user->DSG,X,1.0,(MatSORType)(SOR_ZERO_INITIAL_GUESS | SOR_LOCAL_SYMMETRIC_SWEEP),0.0,1,1,Y)); 522 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"DSG not formed"); 523 PetscFunctionReturn(0); 524 } 525 526 PetscErrorCode StateMatInvMult(Mat J_shell, Vec X, Vec Y) 527 { 528 AppCtx *user; 529 PetscInt its,i; 530 531 PetscFunctionBegin; 532 PetscCall(MatShellGetContext(J_shell,&user)); 533 534 if (Y == user->ytrue) { 535 /* First solve is done with true solution to set up problem */ 536 PetscCall(KSPSetTolerances(user->solver,1e-8,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT)); 537 } else { 538 PetscCall(KSPSetTolerances(user->solver,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT)); 539 } 540 541 PetscCall(Scatter_i(X,user->yi,user->yi_scatter,user->nt)); 542 PetscCall(KSPSolve(user->solver,user->yi[0],user->yiwork[0])); 543 PetscCall(KSPGetIterationNumber(user->solver,&its)); 544 user->ksp_its = user->ksp_its + its; 545 546 for (i=1; i<user->nt; i++) { 547 PetscCall(VecAXPY(user->yi[i],1.0,user->yiwork[i-1])); 548 PetscCall(KSPSolve(user->solver,user->yi[i],user->yiwork[i])); 549 PetscCall(KSPGetIterationNumber(user->solver,&its)); 550 user->ksp_its = user->ksp_its + its; 551 } 552 PetscCall(Gather_i(Y,user->yiwork,user->yi_scatter,user->nt)); 553 PetscFunctionReturn(0); 554 } 555 556 PetscErrorCode StateMatInvTransposeMult(Mat J_shell, Vec X, Vec Y) 557 { 558 AppCtx *user; 559 PetscInt its,i; 560 561 PetscFunctionBegin; 562 PetscCall(MatShellGetContext(J_shell,&user)); 563 564 PetscCall(Scatter_i(X,user->yi,user->yi_scatter,user->nt)); 565 566 i = user->nt - 1; 567 PetscCall(KSPSolve(user->solver,user->yi[i],user->yiwork[i])); 568 569 PetscCall(KSPGetIterationNumber(user->solver,&its)); 570 user->ksp_its = user->ksp_its + its; 571 572 for (i=user->nt-2; i>=0; i--) { 573 PetscCall(VecAXPY(user->yi[i],1.0,user->yiwork[i+1])); 574 PetscCall(KSPSolve(user->solver,user->yi[i],user->yiwork[i])); 575 576 PetscCall(KSPGetIterationNumber(user->solver,&its)); 577 user->ksp_its = user->ksp_its + its; 578 579 } 580 581 PetscCall(Gather_i(Y,user->yiwork,user->yi_scatter,user->nt)); 582 PetscFunctionReturn(0); 583 } 584 585 PetscErrorCode StateMatDuplicate(Mat J_shell, MatDuplicateOption opt, Mat *new_shell) 586 { 587 AppCtx *user; 588 589 PetscFunctionBegin; 590 PetscCall(MatShellGetContext(J_shell,&user)); 591 592 PetscCall(MatCreateShell(PETSC_COMM_WORLD,PETSC_DETERMINE,PETSC_DETERMINE,user->m,user->m,user,new_shell)); 593 PetscCall(MatShellSetOperation(*new_shell,MATOP_MULT,(void(*)(void))StateMatMult)); 594 PetscCall(MatShellSetOperation(*new_shell,MATOP_DUPLICATE,(void(*)(void))StateMatDuplicate)); 595 PetscCall(MatShellSetOperation(*new_shell,MATOP_MULT_TRANSPOSE,(void(*)(void))StateMatMultTranspose)); 596 PetscCall(MatShellSetOperation(*new_shell,MATOP_GET_DIAGONAL,(void(*)(void))StateMatGetDiagonal)); 597 PetscFunctionReturn(0); 598 } 599 600 PetscErrorCode StateMatGetDiagonal(Mat J_shell, Vec X) 601 { 602 AppCtx *user; 603 604 PetscFunctionBegin; 605 PetscCall(MatShellGetContext(J_shell,&user)); 606 PetscCall(VecCopy(user->js_diag,X)); 607 PetscFunctionReturn(0); 608 609 } 610 611 PetscErrorCode FormConstraints(Tao tao, Vec X, Vec C, void *ptr) 612 { 613 /* con = Ay - q, A = [B 0 0 ... 0; 614 -I B 0 ... 0; 615 0 -I B ... 0; 616 ... ; 617 0 ... -I B] 618 B = ht * Div * Sigma * Grad + eye */ 619 PetscInt i; 620 AppCtx *user = (AppCtx*)ptr; 621 622 PetscFunctionBegin; 623 PetscCall(Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter)); 624 PetscCall(Scatter_i(user->y,user->yi,user->yi_scatter,user->nt)); 625 PetscCall(MatMult(user->JsBlock,user->yi[0],user->yiwork[0])); 626 for (i=1; i<user->nt; i++) { 627 PetscCall(MatMult(user->JsBlock,user->yi[i],user->yiwork[i])); 628 PetscCall(VecAXPY(user->yiwork[i],-1.0,user->yi[i-1])); 629 } 630 PetscCall(Gather_i(C,user->yiwork,user->yi_scatter,user->nt)); 631 PetscCall(VecAXPY(C,-1.0,user->q)); 632 PetscFunctionReturn(0); 633 } 634 635 PetscErrorCode Scatter(Vec x, Vec state, VecScatter s_scat, Vec design, VecScatter d_scat) 636 { 637 PetscFunctionBegin; 638 PetscCall(VecScatterBegin(s_scat,x,state,INSERT_VALUES,SCATTER_FORWARD)); 639 PetscCall(VecScatterEnd(s_scat,x,state,INSERT_VALUES,SCATTER_FORWARD)); 640 PetscCall(VecScatterBegin(d_scat,x,design,INSERT_VALUES,SCATTER_FORWARD)); 641 PetscCall(VecScatterEnd(d_scat,x,design,INSERT_VALUES,SCATTER_FORWARD)); 642 PetscFunctionReturn(0); 643 } 644 645 PetscErrorCode Scatter_i(Vec y, Vec *yi, VecScatter *scat, PetscInt nt) 646 { 647 PetscInt i; 648 649 PetscFunctionBegin; 650 for (i=0; i<nt; i++) { 651 PetscCall(VecScatterBegin(scat[i],y,yi[i],INSERT_VALUES,SCATTER_FORWARD)); 652 PetscCall(VecScatterEnd(scat[i],y,yi[i],INSERT_VALUES,SCATTER_FORWARD)); 653 } 654 PetscFunctionReturn(0); 655 } 656 657 PetscErrorCode Gather(Vec x, Vec state, VecScatter s_scat, Vec design, VecScatter d_scat) 658 { 659 PetscFunctionBegin; 660 PetscCall(VecScatterBegin(s_scat,state,x,INSERT_VALUES,SCATTER_REVERSE)); 661 PetscCall(VecScatterEnd(s_scat,state,x,INSERT_VALUES,SCATTER_REVERSE)); 662 PetscCall(VecScatterBegin(d_scat,design,x,INSERT_VALUES,SCATTER_REVERSE)); 663 PetscCall(VecScatterEnd(d_scat,design,x,INSERT_VALUES,SCATTER_REVERSE)); 664 PetscFunctionReturn(0); 665 } 666 667 PetscErrorCode Gather_i(Vec y, Vec *yi, VecScatter *scat, PetscInt nt) 668 { 669 PetscInt i; 670 671 PetscFunctionBegin; 672 for (i=0; i<nt; i++) { 673 PetscCall(VecScatterBegin(scat[i],yi[i],y,INSERT_VALUES,SCATTER_REVERSE)); 674 PetscCall(VecScatterEnd(scat[i],yi[i],y,INSERT_VALUES,SCATTER_REVERSE)); 675 } 676 PetscFunctionReturn(0); 677 } 678 679 PetscErrorCode ParabolicInitialize(AppCtx *user) 680 { 681 PetscInt m,n,i,j,k,linear_index,istart,iend,iblock,lo,hi,lo2,hi2; 682 Vec XX,YY,ZZ,XXwork,YYwork,ZZwork,UTwork,yi,di,bc; 683 PetscReal *x, *y, *z; 684 PetscReal h,stime; 685 PetscScalar hinv,neg_hinv,half = 0.5,sqrt_beta; 686 PetscInt im,indx1,indx2,indy1,indy2,indz1,indz2,nx,ny,nz; 687 PetscReal xri,yri,zri,xim,yim,zim,dx1,dx2,dy1,dy2,dz1,dz2,Dx,Dy,Dz; 688 PetscScalar v,vx,vy,vz; 689 IS is_from_y,is_to_yi,is_from_d,is_to_di; 690 /* Data locations */ 691 PetscScalar xr[64] = {0.4970, 0.8498, 0.7814, 0.6268, 0.7782, 0.6402, 0.3617, 0.3160, 692 0.3610, 0.5298, 0.6987, 0.3331, 0.7962, 0.5596, 0.3866, 0.6774, 693 0.5407, 0.4518, 0.6702, 0.6061, 0.7580, 0.8997, 0.5198, 0.8326, 694 0.2138, 0.9198, 0.3000, 0.2833, 0.8288, 0.7076, 0.1820, 0.0728, 695 0.8447, 0.2367, 0.3239, 0.6413, 0.3114, 0.4731, 0.1192, 0.9273, 696 0.5724, 0.4331, 0.5136, 0.3547, 0.4413, 0.2602, 0.5698, 0.7278, 697 0.5261, 0.6230, 0.2454, 0.3948, 0.7479, 0.6582, 0.4660, 0.5594, 698 0.7574, 0.1143, 0.5900, 0.1065, 0.4260, 0.3294, 0.8276, 0.0756}; 699 700 PetscScalar yr[64] = {0.7345, 0.9120, 0.9288, 0.7528, 0.4463, 0.4985, 0.2497, 0.6256, 701 0.3425, 0.9026, 0.6983, 0.4230, 0.7140, 0.2970, 0.4474, 0.8792, 702 0.6604, 0.2485, 0.7968, 0.6127, 0.1796, 0.2437, 0.5938, 0.6137, 703 0.3867, 0.5658, 0.4575, 0.1009, 0.0863, 0.3361, 0.0738, 0.3985, 704 0.6602, 0.1437, 0.0934, 0.5983, 0.5950, 0.0763, 0.0768, 0.2288, 705 0.5761, 0.1129, 0.3841, 0.6150, 0.6904, 0.6686, 0.1361, 0.4601, 706 0.4491, 0.3716, 0.1969, 0.6537, 0.6743, 0.6991, 0.4811, 0.5480, 707 0.1684, 0.4569, 0.6889, 0.8437, 0.3015, 0.2854, 0.8199, 0.2658}; 708 709 PetscScalar zr[64] = {0.7668, 0.8573, 0.2654, 0.2719, 0.1060, 0.1311, 0.6232, 0.2295, 710 0.8009, 0.2147, 0.2119, 0.9325, 0.4473, 0.3600, 0.3374, 0.3819, 711 0.4066, 0.5801, 0.1673, 0.0959, 0.4638, 0.8236, 0.8800, 0.2939, 712 0.2028, 0.8262, 0.2706, 0.6276, 0.9085, 0.6443, 0.8241, 0.0712, 713 0.1824, 0.7789, 0.4389, 0.8415, 0.7055, 0.6639, 0.3653, 0.2078, 714 0.1987, 0.2297, 0.4321, 0.8115, 0.4915, 0.7764, 0.4657, 0.4627, 715 0.4569, 0.4232, 0.8514, 0.0674, 0.3227, 0.1055, 0.6690, 0.6313, 716 0.9226, 0.5461, 0.4126, 0.2364, 0.6096, 0.7042, 0.3914, 0.0711}; 717 718 PetscFunctionBegin; 719 PetscCall(PetscMalloc1(user->mx,&x)); 720 PetscCall(PetscMalloc1(user->mx,&y)); 721 PetscCall(PetscMalloc1(user->mx,&z)); 722 user->jformed = PETSC_FALSE; 723 user->dsg_formed = PETSC_FALSE; 724 725 n = user->mx * user->mx * user->mx; 726 m = 3 * user->mx * user->mx * (user->mx-1); 727 sqrt_beta = PetscSqrtScalar(user->beta); 728 729 user->ksp_its = 0; 730 user->ksp_its_initial = 0; 731 732 stime = (PetscReal)user->nt/user->ns; 733 PetscCall(PetscMalloc1(user->ns,&user->sample_times)); 734 for (i=0; i<user->ns; i++) { 735 user->sample_times[i] = (PetscInt)(stime*((PetscReal)i+1.0)-0.5); 736 } 737 738 PetscCall(VecCreate(PETSC_COMM_WORLD,&XX)); 739 PetscCall(VecCreate(PETSC_COMM_WORLD,&user->q)); 740 PetscCall(VecSetSizes(XX,PETSC_DECIDE,n)); 741 PetscCall(VecSetSizes(user->q,PETSC_DECIDE,n*user->nt)); 742 PetscCall(VecSetFromOptions(XX)); 743 PetscCall(VecSetFromOptions(user->q)); 744 745 PetscCall(VecDuplicate(XX,&YY)); 746 PetscCall(VecDuplicate(XX,&ZZ)); 747 PetscCall(VecDuplicate(XX,&XXwork)); 748 PetscCall(VecDuplicate(XX,&YYwork)); 749 PetscCall(VecDuplicate(XX,&ZZwork)); 750 PetscCall(VecDuplicate(XX,&UTwork)); 751 PetscCall(VecDuplicate(XX,&user->utrue)); 752 PetscCall(VecDuplicate(XX,&bc)); 753 754 /* Generate 3D grid, and collect ns (1<=ns<=8) right-hand-side vectors into user->q */ 755 h = 1.0/user->mx; 756 hinv = user->mx; 757 neg_hinv = -hinv; 758 759 PetscCall(VecGetOwnershipRange(XX,&istart,&iend)); 760 for (linear_index=istart; linear_index<iend; linear_index++) { 761 i = linear_index % user->mx; 762 j = ((linear_index-i)/user->mx) % user->mx; 763 k = ((linear_index-i)/user->mx-j) / user->mx; 764 vx = h*(i+0.5); 765 vy = h*(j+0.5); 766 vz = h*(k+0.5); 767 PetscCall(VecSetValues(XX,1,&linear_index,&vx,INSERT_VALUES)); 768 PetscCall(VecSetValues(YY,1,&linear_index,&vy,INSERT_VALUES)); 769 PetscCall(VecSetValues(ZZ,1,&linear_index,&vz,INSERT_VALUES)); 770 if ((vx<0.6) && (vx>0.4) && (vy<0.6) && (vy>0.4) && (vy<0.6) && (vz<0.6) && (vz>0.4)) { 771 v = 1000.0; 772 PetscCall(VecSetValues(bc,1,&linear_index,&v,INSERT_VALUES)); 773 } 774 } 775 776 PetscCall(VecAssemblyBegin(XX)); 777 PetscCall(VecAssemblyEnd(XX)); 778 PetscCall(VecAssemblyBegin(YY)); 779 PetscCall(VecAssemblyEnd(YY)); 780 PetscCall(VecAssemblyBegin(ZZ)); 781 PetscCall(VecAssemblyEnd(ZZ)); 782 PetscCall(VecAssemblyBegin(bc)); 783 PetscCall(VecAssemblyEnd(bc)); 784 785 /* Compute true parameter function 786 ut = 0.5 + exp(-10*((x-0.5)^2+(y-0.5)^2+(z-0.5)^2)) */ 787 PetscCall(VecCopy(XX,XXwork)); 788 PetscCall(VecCopy(YY,YYwork)); 789 PetscCall(VecCopy(ZZ,ZZwork)); 790 791 PetscCall(VecShift(XXwork,-0.5)); 792 PetscCall(VecShift(YYwork,-0.5)); 793 PetscCall(VecShift(ZZwork,-0.5)); 794 795 PetscCall(VecPointwiseMult(XXwork,XXwork,XXwork)); 796 PetscCall(VecPointwiseMult(YYwork,YYwork,YYwork)); 797 PetscCall(VecPointwiseMult(ZZwork,ZZwork,ZZwork)); 798 799 PetscCall(VecCopy(XXwork,user->utrue)); 800 PetscCall(VecAXPY(user->utrue,1.0,YYwork)); 801 PetscCall(VecAXPY(user->utrue,1.0,ZZwork)); 802 PetscCall(VecScale(user->utrue,-10.0)); 803 PetscCall(VecExp(user->utrue)); 804 PetscCall(VecShift(user->utrue,0.5)); 805 806 PetscCall(VecDestroy(&XX)); 807 PetscCall(VecDestroy(&YY)); 808 PetscCall(VecDestroy(&ZZ)); 809 PetscCall(VecDestroy(&XXwork)); 810 PetscCall(VecDestroy(&YYwork)); 811 PetscCall(VecDestroy(&ZZwork)); 812 PetscCall(VecDestroy(&UTwork)); 813 814 /* Initial guess and reference model */ 815 PetscCall(VecDuplicate(user->utrue,&user->ur)); 816 PetscCall(VecSet(user->ur,0.0)); 817 818 /* Generate Grad matrix */ 819 PetscCall(MatCreate(PETSC_COMM_WORLD,&user->Grad)); 820 PetscCall(MatSetSizes(user->Grad,PETSC_DECIDE,PETSC_DECIDE,m,n)); 821 PetscCall(MatSetFromOptions(user->Grad)); 822 PetscCall(MatMPIAIJSetPreallocation(user->Grad,2,NULL,2,NULL)); 823 PetscCall(MatSeqAIJSetPreallocation(user->Grad,2,NULL)); 824 PetscCall(MatGetOwnershipRange(user->Grad,&istart,&iend)); 825 826 for (i=istart; i<iend; i++) { 827 if (i<m/3) { 828 iblock = i / (user->mx-1); 829 j = iblock*user->mx + (i % (user->mx-1)); 830 PetscCall(MatSetValues(user->Grad,1,&i,1,&j,&neg_hinv,INSERT_VALUES)); 831 j = j+1; 832 PetscCall(MatSetValues(user->Grad,1,&i,1,&j,&hinv,INSERT_VALUES)); 833 } 834 if (i>=m/3 && i<2*m/3) { 835 iblock = (i-m/3) / (user->mx*(user->mx-1)); 836 j = iblock*user->mx*user->mx + ((i-m/3) % (user->mx*(user->mx-1))); 837 PetscCall(MatSetValues(user->Grad,1,&i,1,&j,&neg_hinv,INSERT_VALUES)); 838 j = j + user->mx; 839 PetscCall(MatSetValues(user->Grad,1,&i,1,&j,&hinv,INSERT_VALUES)); 840 } 841 if (i>=2*m/3) { 842 j = i-2*m/3; 843 PetscCall(MatSetValues(user->Grad,1,&i,1,&j,&neg_hinv,INSERT_VALUES)); 844 j = j + user->mx*user->mx; 845 PetscCall(MatSetValues(user->Grad,1,&i,1,&j,&hinv,INSERT_VALUES)); 846 } 847 } 848 849 PetscCall(MatAssemblyBegin(user->Grad,MAT_FINAL_ASSEMBLY)); 850 PetscCall(MatAssemblyEnd(user->Grad,MAT_FINAL_ASSEMBLY)); 851 852 /* Generate arithmetic averaging matrix Av */ 853 PetscCall(MatCreate(PETSC_COMM_WORLD,&user->Av)); 854 PetscCall(MatSetSizes(user->Av,PETSC_DECIDE,PETSC_DECIDE,m,n)); 855 PetscCall(MatSetFromOptions(user->Av)); 856 PetscCall(MatMPIAIJSetPreallocation(user->Av,2,NULL,2,NULL)); 857 PetscCall(MatSeqAIJSetPreallocation(user->Av,2,NULL)); 858 PetscCall(MatGetOwnershipRange(user->Av,&istart,&iend)); 859 860 for (i=istart; i<iend; i++) { 861 if (i<m/3) { 862 iblock = i / (user->mx-1); 863 j = iblock*user->mx + (i % (user->mx-1)); 864 PetscCall(MatSetValues(user->Av,1,&i,1,&j,&half,INSERT_VALUES)); 865 j = j+1; 866 PetscCall(MatSetValues(user->Av,1,&i,1,&j,&half,INSERT_VALUES)); 867 } 868 if (i>=m/3 && i<2*m/3) { 869 iblock = (i-m/3) / (user->mx*(user->mx-1)); 870 j = iblock*user->mx*user->mx + ((i-m/3) % (user->mx*(user->mx-1))); 871 PetscCall(MatSetValues(user->Av,1,&i,1,&j,&half,INSERT_VALUES)); 872 j = j + user->mx; 873 PetscCall(MatSetValues(user->Av,1,&i,1,&j,&half,INSERT_VALUES)); 874 } 875 if (i>=2*m/3) { 876 j = i-2*m/3; 877 PetscCall(MatSetValues(user->Av,1,&i,1,&j,&half,INSERT_VALUES)); 878 j = j + user->mx*user->mx; 879 PetscCall(MatSetValues(user->Av,1,&i,1,&j,&half,INSERT_VALUES)); 880 } 881 } 882 883 PetscCall(MatAssemblyBegin(user->Av,MAT_FINAL_ASSEMBLY)); 884 PetscCall(MatAssemblyEnd(user->Av,MAT_FINAL_ASSEMBLY)); 885 886 /* Generate transpose of averaging matrix Av */ 887 PetscCall(MatTranspose(user->Av,MAT_INITIAL_MATRIX,&user->AvT)); 888 889 PetscCall(MatCreate(PETSC_COMM_WORLD,&user->L)); 890 PetscCall(MatSetSizes(user->L,PETSC_DECIDE,PETSC_DECIDE,m+n,n)); 891 PetscCall(MatSetFromOptions(user->L)); 892 PetscCall(MatMPIAIJSetPreallocation(user->L,2,NULL,2,NULL)); 893 PetscCall(MatSeqAIJSetPreallocation(user->L,2,NULL)); 894 PetscCall(MatGetOwnershipRange(user->L,&istart,&iend)); 895 896 for (i=istart; i<iend; i++) { 897 if (i<m/3) { 898 iblock = i / (user->mx-1); 899 j = iblock*user->mx + (i % (user->mx-1)); 900 PetscCall(MatSetValues(user->L,1,&i,1,&j,&neg_hinv,INSERT_VALUES)); 901 j = j+1; 902 PetscCall(MatSetValues(user->L,1,&i,1,&j,&hinv,INSERT_VALUES)); 903 } 904 if (i>=m/3 && i<2*m/3) { 905 iblock = (i-m/3) / (user->mx*(user->mx-1)); 906 j = iblock*user->mx*user->mx + ((i-m/3) % (user->mx*(user->mx-1))); 907 PetscCall(MatSetValues(user->L,1,&i,1,&j,&neg_hinv,INSERT_VALUES)); 908 j = j + user->mx; 909 PetscCall(MatSetValues(user->L,1,&i,1,&j,&hinv,INSERT_VALUES)); 910 } 911 if (i>=2*m/3 && i<m) { 912 j = i-2*m/3; 913 PetscCall(MatSetValues(user->L,1,&i,1,&j,&neg_hinv,INSERT_VALUES)); 914 j = j + user->mx*user->mx; 915 PetscCall(MatSetValues(user->L,1,&i,1,&j,&hinv,INSERT_VALUES)); 916 } 917 if (i>=m) { 918 j = i - m; 919 PetscCall(MatSetValues(user->L,1,&i,1,&j,&sqrt_beta,INSERT_VALUES)); 920 } 921 } 922 923 PetscCall(MatAssemblyBegin(user->L,MAT_FINAL_ASSEMBLY)); 924 PetscCall(MatAssemblyEnd(user->L,MAT_FINAL_ASSEMBLY)); 925 926 PetscCall(MatScale(user->L,PetscPowScalar(h,1.5))); 927 928 /* Generate Div matrix */ 929 PetscCall(MatTranspose(user->Grad,MAT_INITIAL_MATRIX,&user->Div)); 930 931 /* Build work vectors and matrices */ 932 PetscCall(VecCreate(PETSC_COMM_WORLD,&user->S)); 933 PetscCall(VecSetSizes(user->S, PETSC_DECIDE, user->mx*user->mx*(user->mx-1)*3)); 934 PetscCall(VecSetFromOptions(user->S)); 935 936 PetscCall(VecCreate(PETSC_COMM_WORLD,&user->lwork)); 937 PetscCall(VecSetSizes(user->lwork,PETSC_DECIDE,m+user->mx*user->mx*user->mx)); 938 PetscCall(VecSetFromOptions(user->lwork)); 939 940 PetscCall(MatDuplicate(user->Div,MAT_SHARE_NONZERO_PATTERN,&user->Divwork)); 941 PetscCall(MatDuplicate(user->Av,MAT_SHARE_NONZERO_PATTERN,&user->Avwork)); 942 943 PetscCall(VecCreate(PETSC_COMM_WORLD,&user->d)); 944 PetscCall(VecSetSizes(user->d,PETSC_DECIDE,user->ndata*user->nt)); 945 PetscCall(VecSetFromOptions(user->d)); 946 947 PetscCall(VecDuplicate(user->S,&user->Swork)); 948 PetscCall(VecDuplicate(user->S,&user->Av_u)); 949 PetscCall(VecDuplicate(user->S,&user->Twork)); 950 PetscCall(VecDuplicate(user->S,&user->Rwork)); 951 PetscCall(VecDuplicate(user->y,&user->ywork)); 952 PetscCall(VecDuplicate(user->u,&user->uwork)); 953 PetscCall(VecDuplicate(user->u,&user->js_diag)); 954 PetscCall(VecDuplicate(user->c,&user->cwork)); 955 PetscCall(VecDuplicate(user->d,&user->dwork)); 956 957 /* Create matrix-free shell user->Js for computing A*x */ 958 PetscCall(MatCreateShell(PETSC_COMM_WORLD,PETSC_DETERMINE,PETSC_DETERMINE,user->m*user->nt,user->m*user->nt,user,&user->Js)); 959 PetscCall(MatShellSetOperation(user->Js,MATOP_MULT,(void(*)(void))StateMatMult)); 960 PetscCall(MatShellSetOperation(user->Js,MATOP_DUPLICATE,(void(*)(void))StateMatDuplicate)); 961 PetscCall(MatShellSetOperation(user->Js,MATOP_MULT_TRANSPOSE,(void(*)(void))StateMatMultTranspose)); 962 PetscCall(MatShellSetOperation(user->Js,MATOP_GET_DIAGONAL,(void(*)(void))StateMatGetDiagonal)); 963 964 /* Diagonal blocks of user->Js */ 965 PetscCall(MatCreateShell(PETSC_COMM_WORLD,PETSC_DETERMINE,PETSC_DETERMINE,user->m,user->m,user,&user->JsBlock)); 966 PetscCall(MatShellSetOperation(user->JsBlock,MATOP_MULT,(void(*)(void))StateMatBlockMult)); 967 /* Blocks are symmetric */ 968 PetscCall(MatShellSetOperation(user->JsBlock,MATOP_MULT_TRANSPOSE,(void(*)(void))StateMatBlockMult)); 969 970 /* Create a matrix-free shell user->JsBlockPrec for computing (U+D)\D*(L+D)\x, where JsBlock = L+D+U, 971 D is diagonal, L is strictly lower triangular, and U is strictly upper triangular. 972 This is an SSOR preconditioner for user->JsBlock. */ 973 PetscCall(MatCreateShell(PETSC_COMM_WORLD,PETSC_DETERMINE,PETSC_DETERMINE,user->m,user->m,user,&user->JsBlockPrec)); 974 PetscCall(MatShellSetOperation(user->JsBlockPrec,MATOP_MULT,(void(*)(void))StateMatBlockPrecMult)); 975 /* JsBlockPrec is symmetric */ 976 PetscCall(MatShellSetOperation(user->JsBlockPrec,MATOP_MULT_TRANSPOSE,(void(*)(void))StateMatBlockPrecMult)); 977 PetscCall(MatSetOption(user->JsBlockPrec,MAT_SYMMETRY_ETERNAL,PETSC_TRUE)); 978 979 /* Create a matrix-free shell user->Jd for computing B*x */ 980 PetscCall(MatCreateShell(PETSC_COMM_WORLD,PETSC_DETERMINE,PETSC_DETERMINE,user->m*user->nt,user->m,user,&user->Jd)); 981 PetscCall(MatShellSetOperation(user->Jd,MATOP_MULT,(void(*)(void))DesignMatMult)); 982 PetscCall(MatShellSetOperation(user->Jd,MATOP_MULT_TRANSPOSE,(void(*)(void))DesignMatMultTranspose)); 983 984 /* User-defined routines for computing user->Js\x and user->Js^T\x*/ 985 PetscCall(MatCreateShell(PETSC_COMM_WORLD,PETSC_DETERMINE,PETSC_DETERMINE,user->m*user->nt,user->m*user->nt,user,&user->JsInv)); 986 PetscCall(MatShellSetOperation(user->JsInv,MATOP_MULT,(void(*)(void))StateMatInvMult)); 987 PetscCall(MatShellSetOperation(user->JsInv,MATOP_MULT_TRANSPOSE,(void(*)(void))StateMatInvTransposeMult)); 988 989 /* Solver options and tolerances */ 990 PetscCall(KSPCreate(PETSC_COMM_WORLD,&user->solver)); 991 PetscCall(KSPSetType(user->solver,KSPCG)); 992 PetscCall(KSPSetOperators(user->solver,user->JsBlock,user->JsBlockPrec)); 993 PetscCall(KSPSetInitialGuessNonzero(user->solver,PETSC_FALSE)); 994 PetscCall(KSPSetTolerances(user->solver,1e-4,1e-20,1e3,500)); 995 PetscCall(KSPSetFromOptions(user->solver)); 996 PetscCall(KSPGetPC(user->solver,&user->prec)); 997 PetscCall(PCSetType(user->prec,PCSHELL)); 998 999 PetscCall(PCShellSetApply(user->prec,StateMatBlockPrecMult)); 1000 PetscCall(PCShellSetApplyTranspose(user->prec,StateMatBlockPrecMult)); 1001 PetscCall(PCShellSetContext(user->prec,user)); 1002 1003 /* Create scatter from y to y_1,y_2,...,y_nt */ 1004 PetscCall(PetscMalloc1(user->nt*user->m,&user->yi_scatter)); 1005 PetscCall(VecCreate(PETSC_COMM_WORLD,&yi)); 1006 PetscCall(VecSetSizes(yi,PETSC_DECIDE,user->mx*user->mx*user->mx)); 1007 PetscCall(VecSetFromOptions(yi)); 1008 PetscCall(VecDuplicateVecs(yi,user->nt,&user->yi)); 1009 PetscCall(VecDuplicateVecs(yi,user->nt,&user->yiwork)); 1010 1011 PetscCall(VecGetOwnershipRange(user->y,&lo2,&hi2)); 1012 istart = 0; 1013 for (i=0; i<user->nt; i++) { 1014 PetscCall(VecGetOwnershipRange(user->yi[i],&lo,&hi)); 1015 PetscCall(ISCreateStride(PETSC_COMM_SELF,hi-lo,lo,1,&is_to_yi)); 1016 PetscCall(ISCreateStride(PETSC_COMM_SELF,hi-lo,lo2+istart,1,&is_from_y)); 1017 PetscCall(VecScatterCreate(user->y,is_from_y,user->yi[i],is_to_yi,&user->yi_scatter[i])); 1018 istart = istart + hi-lo; 1019 PetscCall(ISDestroy(&is_to_yi)); 1020 PetscCall(ISDestroy(&is_from_y)); 1021 } 1022 PetscCall(VecDestroy(&yi)); 1023 1024 /* Create scatter from d to d_1,d_2,...,d_ns */ 1025 PetscCall(PetscMalloc1(user->ns*user->ndata,&user->di_scatter)); 1026 PetscCall(VecCreate(PETSC_COMM_WORLD,&di)); 1027 PetscCall(VecSetSizes(di,PETSC_DECIDE,user->ndata)); 1028 PetscCall(VecSetFromOptions(di)); 1029 PetscCall(VecDuplicateVecs(di,user->ns,&user->di)); 1030 PetscCall(VecGetOwnershipRange(user->d,&lo2,&hi2)); 1031 istart = 0; 1032 for (i=0; i<user->ns; i++) { 1033 PetscCall(VecGetOwnershipRange(user->di[i],&lo,&hi)); 1034 PetscCall(ISCreateStride(PETSC_COMM_SELF,hi-lo,lo,1,&is_to_di)); 1035 PetscCall(ISCreateStride(PETSC_COMM_SELF,hi-lo,lo2+istart,1,&is_from_d)); 1036 PetscCall(VecScatterCreate(user->d,is_from_d,user->di[i],is_to_di,&user->di_scatter[i])); 1037 istart = istart + hi-lo; 1038 PetscCall(ISDestroy(&is_to_di)); 1039 PetscCall(ISDestroy(&is_from_d)); 1040 } 1041 PetscCall(VecDestroy(&di)); 1042 1043 /* Assemble RHS of forward problem */ 1044 PetscCall(VecCopy(bc,user->yiwork[0])); 1045 for (i=1; i<user->nt; i++) { 1046 PetscCall(VecSet(user->yiwork[i],0.0)); 1047 } 1048 PetscCall(Gather_i(user->q,user->yiwork,user->yi_scatter,user->nt)); 1049 PetscCall(VecDestroy(&bc)); 1050 1051 /* Compute true state function yt given ut */ 1052 PetscCall(VecCreate(PETSC_COMM_WORLD,&user->ytrue)); 1053 PetscCall(VecSetSizes(user->ytrue,PETSC_DECIDE,n*user->nt)); 1054 PetscCall(VecSetFromOptions(user->ytrue)); 1055 1056 /* First compute Av_u = Av*exp(-u) */ 1057 PetscCall(VecSet(user->uwork,0)); 1058 PetscCall(VecAXPY(user->uwork,-1.0,user->utrue)); /* Note: user->utrue */ 1059 PetscCall(VecExp(user->uwork)); 1060 PetscCall(MatMult(user->Av,user->uwork,user->Av_u)); 1061 1062 /* Symbolic DSG = Div * Grad */ 1063 PetscCall(MatProductCreate(user->Div,user->Grad,NULL,&user->DSG)); 1064 PetscCall(MatProductSetType(user->DSG,MATPRODUCT_AB)); 1065 PetscCall(MatProductSetAlgorithm(user->DSG,"default")); 1066 PetscCall(MatProductSetFill(user->DSG,PETSC_DEFAULT)); 1067 PetscCall(MatProductSetFromOptions(user->DSG)); 1068 PetscCall(MatProductSymbolic(user->DSG)); 1069 1070 user->dsg_formed = PETSC_TRUE; 1071 1072 /* Next form DSG = Div*Grad */ 1073 PetscCall(MatCopy(user->Div,user->Divwork,SAME_NONZERO_PATTERN)); 1074 PetscCall(MatDiagonalScale(user->Divwork,NULL,user->Av_u)); 1075 if (user->dsg_formed) { 1076 PetscCall(MatProductNumeric(user->DSG)); 1077 } else { 1078 PetscCall(MatMatMult(user->Div,user->Grad,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&user->DSG)); 1079 user->dsg_formed = PETSC_TRUE; 1080 } 1081 /* B = speye(nx^3) + ht*DSG; */ 1082 PetscCall(MatScale(user->DSG,user->ht)); 1083 PetscCall(MatShift(user->DSG,1.0)); 1084 1085 /* Now solve for ytrue */ 1086 PetscCall(StateMatInvMult(user->Js,user->q,user->ytrue)); 1087 1088 /* Initial guess y0 for state given u0 */ 1089 1090 /* First compute Av_u = Av*exp(-u) */ 1091 PetscCall(VecSet(user->uwork,0)); 1092 PetscCall(VecAXPY(user->uwork,-1.0,user->u)); /* Note: user->u */ 1093 PetscCall(VecExp(user->uwork)); 1094 PetscCall(MatMult(user->Av,user->uwork,user->Av_u)); 1095 1096 /* Next form DSG = Div*S*Grad */ 1097 PetscCall(MatCopy(user->Div,user->Divwork,SAME_NONZERO_PATTERN)); 1098 PetscCall(MatDiagonalScale(user->Divwork,NULL,user->Av_u)); 1099 if (user->dsg_formed) { 1100 PetscCall(MatProductNumeric(user->DSG)); 1101 } else { 1102 PetscCall(MatMatMult(user->Div,user->Grad,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&user->DSG)); 1103 1104 user->dsg_formed = PETSC_TRUE; 1105 } 1106 /* B = speye(nx^3) + ht*DSG; */ 1107 PetscCall(MatScale(user->DSG,user->ht)); 1108 PetscCall(MatShift(user->DSG,1.0)); 1109 1110 /* Now solve for y */ 1111 PetscCall(StateMatInvMult(user->Js,user->q,user->y)); 1112 1113 /* Construct projection matrix Q, a block diagonal matrix consisting of nt copies of Qblock along the diagonal */ 1114 PetscCall(MatCreate(PETSC_COMM_WORLD,&user->Qblock)); 1115 PetscCall(MatSetSizes(user->Qblock,PETSC_DECIDE,PETSC_DECIDE,user->ndata,n)); 1116 PetscCall(MatSetFromOptions(user->Qblock)); 1117 PetscCall(MatMPIAIJSetPreallocation(user->Qblock,8,NULL,8,NULL)); 1118 PetscCall(MatSeqAIJSetPreallocation(user->Qblock,8,NULL)); 1119 1120 for (i=0; i<user->mx; i++) { 1121 x[i] = h*(i+0.5); 1122 y[i] = h*(i+0.5); 1123 z[i] = h*(i+0.5); 1124 } 1125 1126 PetscCall(MatGetOwnershipRange(user->Qblock,&istart,&iend)); 1127 nx = user->mx; ny = user->mx; nz = user->mx; 1128 for (i=istart; i<iend; i++) { 1129 xri = xr[i]; 1130 im = 0; 1131 xim = x[im]; 1132 while (xri>xim && im<nx) { 1133 im = im+1; 1134 xim = x[im]; 1135 } 1136 indx1 = im-1; 1137 indx2 = im; 1138 dx1 = xri - x[indx1]; 1139 dx2 = x[indx2] - xri; 1140 1141 yri = yr[i]; 1142 im = 0; 1143 yim = y[im]; 1144 while (yri>yim && im<ny) { 1145 im = im+1; 1146 yim = y[im]; 1147 } 1148 indy1 = im-1; 1149 indy2 = im; 1150 dy1 = yri - y[indy1]; 1151 dy2 = y[indy2] - yri; 1152 1153 zri = zr[i]; 1154 im = 0; 1155 zim = z[im]; 1156 while (zri>zim && im<nz) { 1157 im = im+1; 1158 zim = z[im]; 1159 } 1160 indz1 = im-1; 1161 indz2 = im; 1162 dz1 = zri - z[indz1]; 1163 dz2 = z[indz2] - zri; 1164 1165 Dx = x[indx2] - x[indx1]; 1166 Dy = y[indy2] - y[indy1]; 1167 Dz = z[indz2] - z[indz1]; 1168 1169 j = indx1 + indy1*nx + indz1*nx*ny; 1170 v = (1-dx1/Dx)*(1-dy1/Dy)*(1-dz1/Dz); 1171 PetscCall(MatSetValues(user->Qblock,1,&i,1,&j,&v,INSERT_VALUES)); 1172 1173 j = indx1 + indy1*nx + indz2*nx*ny; 1174 v = (1-dx1/Dx)*(1-dy1/Dy)*(1-dz2/Dz); 1175 PetscCall(MatSetValues(user->Qblock,1,&i,1,&j,&v,INSERT_VALUES)); 1176 1177 j = indx1 + indy2*nx + indz1*nx*ny; 1178 v = (1-dx1/Dx)*(1-dy2/Dy)*(1-dz1/Dz); 1179 PetscCall(MatSetValues(user->Qblock,1,&i,1,&j,&v,INSERT_VALUES)); 1180 1181 j = indx1 + indy2*nx + indz2*nx*ny; 1182 v = (1-dx1/Dx)*(1-dy2/Dy)*(1-dz2/Dz); 1183 PetscCall(MatSetValues(user->Qblock,1,&i,1,&j,&v,INSERT_VALUES)); 1184 1185 j = indx2 + indy1*nx + indz1*nx*ny; 1186 v = (1-dx2/Dx)*(1-dy1/Dy)*(1-dz1/Dz); 1187 PetscCall(MatSetValues(user->Qblock,1,&i,1,&j,&v,INSERT_VALUES)); 1188 1189 j = indx2 + indy1*nx + indz2*nx*ny; 1190 v = (1-dx2/Dx)*(1-dy1/Dy)*(1-dz2/Dz); 1191 PetscCall(MatSetValues(user->Qblock,1,&i,1,&j,&v,INSERT_VALUES)); 1192 1193 j = indx2 + indy2*nx + indz1*nx*ny; 1194 v = (1-dx2/Dx)*(1-dy2/Dy)*(1-dz1/Dz); 1195 PetscCall(MatSetValues(user->Qblock,1,&i,1,&j,&v,INSERT_VALUES)); 1196 1197 j = indx2 + indy2*nx + indz2*nx*ny; 1198 v = (1-dx2/Dx)*(1-dy2/Dy)*(1-dz2/Dz); 1199 PetscCall(MatSetValues(user->Qblock,1,&i,1,&j,&v,INSERT_VALUES)); 1200 } 1201 PetscCall(MatAssemblyBegin(user->Qblock,MAT_FINAL_ASSEMBLY)); 1202 PetscCall(MatAssemblyEnd(user->Qblock,MAT_FINAL_ASSEMBLY)); 1203 1204 PetscCall(MatTranspose(user->Qblock,MAT_INITIAL_MATRIX,&user->QblockT)); 1205 PetscCall(MatTranspose(user->L,MAT_INITIAL_MATRIX,&user->LT)); 1206 1207 /* Add noise to the measurement data */ 1208 PetscCall(VecSet(user->ywork,1.0)); 1209 PetscCall(VecAYPX(user->ywork,user->noise,user->ytrue)); 1210 PetscCall(Scatter_i(user->ywork,user->yiwork,user->yi_scatter,user->nt)); 1211 for (j=0; j<user->ns; j++) { 1212 i = user->sample_times[j]; 1213 PetscCall(MatMult(user->Qblock,user->yiwork[i],user->di[j])); 1214 } 1215 PetscCall(Gather_i(user->d,user->di,user->di_scatter,user->ns)); 1216 1217 /* Now that initial conditions have been set, let the user pass tolerance options to the KSP solver */ 1218 PetscCall(KSPSetFromOptions(user->solver)); 1219 PetscCall(PetscFree(x)); 1220 PetscCall(PetscFree(y)); 1221 PetscCall(PetscFree(z)); 1222 PetscFunctionReturn(0); 1223 } 1224 1225 PetscErrorCode ParabolicDestroy(AppCtx *user) 1226 { 1227 PetscInt i; 1228 1229 PetscFunctionBegin; 1230 PetscCall(MatDestroy(&user->Qblock)); 1231 PetscCall(MatDestroy(&user->QblockT)); 1232 PetscCall(MatDestroy(&user->Div)); 1233 PetscCall(MatDestroy(&user->Divwork)); 1234 PetscCall(MatDestroy(&user->Grad)); 1235 PetscCall(MatDestroy(&user->Av)); 1236 PetscCall(MatDestroy(&user->Avwork)); 1237 PetscCall(MatDestroy(&user->AvT)); 1238 PetscCall(MatDestroy(&user->DSG)); 1239 PetscCall(MatDestroy(&user->L)); 1240 PetscCall(MatDestroy(&user->LT)); 1241 PetscCall(KSPDestroy(&user->solver)); 1242 PetscCall(MatDestroy(&user->Js)); 1243 PetscCall(MatDestroy(&user->Jd)); 1244 PetscCall(MatDestroy(&user->JsInv)); 1245 PetscCall(MatDestroy(&user->JsBlock)); 1246 PetscCall(MatDestroy(&user->JsBlockPrec)); 1247 PetscCall(VecDestroy(&user->u)); 1248 PetscCall(VecDestroy(&user->uwork)); 1249 PetscCall(VecDestroy(&user->utrue)); 1250 PetscCall(VecDestroy(&user->y)); 1251 PetscCall(VecDestroy(&user->ywork)); 1252 PetscCall(VecDestroy(&user->ytrue)); 1253 PetscCall(VecDestroyVecs(user->nt,&user->yi)); 1254 PetscCall(VecDestroyVecs(user->nt,&user->yiwork)); 1255 PetscCall(VecDestroyVecs(user->ns,&user->di)); 1256 PetscCall(PetscFree(user->yi)); 1257 PetscCall(PetscFree(user->yiwork)); 1258 PetscCall(PetscFree(user->di)); 1259 PetscCall(VecDestroy(&user->c)); 1260 PetscCall(VecDestroy(&user->cwork)); 1261 PetscCall(VecDestroy(&user->ur)); 1262 PetscCall(VecDestroy(&user->q)); 1263 PetscCall(VecDestroy(&user->d)); 1264 PetscCall(VecDestroy(&user->dwork)); 1265 PetscCall(VecDestroy(&user->lwork)); 1266 PetscCall(VecDestroy(&user->S)); 1267 PetscCall(VecDestroy(&user->Swork)); 1268 PetscCall(VecDestroy(&user->Av_u)); 1269 PetscCall(VecDestroy(&user->Twork)); 1270 PetscCall(VecDestroy(&user->Rwork)); 1271 PetscCall(VecDestroy(&user->js_diag)); 1272 PetscCall(ISDestroy(&user->s_is)); 1273 PetscCall(ISDestroy(&user->d_is)); 1274 PetscCall(VecScatterDestroy(&user->state_scatter)); 1275 PetscCall(VecScatterDestroy(&user->design_scatter)); 1276 for (i=0; i<user->nt; i++) { 1277 PetscCall(VecScatterDestroy(&user->yi_scatter[i])); 1278 } 1279 for (i=0; i<user->ns; i++) { 1280 PetscCall(VecScatterDestroy(&user->di_scatter[i])); 1281 } 1282 PetscCall(PetscFree(user->yi_scatter)); 1283 PetscCall(PetscFree(user->di_scatter)); 1284 PetscCall(PetscFree(user->sample_times)); 1285 PetscFunctionReturn(0); 1286 } 1287 1288 PetscErrorCode ParabolicMonitor(Tao tao, void *ptr) 1289 { 1290 Vec X; 1291 PetscReal unorm,ynorm; 1292 AppCtx *user = (AppCtx*)ptr; 1293 1294 PetscFunctionBegin; 1295 PetscCall(TaoGetSolution(tao,&X)); 1296 PetscCall(Scatter(X,user->ywork,user->state_scatter,user->uwork,user->design_scatter)); 1297 PetscCall(VecAXPY(user->ywork,-1.0,user->ytrue)); 1298 PetscCall(VecAXPY(user->uwork,-1.0,user->utrue)); 1299 PetscCall(VecNorm(user->uwork,NORM_2,&unorm)); 1300 PetscCall(VecNorm(user->ywork,NORM_2,&ynorm)); 1301 PetscCall(PetscPrintf(MPI_COMM_WORLD, "||u-ut||=%g ||y-yt||=%g\n",(double)unorm,(double)ynorm)); 1302 PetscFunctionReturn(0); 1303 } 1304 1305 /*TEST 1306 1307 build: 1308 requires: !complex 1309 1310 test: 1311 args: -tao_cmonitor -tao_type lcl -ns 1 -tao_gatol 1.e-4 -ksp_max_it 30 1312 requires: !single 1313 1314 TEST*/ 1315