1 /* Program usage: mpiexec -n 1 minsurf1 [-help] [all TAO options] */ 2 3 /* Include "petsctao.h" so we can use TAO solvers. */ 4 #include <petsctao.h> 5 6 static char help[] = 7 "This example demonstrates use of the TAO package to \n\ 8 solve an unconstrained minimization problem. This example is based on a \n\ 9 problem from the MINPACK-2 test suite. Given a rectangular 2-D domain and \n\ 10 boundary values along the edges of the domain, the objective is to find the\n\ 11 surface with the minimal area that satisfies the boundary conditions.\n\ 12 The command line options are:\n\ 13 -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction\n\ 14 -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction\n\ 15 -start <st>, where <st> =0 for zero vector, and an average of the boundary conditions otherwise \n\n"; 16 17 /*T 18 Concepts: TAO^Solving an unconstrained minimization problem 19 Routines: TaoCreate(); TaoSetType(); 20 Routines: TaoSetInitialVector(); 21 Routines: TaoSetObjectiveAndGradientRoutine(); 22 Routines: TaoSetHessianRoutine(); TaoSetFromOptions(); 23 Routines: TaoGetKSP(); TaoSolve(); 24 Routines: TaoDestroy(); 25 Processors: 1 26 T*/ 27 28 /* 29 User-defined application context - contains data needed by the 30 application-provided call-back routines, FormFunctionGradient() 31 and FormHessian(). 32 */ 33 typedef struct { 34 PetscInt mx, my; /* discretization in x, y directions */ 35 PetscReal *bottom, *top, *left, *right; /* boundary values */ 36 Mat H; 37 } AppCtx; 38 39 /* -------- User-defined Routines --------- */ 40 41 static PetscErrorCode MSA_BoundaryConditions(AppCtx*); 42 static PetscErrorCode MSA_InitialPoint(AppCtx*,Vec); 43 static PetscErrorCode QuadraticH(AppCtx*,Vec,Mat); 44 PetscErrorCode FormFunctionGradient(Tao,Vec,PetscReal*,Vec,void*); 45 PetscErrorCode FormHessian(Tao,Vec,Mat,Mat,void*); 46 47 int main(int argc, char **argv) 48 { 49 PetscErrorCode ierr; /* used to check for functions returning nonzeros */ 50 PetscInt N; /* Size of vector */ 51 PetscMPIInt size; /* Number of processors */ 52 Vec x; /* solution, gradient vectors */ 53 KSP ksp; /* PETSc Krylov subspace method */ 54 PetscBool flg; /* A return value when checking for user options */ 55 Tao tao; /* Tao solver context */ 56 AppCtx user; /* user-defined work context */ 57 58 /* Initialize TAO,PETSc */ 59 ierr = PetscInitialize(&argc, &argv,(char *)0,help);if (ierr) return ierr; 60 61 ierr = MPI_Comm_size(MPI_COMM_WORLD,&size);CHKERRMPI(ierr); 62 if (size >1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"Incorrect number of processors"); 63 64 /* Specify default dimension of the problem */ 65 user.mx = 4; user.my = 4; 66 67 /* Check for any command line arguments that override defaults */ 68 ierr = PetscOptionsGetInt(NULL,NULL,"-mx",&user.mx,&flg);CHKERRQ(ierr); 69 ierr = PetscOptionsGetInt(NULL,NULL,"-my",&user.my,&flg);CHKERRQ(ierr); 70 71 ierr = PetscPrintf(PETSC_COMM_SELF,"\n---- Minimum Surface Area Problem -----\n");CHKERRQ(ierr); 72 ierr = PetscPrintf(PETSC_COMM_SELF,"mx: %D my: %D \n\n",user.mx,user.my);CHKERRQ(ierr); 73 74 /* Calculate any derived values from parameters */ 75 N = user.mx*user.my; 76 77 /* Create TAO solver and set desired solution method */ 78 ierr = TaoCreate(PETSC_COMM_SELF,&tao);CHKERRQ(ierr); 79 ierr = TaoSetType(tao,TAOLMVM);CHKERRQ(ierr); 80 81 /* Initialize minsurf application data structure for use in the function evaluations */ 82 ierr = MSA_BoundaryConditions(&user);CHKERRQ(ierr); /* Application specific routine */ 83 84 /* 85 Create a vector to hold the variables. Compute an initial solution. 86 Set this vector, which will be used by TAO. 87 */ 88 ierr = VecCreateSeq(PETSC_COMM_SELF,N,&x);CHKERRQ(ierr); 89 ierr = MSA_InitialPoint(&user,x);CHKERRQ(ierr); /* Application specific routine */ 90 ierr = TaoSetInitialVector(tao,x);CHKERRQ(ierr); /* A TAO routine */ 91 92 /* Provide TAO routines for function, gradient, and Hessian evaluation */ 93 ierr = TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void *)&user);CHKERRQ(ierr); 94 95 /* Create a matrix data structure to store the Hessian. This structure will be used by TAO */ 96 ierr = MatCreateSeqAIJ(PETSC_COMM_SELF,N,N,7,NULL,&(user.H));CHKERRQ(ierr); 97 ierr = MatSetOption(user.H,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr); 98 ierr = TaoSetHessianRoutine(tao,user.H,user.H,FormHessian,(void *)&user);CHKERRQ(ierr); 99 100 /* Check for any TAO command line options */ 101 ierr = TaoSetFromOptions(tao);CHKERRQ(ierr); 102 103 /* Limit the number of iterations in the KSP linear solver */ 104 ierr = TaoGetKSP(tao,&ksp);CHKERRQ(ierr); 105 if (ksp) { 106 ierr = KSPSetTolerances(ksp,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,user.mx*user.my);CHKERRQ(ierr); 107 } 108 109 /* SOLVE THE APPLICATION */ 110 ierr = TaoSolve(tao);CHKERRQ(ierr); 111 112 ierr = TaoDestroy(&tao);CHKERRQ(ierr); 113 ierr = VecDestroy(&x);CHKERRQ(ierr); 114 ierr = MatDestroy(&user.H);CHKERRQ(ierr); 115 ierr = PetscFree(user.bottom);CHKERRQ(ierr); 116 ierr = PetscFree(user.top);CHKERRQ(ierr); 117 ierr = PetscFree(user.left);CHKERRQ(ierr); 118 ierr = PetscFree(user.right);CHKERRQ(ierr); 119 120 ierr = PetscFinalize(); 121 return ierr; 122 } 123 124 /* -------------------------------------------------------------------- */ 125 126 /* FormFunctionGradient - Evaluates function and corresponding gradient. 127 128 Input Parameters: 129 . tao - the Tao context 130 . X - input vector 131 . userCtx - optional user-defined context, as set by TaoSetFunctionGradient() 132 133 Output Parameters: 134 . fcn - the newly evaluated function 135 . G - vector containing the newly evaluated gradient 136 */ 137 PetscErrorCode FormFunctionGradient(Tao tao,Vec X,PetscReal *fcn,Vec G,void *userCtx) 138 { 139 AppCtx *user = (AppCtx *) userCtx; 140 PetscErrorCode ierr; 141 PetscInt i,j,row; 142 PetscInt mx=user->mx, my=user->my; 143 PetscReal rhx=mx+1, rhy=my+1; 144 PetscReal hx=1.0/(mx+1),hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy, area=0.5*hx*hy, ft=0; 145 PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb; 146 PetscReal df1dxc,df2dxc,df3dxc,df4dxc,df5dxc,df6dxc; 147 PetscReal zero=0.0; 148 PetscScalar *g; 149 const PetscScalar *x; 150 151 PetscFunctionBeginUser; 152 ierr = VecSet(G, zero);CHKERRQ(ierr); 153 154 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); 155 ierr = VecGetArray(G,&g);CHKERRQ(ierr); 156 157 /* Compute function over the locally owned part of the mesh */ 158 for (j=0; j<my; j++) { 159 for (i=0; i< mx; i++) { 160 row=(j)*mx + (i); 161 xc = x[row]; 162 xlt=xrb=xl=xr=xb=xt=xc; 163 if (i==0) { /* left side */ 164 xl = user->left[j+1]; 165 xlt = user->left[j+2]; 166 } else { 167 xl = x[row-1]; 168 } 169 170 if (j==0) { /* bottom side */ 171 xb = user->bottom[i+1]; 172 xrb = user->bottom[i+2]; 173 } else { 174 xb = x[row-mx]; 175 } 176 177 if (i+1 == mx) { /* right side */ 178 xr = user->right[j+1]; 179 xrb = user->right[j]; 180 } else { 181 xr = x[row+1]; 182 } 183 184 if (j+1==0+my) { /* top side */ 185 xt = user->top[i+1]; 186 xlt = user->top[i]; 187 }else { 188 xt = x[row+mx]; 189 } 190 191 if (i>0 && j+1<my) { 192 xlt = x[row-1+mx]; 193 } 194 if (j>0 && i+1<mx) { 195 xrb = x[row+1-mx]; 196 } 197 198 d1 = (xc-xl); 199 d2 = (xc-xr); 200 d3 = (xc-xt); 201 d4 = (xc-xb); 202 d5 = (xr-xrb); 203 d6 = (xrb-xb); 204 d7 = (xlt-xl); 205 d8 = (xt-xlt); 206 207 df1dxc = d1*hydhx; 208 df2dxc = (d1*hydhx + d4*hxdhy); 209 df3dxc = d3*hxdhy; 210 df4dxc = (d2*hydhx + d3*hxdhy); 211 df5dxc = d2*hydhx; 212 df6dxc = d4*hxdhy; 213 214 d1 *= rhx; 215 d2 *= rhx; 216 d3 *= rhy; 217 d4 *= rhy; 218 d5 *= rhy; 219 d6 *= rhx; 220 d7 *= rhy; 221 d8 *= rhx; 222 223 f1 = PetscSqrtScalar(1.0 + d1*d1 + d7*d7); 224 f2 = PetscSqrtScalar(1.0 + d1*d1 + d4*d4); 225 f3 = PetscSqrtScalar(1.0 + d3*d3 + d8*d8); 226 f4 = PetscSqrtScalar(1.0 + d3*d3 + d2*d2); 227 f5 = PetscSqrtScalar(1.0 + d2*d2 + d5*d5); 228 f6 = PetscSqrtScalar(1.0 + d4*d4 + d6*d6); 229 230 ft = ft + (f2 + f4); 231 232 df1dxc /= f1; 233 df2dxc /= f2; 234 df3dxc /= f3; 235 df4dxc /= f4; 236 df5dxc /= f5; 237 df6dxc /= f6; 238 239 g[row] = (df1dxc+df2dxc+df3dxc+df4dxc+df5dxc+df6dxc)/2.0; 240 } 241 } 242 243 for (j=0; j<my; j++) { /* left side */ 244 d3 = (user->left[j+1] - user->left[j+2])*rhy; 245 d2 = (user->left[j+1] - x[j*mx])*rhx; 246 ft = ft+PetscSqrtScalar(1.0 + d3*d3 + d2*d2); 247 } 248 249 for (i=0; i<mx; i++) { /* bottom */ 250 d2 = (user->bottom[i+1]-user->bottom[i+2])*rhx; 251 d3 = (user->bottom[i+1]-x[i])*rhy; 252 ft = ft+PetscSqrtScalar(1.0 + d3*d3 + d2*d2); 253 } 254 255 for (j=0; j< my; j++) { /* right side */ 256 d1 = (x[(j+1)*mx-1]-user->right[j+1])*rhx; 257 d4 = (user->right[j]-user->right[j+1])*rhy; 258 ft = ft+PetscSqrtScalar(1.0 + d1*d1 + d4*d4); 259 } 260 261 for (i=0; i<mx; i++) { /* top side */ 262 d1 = (x[(my-1)*mx + i] - user->top[i+1])*rhy; 263 d4 = (user->top[i+1] - user->top[i])*rhx; 264 ft = ft+PetscSqrtScalar(1.0 + d1*d1 + d4*d4); 265 } 266 267 /* Bottom left corner */ 268 d1 = (user->left[0]-user->left[1])*rhy; 269 d2 = (user->bottom[0]-user->bottom[1])*rhx; 270 ft += PetscSqrtScalar(1.0 + d1*d1 + d2*d2); 271 272 /* Top right corner */ 273 d1 = (user->right[my+1] - user->right[my])*rhy; 274 d2 = (user->top[mx+1] - user->top[mx])*rhx; 275 ft += PetscSqrtScalar(1.0 + d1*d1 + d2*d2); 276 277 (*fcn)=ft*area; 278 279 /* Restore vectors */ 280 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); 281 ierr = VecRestoreArray(G,&g);CHKERRQ(ierr); 282 ierr = PetscLogFlops(67.0*mx*my);CHKERRQ(ierr); 283 PetscFunctionReturn(0); 284 } 285 286 /* ------------------------------------------------------------------- */ 287 /* 288 FormHessian - Evaluates the Hessian matrix. 289 290 Input Parameters: 291 . tao - the Tao context 292 . x - input vector 293 . ptr - optional user-defined context, as set by TaoSetHessian() 294 295 Output Parameters: 296 . H - Hessian matrix 297 . Hpre - optionally different preconditioning matrix 298 . flg - flag indicating matrix structure 299 300 */ 301 PetscErrorCode FormHessian(Tao tao,Vec X,Mat H, Mat Hpre, void *ptr) 302 { 303 PetscErrorCode ierr; 304 AppCtx *user = (AppCtx *) ptr; 305 306 PetscFunctionBeginUser; 307 /* Evaluate the Hessian entries*/ 308 ierr = QuadraticH(user,X,H);CHKERRQ(ierr); 309 PetscFunctionReturn(0); 310 } 311 312 /* ------------------------------------------------------------------- */ 313 /* 314 QuadraticH - Evaluates the Hessian matrix. 315 316 Input Parameters: 317 . user - user-defined context, as set by TaoSetHessian() 318 . X - input vector 319 320 Output Parameter: 321 . H - Hessian matrix 322 */ 323 PetscErrorCode QuadraticH(AppCtx *user, Vec X, Mat Hessian) 324 { 325 PetscErrorCode ierr; 326 PetscInt i,j,k,row; 327 PetscInt mx=user->mx, my=user->my; 328 PetscInt col[7]; 329 PetscReal hx=1.0/(mx+1), hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy; 330 PetscReal rhx=mx+1, rhy=my+1; 331 PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb; 332 PetscReal hl,hr,ht,hb,hc,htl,hbr; 333 const PetscScalar *x; 334 PetscReal v[7]; 335 336 PetscFunctionBeginUser; 337 /* Get pointers to vector data */ 338 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); 339 340 /* Initialize matrix entries to zero */ 341 ierr = MatZeroEntries(Hessian);CHKERRQ(ierr); 342 343 /* Set various matrix options */ 344 ierr = MatSetOption(Hessian,MAT_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);CHKERRQ(ierr); 345 346 /* Compute Hessian over the locally owned part of the mesh */ 347 for (i=0; i< mx; i++) { 348 for (j=0; j<my; j++) { 349 350 row=(j)*mx + (i); 351 352 xc = x[row]; 353 xlt=xrb=xl=xr=xb=xt=xc; 354 355 /* Left side */ 356 if (i==0) { 357 xl= user->left[j+1]; 358 xlt = user->left[j+2]; 359 } else { 360 xl = x[row-1]; 361 } 362 363 if (j==0) { 364 xb=user->bottom[i+1]; 365 xrb = user->bottom[i+2]; 366 } else { 367 xb = x[row-mx]; 368 } 369 370 if (i+1 == mx) { 371 xr=user->right[j+1]; 372 xrb = user->right[j]; 373 } else { 374 xr = x[row+1]; 375 } 376 377 if (j+1==my) { 378 xt=user->top[i+1]; 379 xlt = user->top[i]; 380 }else { 381 xt = x[row+mx]; 382 } 383 384 if (i>0 && j+1<my) { 385 xlt = x[row-1+mx]; 386 } 387 if (j>0 && i+1<mx) { 388 xrb = x[row+1-mx]; 389 } 390 391 d1 = (xc-xl)*rhx; 392 d2 = (xc-xr)*rhx; 393 d3 = (xc-xt)*rhy; 394 d4 = (xc-xb)*rhy; 395 d5 = (xrb-xr)*rhy; 396 d6 = (xrb-xb)*rhx; 397 d7 = (xlt-xl)*rhy; 398 d8 = (xlt-xt)*rhx; 399 400 f1 = PetscSqrtScalar(1.0 + d1*d1 + d7*d7); 401 f2 = PetscSqrtScalar(1.0 + d1*d1 + d4*d4); 402 f3 = PetscSqrtScalar(1.0 + d3*d3 + d8*d8); 403 f4 = PetscSqrtScalar(1.0 + d3*d3 + d2*d2); 404 f5 = PetscSqrtScalar(1.0 + d2*d2 + d5*d5); 405 f6 = PetscSqrtScalar(1.0 + d4*d4 + d6*d6); 406 407 hl = (-hydhx*(1.0+d7*d7)+d1*d7)/(f1*f1*f1)+(-hydhx*(1.0+d4*d4)+d1*d4)/(f2*f2*f2); 408 hr = (-hydhx*(1.0+d5*d5)+d2*d5)/(f5*f5*f5)+(-hydhx*(1.0+d3*d3)+d2*d3)/(f4*f4*f4); 409 ht = (-hxdhy*(1.0+d8*d8)+d3*d8)/(f3*f3*f3)+(-hxdhy*(1.0+d2*d2)+d2*d3)/(f4*f4*f4); 410 hb = (-hxdhy*(1.0+d6*d6)+d4*d6)/(f6*f6*f6)+(-hxdhy*(1.0+d1*d1)+d1*d4)/(f2*f2*f2); 411 412 hbr = -d2*d5/(f5*f5*f5) - d4*d6/(f6*f6*f6); 413 htl = -d1*d7/(f1*f1*f1) - d3*d8/(f3*f3*f3); 414 415 hc = hydhx*(1.0+d7*d7)/(f1*f1*f1) + hxdhy*(1.0+d8*d8)/(f3*f3*f3) + hydhx*(1.0+d5*d5)/(f5*f5*f5) + hxdhy*(1.0+d6*d6)/(f6*f6*f6) + 416 (hxdhy*(1.0+d1*d1)+hydhx*(1.0+d4*d4)-2*d1*d4)/(f2*f2*f2) + (hxdhy*(1.0+d2*d2)+hydhx*(1.0+d3*d3)-2*d2*d3)/(f4*f4*f4); 417 418 hl*=0.5; hr*=0.5; ht*=0.5; hb*=0.5; hbr*=0.5; htl*=0.5; hc*=0.5; 419 420 k=0; 421 if (j>0) { 422 v[k]=hb; col[k]=row - mx; k++; 423 } 424 425 if (j>0 && i < mx -1) { 426 v[k]=hbr; col[k]=row - mx+1; k++; 427 } 428 429 if (i>0) { 430 v[k]= hl; col[k]=row - 1; k++; 431 } 432 433 v[k]= hc; col[k]=row; k++; 434 435 if (i < mx-1) { 436 v[k]= hr; col[k]=row+1; k++; 437 } 438 439 if (i>0 && j < my-1) { 440 v[k]= htl; col[k] = row+mx-1; k++; 441 } 442 443 if (j < my-1) { 444 v[k]= ht; col[k] = row+mx; k++; 445 } 446 447 /* 448 Set matrix values using local numbering, which was defined 449 earlier, in the main routine. 450 */ 451 ierr = MatSetValues(Hessian,1,&row,k,col,v,INSERT_VALUES);CHKERRQ(ierr); 452 } 453 } 454 455 /* Restore vectors */ 456 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); 457 458 /* Assemble the matrix */ 459 ierr = MatAssemblyBegin(Hessian,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 460 ierr = MatAssemblyEnd(Hessian,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 461 462 ierr = PetscLogFlops(199.0*mx*my);CHKERRQ(ierr); 463 PetscFunctionReturn(0); 464 } 465 466 /* ------------------------------------------------------------------- */ 467 /* 468 MSA_BoundaryConditions - Calculates the boundary conditions for 469 the region. 470 471 Input Parameter: 472 . user - user-defined application context 473 474 Output Parameter: 475 . user - user-defined application context 476 */ 477 static PetscErrorCode MSA_BoundaryConditions(AppCtx * user) 478 { 479 PetscErrorCode ierr; 480 PetscInt i,j,k,limit=0; 481 PetscInt maxits=5; 482 PetscInt mx=user->mx,my=user->my; 483 PetscInt bsize=0, lsize=0, tsize=0, rsize=0; 484 PetscReal one=1.0, two=2.0, three=3.0, tol=1e-10; 485 PetscReal fnorm,det,hx,hy,xt=0,yt=0; 486 PetscReal u1,u2,nf1,nf2,njac11,njac12,njac21,njac22; 487 PetscReal b=-0.5, t=0.5, l=-0.5, r=0.5; 488 PetscReal *boundary; 489 490 PetscFunctionBeginUser; 491 bsize=mx+2; lsize=my+2; rsize=my+2; tsize=mx+2; 492 493 ierr = PetscMalloc1(bsize,&user->bottom);CHKERRQ(ierr); 494 ierr = PetscMalloc1(tsize,&user->top);CHKERRQ(ierr); 495 ierr = PetscMalloc1(lsize,&user->left);CHKERRQ(ierr); 496 ierr = PetscMalloc1(rsize,&user->right);CHKERRQ(ierr); 497 498 hx= (r-l)/(mx+1); hy=(t-b)/(my+1); 499 500 for (j=0; j<4; j++) { 501 if (j==0) { 502 yt=b; 503 xt=l; 504 limit=bsize; 505 boundary=user->bottom; 506 } else if (j==1) { 507 yt=t; 508 xt=l; 509 limit=tsize; 510 boundary=user->top; 511 } else if (j==2) { 512 yt=b; 513 xt=l; 514 limit=lsize; 515 boundary=user->left; 516 } else { /* if (j==3) */ 517 yt=b; 518 xt=r; 519 limit=rsize; 520 boundary=user->right; 521 } 522 523 for (i=0; i<limit; i++) { 524 u1=xt; 525 u2=-yt; 526 for (k=0; k<maxits; k++) { 527 nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt; 528 nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt; 529 fnorm=PetscSqrtScalar(nf1*nf1+nf2*nf2); 530 if (fnorm <= tol) break; 531 njac11=one+u2*u2-u1*u1; 532 njac12=two*u1*u2; 533 njac21=-two*u1*u2; 534 njac22=-one - u1*u1 + u2*u2; 535 det = njac11*njac22-njac21*njac12; 536 u1 = u1-(njac22*nf1-njac12*nf2)/det; 537 u2 = u2-(njac11*nf2-njac21*nf1)/det; 538 } 539 540 boundary[i]=u1*u1-u2*u2; 541 if (j==0 || j==1) { 542 xt=xt+hx; 543 } else { /* if (j==2 || j==3) */ 544 yt=yt+hy; 545 } 546 } 547 } 548 PetscFunctionReturn(0); 549 } 550 551 /* ------------------------------------------------------------------- */ 552 /* 553 MSA_InitialPoint - Calculates the initial guess in one of three ways. 554 555 Input Parameters: 556 . user - user-defined application context 557 . X - vector for initial guess 558 559 Output Parameters: 560 . X - newly computed initial guess 561 */ 562 static PetscErrorCode MSA_InitialPoint(AppCtx * user, Vec X) 563 { 564 PetscInt start=-1,i,j; 565 PetscErrorCode ierr; 566 PetscReal zero=0.0; 567 PetscBool flg; 568 569 PetscFunctionBeginUser; 570 ierr = VecSet(X, zero);CHKERRQ(ierr); 571 ierr = PetscOptionsGetInt(NULL,NULL,"-start",&start,&flg);CHKERRQ(ierr); 572 573 if (flg && start==0) { /* The zero vector is reasonable */ 574 ierr = VecSet(X, zero);CHKERRQ(ierr); 575 } else { /* Take an average of the boundary conditions */ 576 PetscInt row; 577 PetscInt mx=user->mx,my=user->my; 578 PetscReal *x; 579 580 /* Get pointers to vector data */ 581 ierr = VecGetArray(X,&x);CHKERRQ(ierr); 582 /* Perform local computations */ 583 for (j=0; j<my; j++) { 584 for (i=0; i< mx; i++) { 585 row=(j)*mx + (i); 586 x[row] = (((j+1)*user->bottom[i+1]+(my-j+1)*user->top[i+1])/(my+2)+ ((i+1)*user->left[j+1]+(mx-i+1)*user->right[j+1])/(mx+2))/2.0; 587 } 588 } 589 /* Restore vectors */ 590 ierr = VecRestoreArray(X,&x);CHKERRQ(ierr); 591 } 592 PetscFunctionReturn(0); 593 } 594 595 /*TEST 596 597 build: 598 requires: !complex 599 600 test: 601 args: -tao_smonitor -tao_type nls -mx 10 -my 8 -tao_gatol 1.e-4 602 requires: !single 603 604 test: 605 suffix: 2 606 args: -tao_smonitor -tao_type bmrm -mx 10 -my 8 -tao_gatol 1.e-3 607 requires: !single 608 609 test: 610 suffix: 3 611 args: -tao_smonitor -tao_type lmvm -mx 10 -my 8 -tao_gatol 1.e-3 612 requires: !single 613 614 test: 615 suffix: 4 616 args: -tao_smonitor -tao_type bntr -mx 10 -my 8 -tao_gatol 1.e-4 617 requires: !single 618 619 test: 620 suffix: 5 621 args: -tao_smonitor -tao_type bntl -mx 10 -my 8 -tao_gatol 1.e-4 622 requires: !single 623 624 test: 625 suffix: 6 626 args: -tao_smonitor -tao_type bnls -mx 10 -my 8 -tao_gatol 1.e-4 627 requires: !single 628 629 test: 630 suffix: 7 631 args: -tao_smonitor -tao_type bntr -mx 10 -my 8 -tao_bnk_max_cg_its 3 -tao_gatol 1.e-4 632 requires: !single 633 634 test: 635 suffix: 8 636 args: -tao_smonitor -tao_type bntl -mx 10 -my 8 -tao_bnk_max_cg_its 3 -tao_gatol 1.e-4 637 requires: !single 638 639 test: 640 suffix: 9 641 args: -tao_smonitor -tao_type bnls -mx 10 -my 8 -tao_bnk_max_cg_its 3 -tao_gatol 1.e-4 642 requires: !single 643 644 test: 645 suffix: 10 646 args: -tao_smonitor -tao_type bnls -mx 10 -my 8 -tao_bnk_max_cg_its 3 -tao_gatol 1.e-4 -tao_mf_hessian 647 648 test: 649 suffix: 11 650 args: -tao_smonitor -tao_type bntr -mx 10 -my 8 -tao_gatol 1.e-4 -tao_mf_hessian 651 requires: !single 652 653 test: 654 suffix: 12 655 args: -tao_smonitor -tao_type bntl -mx 10 -my 8 -tao_gatol 1.e-4 -tao_mf_hessian 656 requires: !single 657 658 TEST*/ 659