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