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: Output Parameter:
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: Output Parameter:
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: Output Parameter:
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: Output Parameter:
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: Output Parameter:
284: . S - The index set containing the indices i where veclow[i] <= v[i] <= vechigh[i]
286: Level: advanced
288: .seealso: `Vec`
289: @*/
290: PetscErrorCode VecWhichBetweenOrEqual(Vec VecLow, Vec V, Vec VecHigh, IS *S)
291: {
292: PetscInt i, n_vm = 0;
293: PetscInt n, low, high;
294: PetscInt *vm = NULL;
295: const PetscScalar *v1, *v2, *vmiddle;
297: PetscFunctionBegin;
301: PetscCheckSameComm(V, 2, VecLow, 1);
302: PetscCheckSameComm(V, 2, VecHigh, 3);
303: VecCheckSameSize(V, 2, VecLow, 1);
304: VecCheckSameSize(V, 2, VecHigh, 3);
306: PetscCall(VecGetOwnershipRange(VecLow, &low, &high));
307: PetscCall(VecGetLocalSize(VecLow, &n));
308: if (n > 0) {
309: PetscCall(VecGetArrayRead(VecLow, &v1));
310: if (VecLow != VecHigh) {
311: PetscCall(VecGetArrayRead(VecHigh, &v2));
312: } else {
313: v2 = v1;
314: }
315: if (V != VecLow && V != VecHigh) {
316: PetscCall(VecGetArrayRead(V, &vmiddle));
317: } else if (V == VecLow) {
318: vmiddle = v1;
319: } else {
320: vmiddle = v2;
321: }
323: PetscCall(PetscMalloc1(n, &vm));
325: for (i = 0; i < n; ++i) {
326: if (PetscRealPart(v1[i]) <= PetscRealPart(vmiddle[i]) && PetscRealPart(vmiddle[i]) <= PetscRealPart(v2[i])) {
327: vm[n_vm] = low + i;
328: ++n_vm;
329: }
330: }
332: PetscCall(VecRestoreArrayRead(VecLow, &v1));
333: if (VecLow != VecHigh) PetscCall(VecRestoreArrayRead(VecHigh, &v2));
334: if (V != VecLow && V != VecHigh) PetscCall(VecRestoreArrayRead(V, &vmiddle));
335: }
336: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)V), n_vm, vm, PETSC_OWN_POINTER, S));
337: PetscFunctionReturn(PETSC_SUCCESS);
338: }
340: /*@
341: VecWhichInactive - Creates an `IS` based on a set of vectors
343: Collective
345: Input Parameters:
346: + VecLow - lower bound
347: . V - Vector to compare
348: . D - Direction to compare
349: . VecHigh - higher bound
350: - Strong - indicator for applying strongly inactive test
352: Output Parameter:
353: . S - The index set containing the indices i where the bound is inactive
355: Level: advanced
357: Notes:
358: Creates an index set containing the indices where one of the following holds\:
359: .vb
360: - VecLow(i) < V(i) < VecHigh(i)
361: - VecLow(i) = V(i) and D(i) <= 0 (< 0 when Strong is true)
362: - VecHigh(i) = V(i) and D(i) >= 0 (> 0 when Strong is true)
363: .ve
365: .seealso: `Vec`
366: @*/
367: PetscErrorCode VecWhichInactive(Vec VecLow, Vec V, Vec D, Vec VecHigh, PetscBool Strong, IS *S)
368: {
369: PetscInt i, n_vm = 0;
370: PetscInt n, low, high;
371: PetscInt *vm = NULL;
372: const PetscScalar *v1, *v2, *v, *d;
374: PetscFunctionBegin;
375: if (!VecLow && !VecHigh) {
376: PetscCall(VecGetOwnershipRange(V, &low, &high));
377: PetscCall(ISCreateStride(PetscObjectComm((PetscObject)V), high - low, low, 1, S));
378: PetscFunctionReturn(PETSC_SUCCESS);
379: }
384: PetscCheckSameComm(V, 2, VecLow, 1);
385: PetscCheckSameComm(V, 2, D, 3);
386: PetscCheckSameComm(V, 2, VecHigh, 4);
387: VecCheckSameSize(V, 2, VecLow, 1);
388: VecCheckSameSize(V, 2, D, 3);
389: VecCheckSameSize(V, 2, VecHigh, 4);
391: PetscCall(VecGetOwnershipRange(VecLow, &low, &high));
392: PetscCall(VecGetLocalSize(VecLow, &n));
393: if (n > 0) {
394: PetscCall(VecGetArrayRead(VecLow, &v1));
395: if (VecLow != VecHigh) {
396: PetscCall(VecGetArrayRead(VecHigh, &v2));
397: } else {
398: v2 = v1;
399: }
400: if (V != VecLow && V != VecHigh) {
401: PetscCall(VecGetArrayRead(V, &v));
402: } else if (V == VecLow) {
403: v = v1;
404: } else {
405: v = v2;
406: }
407: if (D != VecLow && D != VecHigh && D != V) {
408: PetscCall(VecGetArrayRead(D, &d));
409: } else if (D == VecLow) {
410: d = v1;
411: } else if (D == VecHigh) {
412: d = v2;
413: } else {
414: d = v;
415: }
417: PetscCall(PetscMalloc1(n, &vm));
419: if (Strong) {
420: for (i = 0; i < n; ++i) {
421: if (PetscRealPart(v1[i]) < PetscRealPart(v[i]) && PetscRealPart(v[i]) < PetscRealPart(v2[i])) {
422: vm[n_vm] = low + i;
423: ++n_vm;
424: } else if (PetscRealPart(v1[i]) == PetscRealPart(v[i]) && PetscRealPart(d[i]) < 0) {
425: vm[n_vm] = low + i;
426: ++n_vm;
427: } else if (PetscRealPart(v2[i]) == PetscRealPart(v[i]) && PetscRealPart(d[i]) > 0) {
428: vm[n_vm] = low + i;
429: ++n_vm;
430: }
431: }
432: } else {
433: for (i = 0; i < n; ++i) {
434: if (PetscRealPart(v1[i]) < PetscRealPart(v[i]) && PetscRealPart(v[i]) < PetscRealPart(v2[i])) {
435: vm[n_vm] = low + i;
436: ++n_vm;
437: } else if (PetscRealPart(v1[i]) == PetscRealPart(v[i]) && PetscRealPart(d[i]) <= 0) {
438: vm[n_vm] = low + i;
439: ++n_vm;
440: } else if (PetscRealPart(v2[i]) == PetscRealPart(v[i]) && PetscRealPart(d[i]) >= 0) {
441: vm[n_vm] = low + i;
442: ++n_vm;
443: }
444: }
445: }
447: PetscCall(VecRestoreArrayRead(VecLow, &v1));
448: if (VecLow != VecHigh) PetscCall(VecRestoreArrayRead(VecHigh, &v2));
449: if (V != VecLow && V != VecHigh) PetscCall(VecRestoreArrayRead(V, &v));
450: if (D != VecLow && D != VecHigh && D != V) PetscCall(VecRestoreArrayRead(D, &d));
451: }
452: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)V), n_vm, vm, PETSC_OWN_POINTER, S));
453: PetscFunctionReturn(PETSC_SUCCESS);
454: }
456: /*@
457: VecISAXPY - Adds a reduced vector to the appropriate elements of a full-space vector.
458: vfull[is[i]] += alpha*vreduced[i]
460: Logically Collective
462: Input Parameters:
463: + vfull - the full-space vector
464: . is - the index set for the reduced space
465: . alpha - the scalar coefficient
466: - vreduced - the reduced-space vector
468: Output Parameter:
469: . vfull - the sum of the full-space vector and reduced-space vector
471: Level: advanced
473: Notes:
474: The index set identifies entries in the global vector.
475: Negative indices are skipped; indices outside the ownership range of `vfull` will raise an error.
477: .seealso: `VecISCopy()`, `VecISSet()`, `VecAXPY()`
478: @*/
479: PetscErrorCode VecISAXPY(Vec vfull, IS is, PetscScalar alpha, Vec vreduced)
480: {
481: PetscInt nfull, nreduced;
482: PetscBool sorted = PETSC_FALSE;
484: PetscFunctionBegin;
487: PetscCheckSameComm(vfull, 1, is, 2);
490: PetscCall(ISGetSize(is, &nfull));
491: if (!nfull) PetscFunctionReturn(PETSC_SUCCESS);
492: PetscCall(VecGetSize(vfull, &nfull));
493: PetscCall(VecGetSize(vreduced, &nreduced));
494: if (nfull == nreduced) PetscCall(ISGetInfo(is, IS_SORTED, IS_GLOBAL, PETSC_TRUE, &sorted));
495: if (sorted) PetscCall(VecAXPY(vfull, alpha, vreduced));
496: else {
497: PetscScalar *y;
498: const PetscScalar *x;
499: PetscInt i, n, m, rstart, rend;
500: const PetscInt *id;
502: PetscCall(VecGetArray(vfull, &y));
503: PetscCall(VecGetArrayRead(vreduced, &x));
504: PetscCall(ISGetIndices(is, &id));
505: PetscCall(ISGetLocalSize(is, &n));
506: PetscCall(VecGetLocalSize(vreduced, &m));
507: PetscCheck(m == n, PETSC_COMM_SELF, PETSC_ERR_SUP, "IS local length %" PetscInt_FMT " not equal to Vec local length %" PetscInt_FMT, n, m);
508: PetscCall(VecGetOwnershipRange(vfull, &rstart, &rend));
509: y = PetscSafePointerPlusOffset(y, -rstart);
510: for (i = 0; i < n; ++i) {
511: if (id[i] < 0) continue;
512: PetscCheck(id[i] >= rstart && id[i] < rend, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only owned values supported");
513: y[id[i]] += alpha * x[i];
514: }
515: y = PetscSafePointerPlusOffset(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: Logically Collective
528: Input Parameters:
529: + vfull - the full-space vector
530: . is - the index set for the reduced space
531: . mode - the direction of copying, `SCATTER_FORWARD` or `SCATTER_REVERSE`
532: - vreduced - the reduced-space vector
534: Output Parameter:
535: . vfull - the sum of the full-space vector and reduced-space vector
537: Level: advanced
539: Notes:
540: The index set identifies entries in the global vector.
541: Negative indices are skipped; indices outside the ownership range of `vfull` will raise an error.
542: .vb
543: mode == SCATTER_FORWARD: vfull[is[i]] = vreduced[i]
544: mode == SCATTER_REVERSE: vreduced[i] = vfull[is[i]]
545: .ve
547: .seealso: `VecISSet()`, `VecISAXPY()`, `VecCopy()`
548: @*/
549: PetscErrorCode VecISCopy(Vec vfull, IS is, ScatterMode mode, Vec vreduced)
550: {
551: PetscInt nfull, nreduced;
552: PetscBool sorted = PETSC_FALSE;
554: PetscFunctionBegin;
557: PetscCheckSameComm(vfull, 1, is, 2);
560: PetscCall(ISGetSize(is, &nfull));
561: if (!nfull) PetscFunctionReturn(PETSC_SUCCESS);
562: PetscCall(VecGetSize(vfull, &nfull));
563: PetscCall(VecGetSize(vreduced, &nreduced));
564: if (nfull == nreduced) PetscCall(ISGetInfo(is, IS_SORTED, IS_GLOBAL, PETSC_TRUE, &sorted));
565: if (sorted) {
566: if (mode == SCATTER_FORWARD) {
567: PetscCall(VecCopy(vreduced, vfull));
568: } else {
569: PetscCall(VecCopy(vfull, vreduced));
570: }
571: } else {
572: const PetscInt *id;
573: PetscInt i, n, m, rstart, rend;
575: PetscCall(ISGetIndices(is, &id));
576: PetscCall(ISGetLocalSize(is, &n));
577: PetscCall(VecGetLocalSize(vreduced, &m));
578: PetscCall(VecGetOwnershipRange(vfull, &rstart, &rend));
579: PetscCheck(m == n, PETSC_COMM_SELF, PETSC_ERR_SUP, "IS local length %" PetscInt_FMT " not equal to Vec local length %" PetscInt_FMT, n, m);
580: if (mode == SCATTER_FORWARD) {
581: PetscScalar *y;
582: const PetscScalar *x;
584: PetscCall(VecGetArray(vfull, &y));
585: PetscCall(VecGetArrayRead(vreduced, &x));
586: y = PetscSafePointerPlusOffset(y, -rstart);
587: for (i = 0; i < n; ++i) {
588: if (id[i] < 0) continue;
589: PetscCheck(id[i] >= rstart && id[i] < rend, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only owned values supported");
590: y[id[i]] = x[i];
591: }
592: y = PetscSafePointerPlusOffset(y, rstart);
593: PetscCall(VecRestoreArrayRead(vreduced, &x));
594: PetscCall(VecRestoreArray(vfull, &y));
595: } else if (mode == SCATTER_REVERSE) {
596: PetscScalar *x;
597: const PetscScalar *y;
599: PetscCall(VecGetArrayRead(vfull, &y));
600: PetscCall(VecGetArray(vreduced, &x));
601: for (i = 0; i < n; ++i) {
602: if (id[i] < 0) continue;
603: PetscCheck(id[i] >= rstart && id[i] < rend, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only owned values supported");
604: x[i] = y[id[i] - rstart];
605: }
606: PetscCall(VecRestoreArray(vreduced, &x));
607: PetscCall(VecRestoreArrayRead(vfull, &y));
608: } else SETERRQ(PetscObjectComm((PetscObject)vfull), PETSC_ERR_ARG_WRONG, "Only forward or reverse modes are legal");
609: PetscCall(ISRestoreIndices(is, &id));
610: }
611: PetscFunctionReturn(PETSC_SUCCESS);
612: }
614: /*@
615: ISComplementVec - Creates the complement of the index set relative to a layout defined by a `Vec`
617: Collective
619: Input Parameters:
620: + S - a PETSc `IS`
621: - V - the reference vector space
623: Output Parameter:
624: . T - the complement of S
626: Level: advanced
628: .seealso: `IS`, `Vec`, `ISCreateGeneral()`
629: @*/
630: PetscErrorCode ISComplementVec(IS S, Vec V, IS *T)
631: {
632: PetscInt start, end;
634: PetscFunctionBegin;
635: PetscCall(VecGetOwnershipRange(V, &start, &end));
636: PetscCall(ISComplement(S, start, end, T));
637: PetscFunctionReturn(PETSC_SUCCESS);
638: }
640: /*@
641: VecISSet - Sets the elements of a vector, specified by an index set, to a constant
643: Logically Collective
645: Input Parameters:
646: + V - the vector
647: . S - index set for the locations in the vector
648: - c - the constant
650: Level: advanced
652: Notes:
653: The index set identifies entries in the global vector.
654: Negative indices are skipped; indices outside the ownership range of V will raise an error.
656: .seealso: `VecISCopy()`, `VecISAXPY()`, `VecISShift()`, `VecSet()`
657: @*/
658: PetscErrorCode VecISSet(Vec V, IS S, PetscScalar c)
659: {
660: PetscInt nloc, low, high, i;
661: const PetscInt *s;
662: PetscScalar *v;
664: PetscFunctionBegin;
667: PetscCheckSameComm(V, 1, S, 2);
668: PetscCall(ISGetSize(S, &nloc));
669: if (nloc) {
670: PetscCall(VecGetOwnershipRange(V, &low, &high));
671: PetscCall(ISGetLocalSize(S, &nloc));
672: PetscCall(ISGetIndices(S, &s));
673: PetscCall(VecGetArray(V, &v));
674: for (i = 0; i < nloc; ++i) {
675: if (s[i] < 0) continue;
676: PetscCheck(s[i] >= low && s[i] < high, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only owned values supported");
677: v[s[i] - low] = c;
678: }
679: PetscCall(ISRestoreIndices(S, &s));
680: PetscCall(VecRestoreArray(V, &v));
681: }
682: PetscFunctionReturn(PETSC_SUCCESS);
683: }
685: /*@
686: VecISShift - Shifts the elements of a vector, specified by an index set, by a constant
688: Logically Collective
690: Input Parameters:
691: + V - the vector
692: . S - index set for the locations in the vector
693: - c - the constant
695: Level: advanced
697: Notes:
698: The index set identifies entries in the global vector.
699: Negative indices are skipped; indices outside the ownership range of V will raise an error.
701: .seealso: `VecISCopy()`, `VecISAXPY()`, `VecISSet()`, `VecShift()`
702: @*/
703: PetscErrorCode VecISShift(Vec V, IS S, PetscScalar c)
704: {
705: PetscInt nloc, low, high, i;
706: const PetscInt *s;
707: PetscScalar *v;
709: PetscFunctionBegin;
712: PetscCheckSameComm(V, 1, S, 2);
713: PetscCall(ISGetSize(S, &nloc));
714: if (nloc) {
715: PetscCall(VecGetOwnershipRange(V, &low, &high));
716: PetscCall(ISGetLocalSize(S, &nloc));
717: PetscCall(ISGetIndices(S, &s));
718: PetscCall(VecGetArray(V, &v));
719: for (i = 0; i < nloc; ++i) {
720: if (s[i] < 0) continue;
721: PetscCheck(s[i] >= low && s[i] < high, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only owned values supported");
722: v[s[i] - low] += c;
723: }
724: PetscCall(ISRestoreIndices(S, &s));
725: PetscCall(VecRestoreArray(V, &v));
726: }
727: PetscFunctionReturn(PETSC_SUCCESS);
728: }
730: /*@
731: VecBoundGradientProjection - Projects vector according to this definition.
732: If XL[i] < X[i] < XU[i], then GP[i] = G[i];
733: If X[i] <= XL[i], then GP[i] = min(G[i],0);
734: If X[i] >= XU[i], then GP[i] = max(G[i],0);
736: Input Parameters:
737: + G - current gradient vector
738: . X - current solution vector with XL[i] <= X[i] <= XU[i]
739: . XL - lower bounds
740: - XU - upper bounds
742: Output Parameter:
743: . GP - gradient projection vector
745: Level: advanced
747: Note:
748: `GP` may be the same vector as `G`
750: For complex numbers only the real part is used in the bounds.
752: .seealso: `Vec`
753: @*/
754: PetscErrorCode VecBoundGradientProjection(Vec G, Vec X, Vec XL, Vec XU, Vec GP)
755: {
756: PetscInt n, i;
757: const PetscScalar *xptr, *xlptr, *xuptr;
758: PetscScalar *gptr, *gpptr;
759: PetscScalar xval, gpval;
761: /* Project variables at the lower and upper bound */
762: PetscFunctionBegin;
768: if (!XL && !XU) {
769: PetscCall(VecCopy(G, GP));
770: PetscFunctionReturn(PETSC_SUCCESS);
771: }
773: PetscCall(VecGetLocalSize(X, &n));
775: PetscCall(VecGetArrayRead(X, &xptr));
776: PetscCall(VecGetArrayRead(XL, &xlptr));
777: PetscCall(VecGetArrayRead(XU, &xuptr));
778: PetscCall(VecGetArrayPair(G, GP, &gptr, &gpptr));
780: for (i = 0; i < n; ++i) {
781: gpval = gptr[i];
782: xval = xptr[i];
783: if (PetscRealPart(gpval) > 0.0 && PetscRealPart(xval) <= PetscRealPart(xlptr[i])) {
784: gpval = 0.0;
785: } else if (PetscRealPart(gpval) < 0.0 && PetscRealPart(xval) >= PetscRealPart(xuptr[i])) {
786: gpval = 0.0;
787: }
788: gpptr[i] = gpval;
789: }
791: PetscCall(VecRestoreArrayRead(X, &xptr));
792: PetscCall(VecRestoreArrayRead(XL, &xlptr));
793: PetscCall(VecRestoreArrayRead(XU, &xuptr));
794: PetscCall(VecRestoreArrayPair(G, GP, &gptr, &gpptr));
795: PetscFunctionReturn(PETSC_SUCCESS);
796: }
798: /*@
799: VecStepMaxBounded - See below
801: Collective
803: Input Parameters:
804: + X - vector with no negative entries
805: . XL - lower bounds
806: . XU - upper bounds
807: - DX - step direction, can have negative, positive or zero entries
809: Output Parameter:
810: . stepmax - minimum value so that X[i] + stepmax*DX[i] <= XL[i] or XU[i] <= X[i] + stepmax*DX[i]
812: Level: intermediate
814: .seealso: `Vec`
815: @*/
816: PetscErrorCode VecStepMaxBounded(Vec X, Vec DX, Vec XL, Vec XU, PetscReal *stepmax)
817: {
818: PetscInt i, nn;
819: const PetscScalar *xx, *dx, *xl, *xu;
820: PetscReal localmax = 0;
822: PetscFunctionBegin;
828: PetscCall(VecGetArrayRead(X, &xx));
829: PetscCall(VecGetArrayRead(XL, &xl));
830: PetscCall(VecGetArrayRead(XU, &xu));
831: PetscCall(VecGetArrayRead(DX, &dx));
832: PetscCall(VecGetLocalSize(X, &nn));
833: for (i = 0; i < nn; i++) {
834: if (PetscRealPart(dx[i]) > 0) {
835: localmax = PetscMax(localmax, PetscRealPart((xu[i] - xx[i]) / dx[i]));
836: } else if (PetscRealPart(dx[i]) < 0) {
837: localmax = PetscMax(localmax, PetscRealPart((xl[i] - xx[i]) / dx[i]));
838: }
839: }
840: PetscCall(VecRestoreArrayRead(X, &xx));
841: PetscCall(VecRestoreArrayRead(XL, &xl));
842: PetscCall(VecRestoreArrayRead(XU, &xu));
843: PetscCall(VecRestoreArrayRead(DX, &dx));
844: PetscCallMPI(MPIU_Allreduce(&localmax, stepmax, 1, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)X)));
845: PetscFunctionReturn(PETSC_SUCCESS);
846: }
848: /*@
849: VecStepBoundInfo - See below
851: Collective
853: Input Parameters:
854: + X - vector with no negative entries
855: . XL - lower bounds
856: . XU - upper bounds
857: - DX - step direction, can have negative, positive or zero entries
859: Output Parameters:
860: + boundmin - (may be `NULL` this it is not computed) maximum value so that XL[i] <= X[i] + boundmax*DX[i] <= XU[i]
861: . wolfemin - (may be `NULL` this it is not computed)
862: - 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]
864: Level: advanced
866: Note:
867: For complex numbers only compares the real part
869: .seealso: `Vec`
870: @*/
871: PetscErrorCode VecStepBoundInfo(Vec X, Vec DX, Vec XL, Vec XU, PetscReal *boundmin, PetscReal *wolfemin, PetscReal *boundmax)
872: {
873: PetscInt n, i;
874: const PetscScalar *x, *xl, *xu, *dx;
875: PetscReal t;
876: PetscReal localmin = PETSC_INFINITY, localwolfemin = PETSC_INFINITY, localmax = -1;
877: MPI_Comm comm;
879: PetscFunctionBegin;
885: PetscCall(VecGetArrayRead(X, &x));
886: PetscCall(VecGetArrayRead(XL, &xl));
887: PetscCall(VecGetArrayRead(XU, &xu));
888: PetscCall(VecGetArrayRead(DX, &dx));
889: PetscCall(VecGetLocalSize(X, &n));
890: for (i = 0; i < n; ++i) {
891: if (PetscRealPart(dx[i]) > 0 && PetscRealPart(xu[i]) < PETSC_INFINITY) {
892: t = PetscRealPart((xu[i] - x[i]) / dx[i]);
893: localmin = PetscMin(t, localmin);
894: if (localmin > 0) localwolfemin = PetscMin(t, localwolfemin);
895: localmax = PetscMax(t, localmax);
896: } else if (PetscRealPart(dx[i]) < 0 && PetscRealPart(xl[i]) > PETSC_NINFINITY) {
897: t = PetscRealPart((xl[i] - x[i]) / dx[i]);
898: localmin = PetscMin(t, localmin);
899: if (localmin > 0) localwolfemin = PetscMin(t, localwolfemin);
900: localmax = PetscMax(t, localmax);
901: }
902: }
904: PetscCall(VecRestoreArrayRead(X, &x));
905: PetscCall(VecRestoreArrayRead(XL, &xl));
906: PetscCall(VecRestoreArrayRead(XU, &xu));
907: PetscCall(VecRestoreArrayRead(DX, &dx));
908: PetscCall(PetscObjectGetComm((PetscObject)X, &comm));
910: if (boundmin) {
911: PetscCallMPI(MPIU_Allreduce(&localmin, boundmin, 1, MPIU_REAL, MPIU_MIN, comm));
912: PetscCall(PetscInfo(X, "Step Bound Info: Closest Bound: %20.19e\n", (double)*boundmin));
913: }
914: if (wolfemin) {
915: PetscCallMPI(MPIU_Allreduce(&localwolfemin, wolfemin, 1, MPIU_REAL, MPIU_MIN, comm));
916: PetscCall(PetscInfo(X, "Step Bound Info: Wolfe: %20.19e\n", (double)*wolfemin));
917: }
918: if (boundmax) {
919: PetscCallMPI(MPIU_Allreduce(&localmax, boundmax, 1, MPIU_REAL, MPIU_MAX, comm));
920: if (*boundmax < 0) *boundmax = PETSC_INFINITY;
921: PetscCall(PetscInfo(X, "Step Bound Info: Max: %20.19e\n", (double)*boundmax));
922: }
923: PetscFunctionReturn(PETSC_SUCCESS);
924: }
926: /*@
927: VecStepMax - Returns the largest value so that x[i] + step*DX[i] >= 0 for all i
929: Collective
931: Input Parameters:
932: + X - vector with no negative entries
933: - DX - a step direction, can have negative, positive or zero entries
935: Output Parameter:
936: . step - largest value such that x[i] + step*DX[i] >= 0 for all i
938: Level: advanced
940: Note:
941: For complex numbers only compares the real part
943: .seealso: `Vec`
944: @*/
945: PetscErrorCode VecStepMax(Vec X, Vec DX, PetscReal *step)
946: {
947: PetscInt i, nn;
948: PetscReal stepmax = PETSC_INFINITY;
949: const PetscScalar *xx, *dx;
951: PetscFunctionBegin;
955: PetscCall(VecGetLocalSize(X, &nn));
956: PetscCall(VecGetArrayRead(X, &xx));
957: PetscCall(VecGetArrayRead(DX, &dx));
958: for (i = 0; i < nn; ++i) {
959: PetscCheck(PetscRealPart(xx[i]) >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Vector must be positive");
960: if (PetscRealPart(dx[i]) < 0) stepmax = PetscMin(stepmax, PetscRealPart(-xx[i] / dx[i]));
961: }
962: PetscCall(VecRestoreArrayRead(X, &xx));
963: PetscCall(VecRestoreArrayRead(DX, &dx));
964: PetscCallMPI(MPIU_Allreduce(&stepmax, step, 1, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)X)));
965: PetscFunctionReturn(PETSC_SUCCESS);
966: }
968: /*@
969: VecPow - Replaces each component of a vector by $ x_i^p $
971: Logically Collective
973: Input Parameters:
974: + v - the vector
975: - p - the exponent to use on each element
977: Level: intermediate
979: Note:
980: This handles negative values, Inf, and Nan in the expected IEEE floating pointing manner. For example, the square root of a negative real number is Nan
981: and 1/0.0 is Inf.
983: .seealso: `Vec`
984: @*/
985: PetscErrorCode VecPow(Vec v, PetscScalar p)
986: {
987: PetscInt n, i;
988: PetscScalar *v1;
990: PetscFunctionBegin;
993: PetscCall(VecGetArray(v, &v1));
994: PetscCall(VecGetLocalSize(v, &n));
996: if (1.0 == p) {
997: } else if (-1.0 == p) {
998: for (i = 0; i < n; ++i) v1[i] = 1.0 / v1[i];
999: } else if (0.0 == p) {
1000: for (i = 0; i < n; ++i) {
1001: /* Not-a-number left alone
1002: Infinity set to one */
1003: if (!PetscIsNanScalar(v1[i])) v1[i] = 1.0;
1004: }
1005: } else if (0.5 == p) {
1006: for (i = 0; i < n; ++i) { v1[i] = PetscSqrtScalar(v1[i]); }
1007: } else if (-0.5 == p) {
1008: for (i = 0; i < n; ++i) { v1[i] = 1.0 / PetscSqrtScalar(v1[i]); }
1009: } else if (2.0 == p) {
1010: for (i = 0; i < n; ++i) v1[i] *= v1[i];
1011: } else if (-2.0 == p) {
1012: for (i = 0; i < n; ++i) v1[i] = 1.0 / (v1[i] * v1[i]);
1013: } else {
1014: for (i = 0; i < n; ++i) { v1[i] = PetscPowScalar(v1[i], p); }
1015: }
1016: PetscCall(VecRestoreArray(v, &v1));
1017: PetscFunctionReturn(PETSC_SUCCESS);
1018: }
1020: /*@
1021: VecMedian - Computes the componentwise median of three vectors
1022: and stores the result in this vector. Used primarily for projecting
1023: a vector within upper and lower bounds.
1025: Logically Collective
1027: Input Parameters:
1028: + Vec1 - The first vector
1029: . Vec2 - The second vector
1030: - Vec3 - The third vector
1032: Output Parameter:
1033: . VMedian - The median vector (this can be any one of the input vectors)
1035: Level: advanced
1037: .seealso: `Vec`
1038: @*/
1039: PetscErrorCode VecMedian(Vec Vec1, Vec Vec2, Vec Vec3, Vec VMedian)
1040: {
1041: PetscInt i, n, low1, high1;
1042: const PetscScalar *v1, *v2, *v3;
1043: PetscScalar *vmed;
1045: PetscFunctionBegin;
1051: if (!Vec1 && !Vec3) {
1052: PetscCall(VecCopy(Vec2, VMedian));
1053: PetscFunctionReturn(PETSC_SUCCESS);
1054: }
1055: if (Vec1 == Vec2 || Vec1 == Vec3) {
1056: PetscCall(VecCopy(Vec1, VMedian));
1057: PetscFunctionReturn(PETSC_SUCCESS);
1058: }
1059: if (Vec2 == Vec3) {
1060: PetscCall(VecCopy(Vec2, VMedian));
1061: PetscFunctionReturn(PETSC_SUCCESS);
1062: }
1064: /* Assert that Vec1 != Vec2 and Vec2 != Vec3 */
1069: PetscCheckSameType(Vec1, 1, Vec2, 2);
1070: PetscCheckSameType(Vec1, 1, Vec3, 3);
1071: PetscCheckSameType(Vec1, 1, VMedian, 4);
1072: PetscCheckSameComm(Vec1, 1, Vec2, 2);
1073: PetscCheckSameComm(Vec1, 1, Vec3, 3);
1074: PetscCheckSameComm(Vec1, 1, VMedian, 4);
1075: VecCheckSameSize(Vec1, 1, Vec2, 2);
1076: VecCheckSameSize(Vec1, 1, Vec3, 3);
1077: VecCheckSameSize(Vec1, 1, VMedian, 4);
1079: PetscCall(VecGetOwnershipRange(Vec1, &low1, &high1));
1080: PetscCall(VecGetLocalSize(Vec1, &n));
1081: if (n > 0) {
1082: PetscCall(VecGetArray(VMedian, &vmed));
1083: if (Vec1 != VMedian) {
1084: PetscCall(VecGetArrayRead(Vec1, &v1));
1085: } else {
1086: v1 = vmed;
1087: }
1088: if (Vec2 != VMedian) {
1089: PetscCall(VecGetArrayRead(Vec2, &v2));
1090: } else {
1091: v2 = vmed;
1092: }
1093: if (Vec3 != VMedian) {
1094: PetscCall(VecGetArrayRead(Vec3, &v3));
1095: } else {
1096: v3 = vmed;
1097: }
1099: 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])));
1101: PetscCall(VecRestoreArray(VMedian, &vmed));
1102: if (VMedian != Vec1) PetscCall(VecRestoreArrayRead(Vec1, &v1));
1103: if (VMedian != Vec2) PetscCall(VecRestoreArrayRead(Vec2, &v2));
1104: if (VMedian != Vec3) PetscCall(VecRestoreArrayRead(Vec3, &v3));
1105: }
1106: PetscFunctionReturn(PETSC_SUCCESS);
1107: }