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