Actual source code: projection.c

  1: #include <petsc/private/vecimpl.h>

  3: /*@
  4:   VecWhichEqual - Creates an index set containing the indices
  5:              where the vectors Vec1 and Vec2 have identical elements.

  7:   Collective on Vec

  9:   Input Parameters:
 10: . Vec1, Vec2 - the two vectors to compare

 12:   OutputParameter:
 13: . S - The index set containing the indices i where vec1[i] == vec2[i]

 15:   Notes:
 16:     the two vectors must have the same parallel layout

 18:   Level: advanced
 19: @*/
 20: PetscErrorCode VecWhichEqual(Vec Vec1, Vec Vec2, IS *S)
 21: {
 22:   PetscInt           i, n_same = 0;
 23:   PetscInt           n, low, high;
 24:   PetscInt          *same = NULL;
 25:   const PetscScalar *v1, *v2;

 30:   VecCheckSameSize(Vec1, 1, Vec2, 2);

 32:   VecGetOwnershipRange(Vec1, &low, &high);
 33:   VecGetLocalSize(Vec1, &n);
 34:   if (n > 0) {
 35:     if (Vec1 == Vec2) {
 36:       VecGetArrayRead(Vec1, &v1);
 37:       v2 = v1;
 38:     } else {
 39:       VecGetArrayRead(Vec1, &v1);
 40:       VecGetArrayRead(Vec2, &v2);
 41:     }

 43:     PetscMalloc1(n, &same);

 45:     for (i = 0; i < n; ++i) {
 46:       if (v1[i] == v2[i]) {
 47:         same[n_same] = low + i;
 48:         ++n_same;
 49:       }
 50:     }

 52:     if (Vec1 == Vec2) {
 53:       VecRestoreArrayRead(Vec1, &v1);
 54:     } else {
 55:       VecRestoreArrayRead(Vec1, &v1);
 56:       VecRestoreArrayRead(Vec2, &v2);
 57:     }
 58:   }
 59:   ISCreateGeneral(PetscObjectComm((PetscObject)Vec1), n_same, same, PETSC_OWN_POINTER, S);
 60:   return 0;
 61: }

 63: /*@
 64:   VecWhichLessThan - Creates an index set containing the indices
 65:   where the vectors Vec1 < Vec2

 67:   Collective on S

 69:   Input Parameters:
 70: . Vec1, Vec2 - the two vectors to compare

 72:   OutputParameter:
 73: . S - The index set containing the indices i where vec1[i] < vec2[i]

 75:   Notes:
 76:   The two vectors must have the same parallel layout

 78:   For complex numbers this only compares the real part

 80:   Level: advanced
 81: @*/
 82: PetscErrorCode VecWhichLessThan(Vec Vec1, Vec Vec2, IS *S)
 83: {
 84:   PetscInt           i, n_lt = 0;
 85:   PetscInt           n, low, high;
 86:   PetscInt          *lt = NULL;
 87:   const PetscScalar *v1, *v2;

 92:   VecCheckSameSize(Vec1, 1, Vec2, 2);

 94:   VecGetOwnershipRange(Vec1, &low, &high);
 95:   VecGetLocalSize(Vec1, &n);
 96:   if (n > 0) {
 97:     if (Vec1 == Vec2) {
 98:       VecGetArrayRead(Vec1, &v1);
 99:       v2 = v1;
100:     } else {
101:       VecGetArrayRead(Vec1, &v1);
102:       VecGetArrayRead(Vec2, &v2);
103:     }

105:     PetscMalloc1(n, &lt);

107:     for (i = 0; i < n; ++i) {
108:       if (PetscRealPart(v1[i]) < PetscRealPart(v2[i])) {
109:         lt[n_lt] = low + i;
110:         ++n_lt;
111:       }
112:     }

114:     if (Vec1 == Vec2) {
115:       VecRestoreArrayRead(Vec1, &v1);
116:     } else {
117:       VecRestoreArrayRead(Vec1, &v1);
118:       VecRestoreArrayRead(Vec2, &v2);
119:     }
120:   }
121:   ISCreateGeneral(PetscObjectComm((PetscObject)Vec1), n_lt, lt, PETSC_OWN_POINTER, S);
122:   return 0;
123: }

125: /*@
126:   VecWhichGreaterThan - Creates an index set containing the indices
127:   where the vectors Vec1 > Vec2

129:   Collective on S

131:   Input Parameters:
132: . Vec1, Vec2 - the two vectors to compare

134:   OutputParameter:
135: . S - The index set containing the indices i where vec1[i] > vec2[i]

137:   Notes:
138:   The two vectors must have the same parallel layout

140:   For complex numbers this only compares the real part

142:   Level: advanced
143: @*/
144: PetscErrorCode VecWhichGreaterThan(Vec Vec1, Vec Vec2, IS *S)
145: {
146:   PetscInt           i, n_gt = 0;
147:   PetscInt           n, low, high;
148:   PetscInt          *gt = NULL;
149:   const PetscScalar *v1, *v2;

154:   VecCheckSameSize(Vec1, 1, Vec2, 2);

156:   VecGetOwnershipRange(Vec1, &low, &high);
157:   VecGetLocalSize(Vec1, &n);
158:   if (n > 0) {
159:     if (Vec1 == Vec2) {
160:       VecGetArrayRead(Vec1, &v1);
161:       v2 = v1;
162:     } else {
163:       VecGetArrayRead(Vec1, &v1);
164:       VecGetArrayRead(Vec2, &v2);
165:     }

167:     PetscMalloc1(n, &gt);

169:     for (i = 0; i < n; ++i) {
170:       if (PetscRealPart(v1[i]) > PetscRealPart(v2[i])) {
171:         gt[n_gt] = low + i;
172:         ++n_gt;
173:       }
174:     }

176:     if (Vec1 == Vec2) {
177:       VecRestoreArrayRead(Vec1, &v1);
178:     } else {
179:       VecRestoreArrayRead(Vec1, &v1);
180:       VecRestoreArrayRead(Vec2, &v2);
181:     }
182:   }
183:   ISCreateGeneral(PetscObjectComm((PetscObject)Vec1), n_gt, gt, PETSC_OWN_POINTER, S);
184:   return 0;
185: }

187: /*@
188:   VecWhichBetween - Creates an index set containing the indices
189:                where  VecLow < V < VecHigh

191:   Collective on S

193:   Input Parameters:
194: + VecLow - lower bound
195: . V - Vector to compare
196: - VecHigh - higher bound

198:   OutputParameter:
199: . S - The index set containing the indices i where veclow[i] < v[i] < vechigh[i]

201:   Notes:
202:   The vectors must have the same parallel layout

204:   For complex numbers this only compares the real part

206:   Level: advanced
207: @*/
208: PetscErrorCode VecWhichBetween(Vec VecLow, Vec V, Vec VecHigh, IS *S)
209: {
210:   PetscInt           i, n_vm = 0;
211:   PetscInt           n, low, high;
212:   PetscInt          *vm = NULL;
213:   const PetscScalar *v1, *v2, *vmiddle;

220:   VecCheckSameSize(V, 2, VecLow, 1);
221:   VecCheckSameSize(V, 2, VecHigh, 3);

223:   VecGetOwnershipRange(VecLow, &low, &high);
224:   VecGetLocalSize(VecLow, &n);
225:   if (n > 0) {
226:     VecGetArrayRead(VecLow, &v1);
227:     if (VecLow != VecHigh) {
228:       VecGetArrayRead(VecHigh, &v2);
229:     } else {
230:       v2 = v1;
231:     }
232:     if (V != VecLow && V != VecHigh) {
233:       VecGetArrayRead(V, &vmiddle);
234:     } else if (V == VecLow) {
235:       vmiddle = v1;
236:     } else {
237:       vmiddle = v2;
238:     }

240:     PetscMalloc1(n, &vm);

242:     for (i = 0; i < n; ++i) {
243:       if (PetscRealPart(v1[i]) < PetscRealPart(vmiddle[i]) && PetscRealPart(vmiddle[i]) < PetscRealPart(v2[i])) {
244:         vm[n_vm] = low + i;
245:         ++n_vm;
246:       }
247:     }

249:     VecRestoreArrayRead(VecLow, &v1);
250:     if (VecLow != VecHigh) VecRestoreArrayRead(VecHigh, &v2);
251:     if (V != VecLow && V != VecHigh) VecRestoreArrayRead(V, &vmiddle);
252:   }
253:   ISCreateGeneral(PetscObjectComm((PetscObject)V), n_vm, vm, PETSC_OWN_POINTER, S);
254:   return 0;
255: }

257: /*@
258:   VecWhichBetweenOrEqual - Creates an index set containing the indices
259:   where  VecLow <= V <= VecHigh

261:   Collective on S

263:   Input Parameters:
264: + VecLow - lower bound
265: . V - Vector to compare
266: - VecHigh - higher bound

268:   OutputParameter:
269: . S - The index set containing the indices i where veclow[i] <= v[i] <= vechigh[i]

271:   Level: advanced
272: @*/

274: PetscErrorCode VecWhichBetweenOrEqual(Vec VecLow, Vec V, Vec VecHigh, IS *S)
275: {
276:   PetscInt           i, n_vm = 0;
277:   PetscInt           n, low, high;
278:   PetscInt          *vm = NULL;
279:   const PetscScalar *v1, *v2, *vmiddle;

286:   VecCheckSameSize(V, 2, VecLow, 1);
287:   VecCheckSameSize(V, 2, VecHigh, 3);

289:   VecGetOwnershipRange(VecLow, &low, &high);
290:   VecGetLocalSize(VecLow, &n);
291:   if (n > 0) {
292:     VecGetArrayRead(VecLow, &v1);
293:     if (VecLow != VecHigh) {
294:       VecGetArrayRead(VecHigh, &v2);
295:     } else {
296:       v2 = v1;
297:     }
298:     if (V != VecLow && V != VecHigh) {
299:       VecGetArrayRead(V, &vmiddle);
300:     } else if (V == VecLow) {
301:       vmiddle = v1;
302:     } else {
303:       vmiddle = v2;
304:     }

306:     PetscMalloc1(n, &vm);

308:     for (i = 0; i < n; ++i) {
309:       if (PetscRealPart(v1[i]) <= PetscRealPart(vmiddle[i]) && PetscRealPart(vmiddle[i]) <= PetscRealPart(v2[i])) {
310:         vm[n_vm] = low + i;
311:         ++n_vm;
312:       }
313:     }

315:     VecRestoreArrayRead(VecLow, &v1);
316:     if (VecLow != VecHigh) VecRestoreArrayRead(VecHigh, &v2);
317:     if (V != VecLow && V != VecHigh) VecRestoreArrayRead(V, &vmiddle);
318:   }
319:   ISCreateGeneral(PetscObjectComm((PetscObject)V), n_vm, vm, PETSC_OWN_POINTER, S);
320:   return 0;
321: }

323: /*@
324:    VecWhichInactive - Creates an index set containing the indices
325:   where one of the following holds:
326:     a) VecLow(i)  < V(i) < VecHigh(i)
327:     b) VecLow(i)  = V(i) and D(i) <= 0 (< 0 when Strong is true)
328:     c) VecHigh(i) = V(i) and D(i) >= 0 (> 0 when Strong is true)

330:   Collective on S

332:   Input Parameters:
333: + VecLow - lower bound
334: . V - Vector to compare
335: . D - Direction to compare
336: . VecHigh - higher bound
337: - Strong - indicator for applying strongly inactive test

339:   OutputParameter:
340: . S - The index set containing the indices i where the bound is inactive

342:   Level: advanced
343: @*/

345: PetscErrorCode VecWhichInactive(Vec VecLow, Vec V, Vec D, Vec VecHigh, PetscBool Strong, IS *S)
346: {
347:   PetscInt           i, n_vm = 0;
348:   PetscInt           n, low, high;
349:   PetscInt          *vm = NULL;
350:   const PetscScalar *v1, *v2, *v, *d;

352:   if (!VecLow && !VecHigh) {
353:     VecGetOwnershipRange(V, &low, &high);
354:     ISCreateStride(PetscObjectComm((PetscObject)V), high - low, low, 1, S);
355:     return 0;
356:   }
364:   VecCheckSameSize(V, 2, VecLow, 1);
365:   VecCheckSameSize(V, 2, D, 3);
366:   VecCheckSameSize(V, 2, VecHigh, 4);

368:   VecGetOwnershipRange(VecLow, &low, &high);
369:   VecGetLocalSize(VecLow, &n);
370:   if (n > 0) {
371:     VecGetArrayRead(VecLow, &v1);
372:     if (VecLow != VecHigh) {
373:       VecGetArrayRead(VecHigh, &v2);
374:     } else {
375:       v2 = v1;
376:     }
377:     if (V != VecLow && V != VecHigh) {
378:       VecGetArrayRead(V, &v);
379:     } else if (V == VecLow) {
380:       v = v1;
381:     } else {
382:       v = v2;
383:     }
384:     if (D != VecLow && D != VecHigh && D != V) {
385:       VecGetArrayRead(D, &d);
386:     } else if (D == VecLow) {
387:       d = v1;
388:     } else if (D == VecHigh) {
389:       d = v2;
390:     } else {
391:       d = v;
392:     }

394:     PetscMalloc1(n, &vm);

396:     if (Strong) {
397:       for (i = 0; i < n; ++i) {
398:         if (PetscRealPart(v1[i]) < PetscRealPart(v[i]) && PetscRealPart(v[i]) < PetscRealPart(v2[i])) {
399:           vm[n_vm] = low + i;
400:           ++n_vm;
401:         } else if (PetscRealPart(v1[i]) == PetscRealPart(v[i]) && PetscRealPart(d[i]) < 0) {
402:           vm[n_vm] = low + i;
403:           ++n_vm;
404:         } else if (PetscRealPart(v2[i]) == PetscRealPart(v[i]) && PetscRealPart(d[i]) > 0) {
405:           vm[n_vm] = low + i;
406:           ++n_vm;
407:         }
408:       }
409:     } else {
410:       for (i = 0; i < n; ++i) {
411:         if (PetscRealPart(v1[i]) < PetscRealPart(v[i]) && PetscRealPart(v[i]) < PetscRealPart(v2[i])) {
412:           vm[n_vm] = low + i;
413:           ++n_vm;
414:         } else if (PetscRealPart(v1[i]) == PetscRealPart(v[i]) && PetscRealPart(d[i]) <= 0) {
415:           vm[n_vm] = low + i;
416:           ++n_vm;
417:         } else if (PetscRealPart(v2[i]) == PetscRealPart(v[i]) && PetscRealPart(d[i]) >= 0) {
418:           vm[n_vm] = low + i;
419:           ++n_vm;
420:         }
421:       }
422:     }

424:     VecRestoreArrayRead(VecLow, &v1);
425:     if (VecLow != VecHigh) VecRestoreArrayRead(VecHigh, &v2);
426:     if (V != VecLow && V != VecHigh) VecRestoreArrayRead(V, &v);
427:     if (D != VecLow && D != VecHigh && D != V) VecRestoreArrayRead(D, &d);
428:   }
429:   ISCreateGeneral(PetscObjectComm((PetscObject)V), n_vm, vm, PETSC_OWN_POINTER, S);
430:   return 0;
431: }

433: /*@
434:   VecISAXPY - Adds a reduced vector to the appropriate elements of a full-space vector.
435:                   vfull[is[i]] += alpha*vreduced[i]

437:   Input Parameters:
438: + vfull    - the full-space vector
439: . is       - the index set for the reduced space
440: . alpha    - the scalar coefficient
441: - vreduced - the reduced-space vector

443:   Output Parameters:
444: . vfull    - the sum of the full-space vector and reduced-space vector

446:   Notes:
447:     The index set identifies entries in the global vector.
448:     Negative indices are skipped; indices outside the ownership range of vfull will raise an error.

450:   Level: advanced

452: .seealso: `VecISCopy()`, `VecISSet()`, `VecAXPY()`
453: @*/
454: PetscErrorCode VecISAXPY(Vec vfull, IS is, PetscScalar alpha, Vec vreduced)
455: {
456:   PetscInt nfull, nreduced;

461:   VecGetSize(vfull, &nfull);
462:   VecGetSize(vreduced, &nreduced);

464:   if (nfull == nreduced) { /* Also takes care of masked vectors */
465:     VecAXPY(vfull, alpha, vreduced);
466:   } else {
467:     PetscScalar       *y;
468:     const PetscScalar *x;
469:     PetscInt           i, n, m, rstart, rend;
470:     const PetscInt    *id;

472:     VecGetArray(vfull, &y);
473:     VecGetArrayRead(vreduced, &x);
474:     ISGetIndices(is, &id);
475:     ISGetLocalSize(is, &n);
476:     VecGetLocalSize(vreduced, &m);
478:     VecGetOwnershipRange(vfull, &rstart, &rend);
479:     y -= rstart;
480:     if (alpha == 1.0) {
481:       for (i = 0; i < n; ++i) {
482:         if (id[i] < 0) continue;
484:         y[id[i]] += x[i];
485:       }
486:     } else {
487:       for (i = 0; i < n; ++i) {
488:         if (id[i] < 0) continue;
490:         y[id[i]] += alpha * x[i];
491:       }
492:     }
493:     y += rstart;
494:     ISRestoreIndices(is, &id);
495:     VecRestoreArray(vfull, &y);
496:     VecRestoreArrayRead(vreduced, &x);
497:   }
498:   return 0;
499: }

501: /*@
502:   VecISCopy - Copies between a reduced vector and the appropriate elements of a full-space vector.

504:   Input Parameters:
505: + vfull    - the full-space vector
506: . is       - the index set for the reduced space
507: . mode     - the direction of copying, SCATTER_FORWARD or SCATTER_REVERSE
508: - vreduced - the reduced-space vector

510:   Output Parameters:
511: . vfull    - the sum of the full-space vector and reduced-space vector

513:   Notes:
514:     The index set identifies entries in the global vector.
515:     Negative indices are skipped; indices outside the ownership range of vfull will raise an error.

517:     mode == SCATTER_FORWARD: vfull[is[i]] = vreduced[i]
518:     mode == SCATTER_REVERSE: vreduced[i] = vfull[is[i]]

520:   Level: advanced

522: .seealso: `VecISSet()`, `VecISAXPY()`, `VecCopy()`
523: @*/
524: PetscErrorCode VecISCopy(Vec vfull, IS is, ScatterMode mode, Vec vreduced)
525: {
526:   PetscInt nfull, nreduced;

531:   VecGetSize(vfull, &nfull);
532:   VecGetSize(vreduced, &nreduced);

534:   if (nfull == nreduced) { /* Also takes care of masked vectors */
535:     if (mode == SCATTER_FORWARD) {
536:       VecCopy(vreduced, vfull);
537:     } else {
538:       VecCopy(vfull, vreduced);
539:     }
540:   } else {
541:     const PetscInt *id;
542:     PetscInt        i, n, m, rstart, rend;

544:     ISGetIndices(is, &id);
545:     ISGetLocalSize(is, &n);
546:     VecGetLocalSize(vreduced, &m);
547:     VecGetOwnershipRange(vfull, &rstart, &rend);
549:     if (mode == SCATTER_FORWARD) {
550:       PetscScalar       *y;
551:       const PetscScalar *x;

553:       VecGetArray(vfull, &y);
554:       VecGetArrayRead(vreduced, &x);
555:       y -= rstart;
556:       for (i = 0; i < n; ++i) {
557:         if (id[i] < 0) continue;
559:         y[id[i]] = x[i];
560:       }
561:       y += rstart;
562:       VecRestoreArrayRead(vreduced, &x);
563:       VecRestoreArray(vfull, &y);
564:     } else if (mode == SCATTER_REVERSE) {
565:       PetscScalar       *x;
566:       const PetscScalar *y;

568:       VecGetArrayRead(vfull, &y);
569:       VecGetArray(vreduced, &x);
570:       for (i = 0; i < n; ++i) {
571:         if (id[i] < 0) continue;
573:         x[i] = y[id[i] - rstart];
574:       }
575:       VecRestoreArray(vreduced, &x);
576:       VecRestoreArrayRead(vfull, &y);
577:     } else SETERRQ(PetscObjectComm((PetscObject)vfull), PETSC_ERR_ARG_WRONG, "Only forward or reverse modes are legal");
578:     ISRestoreIndices(is, &id);
579:   }
580:   return 0;
581: }

583: /*@
584:    ISComplementVec - Creates the complement of the index set relative to a layout defined by a Vec

586:    Collective on IS

588:    Input Parameters:
589: +  S -  a PETSc IS
590: -  V - the reference vector space

592:    Output Parameter:
593: .  T -  the complement of S

595:    Level: advanced

597: .seealso: `ISCreateGeneral()`
598: @*/
599: PetscErrorCode ISComplementVec(IS S, Vec V, IS *T)
600: {
601:   PetscInt start, end;

603:   VecGetOwnershipRange(V, &start, &end);
604:   ISComplement(S, start, end, T);
605:   return 0;
606: }

608: /*@
609:    VecISSet - Sets the elements of a vector, specified by an index set, to a constant

611:    Input Parameters:
612: +  V - the vector
613: .  S - index set for the locations in the vector
614: -  c - the constant

616:   Notes:
617:     The index set identifies entries in the global vector.
618:     Negative indices are skipped; indices outside the ownership range of V will raise an error.

620:    Level: advanced

622: .seealso: `VecISCopy()`, `VecISAXPY()`, `VecSet()`
623: @*/
624: PetscErrorCode VecISSet(Vec V, IS S, PetscScalar c)
625: {
626:   PetscInt        nloc, low, high, i;
627:   const PetscInt *s;
628:   PetscScalar    *v;

630:   if (!S) return 0; /* simply return with no-op if the index set is NULL */

635:   VecGetOwnershipRange(V, &low, &high);
636:   ISGetLocalSize(S, &nloc);
637:   ISGetIndices(S, &s);
638:   VecGetArray(V, &v);
639:   for (i = 0; i < nloc; ++i) {
640:     if (s[i] < 0) continue;
642:     v[s[i] - low] = c;
643:   }
644:   ISRestoreIndices(S, &s);
645:   VecRestoreArray(V, &v);
646:   return 0;
647: }

649: #if !defined(PETSC_USE_COMPLEX)
650: /*@C
651:   VecBoundGradientProjection - Projects  vector according to this definition.
652:   If XL[i] < X[i] < XU[i], then GP[i] = G[i];
653:   If X[i] <= XL[i], then GP[i] = min(G[i],0);
654:   If X[i] >= XU[i], then GP[i] = max(G[i],0);

656:   Input Parameters:
657: + G - current gradient vector
658: . X - current solution vector with XL[i] <= X[i] <= XU[i]
659: . XL - lower bounds
660: - XU - upper bounds

662:   Output Parameter:
663: . GP - gradient projection vector

665:   Notes:
666:     GP may be the same vector as G

668:   Level: advanced
669: @*/
670: PetscErrorCode VecBoundGradientProjection(Vec G, Vec X, Vec XL, Vec XU, Vec GP)
671: {
672:   PetscInt         n, i;
673:   const PetscReal *xptr, *xlptr, *xuptr;
674:   PetscReal       *gptr, *gpptr;
675:   PetscReal        xval, gpval;

677:   /* Project variables at the lower and upper bound */
683:   if (!XL && !XU) {
684:     VecCopy(G, GP);
685:     return 0;
686:   }

688:   VecGetLocalSize(X, &n);

690:   VecGetArrayRead(X, &xptr);
691:   VecGetArrayRead(XL, &xlptr);
692:   VecGetArrayRead(XU, &xuptr);
693:   VecGetArrayPair(G, GP, &gptr, &gpptr);

695:   for (i = 0; i < n; ++i) {
696:     gpval = gptr[i];
697:     xval  = xptr[i];
698:     if (gpval > 0.0 && xval <= xlptr[i]) {
699:       gpval = 0.0;
700:     } else if (gpval < 0.0 && xval >= xuptr[i]) {
701:       gpval = 0.0;
702:     }
703:     gpptr[i] = gpval;
704:   }

706:   VecRestoreArrayRead(X, &xptr);
707:   VecRestoreArrayRead(XL, &xlptr);
708:   VecRestoreArrayRead(XU, &xuptr);
709:   VecRestoreArrayPair(G, GP, &gptr, &gpptr);
710:   return 0;
711: }
712: #endif

714: /*@
715:      VecStepMaxBounded - See below

717:      Collective on Vec

719:      Input Parameters:
720: +      X  - vector with no negative entries
721: .      XL - lower bounds
722: .      XU - upper bounds
723: -      DX  - step direction, can have negative, positive or zero entries

725:      Output Parameter:
726: .     stepmax -   minimum value so that X[i] + stepmax*DX[i] <= XL[i]  or  XU[i] <= X[i] + stepmax*DX[i]

728:   Level: intermediate

730: @*/
731: PetscErrorCode VecStepMaxBounded(Vec X, Vec DX, Vec XL, Vec XU, PetscReal *stepmax)
732: {
733:   PetscInt           i, nn;
734:   const PetscScalar *xx, *dx, *xl, *xu;
735:   PetscReal          localmax = 0;


742:   VecGetArrayRead(X, &xx);
743:   VecGetArrayRead(XL, &xl);
744:   VecGetArrayRead(XU, &xu);
745:   VecGetArrayRead(DX, &dx);
746:   VecGetLocalSize(X, &nn);
747:   for (i = 0; i < nn; i++) {
748:     if (PetscRealPart(dx[i]) > 0) {
749:       localmax = PetscMax(localmax, PetscRealPart((xu[i] - xx[i]) / dx[i]));
750:     } else if (PetscRealPart(dx[i]) < 0) {
751:       localmax = PetscMax(localmax, PetscRealPart((xl[i] - xx[i]) / dx[i]));
752:     }
753:   }
754:   VecRestoreArrayRead(X, &xx);
755:   VecRestoreArrayRead(XL, &xl);
756:   VecRestoreArrayRead(XU, &xu);
757:   VecRestoreArrayRead(DX, &dx);
758:   MPIU_Allreduce(&localmax, stepmax, 1, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)X));
759:   return 0;
760: }

762: /*@
763:      VecStepBoundInfo - See below

765:      Collective on Vec

767:      Input Parameters:
768: +      X  - vector with no negative entries
769: .      XL - lower bounds
770: .      XU - upper bounds
771: -      DX  - step direction, can have negative, positive or zero entries

773:      Output Parameters:
774: +     boundmin -  (may be NULL this it is not computed) maximum value so that   XL[i] <= X[i] + boundmax*DX[i] <= XU[i]
775: .     wolfemin -  (may be NULL this it is not computed)
776: -     boundmax -   (may be NULL this it is not computed) minimum value so that X[i] + boundmax*DX[i] <= XL[i]  or  XU[i] <= X[i] + boundmax*DX[i]

778:      Notes:
779:     For complex numbers only compares the real part

781:   Level: advanced
782: @*/
783: PetscErrorCode VecStepBoundInfo(Vec X, Vec DX, Vec XL, Vec XU, PetscReal *boundmin, PetscReal *wolfemin, PetscReal *boundmax)
784: {
785:   PetscInt           n, i;
786:   const PetscScalar *x, *xl, *xu, *dx;
787:   PetscReal          t;
788:   PetscReal          localmin = PETSC_INFINITY, localwolfemin = PETSC_INFINITY, localmax = -1;
789:   MPI_Comm           comm;


796:   VecGetArrayRead(X, &x);
797:   VecGetArrayRead(XL, &xl);
798:   VecGetArrayRead(XU, &xu);
799:   VecGetArrayRead(DX, &dx);
800:   VecGetLocalSize(X, &n);
801:   for (i = 0; i < n; ++i) {
802:     if (PetscRealPart(dx[i]) > 0 && PetscRealPart(xu[i]) < PETSC_INFINITY) {
803:       t        = PetscRealPart((xu[i] - x[i]) / dx[i]);
804:       localmin = PetscMin(t, localmin);
805:       if (localmin > 0) localwolfemin = PetscMin(t, localwolfemin);
806:       localmax = PetscMax(t, localmax);
807:     } else if (PetscRealPart(dx[i]) < 0 && PetscRealPart(xl[i]) > PETSC_NINFINITY) {
808:       t        = PetscRealPart((xl[i] - x[i]) / dx[i]);
809:       localmin = PetscMin(t, localmin);
810:       if (localmin > 0) localwolfemin = PetscMin(t, localwolfemin);
811:       localmax = PetscMax(t, localmax);
812:     }
813:   }

815:   VecRestoreArrayRead(X, &x);
816:   VecRestoreArrayRead(XL, &xl);
817:   VecRestoreArrayRead(XU, &xu);
818:   VecRestoreArrayRead(DX, &dx);
819:   PetscObjectGetComm((PetscObject)X, &comm);

821:   if (boundmin) {
822:     MPIU_Allreduce(&localmin, boundmin, 1, MPIU_REAL, MPIU_MIN, comm);
823:     PetscInfo(X, "Step Bound Info: Closest Bound: %20.19e\n", (double)*boundmin);
824:   }
825:   if (wolfemin) {
826:     MPIU_Allreduce(&localwolfemin, wolfemin, 1, MPIU_REAL, MPIU_MIN, comm);
827:     PetscInfo(X, "Step Bound Info: Wolfe: %20.19e\n", (double)*wolfemin);
828:   }
829:   if (boundmax) {
830:     MPIU_Allreduce(&localmax, boundmax, 1, MPIU_REAL, MPIU_MAX, comm);
831:     if (*boundmax < 0) *boundmax = PETSC_INFINITY;
832:     PetscInfo(X, "Step Bound Info: Max: %20.19e\n", (double)*boundmax);
833:   }
834:   return 0;
835: }

837: /*@
838:      VecStepMax - Returns the largest value so that x[i] + step*DX[i] >= 0 for all i

840:      Collective on Vec

842:      Input Parameters:
843: +      X  - vector with no negative entries
844: -      DX  - a step direction, can have negative, positive or zero entries

846:      Output Parameter:
847: .    step - largest value such that x[i] + step*DX[i] >= 0 for all i

849:      Notes:
850:     For complex numbers only compares the real part

852:   Level: advanced
853:  @*/
854: PetscErrorCode VecStepMax(Vec X, Vec DX, PetscReal *step)
855: {
856:   PetscInt           i, nn;
857:   PetscReal          stepmax = PETSC_INFINITY;
858:   const PetscScalar *xx, *dx;


863:   VecGetLocalSize(X, &nn);
864:   VecGetArrayRead(X, &xx);
865:   VecGetArrayRead(DX, &dx);
866:   for (i = 0; i < nn; ++i) {
868:     if (PetscRealPart(dx[i]) < 0) stepmax = PetscMin(stepmax, PetscRealPart(-xx[i] / dx[i]));
869:   }
870:   VecRestoreArrayRead(X, &xx);
871:   VecRestoreArrayRead(DX, &dx);
872:   MPIU_Allreduce(&stepmax, step, 1, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)X));
873:   return 0;
874: }

876: /*@
877:   VecPow - Replaces each component of a vector by x_i^p

879:   Logically Collective on v

881:   Input Parameters:
882: + v - the vector
883: - p - the exponent to use on each element

885:   Level: intermediate

887: @*/
888: PetscErrorCode VecPow(Vec v, PetscScalar p)
889: {
890:   PetscInt     n, i;
891:   PetscScalar *v1;


895:   VecGetArray(v, &v1);
896:   VecGetLocalSize(v, &n);

898:   if (1.0 == p) {
899:   } else if (-1.0 == p) {
900:     for (i = 0; i < n; ++i) v1[i] = 1.0 / v1[i];
901:   } else if (0.0 == p) {
902:     for (i = 0; i < n; ++i) {
903:       /*  Not-a-number left alone
904:           Infinity set to one  */
905:       if (v1[i] == v1[i]) v1[i] = 1.0;
906:     }
907:   } else if (0.5 == p) {
908:     for (i = 0; i < n; ++i) {
909:       if (PetscRealPart(v1[i]) >= 0) {
910:         v1[i] = PetscSqrtScalar(v1[i]);
911:       } else {
912:         v1[i] = PETSC_INFINITY;
913:       }
914:     }
915:   } else if (-0.5 == p) {
916:     for (i = 0; i < n; ++i) {
917:       if (PetscRealPart(v1[i]) >= 0) {
918:         v1[i] = 1.0 / PetscSqrtScalar(v1[i]);
919:       } else {
920:         v1[i] = PETSC_INFINITY;
921:       }
922:     }
923:   } else if (2.0 == p) {
924:     for (i = 0; i < n; ++i) v1[i] *= v1[i];
925:   } else if (-2.0 == p) {
926:     for (i = 0; i < n; ++i) v1[i] = 1.0 / (v1[i] * v1[i]);
927:   } else {
928:     for (i = 0; i < n; ++i) {
929:       if (PetscRealPart(v1[i]) >= 0) {
930:         v1[i] = PetscPowScalar(v1[i], p);
931:       } else {
932:         v1[i] = PETSC_INFINITY;
933:       }
934:     }
935:   }
936:   VecRestoreArray(v, &v1);
937:   return 0;
938: }

940: /*@
941:   VecMedian - Computes the componentwise median of three vectors
942:   and stores the result in this vector.  Used primarily for projecting
943:   a vector within upper and lower bounds.

945:   Logically Collective

947:   Input Parameters:
948: + Vec1 - The first vector
949: . Vec2 - The second vector
950: - Vec3 - The third vector

952:   Output Parameter:
953: . VMedian - The median vector (this can be any one of the input vectors)

955:   Level: advanced
956: @*/
957: PetscErrorCode VecMedian(Vec Vec1, Vec Vec2, Vec Vec3, Vec VMedian)
958: {
959:   PetscInt           i, n, low1, high1;
960:   const PetscScalar *v1, *v2, *v3;
961:   PetscScalar       *vmed;


968:   if (!Vec1 && !Vec3) {
969:     VecCopy(Vec2, VMedian);
970:     return 0;
971:   }
972:   if (Vec1 == Vec2 || Vec1 == Vec3) {
973:     VecCopy(Vec1, VMedian);
974:     return 0;
975:   }
976:   if (Vec2 == Vec3) {
977:     VecCopy(Vec2, VMedian);
978:     return 0;
979:   }

981:   /* Assert that Vec1 != Vec2 and Vec2 != Vec3 */
992:   VecCheckSameSize(Vec1, 1, Vec2, 2);
993:   VecCheckSameSize(Vec1, 1, Vec3, 3);
994:   VecCheckSameSize(Vec1, 1, VMedian, 4);

996:   VecGetOwnershipRange(Vec1, &low1, &high1);
997:   VecGetLocalSize(Vec1, &n);
998:   if (n > 0) {
999:     VecGetArray(VMedian, &vmed);
1000:     if (Vec1 != VMedian) {
1001:       VecGetArrayRead(Vec1, &v1);
1002:     } else {
1003:       v1 = vmed;
1004:     }
1005:     if (Vec2 != VMedian) {
1006:       VecGetArrayRead(Vec2, &v2);
1007:     } else {
1008:       v2 = vmed;
1009:     }
1010:     if (Vec3 != VMedian) {
1011:       VecGetArrayRead(Vec3, &v3);
1012:     } else {
1013:       v3 = vmed;
1014:     }

1016:     for (i = 0; i < n; ++i) vmed[i] = PetscMax(PetscMax(PetscMin(PetscRealPart(v1[i]), PetscRealPart(v2[i])), PetscMin(PetscRealPart(v1[i]), PetscRealPart(v3[i]))), PetscMin(PetscRealPart(v2[i]), PetscRealPart(v3[i])));

1018:     VecRestoreArray(VMedian, &vmed);
1019:     if (VMedian != Vec1) VecRestoreArrayRead(Vec1, &v1);
1020:     if (VMedian != Vec2) VecRestoreArrayRead(Vec2, &v2);
1021:     if (VMedian != Vec3) VecRestoreArrayRead(Vec3, &v3);
1022:   }
1023:   return 0;
1024: }