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]) {same[n_same]=low+i; ++n_same;}
47: }
49: if (Vec1 == Vec2) {
50: VecRestoreArrayRead(Vec1,&v1);
51: } else {
52: VecRestoreArrayRead(Vec1,&v1);
53: VecRestoreArrayRead(Vec2,&v2);
54: }
55: }
56: ISCreateGeneral(PetscObjectComm((PetscObject)Vec1),n_same,same,PETSC_OWN_POINTER,S);
57: return 0;
58: }
60: /*@
61: VecWhichLessThan - Creates an index set containing the indices
62: where the vectors Vec1 < Vec2
64: Collective on S
66: Input Parameters:
67: . Vec1, Vec2 - the two vectors to compare
69: OutputParameter:
70: . S - The index set containing the indices i where vec1[i] < vec2[i]
72: Notes:
73: The two vectors must have the same parallel layout
75: For complex numbers this only compares the real part
77: Level: advanced
78: @*/
79: PetscErrorCode VecWhichLessThan(Vec Vec1, Vec Vec2, IS *S)
80: {
81: PetscInt i,n_lt=0;
82: PetscInt n,low,high;
83: PetscInt *lt=NULL;
84: const PetscScalar *v1,*v2;
89: VecCheckSameSize(Vec1,1,Vec2,2);
91: VecGetOwnershipRange(Vec1,&low,&high);
92: VecGetLocalSize(Vec1,&n);
93: if (n>0) {
94: if (Vec1 == Vec2) {
95: VecGetArrayRead(Vec1,&v1);
96: v2=v1;
97: } else {
98: VecGetArrayRead(Vec1,&v1);
99: VecGetArrayRead(Vec2,&v2);
100: }
102: PetscMalloc1(n,<);
104: for (i=0; i<n; ++i) {
105: if (PetscRealPart(v1[i]) < PetscRealPart(v2[i])) {lt[n_lt]=low+i; ++n_lt;}
106: }
108: if (Vec1 == Vec2) {
109: VecRestoreArrayRead(Vec1,&v1);
110: } else {
111: VecRestoreArrayRead(Vec1,&v1);
112: VecRestoreArrayRead(Vec2,&v2);
113: }
114: }
115: ISCreateGeneral(PetscObjectComm((PetscObject)Vec1),n_lt,lt,PETSC_OWN_POINTER,S);
116: return 0;
117: }
119: /*@
120: VecWhichGreaterThan - Creates an index set containing the indices
121: where the vectors Vec1 > Vec2
123: Collective on S
125: Input Parameters:
126: . Vec1, Vec2 - the two vectors to compare
128: OutputParameter:
129: . S - The index set containing the indices i where vec1[i] > vec2[i]
131: Notes:
132: The two vectors must have the same parallel layout
134: For complex numbers this only compares the real part
136: Level: advanced
137: @*/
138: PetscErrorCode VecWhichGreaterThan(Vec Vec1, Vec Vec2, IS *S)
139: {
140: PetscInt i,n_gt=0;
141: PetscInt n,low,high;
142: PetscInt *gt=NULL;
143: const PetscScalar *v1,*v2;
148: VecCheckSameSize(Vec1,1,Vec2,2);
150: VecGetOwnershipRange(Vec1,&low,&high);
151: VecGetLocalSize(Vec1,&n);
152: if (n>0) {
153: if (Vec1 == Vec2) {
154: VecGetArrayRead(Vec1,&v1);
155: v2=v1;
156: } else {
157: VecGetArrayRead(Vec1,&v1);
158: VecGetArrayRead(Vec2,&v2);
159: }
161: PetscMalloc1(n,>);
163: for (i=0; i<n; ++i) {
164: if (PetscRealPart(v1[i]) > PetscRealPart(v2[i])) {gt[n_gt]=low+i; ++n_gt;}
165: }
167: if (Vec1 == Vec2) {
168: VecRestoreArrayRead(Vec1,&v1);
169: } else {
170: VecRestoreArrayRead(Vec1,&v1);
171: VecRestoreArrayRead(Vec2,&v2);
172: }
173: }
174: ISCreateGeneral(PetscObjectComm((PetscObject)Vec1),n_gt,gt,PETSC_OWN_POINTER,S);
175: return 0;
176: }
178: /*@
179: VecWhichBetween - Creates an index set containing the indices
180: where VecLow < V < VecHigh
182: Collective on S
184: Input Parameters:
185: + VecLow - lower bound
186: . V - Vector to compare
187: - VecHigh - higher bound
189: OutputParameter:
190: . S - The index set containing the indices i where veclow[i] < v[i] < vechigh[i]
192: Notes:
193: The vectors must have the same parallel layout
195: For complex numbers this only compares the real part
197: Level: advanced
198: @*/
199: PetscErrorCode VecWhichBetween(Vec VecLow, Vec V, Vec VecHigh, IS *S)
200: {
202: PetscInt i,n_vm=0;
203: PetscInt n,low,high;
204: PetscInt *vm=NULL;
205: const PetscScalar *v1,*v2,*vmiddle;
212: VecCheckSameSize(V,2,VecLow,1);
213: VecCheckSameSize(V,2,VecHigh,3);
215: VecGetOwnershipRange(VecLow,&low,&high);
216: VecGetLocalSize(VecLow,&n);
217: if (n>0) {
218: VecGetArrayRead(VecLow,&v1);
219: if (VecLow != VecHigh) {
220: VecGetArrayRead(VecHigh,&v2);
221: } else {
222: v2=v1;
223: }
224: if (V != VecLow && V != VecHigh) {
225: VecGetArrayRead(V,&vmiddle);
226: } else if (V==VecLow) {
227: vmiddle=v1;
228: } else {
229: vmiddle=v2;
230: }
232: PetscMalloc1(n,&vm);
234: for (i=0; i<n; ++i) {
235: if (PetscRealPart(v1[i]) < PetscRealPart(vmiddle[i]) && PetscRealPart(vmiddle[i]) < PetscRealPart(v2[i])) {vm[n_vm]=low+i; ++n_vm;}
236: }
238: VecRestoreArrayRead(VecLow,&v1);
239: if (VecLow != VecHigh) {
240: VecRestoreArrayRead(VecHigh,&v2);
241: }
242: if (V != VecLow && V != VecHigh) {
243: VecRestoreArrayRead(V,&vmiddle);
244: }
245: }
246: ISCreateGeneral(PetscObjectComm((PetscObject)V),n_vm,vm,PETSC_OWN_POINTER,S);
247: return 0;
248: }
250: /*@
251: VecWhichBetweenOrEqual - Creates an index set containing the indices
252: where VecLow <= V <= VecHigh
254: Collective on S
256: Input Parameters:
257: + VecLow - lower bound
258: . V - Vector to compare
259: - VecHigh - higher bound
261: OutputParameter:
262: . S - The index set containing the indices i where veclow[i] <= v[i] <= vechigh[i]
264: Level: advanced
265: @*/
267: PetscErrorCode VecWhichBetweenOrEqual(Vec VecLow, Vec V, Vec VecHigh, IS * S)
268: {
269: PetscInt i,n_vm=0;
270: PetscInt n,low,high;
271: PetscInt *vm=NULL;
272: const PetscScalar *v1,*v2,*vmiddle;
279: VecCheckSameSize(V,2,VecLow,1);
280: VecCheckSameSize(V,2,VecHigh,3);
282: VecGetOwnershipRange(VecLow,&low,&high);
283: VecGetLocalSize(VecLow,&n);
284: if (n>0) {
285: VecGetArrayRead(VecLow,&v1);
286: if (VecLow != VecHigh) {
287: VecGetArrayRead(VecHigh,&v2);
288: } else {
289: v2=v1;
290: }
291: if (V != VecLow && V != VecHigh) {
292: VecGetArrayRead(V,&vmiddle);
293: } else if (V==VecLow) {
294: vmiddle=v1;
295: } else {
296: vmiddle =v2;
297: }
299: PetscMalloc1(n,&vm);
301: for (i=0; i<n; ++i) {
302: if (PetscRealPart(v1[i]) <= PetscRealPart(vmiddle[i]) && PetscRealPart(vmiddle[i]) <= PetscRealPart(v2[i])) {vm[n_vm]=low+i; ++n_vm;}
303: }
305: VecRestoreArrayRead(VecLow,&v1);
306: if (VecLow != VecHigh) {
307: VecRestoreArrayRead(VecHigh,&v2);
308: }
309: if (V != VecLow && V != VecHigh) {
310: VecRestoreArrayRead(V,&vmiddle);
311: }
312: }
313: ISCreateGeneral(PetscObjectComm((PetscObject)V),n_vm,vm,PETSC_OWN_POINTER,S);
314: return 0;
315: }
317: /*@
318: VecWhichInactive - Creates an index set containing the indices
319: where one of the following holds:
320: a) VecLow(i) < V(i) < VecHigh(i)
321: b) VecLow(i) = V(i) and D(i) <= 0 (< 0 when Strong is true)
322: c) VecHigh(i) = V(i) and D(i) >= 0 (> 0 when Strong is true)
324: Collective on S
326: Input Parameters:
327: + VecLow - lower bound
328: . V - Vector to compare
329: . D - Direction to compare
330: . VecHigh - higher bound
331: - Strong - indicator for applying strongly inactive test
333: OutputParameter:
334: . S - The index set containing the indices i where the bound is inactive
336: Level: advanced
337: @*/
339: PetscErrorCode VecWhichInactive(Vec VecLow, Vec V, Vec D, Vec VecHigh, PetscBool Strong, IS * S)
340: {
341: PetscInt i,n_vm=0;
342: PetscInt n,low,high;
343: PetscInt *vm=NULL;
344: const PetscScalar *v1,*v2,*v,*d;
353: VecCheckSameSize(V,2,VecLow,1);
354: VecCheckSameSize(V,2,D,3);
355: VecCheckSameSize(V,2,VecHigh,4);
357: VecGetOwnershipRange(VecLow,&low,&high);
358: VecGetLocalSize(VecLow,&n);
359: if (n>0) {
360: VecGetArrayRead(VecLow,&v1);
361: if (VecLow != VecHigh) {
362: VecGetArrayRead(VecHigh,&v2);
363: } else {
364: v2=v1;
365: }
366: if (V != VecLow && V != VecHigh) {
367: VecGetArrayRead(V,&v);
368: } else if (V==VecLow) {
369: v = v1;
370: } else {
371: v = v2;
372: }
373: if (D != VecLow && D != VecHigh && D != V) {
374: VecGetArrayRead(D,&d);
375: } else if (D==VecLow) {
376: d = v1;
377: } else if (D==VecHigh) {
378: d = v2;
379: } else {
380: d = v;
381: }
383: PetscMalloc1(n,&vm);
385: if (Strong) {
386: for (i=0; i<n; ++i) {
387: if (PetscRealPart(v1[i]) < PetscRealPart(v[i]) && PetscRealPart(v[i]) < PetscRealPart(v2[i])) {
388: vm[n_vm]=low+i; ++n_vm;
389: } else if (PetscRealPart(v1[i]) == PetscRealPart(v[i]) && PetscRealPart(d[i]) < 0) {
390: vm[n_vm]=low+i; ++n_vm;
391: } else if (PetscRealPart(v2[i]) == PetscRealPart(v[i]) && PetscRealPart(d[i]) > 0) {
392: vm[n_vm]=low+i; ++n_vm;
393: }
394: }
395: } else {
396: for (i=0; i<n; ++i) {
397: if (PetscRealPart(v1[i]) < PetscRealPart(v[i]) && PetscRealPart(v[i]) < PetscRealPart(v2[i])) {
398: vm[n_vm]=low+i; ++n_vm;
399: } else if (PetscRealPart(v1[i]) == PetscRealPart(v[i]) && PetscRealPart(d[i]) <= 0) {
400: vm[n_vm]=low+i; ++n_vm;
401: } else if (PetscRealPart(v2[i]) == PetscRealPart(v[i]) && PetscRealPart(d[i]) >= 0) {
402: vm[n_vm]=low+i; ++n_vm;
403: }
404: }
405: }
407: VecRestoreArrayRead(VecLow,&v1);
408: if (VecLow != VecHigh) {
409: VecRestoreArrayRead(VecHigh,&v2);
410: }
411: if (V != VecLow && V != VecHigh) {
412: VecRestoreArrayRead(V,&v);
413: }
414: if (D != VecLow && D != VecHigh && D != V) {
415: VecRestoreArrayRead(D,&d);
416: }
417: }
418: ISCreateGeneral(PetscObjectComm((PetscObject)V),n_vm,vm,PETSC_OWN_POINTER,S);
419: return 0;
420: }
422: /*@
423: VecISAXPY - Adds a reduced vector to the appropriate elements of a full-space vector.
424: vfull[is[i]] += alpha*vreduced[i]
426: Input Parameters:
427: + vfull - the full-space vector
428: . is - the index set for the reduced space
429: . alpha - the scalar coefficient
430: - vreduced - the reduced-space vector
432: Output Parameters:
433: . vfull - the sum of the full-space vector and reduced-space vector
435: Notes:
436: The index set identifies entries in the global vector.
437: Negative indices are skipped; indices outside the ownership range of vfull will raise an error.
439: Level: advanced
441: .seealso: VecAXPY(), VecGetOwnershipRange()
442: @*/
443: PetscErrorCode VecISAXPY(Vec vfull, IS is, PetscScalar alpha, Vec vreduced)
444: {
445: PetscInt nfull,nreduced;
450: VecGetSize(vfull,&nfull);
451: VecGetSize(vreduced,&nreduced);
453: if (nfull == nreduced) { /* Also takes care of masked vectors */
454: VecAXPY(vfull,alpha,vreduced);
455: } else {
456: PetscScalar *y;
457: const PetscScalar *x;
458: PetscInt i,n,m,rstart,rend;
459: const PetscInt *id;
461: VecGetArray(vfull,&y);
462: VecGetArrayRead(vreduced,&x);
463: ISGetIndices(is,&id);
464: ISGetLocalSize(is,&n);
465: VecGetLocalSize(vreduced,&m);
467: VecGetOwnershipRange(vfull,&rstart,&rend);
468: y -= rstart;
469: if (alpha == 1.0) {
470: for (i=0; i<n; ++i) {
471: if (id[i] < 0) continue;
473: y[id[i]] += x[i];
474: }
475: } else {
476: for (i=0; i<n; ++i) {
477: if (id[i] < 0) continue;
479: y[id[i]] += alpha*x[i];
480: }
481: }
482: y += rstart;
483: ISRestoreIndices(is,&id);
484: VecRestoreArray(vfull,&y);
485: VecRestoreArrayRead(vreduced,&x);
486: }
487: return 0;
488: }
490: /*@
491: VecISCopy - Copies between a reduced vector and the appropriate elements of a full-space vector.
493: Input Parameters:
494: + vfull - the full-space vector
495: . is - the index set for the reduced space
496: . mode - the direction of copying, SCATTER_FORWARD or SCATTER_REVERSE
497: - vreduced - the reduced-space vector
499: Output Parameters:
500: . vfull - the sum of the full-space vector and reduced-space vector
502: Notes:
503: The index set identifies entries in the global vector.
504: Negative indices are skipped; indices outside the ownership range of vfull will raise an error.
506: mode == SCATTER_FORWARD: vfull[is[i]] = vreduced[i]
507: mode == SCATTER_REVERSE: vreduced[i] = vfull[is[i]]
509: Level: advanced
511: .seealso: VecAXPY(), VecGetOwnershipRange()
512: @*/
513: PetscErrorCode VecISCopy(Vec vfull, IS is, ScatterMode mode, Vec vreduced)
514: {
515: PetscInt nfull, nreduced;
520: VecGetSize(vfull, &nfull);
521: VecGetSize(vreduced, &nreduced);
523: if (nfull == nreduced) { /* Also takes care of masked vectors */
524: if (mode == SCATTER_FORWARD) {
525: VecCopy(vreduced, vfull);
526: } else {
527: VecCopy(vfull, vreduced);
528: }
529: } else {
530: const PetscInt *id;
531: PetscInt i, n, m, rstart, rend;
533: ISGetIndices(is, &id);
534: ISGetLocalSize(is, &n);
535: VecGetLocalSize(vreduced, &m);
536: VecGetOwnershipRange(vfull, &rstart, &rend);
538: if (mode == SCATTER_FORWARD) {
539: PetscScalar *y;
540: const PetscScalar *x;
542: VecGetArray(vfull, &y);
543: VecGetArrayRead(vreduced, &x);
544: y -= rstart;
545: for (i = 0; i < n; ++i) {
546: if (id[i] < 0) continue;
548: y[id[i]] = x[i];
549: }
550: y += rstart;
551: VecRestoreArrayRead(vreduced, &x);
552: VecRestoreArray(vfull, &y);
553: } else if (mode == SCATTER_REVERSE) {
554: PetscScalar *x;
555: const PetscScalar *y;
557: VecGetArrayRead(vfull, &y);
558: VecGetArray(vreduced, &x);
559: for (i = 0; i < n; ++i) {
560: if (id[i] < 0) continue;
562: x[i] = y[id[i]-rstart];
563: }
564: VecRestoreArray(vreduced, &x);
565: VecRestoreArrayRead(vfull, &y);
566: } else SETERRQ(PetscObjectComm((PetscObject) vfull), PETSC_ERR_ARG_WRONG, "Only forward or reverse modes are legal");
567: ISRestoreIndices(is, &id);
568: }
569: return 0;
570: }
572: /*@
573: ISComplementVec - Creates the complement of the index set relative to a layout defined by a Vec
575: Collective on IS
577: Input Parameters:
578: + S - a PETSc IS
579: - V - the reference vector space
581: Output Parameter:
582: . T - the complement of S
584: Level: advanced
586: .seealso: ISCreateGeneral()
587: @*/
588: PetscErrorCode ISComplementVec(IS S, Vec V, IS *T)
589: {
590: PetscInt start, end;
592: VecGetOwnershipRange(V,&start,&end);
593: ISComplement(S,start,end,T);
594: return 0;
595: }
597: /*@
598: VecISSet - Sets the elements of a vector, specified by an index set, to a constant
600: Input Parameters:
601: + V - the vector
602: . S - index set for the locations in the vector
603: - c - the constant
605: Notes:
606: The index set identifies entries in the global vector.
607: Negative indices are skipped; indices outside the ownership range of V will raise an error.
609: Level: advanced
611: .seealso: VecSet(), VecGetOwnershipRange()
612: @*/
613: PetscErrorCode VecISSet(Vec V,IS S, PetscScalar c)
614: {
615: PetscInt nloc,low,high,i;
616: const PetscInt *s;
617: PetscScalar *v;
619: if (!S) return 0; /* simply return with no-op if the index set is NULL */
624: VecGetOwnershipRange(V,&low,&high);
625: ISGetLocalSize(S,&nloc);
626: ISGetIndices(S,&s);
627: VecGetArray(V,&v);
628: for (i=0; i<nloc; ++i) {
629: if (s[i] < 0) continue;
631: v[s[i]-low] = c;
632: }
633: ISRestoreIndices(S,&s);
634: VecRestoreArray(V,&v);
635: return 0;
636: }
638: #if !defined(PETSC_USE_COMPLEX)
639: /*@C
640: VecBoundGradientProjection - Projects vector according to this definition.
641: If XL[i] < X[i] < XU[i], then GP[i] = G[i];
642: If X[i] <= XL[i], then GP[i] = min(G[i],0);
643: If X[i] >= XU[i], then GP[i] = max(G[i],0);
645: Input Parameters:
646: + G - current gradient vector
647: . X - current solution vector with XL[i] <= X[i] <= XU[i]
648: . XL - lower bounds
649: - XU - upper bounds
651: Output Parameter:
652: . GP - gradient projection vector
654: Notes:
655: GP may be the same vector as G
657: Level: advanced
658: @*/
659: PetscErrorCode VecBoundGradientProjection(Vec G, Vec X, Vec XL, Vec XU, Vec GP)
660: {
662: PetscInt n,i;
663: const PetscReal *xptr,*xlptr,*xuptr;
664: PetscReal *gptr,*gpptr;
665: PetscReal xval,gpval;
667: /* Project variables at the lower and upper bound */
674: VecGetLocalSize(X,&n);
676: VecGetArrayRead(X,&xptr);
677: VecGetArrayRead(XL,&xlptr);
678: VecGetArrayRead(XU,&xuptr);
679: VecGetArrayPair(G,GP,&gptr,&gpptr);
681: for (i=0; i<n; ++i) {
682: gpval = gptr[i]; xval = xptr[i];
683: if (gpval>0.0 && xval<=xlptr[i]) {
684: gpval = 0.0;
685: } else if (gpval<0.0 && xval>=xuptr[i]) {
686: gpval = 0.0;
687: }
688: gpptr[i] = gpval;
689: }
691: VecRestoreArrayRead(X,&xptr);
692: VecRestoreArrayRead(XL,&xlptr);
693: VecRestoreArrayRead(XU,&xuptr);
694: VecRestoreArrayPair(G,GP,&gptr,&gpptr);
695: return 0;
696: }
697: #endif
699: /*@
700: VecStepMaxBounded - See below
702: Collective on Vec
704: Input Parameters:
705: + X - vector with no negative entries
706: . XL - lower bounds
707: . XU - upper bounds
708: - DX - step direction, can have negative, positive or zero entries
710: Output Parameter:
711: . stepmax - minimum value so that X[i] + stepmax*DX[i] <= XL[i] or XU[i] <= X[i] + stepmax*DX[i]
713: Level: intermediate
715: @*/
716: PetscErrorCode VecStepMaxBounded(Vec X, Vec DX, Vec XL, Vec XU, PetscReal *stepmax)
717: {
718: PetscInt i,nn;
719: const PetscScalar *xx,*dx,*xl,*xu;
720: PetscReal localmax=0;
727: VecGetArrayRead(X,&xx);
728: VecGetArrayRead(XL,&xl);
729: VecGetArrayRead(XU,&xu);
730: VecGetArrayRead(DX,&dx);
731: VecGetLocalSize(X,&nn);
732: for (i=0;i<nn;i++) {
733: if (PetscRealPart(dx[i]) > 0) {
734: localmax=PetscMax(localmax,PetscRealPart((xu[i]-xx[i])/dx[i]));
735: } else if (PetscRealPart(dx[i])<0) {
736: localmax=PetscMax(localmax,PetscRealPart((xl[i]-xx[i])/dx[i]));
737: }
738: }
739: VecRestoreArrayRead(X,&xx);
740: VecRestoreArrayRead(XL,&xl);
741: VecRestoreArrayRead(XU,&xu);
742: VecRestoreArrayRead(DX,&dx);
743: MPIU_Allreduce(&localmax,stepmax,1,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)X));
744: return 0;
745: }
747: /*@
748: VecStepBoundInfo - See below
750: Collective on Vec
752: Input Parameters:
753: + X - vector with no negative entries
754: . XL - lower bounds
755: . XU - upper bounds
756: - DX - step direction, can have negative, positive or zero entries
758: Output Parameters:
759: + boundmin - (may be NULL this it is not computed) maximum value so that XL[i] <= X[i] + boundmax*DX[i] <= XU[i]
760: . wolfemin - (may be NULL this it is not computed)
761: - 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]
763: Notes:
764: For complex numbers only compares the real part
766: Level: advanced
767: @*/
768: PetscErrorCode VecStepBoundInfo(Vec X, Vec DX, Vec XL, Vec XU, PetscReal *boundmin, PetscReal *wolfemin, PetscReal *boundmax)
769: {
770: PetscInt n,i;
771: const PetscScalar *x,*xl,*xu,*dx;
772: PetscReal t;
773: PetscReal localmin=PETSC_INFINITY,localwolfemin=PETSC_INFINITY,localmax=-1;
774: MPI_Comm comm;
781: VecGetArrayRead(X,&x);
782: VecGetArrayRead(XL,&xl);
783: VecGetArrayRead(XU,&xu);
784: VecGetArrayRead(DX,&dx);
785: VecGetLocalSize(X,&n);
786: for (i=0; i<n; ++i) {
787: if (PetscRealPart(dx[i])>0 && PetscRealPart(xu[i]) < PETSC_INFINITY) {
788: t=PetscRealPart((xu[i]-x[i])/dx[i]);
789: localmin=PetscMin(t,localmin);
790: if (localmin>0) {
791: localwolfemin = PetscMin(t,localwolfemin);
792: }
793: localmax = PetscMax(t,localmax);
794: } else if (PetscRealPart(dx[i])<0 && PetscRealPart(xl[i]) > PETSC_NINFINITY) {
795: t=PetscRealPart((xl[i]-x[i])/dx[i]);
796: localmin = PetscMin(t,localmin);
797: if (localmin>0) {
798: localwolfemin = PetscMin(t,localwolfemin);
799: }
800: localmax = PetscMax(t,localmax);
801: }
802: }
804: VecRestoreArrayRead(X,&x);
805: VecRestoreArrayRead(XL,&xl);
806: VecRestoreArrayRead(XU,&xu);
807: VecRestoreArrayRead(DX,&dx);
808: PetscObjectGetComm((PetscObject)X,&comm);
810: if (boundmin) {
811: MPIU_Allreduce(&localmin,boundmin,1,MPIU_REAL,MPIU_MIN,comm);
812: PetscInfo(X,"Step Bound Info: Closest Bound: %20.19e\n",(double)*boundmin);
813: }
814: if (wolfemin) {
815: MPIU_Allreduce(&localwolfemin,wolfemin,1,MPIU_REAL,MPIU_MIN,comm);
816: PetscInfo(X,"Step Bound Info: Wolfe: %20.19e\n",(double)*wolfemin);
817: }
818: if (boundmax) {
819: MPIU_Allreduce(&localmax,boundmax,1,MPIU_REAL,MPIU_MAX,comm);
820: if (*boundmax < 0) *boundmax=PETSC_INFINITY;
821: PetscInfo(X,"Step Bound Info: Max: %20.19e\n",(double)*boundmax);
822: }
823: return 0;
824: }
826: /*@
827: VecStepMax - Returns the largest value so that x[i] + step*DX[i] >= 0 for all i
829: Collective on Vec
831: Input Parameters:
832: + X - vector with no negative entries
833: - DX - a step direction, can have negative, positive or zero entries
835: Output Parameter:
836: . step - largest value such that x[i] + step*DX[i] >= 0 for all i
838: Notes:
839: For complex numbers only compares the real part
841: Level: advanced
842: @*/
843: PetscErrorCode VecStepMax(Vec X, Vec DX, PetscReal *step)
844: {
845: PetscInt i,nn;
846: PetscReal stepmax=PETSC_INFINITY;
847: const PetscScalar *xx,*dx;
852: VecGetLocalSize(X,&nn);
853: VecGetArrayRead(X,&xx);
854: VecGetArrayRead(DX,&dx);
855: for (i=0;i<nn;++i) {
857: else if (PetscRealPart(dx[i])<0) stepmax=PetscMin(stepmax,PetscRealPart(-xx[i]/dx[i]));
858: }
859: VecRestoreArrayRead(X,&xx);
860: VecRestoreArrayRead(DX,&dx);
861: MPIU_Allreduce(&stepmax,step,1,MPIU_REAL,MPIU_MIN,PetscObjectComm((PetscObject)X));
862: return 0;
863: }
865: /*@
866: VecPow - Replaces each component of a vector by x_i^p
868: Logically Collective on v
870: Input Parameters:
871: + v - the vector
872: - p - the exponent to use on each element
874: Level: intermediate
876: @*/
877: PetscErrorCode VecPow(Vec v, PetscScalar p)
878: {
879: PetscInt n,i;
880: PetscScalar *v1;
884: VecGetArray(v,&v1);
885: VecGetLocalSize(v,&n);
887: if (1.0 == p) {
888: } else if (-1.0 == p) {
889: for (i = 0; i < n; ++i) {
890: v1[i] = 1.0 / v1[i];
891: }
892: } else if (0.0 == p) {
893: for (i = 0; i < n; ++i) {
894: /* Not-a-number left alone
895: Infinity set to one */
896: if (v1[i] == v1[i]) {
897: v1[i] = 1.0;
898: }
899: }
900: } else if (0.5 == p) {
901: for (i = 0; i < n; ++i) {
902: if (PetscRealPart(v1[i]) >= 0) {
903: v1[i] = PetscSqrtScalar(v1[i]);
904: } else {
905: v1[i] = PETSC_INFINITY;
906: }
907: }
908: } else if (-0.5 == p) {
909: for (i = 0; i < n; ++i) {
910: if (PetscRealPart(v1[i]) >= 0) {
911: v1[i] = 1.0 / PetscSqrtScalar(v1[i]);
912: } else {
913: v1[i] = PETSC_INFINITY;
914: }
915: }
916: } else if (2.0 == p) {
917: for (i = 0; i < n; ++i) {
918: v1[i] *= v1[i];
919: }
920: } else if (-2.0 == p) {
921: for (i = 0; i < n; ++i) {
922: v1[i] = 1.0 / (v1[i] * v1[i]);
923: }
924: } else {
925: for (i = 0; i < n; ++i) {
926: if (PetscRealPart(v1[i]) >= 0) {
927: v1[i] = PetscPowScalar(v1[i],p);
928: } else {
929: v1[i] = PETSC_INFINITY;
930: }
931: }
932: }
933: VecRestoreArray(v,&v1);
934: return 0;
935: }
937: /*@
938: VecMedian - Computes the componentwise median of three vectors
939: and stores the result in this vector. Used primarily for projecting
940: a vector within upper and lower bounds.
942: Logically Collective
944: Input Parameters:
945: + Vec1 - The first vector
946: . Vec2 - The second vector
947: - Vec3 - The third vector
949: Output Parameter:
950: . VMedian - The median vector (this can be any one of the input vectors)
952: Level: advanced
953: @*/
954: PetscErrorCode VecMedian(Vec Vec1, Vec Vec2, Vec Vec3, Vec VMedian)
955: {
956: PetscInt i,n,low1,high1;
957: const PetscScalar *v1,*v2,*v3;
958: PetscScalar *vmed;
965: if (Vec1==Vec2 || Vec1==Vec3) {
966: VecCopy(Vec1,VMedian);
967: return 0;
968: }
969: if (Vec2==Vec3) {
970: VecCopy(Vec2,VMedian);
971: return 0;
972: }
974: /* Assert that Vec1 != Vec2 and Vec2 != Vec3 */
985: VecCheckSameSize(Vec1,1,Vec2,2);
986: VecCheckSameSize(Vec1,1,Vec3,3);
987: VecCheckSameSize(Vec1,1,VMedian,4);
989: VecGetOwnershipRange(Vec1,&low1,&high1);
990: VecGetLocalSize(Vec1,&n);
991: if (n>0) {
992: VecGetArray(VMedian,&vmed);
993: if (Vec1 != VMedian) {
994: VecGetArrayRead(Vec1,&v1);
995: } else {
996: v1=vmed;
997: }
998: if (Vec2 != VMedian) {
999: VecGetArrayRead(Vec2,&v2);
1000: } else {
1001: v2=vmed;
1002: }
1003: if (Vec3 != VMedian) {
1004: VecGetArrayRead(Vec3,&v3);
1005: } else {
1006: v3=vmed;
1007: }
1009: for (i=0;i<n;++i) {
1010: 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])));
1011: }
1013: VecRestoreArray(VMedian,&vmed);
1014: if (VMedian != Vec1) {
1015: VecRestoreArrayRead(Vec1,&v1);
1016: }
1017: if (VMedian != Vec2) {
1018: VecRestoreArrayRead(Vec2,&v2);
1019: }
1020: if (VMedian != Vec3) {
1021: VecRestoreArrayRead(Vec3,&v3);
1022: }
1023: }
1024: return 0;
1025: }