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