1 #include "petsc-private/petscimpl.h" 2 #include "tao.h" /*I "tao.h" I*/ 3 #include "tao_util.h" /*I "tao_util.h" I*/ 4 5 #undef __FUNCT__ 6 #define __FUNCT__ "VecPow" 7 /*@ 8 VecPow - Replaces each component of a vector by x_i^p 9 10 Logically Collective on v 11 12 Input Parameter: 13 + v - the vector 14 - p - the exponent to use on each element 15 16 Output Parameter: 17 . v - the vector 18 19 Level: intermediate 20 21 @*/ 22 PetscErrorCode VecPow(Vec v, PetscReal p) 23 { 24 PetscErrorCode ierr; 25 PetscInt n,i; 26 PetscReal *v1; 27 28 PetscFunctionBegin; 29 PetscValidHeaderSpecific(v, VEC_CLASSID, 1); 30 31 ierr = VecGetArray(v, &v1); CHKERRQ(ierr); 32 ierr = VecGetLocalSize(v, &n); CHKERRQ(ierr); 33 34 if (1.0 == p) { 35 } 36 else if (-1.0 == p) { 37 for (i = 0; i < n; ++i){ 38 v1[i] = 1.0 / v1[i]; 39 } 40 } 41 else if (0.0 == p) { 42 for (i = 0; i < n; ++i){ 43 /* Not-a-number left alone 44 Infinity set to one */ 45 if (v1[i] == v1[i]) { 46 v1[i] = 1.0; 47 } 48 } 49 } 50 else if (0.5 == p) { 51 for (i = 0; i < n; ++i) { 52 if (v1[i] >= 0) { 53 v1[i] = PetscSqrtScalar(v1[i]); 54 } 55 else { 56 v1[i] = TAO_INFINITY; 57 } 58 } 59 } 60 else if (-0.5 == p) { 61 for (i = 0; i < n; ++i) { 62 if (v1[i] >= 0) { 63 v1[i] = 1.0 / PetscSqrtScalar(v1[i]); 64 } 65 else { 66 v1[i] = TAO_INFINITY; 67 } 68 } 69 } 70 else if (2.0 == p) { 71 for (i = 0; i < n; ++i){ 72 v1[i] *= v1[i]; 73 } 74 } 75 else if (-2.0 == p) { 76 for (i = 0; i < n; ++i){ 77 v1[i] = 1.0 / (v1[i] * v1[i]); 78 } 79 } 80 else { 81 for (i = 0; i < n; ++i) { 82 if (v1[i] >= 0) { 83 v1[i] = PetscPowScalar(v1[i], p); 84 } 85 else { 86 v1[i] = TAO_INFINITY; 87 } 88 } 89 } 90 91 ierr = VecRestoreArray(v,&v1); CHKERRQ(ierr); 92 PetscFunctionReturn(0); 93 } 94 95 #undef __FUNCT__ 96 #define __FUNCT__ "VecMedian" 97 /*@ 98 VecMedian - Computes the componentwise median of three vectors 99 and stores the result in this vector. Used primarily for projecting 100 a vector within upper and lower bounds. 101 102 Logically Collective 103 104 Input Parameters: 105 . Vec1, Vec2, Vec3 - The three vectors 106 107 Output Parameter: 108 . VMedian - The median vector 109 110 Level: advanced 111 @*/ 112 PetscErrorCode VecMedian(Vec Vec1, Vec Vec2, Vec Vec3, Vec VMedian) 113 { 114 PetscErrorCode ierr; 115 PetscInt i,n,low1,low2,low3,low4,high1,high2,high3,high4; 116 PetscReal *v1,*v2,*v3,*vmed; 117 118 PetscFunctionBegin; 119 PetscValidHeaderSpecific(Vec1,VEC_CLASSID,1); 120 PetscValidHeaderSpecific(Vec2,VEC_CLASSID,2); 121 PetscValidHeaderSpecific(Vec3,VEC_CLASSID,3); 122 PetscValidHeaderSpecific(VMedian,VEC_CLASSID,4); 123 124 if (Vec1==Vec2 || Vec1==Vec3){ 125 ierr=VecCopy(Vec1,VMedian); CHKERRQ(ierr); 126 PetscFunctionReturn(0); 127 } 128 if (Vec2==Vec3){ 129 ierr=VecCopy(Vec2,VMedian); CHKERRQ(ierr); 130 PetscFunctionReturn(0); 131 } 132 133 PetscValidType(Vec1,1); 134 PetscValidType(Vec2,2); 135 PetscValidType(VMedian,4); 136 PetscCheckSameType(Vec1,1,Vec2,2); PetscCheckSameType(Vec1,1,VMedian,4); 137 PetscCheckSameComm(Vec1,1,Vec2,2); PetscCheckSameComm(Vec1,1,VMedian,4); 138 139 ierr = VecGetOwnershipRange(Vec1, &low1, &high1); CHKERRQ(ierr); 140 ierr = VecGetOwnershipRange(Vec2, &low2, &high2); CHKERRQ(ierr); 141 ierr = VecGetOwnershipRange(Vec3, &low3, &high3); CHKERRQ(ierr); 142 ierr = VecGetOwnershipRange(VMedian, &low4, &high4); CHKERRQ(ierr); 143 if ( low1!= low2 || low1!= low3 || low1!= low4 || 144 high1!= high2 || high1!= high3 || high1!= high4) 145 SETERRQ(PETSC_COMM_SELF,1,"InCompatible vector local lengths"); 146 147 ierr = VecGetArray(Vec1,&v1); CHKERRQ(ierr); 148 ierr = VecGetArray(Vec2,&v2); CHKERRQ(ierr); 149 ierr = VecGetArray(Vec3,&v3); CHKERRQ(ierr); 150 151 if ( VMedian != Vec1 && VMedian != Vec2 && VMedian != Vec3){ 152 ierr = VecGetArray(VMedian,&vmed); CHKERRQ(ierr); 153 } else if ( VMedian==Vec1 ){ 154 vmed=v1; 155 } else if ( VMedian==Vec2 ){ 156 vmed=v2; 157 } else { 158 vmed=v3; 159 } 160 161 ierr=VecGetLocalSize(Vec1,&n); CHKERRQ(ierr); 162 163 for (i=0;i<n;i++){ 164 vmed[i]=PetscMax(PetscMax(PetscMin(v1[i],v2[i]),PetscMin(v1[i],v3[i])),PetscMin(v2[i],v3[i])); 165 } 166 167 ierr = VecRestoreArray(Vec1,&v1); CHKERRQ(ierr); 168 ierr = VecRestoreArray(Vec2,&v2); CHKERRQ(ierr); 169 ierr = VecRestoreArray(Vec3,&v2); CHKERRQ(ierr); 170 171 if (VMedian!=Vec1 && VMedian != Vec2 && VMedian != Vec3){ 172 ierr = VecRestoreArray(VMedian,&vmed); CHKERRQ(ierr); 173 } 174 175 PetscFunctionReturn(0); 176 } 177 178 PETSC_STATIC_INLINE PetscReal Fischer(PetscReal a, PetscReal b) 179 { 180 /* Method suggested by Bob Vanderbei */ 181 if (a + b <= 0) { 182 return PetscSqrtScalar(a*a + b*b) - (a + b); 183 } 184 return -2.0*a*b / (PetscSqrtScalar(a*a + b*b) + (a + b)); 185 } 186 187 #undef __FUNCT__ 188 #define __FUNCT__ "VecFischer" 189 /*@ 190 VecFischer - Evaluates the Fischer-Burmeister function for complementarity 191 problems. 192 193 Logically Collective on vectors 194 195 Input Parameters: 196 + X - current point 197 . F - function evaluated at x 198 . L - lower bounds 199 - U - upper bounds 200 201 Output Parameters: 202 . FB - The Fischer-Burmeister function vector 203 204 Notes: 205 The Fischer-Burmeister function is defined as 206 $ phi(a,b) := sqrt(a*a + b*b) - a - b 207 and is used reformulate a complementarity problem as a semismooth 208 system of equations. 209 210 The result of this function is done by cases: 211 + l[i] == -infinity, u[i] == infinity -- fb[i] = -f[i] 212 . l[i] == -infinity, u[i] finite -- fb[i] = phi(u[i]-x[i], -f[i]) 213 . l[i] finite, u[i] == infinity -- fb[i] = phi(x[i]-l[i], f[i]) 214 . l[i] finite < u[i] finite -- fb[i] = phi(x[i]-l[i], phi(u[i]-x[i], -f[u])) 215 - otherwise l[i] == u[i] -- fb[i] = l[i] - x[i] 216 217 Level: developer 218 219 @*/ 220 PetscErrorCode VecFischer(Vec X, Vec F, Vec L, Vec U, Vec FB) 221 { 222 PetscReal *x, *f, *l, *u, *fb; 223 PetscReal xval, fval, lval, uval; 224 PetscErrorCode ierr; 225 PetscInt low[5], high[5], n, i; 226 227 PetscFunctionBegin; 228 PetscValidHeaderSpecific(X, VEC_CLASSID,1); 229 PetscValidHeaderSpecific(F, VEC_CLASSID,2); 230 PetscValidHeaderSpecific(L, VEC_CLASSID,3); 231 PetscValidHeaderSpecific(U, VEC_CLASSID,4); 232 PetscValidHeaderSpecific(FB, VEC_CLASSID,4); 233 234 ierr = VecGetOwnershipRange(X, low, high); CHKERRQ(ierr); 235 ierr = VecGetOwnershipRange(F, low + 1, high + 1); CHKERRQ(ierr); 236 ierr = VecGetOwnershipRange(L, low + 2, high + 2); CHKERRQ(ierr); 237 ierr = VecGetOwnershipRange(U, low + 3, high + 3); CHKERRQ(ierr); 238 ierr = VecGetOwnershipRange(FB, low + 4, high + 4); CHKERRQ(ierr); 239 240 for (i = 1; i < 4; ++i) { 241 if (low[0] != low[i] || high[0] != high[i]) 242 SETERRQ(PETSC_COMM_SELF,1,"Vectors must be identically loaded over processors"); 243 } 244 245 ierr = VecGetArray(X, &x); CHKERRQ(ierr); 246 ierr = VecGetArray(F, &f); CHKERRQ(ierr); 247 ierr = VecGetArray(L, &l); CHKERRQ(ierr); 248 ierr = VecGetArray(U, &u); CHKERRQ(ierr); 249 ierr = VecGetArray(FB, &fb); CHKERRQ(ierr); 250 251 ierr = VecGetLocalSize(X, &n); CHKERRQ(ierr); 252 253 for (i = 0; i < n; ++i) { 254 xval = x[i]; fval = f[i]; 255 lval = l[i]; uval = u[i]; 256 257 if ((lval <= -TAO_INFINITY) && (uval >= TAO_INFINITY)) { 258 fb[i] = -fval; 259 } 260 else if (lval <= -TAO_INFINITY) { 261 fb[i] = -Fischer(uval - xval, -fval); 262 } 263 else if (uval >= TAO_INFINITY) { 264 fb[i] = Fischer(xval - lval, fval); 265 } 266 else if (lval == uval) { 267 fb[i] = lval - xval; 268 } 269 else { 270 fval = Fischer(uval - xval, -fval); 271 fb[i] = Fischer(xval - lval, fval); 272 } 273 } 274 275 ierr = VecRestoreArray(X, &x); CHKERRQ(ierr); 276 ierr = VecRestoreArray(F, &f); CHKERRQ(ierr); 277 ierr = VecRestoreArray(L, &l); CHKERRQ(ierr); 278 ierr = VecRestoreArray(U, &u); CHKERRQ(ierr); 279 ierr = VecRestoreArray(FB, &fb); CHKERRQ(ierr); 280 281 PetscFunctionReturn(0); 282 } 283 284 PETSC_STATIC_INLINE PetscReal SFischer(PetscReal a, PetscReal b, PetscReal c) 285 { 286 /* Method suggested by Bob Vanderbei */ 287 if (a + b <= 0) { 288 return PetscSqrtScalar(a*a + b*b + 2.0*c*c) - (a + b); 289 } 290 return 2.0*(c*c - a*b) / (PetscSqrtScalar(a*a + b*b + 2.0*c*c) + (a + b)); 291 } 292 293 #undef __FUNCT__ 294 #define __FUNCT__ "VecSFischer" 295 /*@ 296 VecSFischer - Evaluates the Smoothed Fischer-Burmeister function for 297 complementarity problems. 298 299 Logically Collective on vectors 300 301 Input Parameters: 302 + X - current point 303 . F - function evaluated at x 304 . L - lower bounds 305 . U - upper bounds 306 - mu - smoothing parameter 307 308 Output Parameters: 309 . FB - The Smoothed Fischer-Burmeister function vector 310 311 Notes: 312 The Smoothed Fischer-Burmeister function is defined as 313 $ phi(a,b) := sqrt(a*a + b*b + 2*mu*mu) - a - b 314 and is used reformulate a complementarity problem as a semismooth 315 system of equations. 316 317 The result of this function is done by cases: 318 + l[i] == -infinity, u[i] == infinity -- fb[i] = -f[i] - 2*mu*x[i] 319 . l[i] == -infinity, u[i] finite -- fb[i] = phi(u[i]-x[i], -f[i], mu) 320 . l[i] finite, u[i] == infinity -- fb[i] = phi(x[i]-l[i], f[i], mu) 321 . l[i] finite < u[i] finite -- fb[i] = phi(x[i]-l[i], phi(u[i]-x[i], -f[u], mu), mu) 322 - otherwise l[i] == u[i] -- fb[i] = l[i] - x[i] 323 324 Level: developer 325 326 .seealso VecFischer() 327 @*/ 328 PetscErrorCode VecSFischer(Vec X, Vec F, Vec L, Vec U, PetscReal mu, Vec FB) 329 { 330 PetscReal *x, *f, *l, *u, *fb; 331 PetscReal xval, fval, lval, uval; 332 PetscErrorCode ierr; 333 PetscInt low[5], high[5], n, i; 334 335 PetscFunctionBegin; 336 PetscValidHeaderSpecific(X, VEC_CLASSID,1); 337 PetscValidHeaderSpecific(F, VEC_CLASSID,2); 338 PetscValidHeaderSpecific(L, VEC_CLASSID,3); 339 PetscValidHeaderSpecific(U, VEC_CLASSID,4); 340 PetscValidHeaderSpecific(FB, VEC_CLASSID,6); 341 342 ierr = VecGetOwnershipRange(X, low, high); CHKERRQ(ierr); 343 ierr = VecGetOwnershipRange(F, low + 1, high + 1); CHKERRQ(ierr); 344 ierr = VecGetOwnershipRange(L, low + 2, high + 2); CHKERRQ(ierr); 345 ierr = VecGetOwnershipRange(U, low + 3, high + 3); CHKERRQ(ierr); 346 ierr = VecGetOwnershipRange(FB, low + 4, high + 4); CHKERRQ(ierr); 347 348 for (i = 1; i < 4; ++i) { 349 if (low[0] != low[i] || high[0] != high[i]) 350 SETERRQ(PETSC_COMM_SELF,1,"Vectors must be identically loaded over processors"); 351 } 352 353 ierr = VecGetArray(X, &x); CHKERRQ(ierr); 354 ierr = VecGetArray(F, &f); CHKERRQ(ierr); 355 ierr = VecGetArray(L, &l); CHKERRQ(ierr); 356 ierr = VecGetArray(U, &u); CHKERRQ(ierr); 357 ierr = VecGetArray(FB, &fb); CHKERRQ(ierr); 358 359 ierr = VecGetLocalSize(X, &n); CHKERRQ(ierr); 360 361 for (i = 0; i < n; ++i) { 362 xval = (*x++); fval = (*f++); 363 lval = (*l++); uval = (*u++); 364 365 if ((lval <= -TAO_INFINITY) && (uval >= TAO_INFINITY)) { 366 (*fb++) = -fval - mu*xval; 367 } 368 else if (lval <= -TAO_INFINITY) { 369 (*fb++) = -SFischer(uval - xval, -fval, mu); 370 } 371 else if (uval >= TAO_INFINITY) { 372 (*fb++) = SFischer(xval - lval, fval, mu); 373 } 374 else if (lval == uval) { 375 (*fb++) = lval - xval; 376 } 377 else { 378 fval = SFischer(uval - xval, -fval, mu); 379 (*fb++) = SFischer(xval - lval, fval, mu); 380 } 381 } 382 x -= n; f -= n; l -=n; u -= n; fb -= n; 383 384 ierr = VecRestoreArray(X, &x); CHKERRQ(ierr); 385 ierr = VecRestoreArray(F, &f); CHKERRQ(ierr); 386 ierr = VecRestoreArray(L, &l); CHKERRQ(ierr); 387 ierr = VecRestoreArray(U, &u); CHKERRQ(ierr); 388 ierr = VecRestoreArray(FB, &fb); CHKERRQ(ierr); 389 PetscFunctionReturn(0); 390 } 391 392 PETSC_STATIC_INLINE PetscReal fischnorm(PetscReal a, PetscReal b) 393 { 394 return PetscSqrtScalar(a*a + b*b); 395 } 396 397 PETSC_STATIC_INLINE PetscReal fischsnorm(PetscReal a, PetscReal b, PetscReal c) 398 { 399 return PetscSqrtScalar(a*a + b*b + 2.0*c*c); 400 } 401 402 #undef __FUNCT__ 403 #define __FUNCT__ "D_Fischer" 404 /*@ 405 D_Fischer - Calculates an element of the B-subdifferential of the 406 Fischer-Burmeister function for complementarity problems. 407 408 Collective on jac 409 410 Input Parameters: 411 + jac - the jacobian of f at X 412 . X - current point 413 . Con - constraints function evaluated at X 414 . XL - lower bounds 415 . XU - upper bounds 416 . t1 - work vector 417 - t2 - work vector 418 419 Output Parameters: 420 + Da - diagonal perturbation component of the result 421 - Db - row scaling component of the result 422 423 Level: developer 424 425 .seealso: VecFischer() 426 @*/ 427 PetscErrorCode D_Fischer(Mat jac, Vec X, Vec Con, Vec XL, Vec XU, 428 Vec T1, Vec T2, Vec Da, Vec Db) 429 { 430 PetscErrorCode ierr; 431 PetscInt i,nn; 432 PetscReal *x,*f,*l,*u,*da,*db,*t1,*t2; 433 PetscReal ai,bi,ci,di,ei; 434 435 PetscFunctionBegin; 436 437 ierr = VecGetLocalSize(X,&nn); CHKERRQ(ierr); 438 439 ierr = VecGetArray(X,&x);CHKERRQ(ierr); 440 ierr = VecGetArray(Con,&f);CHKERRQ(ierr); 441 ierr = VecGetArray(XL,&l);CHKERRQ(ierr); 442 ierr = VecGetArray(XU,&u);CHKERRQ(ierr); 443 ierr = VecGetArray(Da,&da);CHKERRQ(ierr); 444 ierr = VecGetArray(Db,&db);CHKERRQ(ierr); 445 ierr = VecGetArray(T1,&t1);CHKERRQ(ierr); 446 ierr = VecGetArray(T2,&t2);CHKERRQ(ierr); 447 448 for (i = 0; i < nn; i++) { 449 da[i] = 0; 450 db[i] = 0; 451 t1[i] = 0; 452 453 if (PetscAbsReal(f[i]) <= PETSC_MACHINE_EPSILON) { 454 if (l[i] > TAO_NINFINITY && PetscAbsReal(x[i] - l[i]) <= PETSC_MACHINE_EPSILON) { 455 t1[i] = 1; 456 da[i] = 1; 457 } 458 459 if (u[i] < TAO_INFINITY && PetscAbsReal(u[i] - x[i]) <= PETSC_MACHINE_EPSILON) { 460 t1[i] = 1; 461 db[i] = 1; 462 } 463 } 464 } 465 466 ierr = VecRestoreArray(T1,&t1); CHKERRQ(ierr); 467 ierr = VecRestoreArray(T2,&t2); CHKERRQ(ierr); 468 ierr = MatMult(jac,T1,T2); CHKERRQ(ierr); 469 ierr = VecGetArray(T2,&t2); CHKERRQ(ierr); 470 471 for (i = 0; i < nn; i++) { 472 if ((l[i] <= TAO_NINFINITY) && (u[i] >= TAO_INFINITY)) { 473 da[i] = 0; 474 db[i] = -1; 475 } 476 else if (l[i] <= TAO_NINFINITY) { 477 if (db[i] >= 1) { 478 ai = fischnorm(1, t2[i]); 479 480 da[i] = -1/ai - 1; 481 db[i] = -t2[i]/ai - 1; 482 } 483 else { 484 bi = u[i] - x[i]; 485 ai = fischnorm(bi, f[i]); 486 ai = PetscMax(PETSC_MACHINE_EPSILON, ai); 487 488 da[i] = bi / ai - 1; 489 db[i] = -f[i] / ai - 1; 490 } 491 } 492 else if (u[i] >= TAO_INFINITY) { 493 if (da[i] >= 1) { 494 ai = fischnorm(1, t2[i]); 495 496 da[i] = 1 / ai - 1; 497 db[i] = t2[i] / ai - 1; 498 } 499 else { 500 bi = x[i] - l[i]; 501 ai = fischnorm(bi, f[i]); 502 ai = PetscMax(PETSC_MACHINE_EPSILON, ai); 503 504 da[i] = bi / ai - 1; 505 db[i] = f[i] / ai - 1; 506 } 507 } 508 else if (l[i] == u[i]) { 509 da[i] = -1; 510 db[i] = 0; 511 } 512 else { 513 if (db[i] >= 1) { 514 ai = fischnorm(1, t2[i]); 515 516 ci = 1 / ai + 1; 517 di = t2[i] / ai + 1; 518 } 519 else { 520 bi = x[i] - u[i]; 521 ai = fischnorm(bi, f[i]); 522 ai = PetscMax(PETSC_MACHINE_EPSILON, ai); 523 524 ci = bi / ai + 1; 525 di = f[i] / ai + 1; 526 } 527 528 if (da[i] >= 1) { 529 bi = ci + di*t2[i]; 530 ai = fischnorm(1, bi); 531 532 bi = bi / ai - 1; 533 ai = 1 / ai - 1; 534 } 535 else { 536 ei = Fischer(u[i] - x[i], -f[i]); 537 ai = fischnorm(x[i] - l[i], ei); 538 ai = PetscMax(PETSC_MACHINE_EPSILON, ai); 539 540 bi = ei / ai - 1; 541 ai = (x[i] - l[i]) / ai - 1; 542 } 543 544 da[i] = ai + bi*ci; 545 db[i] = bi*di; 546 } 547 } 548 549 ierr = VecRestoreArray(Da,&da); CHKERRQ(ierr); 550 ierr = VecRestoreArray(Db,&db); CHKERRQ(ierr); 551 ierr = VecRestoreArray(X,&x); CHKERRQ(ierr); 552 ierr = VecRestoreArray(Con,&f); CHKERRQ(ierr); 553 ierr = VecRestoreArray(XL,&l); CHKERRQ(ierr); 554 ierr = VecRestoreArray(XU,&u); CHKERRQ(ierr); 555 ierr = VecRestoreArray(T2,&t2); CHKERRQ(ierr); 556 PetscFunctionReturn(0); 557 } 558 559 #undef __FUNCT__ 560 #define __FUNCT__ "D_SFischer" 561 /*@ 562 D_SFischer - Calculates an element of the B-subdifferential of the 563 smoothed Fischer-Burmeister function for complementarity problems. 564 565 Collective on jac 566 567 Input Parameters: 568 + jac - the jacobian of f at X 569 . X - current point 570 . F - constraint function evaluated at X 571 . XL - lower bounds 572 . XU - upper bounds 573 . mu - smoothing parameter 574 . T1 - work vector 575 - T2 - work vector 576 577 Output Parameter: 578 + Da - diagonal perturbation component of the result 579 . Db - row scaling component of the result 580 - Dm - derivative with respect to scaling parameter 581 582 Level: developer 583 584 .seealso D_Fischer() 585 @*/ 586 PetscErrorCode D_SFischer(Mat jac, Vec X, Vec Con, 587 Vec XL, Vec XU, PetscReal mu, 588 Vec T1, Vec T2, 589 Vec Da, Vec Db, Vec Dm) 590 { 591 PetscErrorCode ierr; 592 PetscInt i,nn; 593 PetscReal *x, *f, *l, *u, *da, *db, *dm; 594 PetscReal ai, bi, ci, di, ei, fi; 595 596 PetscFunctionBegin; 597 598 if (PetscAbsReal(mu) <= PETSC_MACHINE_EPSILON) { 599 ierr = VecZeroEntries(Dm); CHKERRQ(ierr); 600 ierr = D_Fischer(jac, X, Con, XL, XU, T1, T2, Da, Db); CHKERRQ(ierr); 601 } 602 else { 603 ierr = VecGetLocalSize(X,&nn); CHKERRQ(ierr); 604 ierr = VecGetArray(X,&x); CHKERRQ(ierr); 605 ierr = VecGetArray(Con,&f); CHKERRQ(ierr); 606 ierr = VecGetArray(XL,&l); CHKERRQ(ierr); 607 ierr = VecGetArray(XU,&u); CHKERRQ(ierr); 608 ierr = VecGetArray(Da,&da); CHKERRQ(ierr); 609 ierr = VecGetArray(Db,&db); CHKERRQ(ierr); 610 ierr = VecGetArray(Dm,&dm); CHKERRQ(ierr); 611 612 for (i = 0; i < nn; ++i) { 613 if ((l[i] <= TAO_NINFINITY) && (u[i] >= TAO_INFINITY)) { 614 da[i] = -mu; 615 db[i] = -1; 616 dm[i] = -x[i]; 617 } 618 else if (l[i] <= TAO_NINFINITY) { 619 bi = u[i] - x[i]; 620 ai = fischsnorm(bi, f[i], mu); 621 ai = PetscMax(PETSC_MACHINE_EPSILON, ai); 622 623 da[i] = bi / ai - 1; 624 db[i] = -f[i] / ai - 1; 625 dm[i] = 2.0 * mu / ai; 626 } 627 else if (u[i] >= TAO_INFINITY) { 628 bi = x[i] - l[i]; 629 ai = fischsnorm(bi, f[i], mu); 630 ai = PetscMax(PETSC_MACHINE_EPSILON, ai); 631 632 da[i] = bi / ai - 1; 633 db[i] = f[i] / ai - 1; 634 dm[i] = 2.0 * mu / ai; 635 } 636 else if (l[i] == u[i]) { 637 da[i] = -1; 638 db[i] = 0; 639 dm[i] = 0; 640 } 641 else { 642 bi = x[i] - u[i]; 643 ai = fischsnorm(bi, f[i], mu); 644 ai = PetscMax(PETSC_MACHINE_EPSILON, ai); 645 646 ci = bi / ai + 1; 647 di = f[i] / ai + 1; 648 fi = 2.0 * mu / ai; 649 650 ei = SFischer(u[i] - x[i], -f[i], mu); 651 ai = fischsnorm(x[i] - l[i], ei, mu); 652 ai = PetscMax(PETSC_MACHINE_EPSILON, ai); 653 654 bi = ei / ai - 1; 655 ei = 2.0 * mu / ei; 656 ai = (x[i] - l[i]) / ai - 1; 657 658 da[i] = ai + bi*ci; 659 db[i] = bi*di; 660 dm[i] = ei + bi*fi; 661 } 662 } 663 664 ierr = VecRestoreArray(X,&x); CHKERRQ(ierr); 665 ierr = VecRestoreArray(Con,&f); CHKERRQ(ierr); 666 ierr = VecRestoreArray(XL,&l); CHKERRQ(ierr); 667 ierr = VecRestoreArray(XU,&u); CHKERRQ(ierr); 668 ierr = VecRestoreArray(Da,&da); CHKERRQ(ierr); 669 ierr = VecRestoreArray(Db,&db); CHKERRQ(ierr); 670 ierr = VecRestoreArray(Dm,&dm); CHKERRQ(ierr); 671 672 } 673 PetscFunctionReturn(0); 674 } 675 676