Actual source code: projection.c

petsc-master 2018-02-21
Report Typos and Errors
  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: the two vectors must have the same parallel layout

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

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

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

 44:     PetscMalloc1(n,&same);

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

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

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

 65:   Collective on S

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

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

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

 76:   For complex numbers this only compares the real part

 78:   Level: advanced
 79: @*/
 80: PetscErrorCode VecWhichLessThan(Vec Vec1, Vec Vec2, IS *S)
 81: {
 82:   PetscErrorCode    ierr;
 83:   PetscInt          i,n_lt=0;
 84:   PetscInt          n,low,high;
 85:   PetscInt          *lt=NULL;
 86:   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])) {lt[n_lt]=low+i; ++n_lt;}
109:     }

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

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

126:   Collective on S

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

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

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

137:   For complex numbers this only compares the real part

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

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

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

166:     PetscMalloc1(n,&gt);

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

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

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

187:   Collective on S

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

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

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

200:   For complex numbers this only compares the real part

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

207:   PetscErrorCode    ierr;
208:   PetscInt          i,n_vm=0;
209:   PetscInt          n,low,high;
210:   PetscInt          *vm=NULL;
211:   const PetscScalar *v1,*v2,*vmiddle;

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

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

239:     PetscMalloc1(n,&vm);

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

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


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: }


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

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

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

448:   Level: advanced

450: .seealso:  VecAXPY()
451: @*/
452: PetscErrorCode VecISAXPY(Vec vfull, IS is, PetscScalar alpha, Vec vreduced)
453: {
454:   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;
470:     const PetscInt    *id;

472:     VecGetArray(vfull,&y);
473:     VecGetArrayRead(vreduced,&x);
474:     ISGetIndices(is,&id);
475:     ISGetLocalSize(is,&n);
476:     VecGetLocalSize(vreduced,&m);
477:     if (m != n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"IS local length not equal to Vec local length");
478:     VecGetOwnershipRange(vfull,&rstart,NULL);
479:     y -= rstart;
480:     if (alpha == 1.0) {
481:       for (i=0; i<n; ++i) {
482:         y[id[i]] += x[i];
483:       }
484:     } else {
485:       for (i=0; i<n; ++i) {
486:         y[id[i]] += alpha*x[i];
487:       }
488:     }
489:     y += rstart;
490:     ISRestoreIndices(is,&id);
491:     VecRestoreArray(vfull,&y);
492:     VecRestoreArrayRead(vreduced,&x);
493:   }
494:   return(0);
495: }

497: /*@
498:   VecISCopy - Copies between a reduced vector and the appropriate elements of a full-space vector.
499:               vfull[is[i]] = vreduced[i]

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

507:   Output Parameters:
508: . vfull    - the sum of the full-space vector and reduced-space vector

510:   Level: advanced

512: .seealso:  VecAXPY()
513: @*/
514: PetscErrorCode VecISCopy(Vec vfull, IS is, ScatterMode mode, Vec vreduced)
515: {
516:   PetscInt       nfull, nreduced;

523:   VecGetSize(vfull, &nfull);
524:   VecGetSize(vreduced, &nreduced);

526:   if (nfull == nreduced) { /* Also takes care of masked vectors */
527:     VecCopy(vreduced, vfull);
528:   } else {
529:     const PetscInt *id;
530:     PetscInt        i, n, m, rstart;

532:     ISGetIndices(is, &id);
533:     ISGetLocalSize(is, &n);
534:     VecGetLocalSize(vreduced, &m);
535:     VecGetOwnershipRange(vfull, &rstart, NULL);
536:     if (m != n) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "IS local length %D not equal to Vec local length %D", n, m);
537:     if (mode == SCATTER_FORWARD) {
538:       PetscScalar       *y;
539:       const PetscScalar *x;

541:       VecGetArray(vfull, &y);
542:       VecGetArrayRead(vreduced, &x);
543:       y   -= rstart;
544:       for (i = 0; i < n; ++i) {
545:         y[id[i]] = x[i];
546:       }
547:       y   += rstart;
548:       VecRestoreArrayRead(vreduced, &x);
549:       VecRestoreArray(vfull, &y);
550:     } else if (mode == SCATTER_REVERSE) {
551:       PetscScalar       *x;
552:       const PetscScalar *y;

554:       VecGetArrayRead(vfull, &y);
555:       VecGetArray(vreduced, &x);
556:       for (i = 0; i < n; ++i) {
557:         x[i] = y[id[i]-rstart];
558:       }
559:       VecRestoreArray(vreduced, &x);
560:       VecRestoreArrayRead(vfull, &y);
561:     } else SETERRQ(PetscObjectComm((PetscObject) vfull), PETSC_ERR_ARG_WRONG, "Only forward or reverse modes are legal");
562:     ISRestoreIndices(is, &id);
563:   }
564:   return(0);
565: }

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

570:    Collective on IS

572:    Input Parameter:
573: +  S -  a PETSc IS
574: -  V - the reference vector space

576:    Output Parameter:
577: .  T -  the complement of S

579: .seealso ISCreateGeneral()

581:    Level: advanced
582: @*/
583: PetscErrorCode ISComplementVec(IS S, Vec V, IS *T)
584: {
586:   PetscInt       start, end;

589:   VecGetOwnershipRange(V,&start,&end);
590:   ISComplement(S,start,end,T);
591:   return(0);
592: }

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

597:    Input Parameter:
598: +  V - the vector
599: .  S -  the locations in the vector
600: -  c - the constant

602: .seealso VecSet()

604:    Level: advanced
605: @*/
606: PetscErrorCode VecISSet(Vec V,IS S, PetscScalar c)
607: {
609:   PetscInt       nloc,low,high,i;
610:   const PetscInt *s;
611:   PetscScalar    *v;


619:   VecGetOwnershipRange(V,&low,&high);
620:   ISGetLocalSize(S,&nloc);
621:   ISGetIndices(S,&s);
622:   VecGetArray(V,&v);
623:   for (i=0; i<nloc; ++i){
624:     v[s[i]-low] = c;
625:   }
626:   ISRestoreIndices(S,&s);
627:   VecRestoreArray(V,&v);
628:   return(0);
629: }

631: #if !defined(PETSC_USE_COMPLEX)
632: /*@C
633:   VecBoundGradientProjection - Projects  vector according to this definition.
634:   If XL[i] < X[i] < XU[i], then GP[i] = G[i];
635:   If X[i] <= XL[i], then GP[i] = min(G[i],0);
636:   If X[i] >= XU[i], then GP[i] = max(G[i],0);

638:   Input Parameters:
639: + G - current gradient vector
640: . X - current solution vector with XL[i] <= X[i] <= XU[i]
641: . XL - lower bounds
642: - XU - upper bounds

644:   Output Parameter:
645: . GP - gradient projection vector

647:   Notes: GP may be the same vector as G

649:   Level: advanced
650: @*/
651: PetscErrorCode VecBoundGradientProjection(Vec G, Vec X, Vec XL, Vec XU, Vec GP)
652: {

654:   PetscErrorCode  ierr;
655:   PetscInt        n,i;
656:   const PetscReal *xptr,*xlptr,*xuptr;
657:   PetscReal       *gptr,*gpptr;
658:   PetscReal       xval,gpval;

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

668:   VecGetLocalSize(X,&n);

670:   VecGetArrayRead(X,&xptr);
671:   VecGetArrayRead(XL,&xlptr);
672:   VecGetArrayRead(XU,&xuptr);
673:   VecGetArrayPair(G,GP,&gptr,&gpptr);

675:   for (i=0; i<n; ++i){
676:     gpval = gptr[i]; xval = xptr[i];

678:     if (gpval>0.0 && xval<=xlptr[i]){
679:       gpval = 0.0;
680:     } else if (gpval<0.0 && xval>=xuptr[i]){
681:       gpval = 0.0;
682:     }
683:     gpptr[i] = gpval;
684:   }

686:   VecRestoreArrayRead(X,&xptr);
687:   VecRestoreArrayRead(XL,&xlptr);
688:   VecRestoreArrayRead(XU,&xuptr);
689:   VecRestoreArrayPair(G,GP,&gptr,&gpptr);
690:   return(0);
691: }
692: #endif

694: /*@
695:      VecStepMaxBounded - See below

697:      Collective on Vec

699:      Input Parameters:
700: +      X  - vector with no negative entries
701: .      XL - lower bounds
702: .      XU - upper bounds
703: -      DX  - step direction, can have negative, positive or zero entries

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

708:   Level: intermediate

710: @*/
711: PetscErrorCode VecStepMaxBounded(Vec X, Vec DX, Vec XL, Vec XU, PetscReal *stepmax)
712: {
713:   PetscErrorCode    ierr;
714:   PetscInt          i,nn;
715:   const PetscScalar *xx,*dx,*xl,*xu;
716:   PetscReal         localmax=0;


724:   VecGetArrayRead(X,&xx);
725:   VecGetArrayRead(XL,&xl);
726:   VecGetArrayRead(XU,&xu);
727:   VecGetArrayRead(DX,&dx);
728:   VecGetLocalSize(X,&nn);
729:   for (i=0;i<nn;i++){
730:     if (PetscRealPart(dx[i]) > 0){
731:       localmax=PetscMax(localmax,PetscRealPart((xu[i]-xx[i])/dx[i]));
732:     } else if (PetscRealPart(dx[i])<0){
733:       localmax=PetscMax(localmax,PetscRealPart((xl[i]-xx[i])/dx[i]));
734:     }
735:   }
736:   VecRestoreArrayRead(X,&xx);
737:   VecRestoreArrayRead(XL,&xl);
738:   VecRestoreArrayRead(XU,&xu);
739:   VecRestoreArrayRead(DX,&dx);
740:   MPIU_Allreduce(&localmax,stepmax,1,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)X));
741:   return(0);
742: }

744: /*@
745:      VecStepBoundInfo - See below

747:      Collective on Vec

749:      Input Parameters:
750: +      X  - vector with no negative entries
751: .      XL - lower bounds
752: .      XU - upper bounds
753: -      DX  - step direction, can have negative, positive or zero entries

755:      Output Parameter:
756: +     boundmin -  (may be NULL this it is not computed) maximum value so that   XL[i] <= X[i] + boundmax*DX[i] <= XU[i]
757: .     wolfemin -  (may be NULL this it is not computed)
758: -     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]

760:      Notes: For complex numbers only compares the real part

762:   Level: advanced
763: @*/
764: PetscErrorCode VecStepBoundInfo(Vec X, Vec DX, Vec XL, Vec XU, PetscReal *boundmin, PetscReal *wolfemin, PetscReal *boundmax)
765: {
766:   PetscErrorCode    ierr;
767:   PetscInt          n,i;
768:   const PetscScalar *x,*xl,*xu,*dx;
769:   PetscReal         t;
770:   PetscReal         localmin=PETSC_INFINITY,localwolfemin=PETSC_INFINITY,localmax=-1;
771:   MPI_Comm          comm;


779:   VecGetArrayRead(X,&x);
780:   VecGetArrayRead(XL,&xl);
781:   VecGetArrayRead(XU,&xu);
782:   VecGetArrayRead(DX,&dx);
783:   VecGetLocalSize(X,&n);
784:   for (i=0; i<n; ++i){
785:     if (PetscRealPart(dx[i])>0 && PetscRealPart(xu[i]) < PETSC_INFINITY) {
786:       t=PetscRealPart((xu[i]-x[i])/dx[i]);
787:       localmin=PetscMin(t,localmin);
788:       if (localmin>0){
789:         localwolfemin = PetscMin(t,localwolfemin);
790:       }
791:       localmax = PetscMax(t,localmax);
792:     } else if (PetscRealPart(dx[i])<0 && PetscRealPart(xl[i]) > PETSC_NINFINITY) {
793:       t=PetscRealPart((xl[i]-x[i])/dx[i]);
794:       localmin = PetscMin(t,localmin);
795:       if (localmin>0){
796:         localwolfemin = PetscMin(t,localwolfemin);
797:       }
798:       localmax = PetscMax(t,localmax);
799:     }
800:   }

802:   VecRestoreArrayRead(X,&x);
803:   VecRestoreArrayRead(XL,&xl);
804:   VecRestoreArrayRead(XU,&xu);
805:   VecRestoreArrayRead(DX,&dx);
806:   ierr=PetscObjectGetComm((PetscObject)X,&comm);

808:   if (boundmin){
809:     MPIU_Allreduce(&localmin,boundmin,1,MPIU_REAL,MPIU_MIN,comm);
810:     PetscInfo1(X,"Step Bound Info: Closest Bound: %20.19e\n",(double)*boundmin);
811:   }
812:   if (wolfemin){
813:     MPIU_Allreduce(&localwolfemin,wolfemin,1,MPIU_REAL,MPIU_MIN,comm);
814:     PetscInfo1(X,"Step Bound Info: Wolfe: %20.19e\n",(double)*wolfemin);
815:   }
816:   if (boundmax) {
817:     MPIU_Allreduce(&localmax,boundmax,1,MPIU_REAL,MPIU_MAX,comm);
818:     if (*boundmax < 0) *boundmax=PETSC_INFINITY;
819:     PetscInfo1(X,"Step Bound Info: Max: %20.19e\n",(double)*boundmax);
820:   }
821:   return(0);
822: }

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

827:      Collective on Vec

829:      Input Parameters:
830: +      X  - vector with no negative entries
831: -      DX  - a step direction, can have negative, positive or zero entries

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

836:      Notes: For complex numbers only compares the real part

838:   Level: advanced
839:  @*/
840: PetscErrorCode VecStepMax(Vec X, Vec DX, PetscReal *step)
841: {
842:   PetscErrorCode    ierr;
843:   PetscInt          i,nn;
844:   PetscReal         stepmax=PETSC_INFINITY;
845:   const PetscScalar *xx,*dx;


851:   VecGetLocalSize(X,&nn);
852:   VecGetArrayRead(X,&xx);
853:   VecGetArrayRead(DX,&dx);
854:   for (i=0;i<nn;++i){
855:     if (PetscRealPart(xx[i]) < 0) SETERRQ(PETSC_COMM_SELF,1,"Vector must be positive");
856:     else if (PetscRealPart(dx[i])<0) stepmax=PetscMin(stepmax,PetscRealPart(-xx[i]/dx[i]));
857:   }
858:   VecRestoreArrayRead(X,&xx);
859:   VecRestoreArrayRead(DX,&dx);
860:   MPIU_Allreduce(&stepmax,step,1,MPIU_REAL,MPIU_MIN,PetscObjectComm((PetscObject)X));
861:   return(0);
862: }

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

867:   Logically Collective on v

869:   Input Parameter:
870: + v - the vector
871: - p - the exponent to use on each element

873:   Output Parameter:
874: . v - the vector

876:   Level: intermediate

878: @*/
879: PetscErrorCode VecPow(Vec v, PetscScalar p)
880: {
882:   PetscInt       n,i;
883:   PetscScalar    *v1;


888:   VecGetArray(v,&v1);
889:   VecGetLocalSize(v,&n);

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

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

946:   Logically Collective

948:   Input Parameters:
949: . Vec1, Vec2, Vec3 - The three vectors

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

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


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

978:   /* Assert that Vec1 != Vec2 and Vec2 != Vec3 */
989:   VecCheckSameSize(Vec1,1,Vec2,2);
990:   VecCheckSameSize(Vec1,1,Vec3,3);
991:   VecCheckSameSize(Vec1,1,VMedian,4);

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

1013:     for (i=0;i<n;++i){
1014:       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])));
1015:     }

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