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