Actual source code: projection.c
petsc-3.8.4 2018-03-24
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,< );
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, > );
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: }