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: PetscInt i, n_same = 0;
23: PetscInt n, low, high;
24: PetscInt *same = NULL;
25: const PetscScalar *v1, *v2;
30: VecCheckSameSize(Vec1, 1, Vec2, 2);
32: VecGetOwnershipRange(Vec1, &low, &high);
33: VecGetLocalSize(Vec1, &n);
34: if (n > 0) {
35: if (Vec1 == Vec2) {
36: VecGetArrayRead(Vec1, &v1);
37: v2 = v1;
38: } else {
39: VecGetArrayRead(Vec1, &v1);
40: VecGetArrayRead(Vec2, &v2);
41: }
43: PetscMalloc1(n, &same);
45: for (i = 0; i < n; ++i) {
46: if (v1[i] == v2[i]) {
47: same[n_same] = low + i;
48: ++n_same;
49: }
50: }
52: if (Vec1 == Vec2) {
53: VecRestoreArrayRead(Vec1, &v1);
54: } else {
55: VecRestoreArrayRead(Vec1, &v1);
56: VecRestoreArrayRead(Vec2, &v2);
57: }
58: }
59: ISCreateGeneral(PetscObjectComm((PetscObject)Vec1), n_same, same, PETSC_OWN_POINTER, S);
60: return 0;
61: }
63: /*@
64: VecWhichLessThan - Creates an index set containing the indices
65: where the vectors Vec1 < Vec2
67: Collective on S
69: Input Parameters:
70: . Vec1, Vec2 - the two vectors to compare
72: OutputParameter:
73: . S - The index set containing the indices i where vec1[i] < vec2[i]
75: Notes:
76: The two vectors must have the same parallel layout
78: For complex numbers this only compares the real part
80: Level: advanced
81: @*/
82: PetscErrorCode VecWhichLessThan(Vec Vec1, Vec Vec2, IS *S)
83: {
84: PetscInt i, n_lt = 0;
85: PetscInt n, low, high;
86: PetscInt *lt = NULL;
87: const PetscScalar *v1, *v2;
92: VecCheckSameSize(Vec1, 1, Vec2, 2);
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: }
105: PetscMalloc1(n, <);
107: for (i = 0; i < n; ++i) {
108: if (PetscRealPart(v1[i]) < PetscRealPart(v2[i])) {
109: lt[n_lt] = low + i;
110: ++n_lt;
111: }
112: }
114: if (Vec1 == Vec2) {
115: VecRestoreArrayRead(Vec1, &v1);
116: } else {
117: VecRestoreArrayRead(Vec1, &v1);
118: VecRestoreArrayRead(Vec2, &v2);
119: }
120: }
121: ISCreateGeneral(PetscObjectComm((PetscObject)Vec1), n_lt, lt, PETSC_OWN_POINTER, S);
122: return 0;
123: }
125: /*@
126: VecWhichGreaterThan - Creates an index set containing the indices
127: where the vectors Vec1 > Vec2
129: Collective on S
131: Input Parameters:
132: . Vec1, Vec2 - the two vectors to compare
134: OutputParameter:
135: . S - The index set containing the indices i where vec1[i] > vec2[i]
137: Notes:
138: The two vectors must have the same parallel layout
140: For complex numbers this only compares the real part
142: Level: advanced
143: @*/
144: PetscErrorCode VecWhichGreaterThan(Vec Vec1, Vec Vec2, IS *S)
145: {
146: PetscInt i, n_gt = 0;
147: PetscInt n, low, high;
148: PetscInt *gt = NULL;
149: 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])) {
171: gt[n_gt] = low + i;
172: ++n_gt;
173: }
174: }
176: if (Vec1 == Vec2) {
177: VecRestoreArrayRead(Vec1, &v1);
178: } else {
179: VecRestoreArrayRead(Vec1, &v1);
180: VecRestoreArrayRead(Vec2, &v2);
181: }
182: }
183: ISCreateGeneral(PetscObjectComm((PetscObject)Vec1), n_gt, gt, PETSC_OWN_POINTER, S);
184: return 0;
185: }
187: /*@
188: VecWhichBetween - Creates an index set containing the indices
189: where VecLow < V < VecHigh
191: Collective on S
193: Input Parameters:
194: + VecLow - lower bound
195: . V - Vector to compare
196: - VecHigh - higher bound
198: OutputParameter:
199: . S - The index set containing the indices i where veclow[i] < v[i] < vechigh[i]
201: Notes:
202: The vectors must have the same parallel layout
204: For complex numbers this only compares the real part
206: Level: advanced
207: @*/
208: PetscErrorCode VecWhichBetween(Vec VecLow, Vec V, Vec VecHigh, IS *S)
209: {
210: PetscInt i, n_vm = 0;
211: PetscInt n, low, high;
212: PetscInt *vm = NULL;
213: 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])) {
244: vm[n_vm] = low + i;
245: ++n_vm;
246: }
247: }
249: VecRestoreArrayRead(VecLow, &v1);
250: if (VecLow != VecHigh) VecRestoreArrayRead(VecHigh, &v2);
251: if (V != VecLow && V != VecHigh) VecRestoreArrayRead(V, &vmiddle);
252: }
253: ISCreateGeneral(PetscObjectComm((PetscObject)V), n_vm, vm, PETSC_OWN_POINTER, S);
254: return 0;
255: }
257: /*@
258: VecWhichBetweenOrEqual - Creates an index set containing the indices
259: where VecLow <= V <= VecHigh
261: Collective on S
263: Input Parameters:
264: + VecLow - lower bound
265: . V - Vector to compare
266: - VecHigh - higher bound
268: OutputParameter:
269: . S - The index set containing the indices i where veclow[i] <= v[i] <= vechigh[i]
271: Level: advanced
272: @*/
274: PetscErrorCode VecWhichBetweenOrEqual(Vec VecLow, Vec V, Vec VecHigh, IS *S)
275: {
276: PetscInt i, n_vm = 0;
277: PetscInt n, low, high;
278: PetscInt *vm = NULL;
279: const PetscScalar *v1, *v2, *vmiddle;
286: VecCheckSameSize(V, 2, VecLow, 1);
287: VecCheckSameSize(V, 2, VecHigh, 3);
289: VecGetOwnershipRange(VecLow, &low, &high);
290: VecGetLocalSize(VecLow, &n);
291: if (n > 0) {
292: VecGetArrayRead(VecLow, &v1);
293: if (VecLow != VecHigh) {
294: VecGetArrayRead(VecHigh, &v2);
295: } else {
296: v2 = v1;
297: }
298: if (V != VecLow && V != VecHigh) {
299: VecGetArrayRead(V, &vmiddle);
300: } else if (V == VecLow) {
301: vmiddle = v1;
302: } else {
303: vmiddle = v2;
304: }
306: PetscMalloc1(n, &vm);
308: for (i = 0; i < n; ++i) {
309: if (PetscRealPart(v1[i]) <= PetscRealPart(vmiddle[i]) && PetscRealPart(vmiddle[i]) <= PetscRealPart(v2[i])) {
310: vm[n_vm] = low + i;
311: ++n_vm;
312: }
313: }
315: VecRestoreArrayRead(VecLow, &v1);
316: if (VecLow != VecHigh) VecRestoreArrayRead(VecHigh, &v2);
317: if (V != VecLow && V != VecHigh) VecRestoreArrayRead(V, &vmiddle);
318: }
319: ISCreateGeneral(PetscObjectComm((PetscObject)V), n_vm, vm, PETSC_OWN_POINTER, S);
320: return 0;
321: }
323: /*@
324: VecWhichInactive - Creates an index set containing the indices
325: where one of the following holds:
326: a) VecLow(i) < V(i) < VecHigh(i)
327: b) VecLow(i) = V(i) and D(i) <= 0 (< 0 when Strong is true)
328: c) VecHigh(i) = V(i) and D(i) >= 0 (> 0 when Strong is true)
330: Collective on S
332: Input Parameters:
333: + VecLow - lower bound
334: . V - Vector to compare
335: . D - Direction to compare
336: . VecHigh - higher bound
337: - Strong - indicator for applying strongly inactive test
339: OutputParameter:
340: . S - The index set containing the indices i where the bound is inactive
342: Level: advanced
343: @*/
345: PetscErrorCode VecWhichInactive(Vec VecLow, Vec V, Vec D, Vec VecHigh, PetscBool Strong, IS *S)
346: {
347: PetscInt i, n_vm = 0;
348: PetscInt n, low, high;
349: PetscInt *vm = NULL;
350: const PetscScalar *v1, *v2, *v, *d;
352: if (!VecLow && !VecHigh) {
353: VecGetOwnershipRange(V, &low, &high);
354: ISCreateStride(PetscObjectComm((PetscObject)V), high - low, low, 1, S);
355: return 0;
356: }
364: VecCheckSameSize(V, 2, VecLow, 1);
365: VecCheckSameSize(V, 2, D, 3);
366: VecCheckSameSize(V, 2, VecHigh, 4);
368: VecGetOwnershipRange(VecLow, &low, &high);
369: VecGetLocalSize(VecLow, &n);
370: if (n > 0) {
371: VecGetArrayRead(VecLow, &v1);
372: if (VecLow != VecHigh) {
373: VecGetArrayRead(VecHigh, &v2);
374: } else {
375: v2 = v1;
376: }
377: if (V != VecLow && V != VecHigh) {
378: VecGetArrayRead(V, &v);
379: } else if (V == VecLow) {
380: v = v1;
381: } else {
382: v = v2;
383: }
384: if (D != VecLow && D != VecHigh && D != V) {
385: VecGetArrayRead(D, &d);
386: } else if (D == VecLow) {
387: d = v1;
388: } else if (D == VecHigh) {
389: d = v2;
390: } else {
391: d = v;
392: }
394: PetscMalloc1(n, &vm);
396: if (Strong) {
397: for (i = 0; i < n; ++i) {
398: if (PetscRealPart(v1[i]) < PetscRealPart(v[i]) && PetscRealPart(v[i]) < PetscRealPart(v2[i])) {
399: vm[n_vm] = low + i;
400: ++n_vm;
401: } else if (PetscRealPart(v1[i]) == PetscRealPart(v[i]) && PetscRealPart(d[i]) < 0) {
402: vm[n_vm] = low + i;
403: ++n_vm;
404: } else if (PetscRealPart(v2[i]) == PetscRealPart(v[i]) && PetscRealPart(d[i]) > 0) {
405: vm[n_vm] = low + i;
406: ++n_vm;
407: }
408: }
409: } else {
410: for (i = 0; i < n; ++i) {
411: if (PetscRealPart(v1[i]) < PetscRealPart(v[i]) && PetscRealPart(v[i]) < PetscRealPart(v2[i])) {
412: vm[n_vm] = low + i;
413: ++n_vm;
414: } else if (PetscRealPart(v1[i]) == PetscRealPart(v[i]) && PetscRealPart(d[i]) <= 0) {
415: vm[n_vm] = low + i;
416: ++n_vm;
417: } else if (PetscRealPart(v2[i]) == PetscRealPart(v[i]) && PetscRealPart(d[i]) >= 0) {
418: vm[n_vm] = low + i;
419: ++n_vm;
420: }
421: }
422: }
424: VecRestoreArrayRead(VecLow, &v1);
425: if (VecLow != VecHigh) VecRestoreArrayRead(VecHigh, &v2);
426: if (V != VecLow && V != VecHigh) VecRestoreArrayRead(V, &v);
427: if (D != VecLow && D != VecHigh && D != V) VecRestoreArrayRead(D, &d);
428: }
429: ISCreateGeneral(PetscObjectComm((PetscObject)V), n_vm, vm, PETSC_OWN_POINTER, S);
430: return 0;
431: }
433: /*@
434: VecISAXPY - Adds a reduced vector to the appropriate elements of a full-space vector.
435: vfull[is[i]] += alpha*vreduced[i]
437: Input Parameters:
438: + vfull - the full-space vector
439: . is - the index set for the reduced space
440: . alpha - the scalar coefficient
441: - vreduced - the reduced-space vector
443: Output Parameters:
444: . vfull - the sum of the full-space vector and reduced-space vector
446: Notes:
447: The index set identifies entries in the global vector.
448: Negative indices are skipped; indices outside the ownership range of vfull will raise an error.
450: Level: advanced
452: .seealso: `VecISCopy()`, `VecISSet()`, `VecAXPY()`
453: @*/
454: PetscErrorCode VecISAXPY(Vec vfull, IS is, PetscScalar alpha, Vec vreduced)
455: {
456: PetscInt nfull, nreduced;
461: VecGetSize(vfull, &nfull);
462: VecGetSize(vreduced, &nreduced);
464: if (nfull == nreduced) { /* Also takes care of masked vectors */
465: VecAXPY(vfull, alpha, vreduced);
466: } else {
467: PetscScalar *y;
468: const PetscScalar *x;
469: PetscInt i, n, m, rstart, rend;
470: const PetscInt *id;
472: VecGetArray(vfull, &y);
473: VecGetArrayRead(vreduced, &x);
474: ISGetIndices(is, &id);
475: ISGetLocalSize(is, &n);
476: VecGetLocalSize(vreduced, &m);
478: VecGetOwnershipRange(vfull, &rstart, &rend);
479: y -= rstart;
480: if (alpha == 1.0) {
481: for (i = 0; i < n; ++i) {
482: if (id[i] < 0) continue;
484: y[id[i]] += x[i];
485: }
486: } else {
487: for (i = 0; i < n; ++i) {
488: if (id[i] < 0) continue;
490: y[id[i]] += alpha * x[i];
491: }
492: }
493: y += rstart;
494: ISRestoreIndices(is, &id);
495: VecRestoreArray(vfull, &y);
496: VecRestoreArrayRead(vreduced, &x);
497: }
498: return 0;
499: }
501: /*@
502: VecISCopy - Copies between a reduced vector and the appropriate elements of a full-space vector.
504: Input Parameters:
505: + vfull - the full-space vector
506: . is - the index set for the reduced space
507: . mode - the direction of copying, SCATTER_FORWARD or SCATTER_REVERSE
508: - vreduced - the reduced-space vector
510: Output Parameters:
511: . vfull - the sum of the full-space vector and reduced-space vector
513: Notes:
514: The index set identifies entries in the global vector.
515: Negative indices are skipped; indices outside the ownership range of vfull will raise an error.
517: mode == SCATTER_FORWARD: vfull[is[i]] = vreduced[i]
518: mode == SCATTER_REVERSE: vreduced[i] = vfull[is[i]]
520: Level: advanced
522: .seealso: `VecISSet()`, `VecISAXPY()`, `VecCopy()`
523: @*/
524: PetscErrorCode VecISCopy(Vec vfull, IS is, ScatterMode mode, Vec vreduced)
525: {
526: PetscInt nfull, nreduced;
531: VecGetSize(vfull, &nfull);
532: VecGetSize(vreduced, &nreduced);
534: if (nfull == nreduced) { /* Also takes care of masked vectors */
535: if (mode == SCATTER_FORWARD) {
536: VecCopy(vreduced, vfull);
537: } else {
538: VecCopy(vfull, vreduced);
539: }
540: } else {
541: const PetscInt *id;
542: PetscInt i, n, m, rstart, rend;
544: ISGetIndices(is, &id);
545: ISGetLocalSize(is, &n);
546: VecGetLocalSize(vreduced, &m);
547: VecGetOwnershipRange(vfull, &rstart, &rend);
549: if (mode == SCATTER_FORWARD) {
550: PetscScalar *y;
551: const PetscScalar *x;
553: VecGetArray(vfull, &y);
554: VecGetArrayRead(vreduced, &x);
555: y -= rstart;
556: for (i = 0; i < n; ++i) {
557: if (id[i] < 0) continue;
559: y[id[i]] = x[i];
560: }
561: y += rstart;
562: VecRestoreArrayRead(vreduced, &x);
563: VecRestoreArray(vfull, &y);
564: } else if (mode == SCATTER_REVERSE) {
565: PetscScalar *x;
566: const PetscScalar *y;
568: VecGetArrayRead(vfull, &y);
569: VecGetArray(vreduced, &x);
570: for (i = 0; i < n; ++i) {
571: if (id[i] < 0) continue;
573: x[i] = y[id[i] - rstart];
574: }
575: VecRestoreArray(vreduced, &x);
576: VecRestoreArrayRead(vfull, &y);
577: } else SETERRQ(PetscObjectComm((PetscObject)vfull), PETSC_ERR_ARG_WRONG, "Only forward or reverse modes are legal");
578: ISRestoreIndices(is, &id);
579: }
580: return 0;
581: }
583: /*@
584: ISComplementVec - Creates the complement of the index set relative to a layout defined by a Vec
586: Collective on IS
588: Input Parameters:
589: + S - a PETSc IS
590: - V - the reference vector space
592: Output Parameter:
593: . T - the complement of S
595: Level: advanced
597: .seealso: `ISCreateGeneral()`
598: @*/
599: PetscErrorCode ISComplementVec(IS S, Vec V, IS *T)
600: {
601: PetscInt start, end;
603: VecGetOwnershipRange(V, &start, &end);
604: ISComplement(S, start, end, T);
605: return 0;
606: }
608: /*@
609: VecISSet - Sets the elements of a vector, specified by an index set, to a constant
611: Input Parameters:
612: + V - the vector
613: . S - index set for the locations in the vector
614: - c - the constant
616: Notes:
617: The index set identifies entries in the global vector.
618: Negative indices are skipped; indices outside the ownership range of V will raise an error.
620: Level: advanced
622: .seealso: `VecISCopy()`, `VecISAXPY()`, `VecSet()`
623: @*/
624: PetscErrorCode VecISSet(Vec V, IS S, PetscScalar c)
625: {
626: PetscInt nloc, low, high, i;
627: const PetscInt *s;
628: PetscScalar *v;
630: if (!S) return 0; /* simply return with no-op if the index set is NULL */
635: VecGetOwnershipRange(V, &low, &high);
636: ISGetLocalSize(S, &nloc);
637: ISGetIndices(S, &s);
638: VecGetArray(V, &v);
639: for (i = 0; i < nloc; ++i) {
640: if (s[i] < 0) continue;
642: v[s[i] - low] = c;
643: }
644: ISRestoreIndices(S, &s);
645: VecRestoreArray(V, &v);
646: return 0;
647: }
649: #if !defined(PETSC_USE_COMPLEX)
650: /*@C
651: VecBoundGradientProjection - Projects vector according to this definition.
652: If XL[i] < X[i] < XU[i], then GP[i] = G[i];
653: If X[i] <= XL[i], then GP[i] = min(G[i],0);
654: If X[i] >= XU[i], then GP[i] = max(G[i],0);
656: Input Parameters:
657: + G - current gradient vector
658: . X - current solution vector with XL[i] <= X[i] <= XU[i]
659: . XL - lower bounds
660: - XU - upper bounds
662: Output Parameter:
663: . GP - gradient projection vector
665: Notes:
666: GP may be the same vector as G
668: Level: advanced
669: @*/
670: PetscErrorCode VecBoundGradientProjection(Vec G, Vec X, Vec XL, Vec XU, Vec GP)
671: {
672: PetscInt n, i;
673: const PetscReal *xptr, *xlptr, *xuptr;
674: PetscReal *gptr, *gpptr;
675: PetscReal xval, gpval;
677: /* Project variables at the lower and upper bound */
683: if (!XL && !XU) {
684: VecCopy(G, GP);
685: return 0;
686: }
688: VecGetLocalSize(X, &n);
690: VecGetArrayRead(X, &xptr);
691: VecGetArrayRead(XL, &xlptr);
692: VecGetArrayRead(XU, &xuptr);
693: VecGetArrayPair(G, GP, &gptr, &gpptr);
695: for (i = 0; i < n; ++i) {
696: gpval = gptr[i];
697: xval = xptr[i];
698: if (gpval > 0.0 && xval <= xlptr[i]) {
699: gpval = 0.0;
700: } else if (gpval < 0.0 && xval >= xuptr[i]) {
701: gpval = 0.0;
702: }
703: gpptr[i] = gpval;
704: }
706: VecRestoreArrayRead(X, &xptr);
707: VecRestoreArrayRead(XL, &xlptr);
708: VecRestoreArrayRead(XU, &xuptr);
709: VecRestoreArrayPair(G, GP, &gptr, &gpptr);
710: return 0;
711: }
712: #endif
714: /*@
715: VecStepMaxBounded - See below
717: Collective on Vec
719: Input Parameters:
720: + X - vector with no negative entries
721: . XL - lower bounds
722: . XU - upper bounds
723: - DX - step direction, can have negative, positive or zero entries
725: Output Parameter:
726: . stepmax - minimum value so that X[i] + stepmax*DX[i] <= XL[i] or XU[i] <= X[i] + stepmax*DX[i]
728: Level: intermediate
730: @*/
731: PetscErrorCode VecStepMaxBounded(Vec X, Vec DX, Vec XL, Vec XU, PetscReal *stepmax)
732: {
733: PetscInt i, nn;
734: const PetscScalar *xx, *dx, *xl, *xu;
735: PetscReal localmax = 0;
742: VecGetArrayRead(X, &xx);
743: VecGetArrayRead(XL, &xl);
744: VecGetArrayRead(XU, &xu);
745: VecGetArrayRead(DX, &dx);
746: VecGetLocalSize(X, &nn);
747: for (i = 0; i < nn; i++) {
748: if (PetscRealPart(dx[i]) > 0) {
749: localmax = PetscMax(localmax, PetscRealPart((xu[i] - xx[i]) / dx[i]));
750: } else if (PetscRealPart(dx[i]) < 0) {
751: localmax = PetscMax(localmax, PetscRealPart((xl[i] - xx[i]) / dx[i]));
752: }
753: }
754: VecRestoreArrayRead(X, &xx);
755: VecRestoreArrayRead(XL, &xl);
756: VecRestoreArrayRead(XU, &xu);
757: VecRestoreArrayRead(DX, &dx);
758: MPIU_Allreduce(&localmax, stepmax, 1, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)X));
759: return 0;
760: }
762: /*@
763: VecStepBoundInfo - See below
765: Collective on Vec
767: Input Parameters:
768: + X - vector with no negative entries
769: . XL - lower bounds
770: . XU - upper bounds
771: - DX - step direction, can have negative, positive or zero entries
773: Output Parameters:
774: + boundmin - (may be NULL this it is not computed) maximum value so that XL[i] <= X[i] + boundmax*DX[i] <= XU[i]
775: . wolfemin - (may be NULL this it is not computed)
776: - 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]
778: Notes:
779: For complex numbers only compares the real part
781: Level: advanced
782: @*/
783: PetscErrorCode VecStepBoundInfo(Vec X, Vec DX, Vec XL, Vec XU, PetscReal *boundmin, PetscReal *wolfemin, PetscReal *boundmax)
784: {
785: PetscInt n, i;
786: const PetscScalar *x, *xl, *xu, *dx;
787: PetscReal t;
788: PetscReal localmin = PETSC_INFINITY, localwolfemin = PETSC_INFINITY, localmax = -1;
789: MPI_Comm comm;
796: VecGetArrayRead(X, &x);
797: VecGetArrayRead(XL, &xl);
798: VecGetArrayRead(XU, &xu);
799: VecGetArrayRead(DX, &dx);
800: VecGetLocalSize(X, &n);
801: for (i = 0; i < n; ++i) {
802: if (PetscRealPart(dx[i]) > 0 && PetscRealPart(xu[i]) < PETSC_INFINITY) {
803: t = PetscRealPart((xu[i] - x[i]) / dx[i]);
804: localmin = PetscMin(t, localmin);
805: if (localmin > 0) localwolfemin = PetscMin(t, localwolfemin);
806: localmax = PetscMax(t, localmax);
807: } else if (PetscRealPart(dx[i]) < 0 && PetscRealPart(xl[i]) > PETSC_NINFINITY) {
808: t = PetscRealPart((xl[i] - x[i]) / dx[i]);
809: localmin = PetscMin(t, localmin);
810: if (localmin > 0) localwolfemin = PetscMin(t, localwolfemin);
811: localmax = PetscMax(t, localmax);
812: }
813: }
815: VecRestoreArrayRead(X, &x);
816: VecRestoreArrayRead(XL, &xl);
817: VecRestoreArrayRead(XU, &xu);
818: VecRestoreArrayRead(DX, &dx);
819: PetscObjectGetComm((PetscObject)X, &comm);
821: if (boundmin) {
822: MPIU_Allreduce(&localmin, boundmin, 1, MPIU_REAL, MPIU_MIN, comm);
823: PetscInfo(X, "Step Bound Info: Closest Bound: %20.19e\n", (double)*boundmin);
824: }
825: if (wolfemin) {
826: MPIU_Allreduce(&localwolfemin, wolfemin, 1, MPIU_REAL, MPIU_MIN, comm);
827: PetscInfo(X, "Step Bound Info: Wolfe: %20.19e\n", (double)*wolfemin);
828: }
829: if (boundmax) {
830: MPIU_Allreduce(&localmax, boundmax, 1, MPIU_REAL, MPIU_MAX, comm);
831: if (*boundmax < 0) *boundmax = PETSC_INFINITY;
832: PetscInfo(X, "Step Bound Info: Max: %20.19e\n", (double)*boundmax);
833: }
834: return 0;
835: }
837: /*@
838: VecStepMax - Returns the largest value so that x[i] + step*DX[i] >= 0 for all i
840: Collective on Vec
842: Input Parameters:
843: + X - vector with no negative entries
844: - DX - a step direction, can have negative, positive or zero entries
846: Output Parameter:
847: . step - largest value such that x[i] + step*DX[i] >= 0 for all i
849: Notes:
850: For complex numbers only compares the real part
852: Level: advanced
853: @*/
854: PetscErrorCode VecStepMax(Vec X, Vec DX, PetscReal *step)
855: {
856: PetscInt i, nn;
857: PetscReal stepmax = PETSC_INFINITY;
858: const PetscScalar *xx, *dx;
863: VecGetLocalSize(X, &nn);
864: VecGetArrayRead(X, &xx);
865: VecGetArrayRead(DX, &dx);
866: for (i = 0; i < nn; ++i) {
868: if (PetscRealPart(dx[i]) < 0) stepmax = PetscMin(stepmax, PetscRealPart(-xx[i] / dx[i]));
869: }
870: VecRestoreArrayRead(X, &xx);
871: VecRestoreArrayRead(DX, &dx);
872: MPIU_Allreduce(&stepmax, step, 1, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)X));
873: return 0;
874: }
876: /*@
877: VecPow - Replaces each component of a vector by x_i^p
879: Logically Collective on v
881: Input Parameters:
882: + v - the vector
883: - p - the exponent to use on each element
885: Level: intermediate
887: @*/
888: PetscErrorCode VecPow(Vec v, PetscScalar p)
889: {
890: PetscInt n, i;
891: PetscScalar *v1;
895: VecGetArray(v, &v1);
896: VecGetLocalSize(v, &n);
898: if (1.0 == p) {
899: } else if (-1.0 == p) {
900: for (i = 0; i < n; ++i) v1[i] = 1.0 / v1[i];
901: } else if (0.0 == p) {
902: for (i = 0; i < n; ++i) {
903: /* Not-a-number left alone
904: Infinity set to one */
905: if (v1[i] == v1[i]) v1[i] = 1.0;
906: }
907: } else if (0.5 == p) {
908: for (i = 0; i < n; ++i) {
909: if (PetscRealPart(v1[i]) >= 0) {
910: v1[i] = PetscSqrtScalar(v1[i]);
911: } else {
912: v1[i] = PETSC_INFINITY;
913: }
914: }
915: } else if (-0.5 == p) {
916: for (i = 0; i < n; ++i) {
917: if (PetscRealPart(v1[i]) >= 0) {
918: v1[i] = 1.0 / PetscSqrtScalar(v1[i]);
919: } else {
920: v1[i] = PETSC_INFINITY;
921: }
922: }
923: } else if (2.0 == p) {
924: for (i = 0; i < n; ++i) v1[i] *= v1[i];
925: } else if (-2.0 == p) {
926: for (i = 0; i < n; ++i) v1[i] = 1.0 / (v1[i] * v1[i]);
927: } else {
928: for (i = 0; i < n; ++i) {
929: if (PetscRealPart(v1[i]) >= 0) {
930: v1[i] = PetscPowScalar(v1[i], p);
931: } else {
932: v1[i] = PETSC_INFINITY;
933: }
934: }
935: }
936: VecRestoreArray(v, &v1);
937: return 0;
938: }
940: /*@
941: VecMedian - Computes the componentwise median of three vectors
942: and stores the result in this vector. Used primarily for projecting
943: a vector within upper and lower bounds.
945: Logically Collective
947: Input Parameters:
948: + Vec1 - The first vector
949: . Vec2 - The second vector
950: - Vec3 - The third vector
952: Output Parameter:
953: . VMedian - The median vector (this can be any one of the input vectors)
955: Level: advanced
956: @*/
957: PetscErrorCode VecMedian(Vec Vec1, Vec Vec2, Vec Vec3, Vec VMedian)
958: {
959: PetscInt i, n, low1, high1;
960: const PetscScalar *v1, *v2, *v3;
961: PetscScalar *vmed;
968: if (!Vec1 && !Vec3) {
969: VecCopy(Vec2, VMedian);
970: return 0;
971: }
972: if (Vec1 == Vec2 || Vec1 == Vec3) {
973: VecCopy(Vec1, VMedian);
974: return 0;
975: }
976: if (Vec2 == Vec3) {
977: VecCopy(Vec2, VMedian);
978: return 0;
979: }
981: /* Assert that Vec1 != Vec2 and Vec2 != Vec3 */
992: VecCheckSameSize(Vec1, 1, Vec2, 2);
993: VecCheckSameSize(Vec1, 1, Vec3, 3);
994: VecCheckSameSize(Vec1, 1, VMedian, 4);
996: VecGetOwnershipRange(Vec1, &low1, &high1);
997: VecGetLocalSize(Vec1, &n);
998: if (n > 0) {
999: VecGetArray(VMedian, &vmed);
1000: if (Vec1 != VMedian) {
1001: VecGetArrayRead(Vec1, &v1);
1002: } else {
1003: v1 = vmed;
1004: }
1005: if (Vec2 != VMedian) {
1006: VecGetArrayRead(Vec2, &v2);
1007: } else {
1008: v2 = vmed;
1009: }
1010: if (Vec3 != VMedian) {
1011: VecGetArrayRead(Vec3, &v3);
1012: } else {
1013: v3 = vmed;
1014: }
1016: for (i = 0; i < n; ++i) 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])));
1018: VecRestoreArray(VMedian, &vmed);
1019: if (VMedian != Vec1) VecRestoreArrayRead(Vec1, &v1);
1020: if (VMedian != Vec2) VecRestoreArrayRead(Vec2, &v2);
1021: if (VMedian != Vec3) VecRestoreArrayRead(Vec3, &v3);
1022: }
1023: return 0;
1024: }