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