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, &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:   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, &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:   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: }