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