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 ierr = PetscInitialize(&argc, &argv, (char*)0,help);if (ierr) return ierr; 128 user.mx = 8; 129 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"parabolic example",NULL);CHKERRQ(ierr); 130 CHKERRQ(PetscOptionsInt("-mx","Number of grid points in each direction","",user.mx,&user.mx,NULL)); 131 user.nt = 8; 132 CHKERRQ(PetscOptionsInt("-nt","Number of time steps","",user.nt,&user.nt,NULL)); 133 user.ndata = 64; 134 CHKERRQ(PetscOptionsInt("-ndata","Numbers of data points per sample","",user.ndata,&user.ndata,NULL)); 135 user.ns = 8; 136 CHKERRQ(PetscOptionsInt("-ns","Number of samples","",user.ns,&user.ns,NULL)); 137 user.alpha = 1.0; 138 CHKERRQ(PetscOptionsReal("-alpha","Regularization parameter","",user.alpha,&user.alpha,NULL)); 139 user.beta = 0.01; 140 CHKERRQ(PetscOptionsReal("-beta","Weight attributed to ||u||^2 in regularization functional","",user.beta,&user.beta,NULL)); 141 user.noise = 0.01; 142 CHKERRQ(PetscOptionsReal("-noise","Amount of noise to add to data","",user.noise,&user.noise,NULL)); 143 CHKERRQ(PetscOptionsInt("-ntests","Number of times to repeat TaoSolve","",ntests,&ntests,NULL)); 144 ierr = PetscOptionsEnd();CHKERRQ(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 CHKERRQ(VecCreate(PETSC_COMM_WORLD,&user.u)); 151 CHKERRQ(VecCreate(PETSC_COMM_WORLD,&user.y)); 152 CHKERRQ(VecCreate(PETSC_COMM_WORLD,&user.c)); 153 CHKERRQ(VecSetSizes(user.u,PETSC_DECIDE,user.n-user.m*user.nt)); 154 CHKERRQ(VecSetSizes(user.y,PETSC_DECIDE,user.m*user.nt)); 155 CHKERRQ(VecSetSizes(user.c,PETSC_DECIDE,user.m*user.nt)); 156 CHKERRQ(VecSetFromOptions(user.u)); 157 CHKERRQ(VecSetFromOptions(user.y)); 158 CHKERRQ(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 CHKERRQ(VecCreate(PETSC_COMM_WORLD,&x)); 169 170 CHKERRQ(VecGetOwnershipRange(user.y,&lo,&hi)); 171 CHKERRQ(VecGetOwnershipRange(user.u,&lo2,&hi2)); 172 173 CHKERRQ(ISCreateStride(PETSC_COMM_SELF,hi-lo,lo,1,&is_allstate)); 174 CHKERRQ(ISCreateStride(PETSC_COMM_SELF,hi-lo,lo+lo2,1,&user.s_is)); 175 176 CHKERRQ(ISCreateStride(PETSC_COMM_SELF,hi2-lo2,lo2,1,&is_alldesign)); 177 CHKERRQ(ISCreateStride(PETSC_COMM_SELF,hi2-lo2,hi+lo2,1,&user.d_is)); 178 179 CHKERRQ(VecSetSizes(x,hi-lo+hi2-lo2,user.n)); 180 CHKERRQ(VecSetFromOptions(x)); 181 182 CHKERRQ(VecScatterCreate(x,user.s_is,user.y,is_allstate,&user.state_scatter)); 183 CHKERRQ(VecScatterCreate(x,user.d_is,user.u,is_alldesign,&user.design_scatter)); 184 CHKERRQ(ISDestroy(&is_alldesign)); 185 CHKERRQ(ISDestroy(&is_allstate)); 186 187 /* Create TAO solver and set desired solution method */ 188 CHKERRQ(TaoCreate(PETSC_COMM_WORLD,&tao)); 189 CHKERRQ(TaoSetType(tao,TAOLCL)); 190 191 /* Set up initial vectors and matrices */ 192 CHKERRQ(ParabolicInitialize(&user)); 193 194 CHKERRQ(Gather(x,user.y,user.state_scatter,user.u,user.design_scatter)); 195 CHKERRQ(VecDuplicate(x,&x0)); 196 CHKERRQ(VecCopy(x,x0)); 197 198 /* Set solution vector with an initial guess */ 199 CHKERRQ(TaoSetSolution(tao,x)); 200 CHKERRQ(TaoSetObjective(tao, FormFunction, &user)); 201 CHKERRQ(TaoSetGradient(tao, NULL, FormGradient, &user)); 202 CHKERRQ(TaoSetConstraintsRoutine(tao, user.c, FormConstraints, &user)); 203 204 CHKERRQ(TaoSetJacobianStateRoutine(tao, user.Js, user.JsBlockPrec, user.JsInv, FormJacobianState, &user)); 205 CHKERRQ(TaoSetJacobianDesignRoutine(tao, user.Jd, FormJacobianDesign, &user)); 206 207 CHKERRQ(TaoSetFromOptions(tao)); 208 CHKERRQ(TaoSetStateDesignIS(tao,user.s_is,user.d_is)); 209 210 /* SOLVE THE APPLICATION */ 211 CHKERRQ(PetscLogStageRegister("Trials",&stages[0])); 212 CHKERRQ(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 CHKERRQ(TaoSolve(tao)); 217 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"KSP Iterations = %D\n",user.ksp_its-ksp_old)); 218 CHKERRQ(VecCopy(x0,x)); 219 CHKERRQ(TaoSetSolution(tao,x)); 220 } 221 CHKERRQ(PetscLogStagePop()); 222 CHKERRQ(PetscBarrier((PetscObject)x)); 223 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"KSP iterations within initialization: ")); 224 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"%D\n",user.ksp_its_initial)); 225 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Total KSP iterations over %D trial(s): ",ntests)); 226 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"%D\n",user.ksp_its)); 227 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"KSP iterations per trial: ")); 228 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"%D\n",(user.ksp_its-user.ksp_its_initial)/ntests)); 229 230 CHKERRQ(TaoDestroy(&tao)); 231 CHKERRQ(VecDestroy(&x)); 232 CHKERRQ(VecDestroy(&x0)); 233 CHKERRQ(ParabolicDestroy(&user)); 234 235 ierr = PetscFinalize(); 236 return ierr; 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 CHKERRQ(Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter)); 252 CHKERRQ(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 CHKERRQ(MatMult(user->Qblock,user->yi[i],user->di[j])); 256 } 257 CHKERRQ(Gather_i(user->dwork,user->di,user->di_scatter,user->ns)); 258 CHKERRQ(VecAXPY(user->dwork,-1.0,user->d)); 259 CHKERRQ(VecDot(user->dwork,user->dwork,&d1)); 260 261 CHKERRQ(VecWAXPY(user->uwork,-1.0,user->ur,user->u)); 262 CHKERRQ(MatMult(user->L,user->uwork,user->lwork)); 263 CHKERRQ(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 CHKERRQ(Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter)); 281 CHKERRQ(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 CHKERRQ(MatMult(user->Qblock,user->yi[i],user->di[j])); 285 } 286 CHKERRQ(Gather_i(user->dwork,user->di,user->di_scatter,user->ns)); 287 CHKERRQ(VecAXPY(user->dwork,-1.0,user->d)); 288 CHKERRQ(Scatter_i(user->dwork,user->di,user->di_scatter,user->ns)); 289 CHKERRQ(VecSet(user->ywork,0.0)); 290 CHKERRQ(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 CHKERRQ(MatMult(user->QblockT,user->di[j],user->yiwork[i])); 294 } 295 CHKERRQ(Gather_i(user->ywork,user->yiwork,user->yi_scatter,user->nt)); 296 297 CHKERRQ(VecWAXPY(user->uwork,-1.0,user->ur,user->u)); 298 CHKERRQ(MatMult(user->L,user->uwork,user->lwork)); 299 CHKERRQ(MatMult(user->LT,user->lwork,user->uwork)); 300 CHKERRQ(VecScale(user->uwork, user->alpha)); 301 CHKERRQ(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 CHKERRQ(Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter)); 313 CHKERRQ(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 CHKERRQ(MatMult(user->Qblock,user->yi[i],user->di[j])); 317 } 318 CHKERRQ(Gather_i(user->dwork,user->di,user->di_scatter,user->ns)); 319 CHKERRQ(VecAXPY(user->dwork,-1.0,user->d)); 320 CHKERRQ(VecDot(user->dwork,user->dwork,&d1)); 321 CHKERRQ(Scatter_i(user->dwork,user->di,user->di_scatter,user->ns)); 322 CHKERRQ(VecSet(user->ywork,0.0)); 323 CHKERRQ(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 CHKERRQ(MatMult(user->QblockT,user->di[j],user->yiwork[i])); 327 } 328 CHKERRQ(Gather_i(user->ywork,user->yiwork,user->yi_scatter,user->nt)); 329 330 CHKERRQ(VecWAXPY(user->uwork,-1.0,user->ur,user->u)); 331 CHKERRQ(MatMult(user->L,user->uwork,user->lwork)); 332 CHKERRQ(VecDot(user->lwork,user->lwork,&d2)); 333 CHKERRQ(MatMult(user->LT,user->lwork,user->uwork)); 334 CHKERRQ(VecScale(user->uwork, user->alpha)); 335 *f = 0.5 * (d1 + user->alpha*d2); 336 337 CHKERRQ(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 CHKERRQ(Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter)); 351 CHKERRQ(VecSet(user->uwork,0)); 352 CHKERRQ(VecAXPY(user->uwork,-1.0,user->u)); 353 CHKERRQ(VecExp(user->uwork)); 354 CHKERRQ(MatMult(user->Av,user->uwork,user->Av_u)); 355 CHKERRQ(VecCopy(user->Av_u,user->Swork)); 356 CHKERRQ(VecReciprocal(user->Swork)); 357 CHKERRQ(MatCopy(user->Div,user->Divwork,SAME_NONZERO_PATTERN)); 358 CHKERRQ(MatDiagonalScale(user->Divwork,NULL,user->Swork)); 359 if (user->dsg_formed) { 360 CHKERRQ(MatProductNumeric(user->DSG)); 361 } else { 362 CHKERRQ(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 CHKERRQ(MatScale(user->DSG,user->ht)); 368 CHKERRQ(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 CHKERRQ(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 CHKERRQ(MatShellGetContext(J_shell,&user)); 390 CHKERRQ(Scatter_i(X,user->yi,user->yi_scatter,user->nt)); 391 CHKERRQ(MatMult(user->JsBlock,user->yi[0],user->yiwork[0])); 392 for (i=1; i<user->nt; i++) { 393 CHKERRQ(MatMult(user->JsBlock,user->yi[i],user->yiwork[i])); 394 CHKERRQ(VecAXPY(user->yiwork[i],-1.0,user->yi[i-1])); 395 } 396 CHKERRQ(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 CHKERRQ(MatShellGetContext(J_shell,&user)); 407 CHKERRQ(Scatter_i(X,user->yi,user->yi_scatter,user->nt)); 408 for (i=0; i<user->nt-1; i++) { 409 CHKERRQ(MatMult(user->JsBlock,user->yi[i],user->yiwork[i])); 410 CHKERRQ(VecAXPY(user->yiwork[i],-1.0,user->yi[i+1])); 411 } 412 i = user->nt-1; 413 CHKERRQ(MatMult(user->JsBlock,user->yi[i],user->yiwork[i])); 414 CHKERRQ(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 CHKERRQ(MatShellGetContext(J_shell,&user)); 424 CHKERRQ(MatMult(user->Grad,X,user->Swork)); 425 CHKERRQ(VecPointwiseDivide(user->Swork,user->Swork,user->Av_u)); 426 CHKERRQ(MatMult(user->Div,user->Swork,Y)); 427 CHKERRQ(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 CHKERRQ(MatShellGetContext(J_shell,&user)); 438 439 /* sdiag(1./v) */ 440 CHKERRQ(VecSet(user->uwork,0)); 441 CHKERRQ(VecAXPY(user->uwork,-1.0,user->u)); 442 CHKERRQ(VecExp(user->uwork)); 443 444 /* sdiag(1./((Av*(1./v)).^2)) */ 445 CHKERRQ(MatMult(user->Av,user->uwork,user->Swork)); 446 CHKERRQ(VecPointwiseMult(user->Swork,user->Swork,user->Swork)); 447 CHKERRQ(VecReciprocal(user->Swork)); 448 449 /* (Av * (sdiag(1./v) * b)) */ 450 CHKERRQ(VecPointwiseMult(user->uwork,user->uwork,X)); 451 CHKERRQ(MatMult(user->Av,user->uwork,user->Twork)); 452 453 /* (sdiag(1./((Av*(1./v)).^2)) * (Av * (sdiag(1./v) * b))) */ 454 CHKERRQ(VecPointwiseMult(user->Swork,user->Twork,user->Swork)); 455 456 CHKERRQ(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 CHKERRQ(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 CHKERRQ(VecPointwiseMult(user->Twork,user->Twork,user->Swork)); 463 CHKERRQ(MatMult(user->Div,user->Twork,user->yiwork[i])); 464 CHKERRQ(VecScale(user->yiwork[i],user->ht)); 465 } 466 CHKERRQ(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 CHKERRQ(MatShellGetContext(J_shell,&user)); 478 479 /* sdiag(1./((Av*(1./v)).^2)) */ 480 CHKERRQ(VecSet(user->uwork,0)); 481 CHKERRQ(VecAXPY(user->uwork,-1.0,user->u)); 482 CHKERRQ(VecExp(user->uwork)); 483 CHKERRQ(MatMult(user->Av,user->uwork,user->Rwork)); 484 CHKERRQ(VecPointwiseMult(user->Rwork,user->Rwork,user->Rwork)); 485 CHKERRQ(VecReciprocal(user->Rwork)); 486 487 CHKERRQ(VecSet(Y,0.0)); 488 CHKERRQ(Scatter_i(user->y,user->yi,user->yi_scatter,user->nt)); 489 CHKERRQ(Scatter_i(X,user->yiwork,user->yi_scatter,user->nt)); 490 for (i=0; i<user->nt; i++) { 491 /* (Div' * b(:,i)) */ 492 CHKERRQ(MatMult(user->Grad,user->yiwork[i],user->Swork)); 493 494 /* sdiag(Grad*y(:,i)) */ 495 CHKERRQ(MatMult(user->Grad,user->yi[i],user->Twork)); 496 497 /* (sdiag(Grad*y(:,i)) * (Div' * b(:,i))) */ 498 CHKERRQ(VecPointwiseMult(user->Twork,user->Swork,user->Twork)); 499 500 /* (sdiag(1./((Av*(1./v)).^2)) * (sdiag(Grad*y(:,i)) * (Div' * b(:,i)))) */ 501 CHKERRQ(VecPointwiseMult(user->Twork,user->Rwork,user->Twork)); 502 503 /* (Av' * (sdiag(1./((Av*(1./v)).^2)) * (sdiag(Grad*y(:,i)) * (Div' * b(:,i))))) */ 504 CHKERRQ(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 CHKERRQ(VecPointwiseMult(user->yiwork[i],user->uwork,user->yiwork[i])); 508 CHKERRQ(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 CHKERRQ(PCShellGetContext(PC_shell,&user)); 519 520 if (user->dsg_formed) { 521 CHKERRQ(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 CHKERRQ(MatShellGetContext(J_shell,&user)); 533 534 if (Y == user->ytrue) { 535 /* First solve is done with true solution to set up problem */ 536 CHKERRQ(KSPSetTolerances(user->solver,1e-8,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT)); 537 } else { 538 CHKERRQ(KSPSetTolerances(user->solver,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT)); 539 } 540 541 CHKERRQ(Scatter_i(X,user->yi,user->yi_scatter,user->nt)); 542 CHKERRQ(KSPSolve(user->solver,user->yi[0],user->yiwork[0])); 543 CHKERRQ(KSPGetIterationNumber(user->solver,&its)); 544 user->ksp_its = user->ksp_its + its; 545 546 for (i=1; i<user->nt; i++) { 547 CHKERRQ(VecAXPY(user->yi[i],1.0,user->yiwork[i-1])); 548 CHKERRQ(KSPSolve(user->solver,user->yi[i],user->yiwork[i])); 549 CHKERRQ(KSPGetIterationNumber(user->solver,&its)); 550 user->ksp_its = user->ksp_its + its; 551 } 552 CHKERRQ(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 CHKERRQ(MatShellGetContext(J_shell,&user)); 563 564 CHKERRQ(Scatter_i(X,user->yi,user->yi_scatter,user->nt)); 565 566 i = user->nt - 1; 567 CHKERRQ(KSPSolve(user->solver,user->yi[i],user->yiwork[i])); 568 569 CHKERRQ(KSPGetIterationNumber(user->solver,&its)); 570 user->ksp_its = user->ksp_its + its; 571 572 for (i=user->nt-2; i>=0; i--) { 573 CHKERRQ(VecAXPY(user->yi[i],1.0,user->yiwork[i+1])); 574 CHKERRQ(KSPSolve(user->solver,user->yi[i],user->yiwork[i])); 575 576 CHKERRQ(KSPGetIterationNumber(user->solver,&its)); 577 user->ksp_its = user->ksp_its + its; 578 579 } 580 581 CHKERRQ(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 CHKERRQ(MatShellGetContext(J_shell,&user)); 591 592 CHKERRQ(MatCreateShell(PETSC_COMM_WORLD,PETSC_DETERMINE,PETSC_DETERMINE,user->m,user->m,user,new_shell)); 593 CHKERRQ(MatShellSetOperation(*new_shell,MATOP_MULT,(void(*)(void))StateMatMult)); 594 CHKERRQ(MatShellSetOperation(*new_shell,MATOP_DUPLICATE,(void(*)(void))StateMatDuplicate)); 595 CHKERRQ(MatShellSetOperation(*new_shell,MATOP_MULT_TRANSPOSE,(void(*)(void))StateMatMultTranspose)); 596 CHKERRQ(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 CHKERRQ(MatShellGetContext(J_shell,&user)); 606 CHKERRQ(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 CHKERRQ(Scatter(X,user->y,user->state_scatter,user->u,user->design_scatter)); 624 CHKERRQ(Scatter_i(user->y,user->yi,user->yi_scatter,user->nt)); 625 CHKERRQ(MatMult(user->JsBlock,user->yi[0],user->yiwork[0])); 626 for (i=1; i<user->nt; i++) { 627 CHKERRQ(MatMult(user->JsBlock,user->yi[i],user->yiwork[i])); 628 CHKERRQ(VecAXPY(user->yiwork[i],-1.0,user->yi[i-1])); 629 } 630 CHKERRQ(Gather_i(C,user->yiwork,user->yi_scatter,user->nt)); 631 CHKERRQ(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 CHKERRQ(VecScatterBegin(s_scat,x,state,INSERT_VALUES,SCATTER_FORWARD)); 639 CHKERRQ(VecScatterEnd(s_scat,x,state,INSERT_VALUES,SCATTER_FORWARD)); 640 CHKERRQ(VecScatterBegin(d_scat,x,design,INSERT_VALUES,SCATTER_FORWARD)); 641 CHKERRQ(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 CHKERRQ(VecScatterBegin(scat[i],y,yi[i],INSERT_VALUES,SCATTER_FORWARD)); 652 CHKERRQ(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 CHKERRQ(VecScatterBegin(s_scat,state,x,INSERT_VALUES,SCATTER_REVERSE)); 661 CHKERRQ(VecScatterEnd(s_scat,state,x,INSERT_VALUES,SCATTER_REVERSE)); 662 CHKERRQ(VecScatterBegin(d_scat,design,x,INSERT_VALUES,SCATTER_REVERSE)); 663 CHKERRQ(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 CHKERRQ(VecScatterBegin(scat[i],yi[i],y,INSERT_VALUES,SCATTER_REVERSE)); 674 CHKERRQ(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 CHKERRQ(PetscMalloc1(user->mx,&x)); 720 CHKERRQ(PetscMalloc1(user->mx,&y)); 721 CHKERRQ(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 CHKERRQ(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 CHKERRQ(VecCreate(PETSC_COMM_WORLD,&XX)); 739 CHKERRQ(VecCreate(PETSC_COMM_WORLD,&user->q)); 740 CHKERRQ(VecSetSizes(XX,PETSC_DECIDE,n)); 741 CHKERRQ(VecSetSizes(user->q,PETSC_DECIDE,n*user->nt)); 742 CHKERRQ(VecSetFromOptions(XX)); 743 CHKERRQ(VecSetFromOptions(user->q)); 744 745 CHKERRQ(VecDuplicate(XX,&YY)); 746 CHKERRQ(VecDuplicate(XX,&ZZ)); 747 CHKERRQ(VecDuplicate(XX,&XXwork)); 748 CHKERRQ(VecDuplicate(XX,&YYwork)); 749 CHKERRQ(VecDuplicate(XX,&ZZwork)); 750 CHKERRQ(VecDuplicate(XX,&UTwork)); 751 CHKERRQ(VecDuplicate(XX,&user->utrue)); 752 CHKERRQ(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 CHKERRQ(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 CHKERRQ(VecSetValues(XX,1,&linear_index,&vx,INSERT_VALUES)); 768 CHKERRQ(VecSetValues(YY,1,&linear_index,&vy,INSERT_VALUES)); 769 CHKERRQ(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 CHKERRQ(VecSetValues(bc,1,&linear_index,&v,INSERT_VALUES)); 773 } 774 } 775 776 CHKERRQ(VecAssemblyBegin(XX)); 777 CHKERRQ(VecAssemblyEnd(XX)); 778 CHKERRQ(VecAssemblyBegin(YY)); 779 CHKERRQ(VecAssemblyEnd(YY)); 780 CHKERRQ(VecAssemblyBegin(ZZ)); 781 CHKERRQ(VecAssemblyEnd(ZZ)); 782 CHKERRQ(VecAssemblyBegin(bc)); 783 CHKERRQ(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 CHKERRQ(VecCopy(XX,XXwork)); 788 CHKERRQ(VecCopy(YY,YYwork)); 789 CHKERRQ(VecCopy(ZZ,ZZwork)); 790 791 CHKERRQ(VecShift(XXwork,-0.5)); 792 CHKERRQ(VecShift(YYwork,-0.5)); 793 CHKERRQ(VecShift(ZZwork,-0.5)); 794 795 CHKERRQ(VecPointwiseMult(XXwork,XXwork,XXwork)); 796 CHKERRQ(VecPointwiseMult(YYwork,YYwork,YYwork)); 797 CHKERRQ(VecPointwiseMult(ZZwork,ZZwork,ZZwork)); 798 799 CHKERRQ(VecCopy(XXwork,user->utrue)); 800 CHKERRQ(VecAXPY(user->utrue,1.0,YYwork)); 801 CHKERRQ(VecAXPY(user->utrue,1.0,ZZwork)); 802 CHKERRQ(VecScale(user->utrue,-10.0)); 803 CHKERRQ(VecExp(user->utrue)); 804 CHKERRQ(VecShift(user->utrue,0.5)); 805 806 CHKERRQ(VecDestroy(&XX)); 807 CHKERRQ(VecDestroy(&YY)); 808 CHKERRQ(VecDestroy(&ZZ)); 809 CHKERRQ(VecDestroy(&XXwork)); 810 CHKERRQ(VecDestroy(&YYwork)); 811 CHKERRQ(VecDestroy(&ZZwork)); 812 CHKERRQ(VecDestroy(&UTwork)); 813 814 /* Initial guess and reference model */ 815 CHKERRQ(VecDuplicate(user->utrue,&user->ur)); 816 CHKERRQ(VecSet(user->ur,0.0)); 817 818 /* Generate Grad matrix */ 819 CHKERRQ(MatCreate(PETSC_COMM_WORLD,&user->Grad)); 820 CHKERRQ(MatSetSizes(user->Grad,PETSC_DECIDE,PETSC_DECIDE,m,n)); 821 CHKERRQ(MatSetFromOptions(user->Grad)); 822 CHKERRQ(MatMPIAIJSetPreallocation(user->Grad,2,NULL,2,NULL)); 823 CHKERRQ(MatSeqAIJSetPreallocation(user->Grad,2,NULL)); 824 CHKERRQ(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 CHKERRQ(MatSetValues(user->Grad,1,&i,1,&j,&neg_hinv,INSERT_VALUES)); 831 j = j+1; 832 CHKERRQ(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 CHKERRQ(MatSetValues(user->Grad,1,&i,1,&j,&neg_hinv,INSERT_VALUES)); 838 j = j + user->mx; 839 CHKERRQ(MatSetValues(user->Grad,1,&i,1,&j,&hinv,INSERT_VALUES)); 840 } 841 if (i>=2*m/3) { 842 j = i-2*m/3; 843 CHKERRQ(MatSetValues(user->Grad,1,&i,1,&j,&neg_hinv,INSERT_VALUES)); 844 j = j + user->mx*user->mx; 845 CHKERRQ(MatSetValues(user->Grad,1,&i,1,&j,&hinv,INSERT_VALUES)); 846 } 847 } 848 849 CHKERRQ(MatAssemblyBegin(user->Grad,MAT_FINAL_ASSEMBLY)); 850 CHKERRQ(MatAssemblyEnd(user->Grad,MAT_FINAL_ASSEMBLY)); 851 852 /* Generate arithmetic averaging matrix Av */ 853 CHKERRQ(MatCreate(PETSC_COMM_WORLD,&user->Av)); 854 CHKERRQ(MatSetSizes(user->Av,PETSC_DECIDE,PETSC_DECIDE,m,n)); 855 CHKERRQ(MatSetFromOptions(user->Av)); 856 CHKERRQ(MatMPIAIJSetPreallocation(user->Av,2,NULL,2,NULL)); 857 CHKERRQ(MatSeqAIJSetPreallocation(user->Av,2,NULL)); 858 CHKERRQ(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 CHKERRQ(MatSetValues(user->Av,1,&i,1,&j,&half,INSERT_VALUES)); 865 j = j+1; 866 CHKERRQ(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 CHKERRQ(MatSetValues(user->Av,1,&i,1,&j,&half,INSERT_VALUES)); 872 j = j + user->mx; 873 CHKERRQ(MatSetValues(user->Av,1,&i,1,&j,&half,INSERT_VALUES)); 874 } 875 if (i>=2*m/3) { 876 j = i-2*m/3; 877 CHKERRQ(MatSetValues(user->Av,1,&i,1,&j,&half,INSERT_VALUES)); 878 j = j + user->mx*user->mx; 879 CHKERRQ(MatSetValues(user->Av,1,&i,1,&j,&half,INSERT_VALUES)); 880 } 881 } 882 883 CHKERRQ(MatAssemblyBegin(user->Av,MAT_FINAL_ASSEMBLY)); 884 CHKERRQ(MatAssemblyEnd(user->Av,MAT_FINAL_ASSEMBLY)); 885 886 /* Generate transpose of averaging matrix Av */ 887 CHKERRQ(MatTranspose(user->Av,MAT_INITIAL_MATRIX,&user->AvT)); 888 889 CHKERRQ(MatCreate(PETSC_COMM_WORLD,&user->L)); 890 CHKERRQ(MatSetSizes(user->L,PETSC_DECIDE,PETSC_DECIDE,m+n,n)); 891 CHKERRQ(MatSetFromOptions(user->L)); 892 CHKERRQ(MatMPIAIJSetPreallocation(user->L,2,NULL,2,NULL)); 893 CHKERRQ(MatSeqAIJSetPreallocation(user->L,2,NULL)); 894 CHKERRQ(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 CHKERRQ(MatSetValues(user->L,1,&i,1,&j,&neg_hinv,INSERT_VALUES)); 901 j = j+1; 902 CHKERRQ(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 CHKERRQ(MatSetValues(user->L,1,&i,1,&j,&neg_hinv,INSERT_VALUES)); 908 j = j + user->mx; 909 CHKERRQ(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 CHKERRQ(MatSetValues(user->L,1,&i,1,&j,&neg_hinv,INSERT_VALUES)); 914 j = j + user->mx*user->mx; 915 CHKERRQ(MatSetValues(user->L,1,&i,1,&j,&hinv,INSERT_VALUES)); 916 } 917 if (i>=m) { 918 j = i - m; 919 CHKERRQ(MatSetValues(user->L,1,&i,1,&j,&sqrt_beta,INSERT_VALUES)); 920 } 921 } 922 923 CHKERRQ(MatAssemblyBegin(user->L,MAT_FINAL_ASSEMBLY)); 924 CHKERRQ(MatAssemblyEnd(user->L,MAT_FINAL_ASSEMBLY)); 925 926 CHKERRQ(MatScale(user->L,PetscPowScalar(h,1.5))); 927 928 /* Generate Div matrix */ 929 CHKERRQ(MatTranspose(user->Grad,MAT_INITIAL_MATRIX,&user->Div)); 930 931 /* Build work vectors and matrices */ 932 CHKERRQ(VecCreate(PETSC_COMM_WORLD,&user->S)); 933 CHKERRQ(VecSetSizes(user->S, PETSC_DECIDE, user->mx*user->mx*(user->mx-1)*3)); 934 CHKERRQ(VecSetFromOptions(user->S)); 935 936 CHKERRQ(VecCreate(PETSC_COMM_WORLD,&user->lwork)); 937 CHKERRQ(VecSetSizes(user->lwork,PETSC_DECIDE,m+user->mx*user->mx*user->mx)); 938 CHKERRQ(VecSetFromOptions(user->lwork)); 939 940 CHKERRQ(MatDuplicate(user->Div,MAT_SHARE_NONZERO_PATTERN,&user->Divwork)); 941 CHKERRQ(MatDuplicate(user->Av,MAT_SHARE_NONZERO_PATTERN,&user->Avwork)); 942 943 CHKERRQ(VecCreate(PETSC_COMM_WORLD,&user->d)); 944 CHKERRQ(VecSetSizes(user->d,PETSC_DECIDE,user->ndata*user->nt)); 945 CHKERRQ(VecSetFromOptions(user->d)); 946 947 CHKERRQ(VecDuplicate(user->S,&user->Swork)); 948 CHKERRQ(VecDuplicate(user->S,&user->Av_u)); 949 CHKERRQ(VecDuplicate(user->S,&user->Twork)); 950 CHKERRQ(VecDuplicate(user->S,&user->Rwork)); 951 CHKERRQ(VecDuplicate(user->y,&user->ywork)); 952 CHKERRQ(VecDuplicate(user->u,&user->uwork)); 953 CHKERRQ(VecDuplicate(user->u,&user->js_diag)); 954 CHKERRQ(VecDuplicate(user->c,&user->cwork)); 955 CHKERRQ(VecDuplicate(user->d,&user->dwork)); 956 957 /* Create matrix-free shell user->Js for computing A*x */ 958 CHKERRQ(MatCreateShell(PETSC_COMM_WORLD,PETSC_DETERMINE,PETSC_DETERMINE,user->m*user->nt,user->m*user->nt,user,&user->Js)); 959 CHKERRQ(MatShellSetOperation(user->Js,MATOP_MULT,(void(*)(void))StateMatMult)); 960 CHKERRQ(MatShellSetOperation(user->Js,MATOP_DUPLICATE,(void(*)(void))StateMatDuplicate)); 961 CHKERRQ(MatShellSetOperation(user->Js,MATOP_MULT_TRANSPOSE,(void(*)(void))StateMatMultTranspose)); 962 CHKERRQ(MatShellSetOperation(user->Js,MATOP_GET_DIAGONAL,(void(*)(void))StateMatGetDiagonal)); 963 964 /* Diagonal blocks of user->Js */ 965 CHKERRQ(MatCreateShell(PETSC_COMM_WORLD,PETSC_DETERMINE,PETSC_DETERMINE,user->m,user->m,user,&user->JsBlock)); 966 CHKERRQ(MatShellSetOperation(user->JsBlock,MATOP_MULT,(void(*)(void))StateMatBlockMult)); 967 /* Blocks are symmetric */ 968 CHKERRQ(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 CHKERRQ(MatCreateShell(PETSC_COMM_WORLD,PETSC_DETERMINE,PETSC_DETERMINE,user->m,user->m,user,&user->JsBlockPrec)); 974 CHKERRQ(MatShellSetOperation(user->JsBlockPrec,MATOP_MULT,(void(*)(void))StateMatBlockPrecMult)); 975 /* JsBlockPrec is symmetric */ 976 CHKERRQ(MatShellSetOperation(user->JsBlockPrec,MATOP_MULT_TRANSPOSE,(void(*)(void))StateMatBlockPrecMult)); 977 CHKERRQ(MatSetOption(user->JsBlockPrec,MAT_SYMMETRY_ETERNAL,PETSC_TRUE)); 978 979 /* Create a matrix-free shell user->Jd for computing B*x */ 980 CHKERRQ(MatCreateShell(PETSC_COMM_WORLD,PETSC_DETERMINE,PETSC_DETERMINE,user->m*user->nt,user->m,user,&user->Jd)); 981 CHKERRQ(MatShellSetOperation(user->Jd,MATOP_MULT,(void(*)(void))DesignMatMult)); 982 CHKERRQ(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 CHKERRQ(MatCreateShell(PETSC_COMM_WORLD,PETSC_DETERMINE,PETSC_DETERMINE,user->m*user->nt,user->m*user->nt,user,&user->JsInv)); 986 CHKERRQ(MatShellSetOperation(user->JsInv,MATOP_MULT,(void(*)(void))StateMatInvMult)); 987 CHKERRQ(MatShellSetOperation(user->JsInv,MATOP_MULT_TRANSPOSE,(void(*)(void))StateMatInvTransposeMult)); 988 989 /* Solver options and tolerances */ 990 CHKERRQ(KSPCreate(PETSC_COMM_WORLD,&user->solver)); 991 CHKERRQ(KSPSetType(user->solver,KSPCG)); 992 CHKERRQ(KSPSetOperators(user->solver,user->JsBlock,user->JsBlockPrec)); 993 CHKERRQ(KSPSetInitialGuessNonzero(user->solver,PETSC_FALSE)); 994 CHKERRQ(KSPSetTolerances(user->solver,1e-4,1e-20,1e3,500)); 995 CHKERRQ(KSPSetFromOptions(user->solver)); 996 CHKERRQ(KSPGetPC(user->solver,&user->prec)); 997 CHKERRQ(PCSetType(user->prec,PCSHELL)); 998 999 CHKERRQ(PCShellSetApply(user->prec,StateMatBlockPrecMult)); 1000 CHKERRQ(PCShellSetApplyTranspose(user->prec,StateMatBlockPrecMult)); 1001 CHKERRQ(PCShellSetContext(user->prec,user)); 1002 1003 /* Create scatter from y to y_1,y_2,...,y_nt */ 1004 CHKERRQ(PetscMalloc1(user->nt*user->m,&user->yi_scatter)); 1005 CHKERRQ(VecCreate(PETSC_COMM_WORLD,&yi)); 1006 CHKERRQ(VecSetSizes(yi,PETSC_DECIDE,user->mx*user->mx*user->mx)); 1007 CHKERRQ(VecSetFromOptions(yi)); 1008 CHKERRQ(VecDuplicateVecs(yi,user->nt,&user->yi)); 1009 CHKERRQ(VecDuplicateVecs(yi,user->nt,&user->yiwork)); 1010 1011 CHKERRQ(VecGetOwnershipRange(user->y,&lo2,&hi2)); 1012 istart = 0; 1013 for (i=0; i<user->nt; i++) { 1014 CHKERRQ(VecGetOwnershipRange(user->yi[i],&lo,&hi)); 1015 CHKERRQ(ISCreateStride(PETSC_COMM_SELF,hi-lo,lo,1,&is_to_yi)); 1016 CHKERRQ(ISCreateStride(PETSC_COMM_SELF,hi-lo,lo2+istart,1,&is_from_y)); 1017 CHKERRQ(VecScatterCreate(user->y,is_from_y,user->yi[i],is_to_yi,&user->yi_scatter[i])); 1018 istart = istart + hi-lo; 1019 CHKERRQ(ISDestroy(&is_to_yi)); 1020 CHKERRQ(ISDestroy(&is_from_y)); 1021 } 1022 CHKERRQ(VecDestroy(&yi)); 1023 1024 /* Create scatter from d to d_1,d_2,...,d_ns */ 1025 CHKERRQ(PetscMalloc1(user->ns*user->ndata,&user->di_scatter)); 1026 CHKERRQ(VecCreate(PETSC_COMM_WORLD,&di)); 1027 CHKERRQ(VecSetSizes(di,PETSC_DECIDE,user->ndata)); 1028 CHKERRQ(VecSetFromOptions(di)); 1029 CHKERRQ(VecDuplicateVecs(di,user->ns,&user->di)); 1030 CHKERRQ(VecGetOwnershipRange(user->d,&lo2,&hi2)); 1031 istart = 0; 1032 for (i=0; i<user->ns; i++) { 1033 CHKERRQ(VecGetOwnershipRange(user->di[i],&lo,&hi)); 1034 CHKERRQ(ISCreateStride(PETSC_COMM_SELF,hi-lo,lo,1,&is_to_di)); 1035 CHKERRQ(ISCreateStride(PETSC_COMM_SELF,hi-lo,lo2+istart,1,&is_from_d)); 1036 CHKERRQ(VecScatterCreate(user->d,is_from_d,user->di[i],is_to_di,&user->di_scatter[i])); 1037 istart = istart + hi-lo; 1038 CHKERRQ(ISDestroy(&is_to_di)); 1039 CHKERRQ(ISDestroy(&is_from_d)); 1040 } 1041 CHKERRQ(VecDestroy(&di)); 1042 1043 /* Assemble RHS of forward problem */ 1044 CHKERRQ(VecCopy(bc,user->yiwork[0])); 1045 for (i=1; i<user->nt; i++) { 1046 CHKERRQ(VecSet(user->yiwork[i],0.0)); 1047 } 1048 CHKERRQ(Gather_i(user->q,user->yiwork,user->yi_scatter,user->nt)); 1049 CHKERRQ(VecDestroy(&bc)); 1050 1051 /* Compute true state function yt given ut */ 1052 CHKERRQ(VecCreate(PETSC_COMM_WORLD,&user->ytrue)); 1053 CHKERRQ(VecSetSizes(user->ytrue,PETSC_DECIDE,n*user->nt)); 1054 CHKERRQ(VecSetFromOptions(user->ytrue)); 1055 1056 /* First compute Av_u = Av*exp(-u) */ 1057 CHKERRQ(VecSet(user->uwork,0)); 1058 CHKERRQ(VecAXPY(user->uwork,-1.0,user->utrue)); /* Note: user->utrue */ 1059 CHKERRQ(VecExp(user->uwork)); 1060 CHKERRQ(MatMult(user->Av,user->uwork,user->Av_u)); 1061 1062 /* Symbolic DSG = Div * Grad */ 1063 CHKERRQ(MatProductCreate(user->Div,user->Grad,NULL,&user->DSG)); 1064 CHKERRQ(MatProductSetType(user->DSG,MATPRODUCT_AB)); 1065 CHKERRQ(MatProductSetAlgorithm(user->DSG,"default")); 1066 CHKERRQ(MatProductSetFill(user->DSG,PETSC_DEFAULT)); 1067 CHKERRQ(MatProductSetFromOptions(user->DSG)); 1068 CHKERRQ(MatProductSymbolic(user->DSG)); 1069 1070 user->dsg_formed = PETSC_TRUE; 1071 1072 /* Next form DSG = Div*Grad */ 1073 CHKERRQ(MatCopy(user->Div,user->Divwork,SAME_NONZERO_PATTERN)); 1074 CHKERRQ(MatDiagonalScale(user->Divwork,NULL,user->Av_u)); 1075 if (user->dsg_formed) { 1076 CHKERRQ(MatProductNumeric(user->DSG)); 1077 } else { 1078 CHKERRQ(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 CHKERRQ(MatScale(user->DSG,user->ht)); 1083 CHKERRQ(MatShift(user->DSG,1.0)); 1084 1085 /* Now solve for ytrue */ 1086 CHKERRQ(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 CHKERRQ(VecSet(user->uwork,0)); 1092 CHKERRQ(VecAXPY(user->uwork,-1.0,user->u)); /* Note: user->u */ 1093 CHKERRQ(VecExp(user->uwork)); 1094 CHKERRQ(MatMult(user->Av,user->uwork,user->Av_u)); 1095 1096 /* Next form DSG = Div*S*Grad */ 1097 CHKERRQ(MatCopy(user->Div,user->Divwork,SAME_NONZERO_PATTERN)); 1098 CHKERRQ(MatDiagonalScale(user->Divwork,NULL,user->Av_u)); 1099 if (user->dsg_formed) { 1100 CHKERRQ(MatProductNumeric(user->DSG)); 1101 } else { 1102 CHKERRQ(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 CHKERRQ(MatScale(user->DSG,user->ht)); 1108 CHKERRQ(MatShift(user->DSG,1.0)); 1109 1110 /* Now solve for y */ 1111 CHKERRQ(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 CHKERRQ(MatCreate(PETSC_COMM_WORLD,&user->Qblock)); 1115 CHKERRQ(MatSetSizes(user->Qblock,PETSC_DECIDE,PETSC_DECIDE,user->ndata,n)); 1116 CHKERRQ(MatSetFromOptions(user->Qblock)); 1117 CHKERRQ(MatMPIAIJSetPreallocation(user->Qblock,8,NULL,8,NULL)); 1118 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(MatSetValues(user->Qblock,1,&i,1,&j,&v,INSERT_VALUES)); 1200 } 1201 CHKERRQ(MatAssemblyBegin(user->Qblock,MAT_FINAL_ASSEMBLY)); 1202 CHKERRQ(MatAssemblyEnd(user->Qblock,MAT_FINAL_ASSEMBLY)); 1203 1204 CHKERRQ(MatTranspose(user->Qblock,MAT_INITIAL_MATRIX,&user->QblockT)); 1205 CHKERRQ(MatTranspose(user->L,MAT_INITIAL_MATRIX,&user->LT)); 1206 1207 /* Add noise to the measurement data */ 1208 CHKERRQ(VecSet(user->ywork,1.0)); 1209 CHKERRQ(VecAYPX(user->ywork,user->noise,user->ytrue)); 1210 CHKERRQ(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 CHKERRQ(MatMult(user->Qblock,user->yiwork[i],user->di[j])); 1214 } 1215 CHKERRQ(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 CHKERRQ(KSPSetFromOptions(user->solver)); 1219 CHKERRQ(PetscFree(x)); 1220 CHKERRQ(PetscFree(y)); 1221 CHKERRQ(PetscFree(z)); 1222 PetscFunctionReturn(0); 1223 } 1224 1225 PetscErrorCode ParabolicDestroy(AppCtx *user) 1226 { 1227 PetscInt i; 1228 1229 PetscFunctionBegin; 1230 CHKERRQ(MatDestroy(&user->Qblock)); 1231 CHKERRQ(MatDestroy(&user->QblockT)); 1232 CHKERRQ(MatDestroy(&user->Div)); 1233 CHKERRQ(MatDestroy(&user->Divwork)); 1234 CHKERRQ(MatDestroy(&user->Grad)); 1235 CHKERRQ(MatDestroy(&user->Av)); 1236 CHKERRQ(MatDestroy(&user->Avwork)); 1237 CHKERRQ(MatDestroy(&user->AvT)); 1238 CHKERRQ(MatDestroy(&user->DSG)); 1239 CHKERRQ(MatDestroy(&user->L)); 1240 CHKERRQ(MatDestroy(&user->LT)); 1241 CHKERRQ(KSPDestroy(&user->solver)); 1242 CHKERRQ(MatDestroy(&user->Js)); 1243 CHKERRQ(MatDestroy(&user->Jd)); 1244 CHKERRQ(MatDestroy(&user->JsInv)); 1245 CHKERRQ(MatDestroy(&user->JsBlock)); 1246 CHKERRQ(MatDestroy(&user->JsBlockPrec)); 1247 CHKERRQ(VecDestroy(&user->u)); 1248 CHKERRQ(VecDestroy(&user->uwork)); 1249 CHKERRQ(VecDestroy(&user->utrue)); 1250 CHKERRQ(VecDestroy(&user->y)); 1251 CHKERRQ(VecDestroy(&user->ywork)); 1252 CHKERRQ(VecDestroy(&user->ytrue)); 1253 CHKERRQ(VecDestroyVecs(user->nt,&user->yi)); 1254 CHKERRQ(VecDestroyVecs(user->nt,&user->yiwork)); 1255 CHKERRQ(VecDestroyVecs(user->ns,&user->di)); 1256 CHKERRQ(PetscFree(user->yi)); 1257 CHKERRQ(PetscFree(user->yiwork)); 1258 CHKERRQ(PetscFree(user->di)); 1259 CHKERRQ(VecDestroy(&user->c)); 1260 CHKERRQ(VecDestroy(&user->cwork)); 1261 CHKERRQ(VecDestroy(&user->ur)); 1262 CHKERRQ(VecDestroy(&user->q)); 1263 CHKERRQ(VecDestroy(&user->d)); 1264 CHKERRQ(VecDestroy(&user->dwork)); 1265 CHKERRQ(VecDestroy(&user->lwork)); 1266 CHKERRQ(VecDestroy(&user->S)); 1267 CHKERRQ(VecDestroy(&user->Swork)); 1268 CHKERRQ(VecDestroy(&user->Av_u)); 1269 CHKERRQ(VecDestroy(&user->Twork)); 1270 CHKERRQ(VecDestroy(&user->Rwork)); 1271 CHKERRQ(VecDestroy(&user->js_diag)); 1272 CHKERRQ(ISDestroy(&user->s_is)); 1273 CHKERRQ(ISDestroy(&user->d_is)); 1274 CHKERRQ(VecScatterDestroy(&user->state_scatter)); 1275 CHKERRQ(VecScatterDestroy(&user->design_scatter)); 1276 for (i=0; i<user->nt; i++) { 1277 CHKERRQ(VecScatterDestroy(&user->yi_scatter[i])); 1278 } 1279 for (i=0; i<user->ns; i++) { 1280 CHKERRQ(VecScatterDestroy(&user->di_scatter[i])); 1281 } 1282 CHKERRQ(PetscFree(user->yi_scatter)); 1283 CHKERRQ(PetscFree(user->di_scatter)); 1284 CHKERRQ(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 CHKERRQ(TaoGetSolution(tao,&X)); 1296 CHKERRQ(Scatter(X,user->ywork,user->state_scatter,user->uwork,user->design_scatter)); 1297 CHKERRQ(VecAXPY(user->ywork,-1.0,user->ytrue)); 1298 CHKERRQ(VecAXPY(user->uwork,-1.0,user->utrue)); 1299 CHKERRQ(VecNorm(user->uwork,NORM_2,&unorm)); 1300 CHKERRQ(VecNorm(user->ywork,NORM_2,&ynorm)); 1301 CHKERRQ(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