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
9: Input Parameters:
10: + Vec1 - the first vector to compare
11: - Vec2 - the second two vector to compare
13: OutputParameter:
14: . S - The index set containing the indices i where vec1[i] == vec2[i]
16: Level: advanced
18: Note:
19: The two vectors must have the same parallel layout
21: .seealso: `Vec`
22: @*/
23: PetscErrorCode VecWhichEqual(Vec Vec1, Vec Vec2, IS *S)
24: {
25: PetscInt i, n_same = 0;
26: PetscInt n, low, high;
27: PetscInt *same = NULL;
28: const PetscScalar *v1, *v2;
30: PetscFunctionBegin;
33: PetscCheckSameComm(Vec1, 1, Vec2, 2);
34: VecCheckSameSize(Vec1, 1, Vec2, 2);
36: PetscCall(VecGetOwnershipRange(Vec1, &low, &high));
37: PetscCall(VecGetLocalSize(Vec1, &n));
38: if (n > 0) {
39: if (Vec1 == Vec2) {
40: PetscCall(VecGetArrayRead(Vec1, &v1));
41: v2 = v1;
42: } else {
43: PetscCall(VecGetArrayRead(Vec1, &v1));
44: PetscCall(VecGetArrayRead(Vec2, &v2));
45: }
47: PetscCall(PetscMalloc1(n, &same));
49: for (i = 0; i < n; ++i) {
50: if (v1[i] == v2[i]) {
51: same[n_same] = low + i;
52: ++n_same;
53: }
54: }
56: if (Vec1 == Vec2) {
57: PetscCall(VecRestoreArrayRead(Vec1, &v1));
58: } else {
59: PetscCall(VecRestoreArrayRead(Vec1, &v1));
60: PetscCall(VecRestoreArrayRead(Vec2, &v2));
61: }
62: }
63: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)Vec1), n_same, same, PETSC_OWN_POINTER, S));
64: PetscFunctionReturn(PETSC_SUCCESS);
65: }
67: /*@
68: VecWhichLessThan - Creates an index set containing the indices
69: where the vectors `Vec1` < `Vec2`
71: Collective
73: Input Parameters:
74: + Vec1 - the first vector to compare
75: - Vec2 - the second vector to compare
77: OutputParameter:
78: . S - The index set containing the indices i where vec1[i] < vec2[i]
80: Level: advanced
82: Notes:
83: The two vectors must have the same parallel layout
85: For complex numbers this only compares the real part
87: .seealso: `Vec`
88: @*/
89: PetscErrorCode VecWhichLessThan(Vec Vec1, Vec Vec2, IS *S)
90: {
91: PetscInt i, n_lt = 0;
92: PetscInt n, low, high;
93: PetscInt *lt = NULL;
94: const PetscScalar *v1, *v2;
96: PetscFunctionBegin;
99: PetscCheckSameComm(Vec1, 1, Vec2, 2);
100: VecCheckSameSize(Vec1, 1, Vec2, 2);
102: PetscCall(VecGetOwnershipRange(Vec1, &low, &high));
103: PetscCall(VecGetLocalSize(Vec1, &n));
104: if (n > 0) {
105: if (Vec1 == Vec2) {
106: PetscCall(VecGetArrayRead(Vec1, &v1));
107: v2 = v1;
108: } else {
109: PetscCall(VecGetArrayRead(Vec1, &v1));
110: PetscCall(VecGetArrayRead(Vec2, &v2));
111: }
113: PetscCall(PetscMalloc1(n, <));
115: for (i = 0; i < n; ++i) {
116: if (PetscRealPart(v1[i]) < PetscRealPart(v2[i])) {
117: lt[n_lt] = low + i;
118: ++n_lt;
119: }
120: }
122: if (Vec1 == Vec2) {
123: PetscCall(VecRestoreArrayRead(Vec1, &v1));
124: } else {
125: PetscCall(VecRestoreArrayRead(Vec1, &v1));
126: PetscCall(VecRestoreArrayRead(Vec2, &v2));
127: }
128: }
129: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)Vec1), n_lt, lt, PETSC_OWN_POINTER, S));
130: PetscFunctionReturn(PETSC_SUCCESS);
131: }
133: /*@
134: VecWhichGreaterThan - Creates an index set containing the indices
135: where the vectors `Vec1` > `Vec2`
137: Collective
139: Input Parameters:
140: + Vec1 - the first vector to compare
141: - Vec2 - the second vector to compare
143: OutputParameter:
144: . S - The index set containing the indices i where vec1[i] > vec2[i]
146: Level: advanced
148: Notes:
149: The two vectors must have the same parallel layout
151: For complex numbers this only compares the real part
153: .seealso: `Vec`
154: @*/
155: PetscErrorCode VecWhichGreaterThan(Vec Vec1, Vec Vec2, IS *S)
156: {
157: PetscInt i, n_gt = 0;
158: PetscInt n, low, high;
159: PetscInt *gt = NULL;
160: const PetscScalar *v1, *v2;
162: PetscFunctionBegin;
165: PetscCheckSameComm(Vec1, 1, Vec2, 2);
166: VecCheckSameSize(Vec1, 1, Vec2, 2);
168: PetscCall(VecGetOwnershipRange(Vec1, &low, &high));
169: PetscCall(VecGetLocalSize(Vec1, &n));
170: if (n > 0) {
171: if (Vec1 == Vec2) {
172: PetscCall(VecGetArrayRead(Vec1, &v1));
173: v2 = v1;
174: } else {
175: PetscCall(VecGetArrayRead(Vec1, &v1));
176: PetscCall(VecGetArrayRead(Vec2, &v2));
177: }
179: PetscCall(PetscMalloc1(n, >));
181: for (i = 0; i < n; ++i) {
182: if (PetscRealPart(v1[i]) > PetscRealPart(v2[i])) {
183: gt[n_gt] = low + i;
184: ++n_gt;
185: }
186: }
188: if (Vec1 == Vec2) {
189: PetscCall(VecRestoreArrayRead(Vec1, &v1));
190: } else {
191: PetscCall(VecRestoreArrayRead(Vec1, &v1));
192: PetscCall(VecRestoreArrayRead(Vec2, &v2));
193: }
194: }
195: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)Vec1), n_gt, gt, PETSC_OWN_POINTER, S));
196: PetscFunctionReturn(PETSC_SUCCESS);
197: }
199: /*@
200: VecWhichBetween - Creates an index set containing the indices
201: where `VecLow` < `V` < `VecHigh`
203: Collective
205: Input Parameters:
206: + VecLow - lower bound
207: . V - Vector to compare
208: - VecHigh - higher bound
210: OutputParameter:
211: . S - The index set containing the indices i where veclow[i] < v[i] < vechigh[i]
213: Level: advanced
215: Notes:
216: The vectors must have the same parallel layout
218: For complex numbers this only compares the real part
220: .seealso: `Vec`
221: @*/
222: PetscErrorCode VecWhichBetween(Vec VecLow, Vec V, Vec VecHigh, IS *S)
223: {
224: PetscInt i, n_vm = 0;
225: PetscInt n, low, high;
226: PetscInt *vm = NULL;
227: const PetscScalar *v1, *v2, *vmiddle;
229: PetscFunctionBegin;
233: PetscCheckSameComm(V, 2, VecLow, 1);
234: PetscCheckSameComm(V, 2, VecHigh, 3);
235: VecCheckSameSize(V, 2, VecLow, 1);
236: VecCheckSameSize(V, 2, VecHigh, 3);
238: PetscCall(VecGetOwnershipRange(VecLow, &low, &high));
239: PetscCall(VecGetLocalSize(VecLow, &n));
240: if (n > 0) {
241: PetscCall(VecGetArrayRead(VecLow, &v1));
242: if (VecLow != VecHigh) {
243: PetscCall(VecGetArrayRead(VecHigh, &v2));
244: } else {
245: v2 = v1;
246: }
247: if (V != VecLow && V != VecHigh) {
248: PetscCall(VecGetArrayRead(V, &vmiddle));
249: } else if (V == VecLow) {
250: vmiddle = v1;
251: } else {
252: vmiddle = v2;
253: }
255: PetscCall(PetscMalloc1(n, &vm));
257: for (i = 0; i < n; ++i) {
258: if (PetscRealPart(v1[i]) < PetscRealPart(vmiddle[i]) && PetscRealPart(vmiddle[i]) < PetscRealPart(v2[i])) {
259: vm[n_vm] = low + i;
260: ++n_vm;
261: }
262: }
264: PetscCall(VecRestoreArrayRead(VecLow, &v1));
265: if (VecLow != VecHigh) PetscCall(VecRestoreArrayRead(VecHigh, &v2));
266: if (V != VecLow && V != VecHigh) PetscCall(VecRestoreArrayRead(V, &vmiddle));
267: }
268: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)V), n_vm, vm, PETSC_OWN_POINTER, S));
269: PetscFunctionReturn(PETSC_SUCCESS);
270: }
272: /*@
273: VecWhichBetweenOrEqual - Creates an index set containing the indices
274: where `VecLow` <= `V` <= `VecHigh`
276: Collective
278: Input Parameters:
279: + VecLow - lower bound
280: . V - Vector to compare
281: - VecHigh - higher bound
283: OutputParameter:
284: . S - The index set containing the indices i where veclow[i] <= v[i] <= vechigh[i]
286: Level: advanced
288: .seealso: `Vec`
289: @*/
291: PetscErrorCode VecWhichBetweenOrEqual(Vec VecLow, Vec V, Vec VecHigh, IS *S)
292: {
293: PetscInt i, n_vm = 0;
294: PetscInt n, low, high;
295: PetscInt *vm = NULL;
296: const PetscScalar *v1, *v2, *vmiddle;
298: PetscFunctionBegin;
302: PetscCheckSameComm(V, 2, VecLow, 1);
303: PetscCheckSameComm(V, 2, VecHigh, 3);
304: VecCheckSameSize(V, 2, VecLow, 1);
305: VecCheckSameSize(V, 2, VecHigh, 3);
307: PetscCall(VecGetOwnershipRange(VecLow, &low, &high));
308: PetscCall(VecGetLocalSize(VecLow, &n));
309: if (n > 0) {
310: PetscCall(VecGetArrayRead(VecLow, &v1));
311: if (VecLow != VecHigh) {
312: PetscCall(VecGetArrayRead(VecHigh, &v2));
313: } else {
314: v2 = v1;
315: }
316: if (V != VecLow && V != VecHigh) {
317: PetscCall(VecGetArrayRead(V, &vmiddle));
318: } else if (V == VecLow) {
319: vmiddle = v1;
320: } else {
321: vmiddle = v2;
322: }
324: PetscCall(PetscMalloc1(n, &vm));
326: for (i = 0; i < n; ++i) {
327: if (PetscRealPart(v1[i]) <= PetscRealPart(vmiddle[i]) && PetscRealPart(vmiddle[i]) <= PetscRealPart(v2[i])) {
328: vm[n_vm] = low + i;
329: ++n_vm;
330: }
331: }
333: PetscCall(VecRestoreArrayRead(VecLow, &v1));
334: if (VecLow != VecHigh) PetscCall(VecRestoreArrayRead(VecHigh, &v2));
335: if (V != VecLow && V != VecHigh) PetscCall(VecRestoreArrayRead(V, &vmiddle));
336: }
337: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)V), n_vm, vm, PETSC_OWN_POINTER, S));
338: PetscFunctionReturn(PETSC_SUCCESS);
339: }
341: /*@
342: VecWhichInactive - Creates an index set containing the indices
343: where one of the following holds:
344: a) VecLow(i) < V(i) < VecHigh(i)
345: b) VecLow(i) = V(i) and D(i) <= 0 (< 0 when Strong is true)
346: c) VecHigh(i) = V(i) and D(i) >= 0 (> 0 when Strong is true)
348: Collective
350: Input Parameters:
351: + VecLow - lower bound
352: . V - Vector to compare
353: . D - Direction to compare
354: . VecHigh - higher bound
355: - Strong - indicator for applying strongly inactive test
357: OutputParameter:
358: . S - The index set containing the indices i where the bound is inactive
360: Level: advanced
362: .seealso: `Vec`
363: @*/
365: PetscErrorCode VecWhichInactive(Vec VecLow, Vec V, Vec D, Vec VecHigh, PetscBool Strong, IS *S)
366: {
367: PetscInt i, n_vm = 0;
368: PetscInt n, low, high;
369: PetscInt *vm = NULL;
370: const PetscScalar *v1, *v2, *v, *d;
372: PetscFunctionBegin;
373: if (!VecLow && !VecHigh) {
374: PetscCall(VecGetOwnershipRange(V, &low, &high));
375: PetscCall(ISCreateStride(PetscObjectComm((PetscObject)V), high - low, low, 1, S));
376: PetscFunctionReturn(PETSC_SUCCESS);
377: }
382: PetscCheckSameComm(V, 2, VecLow, 1);
383: PetscCheckSameComm(V, 2, D, 3);
384: PetscCheckSameComm(V, 2, VecHigh, 4);
385: VecCheckSameSize(V, 2, VecLow, 1);
386: VecCheckSameSize(V, 2, D, 3);
387: VecCheckSameSize(V, 2, VecHigh, 4);
389: PetscCall(VecGetOwnershipRange(VecLow, &low, &high));
390: PetscCall(VecGetLocalSize(VecLow, &n));
391: if (n > 0) {
392: PetscCall(VecGetArrayRead(VecLow, &v1));
393: if (VecLow != VecHigh) {
394: PetscCall(VecGetArrayRead(VecHigh, &v2));
395: } else {
396: v2 = v1;
397: }
398: if (V != VecLow && V != VecHigh) {
399: PetscCall(VecGetArrayRead(V, &v));
400: } else if (V == VecLow) {
401: v = v1;
402: } else {
403: v = v2;
404: }
405: if (D != VecLow && D != VecHigh && D != V) {
406: PetscCall(VecGetArrayRead(D, &d));
407: } else if (D == VecLow) {
408: d = v1;
409: } else if (D == VecHigh) {
410: d = v2;
411: } else {
412: d = v;
413: }
415: PetscCall(PetscMalloc1(n, &vm));
417: if (Strong) {
418: for (i = 0; i < n; ++i) {
419: if (PetscRealPart(v1[i]) < PetscRealPart(v[i]) && PetscRealPart(v[i]) < PetscRealPart(v2[i])) {
420: vm[n_vm] = low + i;
421: ++n_vm;
422: } else if (PetscRealPart(v1[i]) == PetscRealPart(v[i]) && PetscRealPart(d[i]) < 0) {
423: vm[n_vm] = low + i;
424: ++n_vm;
425: } else if (PetscRealPart(v2[i]) == PetscRealPart(v[i]) && PetscRealPart(d[i]) > 0) {
426: vm[n_vm] = low + i;
427: ++n_vm;
428: }
429: }
430: } else {
431: for (i = 0; i < n; ++i) {
432: if (PetscRealPart(v1[i]) < PetscRealPart(v[i]) && PetscRealPart(v[i]) < PetscRealPart(v2[i])) {
433: vm[n_vm] = low + i;
434: ++n_vm;
435: } else if (PetscRealPart(v1[i]) == PetscRealPart(v[i]) && PetscRealPart(d[i]) <= 0) {
436: vm[n_vm] = low + i;
437: ++n_vm;
438: } else if (PetscRealPart(v2[i]) == PetscRealPart(v[i]) && PetscRealPart(d[i]) >= 0) {
439: vm[n_vm] = low + i;
440: ++n_vm;
441: }
442: }
443: }
445: PetscCall(VecRestoreArrayRead(VecLow, &v1));
446: if (VecLow != VecHigh) PetscCall(VecRestoreArrayRead(VecHigh, &v2));
447: if (V != VecLow && V != VecHigh) PetscCall(VecRestoreArrayRead(V, &v));
448: if (D != VecLow && D != VecHigh && D != V) PetscCall(VecRestoreArrayRead(D, &d));
449: }
450: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)V), n_vm, vm, PETSC_OWN_POINTER, S));
451: PetscFunctionReturn(PETSC_SUCCESS);
452: }
454: /*@
455: VecISAXPY - Adds a reduced vector to the appropriate elements of a full-space vector.
456: vfull[is[i]] += alpha*vreduced[i]
458: Input Parameters:
459: + vfull - the full-space vector
460: . is - the index set for the reduced space
461: . alpha - the scalar coefficient
462: - vreduced - the reduced-space vector
464: Output Parameter:
465: . vfull - the sum of the full-space vector and reduced-space vector
467: Level: advanced
469: Notes:
470: The index set identifies entries in the global vector.
471: Negative indices are skipped; indices outside the ownership range of `vfull` will raise an error.
473: .seealso: `VecISCopy()`, `VecISSet()`, `VecAXPY()`
474: @*/
475: PetscErrorCode VecISAXPY(Vec vfull, IS is, PetscScalar alpha, Vec vreduced)
476: {
477: PetscInt nfull, nreduced;
478: PetscBool sorted = PETSC_FALSE;
480: PetscFunctionBegin;
484: PetscCall(VecGetSize(vfull, &nfull));
485: PetscCall(VecGetSize(vreduced, &nreduced));
486: if (nfull == nreduced) PetscCall(ISGetInfo(is, IS_SORTED, IS_GLOBAL, PETSC_TRUE, &sorted));
487: if (sorted) PetscCall(VecAXPY(vfull, alpha, vreduced));
488: else {
489: PetscScalar *y;
490: const PetscScalar *x;
491: PetscInt i, n, m, rstart, rend;
492: const PetscInt *id;
494: PetscCall(VecGetArray(vfull, &y));
495: PetscCall(VecGetArrayRead(vreduced, &x));
496: PetscCall(ISGetIndices(is, &id));
497: PetscCall(ISGetLocalSize(is, &n));
498: PetscCall(VecGetLocalSize(vreduced, &m));
499: PetscCheck(m == n, PETSC_COMM_SELF, PETSC_ERR_SUP, "IS local length %" PetscInt_FMT " not equal to Vec local length %" PetscInt_FMT, n, m);
500: PetscCall(VecGetOwnershipRange(vfull, &rstart, &rend));
501: y -= rstart;
502: if (alpha == 1.0) {
503: for (i = 0; i < n; ++i) {
504: if (id[i] < 0) continue;
505: PetscCheck(id[i] >= rstart && id[i] < rend, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only owned values supported");
506: y[id[i]] += x[i];
507: }
508: } else {
509: for (i = 0; i < n; ++i) {
510: if (id[i] < 0) continue;
511: PetscCheck(id[i] >= rstart && id[i] < rend, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only owned values supported");
512: y[id[i]] += alpha * x[i];
513: }
514: }
515: y += rstart;
516: PetscCall(ISRestoreIndices(is, &id));
517: PetscCall(VecRestoreArray(vfull, &y));
518: PetscCall(VecRestoreArrayRead(vreduced, &x));
519: }
520: PetscFunctionReturn(PETSC_SUCCESS);
521: }
523: /*@
524: VecISCopy - Copies between a reduced vector and the appropriate elements of a full-space vector.
526: Input Parameters:
527: + vfull - the full-space vector
528: . is - the index set for the reduced space
529: . mode - the direction of copying, `SCATTER_FORWARD` or `SCATTER_REVERSE`
530: - vreduced - the reduced-space vector
532: Output Parameter:
533: . vfull - the sum of the full-space vector and reduced-space vector
535: Level: advanced
537: Notes:
538: The index set identifies entries in the global vector.
539: Negative indices are skipped; indices outside the ownership range of `vfull` will raise an error.
540: .vb
541: mode == SCATTER_FORWARD: vfull[is[i]] = vreduced[i]
542: mode == SCATTER_REVERSE: vreduced[i] = vfull[is[i]]
543: .ve
545: .seealso: `VecISSet()`, `VecISAXPY()`, `VecCopy()`
546: @*/
547: PetscErrorCode VecISCopy(Vec vfull, IS is, ScatterMode mode, Vec vreduced)
548: {
549: PetscInt nfull, nreduced;
550: PetscBool sorted = PETSC_FALSE;
552: PetscFunctionBegin;
556: PetscCall(VecGetSize(vfull, &nfull));
557: PetscCall(VecGetSize(vreduced, &nreduced));
558: if (nfull == nreduced) PetscCall(ISGetInfo(is, IS_SORTED, IS_GLOBAL, PETSC_TRUE, &sorted));
559: if (sorted) {
560: if (mode == SCATTER_FORWARD) {
561: PetscCall(VecCopy(vreduced, vfull));
562: } else {
563: PetscCall(VecCopy(vfull, vreduced));
564: }
565: } else {
566: const PetscInt *id;
567: PetscInt i, n, m, rstart, rend;
569: PetscCall(ISGetIndices(is, &id));
570: PetscCall(ISGetLocalSize(is, &n));
571: PetscCall(VecGetLocalSize(vreduced, &m));
572: PetscCall(VecGetOwnershipRange(vfull, &rstart, &rend));
573: PetscCheck(m == n, PETSC_COMM_SELF, PETSC_ERR_SUP, "IS local length %" PetscInt_FMT " not equal to Vec local length %" PetscInt_FMT, n, m);
574: if (mode == SCATTER_FORWARD) {
575: PetscScalar *y;
576: const PetscScalar *x;
578: PetscCall(VecGetArray(vfull, &y));
579: PetscCall(VecGetArrayRead(vreduced, &x));
580: y -= rstart;
581: for (i = 0; i < n; ++i) {
582: if (id[i] < 0) continue;
583: PetscCheck(id[i] >= rstart && id[i] < rend, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only owned values supported");
584: y[id[i]] = x[i];
585: }
586: y += rstart;
587: PetscCall(VecRestoreArrayRead(vreduced, &x));
588: PetscCall(VecRestoreArray(vfull, &y));
589: } else if (mode == SCATTER_REVERSE) {
590: PetscScalar *x;
591: const PetscScalar *y;
593: PetscCall(VecGetArrayRead(vfull, &y));
594: PetscCall(VecGetArray(vreduced, &x));
595: for (i = 0; i < n; ++i) {
596: if (id[i] < 0) continue;
597: PetscCheck(id[i] >= rstart && id[i] < rend, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only owned values supported");
598: x[i] = y[id[i] - rstart];
599: }
600: PetscCall(VecRestoreArray(vreduced, &x));
601: PetscCall(VecRestoreArrayRead(vfull, &y));
602: } else SETERRQ(PetscObjectComm((PetscObject)vfull), PETSC_ERR_ARG_WRONG, "Only forward or reverse modes are legal");
603: PetscCall(ISRestoreIndices(is, &id));
604: }
605: PetscFunctionReturn(PETSC_SUCCESS);
606: }
608: /*@
609: ISComplementVec - Creates the complement of the index set relative to a layout defined by a `Vec`
611: Collective
613: Input Parameters:
614: + S - a PETSc `IS`
615: - V - the reference vector space
617: Output Parameter:
618: . T - the complement of S
620: Level: advanced
622: .seealso: `IS`, `Vec`, `ISCreateGeneral()`
623: @*/
624: PetscErrorCode ISComplementVec(IS S, Vec V, IS *T)
625: {
626: PetscInt start, end;
628: PetscFunctionBegin;
629: PetscCall(VecGetOwnershipRange(V, &start, &end));
630: PetscCall(ISComplement(S, start, end, T));
631: PetscFunctionReturn(PETSC_SUCCESS);
632: }
634: /*@
635: VecISSet - Sets the elements of a vector, specified by an index set, to a constant
637: Input Parameters:
638: + V - the vector
639: . S - index set for the locations in the vector
640: - c - the constant
642: Level: advanced
644: Notes:
645: The index set identifies entries in the global vector.
646: Negative indices are skipped; indices outside the ownership range of V will raise an error.
648: .seealso: `VecISCopy()`, `VecISAXPY()`, `VecSet()`
649: @*/
650: PetscErrorCode VecISSet(Vec V, IS S, PetscScalar c)
651: {
652: PetscInt nloc, low, high, i;
653: const PetscInt *s;
654: PetscScalar *v;
656: PetscFunctionBegin;
657: if (!S) PetscFunctionReturn(PETSC_SUCCESS); /* simply return with no-op if the index set is NULL */
662: PetscCall(VecGetOwnershipRange(V, &low, &high));
663: PetscCall(ISGetLocalSize(S, &nloc));
664: PetscCall(ISGetIndices(S, &s));
665: PetscCall(VecGetArray(V, &v));
666: for (i = 0; i < nloc; ++i) {
667: if (s[i] < 0) continue;
668: PetscCheck(s[i] >= low && s[i] < high, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only owned values supported");
669: v[s[i] - low] = c;
670: }
671: PetscCall(ISRestoreIndices(S, &s));
672: PetscCall(VecRestoreArray(V, &v));
673: PetscFunctionReturn(PETSC_SUCCESS);
674: }
676: #if !defined(PETSC_USE_COMPLEX)
677: /*@C
678: VecBoundGradientProjection - Projects vector according to this definition.
679: If XL[i] < X[i] < XU[i], then GP[i] = G[i];
680: If X[i] <= XL[i], then GP[i] = min(G[i],0);
681: If X[i] >= XU[i], then GP[i] = max(G[i],0);
683: Input Parameters:
684: + G - current gradient vector
685: . X - current solution vector with XL[i] <= X[i] <= XU[i]
686: . XL - lower bounds
687: - XU - upper bounds
689: Output Parameter:
690: . GP - gradient projection vector
692: Level: advanced
694: Note:
695: `GP` may be the same vector as `G`
697: .seealso: `Vec`
698: @*/
699: PetscErrorCode VecBoundGradientProjection(Vec G, Vec X, Vec XL, Vec XU, Vec GP)
700: {
701: PetscInt n, i;
702: const PetscReal *xptr, *xlptr, *xuptr;
703: PetscReal *gptr, *gpptr;
704: PetscReal xval, gpval;
706: /* Project variables at the lower and upper bound */
707: PetscFunctionBegin;
713: if (!XL && !XU) {
714: PetscCall(VecCopy(G, GP));
715: PetscFunctionReturn(PETSC_SUCCESS);
716: }
718: PetscCall(VecGetLocalSize(X, &n));
720: PetscCall(VecGetArrayRead(X, &xptr));
721: PetscCall(VecGetArrayRead(XL, &xlptr));
722: PetscCall(VecGetArrayRead(XU, &xuptr));
723: PetscCall(VecGetArrayPair(G, GP, &gptr, &gpptr));
725: for (i = 0; i < n; ++i) {
726: gpval = gptr[i];
727: xval = xptr[i];
728: if (gpval > 0.0 && xval <= xlptr[i]) {
729: gpval = 0.0;
730: } else if (gpval < 0.0 && xval >= xuptr[i]) {
731: gpval = 0.0;
732: }
733: gpptr[i] = gpval;
734: }
736: PetscCall(VecRestoreArrayRead(X, &xptr));
737: PetscCall(VecRestoreArrayRead(XL, &xlptr));
738: PetscCall(VecRestoreArrayRead(XU, &xuptr));
739: PetscCall(VecRestoreArrayPair(G, GP, &gptr, &gpptr));
740: PetscFunctionReturn(PETSC_SUCCESS);
741: }
742: #endif
744: /*@
745: VecStepMaxBounded - See below
747: Collective
749: Input Parameters:
750: + X - vector with no negative entries
751: . XL - lower bounds
752: . XU - upper bounds
753: - DX - step direction, can have negative, positive or zero entries
755: Output Parameter:
756: . stepmax - minimum value so that X[i] + stepmax*DX[i] <= XL[i] or XU[i] <= X[i] + stepmax*DX[i]
758: Level: intermediate
760: .seealso: `Vec`
761: @*/
762: PetscErrorCode VecStepMaxBounded(Vec X, Vec DX, Vec XL, Vec XU, PetscReal *stepmax)
763: {
764: PetscInt i, nn;
765: const PetscScalar *xx, *dx, *xl, *xu;
766: PetscReal localmax = 0;
768: PetscFunctionBegin;
774: PetscCall(VecGetArrayRead(X, &xx));
775: PetscCall(VecGetArrayRead(XL, &xl));
776: PetscCall(VecGetArrayRead(XU, &xu));
777: PetscCall(VecGetArrayRead(DX, &dx));
778: PetscCall(VecGetLocalSize(X, &nn));
779: for (i = 0; i < nn; i++) {
780: if (PetscRealPart(dx[i]) > 0) {
781: localmax = PetscMax(localmax, PetscRealPart((xu[i] - xx[i]) / dx[i]));
782: } else if (PetscRealPart(dx[i]) < 0) {
783: localmax = PetscMax(localmax, PetscRealPart((xl[i] - xx[i]) / dx[i]));
784: }
785: }
786: PetscCall(VecRestoreArrayRead(X, &xx));
787: PetscCall(VecRestoreArrayRead(XL, &xl));
788: PetscCall(VecRestoreArrayRead(XU, &xu));
789: PetscCall(VecRestoreArrayRead(DX, &dx));
790: PetscCall(MPIU_Allreduce(&localmax, stepmax, 1, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)X)));
791: PetscFunctionReturn(PETSC_SUCCESS);
792: }
794: /*@
795: VecStepBoundInfo - See below
797: Collective
799: Input Parameters:
800: + X - vector with no negative entries
801: . XL - lower bounds
802: . XU - upper bounds
803: - DX - step direction, can have negative, positive or zero entries
805: Output Parameters:
806: + boundmin - (may be `NULL` this it is not computed) maximum value so that XL[i] <= X[i] + boundmax*DX[i] <= XU[i]
807: . wolfemin - (may be `NULL` this it is not computed)
808: - 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]
810: Level: advanced
812: Note:
813: For complex numbers only compares the real part
815: .seealso: `Vec`
816: @*/
817: PetscErrorCode VecStepBoundInfo(Vec X, Vec DX, Vec XL, Vec XU, PetscReal *boundmin, PetscReal *wolfemin, PetscReal *boundmax)
818: {
819: PetscInt n, i;
820: const PetscScalar *x, *xl, *xu, *dx;
821: PetscReal t;
822: PetscReal localmin = PETSC_INFINITY, localwolfemin = PETSC_INFINITY, localmax = -1;
823: MPI_Comm comm;
825: PetscFunctionBegin;
831: PetscCall(VecGetArrayRead(X, &x));
832: PetscCall(VecGetArrayRead(XL, &xl));
833: PetscCall(VecGetArrayRead(XU, &xu));
834: PetscCall(VecGetArrayRead(DX, &dx));
835: PetscCall(VecGetLocalSize(X, &n));
836: for (i = 0; i < n; ++i) {
837: if (PetscRealPart(dx[i]) > 0 && PetscRealPart(xu[i]) < PETSC_INFINITY) {
838: t = PetscRealPart((xu[i] - x[i]) / dx[i]);
839: localmin = PetscMin(t, localmin);
840: if (localmin > 0) localwolfemin = PetscMin(t, localwolfemin);
841: localmax = PetscMax(t, localmax);
842: } else if (PetscRealPart(dx[i]) < 0 && PetscRealPart(xl[i]) > PETSC_NINFINITY) {
843: t = PetscRealPart((xl[i] - x[i]) / dx[i]);
844: localmin = PetscMin(t, localmin);
845: if (localmin > 0) localwolfemin = PetscMin(t, localwolfemin);
846: localmax = PetscMax(t, localmax);
847: }
848: }
850: PetscCall(VecRestoreArrayRead(X, &x));
851: PetscCall(VecRestoreArrayRead(XL, &xl));
852: PetscCall(VecRestoreArrayRead(XU, &xu));
853: PetscCall(VecRestoreArrayRead(DX, &dx));
854: PetscCall(PetscObjectGetComm((PetscObject)X, &comm));
856: if (boundmin) {
857: PetscCall(MPIU_Allreduce(&localmin, boundmin, 1, MPIU_REAL, MPIU_MIN, comm));
858: PetscCall(PetscInfo(X, "Step Bound Info: Closest Bound: %20.19e\n", (double)*boundmin));
859: }
860: if (wolfemin) {
861: PetscCall(MPIU_Allreduce(&localwolfemin, wolfemin, 1, MPIU_REAL, MPIU_MIN, comm));
862: PetscCall(PetscInfo(X, "Step Bound Info: Wolfe: %20.19e\n", (double)*wolfemin));
863: }
864: if (boundmax) {
865: PetscCall(MPIU_Allreduce(&localmax, boundmax, 1, MPIU_REAL, MPIU_MAX, comm));
866: if (*boundmax < 0) *boundmax = PETSC_INFINITY;
867: PetscCall(PetscInfo(X, "Step Bound Info: Max: %20.19e\n", (double)*boundmax));
868: }
869: PetscFunctionReturn(PETSC_SUCCESS);
870: }
872: /*@
873: VecStepMax - Returns the largest value so that x[i] + step*DX[i] >= 0 for all i
875: Collective
877: Input Parameters:
878: + X - vector with no negative entries
879: - DX - a step direction, can have negative, positive or zero entries
881: Output Parameter:
882: . step - largest value such that x[i] + step*DX[i] >= 0 for all i
884: Level: advanced
886: Note:
887: For complex numbers only compares the real part
889: .seealso: `Vec`
890: @*/
891: PetscErrorCode VecStepMax(Vec X, Vec DX, PetscReal *step)
892: {
893: PetscInt i, nn;
894: PetscReal stepmax = PETSC_INFINITY;
895: const PetscScalar *xx, *dx;
897: PetscFunctionBegin;
901: PetscCall(VecGetLocalSize(X, &nn));
902: PetscCall(VecGetArrayRead(X, &xx));
903: PetscCall(VecGetArrayRead(DX, &dx));
904: for (i = 0; i < nn; ++i) {
905: PetscCheck(PetscRealPart(xx[i]) >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Vector must be positive");
906: if (PetscRealPart(dx[i]) < 0) stepmax = PetscMin(stepmax, PetscRealPart(-xx[i] / dx[i]));
907: }
908: PetscCall(VecRestoreArrayRead(X, &xx));
909: PetscCall(VecRestoreArrayRead(DX, &dx));
910: PetscCall(MPIU_Allreduce(&stepmax, step, 1, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)X)));
911: PetscFunctionReturn(PETSC_SUCCESS);
912: }
914: /*@
915: VecPow - Replaces each component of a vector by x_i^p
917: Logically Collective
919: Input Parameters:
920: + v - the vector
921: - p - the exponent to use on each element
923: Level: intermediate
925: .seealso: `Vec`
926: @*/
927: PetscErrorCode VecPow(Vec v, PetscScalar p)
928: {
929: PetscInt n, i;
930: PetscScalar *v1;
932: PetscFunctionBegin;
935: PetscCall(VecGetArray(v, &v1));
936: PetscCall(VecGetLocalSize(v, &n));
938: if (1.0 == p) {
939: } else if (-1.0 == p) {
940: for (i = 0; i < n; ++i) v1[i] = 1.0 / v1[i];
941: } else if (0.0 == p) {
942: for (i = 0; i < n; ++i) {
943: /* Not-a-number left alone
944: Infinity set to one */
945: if (v1[i] == v1[i]) v1[i] = 1.0;
946: }
947: } else if (0.5 == p) {
948: for (i = 0; i < n; ++i) {
949: if (PetscRealPart(v1[i]) >= 0) {
950: v1[i] = PetscSqrtScalar(v1[i]);
951: } else {
952: v1[i] = PETSC_INFINITY;
953: }
954: }
955: } else if (-0.5 == p) {
956: for (i = 0; i < n; ++i) {
957: if (PetscRealPart(v1[i]) >= 0) {
958: v1[i] = 1.0 / PetscSqrtScalar(v1[i]);
959: } else {
960: v1[i] = PETSC_INFINITY;
961: }
962: }
963: } else if (2.0 == p) {
964: for (i = 0; i < n; ++i) v1[i] *= v1[i];
965: } else if (-2.0 == p) {
966: for (i = 0; i < n; ++i) v1[i] = 1.0 / (v1[i] * v1[i]);
967: } else {
968: for (i = 0; i < n; ++i) {
969: if (PetscRealPart(v1[i]) >= 0) {
970: v1[i] = PetscPowScalar(v1[i], p);
971: } else {
972: v1[i] = PETSC_INFINITY;
973: }
974: }
975: }
976: PetscCall(VecRestoreArray(v, &v1));
977: PetscFunctionReturn(PETSC_SUCCESS);
978: }
980: /*@
981: VecMedian - Computes the componentwise median of three vectors
982: and stores the result in this vector. Used primarily for projecting
983: a vector within upper and lower bounds.
985: Logically Collective
987: Input Parameters:
988: + Vec1 - The first vector
989: . Vec2 - The second vector
990: - Vec3 - The third vector
992: Output Parameter:
993: . VMedian - The median vector (this can be any one of the input vectors)
995: Level: advanced
997: .seealso: `Vec`
998: @*/
999: PetscErrorCode VecMedian(Vec Vec1, Vec Vec2, Vec Vec3, Vec VMedian)
1000: {
1001: PetscInt i, n, low1, high1;
1002: const PetscScalar *v1, *v2, *v3;
1003: PetscScalar *vmed;
1005: PetscFunctionBegin;
1011: if (!Vec1 && !Vec3) {
1012: PetscCall(VecCopy(Vec2, VMedian));
1013: PetscFunctionReturn(PETSC_SUCCESS);
1014: }
1015: if (Vec1 == Vec2 || Vec1 == Vec3) {
1016: PetscCall(VecCopy(Vec1, VMedian));
1017: PetscFunctionReturn(PETSC_SUCCESS);
1018: }
1019: if (Vec2 == Vec3) {
1020: PetscCall(VecCopy(Vec2, VMedian));
1021: PetscFunctionReturn(PETSC_SUCCESS);
1022: }
1024: /* Assert that Vec1 != Vec2 and Vec2 != Vec3 */
1029: PetscCheckSameType(Vec1, 1, Vec2, 2);
1030: PetscCheckSameType(Vec1, 1, Vec3, 3);
1031: PetscCheckSameType(Vec1, 1, VMedian, 4);
1032: PetscCheckSameComm(Vec1, 1, Vec2, 2);
1033: PetscCheckSameComm(Vec1, 1, Vec3, 3);
1034: PetscCheckSameComm(Vec1, 1, VMedian, 4);
1035: VecCheckSameSize(Vec1, 1, Vec2, 2);
1036: VecCheckSameSize(Vec1, 1, Vec3, 3);
1037: VecCheckSameSize(Vec1, 1, VMedian, 4);
1039: PetscCall(VecGetOwnershipRange(Vec1, &low1, &high1));
1040: PetscCall(VecGetLocalSize(Vec1, &n));
1041: if (n > 0) {
1042: PetscCall(VecGetArray(VMedian, &vmed));
1043: if (Vec1 != VMedian) {
1044: PetscCall(VecGetArrayRead(Vec1, &v1));
1045: } else {
1046: v1 = vmed;
1047: }
1048: if (Vec2 != VMedian) {
1049: PetscCall(VecGetArrayRead(Vec2, &v2));
1050: } else {
1051: v2 = vmed;
1052: }
1053: if (Vec3 != VMedian) {
1054: PetscCall(VecGetArrayRead(Vec3, &v3));
1055: } else {
1056: v3 = vmed;
1057: }
1059: 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])));
1061: PetscCall(VecRestoreArray(VMedian, &vmed));
1062: if (VMedian != Vec1) PetscCall(VecRestoreArrayRead(Vec1, &v1));
1063: if (VMedian != Vec2) PetscCall(VecRestoreArrayRead(Vec2, &v2));
1064: if (VMedian != Vec3) PetscCall(VecRestoreArrayRead(Vec3, &v3));
1065: }
1066: PetscFunctionReturn(PETSC_SUCCESS);
1067: }