1 /* Program usage: mpiexec -n <proc> minsurf2 [-help] [all TAO options] */ 2 3 /* 4 Include "petsctao.h" so we can use TAO solvers. 5 petscdmda.h for distributed array 6 */ 7 #include <petsctao.h> 8 #include <petscdmda.h> 9 10 static char help[] = 11 "This example demonstrates use of the TAO package to \n\ 12 solve an unconstrained minimization problem. This example is based on a \n\ 13 problem from the MINPACK-2 test suite. Given a rectangular 2-D domain and \n\ 14 boundary values along the edges of the domain, the objective is to find the\n\ 15 surface with the minimal area that satisfies the boundary conditions.\n\ 16 The command line options are:\n\ 17 -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction\n\ 18 -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction\n\ 19 -start <st>, where <st> =0 for zero vector, <st> >0 for random start, and <st> <0 \n\ 20 for an average of the boundary conditions\n\n"; 21 22 /*T 23 Concepts: TAO^Solving an unconstrained minimization problem 24 Routines: TaoCreate(); TaoSetType(); 25 Routines: TaoSetSolution(); 26 Routines: TaoSetObjectiveAndGradient(); 27 Routines: TaoSetHessian(); TaoSetFromOptions(); 28 Routines: TaoSetMonitor(); 29 Routines: TaoSolve(); TaoView(); 30 Routines: TaoDestroy(); 31 Processors: n 32 T*/ 33 34 /* 35 User-defined application context - contains data needed by the 36 application-provided call-back routines, FormFunctionGradient() 37 and FormHessian(). 38 */ 39 typedef struct { 40 PetscInt mx, my; /* discretization in x, y directions */ 41 PetscReal *bottom, *top, *left, *right; /* boundary values */ 42 DM dm; /* distributed array data structure */ 43 Mat H; /* Hessian */ 44 } AppCtx; 45 46 /* -------- User-defined Routines --------- */ 47 48 static PetscErrorCode MSA_BoundaryConditions(AppCtx*); 49 static PetscErrorCode MSA_InitialPoint(AppCtx*,Vec); 50 PetscErrorCode QuadraticH(AppCtx*,Vec,Mat); 51 PetscErrorCode FormFunctionGradient(Tao,Vec,PetscReal *,Vec,void*); 52 PetscErrorCode FormGradient(Tao,Vec,Vec,void*); 53 PetscErrorCode FormHessian(Tao,Vec,Mat,Mat,void*); 54 PetscErrorCode My_Monitor(Tao, void *); 55 56 int main(int argc, char **argv) 57 { 58 PetscInt Nx, Ny; /* number of processors in x- and y- directions */ 59 Vec x; /* solution, gradient vectors */ 60 PetscBool flg, viewmat; /* flags */ 61 PetscBool fddefault, fdcoloring; /* flags */ 62 Tao tao; /* TAO solver context */ 63 AppCtx user; /* user-defined work context */ 64 ISColoring iscoloring; 65 MatFDColoring matfdcoloring; 66 67 /* Initialize TAO */ 68 PetscCall(PetscInitialize(&argc, &argv,(char *)0,help)); 69 70 /* Specify dimension of the problem */ 71 user.mx = 10; user.my = 10; 72 73 /* Check for any command line arguments that override defaults */ 74 PetscCall(PetscOptionsGetInt(NULL,NULL,"-mx",&user.mx,&flg)); 75 PetscCall(PetscOptionsGetInt(NULL,NULL,"-my",&user.my,&flg)); 76 77 PetscCall(PetscPrintf(MPI_COMM_WORLD,"\n---- Minimum Surface Area Problem -----\n")); 78 PetscCall(PetscPrintf(MPI_COMM_WORLD,"mx: %D my: %D \n\n",user.mx,user.my)); 79 80 /* Let PETSc determine the vector distribution */ 81 Nx = PETSC_DECIDE; Ny = PETSC_DECIDE; 82 83 /* Create distributed array (DM) to manage parallel grid and vectors */ 84 PetscCall(DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,user.mx, user.my,Nx,Ny,1,1,NULL,NULL,&user.dm)); 85 PetscCall(DMSetFromOptions(user.dm)); 86 PetscCall(DMSetUp(user.dm)); 87 88 /* Create TAO solver and set desired solution method.*/ 89 PetscCall(TaoCreate(PETSC_COMM_WORLD,&tao)); 90 PetscCall(TaoSetType(tao,TAOCG)); 91 92 /* 93 Extract global vector from DA for the vector of variables -- PETSC routine 94 Compute the initial solution -- application specific, see below 95 Set this vector for use by TAO -- TAO routine 96 */ 97 PetscCall(DMCreateGlobalVector(user.dm,&x)); 98 PetscCall(MSA_BoundaryConditions(&user)); 99 PetscCall(MSA_InitialPoint(&user,x)); 100 PetscCall(TaoSetSolution(tao,x)); 101 102 /* 103 Initialize the Application context for use in function evaluations -- application specific, see below. 104 Set routines for function and gradient evaluation 105 */ 106 PetscCall(TaoSetObjectiveAndGradient(tao,NULL,FormFunctionGradient,(void *)&user)); 107 108 /* 109 Given the command line arguments, calculate the hessian with either the user- 110 provided function FormHessian, or the default finite-difference driven Hessian 111 functions 112 */ 113 PetscCall(PetscOptionsHasName(NULL,NULL,"-tao_fddefault",&fddefault)); 114 PetscCall(PetscOptionsHasName(NULL,NULL,"-tao_fdcoloring",&fdcoloring)); 115 116 /* 117 Create a matrix data structure to store the Hessian and set 118 the Hessian evalution routine. 119 Set the matrix structure to be used for Hessian evalutions 120 */ 121 PetscCall(DMCreateMatrix(user.dm,&user.H)); 122 PetscCall(MatSetOption(user.H,MAT_SYMMETRIC,PETSC_TRUE)); 123 124 if (fdcoloring) { 125 PetscCall(DMCreateColoring(user.dm,IS_COLORING_GLOBAL,&iscoloring)); 126 PetscCall(MatFDColoringCreate(user.H,iscoloring,&matfdcoloring)); 127 PetscCall(ISColoringDestroy(&iscoloring)); 128 PetscCall(MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))FormGradient,(void*)&user)); 129 PetscCall(MatFDColoringSetFromOptions(matfdcoloring)); 130 PetscCall(TaoSetHessian(tao,user.H,user.H,TaoDefaultComputeHessianColor,(void *)matfdcoloring)); 131 } else if (fddefault) { 132 PetscCall(TaoSetHessian(tao,user.H,user.H,TaoDefaultComputeHessian,(void *)NULL)); 133 } else { 134 PetscCall(TaoSetHessian(tao,user.H,user.H,FormHessian,(void *)&user)); 135 } 136 137 /* 138 If my_monitor option is in command line, then use the user-provided 139 monitoring function 140 */ 141 PetscCall(PetscOptionsHasName(NULL,NULL,"-my_monitor",&viewmat)); 142 if (viewmat) { 143 PetscCall(TaoSetMonitor(tao,My_Monitor,NULL,NULL)); 144 } 145 146 /* Check for any tao command line options */ 147 PetscCall(TaoSetFromOptions(tao)); 148 149 /* SOLVE THE APPLICATION */ 150 PetscCall(TaoSolve(tao)); 151 152 PetscCall(TaoView(tao,PETSC_VIEWER_STDOUT_WORLD)); 153 154 /* Free TAO data structures */ 155 PetscCall(TaoDestroy(&tao)); 156 157 /* Free PETSc data structures */ 158 PetscCall(VecDestroy(&x)); 159 PetscCall(MatDestroy(&user.H)); 160 if (fdcoloring) { 161 PetscCall(MatFDColoringDestroy(&matfdcoloring)); 162 } 163 PetscCall(PetscFree(user.bottom)); 164 PetscCall(PetscFree(user.top)); 165 PetscCall(PetscFree(user.left)); 166 PetscCall(PetscFree(user.right)); 167 PetscCall(DMDestroy(&user.dm)); 168 PetscCall(PetscFinalize()); 169 return 0; 170 } 171 172 PetscErrorCode FormGradient(Tao tao, Vec X, Vec G,void *userCtx) 173 { 174 PetscReal fcn; 175 176 PetscFunctionBegin; 177 PetscCall(FormFunctionGradient(tao,X,&fcn,G,userCtx)); 178 PetscFunctionReturn(0); 179 } 180 181 /* -------------------------------------------------------------------- */ 182 /* FormFunctionGradient - Evaluates the function and corresponding gradient. 183 184 Input Parameters: 185 . tao - the Tao context 186 . XX - input vector 187 . userCtx - optional user-defined context, as set by TaoSetObjectiveAndGradient() 188 189 Output Parameters: 190 . fcn - the newly evaluated function 191 . GG - vector containing the newly evaluated gradient 192 */ 193 PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *fcn,Vec G,void *userCtx) 194 { 195 AppCtx *user = (AppCtx *) userCtx; 196 PetscInt i,j; 197 PetscInt mx=user->mx, my=user->my; 198 PetscInt xs,xm,gxs,gxm,ys,ym,gys,gym; 199 PetscReal ft=0.0; 200 PetscReal hx=1.0/(mx+1),hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy, area=0.5*hx*hy; 201 PetscReal rhx=mx+1, rhy=my+1; 202 PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb; 203 PetscReal df1dxc,df2dxc,df3dxc,df4dxc,df5dxc,df6dxc; 204 PetscReal **g, **x; 205 Vec localX; 206 207 PetscFunctionBegin; 208 /* Get local mesh boundaries */ 209 PetscCall(DMGetLocalVector(user->dm,&localX)); 210 PetscCall(DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL)); 211 PetscCall(DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL)); 212 213 /* Scatter ghost points to local vector */ 214 PetscCall(DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX)); 215 PetscCall(DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX)); 216 217 /* Get pointers to vector data */ 218 PetscCall(DMDAVecGetArray(user->dm,localX,(void**)&x)); 219 PetscCall(DMDAVecGetArray(user->dm,G,(void**)&g)); 220 221 /* Compute function and gradient over the locally owned part of the mesh */ 222 for (j=ys; j<ys+ym; j++) { 223 for (i=xs; i< xs+xm; i++) { 224 225 xc = x[j][i]; 226 xlt=xrb=xl=xr=xb=xt=xc; 227 228 if (i==0) { /* left side */ 229 xl= user->left[j-ys+1]; 230 xlt = user->left[j-ys+2]; 231 } else { 232 xl = x[j][i-1]; 233 } 234 235 if (j==0) { /* bottom side */ 236 xb=user->bottom[i-xs+1]; 237 xrb =user->bottom[i-xs+2]; 238 } else { 239 xb = x[j-1][i]; 240 } 241 242 if (i+1 == gxs+gxm) { /* right side */ 243 xr=user->right[j-ys+1]; 244 xrb = user->right[j-ys]; 245 } else { 246 xr = x[j][i+1]; 247 } 248 249 if (j+1==gys+gym) { /* top side */ 250 xt=user->top[i-xs+1]; 251 xlt = user->top[i-xs]; 252 }else { 253 xt = x[j+1][i]; 254 } 255 256 if (i>gxs && j+1<gys+gym) { 257 xlt = x[j+1][i-1]; 258 } 259 if (j>gys && i+1<gxs+gxm) { 260 xrb = x[j-1][i+1]; 261 } 262 263 d1 = (xc-xl); 264 d2 = (xc-xr); 265 d3 = (xc-xt); 266 d4 = (xc-xb); 267 d5 = (xr-xrb); 268 d6 = (xrb-xb); 269 d7 = (xlt-xl); 270 d8 = (xt-xlt); 271 272 df1dxc = d1*hydhx; 273 df2dxc = (d1*hydhx + d4*hxdhy); 274 df3dxc = d3*hxdhy; 275 df4dxc = (d2*hydhx + d3*hxdhy); 276 df5dxc = d2*hydhx; 277 df6dxc = d4*hxdhy; 278 279 d1 *= rhx; 280 d2 *= rhx; 281 d3 *= rhy; 282 d4 *= rhy; 283 d5 *= rhy; 284 d6 *= rhx; 285 d7 *= rhy; 286 d8 *= rhx; 287 288 f1 = PetscSqrtReal(1.0 + d1*d1 + d7*d7); 289 f2 = PetscSqrtReal(1.0 + d1*d1 + d4*d4); 290 f3 = PetscSqrtReal(1.0 + d3*d3 + d8*d8); 291 f4 = PetscSqrtReal(1.0 + d3*d3 + d2*d2); 292 f5 = PetscSqrtReal(1.0 + d2*d2 + d5*d5); 293 f6 = PetscSqrtReal(1.0 + d4*d4 + d6*d6); 294 295 ft = ft + (f2 + f4); 296 297 df1dxc /= f1; 298 df2dxc /= f2; 299 df3dxc /= f3; 300 df4dxc /= f4; 301 df5dxc /= f5; 302 df6dxc /= f6; 303 304 g[j][i] = (df1dxc+df2dxc+df3dxc+df4dxc+df5dxc+df6dxc) * 0.5; 305 306 } 307 } 308 309 /* Compute triangular areas along the border of the domain. */ 310 if (xs==0) { /* left side */ 311 for (j=ys; j<ys+ym; j++) { 312 d3=(user->left[j-ys+1] - user->left[j-ys+2])*rhy; 313 d2=(user->left[j-ys+1] - x[j][0]) *rhx; 314 ft = ft+PetscSqrtReal(1.0 + d3*d3 + d2*d2); 315 } 316 } 317 if (ys==0) { /* bottom side */ 318 for (i=xs; i<xs+xm; i++) { 319 d2=(user->bottom[i+1-xs]-user->bottom[i-xs+2])*rhx; 320 d3=(user->bottom[i-xs+1]-x[0][i])*rhy; 321 ft = ft+PetscSqrtReal(1.0 + d3*d3 + d2*d2); 322 } 323 } 324 325 if (xs+xm==mx) { /* right side */ 326 for (j=ys; j< ys+ym; j++) { 327 d1=(x[j][mx-1] - user->right[j-ys+1])*rhx; 328 d4=(user->right[j-ys]-user->right[j-ys+1])*rhy; 329 ft = ft+PetscSqrtReal(1.0 + d1*d1 + d4*d4); 330 } 331 } 332 if (ys+ym==my) { /* top side */ 333 for (i=xs; i<xs+xm; i++) { 334 d1=(x[my-1][i] - user->top[i-xs+1])*rhy; 335 d4=(user->top[i-xs+1] - user->top[i-xs])*rhx; 336 ft = ft+PetscSqrtReal(1.0 + d1*d1 + d4*d4); 337 } 338 } 339 340 if (ys==0 && xs==0) { 341 d1=(user->left[0]-user->left[1])/hy; 342 d2=(user->bottom[0]-user->bottom[1])*rhx; 343 ft +=PetscSqrtReal(1.0 + d1*d1 + d2*d2); 344 } 345 if (ys+ym == my && xs+xm == mx) { 346 d1=(user->right[ym+1] - user->right[ym])*rhy; 347 d2=(user->top[xm+1] - user->top[xm])*rhx; 348 ft +=PetscSqrtReal(1.0 + d1*d1 + d2*d2); 349 } 350 351 ft=ft*area; 352 PetscCallMPI(MPI_Allreduce(&ft,fcn,1,MPIU_REAL,MPIU_SUM,MPI_COMM_WORLD)); 353 354 /* Restore vectors */ 355 PetscCall(DMDAVecRestoreArray(user->dm,localX,(void**)&x)); 356 PetscCall(DMDAVecRestoreArray(user->dm,G,(void**)&g)); 357 358 /* Scatter values to global vector */ 359 PetscCall(DMRestoreLocalVector(user->dm,&localX)); 360 PetscCall(PetscLogFlops(67.0*xm*ym)); 361 PetscFunctionReturn(0); 362 } 363 364 /* ------------------------------------------------------------------- */ 365 /* 366 FormHessian - Evaluates Hessian matrix. 367 368 Input Parameters: 369 . tao - the Tao context 370 . x - input vector 371 . ptr - optional user-defined context, as set by TaoSetHessian() 372 373 Output Parameters: 374 . H - Hessian matrix 375 . Hpre - optionally different preconditioning matrix 376 . flg - flag indicating matrix structure 377 378 */ 379 PetscErrorCode FormHessian(Tao tao,Vec X,Mat H, Mat Hpre, void *ptr) 380 { 381 AppCtx *user = (AppCtx *) ptr; 382 383 PetscFunctionBegin; 384 /* Evaluate the Hessian entries*/ 385 PetscCall(QuadraticH(user,X,H)); 386 PetscFunctionReturn(0); 387 } 388 389 /* ------------------------------------------------------------------- */ 390 /* 391 QuadraticH - Evaluates Hessian matrix. 392 393 Input Parameters: 394 . user - user-defined context, as set by TaoSetHessian() 395 . X - input vector 396 397 Output Parameter: 398 . H - Hessian matrix 399 */ 400 PetscErrorCode QuadraticH(AppCtx *user, Vec X, Mat Hessian) 401 { 402 PetscInt i,j,k; 403 PetscInt mx=user->mx, my=user->my; 404 PetscInt xs,xm,gxs,gxm,ys,ym,gys,gym; 405 PetscReal hx=1.0/(mx+1), hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy; 406 PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb; 407 PetscReal hl,hr,ht,hb,hc,htl,hbr; 408 PetscReal **x, v[7]; 409 MatStencil col[7],row; 410 Vec localX; 411 PetscBool assembled; 412 413 PetscFunctionBegin; 414 /* Get local mesh boundaries */ 415 PetscCall(DMGetLocalVector(user->dm,&localX)); 416 417 PetscCall(DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL)); 418 PetscCall(DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL)); 419 420 /* Scatter ghost points to local vector */ 421 PetscCall(DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX)); 422 PetscCall(DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX)); 423 424 /* Get pointers to vector data */ 425 PetscCall(DMDAVecGetArray(user->dm,localX,(void**)&x)); 426 427 /* Initialize matrix entries to zero */ 428 PetscCall(MatAssembled(Hessian,&assembled)); 429 if (assembled) PetscCall(MatZeroEntries(Hessian)); 430 431 /* Set various matrix options */ 432 PetscCall(MatSetOption(Hessian,MAT_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE)); 433 434 /* Compute Hessian over the locally owned part of the mesh */ 435 436 for (j=ys; j<ys+ym; j++) { 437 438 for (i=xs; i< xs+xm; i++) { 439 440 xc = x[j][i]; 441 xlt=xrb=xl=xr=xb=xt=xc; 442 443 /* Left side */ 444 if (i==0) { 445 xl = user->left[j-ys+1]; 446 xlt = user->left[j-ys+2]; 447 } else { 448 xl = x[j][i-1]; 449 } 450 451 if (j==0) { 452 xb = user->bottom[i-xs+1]; 453 xrb = user->bottom[i-xs+2]; 454 } else { 455 xb = x[j-1][i]; 456 } 457 458 if (i+1 == mx) { 459 xr = user->right[j-ys+1]; 460 xrb = user->right[j-ys]; 461 } else { 462 xr = x[j][i+1]; 463 } 464 465 if (j+1==my) { 466 xt = user->top[i-xs+1]; 467 xlt = user->top[i-xs]; 468 }else { 469 xt = x[j+1][i]; 470 } 471 472 if (i>0 && j+1<my) { 473 xlt = x[j+1][i-1]; 474 } 475 if (j>0 && i+1<mx) { 476 xrb = x[j-1][i+1]; 477 } 478 479 d1 = (xc-xl)/hx; 480 d2 = (xc-xr)/hx; 481 d3 = (xc-xt)/hy; 482 d4 = (xc-xb)/hy; 483 d5 = (xrb-xr)/hy; 484 d6 = (xrb-xb)/hx; 485 d7 = (xlt-xl)/hy; 486 d8 = (xlt-xt)/hx; 487 488 f1 = PetscSqrtReal(1.0 + d1*d1 + d7*d7); 489 f2 = PetscSqrtReal(1.0 + d1*d1 + d4*d4); 490 f3 = PetscSqrtReal(1.0 + d3*d3 + d8*d8); 491 f4 = PetscSqrtReal(1.0 + d3*d3 + d2*d2); 492 f5 = PetscSqrtReal(1.0 + d2*d2 + d5*d5); 493 f6 = PetscSqrtReal(1.0 + d4*d4 + d6*d6); 494 495 hl = (-hydhx*(1.0+d7*d7)+d1*d7)/(f1*f1*f1)+ 496 (-hydhx*(1.0+d4*d4)+d1*d4)/(f2*f2*f2); 497 hr = (-hydhx*(1.0+d5*d5)+d2*d5)/(f5*f5*f5)+ 498 (-hydhx*(1.0+d3*d3)+d2*d3)/(f4*f4*f4); 499 ht = (-hxdhy*(1.0+d8*d8)+d3*d8)/(f3*f3*f3)+ 500 (-hxdhy*(1.0+d2*d2)+d2*d3)/(f4*f4*f4); 501 hb = (-hxdhy*(1.0+d6*d6)+d4*d6)/(f6*f6*f6)+ 502 (-hxdhy*(1.0+d1*d1)+d1*d4)/(f2*f2*f2); 503 504 hbr = -d2*d5/(f5*f5*f5) - d4*d6/(f6*f6*f6); 505 htl = -d1*d7/(f1*f1*f1) - d3*d8/(f3*f3*f3); 506 507 hc = hydhx*(1.0+d7*d7)/(f1*f1*f1) + hxdhy*(1.0+d8*d8)/(f3*f3*f3) + 508 hydhx*(1.0+d5*d5)/(f5*f5*f5) + hxdhy*(1.0+d6*d6)/(f6*f6*f6) + 509 (hxdhy*(1.0+d1*d1)+hydhx*(1.0+d4*d4)-2*d1*d4)/(f2*f2*f2) + 510 (hxdhy*(1.0+d2*d2)+hydhx*(1.0+d3*d3)-2*d2*d3)/(f4*f4*f4); 511 512 hl/=2.0; hr/=2.0; ht/=2.0; hb/=2.0; hbr/=2.0; htl/=2.0; hc/=2.0; 513 514 row.j = j; row.i = i; 515 k=0; 516 if (j>0) { 517 v[k]=hb; 518 col[k].j = j - 1; col[k].i = i; 519 k++; 520 } 521 522 if (j>0 && i < mx -1) { 523 v[k]=hbr; 524 col[k].j = j - 1; col[k].i = i+1; 525 k++; 526 } 527 528 if (i>0) { 529 v[k]= hl; 530 col[k].j = j; col[k].i = i-1; 531 k++; 532 } 533 534 v[k]= hc; 535 col[k].j = j; col[k].i = i; 536 k++; 537 538 if (i < mx-1) { 539 v[k]= hr; 540 col[k].j = j; col[k].i = i+1; 541 k++; 542 } 543 544 if (i>0 && j < my-1) { 545 v[k]= htl; 546 col[k].j = j+1; col[k].i = i-1; 547 k++; 548 } 549 550 if (j < my-1) { 551 v[k]= ht; 552 col[k].j = j+1; col[k].i = i; 553 k++; 554 } 555 556 /* 557 Set matrix values using local numbering, which was defined 558 earlier, in the main routine. 559 */ 560 PetscCall(MatSetValuesStencil(Hessian,1,&row,k,col,v,INSERT_VALUES)); 561 } 562 } 563 564 PetscCall(DMDAVecRestoreArray(user->dm,localX,(void**)&x)); 565 PetscCall(DMRestoreLocalVector(user->dm,&localX)); 566 567 PetscCall(MatAssemblyBegin(Hessian,MAT_FINAL_ASSEMBLY)); 568 PetscCall(MatAssemblyEnd(Hessian,MAT_FINAL_ASSEMBLY)); 569 570 PetscCall(PetscLogFlops(199.0*xm*ym)); 571 PetscFunctionReturn(0); 572 } 573 574 /* ------------------------------------------------------------------- */ 575 /* 576 MSA_BoundaryConditions - Calculates the boundary conditions for 577 the region. 578 579 Input Parameter: 580 . user - user-defined application context 581 582 Output Parameter: 583 . user - user-defined application context 584 */ 585 static PetscErrorCode MSA_BoundaryConditions(AppCtx * user) 586 { 587 PetscInt i,j,k,limit=0,maxits=5; 588 PetscInt xs,ys,xm,ym,gxs,gys,gxm,gym; 589 PetscInt mx=user->mx,my=user->my; 590 PetscInt bsize=0, lsize=0, tsize=0, rsize=0; 591 PetscReal one=1.0, two=2.0, three=3.0, tol=1e-10; 592 PetscReal fnorm,det,hx,hy,xt=0,yt=0; 593 PetscReal u1,u2,nf1,nf2,njac11,njac12,njac21,njac22; 594 PetscReal b=-0.5, t=0.5, l=-0.5, r=0.5; 595 PetscReal *boundary; 596 PetscBool flg; 597 598 PetscFunctionBegin; 599 /* Get local mesh boundaries */ 600 PetscCall(DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL)); 601 PetscCall(DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL)); 602 603 bsize=xm+2; 604 lsize=ym+2; 605 rsize=ym+2; 606 tsize=xm+2; 607 608 PetscCall(PetscMalloc1(bsize,&user->bottom)); 609 PetscCall(PetscMalloc1(tsize,&user->top)); 610 PetscCall(PetscMalloc1(lsize,&user->left)); 611 PetscCall(PetscMalloc1(rsize,&user->right)); 612 613 hx= (r-l)/(mx+1); hy=(t-b)/(my+1); 614 615 for (j=0; j<4; j++) { 616 if (j==0) { 617 yt=b; 618 xt=l+hx*xs; 619 limit=bsize; 620 boundary=user->bottom; 621 } else if (j==1) { 622 yt=t; 623 xt=l+hx*xs; 624 limit=tsize; 625 boundary=user->top; 626 } else if (j==2) { 627 yt=b+hy*ys; 628 xt=l; 629 limit=lsize; 630 boundary=user->left; 631 } else { /* if (j==3) */ 632 yt=b+hy*ys; 633 xt=r; 634 limit=rsize; 635 boundary=user->right; 636 } 637 638 for (i=0; i<limit; i++) { 639 u1=xt; 640 u2=-yt; 641 for (k=0; k<maxits; k++) { 642 nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt; 643 nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt; 644 fnorm=PetscSqrtReal(nf1*nf1+nf2*nf2); 645 if (fnorm <= tol) break; 646 njac11=one+u2*u2-u1*u1; 647 njac12=two*u1*u2; 648 njac21=-two*u1*u2; 649 njac22=-one - u1*u1 + u2*u2; 650 det = njac11*njac22-njac21*njac12; 651 u1 = u1-(njac22*nf1-njac12*nf2)/det; 652 u2 = u2-(njac11*nf2-njac21*nf1)/det; 653 } 654 655 boundary[i]=u1*u1-u2*u2; 656 if (j==0 || j==1) { 657 xt=xt+hx; 658 } else { /* if (j==2 || j==3) */ 659 yt=yt+hy; 660 } 661 662 } 663 664 } 665 666 /* Scale the boundary if desired */ 667 if (1==1) { 668 PetscReal scl = 1.0; 669 670 PetscCall(PetscOptionsGetReal(NULL,NULL,"-bottom",&scl,&flg)); 671 if (flg) { 672 for (i=0;i<bsize;i++) user->bottom[i]*=scl; 673 } 674 675 PetscCall(PetscOptionsGetReal(NULL,NULL,"-top",&scl,&flg)); 676 if (flg) { 677 for (i=0;i<tsize;i++) user->top[i]*=scl; 678 } 679 680 PetscCall(PetscOptionsGetReal(NULL,NULL,"-right",&scl,&flg)); 681 if (flg) { 682 for (i=0;i<rsize;i++) user->right[i]*=scl; 683 } 684 685 PetscCall(PetscOptionsGetReal(NULL,NULL,"-left",&scl,&flg)); 686 if (flg) { 687 for (i=0;i<lsize;i++) user->left[i]*=scl; 688 } 689 } 690 PetscFunctionReturn(0); 691 } 692 693 /* ------------------------------------------------------------------- */ 694 /* 695 MSA_InitialPoint - Calculates the initial guess in one of three ways. 696 697 Input Parameters: 698 . user - user-defined application context 699 . X - vector for initial guess 700 701 Output Parameters: 702 . X - newly computed initial guess 703 */ 704 static PetscErrorCode MSA_InitialPoint(AppCtx * user, Vec X) 705 { 706 PetscInt start2=-1,i,j; 707 PetscReal start1=0; 708 PetscBool flg1,flg2; 709 710 PetscFunctionBegin; 711 PetscCall(PetscOptionsGetReal(NULL,NULL,"-start",&start1,&flg1)); 712 PetscCall(PetscOptionsGetInt(NULL,NULL,"-random",&start2,&flg2)); 713 714 if (flg1) { /* The zero vector is reasonable */ 715 716 PetscCall(VecSet(X, start1)); 717 718 } else if (flg2 && start2>0) { /* Try a random start between -0.5 and 0.5 */ 719 PetscRandom rctx; PetscReal np5=-0.5; 720 721 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD,&rctx)); 722 PetscCall(VecSetRandom(X, rctx)); 723 PetscCall(PetscRandomDestroy(&rctx)); 724 PetscCall(VecShift(X, np5)); 725 726 } else { /* Take an average of the boundary conditions */ 727 PetscInt xs,xm,ys,ym; 728 PetscInt mx=user->mx,my=user->my; 729 PetscReal **x; 730 731 /* Get local mesh boundaries */ 732 PetscCall(DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL)); 733 734 /* Get pointers to vector data */ 735 PetscCall(DMDAVecGetArray(user->dm,X,(void**)&x)); 736 737 /* Perform local computations */ 738 for (j=ys; j<ys+ym; j++) { 739 for (i=xs; i< xs+xm; i++) { 740 x[j][i] = (((j+1)*user->bottom[i-xs+1]+(my-j+1)*user->top[i-xs+1])/(my+2)+((i+1)*user->left[j-ys+1]+(mx-i+1)*user->right[j-ys+1])/(mx+2))/2.0; 741 } 742 } 743 PetscCall(DMDAVecRestoreArray(user->dm,X,(void**)&x)); 744 PetscCall(PetscLogFlops(9.0*xm*ym)); 745 } 746 PetscFunctionReturn(0); 747 } 748 749 /*-----------------------------------------------------------------------*/ 750 PetscErrorCode My_Monitor(Tao tao, void *ctx) 751 { 752 Vec X; 753 754 PetscFunctionBegin; 755 PetscCall(TaoGetSolution(tao,&X)); 756 PetscCall(VecView(X,PETSC_VIEWER_STDOUT_WORLD)); 757 PetscFunctionReturn(0); 758 } 759 760 /*TEST 761 762 build: 763 requires: !complex 764 765 test: 766 args: -tao_smonitor -tao_type lmvm -mx 10 -my 8 -tao_gatol 1.e-3 767 requires: !single 768 769 test: 770 suffix: 2 771 nsize: 2 772 args: -tao_smonitor -tao_type nls -tao_nls_ksp_max_it 15 -tao_gatol 1.e-4 773 filter: grep -v "nls ksp" 774 requires: !single 775 776 test: 777 suffix: 3 778 nsize: 3 779 args: -tao_smonitor -tao_type cg -tao_cg_type fr -mx 10 -my 10 -tao_gatol 1.e-3 780 requires: !single 781 782 test: 783 suffix: 5 784 nsize: 2 785 args: -tao_smonitor -tao_type bmrm -mx 10 -my 8 -tao_gatol 1.e-3 786 requires: !single 787 788 TEST*/ 789