1 #include <petsc-private/petscimpl.h> 2 #include <petsctao.h> /*I "petsctao.h" I*/ 3 4 5 PETSC_STATIC_INLINE PetscReal Fischer(PetscReal a, PetscReal b) 6 { 7 /* Method suggested by Bob Vanderbei */ 8 if (a + b <= 0) { 9 return PetscSqrtScalar(a*a + b*b) - (a + b); 10 } 11 return -2.0*a*b / (PetscSqrtScalar(a*a + b*b) + (a + b)); 12 } 13 14 #undef __FUNCT__ 15 #define __FUNCT__ "VecFischer" 16 /*@ 17 VecFischer - Evaluates the Fischer-Burmeister function for complementarity 18 problems. 19 20 Logically Collective on vectors 21 22 Input Parameters: 23 + X - current point 24 . F - function evaluated at x 25 . L - lower bounds 26 - U - upper bounds 27 28 Output Parameters: 29 . FB - The Fischer-Burmeister function vector 30 31 Notes: 32 The Fischer-Burmeister function is defined as 33 $ phi(a,b) := sqrt(a*a + b*b) - a - b 34 and is used reformulate a complementarity problem as a semismooth 35 system of equations. 36 37 The result of this function is done by cases: 38 + l[i] == -infinity, u[i] == infinity -- fb[i] = -f[i] 39 . l[i] == -infinity, u[i] finite -- fb[i] = phi(u[i]-x[i], -f[i]) 40 . l[i] finite, u[i] == infinity -- fb[i] = phi(x[i]-l[i], f[i]) 41 . l[i] finite < u[i] finite -- fb[i] = phi(x[i]-l[i], phi(u[i]-x[i], -f[u])) 42 - otherwise l[i] == u[i] -- fb[i] = l[i] - x[i] 43 44 Level: developer 45 46 @*/ 47 PetscErrorCode VecFischer(Vec X, Vec F, Vec L, Vec U, Vec FB) 48 { 49 PetscReal *x, *f, *l, *u, *fb; 50 PetscReal xval, fval, lval, uval; 51 PetscErrorCode ierr; 52 PetscInt low[5], high[5], n, i; 53 54 PetscFunctionBegin; 55 PetscValidHeaderSpecific(X, VEC_CLASSID,1); 56 PetscValidHeaderSpecific(F, VEC_CLASSID,2); 57 PetscValidHeaderSpecific(L, VEC_CLASSID,3); 58 PetscValidHeaderSpecific(U, VEC_CLASSID,4); 59 PetscValidHeaderSpecific(FB, VEC_CLASSID,4); 60 61 ierr = VecGetOwnershipRange(X, low, high);CHKERRQ(ierr); 62 ierr = VecGetOwnershipRange(F, low + 1, high + 1);CHKERRQ(ierr); 63 ierr = VecGetOwnershipRange(L, low + 2, high + 2);CHKERRQ(ierr); 64 ierr = VecGetOwnershipRange(U, low + 3, high + 3);CHKERRQ(ierr); 65 ierr = VecGetOwnershipRange(FB, low + 4, high + 4);CHKERRQ(ierr); 66 67 for (i = 1; i < 4; ++i) { 68 if (low[0] != low[i] || high[0] != high[i]) SETERRQ(PETSC_COMM_SELF,1,"Vectors must be identically loaded over processors"); 69 } 70 71 ierr = VecGetArray(X, &x);CHKERRQ(ierr); 72 ierr = VecGetArray(F, &f);CHKERRQ(ierr); 73 ierr = VecGetArray(L, &l);CHKERRQ(ierr); 74 ierr = VecGetArray(U, &u);CHKERRQ(ierr); 75 ierr = VecGetArray(FB, &fb);CHKERRQ(ierr); 76 77 ierr = VecGetLocalSize(X, &n);CHKERRQ(ierr); 78 79 for (i = 0; i < n; ++i) { 80 xval = x[i]; fval = f[i]; 81 lval = l[i]; uval = u[i]; 82 83 if ((lval <= -PETSC_INFINITY) && (uval >= PETSC_INFINITY)) { 84 fb[i] = -fval; 85 } else if (lval <= -PETSC_INFINITY) { 86 fb[i] = -Fischer(uval - xval, -fval); 87 } else if (uval >= PETSC_INFINITY) { 88 fb[i] = Fischer(xval - lval, fval); 89 } else if (lval == uval) { 90 fb[i] = lval - xval; 91 } else { 92 fval = Fischer(uval - xval, -fval); 93 fb[i] = Fischer(xval - lval, fval); 94 } 95 } 96 97 ierr = VecRestoreArray(X, &x);CHKERRQ(ierr); 98 ierr = VecRestoreArray(F, &f);CHKERRQ(ierr); 99 ierr = VecRestoreArray(L, &l);CHKERRQ(ierr); 100 ierr = VecRestoreArray(U, &u);CHKERRQ(ierr); 101 ierr = VecRestoreArray(FB, &fb);CHKERRQ(ierr); 102 PetscFunctionReturn(0); 103 } 104 105 PETSC_STATIC_INLINE PetscReal SFischer(PetscReal a, PetscReal b, PetscReal c) 106 { 107 /* Method suggested by Bob Vanderbei */ 108 if (a + b <= 0) { 109 return PetscSqrtScalar(a*a + b*b + 2.0*c*c) - (a + b); 110 } 111 return 2.0*(c*c - a*b) / (PetscSqrtScalar(a*a + b*b + 2.0*c*c) + (a + b)); 112 } 113 114 #undef __FUNCT__ 115 #define __FUNCT__ "VecSFischer" 116 /*@ 117 VecSFischer - Evaluates the Smoothed Fischer-Burmeister function for 118 complementarity problems. 119 120 Logically Collective on vectors 121 122 Input Parameters: 123 + X - current point 124 . F - function evaluated at x 125 . L - lower bounds 126 . U - upper bounds 127 - mu - smoothing parameter 128 129 Output Parameters: 130 . FB - The Smoothed Fischer-Burmeister function vector 131 132 Notes: 133 The Smoothed Fischer-Burmeister function is defined as 134 $ phi(a,b) := sqrt(a*a + b*b + 2*mu*mu) - a - b 135 and is used reformulate a complementarity problem as a semismooth 136 system of equations. 137 138 The result of this function is done by cases: 139 + l[i] == -infinity, u[i] == infinity -- fb[i] = -f[i] - 2*mu*x[i] 140 . l[i] == -infinity, u[i] finite -- fb[i] = phi(u[i]-x[i], -f[i], mu) 141 . l[i] finite, u[i] == infinity -- fb[i] = phi(x[i]-l[i], f[i], mu) 142 . l[i] finite < u[i] finite -- fb[i] = phi(x[i]-l[i], phi(u[i]-x[i], -f[u], mu), mu) 143 - otherwise l[i] == u[i] -- fb[i] = l[i] - x[i] 144 145 Level: developer 146 147 .seealso VecFischer() 148 @*/ 149 PetscErrorCode VecSFischer(Vec X, Vec F, Vec L, Vec U, PetscReal mu, Vec FB) 150 { 151 PetscReal *x, *f, *l, *u, *fb; 152 PetscReal xval, fval, lval, uval; 153 PetscErrorCode ierr; 154 PetscInt low[5], high[5], n, i; 155 156 PetscFunctionBegin; 157 PetscValidHeaderSpecific(X, VEC_CLASSID,1); 158 PetscValidHeaderSpecific(F, VEC_CLASSID,2); 159 PetscValidHeaderSpecific(L, VEC_CLASSID,3); 160 PetscValidHeaderSpecific(U, VEC_CLASSID,4); 161 PetscValidHeaderSpecific(FB, VEC_CLASSID,6); 162 163 ierr = VecGetOwnershipRange(X, low, high);CHKERRQ(ierr); 164 ierr = VecGetOwnershipRange(F, low + 1, high + 1);CHKERRQ(ierr); 165 ierr = VecGetOwnershipRange(L, low + 2, high + 2);CHKERRQ(ierr); 166 ierr = VecGetOwnershipRange(U, low + 3, high + 3);CHKERRQ(ierr); 167 ierr = VecGetOwnershipRange(FB, low + 4, high + 4);CHKERRQ(ierr); 168 169 for (i = 1; i < 4; ++i) { 170 if (low[0] != low[i] || high[0] != high[i]) SETERRQ(PETSC_COMM_SELF,1,"Vectors must be identically loaded over processors"); 171 } 172 173 ierr = VecGetArray(X, &x);CHKERRQ(ierr); 174 ierr = VecGetArray(F, &f);CHKERRQ(ierr); 175 ierr = VecGetArray(L, &l);CHKERRQ(ierr); 176 ierr = VecGetArray(U, &u);CHKERRQ(ierr); 177 ierr = VecGetArray(FB, &fb);CHKERRQ(ierr); 178 179 ierr = VecGetLocalSize(X, &n);CHKERRQ(ierr); 180 181 for (i = 0; i < n; ++i) { 182 xval = (*x++); fval = (*f++); 183 lval = (*l++); uval = (*u++); 184 185 if ((lval <= -PETSC_INFINITY) && (uval >= PETSC_INFINITY)) { 186 (*fb++) = -fval - mu*xval; 187 } else if (lval <= -PETSC_INFINITY) { 188 (*fb++) = -SFischer(uval - xval, -fval, mu); 189 } else if (uval >= PETSC_INFINITY) { 190 (*fb++) = SFischer(xval - lval, fval, mu); 191 } else if (lval == uval) { 192 (*fb++) = lval - xval; 193 } else { 194 fval = SFischer(uval - xval, -fval, mu); 195 (*fb++) = SFischer(xval - lval, fval, mu); 196 } 197 } 198 x -= n; f -= n; l -=n; u -= n; fb -= n; 199 200 ierr = VecRestoreArray(X, &x);CHKERRQ(ierr); 201 ierr = VecRestoreArray(F, &f);CHKERRQ(ierr); 202 ierr = VecRestoreArray(L, &l);CHKERRQ(ierr); 203 ierr = VecRestoreArray(U, &u);CHKERRQ(ierr); 204 ierr = VecRestoreArray(FB, &fb);CHKERRQ(ierr); 205 PetscFunctionReturn(0); 206 } 207 208 PETSC_STATIC_INLINE PetscReal fischnorm(PetscReal a, PetscReal b) 209 { 210 return PetscSqrtScalar(a*a + b*b); 211 } 212 213 PETSC_STATIC_INLINE PetscReal fischsnorm(PetscReal a, PetscReal b, PetscReal c) 214 { 215 return PetscSqrtScalar(a*a + b*b + 2.0*c*c); 216 } 217 218 #undef __FUNCT__ 219 #define __FUNCT__ "D_Fischer" 220 /*@ 221 D_Fischer - Calculates an element of the B-subdifferential of the 222 Fischer-Burmeister function for complementarity problems. 223 224 Collective on jac 225 226 Input Parameters: 227 + jac - the jacobian of f at X 228 . X - current point 229 . Con - constraints function evaluated at X 230 . XL - lower bounds 231 . XU - upper bounds 232 . t1 - work vector 233 - t2 - work vector 234 235 Output Parameters: 236 + Da - diagonal perturbation component of the result 237 - Db - row scaling component of the result 238 239 Level: developer 240 241 .seealso: VecFischer() 242 @*/ 243 PetscErrorCode D_Fischer(Mat jac, Vec X, Vec Con, Vec XL, Vec XU, Vec T1, Vec T2, Vec Da, Vec Db) 244 { 245 PetscErrorCode ierr; 246 PetscInt i,nn; 247 PetscReal *x,*f,*l,*u,*da,*db,*t1,*t2; 248 PetscReal ai,bi,ci,di,ei; 249 250 PetscFunctionBegin; 251 ierr = VecGetLocalSize(X,&nn);CHKERRQ(ierr); 252 ierr = VecGetArray(X,&x);CHKERRQ(ierr); 253 ierr = VecGetArray(Con,&f);CHKERRQ(ierr); 254 ierr = VecGetArray(XL,&l);CHKERRQ(ierr); 255 ierr = VecGetArray(XU,&u);CHKERRQ(ierr); 256 ierr = VecGetArray(Da,&da);CHKERRQ(ierr); 257 ierr = VecGetArray(Db,&db);CHKERRQ(ierr); 258 ierr = VecGetArray(T1,&t1);CHKERRQ(ierr); 259 ierr = VecGetArray(T2,&t2);CHKERRQ(ierr); 260 261 for (i = 0; i < nn; i++) { 262 da[i] = 0; 263 db[i] = 0; 264 t1[i] = 0; 265 266 if (PetscAbsReal(f[i]) <= PETSC_MACHINE_EPSILON) { 267 if (l[i] > PETSC_NINFINITY && PetscAbsReal(x[i] - l[i]) <= PETSC_MACHINE_EPSILON) { 268 t1[i] = 1; 269 da[i] = 1; 270 } 271 272 if (u[i] < PETSC_INFINITY && PetscAbsReal(u[i] - x[i]) <= PETSC_MACHINE_EPSILON) { 273 t1[i] = 1; 274 db[i] = 1; 275 } 276 } 277 } 278 279 ierr = VecRestoreArray(T1,&t1);CHKERRQ(ierr); 280 ierr = VecRestoreArray(T2,&t2);CHKERRQ(ierr); 281 ierr = MatMult(jac,T1,T2);CHKERRQ(ierr); 282 ierr = VecGetArray(T2,&t2);CHKERRQ(ierr); 283 284 for (i = 0; i < nn; i++) { 285 if ((l[i] <= PETSC_NINFINITY) && (u[i] >= PETSC_INFINITY)) { 286 da[i] = 0; 287 db[i] = -1; 288 } else if (l[i] <= PETSC_NINFINITY) { 289 if (db[i] >= 1) { 290 ai = fischnorm(1, t2[i]); 291 292 da[i] = -1/ai - 1; 293 db[i] = -t2[i]/ai - 1; 294 } else { 295 bi = u[i] - x[i]; 296 ai = fischnorm(bi, f[i]); 297 ai = PetscMax(PETSC_MACHINE_EPSILON, ai); 298 299 da[i] = bi / ai - 1; 300 db[i] = -f[i] / ai - 1; 301 } 302 } else if (u[i] >= PETSC_INFINITY) { 303 if (da[i] >= 1) { 304 ai = fischnorm(1, t2[i]); 305 306 da[i] = 1 / ai - 1; 307 db[i] = t2[i] / ai - 1; 308 } else { 309 bi = x[i] - l[i]; 310 ai = fischnorm(bi, f[i]); 311 ai = PetscMax(PETSC_MACHINE_EPSILON, ai); 312 313 da[i] = bi / ai - 1; 314 db[i] = f[i] / ai - 1; 315 } 316 } else if (l[i] == u[i]) { 317 da[i] = -1; 318 db[i] = 0; 319 } else { 320 if (db[i] >= 1) { 321 ai = fischnorm(1, t2[i]); 322 323 ci = 1 / ai + 1; 324 di = t2[i] / ai + 1; 325 } else { 326 bi = x[i] - u[i]; 327 ai = fischnorm(bi, f[i]); 328 ai = PetscMax(PETSC_MACHINE_EPSILON, ai); 329 330 ci = bi / ai + 1; 331 di = f[i] / ai + 1; 332 } 333 334 if (da[i] >= 1) { 335 bi = ci + di*t2[i]; 336 ai = fischnorm(1, bi); 337 338 bi = bi / ai - 1; 339 ai = 1 / ai - 1; 340 } else { 341 ei = Fischer(u[i] - x[i], -f[i]); 342 ai = fischnorm(x[i] - l[i], ei); 343 ai = PetscMax(PETSC_MACHINE_EPSILON, ai); 344 345 bi = ei / ai - 1; 346 ai = (x[i] - l[i]) / ai - 1; 347 } 348 349 da[i] = ai + bi*ci; 350 db[i] = bi*di; 351 } 352 } 353 354 ierr = VecRestoreArray(Da,&da);CHKERRQ(ierr); 355 ierr = VecRestoreArray(Db,&db);CHKERRQ(ierr); 356 ierr = VecRestoreArray(X,&x);CHKERRQ(ierr); 357 ierr = VecRestoreArray(Con,&f);CHKERRQ(ierr); 358 ierr = VecRestoreArray(XL,&l);CHKERRQ(ierr); 359 ierr = VecRestoreArray(XU,&u);CHKERRQ(ierr); 360 ierr = VecRestoreArray(T2,&t2);CHKERRQ(ierr); 361 PetscFunctionReturn(0); 362 } 363 364 #undef __FUNCT__ 365 #define __FUNCT__ "D_SFischer" 366 /*@ 367 D_SFischer - Calculates an element of the B-subdifferential of the 368 smoothed Fischer-Burmeister function for complementarity problems. 369 370 Collective on jac 371 372 Input Parameters: 373 + jac - the jacobian of f at X 374 . X - current point 375 . F - constraint function evaluated at X 376 . XL - lower bounds 377 . XU - upper bounds 378 . mu - smoothing parameter 379 . T1 - work vector 380 - T2 - work vector 381 382 Output Parameter: 383 + Da - diagonal perturbation component of the result 384 . Db - row scaling component of the result 385 - Dm - derivative with respect to scaling parameter 386 387 Level: developer 388 389 .seealso D_Fischer() 390 @*/ 391 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) 392 { 393 PetscErrorCode ierr; 394 PetscInt i,nn; 395 PetscReal *x, *f, *l, *u, *da, *db, *dm; 396 PetscReal ai, bi, ci, di, ei, fi; 397 398 PetscFunctionBegin; 399 if (PetscAbsReal(mu) <= PETSC_MACHINE_EPSILON) { 400 ierr = VecZeroEntries(Dm);CHKERRQ(ierr); 401 ierr = D_Fischer(jac, X, Con, XL, XU, T1, T2, Da, Db);CHKERRQ(ierr); 402 } else { 403 ierr = VecGetLocalSize(X,&nn);CHKERRQ(ierr); 404 ierr = VecGetArray(X,&x);CHKERRQ(ierr); 405 ierr = VecGetArray(Con,&f);CHKERRQ(ierr); 406 ierr = VecGetArray(XL,&l);CHKERRQ(ierr); 407 ierr = VecGetArray(XU,&u);CHKERRQ(ierr); 408 ierr = VecGetArray(Da,&da);CHKERRQ(ierr); 409 ierr = VecGetArray(Db,&db);CHKERRQ(ierr); 410 ierr = VecGetArray(Dm,&dm);CHKERRQ(ierr); 411 412 for (i = 0; i < nn; ++i) { 413 if ((l[i] <= PETSC_NINFINITY) && (u[i] >= PETSC_INFINITY)) { 414 da[i] = -mu; 415 db[i] = -1; 416 dm[i] = -x[i]; 417 } else if (l[i] <= PETSC_NINFINITY) { 418 bi = u[i] - x[i]; 419 ai = fischsnorm(bi, f[i], mu); 420 ai = PetscMax(PETSC_MACHINE_EPSILON, ai); 421 422 da[i] = bi / ai - 1; 423 db[i] = -f[i] / ai - 1; 424 dm[i] = 2.0 * mu / ai; 425 } else if (u[i] >= PETSC_INFINITY) { 426 bi = x[i] - l[i]; 427 ai = fischsnorm(bi, f[i], mu); 428 ai = PetscMax(PETSC_MACHINE_EPSILON, ai); 429 430 da[i] = bi / ai - 1; 431 db[i] = f[i] / ai - 1; 432 dm[i] = 2.0 * mu / ai; 433 } else if (l[i] == u[i]) { 434 da[i] = -1; 435 db[i] = 0; 436 dm[i] = 0; 437 } else { 438 bi = x[i] - u[i]; 439 ai = fischsnorm(bi, f[i], mu); 440 ai = PetscMax(PETSC_MACHINE_EPSILON, ai); 441 442 ci = bi / ai + 1; 443 di = f[i] / ai + 1; 444 fi = 2.0 * mu / ai; 445 446 ei = SFischer(u[i] - x[i], -f[i], mu); 447 ai = fischsnorm(x[i] - l[i], ei, mu); 448 ai = PetscMax(PETSC_MACHINE_EPSILON, ai); 449 450 bi = ei / ai - 1; 451 ei = 2.0 * mu / ei; 452 ai = (x[i] - l[i]) / ai - 1; 453 454 da[i] = ai + bi*ci; 455 db[i] = bi*di; 456 dm[i] = ei + bi*fi; 457 } 458 } 459 460 ierr = VecRestoreArray(X,&x);CHKERRQ(ierr); 461 ierr = VecRestoreArray(Con,&f);CHKERRQ(ierr); 462 ierr = VecRestoreArray(XL,&l);CHKERRQ(ierr); 463 ierr = VecRestoreArray(XU,&u);CHKERRQ(ierr); 464 ierr = VecRestoreArray(Da,&da);CHKERRQ(ierr); 465 ierr = VecRestoreArray(Db,&db);CHKERRQ(ierr); 466 ierr = VecRestoreArray(Dm,&dm);CHKERRQ(ierr); 467 } 468 PetscFunctionReturn(0); 469 } 470 471