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