Actual source code: projection.c

petsc-3.8.3 2017-12-09
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;
 84:   PetscInt          n,low,high,n_lt=0;
 85:   PetscInt          *lt = NULL;
 86:   const PetscScalar *v1,*v2;

 92:   VecCheckSameSize(Vec1,1,Vec2,2);
 93: 
 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:     }
104:     PetscMalloc1(n,&lt );

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

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

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

125:   Collective on S

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

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

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

136:   For complex numbers this only compares the real part

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

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

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

164:     PetscMalloc1(n, &gt );

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

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

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

185:   Collective on S

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

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

195:   Notes:
196:   The vectors must have the same parallel layout

198:   For complex numbers this only compares the real part

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

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

213:   VecCheckSameSize(V,2,VecLow,1);
214:   VecCheckSameSize(V,2,VecHigh,3);

216:   VecGetOwnershipRange(VecLow, &low, &high);
217:   VecGetLocalSize(VecLow,&n);
218:   if (n>0){
219:     VecGetArrayRead(VecLow,&v1);
220:     if (VecLow != VecHigh){
221:       VecGetArrayRead(VecHigh,&v2);
222:     } else {
223:       v2=v1;
224:     }
225:     if (V != VecLow && V != VecHigh){
226:       VecGetArrayRead(V,&vmiddle);
227:     } else if ( V==VecLow ){
228:       vmiddle=v1;
229:     } else {
230:       vmiddle =v2;
231:     }
232:     PetscMalloc1(n, &vm );
233:     for (i=0; i<n; i++){
234:       if (PetscRealPart(v1[i]) < PetscRealPart(vmiddle[i]) && PetscRealPart(vmiddle[i]) < PetscRealPart(v2[i])) {vm[n_vm]=low+i; n_vm++;}
235:     }
236:     VecRestoreArrayRead(VecLow,&v1);
237:     if (VecLow != VecHigh){
238:       VecRestoreArrayRead(VecHigh,&v2);
239:     }
240:     if (V != VecLow && V != VecHigh){
241:       VecRestoreArrayRead(V,&vmiddle);
242:     }
243:   }
244:   ISCreateGeneral(PetscObjectComm((PetscObject)V),n_vm,vm,PETSC_OWN_POINTER,S);
245:   return(0);
246: }


249: /*@
250:   VecWhichBetweenOrEqual - Creates an index set containing the indices
251:   where  VecLow <= V <= VecHigh

253:   Collective on S

255:   Input Parameters:
256: + VecLow - lower bound
257: . V - Vector to compare
258: - VecHigh - higher bound

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

263:   Level: advanced
264: @*/

266: PetscErrorCode VecWhichBetweenOrEqual(Vec VecLow, Vec V, Vec VecHigh, IS * S)
267: {
269:   PetscInt       n,low,high,n_vm=0,i;
270:   PetscInt       *vm = NULL;
271:   PetscScalar    *v1,*v2,*vmiddle;

276:   VecCheckSameSize(V,2,VecLow,1);
277:   VecCheckSameSize(V,2,VecHigh,3);

279:   VecGetOwnershipRange(VecLow, &low, &high);
280:   VecGetLocalSize(VecLow,&n);

282:   if (n>0){
283:     VecGetArray(VecLow,&v1);
284:     if (VecLow != VecHigh){
285:       VecGetArray(VecHigh,&v2);
286:     } else {
287:       v2=v1;
288:     }
289:     if ( V != VecLow && V != VecHigh){
290:       VecGetArray(V,&vmiddle);
291:     } else if ( V==VecLow ){
292:       vmiddle=v1;
293:     } else {
294:       vmiddle =v2;
295:     }

297:     PetscMalloc1(n, &vm );

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

303:     VecRestoreArray(VecLow,&v1);
304:     if (VecLow != VecHigh){
305:       VecRestoreArray(VecHigh,&v2);
306:     }
307:     if ( V != VecLow && V != VecHigh){
308:       VecRestoreArray(V,&vmiddle);
309:     }
310:   }
311:   ISCreateGeneral(PetscObjectComm((PetscObject)V),n_vm,vm,PETSC_OWN_POINTER,S);
312:   return(0);
313: }

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

319:   Input Parameters:
320: + vfull - the full-space vector
321: . vreduced - the reduced-space vector
322: - is - the index set for the reduced space

324:   Output Parameters:
325: . vfull - the sum of the full-space vector and reduced-space vector

327:   Level: advanced

329: .seealso:  VecAXPY()
330: @*/
331: PetscErrorCode VecISAXPY(Vec vfull, IS is, PetscScalar alpha,Vec vreduced)
332: {
333:   PetscInt       nfull,nreduced;

340:   VecGetSize(vfull,&nfull);
341:   VecGetSize(vreduced,&nreduced);

343:   if (nfull == nreduced) { /* Also takes care of masked vectors */
344:     VecAXPY(vfull,alpha,vreduced);
345:   } else {
346:     PetscScalar      *y;
347:     const PetscScalar *x;
348:     PetscInt          i,n,m,rstart;
349:     const PetscInt    *id;

351:     VecGetArray(vfull,&y);
352:     VecGetArrayRead(vreduced,&x);
353:     ISGetIndices(is,&id);
354:     ISGetLocalSize(is,&n);
355:     VecGetLocalSize(vreduced,&m);
356:     if (m != n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"IS local length not equal to Vec local length");
357:     VecGetOwnershipRange(vfull,&rstart,NULL);
358:     y -= rstart;
359:     if (alpha == 1.0) {
360:       for (i=0; i<n; i++) {
361:         y[id[i]] += x[i];
362:       }
363:     } else {
364:       for (i=0; i<n; i++) {
365:         y[id[i]] += alpha*x[i];
366:       }
367:     }
368:     y += rstart;
369:     ISRestoreIndices(is,&id);
370:     VecRestoreArray(vfull,&y);
371:     VecRestoreArrayRead(vreduced,&x);
372:   }
373:   return(0);
374: }

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

379:    Collective on IS

381:    Input Parameter:
382: +  S -  a PETSc IS
383: -  V - the reference vector space

385:    Output Parameter:
386: .  T -  the complement of S

388: .seealso ISCreateGeneral()

390:    Level: advanced
391: @*/
392: PetscErrorCode ISComplementVec(IS S, Vec V, IS *T)
393: {
395:   PetscInt       start, end;

398:   VecGetOwnershipRange(V,&start,&end);
399:   ISComplement(S,start,end,T);
400:   return(0);
401: }

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

406:    Input Parameter:
407: +  V - the vector
408: .  S -  the locations in the vector
409: -  c - the constant

411: .seealso VecSet()

413:    Level: advanced
414: @*/
415: PetscErrorCode VecISSet(Vec V,IS S, PetscScalar c)
416: {
418:   PetscInt       nloc,low,high,i;
419:   const PetscInt *s;
420:   PetscScalar    *v;


428:   VecGetOwnershipRange(V, &low, &high);
429:   ISGetLocalSize(S,&nloc);
430:   ISGetIndices(S, &s);
431:   VecGetArray(V,&v);
432:   for (i=0; i<nloc; i++){
433:     v[s[i]-low] = c;
434:   }
435:   ISRestoreIndices(S, &s);
436:   VecRestoreArray(V,&v);
437:   return(0);
438: }

440: #if !defined(PETSC_USE_COMPLEX)
441: /*@C
442:   VecBoundGradientProjection - Projects  vector according to this definition.
443:   If XL[i] < X[i] < XU[i], then GP[i] = G[i];
444:   If X[i]<=XL[i], then GP[i] = min(G[i],0);
445:   If X[i]>=XU[i], then GP[i] = max(G[i],0);

447:   Input Parameters:
448: + G - current gradient vector
449: . X - current solution vector
450: . XL - lower bounds
451: - XU - upper bounds

453:   Output Parameter:
454: . GP - gradient projection vector

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

458:   Level: advanced
459: C@*/
460: PetscErrorCode VecBoundGradientProjection(Vec G, Vec X, Vec XL, Vec XU, Vec GP)
461: {

463:   PetscErrorCode  ierr;
464:   PetscInt        n,i;
465:   const PetscReal *xptr,*xlptr,*xuptr;
466:   PetscReal       *gptr,*gpptr;
467:   PetscReal       xval,gpval;

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

477:   VecGetLocalSize(X,&n);

479:   VecGetArrayRead(X,&xptr);
480:   VecGetArrayRead(XL,&xlptr);
481:   VecGetArrayRead(XU,&xuptr);
482:   VecGetArrayPair(G,GP,&gptr,&gpptr);

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

487:     if (gpval>0 && xval<=xlptr[i]){
488:       gpval = 0;
489:     } else if (gpval<0 && xval>=xuptr[i]){
490:       gpval = 0;
491:     }
492:     gpptr[i] = gpval;
493:   }

495:   VecRestoreArrayRead(X,&xptr);
496:   VecRestoreArrayRead(XL,&xlptr);
497:   VecRestoreArrayRead(XU,&xuptr);
498:   VecRestoreArray(G,&gptr);
499:   VecRestoreArrayPair(G,GP,&gptr,&gpptr);
500:   return(0);
501: }
502: #endif

504: /*@
505:      VecStepMaxBounded - See below

507:      Collective on Vec

509:      Input Parameters:
510: +      X  - vector with no negative entries
511: .      XL - lower bounds
512: .      XU - upper bounds
513: -      DX  - step direction, can have negative, positive or zero entries

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

518: @*/
519: PetscErrorCode VecStepMaxBounded(Vec X, Vec DX, Vec XL, Vec XU, PetscReal *stepmax)
520: {
521:   PetscErrorCode    ierr;
522:   PetscInt          i,nn;
523:   const PetscScalar *xx,*dx,*xl,*xu;
524:   PetscReal         localmax=0;


532:   VecGetArrayRead(X,&xx);
533:   VecGetArrayRead(XL,&xl);
534:   VecGetArrayRead(XU,&xu);
535:   VecGetArrayRead(DX,&dx);
536:   VecGetLocalSize(X,&nn);
537:   for (i=0;i<nn;i++){
538:     if (PetscRealPart(dx[i]) > 0){
539:       localmax=PetscMax(localmax,PetscRealPart((xu[i]-xx[i])/dx[i]));
540:     } else if (PetscRealPart(dx[i])<0){
541:       localmax=PetscMax(localmax,PetscRealPart((xl[i]-xx[i])/dx[i]));
542:     }
543:   }
544:   VecRestoreArrayRead(X,&xx);
545:   VecRestoreArrayRead(XL,&xl);
546:   VecRestoreArrayRead(XU,&xu);
547:   VecRestoreArrayRead(DX,&dx);
548:   MPIU_Allreduce(&localmax,stepmax,1,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)X));
549:   return(0);
550: }

552: /*@
553:      VecStepBoundInfo - See below

555:      Collective on Vec

557:      Input Parameters:
558: +      X  - vector with no negative entries
559: .      XL - lower bounds
560: .      XU - upper bounds
561: -      DX  - step direction, can have negative, positive or zero entries

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

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

570:   Level: advanced
571: @*/
572: PetscErrorCode VecStepBoundInfo(Vec X, Vec DX, Vec XL, Vec XU, PetscReal *boundmin, PetscReal *wolfemin, PetscReal *boundmax)
573: {
574:   PetscErrorCode    ierr;
575:   PetscInt          n,i;
576:   const PetscScalar *x,*xl,*xu,*dx;
577:   PetscReal         t;
578:   PetscReal         localmin=PETSC_INFINITY,localwolfemin=PETSC_INFINITY,localmax=-1;
579:   MPI_Comm          comm;


587:   VecGetArrayRead(X,&x);
588:   VecGetArrayRead(XL,&xl);
589:   VecGetArrayRead(XU,&xu);
590:   VecGetArrayRead(DX,&dx);
591:   VecGetLocalSize(X,&n);
592:   for (i=0;i<n;i++){
593:     if (PetscRealPart(dx[i])>0 && PetscRealPart(xu[i]) < PETSC_INFINITY) {
594:       t=PetscRealPart((xu[i]-x[i])/dx[i]);
595:       localmin=PetscMin(t,localmin);
596:       if (localmin>0){
597:         localwolfemin = PetscMin(t,localwolfemin);
598:       }
599:       localmax = PetscMax(t,localmax);
600:     } else if (PetscRealPart(dx[i])<0 && PetscRealPart(xl[i]) > PETSC_NINFINITY) {
601:       t=PetscRealPart((xl[i]-x[i])/dx[i]);
602:       localmin = PetscMin(t,localmin);
603:       if (localmin>0){
604:         localwolfemin = PetscMin(t,localwolfemin);
605:       }
606:       localmax = PetscMax(t,localmax);
607:     }
608:   }

610:   VecRestoreArrayRead(X,&x);
611:   VecRestoreArrayRead(XL,&xl);
612:   VecRestoreArrayRead(XU,&xu);
613:   VecRestoreArrayRead(DX,&dx);
614:   ierr=PetscObjectGetComm((PetscObject)X,&comm);

616:   if (boundmin){
617:     MPIU_Allreduce(&localmin,boundmin,1,MPIU_REAL,MPIU_MIN,comm);
618:     PetscInfo1(X,"Step Bound Info: Closest Bound: %g \n",(double)*boundmin);
619:   }
620:   if (wolfemin){
621:     MPIU_Allreduce(&localwolfemin,wolfemin,1,MPIU_REAL,MPIU_MIN,comm);
622:     PetscInfo1(X,"Step Bound Info: Wolfe: %g \n",(double)*wolfemin);
623:   }
624:   if (boundmax) {
625:     MPIU_Allreduce(&localmax,boundmax,1,MPIU_REAL,MPIU_MAX,comm);
626:     if (*boundmax < 0) *boundmax=PETSC_INFINITY;
627:     PetscInfo1(X,"Step Bound Info: Max: %g \n",(double)*boundmax);
628:   }
629:   return(0);
630: }

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

635:      Collective on Vec

637:      Input Parameters:
638: +      X  - vector with no negative entries
639: -      DX  - a step direction, can have negative, positive or zero entries

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

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

646:   Level: advanced
647:  @*/
648: PetscErrorCode VecStepMax(Vec X, Vec DX, PetscReal *step)
649: {
650:   PetscErrorCode    ierr;
651:   PetscInt          i, nn;
652:   PetscReal         stepmax=PETSC_INFINITY;
653:   const PetscScalar *xx, *dx;


659:   VecGetLocalSize(X,&nn);
660:   VecGetArrayRead(X,&xx);
661:   VecGetArrayRead(DX,&dx);
662:   for (i=0;i<nn;i++){
663:     if (PetscRealPart(xx[i]) < 0) SETERRQ(PETSC_COMM_SELF,1,"Vector must be positive");
664:     else if (PetscRealPart(dx[i])<0) stepmax=PetscMin(stepmax,PetscRealPart(-xx[i]/dx[i]));
665:   }
666:   VecRestoreArrayRead(X,&xx);
667:   VecRestoreArrayRead(DX,&dx);
668:   MPIU_Allreduce(&stepmax,step,1,MPIU_REAL,MPIU_MIN,PetscObjectComm((PetscObject)X));
669:   return(0);
670: }

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

675:   Logically Collective on v

677:   Input Parameter:
678: + v - the vector
679: - p - the exponent to use on each element

681:   Output Parameter:
682: . v - the vector

684:   Level: intermediate

686: @*/
687: PetscErrorCode VecPow(Vec v, PetscScalar p)
688: {
690:   PetscInt       n,i;
691:   PetscScalar    *v1;


696:   VecGetArray(v, &v1);
697:   VecGetLocalSize(v, &n);

699:   if (1.0 == p) {
700:   } else if (-1.0 == p) {
701:     for (i = 0; i < n; ++i){
702:       v1[i] = 1.0 / v1[i];
703:     }
704:   } else if (0.0 == p) {
705:     for (i = 0; i < n; ++i){
706:       /*  Not-a-number left alone
707:           Infinity set to one  */
708:       if (v1[i] == v1[i]) {
709:         v1[i] = 1.0;
710:       }
711:     }
712:   } else if (0.5 == p) {
713:     for (i = 0; i < n; ++i) {
714:       if (PetscRealPart(v1[i]) >= 0) {
715:         v1[i] = PetscSqrtScalar(v1[i]);
716:       } else {
717:         v1[i] = PETSC_INFINITY;
718:       }
719:     }
720:   } else if (-0.5 == p) {
721:     for (i = 0; i < n; ++i) {
722:       if (PetscRealPart(v1[i]) >= 0) {
723:         v1[i] = 1.0 / PetscSqrtScalar(v1[i]);
724:       } else {
725:         v1[i] = PETSC_INFINITY;
726:       }
727:     }
728:   } else if (2.0 == p) {
729:     for (i = 0; i < n; ++i){
730:       v1[i] *= v1[i];
731:     }
732:   } else if (-2.0 == p) {
733:     for (i = 0; i < n; ++i){
734:       v1[i] = 1.0 / (v1[i] * v1[i]);
735:     }
736:   } else {
737:     for (i = 0; i < n; ++i) {
738:       if (PetscRealPart(v1[i]) >= 0) {
739:         v1[i] = PetscPowScalar(v1[i], p);
740:       } else {
741:         v1[i] = PETSC_INFINITY;
742:       }
743:     }
744:   }
745:   VecRestoreArray(v,&v1);
746:   return(0);
747: }

749: /*@
750:   VecMedian - Computes the componentwise median of three vectors
751:   and stores the result in this vector.  Used primarily for projecting
752:   a vector within upper and lower bounds.

754:   Logically Collective

756:   Input Parameters:
757: . Vec1, Vec2, Vec3 - The three vectors

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

762:   Developers Note: Should VMedian be allow to be one of the input vectors?

764:   Level: advanced
765: @*/
766: PetscErrorCode VecMedian(Vec Vec1, Vec Vec2, Vec Vec3, Vec VMedian)
767: {
769:   PetscInt       i,n,low1,high1;
770:   PetscScalar    *v1,*v2,*v3,*vmed;


778:   if (Vec1==Vec2 || Vec1==Vec3){
779:     VecCopy(Vec1,VMedian);
780:     return(0);
781:   }
782:   if (Vec2==Vec3){
783:     VecCopy(Vec2,VMedian);
784:     return(0);
785:   }

794:   VecCheckSameSize(Vec1,1,Vec2,2);
795:   VecCheckSameSize(Vec1,1,VMedian,4);

797:   VecGetOwnershipRange(Vec1, &low1, &high1);
798:   VecGetArray(Vec1,&v1);
799:   VecGetArray(Vec2,&v2);
800:   VecGetArray(Vec3,&v3);

802:   if ( VMedian != Vec1 && VMedian != Vec2 && VMedian != Vec3){
803:     VecGetArray(VMedian,&vmed);
804:   } else if ( VMedian==Vec1 ){
805:     vmed=v1;
806:   } else if ( VMedian==Vec2 ){
807:     vmed=v2;
808:   } else {
809:     vmed=v3;
810:   }

812:   VecGetLocalSize(Vec1,&n);
813:   for (i=0;i<n;i++){
814:     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])));
815:   }
816:   VecRestoreArray(Vec1,&v1);
817:   VecRestoreArray(Vec2,&v2);
818:   VecRestoreArray(Vec3,&v3);

820:   if (VMedian!=Vec1 && VMedian != Vec2 && VMedian != Vec3){
821:     VecRestoreArray(VMedian,&vmed);
822:   }
823:   return(0);
824: }