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:   PetscErrorCode    ierr;
 23:   PetscInt          i,n_same=0;
 24:   PetscInt          n,low,high;
 25:   PetscInt          *same=NULL;
 26:   const PetscScalar *v1,*v2;

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

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

 45:     PetscMalloc1(n,&same);

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

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

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

 66:   Collective on S

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

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

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

 77:   For complex numbers this only compares the real part

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

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

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

106:     PetscMalloc1(n,&lt);

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

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

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

127:   Collective on S

129:   Input Parameters:
130: . Vec1, Vec2 - the two vectors to compare

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

135:   Notes:
136:   The two vectors must have the same parallel layout

138:   For complex numbers this only compares the real part

140:   Level: advanced
141: @*/
142: PetscErrorCode VecWhichGreaterThan(Vec Vec1, Vec Vec2, IS *S)
143: {
144:   PetscErrorCode    ierr;
145:   PetscInt          i,n_gt=0;
146:   PetscInt          n,low,high;
147:   PetscInt          *gt=NULL;
148:   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])) {gt[n_gt]=low+i; ++n_gt;}
171:     }

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

184: /*@
185:   VecWhichBetween - Creates an index set containing the indices
186:                where  VecLow < V < VecHigh

188:   Collective on S

190:   Input Parameters:
191: + VecLow - lower bound
192: . V - Vector to compare
193: - VecHigh - higher bound

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

198:   Notes:
199:   The vectors must have the same parallel layout

201:   For complex numbers this only compares the real part

203:   Level: advanced
204: @*/
205: PetscErrorCode VecWhichBetween(Vec VecLow, Vec V, Vec VecHigh, IS *S)
206: {

208:   PetscErrorCode    ierr;
209:   PetscInt          i,n_vm=0;
210:   PetscInt          n,low,high;
211:   PetscInt          *vm=NULL;
212:   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])) {vm[n_vm]=low+i; ++n_vm;}
244:     }

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

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

262:   Collective on S

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

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

272:   Level: advanced
273: @*/

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

289:   VecCheckSameSize(V,2,VecLow,1);
290:   VecCheckSameSize(V,2,VecHigh,3);

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

309:     PetscMalloc1(n,&vm);

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

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

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

334:   Collective on S

336:   Input Parameters:
337: + VecLow - lower bound
338: . V - Vector to compare
339: . D - Direction to compare
340: . VecHigh - higher bound
341: - Strong - indicator for applying strongly inactive test

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

346:   Level: advanced
347: @*/

349: PetscErrorCode VecWhichInactive(Vec VecLow, Vec V, Vec D, Vec VecHigh, PetscBool Strong, IS * S)
350: {
351:   PetscErrorCode    ierr;
352:   PetscInt          i,n_vm=0;
353:   PetscInt          n,low,high;
354:   PetscInt          *vm=NULL;
355:   const PetscScalar *v1,*v2,*v,*d;

365:   VecCheckSameSize(V,2,VecLow,1);
366:   VecCheckSameSize(V,2,D,3);
367:   VecCheckSameSize(V,2,VecHigh,4);

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

395:     PetscMalloc1(n,&vm);

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

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

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

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

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

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

451:   Level: advanced

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

464:   VecGetSize(vfull,&nfull);
465:   VecGetSize(vreduced,&nreduced);

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

475:     VecGetArray(vfull,&y);
476:     VecGetArrayRead(vreduced,&x);
477:     ISGetIndices(is,&id);
478:     ISGetLocalSize(is,&n);
479:     VecGetLocalSize(vreduced,&m);
480:     if (m != n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"IS local length not equal to Vec local length");
481:     VecGetOwnershipRange(vfull,&rstart,&rend);
482:     y   -= rstart;
483:     if (alpha == 1.0) {
484:       for (i=0; i<n; ++i) {
485:         if (id[i] < 0) continue;
486:         if (id[i] < rstart || id[i] >= rend) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only owned values supported");
487:         y[id[i]] += x[i];
488:       }
489:     } else {
490:       for (i=0; i<n; ++i) {
491:         if (id[i] < 0) continue;
492:         if (id[i] < rstart || id[i] >= rend) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only owned values supported");
493:         y[id[i]] += alpha*x[i];
494:       }
495:     }
496:     y += rstart;
497:     ISRestoreIndices(is,&id);
498:     VecRestoreArray(vfull,&y);
499:     VecRestoreArrayRead(vreduced,&x);
500:   }
501:   return(0);
502: }

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

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

513:   Output Parameters:
514: . vfull    - the sum of the full-space vector and reduced-space vector

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

520:     mode == SCATTER_FORWARD: vfull[is[i]] = vreduced[i]
521:     mode == SCATTER_REVERSE: vreduced[i] = vfull[is[i]]

523:   Level: advanced

525: .seealso:  VecAXPY(), VecGetOwnershipRange()
526: @*/
527: PetscErrorCode VecISCopy(Vec vfull, IS is, ScatterMode mode, Vec vreduced)
528: {
529:   PetscInt       nfull, nreduced;

536:   VecGetSize(vfull, &nfull);
537:   VecGetSize(vreduced, &nreduced);

539:   if (nfull == nreduced) { /* Also takes care of masked vectors */
540:     if (mode == SCATTER_FORWARD) {
541:       VecCopy(vreduced, vfull);
542:     } else {
543:       VecCopy(vfull, vreduced);
544:     }
545:   } else {
546:     const PetscInt *id;
547:     PetscInt        i, n, m, rstart, rend;

549:     ISGetIndices(is, &id);
550:     ISGetLocalSize(is, &n);
551:     VecGetLocalSize(vreduced, &m);
552:     VecGetOwnershipRange(vfull, &rstart, &rend);
553:     if (m != n) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "IS local length %D not equal to Vec local length %D", n, m);
554:     if (mode == SCATTER_FORWARD) {
555:       PetscScalar       *y;
556:       const PetscScalar *x;

558:       VecGetArray(vfull, &y);
559:       VecGetArrayRead(vreduced, &x);
560:       y   -= rstart;
561:       for (i = 0; i < n; ++i) {
562:         if (id[i] < 0) continue;
563:         if (id[i] < rstart || id[i] >= rend) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only owned values supported");
564:         y[id[i]] = x[i];
565:       }
566:       y   += rstart;
567:       VecRestoreArrayRead(vreduced, &x);
568:       VecRestoreArray(vfull, &y);
569:     } else if (mode == SCATTER_REVERSE) {
570:       PetscScalar       *x;
571:       const PetscScalar *y;

573:       VecGetArrayRead(vfull, &y);
574:       VecGetArray(vreduced, &x);
575:       for (i = 0; i < n; ++i) {
576:         if (id[i] < 0) continue;
577:         if (id[i] < rstart || id[i] >= rend) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only owned values supported");
578:         x[i] = y[id[i]-rstart];
579:       }
580:       VecRestoreArray(vreduced, &x);
581:       VecRestoreArrayRead(vfull, &y);
582:     } else SETERRQ(PetscObjectComm((PetscObject) vfull), PETSC_ERR_ARG_WRONG, "Only forward or reverse modes are legal");
583:     ISRestoreIndices(is, &id);
584:   }
585:   return(0);
586: }

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

591:    Collective on IS

593:    Input Parameters:
594: +  S -  a PETSc IS
595: -  V - the reference vector space

597:    Output Parameter:
598: .  T -  the complement of S

600:    Level: advanced

602: .seealso: ISCreateGeneral()
603: @*/
604: PetscErrorCode ISComplementVec(IS S, Vec V, IS *T)
605: {
607:   PetscInt       start, end;

610:   VecGetOwnershipRange(V,&start,&end);
611:   ISComplement(S,start,end,T);
612:   return(0);
613: }

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

618:    Input Parameters:
619: +  V - the vector
620: .  S - index set for the locations in the vector
621: -  c - the constant

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

627:    Level: advanced

629: .seealso: VecSet(), VecGetOwnershipRange()
630: @*/
631: PetscErrorCode VecISSet(Vec V,IS S, PetscScalar c)
632: {
634:   PetscInt       nloc,low,high,i;
635:   const PetscInt *s;
636:   PetscScalar    *v;

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

644:   VecGetOwnershipRange(V,&low,&high);
645:   ISGetLocalSize(S,&nloc);
646:   ISGetIndices(S,&s);
647:   VecGetArray(V,&v);
648:   for (i=0; i<nloc; ++i) {
649:     if (s[i] < 0) continue;
650:     if (s[i] < low || s[i] >= high) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only owned values supported");
651:     v[s[i]-low] = c;
652:   }
653:   ISRestoreIndices(S,&s);
654:   VecRestoreArray(V,&v);
655:   return(0);
656: }

658: #if !defined(PETSC_USE_COMPLEX)
659: /*@C
660:   VecBoundGradientProjection - Projects  vector according to this definition.
661:   If XL[i] < X[i] < XU[i], then GP[i] = G[i];
662:   If X[i] <= XL[i], then GP[i] = min(G[i],0);
663:   If X[i] >= XU[i], then GP[i] = max(G[i],0);

665:   Input Parameters:
666: + G - current gradient vector
667: . X - current solution vector with XL[i] <= X[i] <= XU[i]
668: . XL - lower bounds
669: - XU - upper bounds

671:   Output Parameter:
672: . GP - gradient projection vector

674:   Notes:
675:     GP may be the same vector as G

677:   Level: advanced
678: @*/
679: PetscErrorCode VecBoundGradientProjection(Vec G, Vec X, Vec XL, Vec XU, Vec GP)
680: {

682:   PetscErrorCode  ierr;
683:   PetscInt        n,i;
684:   const PetscReal *xptr,*xlptr,*xuptr;
685:   PetscReal       *gptr,*gpptr;
686:   PetscReal       xval,gpval;

688:   /* Project variables at the lower and upper bound */

696:   VecGetLocalSize(X,&n);

698:   VecGetArrayRead(X,&xptr);
699:   VecGetArrayRead(XL,&xlptr);
700:   VecGetArrayRead(XU,&xuptr);
701:   VecGetArrayPair(G,GP,&gptr,&gpptr);

703:   for (i=0; i<n; ++i) {
704:     gpval = gptr[i]; xval = xptr[i];
705:     if (gpval>0.0 && xval<=xlptr[i]) {
706:       gpval = 0.0;
707:     } else if (gpval<0.0 && xval>=xuptr[i]) {
708:       gpval = 0.0;
709:     }
710:     gpptr[i] = gpval;
711:   }

713:   VecRestoreArrayRead(X,&xptr);
714:   VecRestoreArrayRead(XL,&xlptr);
715:   VecRestoreArrayRead(XU,&xuptr);
716:   VecRestoreArrayPair(G,GP,&gptr,&gpptr);
717:   return(0);
718: }
719: #endif

721: /*@
722:      VecStepMaxBounded - See below

724:      Collective on Vec

726:      Input Parameters:
727: +      X  - vector with no negative entries
728: .      XL - lower bounds
729: .      XU - upper bounds
730: -      DX  - step direction, can have negative, positive or zero entries

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

735:   Level: intermediate

737: @*/
738: PetscErrorCode VecStepMaxBounded(Vec X, Vec DX, Vec XL, Vec XU, PetscReal *stepmax)
739: {
740:   PetscErrorCode    ierr;
741:   PetscInt          i,nn;
742:   const PetscScalar *xx,*dx,*xl,*xu;
743:   PetscReal         localmax=0;


751:   VecGetArrayRead(X,&xx);
752:   VecGetArrayRead(XL,&xl);
753:   VecGetArrayRead(XU,&xu);
754:   VecGetArrayRead(DX,&dx);
755:   VecGetLocalSize(X,&nn);
756:   for (i=0;i<nn;i++) {
757:     if (PetscRealPart(dx[i]) > 0) {
758:       localmax=PetscMax(localmax,PetscRealPart((xu[i]-xx[i])/dx[i]));
759:     } else if (PetscRealPart(dx[i])<0) {
760:       localmax=PetscMax(localmax,PetscRealPart((xl[i]-xx[i])/dx[i]));
761:     }
762:   }
763:   VecRestoreArrayRead(X,&xx);
764:   VecRestoreArrayRead(XL,&xl);
765:   VecRestoreArrayRead(XU,&xu);
766:   VecRestoreArrayRead(DX,&dx);
767:   MPIU_Allreduce(&localmax,stepmax,1,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)X));
768:   return(0);
769: }

771: /*@
772:      VecStepBoundInfo - See below

774:      Collective on Vec

776:      Input Parameters:
777: +      X  - vector with no negative entries
778: .      XL - lower bounds
779: .      XU - upper bounds
780: -      DX  - step direction, can have negative, positive or zero entries

782:      Output Parameters:
783: +     boundmin -  (may be NULL this it is not computed) maximum value so that   XL[i] <= X[i] + boundmax*DX[i] <= XU[i]
784: .     wolfemin -  (may be NULL this it is not computed)
785: -     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]

787:      Notes:
788:     For complex numbers only compares the real part

790:   Level: advanced
791: @*/
792: PetscErrorCode VecStepBoundInfo(Vec X, Vec DX, Vec XL, Vec XU, PetscReal *boundmin, PetscReal *wolfemin, PetscReal *boundmax)
793: {
794:   PetscErrorCode    ierr;
795:   PetscInt          n,i;
796:   const PetscScalar *x,*xl,*xu,*dx;
797:   PetscReal         t;
798:   PetscReal         localmin=PETSC_INFINITY,localwolfemin=PETSC_INFINITY,localmax=-1;
799:   MPI_Comm          comm;


807:   VecGetArrayRead(X,&x);
808:   VecGetArrayRead(XL,&xl);
809:   VecGetArrayRead(XU,&xu);
810:   VecGetArrayRead(DX,&dx);
811:   VecGetLocalSize(X,&n);
812:   for (i=0; i<n; ++i) {
813:     if (PetscRealPart(dx[i])>0 && PetscRealPart(xu[i]) < PETSC_INFINITY) {
814:       t=PetscRealPart((xu[i]-x[i])/dx[i]);
815:       localmin=PetscMin(t,localmin);
816:       if (localmin>0) {
817:         localwolfemin = PetscMin(t,localwolfemin);
818:       }
819:       localmax = PetscMax(t,localmax);
820:     } else if (PetscRealPart(dx[i])<0 && PetscRealPart(xl[i]) > PETSC_NINFINITY) {
821:       t=PetscRealPart((xl[i]-x[i])/dx[i]);
822:       localmin = PetscMin(t,localmin);
823:       if (localmin>0) {
824:         localwolfemin = PetscMin(t,localwolfemin);
825:       }
826:       localmax = PetscMax(t,localmax);
827:     }
828:   }

830:   VecRestoreArrayRead(X,&x);
831:   VecRestoreArrayRead(XL,&xl);
832:   VecRestoreArrayRead(XU,&xu);
833:   VecRestoreArrayRead(DX,&dx);
834:   PetscObjectGetComm((PetscObject)X,&comm);

836:   if (boundmin) {
837:     MPIU_Allreduce(&localmin,boundmin,1,MPIU_REAL,MPIU_MIN,comm);
838:     PetscInfo1(X,"Step Bound Info: Closest Bound: %20.19e\n",(double)*boundmin);
839:   }
840:   if (wolfemin) {
841:     MPIU_Allreduce(&localwolfemin,wolfemin,1,MPIU_REAL,MPIU_MIN,comm);
842:     PetscInfo1(X,"Step Bound Info: Wolfe: %20.19e\n",(double)*wolfemin);
843:   }
844:   if (boundmax) {
845:     MPIU_Allreduce(&localmax,boundmax,1,MPIU_REAL,MPIU_MAX,comm);
846:     if (*boundmax < 0) *boundmax=PETSC_INFINITY;
847:     PetscInfo1(X,"Step Bound Info: Max: %20.19e\n",(double)*boundmax);
848:   }
849:   return(0);
850: }

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

855:      Collective on Vec

857:      Input Parameters:
858: +      X  - vector with no negative entries
859: -      DX  - a step direction, can have negative, positive or zero entries

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

864:      Notes:
865:     For complex numbers only compares the real part

867:   Level: advanced
868:  @*/
869: PetscErrorCode VecStepMax(Vec X, Vec DX, PetscReal *step)
870: {
871:   PetscErrorCode    ierr;
872:   PetscInt          i,nn;
873:   PetscReal         stepmax=PETSC_INFINITY;
874:   const PetscScalar *xx,*dx;


880:   VecGetLocalSize(X,&nn);
881:   VecGetArrayRead(X,&xx);
882:   VecGetArrayRead(DX,&dx);
883:   for (i=0;i<nn;++i) {
884:     if (PetscRealPart(xx[i]) < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Vector must be positive");
885:     else if (PetscRealPart(dx[i])<0) stepmax=PetscMin(stepmax,PetscRealPart(-xx[i]/dx[i]));
886:   }
887:   VecRestoreArrayRead(X,&xx);
888:   VecRestoreArrayRead(DX,&dx);
889:   MPIU_Allreduce(&stepmax,step,1,MPIU_REAL,MPIU_MIN,PetscObjectComm((PetscObject)X));
890:   return(0);
891: }

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

896:   Logically Collective on v

898:   Input Parameters:
899: + v - the vector
900: - p - the exponent to use on each element

902:   Level: intermediate

904: @*/
905: PetscErrorCode VecPow(Vec v, PetscScalar p)
906: {
908:   PetscInt       n,i;
909:   PetscScalar    *v1;


914:   VecGetArray(v,&v1);
915:   VecGetLocalSize(v,&n);

917:   if (1.0 == p) {
918:   } else if (-1.0 == p) {
919:     for (i = 0; i < n; ++i) {
920:       v1[i] = 1.0 / v1[i];
921:     }
922:   } else if (0.0 == p) {
923:     for (i = 0; i < n; ++i) {
924:       /*  Not-a-number left alone
925:           Infinity set to one  */
926:       if (v1[i] == v1[i]) {
927:         v1[i] = 1.0;
928:       }
929:     }
930:   } else if (0.5 == p) {
931:     for (i = 0; i < n; ++i) {
932:       if (PetscRealPart(v1[i]) >= 0) {
933:         v1[i] = PetscSqrtScalar(v1[i]);
934:       } else {
935:         v1[i] = PETSC_INFINITY;
936:       }
937:     }
938:   } else if (-0.5 == p) {
939:     for (i = 0; i < n; ++i) {
940:       if (PetscRealPart(v1[i]) >= 0) {
941:         v1[i] = 1.0 / PetscSqrtScalar(v1[i]);
942:       } else {
943:         v1[i] = PETSC_INFINITY;
944:       }
945:     }
946:   } else if (2.0 == p) {
947:     for (i = 0; i < n; ++i) {
948:       v1[i] *= v1[i];
949:     }
950:   } else if (-2.0 == p) {
951:     for (i = 0; i < n; ++i) {
952:       v1[i] = 1.0 / (v1[i] * v1[i]);
953:     }
954:   } else {
955:     for (i = 0; i < n; ++i) {
956:       if (PetscRealPart(v1[i]) >= 0) {
957:         v1[i] = PetscPowScalar(v1[i],p);
958:       } else {
959:         v1[i] = PETSC_INFINITY;
960:       }
961:     }
962:   }
963:   VecRestoreArray(v,&v1);
964:   return(0);
965: }

967: /*@
968:   VecMedian - Computes the componentwise median of three vectors
969:   and stores the result in this vector.  Used primarily for projecting
970:   a vector within upper and lower bounds.

972:   Logically Collective

974:   Input Parameters:
975: + Vec1 - The first vector
976: . Vec2 - The second vector
977: - Vec3 - The third vector

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

982:   Level: advanced
983: @*/
984: PetscErrorCode VecMedian(Vec Vec1, Vec Vec2, Vec Vec3, Vec VMedian)
985: {
986:   PetscErrorCode    ierr;
987:   PetscInt          i,n,low1,high1;
988:   const PetscScalar *v1,*v2,*v3;
989:   PetscScalar       *vmed;


997:   if (Vec1==Vec2 || Vec1==Vec3) {
998:     VecCopy(Vec1,VMedian);
999:     return(0);
1000:   }
1001:   if (Vec2==Vec3) {
1002:     VecCopy(Vec2,VMedian);
1003:     return(0);
1004:   }

1006:   /* Assert that Vec1 != Vec2 and Vec2 != Vec3 */
1017:   VecCheckSameSize(Vec1,1,Vec2,2);
1018:   VecCheckSameSize(Vec1,1,Vec3,3);
1019:   VecCheckSameSize(Vec1,1,VMedian,4);

1021:   VecGetOwnershipRange(Vec1,&low1,&high1);
1022:   VecGetLocalSize(Vec1,&n);
1023:   if (n>0) {
1024:     VecGetArray(VMedian,&vmed);
1025:     if (Vec1 != VMedian) {
1026:       VecGetArrayRead(Vec1,&v1);
1027:     } else {
1028:       v1=vmed;
1029:     }
1030:     if (Vec2 != VMedian) {
1031:       VecGetArrayRead(Vec2,&v2);
1032:     } else {
1033:       v2=vmed;
1034:     }
1035:     if (Vec3 != VMedian) {
1036:       VecGetArrayRead(Vec3,&v3);
1037:     } else {
1038:       v3=vmed;
1039:     }

1041:     for (i=0;i<n;++i) {
1042:       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])));
1043:     }

1045:     VecRestoreArray(VMedian,&vmed);
1046:     if (VMedian != Vec1) {
1047:       VecRestoreArrayRead(Vec1,&v1);
1048:     }
1049:     if (VMedian != Vec2) {
1050:       VecRestoreArrayRead(Vec2,&v2);
1051:     }
1052:     if (VMedian != Vec3) {
1053:       VecRestoreArrayRead(Vec3,&v3);
1054:     }
1055:   }
1056:   return(0);
1057: }