Actual source code: projection.c
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,<);
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,>);
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: }
258: /*@
259: VecWhichBetweenOrEqual - Creates an index set containing the indices
260: where VecLow <= V <= VecHigh
262: Collective on S
264: Input Parameters:
265: + VecLow - lower bound
266: . V - Vector to compare
267: - VecHigh - higher bound
269: OutputParameter:
270: . S - The index set containing the indices i where veclow[i] <= v[i] <= vechigh[i]
272: Level: advanced
273: @*/
275: PetscErrorCode VecWhichBetweenOrEqual(Vec VecLow, Vec V, Vec VecHigh, IS * S)
276: {
277: PetscErrorCode ierr;
278: PetscInt i,n_vm=0;
279: PetscInt n,low,high;
280: PetscInt *vm=NULL;
281: const PetscScalar *v1,*v2,*vmiddle;
289: VecCheckSameSize(V,2,VecLow,1);
290: VecCheckSameSize(V,2,VecHigh,3);
292: VecGetOwnershipRange(VecLow,&low,&high);
293: VecGetLocalSize(VecLow,&n);
294: if (n>0) {
295: VecGetArrayRead(VecLow,&v1);
296: if (VecLow != VecHigh) {
297: VecGetArrayRead(VecHigh,&v2);
298: } else {
299: v2=v1;
300: }
301: if (V != VecLow && V != VecHigh) {
302: VecGetArrayRead(V,&vmiddle);
303: } else if (V==VecLow) {
304: vmiddle=v1;
305: } else {
306: vmiddle =v2;
307: }
309: PetscMalloc1(n,&vm);
311: for (i=0; i<n; ++i) {
312: if (PetscRealPart(v1[i]) <= PetscRealPart(vmiddle[i]) && PetscRealPart(vmiddle[i]) <= PetscRealPart(v2[i])) {vm[n_vm]=low+i; ++n_vm;}
313: }
315: VecRestoreArrayRead(VecLow,&v1);
316: if (VecLow != VecHigh) {
317: VecRestoreArrayRead(VecHigh,&v2);
318: }
319: if (V != VecLow && V != VecHigh) {
320: VecRestoreArrayRead(V,&vmiddle);
321: }
322: }
323: ISCreateGeneral(PetscObjectComm((PetscObject)V),n_vm,vm,PETSC_OWN_POINTER,S);
324: return(0);
325: }
327: /*@
328: VecWhichInactive - Creates an index set containing the indices
329: where one of the following holds:
330: a) VecLow(i) < V(i) < VecHigh(i)
331: b) VecLow(i) = V(i) and D(i) <= 0 (< 0 when Strong is true)
332: c) VecHigh(i) = V(i) and D(i) >= 0 (> 0 when Strong is true)
334: Collective on S
336: Input Parameters:
337: + VecLow - lower bound
338: . V - Vector to compare
339: . D - Direction to compare
340: . VecHigh - higher bound
341: - Strong - indicator for applying strongly inactive test
343: OutputParameter:
344: . S - The index set containing the indices i where the bound is inactive
346: Level: advanced
347: @*/
349: PetscErrorCode VecWhichInactive(Vec VecLow, Vec V, Vec D, Vec VecHigh, PetscBool Strong, IS * S)
350: {
351: PetscErrorCode ierr;
352: PetscInt i,n_vm=0;
353: PetscInt n,low,high;
354: PetscInt *vm=NULL;
355: const PetscScalar *v1,*v2,*v,*d;
365: VecCheckSameSize(V,2,VecLow,1);
366: VecCheckSameSize(V,2,D,3);
367: VecCheckSameSize(V,2,VecHigh,4);
369: VecGetOwnershipRange(VecLow,&low,&high);
370: VecGetLocalSize(VecLow,&n);
371: if (n>0) {
372: VecGetArrayRead(VecLow,&v1);
373: if (VecLow != VecHigh) {
374: VecGetArrayRead(VecHigh,&v2);
375: } else {
376: v2=v1;
377: }
378: if (V != VecLow && V != VecHigh) {
379: VecGetArrayRead(V,&v);
380: } else if (V==VecLow) {
381: v = v1;
382: } else {
383: v = v2;
384: }
385: if (D != VecLow && D != VecHigh && D != V) {
386: VecGetArrayRead(D,&d);
387: } else if (D==VecLow) {
388: d = v1;
389: } else if (D==VecHigh) {
390: d = v2;
391: } else {
392: d = v;
393: }
395: PetscMalloc1(n,&vm);
397: if (Strong) {
398: for (i=0; i<n; ++i) {
399: if (PetscRealPart(v1[i]) < PetscRealPart(v[i]) && PetscRealPart(v[i]) < PetscRealPart(v2[i])) {
400: vm[n_vm]=low+i; ++n_vm;
401: } else if (PetscRealPart(v1[i]) == PetscRealPart(v[i]) && PetscRealPart(d[i]) < 0) {
402: vm[n_vm]=low+i; ++n_vm;
403: } else if (PetscRealPart(v2[i]) == PetscRealPart(v[i]) && PetscRealPart(d[i]) > 0) {
404: vm[n_vm]=low+i; ++n_vm;
405: }
406: }
407: } else {
408: for (i=0; i<n; ++i) {
409: if (PetscRealPart(v1[i]) < PetscRealPart(v[i]) && PetscRealPart(v[i]) < PetscRealPart(v2[i])) {
410: vm[n_vm]=low+i; ++n_vm;
411: } else if (PetscRealPart(v1[i]) == PetscRealPart(v[i]) && PetscRealPart(d[i]) <= 0) {
412: vm[n_vm]=low+i; ++n_vm;
413: } else if (PetscRealPart(v2[i]) == PetscRealPart(v[i]) && PetscRealPart(d[i]) >= 0) {
414: vm[n_vm]=low+i; ++n_vm;
415: }
416: }
417: }
419: VecRestoreArrayRead(VecLow,&v1);
420: if (VecLow != VecHigh) {
421: VecRestoreArrayRead(VecHigh,&v2);
422: }
423: if (V != VecLow && V != VecHigh) {
424: VecRestoreArrayRead(V,&v);
425: }
426: if (D != VecLow && D != VecHigh && D != V) {
427: VecRestoreArrayRead(D,&d);
428: }
429: }
430: ISCreateGeneral(PetscObjectComm((PetscObject)V),n_vm,vm,PETSC_OWN_POINTER,S);
431: return(0);
432: }
434: /*@
435: VecISAXPY - Adds a reduced vector to the appropriate elements of a full-space vector.
436: vfull[is[i]] += alpha*vreduced[i]
438: Input Parameters:
439: + vfull - the full-space vector
440: . is - the index set for the reduced space
441: . alpha - the scalar coefficient
442: - vreduced - the reduced-space vector
444: Output Parameters:
445: . vfull - the sum of the full-space vector and reduced-space vector
447: Notes:
448: The index set identifies entries in the global vector.
449: Negative indices are skipped; indices outside the ownership range of vfull will raise an error.
451: Level: advanced
453: .seealso: VecAXPY(), VecGetOwnershipRange()
454: @*/
455: PetscErrorCode VecISAXPY(Vec vfull, IS is, PetscScalar alpha, Vec vreduced)
456: {
457: PetscInt nfull,nreduced;
464: VecGetSize(vfull,&nfull);
465: VecGetSize(vreduced,&nreduced);
467: if (nfull == nreduced) { /* Also takes care of masked vectors */
468: VecAXPY(vfull,alpha,vreduced);
469: } else {
470: PetscScalar *y;
471: const PetscScalar *x;
472: PetscInt i,n,m,rstart,rend;
473: const PetscInt *id;
475: VecGetArray(vfull,&y);
476: VecGetArrayRead(vreduced,&x);
477: ISGetIndices(is,&id);
478: ISGetLocalSize(is,&n);
479: VecGetLocalSize(vreduced,&m);
480: if (m != n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"IS local length not equal to Vec local length");
481: VecGetOwnershipRange(vfull,&rstart,&rend);
482: y -= rstart;
483: if (alpha == 1.0) {
484: for (i=0; i<n; ++i) {
485: if (id[i] < 0) continue;
486: if (id[i] < rstart || id[i] >= rend) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only owned values supported");
487: y[id[i]] += x[i];
488: }
489: } else {
490: for (i=0; i<n; ++i) {
491: if (id[i] < 0) continue;
492: if (id[i] < rstart || id[i] >= rend) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only owned values supported");
493: y[id[i]] += alpha*x[i];
494: }
495: }
496: y += rstart;
497: ISRestoreIndices(is,&id);
498: VecRestoreArray(vfull,&y);
499: VecRestoreArrayRead(vreduced,&x);
500: }
501: return(0);
502: }
504: /*@
505: VecISCopy - Copies between a reduced vector and the appropriate elements of a full-space vector.
507: Input Parameters:
508: + vfull - the full-space vector
509: . is - the index set for the reduced space
510: . mode - the direction of copying, SCATTER_FORWARD or SCATTER_REVERSE
511: - vreduced - the reduced-space vector
513: Output Parameters:
514: . vfull - the sum of the full-space vector and reduced-space vector
516: Notes:
517: The index set identifies entries in the global vector.
518: Negative indices are skipped; indices outside the ownership range of vfull will raise an error.
520: mode == SCATTER_FORWARD: vfull[is[i]] = vreduced[i]
521: mode == SCATTER_REVERSE: vreduced[i] = vfull[is[i]]
523: Level: advanced
525: .seealso: VecAXPY(), VecGetOwnershipRange()
526: @*/
527: PetscErrorCode VecISCopy(Vec vfull, IS is, ScatterMode mode, Vec vreduced)
528: {
529: PetscInt nfull, nreduced;
536: VecGetSize(vfull, &nfull);
537: VecGetSize(vreduced, &nreduced);
539: if (nfull == nreduced) { /* Also takes care of masked vectors */
540: if (mode == SCATTER_FORWARD) {
541: VecCopy(vreduced, vfull);
542: } else {
543: VecCopy(vfull, vreduced);
544: }
545: } else {
546: const PetscInt *id;
547: PetscInt i, n, m, rstart, rend;
549: ISGetIndices(is, &id);
550: ISGetLocalSize(is, &n);
551: VecGetLocalSize(vreduced, &m);
552: VecGetOwnershipRange(vfull, &rstart, &rend);
553: if (m != n) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "IS local length %D not equal to Vec local length %D", n, m);
554: if (mode == SCATTER_FORWARD) {
555: PetscScalar *y;
556: const PetscScalar *x;
558: VecGetArray(vfull, &y);
559: VecGetArrayRead(vreduced, &x);
560: y -= rstart;
561: for (i = 0; i < n; ++i) {
562: if (id[i] < 0) continue;
563: if (id[i] < rstart || id[i] >= rend) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only owned values supported");
564: y[id[i]] = x[i];
565: }
566: y += rstart;
567: VecRestoreArrayRead(vreduced, &x);
568: VecRestoreArray(vfull, &y);
569: } else if (mode == SCATTER_REVERSE) {
570: PetscScalar *x;
571: const PetscScalar *y;
573: VecGetArrayRead(vfull, &y);
574: VecGetArray(vreduced, &x);
575: for (i = 0; i < n; ++i) {
576: if (id[i] < 0) continue;
577: if (id[i] < rstart || id[i] >= rend) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only owned values supported");
578: x[i] = y[id[i]-rstart];
579: }
580: VecRestoreArray(vreduced, &x);
581: VecRestoreArrayRead(vfull, &y);
582: } else SETERRQ(PetscObjectComm((PetscObject) vfull), PETSC_ERR_ARG_WRONG, "Only forward or reverse modes are legal");
583: ISRestoreIndices(is, &id);
584: }
585: return(0);
586: }
588: /*@
589: ISComplementVec - Creates the complement of the index set relative to a layout defined by a Vec
591: Collective on IS
593: Input Parameters:
594: + S - a PETSc IS
595: - V - the reference vector space
597: Output Parameter:
598: . T - the complement of S
600: Level: advanced
602: .seealso: ISCreateGeneral()
603: @*/
604: PetscErrorCode ISComplementVec(IS S, Vec V, IS *T)
605: {
607: PetscInt start, end;
610: VecGetOwnershipRange(V,&start,&end);
611: ISComplement(S,start,end,T);
612: return(0);
613: }
615: /*@
616: VecISSet - Sets the elements of a vector, specified by an index set, to a constant
618: Input Parameters:
619: + V - the vector
620: . S - index set for the locations in the vector
621: - c - the constant
623: Notes:
624: The index set identifies entries in the global vector.
625: Negative indices are skipped; indices outside the ownership range of V will raise an error.
627: Level: advanced
629: .seealso: VecSet(), VecGetOwnershipRange()
630: @*/
631: PetscErrorCode VecISSet(Vec V,IS S, PetscScalar c)
632: {
634: PetscInt nloc,low,high,i;
635: const PetscInt *s;
636: PetscScalar *v;
639: if (!S) return(0); /* simply return with no-op if the index set is NULL */
644: VecGetOwnershipRange(V,&low,&high);
645: ISGetLocalSize(S,&nloc);
646: ISGetIndices(S,&s);
647: VecGetArray(V,&v);
648: for (i=0; i<nloc; ++i) {
649: if (s[i] < 0) continue;
650: if (s[i] < low || s[i] >= high) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only owned values supported");
651: v[s[i]-low] = c;
652: }
653: ISRestoreIndices(S,&s);
654: VecRestoreArray(V,&v);
655: return(0);
656: }
658: #if !defined(PETSC_USE_COMPLEX)
659: /*@C
660: VecBoundGradientProjection - Projects vector according to this definition.
661: If XL[i] < X[i] < XU[i], then GP[i] = G[i];
662: If X[i] <= XL[i], then GP[i] = min(G[i],0);
663: If X[i] >= XU[i], then GP[i] = max(G[i],0);
665: Input Parameters:
666: + G - current gradient vector
667: . X - current solution vector with XL[i] <= X[i] <= XU[i]
668: . XL - lower bounds
669: - XU - upper bounds
671: Output Parameter:
672: . GP - gradient projection vector
674: Notes:
675: GP may be the same vector as G
677: Level: advanced
678: @*/
679: PetscErrorCode VecBoundGradientProjection(Vec G, Vec X, Vec XL, Vec XU, Vec GP)
680: {
682: PetscErrorCode ierr;
683: PetscInt n,i;
684: const PetscReal *xptr,*xlptr,*xuptr;
685: PetscReal *gptr,*gpptr;
686: PetscReal xval,gpval;
688: /* Project variables at the lower and upper bound */
696: VecGetLocalSize(X,&n);
698: VecGetArrayRead(X,&xptr);
699: VecGetArrayRead(XL,&xlptr);
700: VecGetArrayRead(XU,&xuptr);
701: VecGetArrayPair(G,GP,&gptr,&gpptr);
703: for (i=0; i<n; ++i) {
704: gpval = gptr[i]; xval = xptr[i];
705: if (gpval>0.0 && xval<=xlptr[i]) {
706: gpval = 0.0;
707: } else if (gpval<0.0 && xval>=xuptr[i]) {
708: gpval = 0.0;
709: }
710: gpptr[i] = gpval;
711: }
713: VecRestoreArrayRead(X,&xptr);
714: VecRestoreArrayRead(XL,&xlptr);
715: VecRestoreArrayRead(XU,&xuptr);
716: VecRestoreArrayPair(G,GP,&gptr,&gpptr);
717: return(0);
718: }
719: #endif
721: /*@
722: VecStepMaxBounded - See below
724: Collective on Vec
726: Input Parameters:
727: + X - vector with no negative entries
728: . XL - lower bounds
729: . XU - upper bounds
730: - DX - step direction, can have negative, positive or zero entries
732: Output Parameter:
733: . stepmax - minimum value so that X[i] + stepmax*DX[i] <= XL[i] or XU[i] <= X[i] + stepmax*DX[i]
735: Level: intermediate
737: @*/
738: PetscErrorCode VecStepMaxBounded(Vec X, Vec DX, Vec XL, Vec XU, PetscReal *stepmax)
739: {
740: PetscErrorCode ierr;
741: PetscInt i,nn;
742: const PetscScalar *xx,*dx,*xl,*xu;
743: PetscReal localmax=0;
751: VecGetArrayRead(X,&xx);
752: VecGetArrayRead(XL,&xl);
753: VecGetArrayRead(XU,&xu);
754: VecGetArrayRead(DX,&dx);
755: VecGetLocalSize(X,&nn);
756: for (i=0;i<nn;i++) {
757: if (PetscRealPart(dx[i]) > 0) {
758: localmax=PetscMax(localmax,PetscRealPart((xu[i]-xx[i])/dx[i]));
759: } else if (PetscRealPart(dx[i])<0) {
760: localmax=PetscMax(localmax,PetscRealPart((xl[i]-xx[i])/dx[i]));
761: }
762: }
763: VecRestoreArrayRead(X,&xx);
764: VecRestoreArrayRead(XL,&xl);
765: VecRestoreArrayRead(XU,&xu);
766: VecRestoreArrayRead(DX,&dx);
767: MPIU_Allreduce(&localmax,stepmax,1,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)X));
768: return(0);
769: }
771: /*@
772: VecStepBoundInfo - See below
774: Collective on Vec
776: Input Parameters:
777: + X - vector with no negative entries
778: . XL - lower bounds
779: . XU - upper bounds
780: - DX - step direction, can have negative, positive or zero entries
782: Output Parameters:
783: + boundmin - (may be NULL this it is not computed) maximum value so that XL[i] <= X[i] + boundmax*DX[i] <= XU[i]
784: . wolfemin - (may be NULL this it is not computed)
785: - 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]
787: Notes:
788: For complex numbers only compares the real part
790: Level: advanced
791: @*/
792: PetscErrorCode VecStepBoundInfo(Vec X, Vec DX, Vec XL, Vec XU, PetscReal *boundmin, PetscReal *wolfemin, PetscReal *boundmax)
793: {
794: PetscErrorCode ierr;
795: PetscInt n,i;
796: const PetscScalar *x,*xl,*xu,*dx;
797: PetscReal t;
798: PetscReal localmin=PETSC_INFINITY,localwolfemin=PETSC_INFINITY,localmax=-1;
799: MPI_Comm comm;
807: VecGetArrayRead(X,&x);
808: VecGetArrayRead(XL,&xl);
809: VecGetArrayRead(XU,&xu);
810: VecGetArrayRead(DX,&dx);
811: VecGetLocalSize(X,&n);
812: for (i=0; i<n; ++i) {
813: if (PetscRealPart(dx[i])>0 && PetscRealPart(xu[i]) < PETSC_INFINITY) {
814: t=PetscRealPart((xu[i]-x[i])/dx[i]);
815: localmin=PetscMin(t,localmin);
816: if (localmin>0) {
817: localwolfemin = PetscMin(t,localwolfemin);
818: }
819: localmax = PetscMax(t,localmax);
820: } else if (PetscRealPart(dx[i])<0 && PetscRealPart(xl[i]) > PETSC_NINFINITY) {
821: t=PetscRealPart((xl[i]-x[i])/dx[i]);
822: localmin = PetscMin(t,localmin);
823: if (localmin>0) {
824: localwolfemin = PetscMin(t,localwolfemin);
825: }
826: localmax = PetscMax(t,localmax);
827: }
828: }
830: VecRestoreArrayRead(X,&x);
831: VecRestoreArrayRead(XL,&xl);
832: VecRestoreArrayRead(XU,&xu);
833: VecRestoreArrayRead(DX,&dx);
834: PetscObjectGetComm((PetscObject)X,&comm);
836: if (boundmin) {
837: MPIU_Allreduce(&localmin,boundmin,1,MPIU_REAL,MPIU_MIN,comm);
838: PetscInfo1(X,"Step Bound Info: Closest Bound: %20.19e\n",(double)*boundmin);
839: }
840: if (wolfemin) {
841: MPIU_Allreduce(&localwolfemin,wolfemin,1,MPIU_REAL,MPIU_MIN,comm);
842: PetscInfo1(X,"Step Bound Info: Wolfe: %20.19e\n",(double)*wolfemin);
843: }
844: if (boundmax) {
845: MPIU_Allreduce(&localmax,boundmax,1,MPIU_REAL,MPIU_MAX,comm);
846: if (*boundmax < 0) *boundmax=PETSC_INFINITY;
847: PetscInfo1(X,"Step Bound Info: Max: %20.19e\n",(double)*boundmax);
848: }
849: return(0);
850: }
852: /*@
853: VecStepMax - Returns the largest value so that x[i] + step*DX[i] >= 0 for all i
855: Collective on Vec
857: Input Parameters:
858: + X - vector with no negative entries
859: - DX - a step direction, can have negative, positive or zero entries
861: Output Parameter:
862: . step - largest value such that x[i] + step*DX[i] >= 0 for all i
864: Notes:
865: For complex numbers only compares the real part
867: Level: advanced
868: @*/
869: PetscErrorCode VecStepMax(Vec X, Vec DX, PetscReal *step)
870: {
871: PetscErrorCode ierr;
872: PetscInt i,nn;
873: PetscReal stepmax=PETSC_INFINITY;
874: const PetscScalar *xx,*dx;
880: VecGetLocalSize(X,&nn);
881: VecGetArrayRead(X,&xx);
882: VecGetArrayRead(DX,&dx);
883: for (i=0;i<nn;++i) {
884: if (PetscRealPart(xx[i]) < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Vector must be positive");
885: else if (PetscRealPart(dx[i])<0) stepmax=PetscMin(stepmax,PetscRealPart(-xx[i]/dx[i]));
886: }
887: VecRestoreArrayRead(X,&xx);
888: VecRestoreArrayRead(DX,&dx);
889: MPIU_Allreduce(&stepmax,step,1,MPIU_REAL,MPIU_MIN,PetscObjectComm((PetscObject)X));
890: return(0);
891: }
893: /*@
894: VecPow - Replaces each component of a vector by x_i^p
896: Logically Collective on v
898: Input Parameters:
899: + v - the vector
900: - p - the exponent to use on each element
902: Level: intermediate
904: @*/
905: PetscErrorCode VecPow(Vec v, PetscScalar p)
906: {
908: PetscInt n,i;
909: PetscScalar *v1;
914: VecGetArray(v,&v1);
915: VecGetLocalSize(v,&n);
917: if (1.0 == p) {
918: } else if (-1.0 == p) {
919: for (i = 0; i < n; ++i) {
920: v1[i] = 1.0 / v1[i];
921: }
922: } else if (0.0 == p) {
923: for (i = 0; i < n; ++i) {
924: /* Not-a-number left alone
925: Infinity set to one */
926: if (v1[i] == v1[i]) {
927: v1[i] = 1.0;
928: }
929: }
930: } else if (0.5 == p) {
931: for (i = 0; i < n; ++i) {
932: if (PetscRealPart(v1[i]) >= 0) {
933: v1[i] = PetscSqrtScalar(v1[i]);
934: } else {
935: v1[i] = PETSC_INFINITY;
936: }
937: }
938: } else if (-0.5 == p) {
939: for (i = 0; i < n; ++i) {
940: if (PetscRealPart(v1[i]) >= 0) {
941: v1[i] = 1.0 / PetscSqrtScalar(v1[i]);
942: } else {
943: v1[i] = PETSC_INFINITY;
944: }
945: }
946: } else if (2.0 == p) {
947: for (i = 0; i < n; ++i) {
948: v1[i] *= v1[i];
949: }
950: } else if (-2.0 == p) {
951: for (i = 0; i < n; ++i) {
952: v1[i] = 1.0 / (v1[i] * v1[i]);
953: }
954: } else {
955: for (i = 0; i < n; ++i) {
956: if (PetscRealPart(v1[i]) >= 0) {
957: v1[i] = PetscPowScalar(v1[i],p);
958: } else {
959: v1[i] = PETSC_INFINITY;
960: }
961: }
962: }
963: VecRestoreArray(v,&v1);
964: return(0);
965: }
967: /*@
968: VecMedian - Computes the componentwise median of three vectors
969: and stores the result in this vector. Used primarily for projecting
970: a vector within upper and lower bounds.
972: Logically Collective
974: Input Parameters:
975: + Vec1 - The first vector
976: . Vec2 - The second vector
977: - Vec3 - The third vector
979: Output Parameter:
980: . VMedian - The median vector (this can be any one of the input vectors)
982: Level: advanced
983: @*/
984: PetscErrorCode VecMedian(Vec Vec1, Vec Vec2, Vec Vec3, Vec VMedian)
985: {
986: PetscErrorCode ierr;
987: PetscInt i,n,low1,high1;
988: const PetscScalar *v1,*v2,*v3;
989: PetscScalar *vmed;
997: if (Vec1==Vec2 || Vec1==Vec3) {
998: VecCopy(Vec1,VMedian);
999: return(0);
1000: }
1001: if (Vec2==Vec3) {
1002: VecCopy(Vec2,VMedian);
1003: return(0);
1004: }
1006: /* Assert that Vec1 != Vec2 and Vec2 != Vec3 */
1017: VecCheckSameSize(Vec1,1,Vec2,2);
1018: VecCheckSameSize(Vec1,1,Vec3,3);
1019: VecCheckSameSize(Vec1,1,VMedian,4);
1021: VecGetOwnershipRange(Vec1,&low1,&high1);
1022: VecGetLocalSize(Vec1,&n);
1023: if (n>0) {
1024: VecGetArray(VMedian,&vmed);
1025: if (Vec1 != VMedian) {
1026: VecGetArrayRead(Vec1,&v1);
1027: } else {
1028: v1=vmed;
1029: }
1030: if (Vec2 != VMedian) {
1031: VecGetArrayRead(Vec2,&v2);
1032: } else {
1033: v2=vmed;
1034: }
1035: if (Vec3 != VMedian) {
1036: VecGetArrayRead(Vec3,&v3);
1037: } else {
1038: v3=vmed;
1039: }
1041: for (i=0;i<n;++i) {
1042: 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])));
1043: }
1045: VecRestoreArray(VMedian,&vmed);
1046: if (VMedian != Vec1) {
1047: VecRestoreArrayRead(Vec1,&v1);
1048: }
1049: if (VMedian != Vec2) {
1050: VecRestoreArrayRead(Vec2,&v2);
1051: }
1052: if (VMedian != Vec3) {
1053: VecRestoreArrayRead(Vec3,&v3);
1054: }
1055: }
1056: return(0);
1057: }