1 #include <petsc/private/petscimpl.h> 2 #include <petsctao.h> /*I "petsctao.h" I*/ 3 #include <petscsys.h> 4 5 static inline PetscReal Fischer(PetscReal a, PetscReal b) 6 { 7 /* Method suggested by Bob Vanderbei */ 8 if (a + b <= 0) return PetscSqrtReal(a * a + b * b) - (a + b); 9 return -2.0 * a * b / (PetscSqrtReal(a * a + b * b) + (a + b)); 10 } 11 12 /*@ 13 VecFischer - Evaluates the Fischer-Burmeister function for complementarity 14 problems. 15 16 Logically Collective 17 18 Input Parameters: 19 + X - current point 20 . F - function evaluated at x 21 . L - lower bounds 22 - U - upper bounds 23 24 Output Parameter: 25 . FB - The Fischer-Burmeister function vector 26 27 Level: developer 28 29 Notes: 30 The Fischer-Burmeister function is defined as 31 32 $$ 33 \phi(a,b) := \sqrt{a^2 + b^2} - a - b 34 $$ 35 36 and is used reformulate a complementarity problem as a semismooth 37 system of equations. 38 39 The result of this function is done by cases\: 40 + l[i] == -infinity, u[i] == infinity -- fb[i] = -f[i] 41 . l[i] == -infinity, u[i] finite -- fb[i] = phi(u[i]-x[i], -f[i]) 42 . l[i] finite, u[i] == infinity -- fb[i] = phi(x[i]-l[i], f[i]) 43 . l[i] finite < u[i] finite -- fb[i] = phi(x[i]-l[i], phi(u[i]-x[i], -f[u])) 44 - otherwise l[i] == u[i] -- fb[i] = l[i] - x[i] 45 46 .seealso: `Vec`, `VecSFischer()`, `MatDFischer()`, `MatDSFischer()` 47 @*/ 48 PetscErrorCode VecFischer(Vec X, Vec F, Vec L, Vec U, Vec FB) 49 { 50 const PetscScalar *x, *f, *l, *u; 51 PetscScalar *fb; 52 PetscReal xval, fval, lval, uval; 53 PetscInt low[5], high[5], n, i; 54 55 PetscFunctionBegin; 56 PetscValidHeaderSpecific(X, VEC_CLASSID, 1); 57 PetscValidHeaderSpecific(F, VEC_CLASSID, 2); 58 if (L) PetscValidHeaderSpecific(L, VEC_CLASSID, 3); 59 if (U) PetscValidHeaderSpecific(U, VEC_CLASSID, 4); 60 PetscValidHeaderSpecific(FB, VEC_CLASSID, 5); 61 62 if (!L && !U) { 63 PetscCall(VecAXPBY(FB, -1.0, 0.0, F)); 64 PetscFunctionReturn(PETSC_SUCCESS); 65 } 66 67 PetscCall(VecGetOwnershipRange(X, low, high)); 68 PetscCall(VecGetOwnershipRange(F, low + 1, high + 1)); 69 PetscCall(VecGetOwnershipRange(L, low + 2, high + 2)); 70 PetscCall(VecGetOwnershipRange(U, low + 3, high + 3)); 71 PetscCall(VecGetOwnershipRange(FB, low + 4, high + 4)); 72 73 for (i = 1; i < 4; ++i) PetscCheck(low[0] == low[i] && high[0] == high[i], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Vectors must be identically loaded over processors"); 74 75 PetscCall(VecGetArrayRead(X, &x)); 76 PetscCall(VecGetArrayRead(F, &f)); 77 PetscCall(VecGetArrayRead(L, &l)); 78 PetscCall(VecGetArrayRead(U, &u)); 79 PetscCall(VecGetArray(FB, &fb)); 80 81 PetscCall(VecGetLocalSize(X, &n)); 82 83 for (i = 0; i < n; ++i) { 84 xval = PetscRealPart(x[i]); 85 fval = PetscRealPart(f[i]); 86 lval = PetscRealPart(l[i]); 87 uval = PetscRealPart(u[i]); 88 89 if (lval <= -PETSC_INFINITY && uval >= PETSC_INFINITY) { 90 fb[i] = -fval; 91 } else if (lval <= -PETSC_INFINITY) { 92 fb[i] = -Fischer(uval - xval, -fval); 93 } else if (uval >= PETSC_INFINITY) { 94 fb[i] = Fischer(xval - lval, fval); 95 } else if (lval == uval) { 96 fb[i] = lval - xval; 97 } else { 98 fval = Fischer(uval - xval, -fval); 99 fb[i] = Fischer(xval - lval, fval); 100 } 101 } 102 103 PetscCall(VecRestoreArrayRead(X, &x)); 104 PetscCall(VecRestoreArrayRead(F, &f)); 105 PetscCall(VecRestoreArrayRead(L, &l)); 106 PetscCall(VecRestoreArrayRead(U, &u)); 107 PetscCall(VecRestoreArray(FB, &fb)); 108 PetscFunctionReturn(PETSC_SUCCESS); 109 } 110 111 static inline PetscReal SFischer(PetscReal a, PetscReal b, PetscReal c) 112 { 113 /* Method suggested by Bob Vanderbei */ 114 if (a + b <= 0) return PetscSqrtReal(a * a + b * b + 2.0 * c * c) - (a + b); 115 return 2.0 * (c * c - a * b) / (PetscSqrtReal(a * a + b * b + 2.0 * c * c) + (a + b)); 116 } 117 118 /*@ 119 VecSFischer - Evaluates the Smoothed Fischer-Burmeister function for 120 complementarity problems. 121 122 Logically Collective 123 124 Input Parameters: 125 + X - current point 126 . F - function evaluated at x 127 . L - lower bounds 128 . U - upper bounds 129 - mu - smoothing parameter 130 131 Output Parameter: 132 . FB - The Smoothed Fischer-Burmeister function vector 133 134 Notes: 135 The Smoothed Fischer-Burmeister function is defined as 136 $ phi(a,b) := sqrt(a*a + b*b + 2*mu*mu) - a - b 137 and is used reformulate a complementarity problem as a semismooth 138 system of equations. 139 140 The result of this function is done by cases\: 141 + l[i] == -infinity, u[i] == infinity -- fb[i] = -f[i] - 2*mu*x[i] 142 . l[i] == -infinity, u[i] finite -- fb[i] = phi(u[i]-x[i], -f[i], mu) 143 . l[i] finite, u[i] == infinity -- fb[i] = phi(x[i]-l[i], f[i], mu) 144 . l[i] finite < u[i] finite -- fb[i] = phi(x[i]-l[i], phi(u[i]-x[i], -f[u], mu), mu) 145 - otherwise l[i] == u[i] -- fb[i] = l[i] - x[i] 146 147 Level: developer 148 149 .seealso: `Vec`, `VecFischer()`, `MatDFischer()`, `MatDSFischer()` 150 @*/ 151 PetscErrorCode VecSFischer(Vec X, Vec F, Vec L, Vec U, PetscReal mu, Vec FB) 152 { 153 const PetscScalar *x, *f, *l, *u; 154 PetscScalar *fb; 155 PetscReal xval, fval, lval, uval; 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 PetscCall(VecGetOwnershipRange(X, low, high)); 166 PetscCall(VecGetOwnershipRange(F, low + 1, high + 1)); 167 PetscCall(VecGetOwnershipRange(L, low + 2, high + 2)); 168 PetscCall(VecGetOwnershipRange(U, low + 3, high + 3)); 169 PetscCall(VecGetOwnershipRange(FB, low + 4, high + 4)); 170 171 for (i = 1; i < 4; ++i) PetscCheck(low[0] == low[i] && high[0] == high[i], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Vectors must be identically loaded over processors"); 172 173 PetscCall(VecGetArrayRead(X, &x)); 174 PetscCall(VecGetArrayRead(F, &f)); 175 PetscCall(VecGetArrayRead(L, &l)); 176 PetscCall(VecGetArrayRead(U, &u)); 177 PetscCall(VecGetArray(FB, &fb)); 178 179 PetscCall(VecGetLocalSize(X, &n)); 180 181 for (i = 0; i < n; ++i) { 182 xval = PetscRealPart(*x++); 183 fval = PetscRealPart(*f++); 184 lval = PetscRealPart(*l++); 185 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; 201 f -= n; 202 l -= n; 203 u -= n; 204 fb -= n; 205 206 PetscCall(VecRestoreArrayRead(X, &x)); 207 PetscCall(VecRestoreArrayRead(F, &f)); 208 PetscCall(VecRestoreArrayRead(L, &l)); 209 PetscCall(VecRestoreArrayRead(U, &u)); 210 PetscCall(VecRestoreArray(FB, &fb)); 211 PetscFunctionReturn(PETSC_SUCCESS); 212 } 213 214 static inline PetscReal fischnorm(PetscReal a, PetscReal b) 215 { 216 return PetscSqrtReal(a * a + b * b); 217 } 218 219 static inline PetscReal fischsnorm(PetscReal a, PetscReal b, PetscReal c) 220 { 221 return PetscSqrtReal(a * a + b * b + 2.0 * c * c); 222 } 223 224 /*@ 225 MatDFischer - Calculates an element of the B-subdifferential of the 226 Fischer-Burmeister function for complementarity problems. 227 228 Collective 229 230 Input Parameters: 231 + jac - the jacobian of `f` at `X` 232 . X - current point 233 . Con - constraints function evaluated at `X` 234 . XL - lower bounds 235 . XU - upper bounds 236 . T1 - work vector 237 - T2 - work vector 238 239 Output Parameters: 240 + Da - diagonal perturbation component of the result 241 - Db - row scaling component of the result 242 243 Level: developer 244 245 .seealso: `Mat`, `VecFischer()`, `VecSFischer()`, `MatDSFischer()` 246 @*/ 247 PetscErrorCode MatDFischer(Mat jac, Vec X, Vec Con, Vec XL, Vec XU, Vec T1, Vec T2, Vec Da, Vec Db) 248 { 249 PetscInt i, nn; 250 const PetscScalar *x, *f, *l, *u, *t2; 251 PetscScalar *da, *db, *t1; 252 PetscReal ai, bi, ci, di, ei; 253 254 PetscFunctionBegin; 255 PetscCall(VecGetLocalSize(X, &nn)); 256 PetscCall(VecGetArrayRead(X, &x)); 257 PetscCall(VecGetArrayRead(Con, &f)); 258 PetscCall(VecGetArrayRead(XL, &l)); 259 PetscCall(VecGetArrayRead(XU, &u)); 260 PetscCall(VecGetArray(Da, &da)); 261 PetscCall(VecGetArray(Db, &db)); 262 PetscCall(VecGetArray(T1, &t1)); 263 PetscCall(VecGetArrayRead(T2, &t2)); 264 265 for (i = 0; i < nn; i++) { 266 da[i] = 0.0; 267 db[i] = 0.0; 268 t1[i] = 0.0; 269 270 if (PetscAbsScalar(f[i]) <= PETSC_MACHINE_EPSILON) { 271 if (PetscRealPart(l[i]) > PETSC_NINFINITY && PetscAbsScalar(x[i] - l[i]) <= PETSC_MACHINE_EPSILON) { 272 t1[i] = 1.0; 273 da[i] = 1.0; 274 } 275 276 if (PetscRealPart(u[i]) < PETSC_INFINITY && PetscAbsScalar(u[i] - x[i]) <= PETSC_MACHINE_EPSILON) { 277 t1[i] = 1.0; 278 db[i] = 1.0; 279 } 280 } 281 } 282 283 PetscCall(VecRestoreArray(T1, &t1)); 284 PetscCall(VecRestoreArrayRead(T2, &t2)); 285 PetscCall(MatMult(jac, T1, T2)); 286 PetscCall(VecGetArrayRead(T2, &t2)); 287 288 for (i = 0; i < nn; i++) { 289 if ((PetscRealPart(l[i]) <= PETSC_NINFINITY) && (PetscRealPart(u[i]) >= PETSC_INFINITY)) { 290 da[i] = 0.0; 291 db[i] = -1.0; 292 } else if (PetscRealPart(l[i]) <= PETSC_NINFINITY) { 293 if (PetscRealPart(db[i]) >= 1) { 294 ai = fischnorm(1.0, PetscRealPart(t2[i])); 295 296 da[i] = -1.0 / ai - 1.0; 297 db[i] = -t2[i] / ai - 1.0; 298 } else { 299 bi = PetscRealPart(u[i]) - PetscRealPart(x[i]); 300 ai = fischnorm(bi, PetscRealPart(f[i])); 301 ai = PetscMax(PETSC_MACHINE_EPSILON, ai); 302 303 da[i] = bi / ai - 1.0; 304 db[i] = -f[i] / ai - 1.0; 305 } 306 } else if (PetscRealPart(u[i]) >= PETSC_INFINITY) { 307 if (PetscRealPart(da[i]) >= 1) { 308 ai = fischnorm(1.0, PetscRealPart(t2[i])); 309 310 da[i] = 1.0 / ai - 1.0; 311 db[i] = t2[i] / ai - 1.0; 312 } else { 313 bi = PetscRealPart(x[i]) - PetscRealPart(l[i]); 314 ai = fischnorm(bi, PetscRealPart(f[i])); 315 ai = PetscMax(PETSC_MACHINE_EPSILON, ai); 316 317 da[i] = bi / ai - 1.0; 318 db[i] = f[i] / ai - 1.0; 319 } 320 } else if (PetscRealPart(l[i]) == PetscRealPart(u[i])) { 321 da[i] = -1.0; 322 db[i] = 0.0; 323 } else { 324 if (PetscRealPart(db[i]) >= 1) { 325 ai = fischnorm(1.0, PetscRealPart(t2[i])); 326 327 ci = 1.0 / ai + 1.0; 328 di = PetscRealPart(t2[i]) / ai + 1.0; 329 } else { 330 bi = PetscRealPart(x[i]) - PetscRealPart(u[i]); 331 ai = fischnorm(bi, PetscRealPart(f[i])); 332 ai = PetscMax(PETSC_MACHINE_EPSILON, ai); 333 334 ci = bi / ai + 1.0; 335 di = PetscRealPart(f[i]) / ai + 1.0; 336 } 337 338 if (PetscRealPart(da[i]) >= 1) { 339 bi = ci + di * PetscRealPart(t2[i]); 340 ai = fischnorm(1.0, bi); 341 342 bi = bi / ai - 1.0; 343 ai = 1.0 / ai - 1.0; 344 } else { 345 ei = Fischer(PetscRealPart(u[i]) - PetscRealPart(x[i]), -PetscRealPart(f[i])); 346 ai = fischnorm(PetscRealPart(x[i]) - PetscRealPart(l[i]), ei); 347 ai = PetscMax(PETSC_MACHINE_EPSILON, ai); 348 349 bi = ei / ai - 1.0; 350 ai = (PetscRealPart(x[i]) - PetscRealPart(l[i])) / ai - 1.0; 351 } 352 353 da[i] = ai + bi * ci; 354 db[i] = bi * di; 355 } 356 } 357 358 PetscCall(VecRestoreArray(Da, &da)); 359 PetscCall(VecRestoreArray(Db, &db)); 360 PetscCall(VecRestoreArrayRead(X, &x)); 361 PetscCall(VecRestoreArrayRead(Con, &f)); 362 PetscCall(VecRestoreArrayRead(XL, &l)); 363 PetscCall(VecRestoreArrayRead(XU, &u)); 364 PetscCall(VecRestoreArrayRead(T2, &t2)); 365 PetscFunctionReturn(PETSC_SUCCESS); 366 } 367 368 /*@ 369 MatDSFischer - Calculates an element of the B-subdifferential of the 370 smoothed Fischer-Burmeister function for complementarity problems. 371 372 Collective 373 374 Input Parameters: 375 + jac - the jacobian of f at X 376 . X - current point 377 . Con - constraint function evaluated at X 378 . XL - lower bounds 379 . XU - upper bounds 380 . mu - smoothing parameter 381 . T1 - work vector 382 - T2 - work vector 383 384 Output Parameters: 385 + Da - diagonal perturbation component of the result 386 . Db - row scaling component of the result 387 - Dm - derivative with respect to scaling parameter 388 389 Level: developer 390 391 .seealso: `Mat`, `VecFischer()`, `VecDFischer()`, `MatDFischer()` 392 @*/ 393 PetscErrorCode MatDSFischer(Mat jac, Vec X, Vec Con, Vec XL, Vec XU, PetscReal mu, Vec T1, Vec T2, Vec Da, Vec Db, Vec Dm) 394 { 395 PetscInt i, nn; 396 const PetscScalar *x, *f, *l, *u; 397 PetscScalar *da, *db, *dm; 398 PetscReal ai, bi, ci, di, ei, fi; 399 400 PetscFunctionBegin; 401 if (PetscAbsReal(mu) <= PETSC_MACHINE_EPSILON) { 402 PetscCall(VecZeroEntries(Dm)); 403 PetscCall(MatDFischer(jac, X, Con, XL, XU, T1, T2, Da, Db)); 404 } else { 405 PetscCall(VecGetLocalSize(X, &nn)); 406 PetscCall(VecGetArrayRead(X, &x)); 407 PetscCall(VecGetArrayRead(Con, &f)); 408 PetscCall(VecGetArrayRead(XL, &l)); 409 PetscCall(VecGetArrayRead(XU, &u)); 410 PetscCall(VecGetArray(Da, &da)); 411 PetscCall(VecGetArray(Db, &db)); 412 PetscCall(VecGetArray(Dm, &dm)); 413 414 for (i = 0; i < nn; ++i) { 415 if ((PetscRealPart(l[i]) <= PETSC_NINFINITY) && (PetscRealPart(u[i]) >= PETSC_INFINITY)) { 416 da[i] = -mu; 417 db[i] = -1.0; 418 dm[i] = -x[i]; 419 } else if (PetscRealPart(l[i]) <= PETSC_NINFINITY) { 420 bi = PetscRealPart(u[i]) - PetscRealPart(x[i]); 421 ai = fischsnorm(bi, PetscRealPart(f[i]), mu); 422 ai = PetscMax(PETSC_MACHINE_EPSILON, ai); 423 424 da[i] = bi / ai - 1.0; 425 db[i] = -PetscRealPart(f[i]) / ai - 1.0; 426 dm[i] = 2.0 * mu / ai; 427 } else if (PetscRealPart(u[i]) >= PETSC_INFINITY) { 428 bi = PetscRealPart(x[i]) - PetscRealPart(l[i]); 429 ai = fischsnorm(bi, PetscRealPart(f[i]), mu); 430 ai = PetscMax(PETSC_MACHINE_EPSILON, ai); 431 432 da[i] = bi / ai - 1.0; 433 db[i] = PetscRealPart(f[i]) / ai - 1.0; 434 dm[i] = 2.0 * mu / ai; 435 } else if (PetscRealPart(l[i]) == PetscRealPart(u[i])) { 436 da[i] = -1.0; 437 db[i] = 0.0; 438 dm[i] = 0.0; 439 } else { 440 bi = PetscRealPart(x[i]) - PetscRealPart(u[i]); 441 ai = fischsnorm(bi, PetscRealPart(f[i]), mu); 442 ai = PetscMax(PETSC_MACHINE_EPSILON, ai); 443 444 ci = bi / ai + 1.0; 445 di = PetscRealPart(f[i]) / ai + 1.0; 446 fi = 2.0 * mu / ai; 447 448 ei = SFischer(PetscRealPart(u[i]) - PetscRealPart(x[i]), -PetscRealPart(f[i]), mu); 449 ai = fischsnorm(PetscRealPart(x[i]) - PetscRealPart(l[i]), ei, mu); 450 ai = PetscMax(PETSC_MACHINE_EPSILON, ai); 451 452 bi = ei / ai - 1.0; 453 ei = 2.0 * mu / ei; 454 ai = (PetscRealPart(x[i]) - PetscRealPart(l[i])) / ai - 1.0; 455 456 da[i] = ai + bi * ci; 457 db[i] = bi * di; 458 dm[i] = ei + bi * fi; 459 } 460 } 461 462 PetscCall(VecRestoreArrayRead(X, &x)); 463 PetscCall(VecRestoreArrayRead(Con, &f)); 464 PetscCall(VecRestoreArrayRead(XL, &l)); 465 PetscCall(VecRestoreArrayRead(XU, &u)); 466 PetscCall(VecRestoreArray(Da, &da)); 467 PetscCall(VecRestoreArray(Db, &db)); 468 PetscCall(VecRestoreArray(Dm, &dm)); 469 } 470 PetscFunctionReturn(PETSC_SUCCESS); 471 } 472 473 static inline PetscReal ST_InternalPN(PetscScalar in, PetscReal lb, PetscReal ub) 474 { 475 return PetscMax(0, PetscRealPart(in) - ub) - PetscMax(0, -PetscRealPart(in) - PetscAbsReal(lb)); 476 } 477 478 static inline PetscReal ST_InternalNN(PetscScalar in, PetscReal lb, PetscReal ub) 479 { 480 return PetscMax(0, PetscRealPart(in) + PetscAbsReal(ub)) - PetscMax(0, -PetscRealPart(in) - PetscAbsReal(lb)); 481 } 482 483 static inline PetscReal ST_InternalPP(PetscScalar in, PetscReal lb, PetscReal ub) 484 { 485 return PetscMax(0, PetscRealPart(in) - ub) + PetscMin(0, PetscRealPart(in) - lb); 486 } 487 488 /*@ 489 TaoSoftThreshold - Calculates soft thresholding routine with input vector 490 and given lower and upper bound and returns it to output vector. 491 492 Input Parameters: 493 + in - input vector to be thresholded 494 . lb - lower bound 495 - ub - upper bound 496 497 Output Parameter: 498 . out - Soft thresholded output vector 499 500 Notes: 501 Soft thresholding is defined as 502 \[ S(input,lb,ub) = 503 \begin{cases} 504 input - ub \text{input > ub} \\ 505 0 \text{lb =< input <= ub} \\ 506 input + lb \text{input < lb} \\ 507 \] 508 509 Level: developer 510 511 .seealso: `Tao`, `Vec` 512 @*/ 513 PetscErrorCode TaoSoftThreshold(Vec in, PetscReal lb, PetscReal ub, Vec out) 514 { 515 PetscInt i, nlocal, mlocal; 516 PetscScalar *inarray, *outarray; 517 518 PetscFunctionBegin; 519 PetscCall(VecGetArrayPair(in, out, &inarray, &outarray)); 520 PetscCall(VecGetLocalSize(in, &nlocal)); 521 PetscCall(VecGetLocalSize(out, &mlocal)); 522 523 PetscCheck(nlocal == mlocal, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Input and output vectors need to be of same size."); 524 if (lb == ub) { 525 PetscCall(VecRestoreArrayPair(in, out, &inarray, &outarray)); 526 PetscFunctionReturn(PETSC_SUCCESS); 527 } 528 PetscCheck(lb < ub, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Lower bound needs to be lower than upper bound."); 529 530 if (ub >= 0 && lb < 0) { 531 for (i = 0; i < nlocal; i++) outarray[i] = ST_InternalPN(inarray[i], lb, ub); 532 } else if (ub < 0 && lb < 0) { 533 for (i = 0; i < nlocal; i++) outarray[i] = ST_InternalNN(inarray[i], lb, ub); 534 } else { 535 for (i = 0; i < nlocal; i++) outarray[i] = ST_InternalPP(inarray[i], lb, ub); 536 } 537 538 PetscCall(VecRestoreArrayPair(in, out, &inarray, &outarray)); 539 PetscFunctionReturn(PETSC_SUCCESS); 540 } 541