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