Actual source code: projection.c

petsc-master 2018-06-18
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:
 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: }


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

263:   Collective on S

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

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

273:   Level: advanced
274: @*/

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

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

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

310:     PetscMalloc1(n,&vm);

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

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

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

335:   Collective on S

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

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

347:   Level: advanced
348: @*/

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

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

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

396:     PetscMalloc1(n,&vm);

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

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


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

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

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

449:   Level: advanced

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

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

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

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

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

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

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

511:   Level: advanced

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

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

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

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

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

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

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

571:    Collective on IS

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

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

580: .seealso ISCreateGeneral()

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

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

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

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

603: .seealso VecSet()

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

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

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

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

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

646:   Output Parameter:
647: . GP - gradient projection vector

649:   Notes:
650:     GP may be the same vector as G

652:   Level: advanced
653: @*/
654: PetscErrorCode VecBoundGradientProjection(Vec G, Vec X, Vec XL, Vec XU, Vec GP)
655: {

657:   PetscErrorCode  ierr;
658:   PetscInt        n,i;
659:   const PetscReal *xptr,*xlptr,*xuptr;
660:   PetscReal       *gptr,*gpptr;
661:   PetscReal       xval,gpval;

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

671:   VecGetLocalSize(X,&n);

673:   VecGetArrayRead(X,&xptr);
674:   VecGetArrayRead(XL,&xlptr);
675:   VecGetArrayRead(XU,&xuptr);
676:   VecGetArrayPair(G,GP,&gptr,&gpptr);

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

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

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

697: /*@
698:      VecStepMaxBounded - See below

700:      Collective on Vec

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

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

711:   Level: intermediate

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


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

747: /*@
748:      VecStepBoundInfo - See below

750:      Collective on Vec

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

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

763:      Notes:
764:     For complex numbers only compares the real part

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


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

806:   VecRestoreArrayRead(X,&x);
807:   VecRestoreArrayRead(XL,&xl);
808:   VecRestoreArrayRead(XU,&xu);
809:   VecRestoreArrayRead(DX,&dx);
810:   ierr=PetscObjectGetComm((PetscObject)X,&comm);

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

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

831:      Collective on Vec

833:      Input Parameters:
834: +      X  - vector with no negative entries
835: -      DX  - a step direction, can have negative, positive or zero entries

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

840:      Notes:
841:     For complex numbers only compares the real part

843:   Level: advanced
844:  @*/
845: PetscErrorCode VecStepMax(Vec X, Vec DX, PetscReal *step)
846: {
847:   PetscErrorCode    ierr;
848:   PetscInt          i,nn;
849:   PetscReal         stepmax=PETSC_INFINITY;
850:   const PetscScalar *xx,*dx;


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

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

872:   Logically Collective on v

874:   Input Parameter:
875: + v - the vector
876: - p - the exponent to use on each element

878:   Output Parameter:
879: . v - the vector

881:   Level: intermediate

883: @*/
884: PetscErrorCode VecPow(Vec v, PetscScalar p)
885: {
887:   PetscInt       n,i;
888:   PetscScalar    *v1;


893:   VecGetArray(v,&v1);
894:   VecGetLocalSize(v,&n);

896:   if (1.0 == p) {
897:   } else if (-1.0 == p) {
898:     for (i = 0; i < n; ++i){
899:       v1[i] = 1.0 / v1[i];
900:     }
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]) {
906:         v1[i] = 1.0;
907:       }
908:     }
909:   } else if (0.5 == p) {
910:     for (i = 0; i < n; ++i) {
911:       if (PetscRealPart(v1[i]) >= 0) {
912:         v1[i] = PetscSqrtScalar(v1[i]);
913:       } else {
914:         v1[i] = PETSC_INFINITY;
915:       }
916:     }
917:   } else if (-0.5 == p) {
918:     for (i = 0; i < n; ++i) {
919:       if (PetscRealPart(v1[i]) >= 0) {
920:         v1[i] = 1.0 / PetscSqrtScalar(v1[i]);
921:       } else {
922:         v1[i] = PETSC_INFINITY;
923:       }
924:     }
925:   } else if (2.0 == p) {
926:     for (i = 0; i < n; ++i){
927:       v1[i] *= v1[i];
928:     }
929:   } else if (-2.0 == p) {
930:     for (i = 0; i < n; ++i){
931:       v1[i] = 1.0 / (v1[i] * v1[i]);
932:     }
933:   } else {
934:     for (i = 0; i < n; ++i) {
935:       if (PetscRealPart(v1[i]) >= 0) {
936:         v1[i] = PetscPowScalar(v1[i],p);
937:       } else {
938:         v1[i] = PETSC_INFINITY;
939:       }
940:     }
941:   }
942:   VecRestoreArray(v,&v1);
943:   return(0);
944: }

946: /*@
947:   VecMedian - Computes the componentwise median of three vectors
948:   and stores the result in this vector.  Used primarily for projecting
949:   a vector within upper and lower bounds.

951:   Logically Collective

953:   Input Parameters:
954: . Vec1, Vec2, Vec3 - The three vectors

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

959:   Level: advanced
960: @*/
961: PetscErrorCode VecMedian(Vec Vec1, Vec Vec2, Vec Vec3, Vec VMedian)
962: {
963:   PetscErrorCode    ierr;
964:   PetscInt          i,n,low1,high1;
965:   const PetscScalar *v1,*v2,*v3;
966:   PetscScalar       *vmed;


974:   if (Vec1==Vec2 || Vec1==Vec3){
975:     VecCopy(Vec1,VMedian);
976:     return(0);
977:   }
978:   if (Vec2==Vec3){
979:     VecCopy(Vec2,VMedian);
980:     return(0);
981:   }

983:   /* Assert that Vec1 != Vec2 and Vec2 != Vec3 */
994:   VecCheckSameSize(Vec1,1,Vec2,2);
995:   VecCheckSameSize(Vec1,1,Vec3,3);
996:   VecCheckSameSize(Vec1,1,VMedian,4);

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

1018:     for (i=0;i<n;++i){
1019:       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])));
1020:     }

1022:     VecRestoreArray(VMedian,&vmed);
1023:     if (VMedian != Vec1){
1024:       VecRestoreArrayRead(Vec1,&v1);
1025:     }
1026:     if (VMedian != Vec2){
1027:       VecRestoreArrayRead(Vec2,&v2);
1028:     }
1029:     if (VMedian != Vec3){
1030:       VecRestoreArrayRead(Vec3,&v3);
1031:     }
1032:   }
1033:   return(0);
1034: }