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 /* Get pointers to vector data */ 337 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); 338 339 /* Initialize matrix entries to zero */ 340 ierr = MatZeroEntries(Hessian);CHKERRQ(ierr); 341 342 /* Set various matrix options */ 343 ierr = MatSetOption(Hessian,MAT_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);CHKERRQ(ierr); 344 345 /* Compute Hessian over the locally owned part of the mesh */ 346 for (i=0; i< mx; i++){ 347 for (j=0; j<my; j++){ 348 349 row=(j)*mx + (i); 350 351 xc = x[row]; 352 xlt=xrb=xl=xr=xb=xt=xc; 353 354 /* Left side */ 355 if (i==0){ 356 xl= user->left[j+1]; 357 xlt = user->left[j+2]; 358 } else { 359 xl = x[row-1]; 360 } 361 362 if (j==0){ 363 xb=user->bottom[i+1]; 364 xrb = user->bottom[i+2]; 365 } else { 366 xb = x[row-mx]; 367 } 368 369 if (i+1 == mx){ 370 xr=user->right[j+1]; 371 xrb = user->right[j]; 372 } else { 373 xr = x[row+1]; 374 } 375 376 if (j+1==my){ 377 xt=user->top[i+1]; 378 xlt = user->top[i]; 379 }else { 380 xt = x[row+mx]; 381 } 382 383 if (i>0 && j+1<my){ 384 xlt = x[row-1+mx]; 385 } 386 if (j>0 && i+1<mx){ 387 xrb = x[row+1-mx]; 388 } 389 390 d1 = (xc-xl)*rhx; 391 d2 = (xc-xr)*rhx; 392 d3 = (xc-xt)*rhy; 393 d4 = (xc-xb)*rhy; 394 d5 = (xrb-xr)*rhy; 395 d6 = (xrb-xb)*rhx; 396 d7 = (xlt-xl)*rhy; 397 d8 = (xlt-xt)*rhx; 398 399 f1 = PetscSqrtScalar(1.0 + d1*d1 + d7*d7); 400 f2 = PetscSqrtScalar(1.0 + d1*d1 + d4*d4); 401 f3 = PetscSqrtScalar(1.0 + d3*d3 + d8*d8); 402 f4 = PetscSqrtScalar(1.0 + d3*d3 + d2*d2); 403 f5 = PetscSqrtScalar(1.0 + d2*d2 + d5*d5); 404 f6 = PetscSqrtScalar(1.0 + d4*d4 + d6*d6); 405 406 hl = (-hydhx*(1.0+d7*d7)+d1*d7)/(f1*f1*f1)+(-hydhx*(1.0+d4*d4)+d1*d4)/(f2*f2*f2); 407 hr = (-hydhx*(1.0+d5*d5)+d2*d5)/(f5*f5*f5)+(-hydhx*(1.0+d3*d3)+d2*d3)/(f4*f4*f4); 408 ht = (-hxdhy*(1.0+d8*d8)+d3*d8)/(f3*f3*f3)+(-hxdhy*(1.0+d2*d2)+d2*d3)/(f4*f4*f4); 409 hb = (-hxdhy*(1.0+d6*d6)+d4*d6)/(f6*f6*f6)+(-hxdhy*(1.0+d1*d1)+d1*d4)/(f2*f2*f2); 410 411 hbr = -d2*d5/(f5*f5*f5) - d4*d6/(f6*f6*f6); 412 htl = -d1*d7/(f1*f1*f1) - d3*d8/(f3*f3*f3); 413 414 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) + 415 (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); 416 417 hl*=0.5; hr*=0.5; ht*=0.5; hb*=0.5; hbr*=0.5; htl*=0.5; hc*=0.5; 418 419 k=0; 420 if (j>0){ 421 v[k]=hb; col[k]=row - mx; k++; 422 } 423 424 if (j>0 && i < mx -1){ 425 v[k]=hbr; col[k]=row - mx+1; k++; 426 } 427 428 if (i>0){ 429 v[k]= hl; col[k]=row - 1; k++; 430 } 431 432 v[k]= hc; col[k]=row; k++; 433 434 if (i < mx-1){ 435 v[k]= hr; col[k]=row+1; k++; 436 } 437 438 if (i>0 && j < my-1){ 439 v[k]= htl; col[k] = row+mx-1; k++; 440 } 441 442 if (j < my-1){ 443 v[k]= ht; col[k] = row+mx; k++; 444 } 445 446 /* 447 Set matrix values using local numbering, which was defined 448 earlier, in the main routine. 449 */ 450 ierr = MatSetValues(Hessian,1,&row,k,col,v,INSERT_VALUES);CHKERRQ(ierr); 451 } 452 } 453 454 /* Restore vectors */ 455 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); 456 457 /* Assemble the matrix */ 458 ierr = MatAssemblyBegin(Hessian,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 459 ierr = MatAssemblyEnd(Hessian,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 460 461 ierr = PetscLogFlops(199.0*mx*my);CHKERRQ(ierr); 462 PetscFunctionReturn(0); 463 } 464 465 /* ------------------------------------------------------------------- */ 466 /* 467 MSA_BoundaryConditions - Calculates the boundary conditions for 468 the region. 469 470 Input Parameter: 471 . user - user-defined application context 472 473 Output Parameter: 474 . user - user-defined application context 475 */ 476 static PetscErrorCode MSA_BoundaryConditions(AppCtx * user) 477 { 478 PetscErrorCode ierr; 479 PetscInt i,j,k,limit=0; 480 PetscInt maxits=5; 481 PetscInt mx=user->mx,my=user->my; 482 PetscInt bsize=0, lsize=0, tsize=0, rsize=0; 483 PetscReal one=1.0, two=2.0, three=3.0, tol=1e-10; 484 PetscReal fnorm,det,hx,hy,xt=0,yt=0; 485 PetscReal u1,u2,nf1,nf2,njac11,njac12,njac21,njac22; 486 PetscReal b=-0.5, t=0.5, l=-0.5, r=0.5; 487 PetscReal *boundary; 488 489 PetscFunctionBeginUser; 490 bsize=mx+2; lsize=my+2; rsize=my+2; tsize=mx+2; 491 492 ierr = PetscMalloc1(bsize,&user->bottom);CHKERRQ(ierr); 493 ierr = PetscMalloc1(tsize,&user->top);CHKERRQ(ierr); 494 ierr = PetscMalloc1(lsize,&user->left);CHKERRQ(ierr); 495 ierr = PetscMalloc1(rsize,&user->right);CHKERRQ(ierr); 496 497 hx= (r-l)/(mx+1); hy=(t-b)/(my+1); 498 499 for (j=0; j<4; j++){ 500 if (j==0){ 501 yt=b; 502 xt=l; 503 limit=bsize; 504 boundary=user->bottom; 505 } else if (j==1){ 506 yt=t; 507 xt=l; 508 limit=tsize; 509 boundary=user->top; 510 } else if (j==2){ 511 yt=b; 512 xt=l; 513 limit=lsize; 514 boundary=user->left; 515 } else { /* if (j==3) */ 516 yt=b; 517 xt=r; 518 limit=rsize; 519 boundary=user->right; 520 } 521 522 for (i=0; i<limit; i++){ 523 u1=xt; 524 u2=-yt; 525 for (k=0; k<maxits; k++){ 526 nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt; 527 nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt; 528 fnorm=PetscSqrtScalar(nf1*nf1+nf2*nf2); 529 if (fnorm <= tol) break; 530 njac11=one+u2*u2-u1*u1; 531 njac12=two*u1*u2; 532 njac21=-two*u1*u2; 533 njac22=-one - u1*u1 + u2*u2; 534 det = njac11*njac22-njac21*njac12; 535 u1 = u1-(njac22*nf1-njac12*nf2)/det; 536 u2 = u2-(njac11*nf2-njac21*nf1)/det; 537 } 538 539 boundary[i]=u1*u1-u2*u2; 540 if (j==0 || j==1) { 541 xt=xt+hx; 542 } else { /* if (j==2 || j==3) */ 543 yt=yt+hy; 544 } 545 } 546 } 547 PetscFunctionReturn(0); 548 } 549 550 /* ------------------------------------------------------------------- */ 551 /* 552 MSA_InitialPoint - Calculates the initial guess in one of three ways. 553 554 Input Parameters: 555 . user - user-defined application context 556 . X - vector for initial guess 557 558 Output Parameters: 559 . X - newly computed initial guess 560 */ 561 static PetscErrorCode MSA_InitialPoint(AppCtx * user, Vec X) 562 { 563 PetscInt start=-1,i,j; 564 PetscErrorCode ierr; 565 PetscReal zero=0.0; 566 PetscBool flg; 567 568 PetscFunctionBeginUser; 569 ierr = VecSet(X, zero);CHKERRQ(ierr); 570 ierr = PetscOptionsGetInt(NULL,NULL,"-start",&start,&flg);CHKERRQ(ierr); 571 572 if (flg && start==0){ /* The zero vector is reasonable */ 573 ierr = VecSet(X, zero);CHKERRQ(ierr); 574 } else { /* Take an average of the boundary conditions */ 575 PetscInt row; 576 PetscInt mx=user->mx,my=user->my; 577 PetscReal *x; 578 579 /* Get pointers to vector data */ 580 ierr = VecGetArray(X,&x);CHKERRQ(ierr); 581 /* Perform local computations */ 582 for (j=0; j<my; j++){ 583 for (i=0; i< mx; i++){ 584 row=(j)*mx + (i); 585 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; 586 } 587 } 588 /* Restore vectors */ 589 ierr = VecRestoreArray(X,&x);CHKERRQ(ierr); 590 } 591 PetscFunctionReturn(0); 592 } 593 594 /*TEST 595 596 build: 597 requires: !complex 598 599 test: 600 args: -tao_smonitor -tao_type nls -mx 10 -my 8 -tao_gatol 1.e-4 601 requires: !single 602 603 test: 604 suffix: 2 605 args: -tao_smonitor -tao_type bmrm -mx 10 -my 8 -tao_gatol 1.e-3 606 requires: !single 607 608 test: 609 suffix: 3 610 args: -tao_smonitor -tao_type lmvm -mx 10 -my 8 -tao_gatol 1.e-3 611 requires: !single 612 613 test: 614 suffix: 4 615 args: -tao_smonitor -tao_type bntr -mx 10 -my 8 -tao_gatol 1.e-4 616 requires: !single 617 618 test: 619 suffix: 5 620 args: -tao_smonitor -tao_type bntl -mx 10 -my 8 -tao_gatol 1.e-4 621 requires: !single 622 623 test: 624 suffix: 6 625 args: -tao_smonitor -tao_type bnls -mx 10 -my 8 -tao_gatol 1.e-4 626 requires: !single 627 628 test: 629 suffix: 7 630 args: -tao_smonitor -tao_type bntr -mx 10 -my 8 -tao_bnk_max_cg_its 3 -tao_gatol 1.e-4 631 requires: !single 632 633 test: 634 suffix: 8 635 args: -tao_smonitor -tao_type bntl -mx 10 -my 8 -tao_bnk_max_cg_its 3 -tao_gatol 1.e-4 636 requires: !single 637 638 test: 639 suffix: 9 640 args: -tao_smonitor -tao_type bnls -mx 10 -my 8 -tao_bnk_max_cg_its 3 -tao_gatol 1.e-4 641 requires: !single 642 643 TEST*/ 644