1 /* Program usage: mpiexec -n 1 eptorsion1 [-help] [all TAO options] */ 2 3 /* ---------------------------------------------------------------------- 4 5 Elastic-plastic torsion problem. 6 7 The elastic plastic torsion problem arises from the determination 8 of the stress field on an infinitely long cylindrical bar, which is 9 equivalent to the solution of the following problem: 10 11 min{ .5 * integral(||gradient(v(x))||^2 dx) - C * integral(v(x) dx)} 12 13 where C is the torsion angle per unit length. 14 15 The multiprocessor version of this code is eptorsion2.c. 16 17 ---------------------------------------------------------------------- */ 18 19 /* 20 Include "petsctao.h" so that we can use TAO solvers. Note that this 21 file automatically includes files for lower-level support, such as those 22 provided by the PETSc library: 23 petsc.h - base PETSc routines petscvec.h - vectors 24 petscsys.h - system routines petscmat.h - matrices 25 petscis.h - index sets petscksp.h - Krylov subspace methods 26 petscviewer.h - viewers petscpc.h - preconditioners 27 */ 28 29 #include <petsctao.h> 30 31 static char help[]= 32 "Demonstrates use of the TAO package to solve \n\ 33 unconstrained minimization problems on a single processor. This example \n\ 34 is based on the Elastic-Plastic Torsion (dept) problem from the MINPACK-2 \n\ 35 test suite.\n\ 36 The command line options are:\n\ 37 -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction\n\ 38 -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction\n\ 39 -par <param>, where <param> = angle of twist per unit length\n\n"; 40 41 /*T 42 Concepts: TAO^Solving an unconstrained minimization problem 43 Routines: TaoCreate(); TaoSetType(); 44 Routines: TaoSetSolution(); 45 Routines: TaoSetObjectiveAndGradient(); 46 Routines: TaoSetHessian(); TaoSetFromOptions(); 47 Routines: TaoGetKSP(); TaoSolve(); 48 Routines: TaoDestroy(); 49 Processors: 1 50 T*/ 51 52 /* 53 User-defined application context - contains data needed by the 54 application-provided call-back routines, FormFunction(), 55 FormGradient(), and FormHessian(). 56 */ 57 58 typedef struct { 59 PetscReal param; /* nonlinearity parameter */ 60 PetscInt mx, my; /* discretization in x- and y-directions */ 61 PetscInt ndim; /* problem dimension */ 62 Vec s, y, xvec; /* work space for computing Hessian */ 63 PetscReal hx, hy; /* mesh spacing in x- and y-directions */ 64 } AppCtx; 65 66 /* -------- User-defined Routines --------- */ 67 68 PetscErrorCode FormInitialGuess(AppCtx*,Vec); 69 PetscErrorCode FormFunction(Tao,Vec,PetscReal*,void*); 70 PetscErrorCode FormGradient(Tao,Vec,Vec,void*); 71 PetscErrorCode FormHessian(Tao,Vec,Mat,Mat, void*); 72 PetscErrorCode HessianProductMat(Mat,Vec,Vec); 73 PetscErrorCode HessianProduct(void*,Vec,Vec); 74 PetscErrorCode MatrixFreeHessian(Tao,Vec,Mat,Mat,void*); 75 PetscErrorCode FormFunctionGradient(Tao,Vec,PetscReal *,Vec,void *); 76 77 PetscErrorCode main(int argc,char **argv) 78 { 79 PetscInt mx=10; /* discretization in x-direction */ 80 PetscInt my=10; /* discretization in y-direction */ 81 Vec x; /* solution, gradient vectors */ 82 PetscBool flg; /* A return value when checking for use options */ 83 Tao tao; /* Tao solver context */ 84 Mat H; /* Hessian matrix */ 85 AppCtx user; /* application context */ 86 PetscMPIInt size; /* number of processes */ 87 PetscReal one=1.0; 88 89 PetscBool test_lmvm = PETSC_FALSE; 90 KSP ksp; 91 PC pc; 92 Mat M; 93 Vec in, out, out2; 94 PetscReal mult_solve_dist; 95 96 /* Initialize TAO,PETSc */ 97 PetscCall(PetscInitialize(&argc,&argv,(char *)0,help)); 98 PetscCallMPI(MPI_Comm_size(MPI_COMM_WORLD,&size)); 99 PetscCheck(size == 1,PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"Incorrect number of processors"); 100 101 /* Specify default parameters for the problem, check for command-line overrides */ 102 user.param = 5.0; 103 PetscCall(PetscOptionsGetInt(NULL,NULL,"-my",&my,&flg)); 104 PetscCall(PetscOptionsGetInt(NULL,NULL,"-mx",&mx,&flg)); 105 PetscCall(PetscOptionsGetReal(NULL,NULL,"-par",&user.param,&flg)); 106 PetscCall(PetscOptionsGetBool(NULL,NULL,"-test_lmvm",&test_lmvm,&flg)); 107 108 PetscCall(PetscPrintf(PETSC_COMM_SELF,"\n---- Elastic-Plastic Torsion Problem -----\n")); 109 PetscCall(PetscPrintf(PETSC_COMM_SELF,"mx: %D my: %D \n\n",mx,my)); 110 user.ndim = mx * my; user.mx = mx; user.my = my; 111 user.hx = one/(mx+1); user.hy = one/(my+1); 112 113 /* Allocate vectors */ 114 PetscCall(VecCreateSeq(PETSC_COMM_SELF,user.ndim,&user.y)); 115 PetscCall(VecDuplicate(user.y,&user.s)); 116 PetscCall(VecDuplicate(user.y,&x)); 117 118 /* Create TAO solver and set desired solution method */ 119 PetscCall(TaoCreate(PETSC_COMM_SELF,&tao)); 120 PetscCall(TaoSetType(tao,TAOLMVM)); 121 122 /* Set solution vector with an initial guess */ 123 PetscCall(FormInitialGuess(&user,x)); 124 PetscCall(TaoSetSolution(tao,x)); 125 126 /* Set routine for function and gradient evaluation */ 127 PetscCall(TaoSetObjectiveAndGradient(tao,NULL,FormFunctionGradient,(void *)&user)); 128 129 /* From command line options, determine if using matrix-free hessian */ 130 PetscCall(PetscOptionsHasName(NULL,NULL,"-my_tao_mf",&flg)); 131 if (flg) { 132 PetscCall(MatCreateShell(PETSC_COMM_SELF,user.ndim,user.ndim,user.ndim,user.ndim,(void*)&user,&H)); 133 PetscCall(MatShellSetOperation(H,MATOP_MULT,(void(*)(void))HessianProductMat)); 134 PetscCall(MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE)); 135 136 PetscCall(TaoSetHessian(tao,H,H,MatrixFreeHessian,(void *)&user)); 137 } else { 138 PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF,user.ndim,user.ndim,5,NULL,&H)); 139 PetscCall(MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE)); 140 PetscCall(TaoSetHessian(tao,H,H,FormHessian,(void *)&user)); 141 } 142 143 /* Test the LMVM matrix */ 144 if (test_lmvm) { 145 PetscCall(PetscOptionsSetValue(NULL, "-tao_type", "bntr")); 146 PetscCall(PetscOptionsSetValue(NULL, "-tao_bnk_pc_type", "lmvm")); 147 } 148 149 /* Check for any TAO command line options */ 150 PetscCall(TaoSetFromOptions(tao)); 151 152 /* SOLVE THE APPLICATION */ 153 PetscCall(TaoSolve(tao)); 154 155 /* Test the LMVM matrix */ 156 if (test_lmvm) { 157 PetscCall(TaoGetKSP(tao, &ksp)); 158 PetscCall(KSPGetPC(ksp, &pc)); 159 PetscCall(PCLMVMGetMatLMVM(pc, &M)); 160 PetscCall(VecDuplicate(x, &in)); 161 PetscCall(VecDuplicate(x, &out)); 162 PetscCall(VecDuplicate(x, &out2)); 163 PetscCall(VecSet(in, 5.0)); 164 PetscCall(MatMult(M, in, out)); 165 PetscCall(MatSolve(M, out, out2)); 166 PetscCall(VecAXPY(out2, -1.0, in)); 167 PetscCall(VecNorm(out2, NORM_2, &mult_solve_dist)); 168 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)tao), "error between MatMult and MatSolve: %e\n", mult_solve_dist)); 169 PetscCall(VecDestroy(&in)); 170 PetscCall(VecDestroy(&out)); 171 PetscCall(VecDestroy(&out2)); 172 } 173 174 PetscCall(TaoDestroy(&tao)); 175 PetscCall(VecDestroy(&user.s)); 176 PetscCall(VecDestroy(&user.y)); 177 PetscCall(VecDestroy(&x)); 178 PetscCall(MatDestroy(&H)); 179 180 PetscCall(PetscFinalize()); 181 return 0; 182 } 183 184 /* ------------------------------------------------------------------- */ 185 /* 186 FormInitialGuess - Computes an initial approximation to the solution. 187 188 Input Parameters: 189 . user - user-defined application context 190 . X - vector 191 192 Output Parameters: 193 . X - vector 194 */ 195 PetscErrorCode FormInitialGuess(AppCtx *user,Vec X) 196 { 197 PetscReal hx = user->hx, hy = user->hy, temp; 198 PetscReal val; 199 PetscInt i, j, k, nx = user->mx, ny = user->my; 200 201 /* Compute initial guess */ 202 PetscFunctionBeginUser; 203 for (j=0; j<ny; j++) { 204 temp = PetscMin(j+1,ny-j)*hy; 205 for (i=0; i<nx; i++) { 206 k = nx*j + i; 207 val = PetscMin((PetscMin(i+1,nx-i))*hx,temp); 208 PetscCall(VecSetValues(X,1,&k,&val,ADD_VALUES)); 209 } 210 } 211 PetscCall(VecAssemblyBegin(X)); 212 PetscCall(VecAssemblyEnd(X)); 213 PetscFunctionReturn(0); 214 } 215 216 /* ------------------------------------------------------------------- */ 217 /* 218 FormFunctionGradient - Evaluates the function and corresponding gradient. 219 220 Input Parameters: 221 tao - the Tao context 222 X - the input vector 223 ptr - optional user-defined context, as set by TaoSetFunction() 224 225 Output Parameters: 226 f - the newly evaluated function 227 G - the newly evaluated gradient 228 */ 229 PetscErrorCode FormFunctionGradient(Tao tao,Vec X,PetscReal *f,Vec G,void *ptr) 230 { 231 PetscFunctionBeginUser; 232 PetscCall(FormFunction(tao,X,f,ptr)); 233 PetscCall(FormGradient(tao,X,G,ptr)); 234 PetscFunctionReturn(0); 235 } 236 237 /* ------------------------------------------------------------------- */ 238 /* 239 FormFunction - Evaluates the function, f(X). 240 241 Input Parameters: 242 . tao - the Tao context 243 . X - the input vector 244 . ptr - optional user-defined context, as set by TaoSetFunction() 245 246 Output Parameters: 247 . f - the newly evaluated function 248 */ 249 PetscErrorCode FormFunction(Tao tao,Vec X,PetscReal *f,void *ptr) 250 { 251 AppCtx *user = (AppCtx *) ptr; 252 PetscReal hx = user->hx, hy = user->hy, area, three = 3.0, p5 = 0.5; 253 PetscReal zero = 0.0, vb, vl, vr, vt, dvdx, dvdy, flin = 0.0, fquad = 0.0; 254 PetscReal v, cdiv3 = user->param/three; 255 const PetscScalar *x; 256 PetscInt nx = user->mx, ny = user->my, i, j, k; 257 258 PetscFunctionBeginUser; 259 /* Get pointer to vector data */ 260 PetscCall(VecGetArrayRead(X,&x)); 261 262 /* Compute function contributions over the lower triangular elements */ 263 for (j=-1; j<ny; j++) { 264 for (i=-1; i<nx; i++) { 265 k = nx*j + i; 266 v = zero; 267 vr = zero; 268 vt = zero; 269 if (i >= 0 && j >= 0) v = x[k]; 270 if (i < nx-1 && j > -1) vr = x[k+1]; 271 if (i > -1 && j < ny-1) vt = x[k+nx]; 272 dvdx = (vr-v)/hx; 273 dvdy = (vt-v)/hy; 274 fquad += dvdx*dvdx + dvdy*dvdy; 275 flin -= cdiv3*(v+vr+vt); 276 } 277 } 278 279 /* Compute function contributions over the upper triangular elements */ 280 for (j=0; j<=ny; j++) { 281 for (i=0; i<=nx; i++) { 282 k = nx*j + i; 283 vb = zero; 284 vl = zero; 285 v = zero; 286 if (i < nx && j > 0) vb = x[k-nx]; 287 if (i > 0 && j < ny) vl = x[k-1]; 288 if (i < nx && j < ny) v = x[k]; 289 dvdx = (v-vl)/hx; 290 dvdy = (v-vb)/hy; 291 fquad = fquad + dvdx*dvdx + dvdy*dvdy; 292 flin = flin - cdiv3*(vb+vl+v); 293 } 294 } 295 296 /* Restore vector */ 297 PetscCall(VecRestoreArrayRead(X,&x)); 298 299 /* Scale the function */ 300 area = p5*hx*hy; 301 *f = area*(p5*fquad+flin); 302 303 PetscCall(PetscLogFlops(24.0*nx*ny)); 304 PetscFunctionReturn(0); 305 } 306 307 /* ------------------------------------------------------------------- */ 308 /* 309 FormGradient - Evaluates the gradient, G(X). 310 311 Input Parameters: 312 . tao - the Tao context 313 . X - input vector 314 . ptr - optional user-defined context 315 316 Output Parameters: 317 . G - vector containing the newly evaluated gradient 318 */ 319 PetscErrorCode FormGradient(Tao tao,Vec X,Vec G,void *ptr) 320 { 321 AppCtx *user = (AppCtx *) ptr; 322 PetscReal zero=0.0, p5=0.5, three = 3.0, area, val; 323 PetscInt nx = user->mx, ny = user->my, ind, i, j, k; 324 PetscReal hx = user->hx, hy = user->hy; 325 PetscReal vb, vl, vr, vt, dvdx, dvdy; 326 PetscReal v, cdiv3 = user->param/three; 327 const PetscScalar *x; 328 329 PetscFunctionBeginUser; 330 /* Initialize gradient to zero */ 331 PetscCall(VecSet(G, zero)); 332 333 /* Get pointer to vector data */ 334 PetscCall(VecGetArrayRead(X,&x)); 335 336 /* Compute gradient contributions over the lower triangular elements */ 337 for (j=-1; j<ny; j++) { 338 for (i=-1; i<nx; i++) { 339 k = nx*j + i; 340 v = zero; 341 vr = zero; 342 vt = zero; 343 if (i >= 0 && j >= 0) v = x[k]; 344 if (i < nx-1 && j > -1) vr = x[k+1]; 345 if (i > -1 && j < ny-1) vt = x[k+nx]; 346 dvdx = (vr-v)/hx; 347 dvdy = (vt-v)/hy; 348 if (i != -1 && j != -1) { 349 ind = k; val = - dvdx/hx - dvdy/hy - cdiv3; 350 PetscCall(VecSetValues(G,1,&ind,&val,ADD_VALUES)); 351 } 352 if (i != nx-1 && j != -1) { 353 ind = k+1; val = dvdx/hx - cdiv3; 354 PetscCall(VecSetValues(G,1,&ind,&val,ADD_VALUES)); 355 } 356 if (i != -1 && j != ny-1) { 357 ind = k+nx; val = dvdy/hy - cdiv3; 358 PetscCall(VecSetValues(G,1,&ind,&val,ADD_VALUES)); 359 } 360 } 361 } 362 363 /* Compute gradient contributions over the upper triangular elements */ 364 for (j=0; j<=ny; j++) { 365 for (i=0; i<=nx; i++) { 366 k = nx*j + i; 367 vb = zero; 368 vl = zero; 369 v = zero; 370 if (i < nx && j > 0) vb = x[k-nx]; 371 if (i > 0 && j < ny) vl = x[k-1]; 372 if (i < nx && j < ny) v = x[k]; 373 dvdx = (v-vl)/hx; 374 dvdy = (v-vb)/hy; 375 if (i != nx && j != 0) { 376 ind = k-nx; val = - dvdy/hy - cdiv3; 377 PetscCall(VecSetValues(G,1,&ind,&val,ADD_VALUES)); 378 } 379 if (i != 0 && j != ny) { 380 ind = k-1; val = - dvdx/hx - cdiv3; 381 PetscCall(VecSetValues(G,1,&ind,&val,ADD_VALUES)); 382 } 383 if (i != nx && j != ny) { 384 ind = k; val = dvdx/hx + dvdy/hy - cdiv3; 385 PetscCall(VecSetValues(G,1,&ind,&val,ADD_VALUES)); 386 } 387 } 388 } 389 PetscCall(VecRestoreArrayRead(X,&x)); 390 391 /* Assemble gradient vector */ 392 PetscCall(VecAssemblyBegin(G)); 393 PetscCall(VecAssemblyEnd(G)); 394 395 /* Scale the gradient */ 396 area = p5*hx*hy; 397 PetscCall(VecScale(G, area)); 398 PetscCall(PetscLogFlops(24.0*nx*ny)); 399 PetscFunctionReturn(0); 400 } 401 402 /* ------------------------------------------------------------------- */ 403 /* 404 FormHessian - Forms the Hessian matrix. 405 406 Input Parameters: 407 . tao - the Tao context 408 . X - the input vector 409 . ptr - optional user-defined context, as set by TaoSetHessian() 410 411 Output Parameters: 412 . H - Hessian matrix 413 . PrecH - optionally different preconditioning Hessian 414 . flag - flag indicating matrix structure 415 416 Notes: 417 This routine is intended simply as an example of the interface 418 to a Hessian evaluation routine. Since this example compute the 419 Hessian a column at a time, it is not particularly efficient and 420 is not recommended. 421 */ 422 PetscErrorCode FormHessian(Tao tao,Vec X,Mat H,Mat Hpre, void *ptr) 423 { 424 AppCtx *user = (AppCtx *) ptr; 425 PetscInt i,j, ndim = user->ndim; 426 PetscReal *y, zero = 0.0, one = 1.0; 427 PetscBool assembled; 428 429 PetscFunctionBeginUser; 430 user->xvec = X; 431 432 /* Initialize Hessian entries and work vector to zero */ 433 PetscCall(MatAssembled(H,&assembled)); 434 if (assembled)PetscCall(MatZeroEntries(H)); 435 436 PetscCall(VecSet(user->s, zero)); 437 438 /* Loop over matrix columns to compute entries of the Hessian */ 439 for (j=0; j<ndim; j++) { 440 PetscCall(VecSetValues(user->s,1,&j,&one,INSERT_VALUES)); 441 PetscCall(VecAssemblyBegin(user->s)); 442 PetscCall(VecAssemblyEnd(user->s)); 443 444 PetscCall(HessianProduct(ptr,user->s,user->y)); 445 446 PetscCall(VecSetValues(user->s,1,&j,&zero,INSERT_VALUES)); 447 PetscCall(VecAssemblyBegin(user->s)); 448 PetscCall(VecAssemblyEnd(user->s)); 449 450 PetscCall(VecGetArray(user->y,&y)); 451 for (i=0; i<ndim; i++) { 452 if (y[i] != zero) { 453 PetscCall(MatSetValues(H,1,&i,1,&j,&y[i],ADD_VALUES)); 454 } 455 } 456 PetscCall(VecRestoreArray(user->y,&y)); 457 } 458 PetscCall(MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY)); 459 PetscCall(MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY)); 460 PetscFunctionReturn(0); 461 } 462 463 /* ------------------------------------------------------------------- */ 464 /* 465 MatrixFreeHessian - Sets a pointer for use in computing Hessian-vector 466 products. 467 468 Input Parameters: 469 . tao - the Tao context 470 . X - the input vector 471 . ptr - optional user-defined context, as set by TaoSetHessian() 472 473 Output Parameters: 474 . H - Hessian matrix 475 . PrecH - optionally different preconditioning Hessian 476 . flag - flag indicating matrix structure 477 */ 478 PetscErrorCode MatrixFreeHessian(Tao tao,Vec X,Mat H,Mat PrecH, void *ptr) 479 { 480 AppCtx *user = (AppCtx *) ptr; 481 482 /* Sets location of vector for use in computing matrix-vector products of the form H(X)*y */ 483 PetscFunctionBeginUser; 484 user->xvec = X; 485 PetscFunctionReturn(0); 486 } 487 488 /* ------------------------------------------------------------------- */ 489 /* 490 HessianProductMat - Computes the matrix-vector product 491 y = mat*svec. 492 493 Input Parameters: 494 . mat - input matrix 495 . svec - input vector 496 497 Output Parameters: 498 . y - solution vector 499 */ 500 PetscErrorCode HessianProductMat(Mat mat,Vec svec,Vec y) 501 { 502 void *ptr; 503 504 PetscFunctionBeginUser; 505 PetscCall(MatShellGetContext(mat,&ptr)); 506 PetscCall(HessianProduct(ptr,svec,y)); 507 PetscFunctionReturn(0); 508 } 509 510 /* ------------------------------------------------------------------- */ 511 /* 512 Hessian Product - Computes the matrix-vector product: 513 y = f''(x)*svec. 514 515 Input Parameters: 516 . ptr - pointer to the user-defined context 517 . svec - input vector 518 519 Output Parameters: 520 . y - product vector 521 */ 522 PetscErrorCode HessianProduct(void *ptr,Vec svec,Vec y) 523 { 524 AppCtx *user = (AppCtx *)ptr; 525 PetscReal p5 = 0.5, zero = 0.0, one = 1.0, hx, hy, val, area; 526 const PetscScalar *x, *s; 527 PetscReal v, vb, vl, vr, vt, hxhx, hyhy; 528 PetscInt nx, ny, i, j, k, ind; 529 530 PetscFunctionBeginUser; 531 nx = user->mx; 532 ny = user->my; 533 hx = user->hx; 534 hy = user->hy; 535 hxhx = one/(hx*hx); 536 hyhy = one/(hy*hy); 537 538 /* Get pointers to vector data */ 539 PetscCall(VecGetArrayRead(user->xvec,&x)); 540 PetscCall(VecGetArrayRead(svec,&s)); 541 542 /* Initialize product vector to zero */ 543 PetscCall(VecSet(y, zero)); 544 545 /* Compute f''(x)*s over the lower triangular elements */ 546 for (j=-1; j<ny; j++) { 547 for (i=-1; i<nx; i++) { 548 k = nx*j + i; 549 v = zero; 550 vr = zero; 551 vt = zero; 552 if (i != -1 && j != -1) v = s[k]; 553 if (i != nx-1 && j != -1) { 554 vr = s[k+1]; 555 ind = k+1; val = hxhx*(vr-v); 556 PetscCall(VecSetValues(y,1,&ind,&val,ADD_VALUES)); 557 } 558 if (i != -1 && j != ny-1) { 559 vt = s[k+nx]; 560 ind = k+nx; val = hyhy*(vt-v); 561 PetscCall(VecSetValues(y,1,&ind,&val,ADD_VALUES)); 562 } 563 if (i != -1 && j != -1) { 564 ind = k; val = hxhx*(v-vr) + hyhy*(v-vt); 565 PetscCall(VecSetValues(y,1,&ind,&val,ADD_VALUES)); 566 } 567 } 568 } 569 570 /* Compute f''(x)*s over the upper triangular elements */ 571 for (j=0; j<=ny; j++) { 572 for (i=0; i<=nx; i++) { 573 k = nx*j + i; 574 v = zero; 575 vl = zero; 576 vb = zero; 577 if (i != nx && j != ny) v = s[k]; 578 if (i != nx && j != 0) { 579 vb = s[k-nx]; 580 ind = k-nx; val = hyhy*(vb-v); 581 PetscCall(VecSetValues(y,1,&ind,&val,ADD_VALUES)); 582 } 583 if (i != 0 && j != ny) { 584 vl = s[k-1]; 585 ind = k-1; val = hxhx*(vl-v); 586 PetscCall(VecSetValues(y,1,&ind,&val,ADD_VALUES)); 587 } 588 if (i != nx && j != ny) { 589 ind = k; val = hxhx*(v-vl) + hyhy*(v-vb); 590 PetscCall(VecSetValues(y,1,&ind,&val,ADD_VALUES)); 591 } 592 } 593 } 594 /* Restore vector data */ 595 PetscCall(VecRestoreArrayRead(svec,&s)); 596 PetscCall(VecRestoreArrayRead(user->xvec,&x)); 597 598 /* Assemble vector */ 599 PetscCall(VecAssemblyBegin(y)); 600 PetscCall(VecAssemblyEnd(y)); 601 602 /* Scale resulting vector by area */ 603 area = p5*hx*hy; 604 PetscCall(VecScale(y, area)); 605 PetscCall(PetscLogFlops(18.0*nx*ny)); 606 PetscFunctionReturn(0); 607 } 608 609 /*TEST 610 611 build: 612 requires: !complex 613 614 test: 615 suffix: 1 616 args: -tao_smonitor -tao_type ntl -tao_gatol 1.e-4 617 618 test: 619 suffix: 2 620 args: -tao_smonitor -tao_type ntr -tao_gatol 1.e-4 621 622 test: 623 suffix: 3 624 args: -tao_smonitor -tao_type bntr -tao_gatol 1.e-4 -my_tao_mf -tao_test_hessian 625 626 test: 627 suffix: 4 628 args: -tao_smonitor -tao_gatol 1e-3 -tao_type bqnls 629 630 test: 631 suffix: 5 632 args: -tao_smonitor -tao_gatol 1e-3 -tao_type blmvm 633 634 test: 635 suffix: 6 636 args: -tao_smonitor -tao_gatol 1e-3 -tao_type bqnktr -tao_bqnk_mat_type lmvmsr1 637 638 TEST*/ 639