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