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