Actual source code: projection.c

petsc-3.5.2 2014-09-08
Report Typos and Errors
  1: #include <petsc-private/vecimpl.h>    /*I   "petscvec.h"  I*/

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

  9:   Collective on Vec

 11:   Input Parameters:
 12: . Vec1, Vec2 - the two vectors to compare

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

 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,low2,high2;
 24:   PetscInt        *same = NULL;
 25:   PetscScalar     *v1,*v2;
 26:   MPI_Comm        comm;


 33:   VecGetOwnershipRange(Vec1, &low, &high);
 34:   VecGetOwnershipRange(Vec2, &low2, &high2);
 35:   if ( low != low2 || high != high2 ) SETERRQ(PETSC_COMM_SELF,1,"Vectors must have identical layout");

 37:   VecGetLocalSize(Vec1,&n);
 38:   if (n>0){
 39:     if (Vec1 == Vec2){
 40:       VecGetArray(Vec1,&v1);
 41:       v2=v1;
 42:     } else {
 43:       VecGetArray(Vec1,&v1);
 44:       VecGetArray(Vec2,&v2);
 45:     }

 47:     PetscMalloc1( n,&same );

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

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

 67: /*@
 68:   VecWhichLessThan - Creates an index set containing the indices
 69:   where the vectors Vec1 < Vec2

 71:   Collective on S

 73:   Input Parameters:
 74: . Vec1, Vec2 - the two vectors to compare

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

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


 95:   VecGetOwnershipRange(Vec1, &low, &high);
 96:   VecGetOwnershipRange(Vec2, &low2, &high2);
 97:   if ( low != low2 || high != high2 ) SETERRQ(PETSC_COMM_SELF,1,"Vectors must haveidentical layout");

 99:   VecGetLocalSize(Vec1,&n);
100:   if (n>0){
101:     if (Vec1 == Vec2){
102:       VecGetArray(Vec1,&v1);
103:       v2=v1;
104:     } else {
105:       VecGetArray(Vec1,&v1);
106:       VecGetArray(Vec2,&v2);
107:     }
108:     PetscMalloc1(n,&lt );

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

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

128: /*@
129:   VecWhichGreaterThan - Creates an index set containing the indices
130:   where the vectors Vec1 > Vec2

132:   Collective on S

134:   Input Parameters:
135: . Vec1, Vec2 - the two vectors to compare

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

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


155:   VecGetOwnershipRange(Vec1, &low, &high);
156:   VecGetOwnershipRange(Vec2, &low2, &high2);
157:   if ( low != low2 || high != high2 ) SETERRQ(PETSC_COMM_SELF,1,"Vectors must be have identical layout");

159:   VecGetLocalSize(Vec1,&n);

161:   if (n>0){

163:     if (Vec1 == Vec2){
164:       VecGetArray(Vec1,&v1);
165:       v2=v1;
166:     } else {
167:       VecGetArray(Vec1,&v1);
168:       VecGetArray(Vec2,&v2);
169:     }

171:     PetscMalloc1(n, &gt );

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

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

191: /*@
192:   VecWhichBetween - Creates an index set containing the indices
193:                where  VecLow < V < VecHigh

195:   Collective on S

197:   Input Parameters:
198: + VecLow - lower bound
199: . V - Vector to compare
200: - VecHigh - higher bound

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

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

211:   PetscInt       n,low,high,low2,high2,low3,high3,n_vm=0;
212:   PetscInt       *vm,i;
213:   PetscScalar    *v1,*v2,*vmiddle;
214:   MPI_Comm       comm;


219:   VecGetOwnershipRange(VecLow, &low, &high);
220:   VecGetOwnershipRange(VecHigh, &low2, &high2);
221:   VecGetOwnershipRange(V, &low3, &high3);
222:   if ( low!=low2 || high!=high2 || low!=low3 || high!=high3) SETERRQ(PETSC_COMM_SELF,1,"Vectors must have identical layout");

224:   VecGetLocalSize(VecLow,&n);
225:   if (n>0){
226:     VecGetArray(VecLow,&v1);
227:     if (VecLow != VecHigh){
228:       VecGetArray(VecHigh,&v2);
229:     } else {
230:       v2=v1;
231:     }
232:     if ( V != VecLow && V != VecHigh){
233:       VecGetArray(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:     VecRestoreArray(VecLow,&v1);
247:     if (VecLow != VecHigh){
248:       VecRestoreArray(VecHigh,&v2);
249:     }
250:     if ( V != VecLow && V != VecHigh){
251:       VecRestoreArray(V,&vmiddle);
252:     }
253:   }
254:   PetscObjectGetComm((PetscObject)V,&comm);
255:   ISCreateGeneral(comm,n_vm,vm,PETSC_OWN_POINTER,S);
256:   return(0);
257: }


262: /*@
263:   VecWhichBetweenOrEqual - Creates an index set containing the indices
264:   where  VecLow <= V <= VecHigh

266:   Collective on S

268:   Input Parameters:
269: + VecLow - lower bound
270: . V - Vector to compare
271: - VecHigh - higher bound

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

276:   Level: advanced
277: @*/

279: PetscErrorCode VecWhichBetweenOrEqual(Vec VecLow, Vec V, Vec VecHigh, IS * S)
280: {
282:   PetscInt       n,low,high,low2,high2,low3,high3,n_vm=0,i;
283:   PetscInt       *vm = NULL;
284:   PetscScalar    *v1,*v2,*vmiddle;
285:   MPI_Comm       comm;


290:   VecGetOwnershipRange(VecLow, &low, &high);
291:   VecGetOwnershipRange(VecHigh, &low2, &high2);
292:   VecGetOwnershipRange(V, &low3, &high3);
293:   if ( low!=low2 || high!=high2 || low!=low3 || high!=high3 ) SETERRQ(PETSC_COMM_SELF,1,"Vectors must have identical layout");

295:   VecGetLocalSize(VecLow,&n);

297:   if (n>0){
298:     VecGetArray(VecLow,&v1);
299:     if (VecLow != VecHigh){
300:       VecGetArray(VecHigh,&v2);
301:     } else {
302:       v2=v1;
303:     }
304:     if ( V != VecLow && V != VecHigh){
305:       VecGetArray(V,&vmiddle);
306:     } else if ( V==VecLow ){
307:       vmiddle=v1;
308:     } else {
309:       vmiddle =v2;
310:     }

312:     PetscMalloc1(n, &vm );

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

318:     VecRestoreArray(VecLow,&v1);
319:     if (VecLow != VecHigh){
320:       VecRestoreArray(VecHigh,&v2);
321:     }
322:     if ( V != VecLow && V != VecHigh){
323:       VecRestoreArray(V,&vmiddle);
324:     }
325:   }
326:   PetscObjectGetComm((PetscObject)V,&comm);
327:   ISCreateGeneral(comm,n_vm,vm,PETSC_OWN_POINTER,S);
328:   return(0);
329: }

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

337:   Input Parameters:
338: + vfull - the full-space vector
339: . vreduced - the reduced-space vector
340: - is - the index set for the reduced space

342:   Output Parameters:
343: . vfull - the sum of the full-space vector and reduced-space vector

345:   Level: advanced

347: .seealso:  VecAXPY()
348: @*/
349: PetscErrorCode VecISAXPY(Vec vfull, IS is, PetscScalar alpha,Vec vreduced)
350: {
351:   PetscInt       nfull,nreduced;
352:   MPI_Comm       comm;

359:   VecGetSize(vfull,&nfull);
360:   VecGetSize(vreduced,&nreduced);

362:   if (nfull == nreduced) { /* Also takes care of masked vectors */
363:     VecAXPY(vfull,alpha,vreduced);
364:   } else {
365:     PetscScalar      *y;
366:     const PetscScalar *x;
367:     PetscInt          i,n,m,rstart;
368:     const PetscInt    *id;

370:     PetscObjectGetComm((PetscObject)vfull,&comm);
371:     VecGetArray(vfull,&y);
372:     VecGetArrayRead(vreduced,&x);
373:     ISGetIndices(is,&id);
374:     ISGetLocalSize(is,&n);
375:     VecGetLocalSize(vreduced,&m);
376:     if (m != n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"IS local length not equal to Vec local length");
377:     VecGetOwnershipRange(vfull,&rstart,NULL);
378:     y -= rstart;
379:     if (alpha == 1.0) {
380:       for (i=0; i<n; i++) {
381:         y[id[i]] += x[i];
382:       }
383:     } else {
384:       for (i=0; i<n; i++) {
385:         y[id[i]] += alpha*x[i];
386:       }
387:     }
388:     y += rstart;
389:     ISRestoreIndices(is,&id);
390:     VecRestoreArray(vfull,&y);
391:     VecRestoreArrayRead(vreduced,&x);
392:   }
393:   return(0);
394: }

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

401:    Collective on IS

403:    Input Parameter:
404: +  S -  a PETSc IS
405: -  V - the reference vector space

407:    Output Parameter:
408: .  T -  the complement of S

410: .seealso ISCreateGeneral()

412:    Level: advanced
413: @*/
414: PetscErrorCode ISComplementVec(IS S, Vec V, IS *T)
415: {
417:   PetscInt       start, end;

420:   VecGetOwnershipRange(V,&start,&end);
421:   ISComplement(S,start,end,T);
422:   return(0);
423: }

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

430:    Input Parameter:
431: +  V - the vector
432: .  S -  the locations in the vector
433: -  c - the constant

435: .seealso VecSet()

437:    Level: advanced
438: @*/
439: PetscErrorCode VecISSet(Vec V,IS S, PetscScalar c)
440: {
442:   PetscInt       nloc,low,high,i;
443:   const PetscInt *s;
444:   PetscScalar    *v;


452:   VecGetOwnershipRange(V, &low, &high);
453:   ISGetLocalSize(S,&nloc);
454:   ISGetIndices(S, &s);
455:   VecGetArray(V,&v);
456:   for (i=0; i<nloc; i++){
457:     v[s[i]-low] = c;
458:   }
459:   ISRestoreIndices(S, &s);
460:   VecRestoreArray(V,&v);
461:   return(0);
462: }

464: #if !defined(PETSC_USE_COMPLEX)
467: /*@C
468:   VecBoundGradientProjection - Projects  vector according to this definition.
469:   If XL[i] < X[i] < XU[i], then GP[i] = G[i];
470:   If X[i]<=XL[i], then GP[i] = min(G[i],0);
471:   If X[i]>=XU[i], then GP[i] = max(G[i],0);

473:   Input Parameters:
474: + G - current gradient vector
475: . X - current solution vector
476: . XL - lower bounds
477: - XU - upper bounds

479:   Output Parameter:
480: . GP - gradient projection vector

482:   Level: advanced
483: C@*/
484: PetscErrorCode VecBoundGradientProjection(Vec G, Vec X, Vec XL, Vec XU, Vec GP)
485: {

488:   PetscInt       n,i;
489:   PetscReal      *xptr,*xlptr,*xuptr,*gptr,*gpptr;
490:   PetscReal      xval,gpval;

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

500:   VecGetLocalSize(X,&n);

502:   ierr=VecGetArray(X,&xptr);
503:   ierr=VecGetArray(XL,&xlptr);
504:   ierr=VecGetArray(XU,&xuptr);
505:   ierr=VecGetArray(G,&gptr);
506:   if (G!=GP){
507:     ierr=VecGetArray(GP,&gpptr);
508:   } else { gpptr=gptr; }

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

513:     if (gpval>0 && xval<=xlptr[i]){
514:       gpval = 0;
515:     } else if (gpval<0 && xval>=xuptr[i]){
516:       gpval = 0;
517:     }
518:     gpptr[i] = gpval;
519:   }

521:   ierr=VecRestoreArray(X,&xptr);
522:   ierr=VecRestoreArray(XL,&xlptr);
523:   ierr=VecRestoreArray(XU,&xuptr);
524:   ierr=VecRestoreArray(G,&gptr);
525:   if (G!=GP){
526:     ierr=VecRestoreArray(GP,&gpptr);
527:   }
528:   return(0);
529: }
530: #endif

534: /*@
535:      VecStepMaxBounded - See below

537:      Collective on Vec

539:      Input Parameters:
540: +      X  - vector with no negative entries
541: .      XL - lower bounds
542: .      XU - upper bounds
543: -      DX  - step direction, can have negative, positive or zero entries

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

548: @*/
549: PetscErrorCode VecStepMaxBounded(Vec X, Vec DX, Vec XL, Vec XU, PetscReal *stepmax)
550: {
552:   PetscInt       i,nn;
553:   PetscScalar    *xx,*dx,*xl,*xu;
554:   PetscReal      localmax=0;
555:   MPI_Comm       comm;


563:   VecGetArray(X,&xx);
564:   VecGetArray(XL,&xl);
565:   VecGetArray(XU,&xu);
566:   VecGetArray(DX,&dx);
567:   VecGetLocalSize(X,&nn);
568:   for (i=0;i<nn;i++){
569:     if (PetscRealPart(dx[i]) > 0){
570:       localmax=PetscMax(localmax,PetscRealPart((xu[i]-xx[i])/dx[i]));
571:     } else if (PetscRealPart(dx[i])<0){
572:       localmax=PetscMax(localmax,PetscRealPart((xl[i]-xx[i])/dx[i]));
573:     }
574:   }
575:   VecRestoreArray(X,&xx);
576:   VecRestoreArray(XL,&xl);
577:   VecRestoreArray(XU,&xu);
578:   VecRestoreArray(DX,&dx);
579:   PetscObjectGetComm((PetscObject)X,&comm);
580:   MPI_Allreduce(&localmax,stepmax,1,MPIU_REAL,MPIU_MAX,comm);
581:   return(0);
582: }

586: /*@
587:      VecStepBoundInfo - See below

589:      Collective on Vec

591:      Input Parameters:
592: +      X  - vector with no negative entries
593: .      XL - lower bounds
594: .      XU - upper bounds
595: -      DX  - step direction, can have negative, positive or zero entries

597:      Output Parameter:
598: +     boundmin -  maximum value so that   XL[i] <= X[i] + boundmax*DX[i] <= XU[i]
599: .     wolfemin -
600: -     boundmax -   minimum value so that X[i] + boundmax*DX[i] <= XL[i]  or  XU[i] <= X[i] + boundmax*DX[i]

602:   Level: advanced
603: @*/
604: PetscErrorCode VecStepBoundInfo(Vec X, Vec DX, Vec XL, Vec XU, PetscReal *boundmin, PetscReal *wolfemin, PetscReal *boundmax)
605: {
607:   PetscInt       n,i;
608:   PetscScalar    *x,*xl,*xu,*dx;
609:   PetscReal      t;
610:   PetscReal      localmin=PETSC_INFINITY,localwolfemin=PETSC_INFINITY,localmax=0;
611:   MPI_Comm       comm;


619:   ierr=VecGetArray(X,&x);
620:   ierr=VecGetArray(XL,&xl);
621:   ierr=VecGetArray(XU,&xu);
622:   ierr=VecGetArray(DX,&dx);
623:   VecGetLocalSize(X,&n);
624:   for (i=0;i<n;i++){
625:     if (PetscRealPart(dx[i])>0){
626:       t=PetscRealPart((xu[i]-x[i])/dx[i]);
627:       localmin=PetscMin(t,localmin);
628:       if (localmin>0){
629:         localwolfemin = PetscMin(t,localwolfemin);
630:       }
631:       localmax = PetscMax(t,localmax);
632:     } else if (PetscRealPart(dx[i])<0){
633:       t=PetscRealPart((xl[i]-x[i])/dx[i]);
634:       localmin = PetscMin(t,localmin);
635:       if (localmin>0){
636:         localwolfemin = PetscMin(t,localwolfemin);
637:       }
638:       localmax = PetscMax(t,localmax);
639:     }
640:   }
641:   ierr=VecRestoreArray(X,&x);
642:   ierr=VecRestoreArray(XL,&xl);
643:   ierr=VecRestoreArray(XU,&xu);
644:   ierr=VecRestoreArray(DX,&dx);
645:   ierr=PetscObjectGetComm((PetscObject)X,&comm);

647:   if (boundmin){
648:     MPI_Allreduce(&localmin,boundmin,1,MPIU_REAL,MPIU_MIN,comm);
649:     PetscInfo1(X,"Step Bound Info: Closest Bound: %g \n",(double)*boundmin);
650:   }
651:   if (wolfemin){
652:     MPI_Allreduce(&localwolfemin,wolfemin,1,MPIU_REAL,MPIU_MIN,comm);
653:     PetscInfo1(X,"Step Bound Info: Wolfe: %g \n",(double)*wolfemin);
654:   }
655:   if (boundmax) {
656:     MPI_Allreduce(&localmax,boundmax,1,MPIU_REAL,MPIU_MAX,comm);
657:     PetscInfo1(X,"Step Bound Info: Max: %g \n",(double)*boundmax);
658:   }
659:   return(0);
660: }

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

667:      Collective on Vec

669:      Input Parameters:
670: +      X  - vector with no negative entries
671: -      DX  - a step direction, can have negative, positive or zero entries

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

676:   Level: advanced
677:  @*/
678: PetscErrorCode VecStepMax(Vec X, Vec DX, PetscReal *step)
679: {
681:   PetscInt       i, nn;
682:   PetscReal      stepmax=PETSC_INFINITY;
683:   PetscScalar    *xx, *dx;
684:   MPI_Comm       comm;


690:   VecGetLocalSize(X,&nn);
691:   VecGetArray(X,&xx);
692:   VecGetArray(DX,&dx);
693:   for (i=0;i<nn;i++){
694:     if (PetscRealPart(xx[i]) < 0) SETERRQ(PETSC_COMM_SELF,1,"Vector must be positive");
695:     else if (PetscRealPart(dx[i])<0) stepmax=PetscMin(stepmax,PetscRealPart(-xx[i]/dx[i]));
696:   }
697:   VecRestoreArray(X,&xx);
698:   VecRestoreArray(DX,&dx);
699:   PetscObjectGetComm((PetscObject)X,&comm);
700:   MPI_Allreduce(&stepmax,step,1,MPIU_REAL,MPIU_MIN,comm);
701:   return(0);
702: }

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

709:   Logically Collective on v

711:   Input Parameter:
712: + v - the vector
713: - p - the exponent to use on each element

715:   Output Parameter:
716: . v - the vector

718:   Level: intermediate

720: @*/
721: PetscErrorCode VecPow(Vec v, PetscScalar p)
722: {
724:   PetscInt       n,i;
725:   PetscScalar    *v1;


730:   VecGetArray(v, &v1);
731:   VecGetLocalSize(v, &n);

733:   if (1.0 == p) {
734:   } else if (-1.0 == p) {
735:     for (i = 0; i < n; ++i){
736:       v1[i] = 1.0 / v1[i];
737:     }
738:   } else if (0.0 == p) {
739:     for (i = 0; i < n; ++i){
740:       /*  Not-a-number left alone
741:           Infinity set to one  */
742:       if (v1[i] == v1[i]) {
743:         v1[i] = 1.0;
744:       }
745:     }
746:   } else if (0.5 == p) {
747:     for (i = 0; i < n; ++i) {
748:       if (PetscRealPart(v1[i]) >= 0) {
749:         v1[i] = PetscSqrtScalar(v1[i]);
750:       } else {
751:         v1[i] = PETSC_INFINITY;
752:       }
753:     }
754:   } else if (-0.5 == p) {
755:     for (i = 0; i < n; ++i) {
756:       if (PetscRealPart(v1[i]) >= 0) {
757:         v1[i] = 1.0 / PetscSqrtScalar(v1[i]);
758:       } else {
759:         v1[i] = PETSC_INFINITY;
760:       }
761:     }
762:   } else if (2.0 == p) {
763:     for (i = 0; i < n; ++i){
764:       v1[i] *= v1[i];
765:     }
766:   } else if (-2.0 == p) {
767:     for (i = 0; i < n; ++i){
768:       v1[i] = 1.0 / (v1[i] * v1[i]);
769:     }
770:   } else {
771:     for (i = 0; i < n; ++i) {
772:       if (PetscRealPart(v1[i]) >= 0) {
773:         v1[i] = PetscPowScalar(v1[i], p);
774:       } else {
775:         v1[i] = PETSC_INFINITY;
776:       }
777:     }
778:   }
779:   VecRestoreArray(v,&v1);
780:   return(0);
781: }

785: /*@
786:   VecMedian - Computes the componentwise median of three vectors
787:   and stores the result in this vector.  Used primarily for projecting
788:   a vector within upper and lower bounds.

790:   Logically Collective

792:   Input Parameters:
793: . Vec1, Vec2, Vec3 - The three vectors

795:   Output Parameter:
796: . VMedian - The median vector

798:   Level: advanced
799: @*/
800: PetscErrorCode VecMedian(Vec Vec1, Vec Vec2, Vec Vec3, Vec VMedian)
801: {
803:   PetscInt       i,n,low1,low2,low3,low4,high1,high2,high3,high4;
804:   PetscScalar    *v1,*v2,*v3,*vmed;


812:   if (Vec1==Vec2 || Vec1==Vec3){
813:     ierr=VecCopy(Vec1,VMedian);
814:     return(0);
815:   }
816:   if (Vec2==Vec3){
817:     ierr=VecCopy(Vec2,VMedian);
818:     return(0);
819:   }


827:   VecGetOwnershipRange(Vec1, &low1, &high1);
828:   VecGetOwnershipRange(Vec2, &low2, &high2);
829:   VecGetOwnershipRange(Vec3, &low3, &high3);
830:   VecGetOwnershipRange(VMedian, &low4, &high4);
831:   if ( low1!= low2 || low1!= low3 || low1!= low4 || high1!= high2 || high1!= high3 || high1!= high4) SETERRQ(PETSC_COMM_SELF,1,"InCompatible vector local lengths");

833:   VecGetArray(Vec1,&v1);
834:   VecGetArray(Vec2,&v2);
835:   VecGetArray(Vec3,&v3);

837:   if ( VMedian != Vec1 && VMedian != Vec2 && VMedian != Vec3){
838:     VecGetArray(VMedian,&vmed);
839:   } else if ( VMedian==Vec1 ){
840:     vmed=v1;
841:   } else if ( VMedian==Vec2 ){
842:     vmed=v2;
843:   } else {
844:     vmed=v3;
845:   }

847:   ierr=VecGetLocalSize(Vec1,&n);

849:   for (i=0;i<n;i++){
850:     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])));
851:   }

853:   VecRestoreArray(Vec1,&v1);
854:   VecRestoreArray(Vec2,&v2);
855:   VecRestoreArray(Vec3,&v2);

857:   if (VMedian!=Vec1 && VMedian != Vec2 && VMedian != Vec3){
858:     VecRestoreArray(VMedian,&vmed);
859:   }
860:   return(0);
861: }