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: }
259: /*@
260: VecWhichBetweenOrEqual - Creates an index set containing the indices
261: where VecLow <= V <= VecHigh
263: Collective on S
265: Input Parameters:
266: + VecLow - lower bound
267: . V - Vector to compare
268: - VecHigh - higher bound
270: OutputParameter:
271: . S - The index set containing the indices i where veclow[i] <= v[i] <= vechigh[i]
273: Level: advanced
274: @*/
276: PetscErrorCode VecWhichBetweenOrEqual(Vec VecLow, Vec V, Vec VecHigh, IS * S)
277: {
278: PetscErrorCode ierr;
279: PetscInt i,n_vm=0;
280: PetscInt n,low,high;
281: PetscInt *vm=NULL;
282: const PetscScalar *v1,*v2,*vmiddle;
290: VecCheckSameSize(V,2,VecLow,1);
291: VecCheckSameSize(V,2,VecHigh,3);
293: VecGetOwnershipRange(VecLow,&low,&high);
294: VecGetLocalSize(VecLow,&n);
295: if (n>0){
296: VecGetArrayRead(VecLow,&v1);
297: if (VecLow != VecHigh){
298: VecGetArrayRead(VecHigh,&v2);
299: } else {
300: v2=v1;
301: }
302: if (V != VecLow && V != VecHigh){
303: VecGetArrayRead(V,&vmiddle);
304: } else if (V==VecLow){
305: vmiddle=v1;
306: } else {
307: vmiddle =v2;
308: }
310: PetscMalloc1(n,&vm);
312: for (i=0; i<n; ++i){
313: if (PetscRealPart(v1[i]) <= PetscRealPart(vmiddle[i]) && PetscRealPart(vmiddle[i]) <= PetscRealPart(v2[i])) {vm[n_vm]=low+i; ++n_vm;}
314: }
316: VecRestoreArrayRead(VecLow,&v1);
317: if (VecLow != VecHigh){
318: VecRestoreArrayRead(VecHigh,&v2);
319: }
320: if (V != VecLow && V != VecHigh){
321: VecRestoreArrayRead(V,&vmiddle);
322: }
323: }
324: ISCreateGeneral(PetscObjectComm((PetscObject)V),n_vm,vm,PETSC_OWN_POINTER,S);
325: return(0);
326: }
328: /*@
329: VecWhichInactive - Creates an index set containing the indices
330: where one of the following holds:
331: a) VecLow(i) < V(i) < VecHigh(i)
332: b) VecLow(i) = V(i) and D(i) <= 0 (< 0 when Strong is true)
333: c) VecHigh(i) = V(i) and D(i) >= 0 (> 0 when Strong is true)
335: Collective on S
337: Input Parameters:
338: + VecLow - lower bound
339: . V - Vector to compare
340: . D - Direction to compare
341: . VecHigh - higher bound
342: - Strong - indicator for applying strongly inactive test
344: OutputParameter:
345: . S - The index set containing the indices i where the bound is inactive
347: Level: advanced
348: @*/
350: PetscErrorCode VecWhichInactive(Vec VecLow, Vec V, Vec D, Vec VecHigh, PetscBool Strong, IS * S)
351: {
352: PetscErrorCode ierr;
353: PetscInt i,n_vm=0;
354: PetscInt n,low,high;
355: PetscInt *vm=NULL;
356: const PetscScalar *v1,*v2,*v,*d;
366: VecCheckSameSize(V,2,VecLow,1);
367: VecCheckSameSize(V,2,D,3);
368: VecCheckSameSize(V,2,VecHigh,4);
370: VecGetOwnershipRange(VecLow,&low,&high);
371: VecGetLocalSize(VecLow,&n);
372: if (n>0){
373: VecGetArrayRead(VecLow,&v1);
374: if (VecLow != VecHigh){
375: VecGetArrayRead(VecHigh,&v2);
376: } else {
377: v2=v1;
378: }
379: if (V != VecLow && V != VecHigh){
380: VecGetArrayRead(V,&v);
381: } else if (V==VecLow){
382: v = v1;
383: } else {
384: v = v2;
385: }
386: if (D != VecLow && D != VecHigh && D != V){
387: VecGetArrayRead(D,&d);
388: } else if (D==VecLow){
389: d = v1;
390: } else if (D==VecHigh){
391: d = v2;
392: } else {
393: d = v;
394: }
396: PetscMalloc1(n,&vm);
398: if (Strong){
399: for (i=0; i<n; ++i) {
400: if (PetscRealPart(v1[i]) < PetscRealPart(v[i]) && PetscRealPart(v[i]) < PetscRealPart(v2[i])){
401: vm[n_vm]=low+i; ++n_vm;
402: } else if (PetscRealPart(v1[i]) == PetscRealPart(v[i]) && PetscRealPart(d[i]) < 0){
403: vm[n_vm]=low+i; ++n_vm;
404: } else if (PetscRealPart(v2[i]) == PetscRealPart(v[i]) && PetscRealPart(d[i]) > 0){
405: vm[n_vm]=low+i; ++n_vm;
406: }
407: }
408: } else {
409: for (i=0; i<n; ++i) {
410: if (PetscRealPart(v1[i]) < PetscRealPart(v[i]) && PetscRealPart(v[i]) < PetscRealPart(v2[i])){
411: vm[n_vm]=low+i; ++n_vm;
412: } else if (PetscRealPart(v1[i]) == PetscRealPart(v[i]) && PetscRealPart(d[i]) <= 0){
413: vm[n_vm]=low+i; ++n_vm;
414: } else if (PetscRealPart(v2[i]) == PetscRealPart(v[i]) && PetscRealPart(d[i]) >= 0){
415: vm[n_vm]=low+i; ++n_vm;
416: }
417: }
418: }
420: VecRestoreArrayRead(VecLow,&v1);
421: if (VecLow != VecHigh){
422: VecRestoreArrayRead(VecHigh,&v2);
423: }
424: if (V != VecLow && V != VecHigh){
425: VecRestoreArrayRead(V,&v);
426: }
427: if (D != VecLow && D != VecHigh && D != V){
428: VecRestoreArrayRead(D,&d);
429: }
430: }
431: ISCreateGeneral(PetscObjectComm((PetscObject)V),n_vm,vm,PETSC_OWN_POINTER,S);
432: return(0);
433: }
436: /*@
437: VecISAXPY - Adds a reduced vector to the appropriate elements of a full-space vector.
438: vfull[is[i]] += alpha*vreduced[i]
440: Input Parameters:
441: + vfull - the full-space vector
442: . is - the index set for the reduced space
443: . alpha - the scalar coefficient
444: - vreduced - the reduced-space vector
446: Output Parameters:
447: . vfull - the sum of the full-space vector and reduced-space vector
449: Notes:
450: The index set identifies entries in the global vector.
451: Negative indices are skipped; indices outside the ownership range of vfull will raise an error.
453: Level: advanced
455: .seealso: VecAXPY(), VecGetOwnershipRange()
456: @*/
457: PetscErrorCode VecISAXPY(Vec vfull, IS is, PetscScalar alpha, Vec vreduced)
458: {
459: PetscInt nfull,nreduced;
466: VecGetSize(vfull,&nfull);
467: VecGetSize(vreduced,&nreduced);
469: if (nfull == nreduced) { /* Also takes care of masked vectors */
470: VecAXPY(vfull,alpha,vreduced);
471: } else {
472: PetscScalar *y;
473: const PetscScalar *x;
474: PetscInt i,n,m,rstart,rend;
475: const PetscInt *id;
477: VecGetArray(vfull,&y);
478: VecGetArrayRead(vreduced,&x);
479: ISGetIndices(is,&id);
480: ISGetLocalSize(is,&n);
481: VecGetLocalSize(vreduced,&m);
482: if (m != n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"IS local length not equal to Vec local length");
483: VecGetOwnershipRange(vfull,&rstart,&rend);
484: y -= rstart;
485: if (alpha == 1.0) {
486: for (i=0; i<n; ++i) {
487: if (id[i] < 0) continue;
488: if (id[i] < rstart || id[i] >= rend) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only owned values supported");
489: y[id[i]] += x[i];
490: }
491: } else {
492: for (i=0; i<n; ++i) {
493: if (id[i] < 0) continue;
494: if (id[i] < rstart || id[i] >= rend) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only owned values supported");
495: y[id[i]] += alpha*x[i];
496: }
497: }
498: y += rstart;
499: ISRestoreIndices(is,&id);
500: VecRestoreArray(vfull,&y);
501: VecRestoreArrayRead(vreduced,&x);
502: }
503: return(0);
504: }
506: /*@
507: VecISCopy - Copies between a reduced vector and the appropriate elements of a full-space vector.
509: Input Parameters:
510: + vfull - the full-space vector
511: . is - the index set for the reduced space
512: . mode - the direction of copying, SCATTER_FORWARD or SCATTER_REVERSE
513: - vreduced - the reduced-space vector
515: Output Parameters:
516: . vfull - the sum of the full-space vector and reduced-space vector
518: Notes:
519: The index set identifies entries in the global vector.
520: Negative indices are skipped; indices outside the ownership range of vfull will raise an error.
522: mode == SCATTER_FORWARD: vfull[is[i]] = vreduced[i]
523: mode == SCATTER_REVERSE: vreduced[i] = vfull[is[i]]
525: Level: advanced
527: .seealso: VecAXPY(), VecGetOwnershipRange()
528: @*/
529: PetscErrorCode VecISCopy(Vec vfull, IS is, ScatterMode mode, Vec vreduced)
530: {
531: PetscInt nfull, nreduced;
538: VecGetSize(vfull, &nfull);
539: VecGetSize(vreduced, &nreduced);
541: if (nfull == nreduced) { /* Also takes care of masked vectors */
542: if (mode == SCATTER_FORWARD) {
543: VecCopy(vreduced, vfull);
544: } else {
545: VecCopy(vfull, vreduced);
546: }
547: } else {
548: const PetscInt *id;
549: PetscInt i, n, m, rstart, rend;
551: ISGetIndices(is, &id);
552: ISGetLocalSize(is, &n);
553: VecGetLocalSize(vreduced, &m);
554: VecGetOwnershipRange(vfull, &rstart, &rend);
555: if (m != n) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "IS local length %D not equal to Vec local length %D", n, m);
556: if (mode == SCATTER_FORWARD) {
557: PetscScalar *y;
558: const PetscScalar *x;
560: VecGetArray(vfull, &y);
561: VecGetArrayRead(vreduced, &x);
562: y -= rstart;
563: for (i = 0; i < n; ++i) {
564: if (id[i] < 0) continue;
565: if (id[i] < rstart || id[i] >= rend) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only owned values supported");
566: y[id[i]] = x[i];
567: }
568: y += rstart;
569: VecRestoreArrayRead(vreduced, &x);
570: VecRestoreArray(vfull, &y);
571: } else if (mode == SCATTER_REVERSE) {
572: PetscScalar *x;
573: const PetscScalar *y;
575: VecGetArrayRead(vfull, &y);
576: VecGetArray(vreduced, &x);
577: for (i = 0; i < n; ++i) {
578: if (id[i] < 0) continue;
579: if (id[i] < rstart || id[i] >= rend) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only owned values supported");
580: x[i] = y[id[i]-rstart];
581: }
582: VecRestoreArray(vreduced, &x);
583: VecRestoreArrayRead(vfull, &y);
584: } else SETERRQ(PetscObjectComm((PetscObject) vfull), PETSC_ERR_ARG_WRONG, "Only forward or reverse modes are legal");
585: ISRestoreIndices(is, &id);
586: }
587: return(0);
588: }
590: /*@
591: ISComplementVec - Creates the complement of the index set relative to a layout defined by a Vec
593: Collective on IS
595: Input Parameter:
596: + S - a PETSc IS
597: - V - the reference vector space
599: Output Parameter:
600: . T - the complement of S
602: Level: advanced
604: .seealso: ISCreateGeneral()
605: @*/
606: PetscErrorCode ISComplementVec(IS S, Vec V, IS *T)
607: {
609: PetscInt start, end;
612: VecGetOwnershipRange(V,&start,&end);
613: ISComplement(S,start,end,T);
614: return(0);
615: }
617: /*@
618: VecISSet - Sets the elements of a vector, specified by an index set, to a constant
620: Input Parameter:
621: + V - the vector
622: . S - index set for the locations in the vector
623: - c - the constant
625: Notes:
626: The index set identifies entries in the global vector.
627: Negative indices are skipped; indices outside the ownership range of V will raise an error.
629: Level: advanced
631: .seealso: VecSet(), VecGetOwnershipRange()
632: @*/
633: PetscErrorCode VecISSet(Vec V,IS S, PetscScalar c)
634: {
636: PetscInt nloc,low,high,i;
637: const PetscInt *s;
638: PetscScalar *v;
641: if (!S) return(0); /* simply return with no-op if the index set is NULL */
646: VecGetOwnershipRange(V,&low,&high);
647: ISGetLocalSize(S,&nloc);
648: ISGetIndices(S,&s);
649: VecGetArray(V,&v);
650: for (i=0; i<nloc; ++i){
651: if (s[i] < 0) continue;
652: if (s[i] < low || s[i] >= high) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only owned values supported");
653: v[s[i]-low] = c;
654: }
655: ISRestoreIndices(S,&s);
656: VecRestoreArray(V,&v);
657: return(0);
658: }
660: #if !defined(PETSC_USE_COMPLEX)
661: /*@C
662: VecBoundGradientProjection - Projects vector according to this definition.
663: If XL[i] < X[i] < XU[i], then GP[i] = G[i];
664: If X[i] <= XL[i], then GP[i] = min(G[i],0);
665: If X[i] >= XU[i], then GP[i] = max(G[i],0);
667: Input Parameters:
668: + G - current gradient vector
669: . X - current solution vector with XL[i] <= X[i] <= XU[i]
670: . XL - lower bounds
671: - XU - upper bounds
673: Output Parameter:
674: . GP - gradient projection vector
676: Notes:
677: GP may be the same vector as G
679: Level: advanced
680: @*/
681: PetscErrorCode VecBoundGradientProjection(Vec G, Vec X, Vec XL, Vec XU, Vec GP)
682: {
684: PetscErrorCode ierr;
685: PetscInt n,i;
686: const PetscReal *xptr,*xlptr,*xuptr;
687: PetscReal *gptr,*gpptr;
688: PetscReal xval,gpval;
690: /* Project variables at the lower and upper bound */
698: VecGetLocalSize(X,&n);
700: VecGetArrayRead(X,&xptr);
701: VecGetArrayRead(XL,&xlptr);
702: VecGetArrayRead(XU,&xuptr);
703: VecGetArrayPair(G,GP,&gptr,&gpptr);
705: for (i=0; i<n; ++i){
706: gpval = gptr[i]; xval = xptr[i];
707: if (gpval>0.0 && xval<=xlptr[i]){
708: gpval = 0.0;
709: } else if (gpval<0.0 && xval>=xuptr[i]){
710: gpval = 0.0;
711: }
712: gpptr[i] = gpval;
713: }
715: VecRestoreArrayRead(X,&xptr);
716: VecRestoreArrayRead(XL,&xlptr);
717: VecRestoreArrayRead(XU,&xuptr);
718: VecRestoreArrayPair(G,GP,&gptr,&gpptr);
719: return(0);
720: }
721: #endif
723: /*@
724: VecStepMaxBounded - See below
726: Collective on Vec
728: Input Parameters:
729: + X - vector with no negative entries
730: . XL - lower bounds
731: . XU - upper bounds
732: - DX - step direction, can have negative, positive or zero entries
734: Output Parameter:
735: . stepmax - minimum value so that X[i] + stepmax*DX[i] <= XL[i] or XU[i] <= X[i] + stepmax*DX[i]
737: Level: intermediate
739: @*/
740: PetscErrorCode VecStepMaxBounded(Vec X, Vec DX, Vec XL, Vec XU, PetscReal *stepmax)
741: {
742: PetscErrorCode ierr;
743: PetscInt i,nn;
744: const PetscScalar *xx,*dx,*xl,*xu;
745: PetscReal localmax=0;
753: VecGetArrayRead(X,&xx);
754: VecGetArrayRead(XL,&xl);
755: VecGetArrayRead(XU,&xu);
756: VecGetArrayRead(DX,&dx);
757: VecGetLocalSize(X,&nn);
758: for (i=0;i<nn;i++){
759: if (PetscRealPart(dx[i]) > 0){
760: localmax=PetscMax(localmax,PetscRealPart((xu[i]-xx[i])/dx[i]));
761: } else if (PetscRealPart(dx[i])<0){
762: localmax=PetscMax(localmax,PetscRealPart((xl[i]-xx[i])/dx[i]));
763: }
764: }
765: VecRestoreArrayRead(X,&xx);
766: VecRestoreArrayRead(XL,&xl);
767: VecRestoreArrayRead(XU,&xu);
768: VecRestoreArrayRead(DX,&dx);
769: MPIU_Allreduce(&localmax,stepmax,1,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)X));
770: return(0);
771: }
773: /*@
774: VecStepBoundInfo - See below
776: Collective on Vec
778: Input Parameters:
779: + X - vector with no negative entries
780: . XL - lower bounds
781: . XU - upper bounds
782: - DX - step direction, can have negative, positive or zero entries
784: Output Parameter:
785: + boundmin - (may be NULL this it is not computed) maximum value so that XL[i] <= X[i] + boundmax*DX[i] <= XU[i]
786: . wolfemin - (may be NULL this it is not computed)
787: - 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]
789: Notes:
790: For complex numbers only compares the real part
792: Level: advanced
793: @*/
794: PetscErrorCode VecStepBoundInfo(Vec X, Vec DX, Vec XL, Vec XU, PetscReal *boundmin, PetscReal *wolfemin, PetscReal *boundmax)
795: {
796: PetscErrorCode ierr;
797: PetscInt n,i;
798: const PetscScalar *x,*xl,*xu,*dx;
799: PetscReal t;
800: PetscReal localmin=PETSC_INFINITY,localwolfemin=PETSC_INFINITY,localmax=-1;
801: MPI_Comm comm;
809: VecGetArrayRead(X,&x);
810: VecGetArrayRead(XL,&xl);
811: VecGetArrayRead(XU,&xu);
812: VecGetArrayRead(DX,&dx);
813: VecGetLocalSize(X,&n);
814: for (i=0; i<n; ++i){
815: if (PetscRealPart(dx[i])>0 && PetscRealPart(xu[i]) < PETSC_INFINITY) {
816: t=PetscRealPart((xu[i]-x[i])/dx[i]);
817: localmin=PetscMin(t,localmin);
818: if (localmin>0){
819: localwolfemin = PetscMin(t,localwolfemin);
820: }
821: localmax = PetscMax(t,localmax);
822: } else if (PetscRealPart(dx[i])<0 && PetscRealPart(xl[i]) > PETSC_NINFINITY) {
823: t=PetscRealPart((xl[i]-x[i])/dx[i]);
824: localmin = PetscMin(t,localmin);
825: if (localmin>0){
826: localwolfemin = PetscMin(t,localwolfemin);
827: }
828: localmax = PetscMax(t,localmax);
829: }
830: }
832: VecRestoreArrayRead(X,&x);
833: VecRestoreArrayRead(XL,&xl);
834: VecRestoreArrayRead(XU,&xu);
835: VecRestoreArrayRead(DX,&dx);
836: PetscObjectGetComm((PetscObject)X,&comm);
838: if (boundmin){
839: MPIU_Allreduce(&localmin,boundmin,1,MPIU_REAL,MPIU_MIN,comm);
840: PetscInfo1(X,"Step Bound Info: Closest Bound: %20.19e\n",(double)*boundmin);
841: }
842: if (wolfemin){
843: MPIU_Allreduce(&localwolfemin,wolfemin,1,MPIU_REAL,MPIU_MIN,comm);
844: PetscInfo1(X,"Step Bound Info: Wolfe: %20.19e\n",(double)*wolfemin);
845: }
846: if (boundmax) {
847: MPIU_Allreduce(&localmax,boundmax,1,MPIU_REAL,MPIU_MAX,comm);
848: if (*boundmax < 0) *boundmax=PETSC_INFINITY;
849: PetscInfo1(X,"Step Bound Info: Max: %20.19e\n",(double)*boundmax);
850: }
851: return(0);
852: }
854: /*@
855: VecStepMax - Returns the largest value so that x[i] + step*DX[i] >= 0 for all i
857: Collective on Vec
859: Input Parameters:
860: + X - vector with no negative entries
861: - DX - a step direction, can have negative, positive or zero entries
863: Output Parameter:
864: . step - largest value such that x[i] + step*DX[i] >= 0 for all i
866: Notes:
867: For complex numbers only compares the real part
869: Level: advanced
870: @*/
871: PetscErrorCode VecStepMax(Vec X, Vec DX, PetscReal *step)
872: {
873: PetscErrorCode ierr;
874: PetscInt i,nn;
875: PetscReal stepmax=PETSC_INFINITY;
876: const PetscScalar *xx,*dx;
882: VecGetLocalSize(X,&nn);
883: VecGetArrayRead(X,&xx);
884: VecGetArrayRead(DX,&dx);
885: for (i=0;i<nn;++i){
886: if (PetscRealPart(xx[i]) < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Vector must be positive");
887: else if (PetscRealPart(dx[i])<0) stepmax=PetscMin(stepmax,PetscRealPart(-xx[i]/dx[i]));
888: }
889: VecRestoreArrayRead(X,&xx);
890: VecRestoreArrayRead(DX,&dx);
891: MPIU_Allreduce(&stepmax,step,1,MPIU_REAL,MPIU_MIN,PetscObjectComm((PetscObject)X));
892: return(0);
893: }
895: /*@
896: VecPow - Replaces each component of a vector by x_i^p
898: Logically Collective on v
900: Input Parameter:
901: + v - the vector
902: - p - the exponent to use on each element
904: Output Parameter:
905: . v - the vector
907: Level: intermediate
909: @*/
910: PetscErrorCode VecPow(Vec v, PetscScalar p)
911: {
913: PetscInt n,i;
914: PetscScalar *v1;
919: VecGetArray(v,&v1);
920: VecGetLocalSize(v,&n);
922: if (1.0 == p) {
923: } else if (-1.0 == p) {
924: for (i = 0; i < n; ++i){
925: v1[i] = 1.0 / v1[i];
926: }
927: } else if (0.0 == p) {
928: for (i = 0; i < n; ++i){
929: /* Not-a-number left alone
930: Infinity set to one */
931: if (v1[i] == v1[i]) {
932: v1[i] = 1.0;
933: }
934: }
935: } else if (0.5 == p) {
936: for (i = 0; i < n; ++i) {
937: if (PetscRealPart(v1[i]) >= 0) {
938: v1[i] = PetscSqrtScalar(v1[i]);
939: } else {
940: v1[i] = PETSC_INFINITY;
941: }
942: }
943: } else if (-0.5 == p) {
944: for (i = 0; i < n; ++i) {
945: if (PetscRealPart(v1[i]) >= 0) {
946: v1[i] = 1.0 / PetscSqrtScalar(v1[i]);
947: } else {
948: v1[i] = PETSC_INFINITY;
949: }
950: }
951: } else if (2.0 == p) {
952: for (i = 0; i < n; ++i){
953: v1[i] *= v1[i];
954: }
955: } else if (-2.0 == p) {
956: for (i = 0; i < n; ++i){
957: v1[i] = 1.0 / (v1[i] * v1[i]);
958: }
959: } else {
960: for (i = 0; i < n; ++i) {
961: if (PetscRealPart(v1[i]) >= 0) {
962: v1[i] = PetscPowScalar(v1[i],p);
963: } else {
964: v1[i] = PETSC_INFINITY;
965: }
966: }
967: }
968: VecRestoreArray(v,&v1);
969: return(0);
970: }
972: /*@
973: VecMedian - Computes the componentwise median of three vectors
974: and stores the result in this vector. Used primarily for projecting
975: a vector within upper and lower bounds.
977: Logically Collective
979: Input Parameters:
980: . Vec1, Vec2, Vec3 - The three vectors
982: Output Parameter:
983: . VMedian - The median vector (this can be any one of the input vectors)
985: Level: advanced
986: @*/
987: PetscErrorCode VecMedian(Vec Vec1, Vec Vec2, Vec Vec3, Vec VMedian)
988: {
989: PetscErrorCode ierr;
990: PetscInt i,n,low1,high1;
991: const PetscScalar *v1,*v2,*v3;
992: PetscScalar *vmed;
1000: if (Vec1==Vec2 || Vec1==Vec3){
1001: VecCopy(Vec1,VMedian);
1002: return(0);
1003: }
1004: if (Vec2==Vec3){
1005: VecCopy(Vec2,VMedian);
1006: return(0);
1007: }
1009: /* Assert that Vec1 != Vec2 and Vec2 != Vec3 */
1020: VecCheckSameSize(Vec1,1,Vec2,2);
1021: VecCheckSameSize(Vec1,1,Vec3,3);
1022: VecCheckSameSize(Vec1,1,VMedian,4);
1024: VecGetOwnershipRange(Vec1,&low1,&high1);
1025: VecGetLocalSize(Vec1,&n);
1026: if (n>0){
1027: VecGetArray(VMedian,&vmed);
1028: if (Vec1 != VMedian){
1029: VecGetArrayRead(Vec1,&v1);
1030: } else {
1031: v1=vmed;
1032: }
1033: if (Vec2 != VMedian){
1034: VecGetArrayRead(Vec2,&v2);
1035: } else {
1036: v2=vmed;
1037: }
1038: if (Vec3 != VMedian){
1039: VecGetArrayRead(Vec3,&v3);
1040: } else {
1041: v3=vmed;
1042: }
1044: for (i=0;i<n;++i){
1045: 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])));
1046: }
1048: VecRestoreArray(VMedian,&vmed);
1049: if (VMedian != Vec1){
1050: VecRestoreArrayRead(Vec1,&v1);
1051: }
1052: if (VMedian != Vec2){
1053: VecRestoreArrayRead(Vec2,&v2);
1054: }
1055: if (VMedian != Vec3){
1056: VecRestoreArrayRead(Vec3,&v3);
1057: }
1058: }
1059: return(0);
1060: }