Actual source code: projection.c
petsc-3.7.5 2017-01-01
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,< );
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, > );
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 = NULL,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: MPIU_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=-1;
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 && PetscRealPart(xu[i]) < PETSC_INFINITY) {
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 && PetscRealPart(xl[i]) > PETSC_NINFINITY) {
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: }
642: ierr=VecRestoreArray(X,&x);
643: ierr=VecRestoreArray(XL,&xl);
644: ierr=VecRestoreArray(XU,&xu);
645: ierr=VecRestoreArray(DX,&dx);
646: ierr=PetscObjectGetComm((PetscObject)X,&comm);
648: if (boundmin){
649: MPIU_Allreduce(&localmin,boundmin,1,MPIU_REAL,MPIU_MIN,comm);
650: PetscInfo1(X,"Step Bound Info: Closest Bound: %g \n",(double)*boundmin);
651: }
652: if (wolfemin){
653: MPIU_Allreduce(&localwolfemin,wolfemin,1,MPIU_REAL,MPIU_MIN,comm);
654: PetscInfo1(X,"Step Bound Info: Wolfe: %g \n",(double)*wolfemin);
655: }
656: if (boundmax) {
657: MPIU_Allreduce(&localmax,boundmax,1,MPIU_REAL,MPIU_MAX,comm);
658: if (*boundmax < 0) *boundmax=PETSC_INFINITY;
659: PetscInfo1(X,"Step Bound Info: Max: %g \n",(double)*boundmax);
660: }
661: return(0);
662: }
666: /*@
667: VecStepMax - Returns the largest value so that x[i] + step*DX[i] >= 0 for all i
669: Collective on Vec
671: Input Parameters:
672: + X - vector with no negative entries
673: - DX - a step direction, can have negative, positive or zero entries
675: Output Parameter:
676: . step - largest value such that x[i] + step*DX[i] >= 0 for all i
678: Level: advanced
679: @*/
680: PetscErrorCode VecStepMax(Vec X, Vec DX, PetscReal *step)
681: {
683: PetscInt i, nn;
684: PetscReal stepmax=PETSC_INFINITY;
685: PetscScalar *xx, *dx;
686: MPI_Comm comm;
692: VecGetLocalSize(X,&nn);
693: VecGetArray(X,&xx);
694: VecGetArray(DX,&dx);
695: for (i=0;i<nn;i++){
696: if (PetscRealPart(xx[i]) < 0) SETERRQ(PETSC_COMM_SELF,1,"Vector must be positive");
697: else if (PetscRealPart(dx[i])<0) stepmax=PetscMin(stepmax,PetscRealPart(-xx[i]/dx[i]));
698: }
699: VecRestoreArray(X,&xx);
700: VecRestoreArray(DX,&dx);
701: PetscObjectGetComm((PetscObject)X,&comm);
702: MPIU_Allreduce(&stepmax,step,1,MPIU_REAL,MPIU_MIN,comm);
703: return(0);
704: }
708: /*@
709: VecPow - Replaces each component of a vector by x_i^p
711: Logically Collective on v
713: Input Parameter:
714: + v - the vector
715: - p - the exponent to use on each element
717: Output Parameter:
718: . v - the vector
720: Level: intermediate
722: @*/
723: PetscErrorCode VecPow(Vec v, PetscScalar p)
724: {
726: PetscInt n,i;
727: PetscScalar *v1;
732: VecGetArray(v, &v1);
733: VecGetLocalSize(v, &n);
735: if (1.0 == p) {
736: } else if (-1.0 == p) {
737: for (i = 0; i < n; ++i){
738: v1[i] = 1.0 / v1[i];
739: }
740: } else if (0.0 == p) {
741: for (i = 0; i < n; ++i){
742: /* Not-a-number left alone
743: Infinity set to one */
744: if (v1[i] == v1[i]) {
745: v1[i] = 1.0;
746: }
747: }
748: } else if (0.5 == p) {
749: for (i = 0; i < n; ++i) {
750: if (PetscRealPart(v1[i]) >= 0) {
751: v1[i] = PetscSqrtScalar(v1[i]);
752: } else {
753: v1[i] = PETSC_INFINITY;
754: }
755: }
756: } else if (-0.5 == p) {
757: for (i = 0; i < n; ++i) {
758: if (PetscRealPart(v1[i]) >= 0) {
759: v1[i] = 1.0 / PetscSqrtScalar(v1[i]);
760: } else {
761: v1[i] = PETSC_INFINITY;
762: }
763: }
764: } else if (2.0 == p) {
765: for (i = 0; i < n; ++i){
766: v1[i] *= v1[i];
767: }
768: } else if (-2.0 == p) {
769: for (i = 0; i < n; ++i){
770: v1[i] = 1.0 / (v1[i] * v1[i]);
771: }
772: } else {
773: for (i = 0; i < n; ++i) {
774: if (PetscRealPart(v1[i]) >= 0) {
775: v1[i] = PetscPowScalar(v1[i], p);
776: } else {
777: v1[i] = PETSC_INFINITY;
778: }
779: }
780: }
781: VecRestoreArray(v,&v1);
782: return(0);
783: }
787: /*@
788: VecMedian - Computes the componentwise median of three vectors
789: and stores the result in this vector. Used primarily for projecting
790: a vector within upper and lower bounds.
792: Logically Collective
794: Input Parameters:
795: . Vec1, Vec2, Vec3 - The three vectors
797: Output Parameter:
798: . VMedian - The median vector
800: Level: advanced
801: @*/
802: PetscErrorCode VecMedian(Vec Vec1, Vec Vec2, Vec Vec3, Vec VMedian)
803: {
805: PetscInt i,n,low1,low2,low3,low4,high1,high2,high3,high4;
806: PetscScalar *v1,*v2,*v3,*vmed;
814: if (Vec1==Vec2 || Vec1==Vec3){
815: ierr=VecCopy(Vec1,VMedian);
816: return(0);
817: }
818: if (Vec2==Vec3){
819: ierr=VecCopy(Vec2,VMedian);
820: return(0);
821: }
829: VecGetOwnershipRange(Vec1, &low1, &high1);
830: VecGetOwnershipRange(Vec2, &low2, &high2);
831: VecGetOwnershipRange(Vec3, &low3, &high3);
832: VecGetOwnershipRange(VMedian, &low4, &high4);
833: if ( low1!= low2 || low1!= low3 || low1!= low4 || high1!= high2 || high1!= high3 || high1!= high4) SETERRQ(PETSC_COMM_SELF,1,"InCompatible vector local lengths");
835: VecGetArray(Vec1,&v1);
836: VecGetArray(Vec2,&v2);
837: VecGetArray(Vec3,&v3);
839: if ( VMedian != Vec1 && VMedian != Vec2 && VMedian != Vec3){
840: VecGetArray(VMedian,&vmed);
841: } else if ( VMedian==Vec1 ){
842: vmed=v1;
843: } else if ( VMedian==Vec2 ){
844: vmed=v2;
845: } else {
846: vmed=v3;
847: }
849: ierr=VecGetLocalSize(Vec1,&n);
851: for (i=0;i<n;i++){
852: 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])));
853: }
855: VecRestoreArray(Vec1,&v1);
856: VecRestoreArray(Vec2,&v2);
857: VecRestoreArray(Vec3,&v3);
859: if (VMedian!=Vec1 && VMedian != Vec2 && VMedian != Vec3){
860: VecRestoreArray(VMedian,&vmed);
861: }
862: return(0);
863: }