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, &lt));

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, &gt));

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;
488:   PetscCall(VecGetSize(vfull, &nfull));
489:   PetscCall(VecGetSize(vreduced, &nreduced));
490:   if (nfull == nreduced) PetscCall(ISGetInfo(is, IS_SORTED, IS_GLOBAL, PETSC_TRUE, &sorted));
491:   if (sorted) PetscCall(VecAXPY(vfull, alpha, vreduced));
492:   else {
493:     PetscScalar       *y;
494:     const PetscScalar *x;
495:     PetscInt           i, n, m, rstart, rend;
496:     const PetscInt    *id;

498:     PetscCall(VecGetArray(vfull, &y));
499:     PetscCall(VecGetArrayRead(vreduced, &x));
500:     PetscCall(ISGetIndices(is, &id));
501:     PetscCall(ISGetLocalSize(is, &n));
502:     PetscCall(VecGetLocalSize(vreduced, &m));
503:     PetscCheck(m == n, PETSC_COMM_SELF, PETSC_ERR_SUP, "IS local length %" PetscInt_FMT " not equal to Vec local length %" PetscInt_FMT, n, m);
504:     PetscCall(VecGetOwnershipRange(vfull, &rstart, &rend));
505:     y = PetscSafePointerPlusOffset(y, -rstart);
506:     for (i = 0; i < n; ++i) {
507:       if (id[i] < 0) continue;
508:       PetscCheck(id[i] >= rstart && id[i] < rend, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only owned values supported");
509:       y[id[i]] += alpha * x[i];
510:     }
511:     y = PetscSafePointerPlusOffset(y, rstart);
512:     PetscCall(ISRestoreIndices(is, &id));
513:     PetscCall(VecRestoreArray(vfull, &y));
514:     PetscCall(VecRestoreArrayRead(vreduced, &x));
515:   }
516:   PetscFunctionReturn(PETSC_SUCCESS);
517: }

519: /*@
520:   VecISCopy - Copies between a reduced vector and the appropriate elements of a full-space vector.

522:   Logically Collective

524:   Input Parameters:
525: + vfull    - the full-space vector
526: . is       - the index set for the reduced space
527: . mode     - the direction of copying, `SCATTER_FORWARD` or `SCATTER_REVERSE`
528: - vreduced - the reduced-space vector

530:   Output Parameter:
531: . vfull - the sum of the full-space vector and reduced-space vector

533:   Level: advanced

535:   Notes:
536:   The index set identifies entries in the global vector.
537:   Negative indices are skipped; indices outside the ownership range of `vfull` will raise an error.
538: .vb
539:     mode == SCATTER_FORWARD: vfull[is[i]] = vreduced[i]
540:     mode == SCATTER_REVERSE: vreduced[i] = vfull[is[i]]
541: .ve

543: .seealso: `VecISSet()`, `VecISAXPY()`, `VecCopy()`
544: @*/
545: PetscErrorCode VecISCopy(Vec vfull, IS is, ScatterMode mode, Vec vreduced)
546: {
547:   PetscInt  nfull, nreduced;
548:   PetscBool sorted = PETSC_FALSE;

550:   PetscFunctionBegin;
555:   PetscCall(VecGetSize(vfull, &nfull));
556:   PetscCall(VecGetSize(vreduced, &nreduced));
557:   if (nfull == nreduced) PetscCall(ISGetInfo(is, IS_SORTED, IS_GLOBAL, PETSC_TRUE, &sorted));
558:   if (sorted) {
559:     if (mode == SCATTER_FORWARD) {
560:       PetscCall(VecCopy(vreduced, vfull));
561:     } else {
562:       PetscCall(VecCopy(vfull, vreduced));
563:     }
564:   } else {
565:     const PetscInt *id;
566:     PetscInt        i, n, m, rstart, rend;

568:     PetscCall(ISGetIndices(is, &id));
569:     PetscCall(ISGetLocalSize(is, &n));
570:     PetscCall(VecGetLocalSize(vreduced, &m));
571:     PetscCall(VecGetOwnershipRange(vfull, &rstart, &rend));
572:     PetscCheck(m == n, PETSC_COMM_SELF, PETSC_ERR_SUP, "IS local length %" PetscInt_FMT " not equal to Vec local length %" PetscInt_FMT, n, m);
573:     if (mode == SCATTER_FORWARD) {
574:       PetscScalar       *y;
575:       const PetscScalar *x;

577:       PetscCall(VecGetArray(vfull, &y));
578:       PetscCall(VecGetArrayRead(vreduced, &x));
579:       y = PetscSafePointerPlusOffset(y, -rstart);
580:       for (i = 0; i < n; ++i) {
581:         if (id[i] < 0) continue;
582:         PetscCheck(id[i] >= rstart && id[i] < rend, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only owned values supported");
583:         y[id[i]] = x[i];
584:       }
585:       y = PetscSafePointerPlusOffset(y, rstart);
586:       PetscCall(VecRestoreArrayRead(vreduced, &x));
587:       PetscCall(VecRestoreArray(vfull, &y));
588:     } else if (mode == SCATTER_REVERSE) {
589:       PetscScalar       *x;
590:       const PetscScalar *y;

592:       PetscCall(VecGetArrayRead(vfull, &y));
593:       PetscCall(VecGetArray(vreduced, &x));
594:       for (i = 0; i < n; ++i) {
595:         if (id[i] < 0) continue;
596:         PetscCheck(id[i] >= rstart && id[i] < rend, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only owned values supported");
597:         x[i] = y[id[i] - rstart];
598:       }
599:       PetscCall(VecRestoreArray(vreduced, &x));
600:       PetscCall(VecRestoreArrayRead(vfull, &y));
601:     } else SETERRQ(PetscObjectComm((PetscObject)vfull), PETSC_ERR_ARG_WRONG, "Only forward or reverse modes are legal");
602:     PetscCall(ISRestoreIndices(is, &id));
603:   }
604:   PetscFunctionReturn(PETSC_SUCCESS);
605: }

607: /*@
608:   ISComplementVec - Creates the complement of the index set relative to a layout defined by a `Vec`

610:   Collective

612:   Input Parameters:
613: + S - a PETSc `IS`
614: - V - the reference vector space

616:   Output Parameter:
617: . T - the complement of S

619:   Level: advanced

621: .seealso: `IS`, `Vec`, `ISCreateGeneral()`
622: @*/
623: PetscErrorCode ISComplementVec(IS S, Vec V, IS *T)
624: {
625:   PetscInt start, end;

627:   PetscFunctionBegin;
628:   PetscCall(VecGetOwnershipRange(V, &start, &end));
629:   PetscCall(ISComplement(S, start, end, T));
630:   PetscFunctionReturn(PETSC_SUCCESS);
631: }

633: /*@
634:   VecISSet - Sets the elements of a vector, specified by an index set, to a constant

636:   Logically Collective

638:   Input Parameters:
639: + V - the vector
640: . S - index set for the locations in the vector
641: - c - the constant

643:   Level: advanced

645:   Notes:
646:   The index set identifies entries in the global vector.
647:   Negative indices are skipped; indices outside the ownership range of V will raise an error.

649: .seealso: `VecISCopy()`, `VecISAXPY()`, `VecISShift()`, `VecSet()`
650: @*/
651: PetscErrorCode VecISSet(Vec V, IS S, PetscScalar c)
652: {
653:   PetscInt        nloc, low, high, i;
654:   const PetscInt *s;
655:   PetscScalar    *v;

657:   PetscFunctionBegin;
661:   PetscCall(VecGetOwnershipRange(V, &low, &high));
662:   PetscCall(ISGetLocalSize(S, &nloc));
663:   PetscCall(ISGetIndices(S, &s));
664:   PetscCall(VecGetArray(V, &v));
665:   for (i = 0; i < nloc; ++i) {
666:     if (s[i] < 0) continue;
667:     PetscCheck(s[i] >= low && s[i] < high, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only owned values supported");
668:     v[s[i] - low] = c;
669:   }
670:   PetscCall(ISRestoreIndices(S, &s));
671:   PetscCall(VecRestoreArray(V, &v));
672:   PetscFunctionReturn(PETSC_SUCCESS);
673: }

675: /*@
676:   VecISShift - Shifts the elements of a vector, specified by an index set, by a constant

678:   Logically Collective

680:   Input Parameters:
681: + V - the vector
682: . S - index set for the locations in the vector
683: - c - the constant

685:   Level: advanced

687:   Notes:
688:   The index set identifies entries in the global vector.
689:   Negative indices are skipped; indices outside the ownership range of V will raise an error.

691: .seealso: `VecISCopy()`, `VecISAXPY()`, `VecISSet()`, `VecShift()`
692: @*/
693: PetscErrorCode VecISShift(Vec V, IS S, PetscScalar c)
694: {
695:   PetscInt        nloc, low, high, i;
696:   const PetscInt *s;
697:   PetscScalar    *v;

699:   PetscFunctionBegin;
703:   PetscCall(VecGetOwnershipRange(V, &low, &high));
704:   PetscCall(ISGetLocalSize(S, &nloc));
705:   PetscCall(ISGetIndices(S, &s));
706:   PetscCall(VecGetArray(V, &v));
707:   for (i = 0; i < nloc; ++i) {
708:     if (s[i] < 0) continue;
709:     PetscCheck(s[i] >= low && s[i] < high, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only owned values supported");
710:     v[s[i] - low] += c;
711:   }
712:   PetscCall(ISRestoreIndices(S, &s));
713:   PetscCall(VecRestoreArray(V, &v));
714:   PetscFunctionReturn(PETSC_SUCCESS);
715: }

717: #if !defined(PETSC_USE_COMPLEX)
718: /*@C
719:   VecBoundGradientProjection - Projects  vector according to this definition.
720:   If XL[i] < X[i] < XU[i], then GP[i] = G[i];
721:   If X[i] <= XL[i], then GP[i] = min(G[i],0);
722:   If X[i] >= XU[i], then GP[i] = max(G[i],0);

724:   Input Parameters:
725: + G  - current gradient vector
726: . X  - current solution vector with XL[i] <= X[i] <= XU[i]
727: . XL - lower bounds
728: - XU - upper bounds

730:   Output Parameter:
731: . GP - gradient projection vector

733:   Level: advanced

735:   Note:
736:   `GP` may be the same vector as `G`

738: .seealso: `Vec`
739: @*/
740: PetscErrorCode VecBoundGradientProjection(Vec G, Vec X, Vec XL, Vec XU, Vec GP)
741: {
742:   PetscInt         n, i;
743:   const PetscReal *xptr, *xlptr, *xuptr;
744:   PetscReal       *gptr, *gpptr;
745:   PetscReal        xval, gpval;

747:   /* Project variables at the lower and upper bound */
748:   PetscFunctionBegin;
754:   if (!XL && !XU) {
755:     PetscCall(VecCopy(G, GP));
756:     PetscFunctionReturn(PETSC_SUCCESS);
757:   }

759:   PetscCall(VecGetLocalSize(X, &n));

761:   PetscCall(VecGetArrayRead(X, &xptr));
762:   PetscCall(VecGetArrayRead(XL, &xlptr));
763:   PetscCall(VecGetArrayRead(XU, &xuptr));
764:   PetscCall(VecGetArrayPair(G, GP, &gptr, &gpptr));

766:   for (i = 0; i < n; ++i) {
767:     gpval = gptr[i];
768:     xval  = xptr[i];
769:     if (gpval > 0.0 && xval <= xlptr[i]) {
770:       gpval = 0.0;
771:     } else if (gpval < 0.0 && xval >= xuptr[i]) {
772:       gpval = 0.0;
773:     }
774:     gpptr[i] = gpval;
775:   }

777:   PetscCall(VecRestoreArrayRead(X, &xptr));
778:   PetscCall(VecRestoreArrayRead(XL, &xlptr));
779:   PetscCall(VecRestoreArrayRead(XU, &xuptr));
780:   PetscCall(VecRestoreArrayPair(G, GP, &gptr, &gpptr));
781:   PetscFunctionReturn(PETSC_SUCCESS);
782: }
783: #endif

785: /*@
786:   VecStepMaxBounded - See below

788:   Collective

790:   Input Parameters:
791: + X  - vector with no negative entries
792: . XL - lower bounds
793: . XU - upper bounds
794: - DX - step direction, can have negative, positive or zero entries

796:   Output Parameter:
797: . stepmax - minimum value so that X[i] + stepmax*DX[i] <= XL[i]  or  XU[i] <= X[i] + stepmax*DX[i]

799:   Level: intermediate

801: .seealso: `Vec`
802: @*/
803: PetscErrorCode VecStepMaxBounded(Vec X, Vec DX, Vec XL, Vec XU, PetscReal *stepmax)
804: {
805:   PetscInt           i, nn;
806:   const PetscScalar *xx, *dx, *xl, *xu;
807:   PetscReal          localmax = 0;

809:   PetscFunctionBegin;

815:   PetscCall(VecGetArrayRead(X, &xx));
816:   PetscCall(VecGetArrayRead(XL, &xl));
817:   PetscCall(VecGetArrayRead(XU, &xu));
818:   PetscCall(VecGetArrayRead(DX, &dx));
819:   PetscCall(VecGetLocalSize(X, &nn));
820:   for (i = 0; i < nn; i++) {
821:     if (PetscRealPart(dx[i]) > 0) {
822:       localmax = PetscMax(localmax, PetscRealPart((xu[i] - xx[i]) / dx[i]));
823:     } else if (PetscRealPart(dx[i]) < 0) {
824:       localmax = PetscMax(localmax, PetscRealPart((xl[i] - xx[i]) / dx[i]));
825:     }
826:   }
827:   PetscCall(VecRestoreArrayRead(X, &xx));
828:   PetscCall(VecRestoreArrayRead(XL, &xl));
829:   PetscCall(VecRestoreArrayRead(XU, &xu));
830:   PetscCall(VecRestoreArrayRead(DX, &dx));
831:   PetscCall(MPIU_Allreduce(&localmax, stepmax, 1, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)X)));
832:   PetscFunctionReturn(PETSC_SUCCESS);
833: }

835: /*@
836:   VecStepBoundInfo - See below

838:   Collective

840:   Input Parameters:
841: + X  - vector with no negative entries
842: . XL - lower bounds
843: . XU - upper bounds
844: - DX - step direction, can have negative, positive or zero entries

846:   Output Parameters:
847: + boundmin - (may be `NULL` this it is not computed) maximum value so that   XL[i] <= X[i] + boundmax*DX[i] <= XU[i]
848: . wolfemin - (may be `NULL` this it is not computed)
849: - 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]

851:   Level: advanced

853:   Note:
854:   For complex numbers only compares the real part

856: .seealso: `Vec`
857: @*/
858: PetscErrorCode VecStepBoundInfo(Vec X, Vec DX, Vec XL, Vec XU, PetscReal *boundmin, PetscReal *wolfemin, PetscReal *boundmax)
859: {
860:   PetscInt           n, i;
861:   const PetscScalar *x, *xl, *xu, *dx;
862:   PetscReal          t;
863:   PetscReal          localmin = PETSC_INFINITY, localwolfemin = PETSC_INFINITY, localmax = -1;
864:   MPI_Comm           comm;

866:   PetscFunctionBegin;

872:   PetscCall(VecGetArrayRead(X, &x));
873:   PetscCall(VecGetArrayRead(XL, &xl));
874:   PetscCall(VecGetArrayRead(XU, &xu));
875:   PetscCall(VecGetArrayRead(DX, &dx));
876:   PetscCall(VecGetLocalSize(X, &n));
877:   for (i = 0; i < n; ++i) {
878:     if (PetscRealPart(dx[i]) > 0 && PetscRealPart(xu[i]) < PETSC_INFINITY) {
879:       t        = PetscRealPart((xu[i] - x[i]) / dx[i]);
880:       localmin = PetscMin(t, localmin);
881:       if (localmin > 0) localwolfemin = PetscMin(t, localwolfemin);
882:       localmax = PetscMax(t, localmax);
883:     } else if (PetscRealPart(dx[i]) < 0 && PetscRealPart(xl[i]) > PETSC_NINFINITY) {
884:       t        = PetscRealPart((xl[i] - x[i]) / dx[i]);
885:       localmin = PetscMin(t, localmin);
886:       if (localmin > 0) localwolfemin = PetscMin(t, localwolfemin);
887:       localmax = PetscMax(t, localmax);
888:     }
889:   }

891:   PetscCall(VecRestoreArrayRead(X, &x));
892:   PetscCall(VecRestoreArrayRead(XL, &xl));
893:   PetscCall(VecRestoreArrayRead(XU, &xu));
894:   PetscCall(VecRestoreArrayRead(DX, &dx));
895:   PetscCall(PetscObjectGetComm((PetscObject)X, &comm));

897:   if (boundmin) {
898:     PetscCall(MPIU_Allreduce(&localmin, boundmin, 1, MPIU_REAL, MPIU_MIN, comm));
899:     PetscCall(PetscInfo(X, "Step Bound Info: Closest Bound: %20.19e\n", (double)*boundmin));
900:   }
901:   if (wolfemin) {
902:     PetscCall(MPIU_Allreduce(&localwolfemin, wolfemin, 1, MPIU_REAL, MPIU_MIN, comm));
903:     PetscCall(PetscInfo(X, "Step Bound Info: Wolfe: %20.19e\n", (double)*wolfemin));
904:   }
905:   if (boundmax) {
906:     PetscCall(MPIU_Allreduce(&localmax, boundmax, 1, MPIU_REAL, MPIU_MAX, comm));
907:     if (*boundmax < 0) *boundmax = PETSC_INFINITY;
908:     PetscCall(PetscInfo(X, "Step Bound Info: Max: %20.19e\n", (double)*boundmax));
909:   }
910:   PetscFunctionReturn(PETSC_SUCCESS);
911: }

913: /*@
914:   VecStepMax - Returns the largest value so that x[i] + step*DX[i] >= 0 for all i

916:   Collective

918:   Input Parameters:
919: + X  - vector with no negative entries
920: - DX - a step direction, can have negative, positive or zero entries

922:   Output Parameter:
923: . step - largest value such that x[i] + step*DX[i] >= 0 for all i

925:   Level: advanced

927:   Note:
928:   For complex numbers only compares the real part

930: .seealso: `Vec`
931: @*/
932: PetscErrorCode VecStepMax(Vec X, Vec DX, PetscReal *step)
933: {
934:   PetscInt           i, nn;
935:   PetscReal          stepmax = PETSC_INFINITY;
936:   const PetscScalar *xx, *dx;

938:   PetscFunctionBegin;

942:   PetscCall(VecGetLocalSize(X, &nn));
943:   PetscCall(VecGetArrayRead(X, &xx));
944:   PetscCall(VecGetArrayRead(DX, &dx));
945:   for (i = 0; i < nn; ++i) {
946:     PetscCheck(PetscRealPart(xx[i]) >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Vector must be positive");
947:     if (PetscRealPart(dx[i]) < 0) stepmax = PetscMin(stepmax, PetscRealPart(-xx[i] / dx[i]));
948:   }
949:   PetscCall(VecRestoreArrayRead(X, &xx));
950:   PetscCall(VecRestoreArrayRead(DX, &dx));
951:   PetscCall(MPIU_Allreduce(&stepmax, step, 1, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)X)));
952:   PetscFunctionReturn(PETSC_SUCCESS);
953: }

955: /*@
956:   VecPow - Replaces each component of a vector by x_i^p

958:   Logically Collective

960:   Input Parameters:
961: + v - the vector
962: - p - the exponent to use on each element

964:   Level: intermediate

966: .seealso: `Vec`
967: @*/
968: PetscErrorCode VecPow(Vec v, PetscScalar p)
969: {
970:   PetscInt     n, i;
971:   PetscScalar *v1;

973:   PetscFunctionBegin;

976:   PetscCall(VecGetArray(v, &v1));
977:   PetscCall(VecGetLocalSize(v, &n));

979:   if (1.0 == p) {
980:   } else if (-1.0 == p) {
981:     for (i = 0; i < n; ++i) v1[i] = 1.0 / v1[i];
982:   } else if (0.0 == p) {
983:     for (i = 0; i < n; ++i) {
984:       /*  Not-a-number left alone
985:           Infinity set to one  */
986:       if (v1[i] == v1[i]) v1[i] = 1.0;
987:     }
988:   } else if (0.5 == p) {
989:     for (i = 0; i < n; ++i) {
990:       if (PetscRealPart(v1[i]) >= 0) {
991:         v1[i] = PetscSqrtScalar(v1[i]);
992:       } else {
993:         v1[i] = PETSC_INFINITY;
994:       }
995:     }
996:   } else if (-0.5 == p) {
997:     for (i = 0; i < n; ++i) {
998:       if (PetscRealPart(v1[i]) >= 0) {
999:         v1[i] = 1.0 / PetscSqrtScalar(v1[i]);
1000:       } else {
1001:         v1[i] = PETSC_INFINITY;
1002:       }
1003:     }
1004:   } else if (2.0 == p) {
1005:     for (i = 0; i < n; ++i) v1[i] *= v1[i];
1006:   } else if (-2.0 == p) {
1007:     for (i = 0; i < n; ++i) v1[i] = 1.0 / (v1[i] * v1[i]);
1008:   } else {
1009:     for (i = 0; i < n; ++i) {
1010:       if (PetscRealPart(v1[i]) >= 0) {
1011:         v1[i] = PetscPowScalar(v1[i], p);
1012:       } else {
1013:         v1[i] = PETSC_INFINITY;
1014:       }
1015:     }
1016:   }
1017:   PetscCall(VecRestoreArray(v, &v1));
1018:   PetscFunctionReturn(PETSC_SUCCESS);
1019: }

1021: /*@
1022:   VecMedian - Computes the componentwise median of three vectors
1023:   and stores the result in this vector.  Used primarily for projecting
1024:   a vector within upper and lower bounds.

1026:   Logically Collective

1028:   Input Parameters:
1029: + Vec1 - The first vector
1030: . Vec2 - The second vector
1031: - Vec3 - The third vector

1033:   Output Parameter:
1034: . VMedian - The median vector (this can be any one of the input vectors)

1036:   Level: advanced

1038: .seealso: `Vec`
1039: @*/
1040: PetscErrorCode VecMedian(Vec Vec1, Vec Vec2, Vec Vec3, Vec VMedian)
1041: {
1042:   PetscInt           i, n, low1, high1;
1043:   const PetscScalar *v1, *v2, *v3;
1044:   PetscScalar       *vmed;

1046:   PetscFunctionBegin;

1052:   if (!Vec1 && !Vec3) {
1053:     PetscCall(VecCopy(Vec2, VMedian));
1054:     PetscFunctionReturn(PETSC_SUCCESS);
1055:   }
1056:   if (Vec1 == Vec2 || Vec1 == Vec3) {
1057:     PetscCall(VecCopy(Vec1, VMedian));
1058:     PetscFunctionReturn(PETSC_SUCCESS);
1059:   }
1060:   if (Vec2 == Vec3) {
1061:     PetscCall(VecCopy(Vec2, VMedian));
1062:     PetscFunctionReturn(PETSC_SUCCESS);
1063:   }

1065:   /* Assert that Vec1 != Vec2 and Vec2 != Vec3 */
1070:   PetscCheckSameType(Vec1, 1, Vec2, 2);
1071:   PetscCheckSameType(Vec1, 1, Vec3, 3);
1072:   PetscCheckSameType(Vec1, 1, VMedian, 4);
1073:   PetscCheckSameComm(Vec1, 1, Vec2, 2);
1074:   PetscCheckSameComm(Vec1, 1, Vec3, 3);
1075:   PetscCheckSameComm(Vec1, 1, VMedian, 4);
1076:   VecCheckSameSize(Vec1, 1, Vec2, 2);
1077:   VecCheckSameSize(Vec1, 1, Vec3, 3);
1078:   VecCheckSameSize(Vec1, 1, VMedian, 4);

1080:   PetscCall(VecGetOwnershipRange(Vec1, &low1, &high1));
1081:   PetscCall(VecGetLocalSize(Vec1, &n));
1082:   if (n > 0) {
1083:     PetscCall(VecGetArray(VMedian, &vmed));
1084:     if (Vec1 != VMedian) {
1085:       PetscCall(VecGetArrayRead(Vec1, &v1));
1086:     } else {
1087:       v1 = vmed;
1088:     }
1089:     if (Vec2 != VMedian) {
1090:       PetscCall(VecGetArrayRead(Vec2, &v2));
1091:     } else {
1092:       v2 = vmed;
1093:     }
1094:     if (Vec3 != VMedian) {
1095:       PetscCall(VecGetArrayRead(Vec3, &v3));
1096:     } else {
1097:       v3 = vmed;
1098:     }

1100:     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])));

1102:     PetscCall(VecRestoreArray(VMedian, &vmed));
1103:     if (VMedian != Vec1) PetscCall(VecRestoreArrayRead(Vec1, &v1));
1104:     if (VMedian != Vec2) PetscCall(VecRestoreArrayRead(Vec2, &v2));
1105:     if (VMedian != Vec3) PetscCall(VecRestoreArrayRead(Vec3, &v3));
1106:   }
1107:   PetscFunctionReturn(PETSC_SUCCESS);
1108: }