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