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