Actual source code: vector.c

  1: /*
  2:      Provides the interface functions for vector operations that do NOT have PetscScalar/PetscReal in the signature
  3:    These are the vector functions the user calls.
  4: */
  5: #include <petsc/private/vecimpl.h>
  6: #include <petsc/private/deviceimpl.h>

  8: /* Logging support */
  9: PetscClassId  VEC_CLASSID;
 10: PetscLogEvent VEC_View, VEC_Max, VEC_Min, VEC_Dot, VEC_MDot, VEC_TDot;
 11: PetscLogEvent VEC_Norm, VEC_Normalize, VEC_Scale, VEC_Shift, VEC_Copy, VEC_Set, VEC_AXPY, VEC_AYPX, VEC_WAXPY;
 12: PetscLogEvent VEC_MTDot, VEC_MAXPY, VEC_Swap, VEC_AssemblyBegin, VEC_ScatterBegin, VEC_ScatterEnd;
 13: PetscLogEvent VEC_AssemblyEnd, VEC_PointwiseMult, VEC_SetValues, VEC_Load, VEC_SetPreallocateCOO, VEC_SetValuesCOO;
 14: PetscLogEvent VEC_SetRandom, VEC_ReduceArithmetic, VEC_ReduceCommunication, VEC_ReduceBegin, VEC_ReduceEnd, VEC_Ops;
 15: PetscLogEvent VEC_DotNorm2, VEC_AXPBYPCZ;
 16: PetscLogEvent VEC_ViennaCLCopyFromGPU, VEC_ViennaCLCopyToGPU;
 17: PetscLogEvent VEC_CUDACopyFromGPU, VEC_CUDACopyToGPU;
 18: PetscLogEvent VEC_HIPCopyFromGPU, VEC_HIPCopyToGPU;

 20: /*@
 21:   VecStashGetInfo - Gets how many values are currently in the vector stash, i.e. need
 22:   to be communicated to other processors during the `VecAssemblyBegin()`/`VecAssemblyEnd()` process

 24:   Not Collective

 26:   Input Parameter:
 27: . vec - the vector

 29:   Output Parameters:
 30: + nstash    - the size of the stash
 31: . reallocs  - the number of additional mallocs incurred in building the stash
 32: . bnstash   - the size of the block stash
 33: - breallocs - the number of additional mallocs incurred in building the block stash (from `VecSetValuesBlocked()`)

 35:   Level: advanced

 37: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecStashSetInitialSize()`, `VecStashView()`
 38: @*/
 39: PetscErrorCode VecStashGetInfo(Vec vec, PetscInt *nstash, PetscInt *reallocs, PetscInt *bnstash, PetscInt *breallocs)
 40: {
 41:   PetscFunctionBegin;
 42:   PetscCall(VecStashGetInfo_Private(&vec->stash, nstash, reallocs));
 43:   PetscCall(VecStashGetInfo_Private(&vec->bstash, bnstash, breallocs));
 44:   PetscFunctionReturn(PETSC_SUCCESS);
 45: }

 47: /*@
 48:   VecSetLocalToGlobalMapping - Sets a local numbering to global numbering used
 49:   by the routine `VecSetValuesLocal()` to allow users to insert vector entries
 50:   using a local (per-processor) numbering.

 52:   Logically Collective

 54:   Input Parameters:
 55: + x       - vector
 56: - mapping - mapping created with `ISLocalToGlobalMappingCreate()` or `ISLocalToGlobalMappingCreateIS()`

 58:   Level: intermediate

 60:   Notes:
 61:   All vectors obtained with `VecDuplicate()` from this vector inherit the same mapping.

 63:   Vectors obtained with `DMCreateGlobaVector()` will often have this attribute attached to the vector so this call is not needed

 65: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecAssemblyEnd()`, `VecSetValues()`, `VecSetValuesLocal()`,
 66:            `VecGetLocalToGlobalMapping()`, `VecSetValuesBlockedLocal()`
 67: @*/
 68: PetscErrorCode VecSetLocalToGlobalMapping(Vec x, ISLocalToGlobalMapping mapping)
 69: {
 70:   PetscFunctionBegin;
 73:   if (x->ops->setlocaltoglobalmapping) PetscUseTypeMethod(x, setlocaltoglobalmapping, mapping);
 74:   else PetscCall(PetscLayoutSetISLocalToGlobalMapping(x->map, mapping));
 75:   PetscFunctionReturn(PETSC_SUCCESS);
 76: }

 78: /*@
 79:   VecGetLocalToGlobalMapping - Gets the local-to-global numbering set by `VecSetLocalToGlobalMapping()`

 81:   Not Collective

 83:   Input Parameter:
 84: . X - the vector

 86:   Output Parameter:
 87: . mapping - the mapping

 89:   Level: advanced

 91: .seealso: [](ch_vectors), `Vec`, `VecSetValuesLocal()`, `VecSetLocalToGlobalMapping()`
 92: @*/
 93: PetscErrorCode VecGetLocalToGlobalMapping(Vec X, ISLocalToGlobalMapping *mapping)
 94: {
 95:   PetscFunctionBegin;
 98:   PetscAssertPointer(mapping, 2);
 99:   if (X->ops->getlocaltoglobalmapping) PetscUseTypeMethod(X, getlocaltoglobalmapping, mapping);
100:   else *mapping = X->map->mapping;
101:   PetscFunctionReturn(PETSC_SUCCESS);
102: }

104: /*@
105:   VecAssemblyBegin - Begins assembling the vector; that is ensuring all the vector's entries are stored on the correct MPI process. This routine should
106:   be called after completing all calls to `VecSetValues()`.

108:   Collective

110:   Input Parameter:
111: . vec - the vector

113:   Level: beginner

115: .seealso: [](ch_vectors), `Vec`, `VecAssemblyEnd()`, `VecSetValues()`
116: @*/
117: PetscErrorCode VecAssemblyBegin(Vec vec)
118: {
119:   PetscFunctionBegin;
122:   PetscCall(VecStashViewFromOptions(vec, NULL, "-vec_view_stash"));
123:   PetscCall(PetscLogEventBegin(VEC_AssemblyBegin, vec, 0, 0, 0));
124:   PetscTryTypeMethod(vec, assemblybegin);
125:   PetscCall(PetscLogEventEnd(VEC_AssemblyBegin, vec, 0, 0, 0));
126:   PetscCall(PetscObjectStateIncrease((PetscObject)vec));
127:   PetscFunctionReturn(PETSC_SUCCESS);
128: }

130: /*@
131:   VecAssemblyEnd - Completes assembling the vector.  This routine should be called after `VecAssemblyBegin()`.

133:   Collective

135:   Input Parameter:
136: . vec - the vector

138:   Options Database Keys:
139: + -vec_view                 - Prints vector in `PETSC_VIEWER_DEFAULT` format
140: . -vec_view ::ascii_matlab  - Prints vector in `PETSC_VIEWER_ASCII_MATLAB` format to stdout
141: . -vec_view matlab:filename - Prints vector in MATLAB .mat file to filename (requires PETSc configured with --with-matlab)
142: . -vec_view draw            - Activates vector viewing using drawing tools
143: . -display <name>           - Sets display name (default is host)
144: . -draw_pause <sec>         - Sets number of seconds to pause after display
145: - -vec_view socket          - Activates vector viewing using a socket

147:   Level: beginner

149: .seealso: [](ch_vectors), `Vec`, `VecAssemblyBegin()`, `VecSetValues()`
150: @*/
151: PetscErrorCode VecAssemblyEnd(Vec vec)
152: {
153:   PetscFunctionBegin;
155:   PetscCall(PetscLogEventBegin(VEC_AssemblyEnd, vec, 0, 0, 0));
157:   PetscTryTypeMethod(vec, assemblyend);
158:   PetscCall(PetscLogEventEnd(VEC_AssemblyEnd, vec, 0, 0, 0));
159:   PetscCall(VecViewFromOptions(vec, NULL, "-vec_view"));
160:   PetscFunctionReturn(PETSC_SUCCESS);
161: }

163: /*@
164:   VecSetPreallocationCOO - set preallocation for a vector using a coordinate format of the entries with global indices

166:   Collective

168:   Input Parameters:
169: + x     - vector being preallocated
170: . ncoo  - number of entries
171: - coo_i - entry indices

173:   Level: beginner

175:   Notes:
176:   This and `VecSetValuesCOO()` provide an alternative API to using `VecSetValues()` to provide vector values.

178:   This API is particularly efficient for use on GPUs.

180:   Entries can be repeated, see `VecSetValuesCOO()`. Negative indices are not allowed unless vector option `VEC_IGNORE_NEGATIVE_INDICES` is set,
181:   in which case they, along with the corresponding entries in `VecSetValuesCOO()`, are ignored. If vector option `VEC_NO_OFF_PROC_ENTRIES` is set,
182:   remote entries are ignored, otherwise, they will be properly added or inserted to the vector.

184:   The array coo_i[] may be freed immediately after calling this function.

186: .seealso: [](ch_vectors), `Vec`, `VecSetValuesCOO()`, `VecSetPreallocationCOOLocal()`
187: @*/
188: PetscErrorCode VecSetPreallocationCOO(Vec x, PetscCount ncoo, const PetscInt coo_i[])
189: {
190:   PetscFunctionBegin;
193:   if (ncoo) PetscAssertPointer(coo_i, 3);
194:   PetscCall(PetscLogEventBegin(VEC_SetPreallocateCOO, x, 0, 0, 0));
195:   PetscCall(PetscLayoutSetUp(x->map));
196:   if (x->ops->setpreallocationcoo) {
197:     PetscUseTypeMethod(x, setpreallocationcoo, ncoo, coo_i);
198:   } else {
199:     IS is_coo_i;
200:     /* The default implementation only supports ncoo within limit of PetscInt */
201:     PetscCheck(ncoo <= PETSC_MAX_INT, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "ncoo %" PetscCount_FMT " overflowed PetscInt; configure --with-64-bit-indices or request support", ncoo);
202:     PetscCall(ISCreateGeneral(PETSC_COMM_SELF, ncoo, coo_i, PETSC_COPY_VALUES, &is_coo_i));
203:     PetscCall(PetscObjectCompose((PetscObject)x, "__PETSc_coo_i", (PetscObject)is_coo_i));
204:     PetscCall(ISDestroy(&is_coo_i));
205:   }
206:   PetscCall(PetscLogEventEnd(VEC_SetPreallocateCOO, x, 0, 0, 0));
207:   PetscFunctionReturn(PETSC_SUCCESS);
208: }

210: /*@
211:   VecSetPreallocationCOOLocal - set preallocation for vectors using a coordinate format of the entries with local indices

213:   Collective

215:   Input Parameters:
216: + x     - vector being preallocated
217: . ncoo  - number of entries
218: - coo_i - row indices (local numbering; may be modified)

220:   Level: beginner

222:   Notes:
223:   This and `VecSetValuesCOO()` provide an alternative API to using `VecSetValuesLocal()` to provide vector values.

225:   This API is particularly efficient for use on GPUs.

227:   The local indices are translated using the local to global mapping, thus `VecSetLocalToGlobalMapping()` must have been
228:   called prior to this function.

230:   The indices coo_i may be modified within this function. They might be translated to corresponding global
231:   indices, but the caller should not rely on them having any specific value after this function returns. The arrays
232:   can be freed or reused immediately after this function returns.

234:   Entries can be repeated. Negative indices and remote indices might be allowed. see `VecSetPreallocationCOO()`.

236: .seealso: [](ch_vectors), `Vec`, `VecSetPreallocationCOO()`, `VecSetValuesCOO()`
237: @*/
238: PetscErrorCode VecSetPreallocationCOOLocal(Vec x, PetscCount ncoo, PetscInt coo_i[])
239: {
240:   ISLocalToGlobalMapping ltog;

242:   PetscFunctionBegin;
245:   if (ncoo) PetscAssertPointer(coo_i, 3);
246:   PetscCheck(ncoo <= PETSC_MAX_INT, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "ncoo %" PetscCount_FMT " overflowed PetscInt; configure --with-64-bit-indices or request support", ncoo);
247:   PetscCall(PetscLayoutSetUp(x->map));
248:   PetscCall(VecGetLocalToGlobalMapping(x, &ltog));
249:   if (ltog) PetscCall(ISLocalToGlobalMappingApply(ltog, ncoo, coo_i, coo_i));
250:   PetscCall(VecSetPreallocationCOO(x, ncoo, coo_i));
251:   PetscFunctionReturn(PETSC_SUCCESS);
252: }

254: /*@
255:   VecSetValuesCOO - set values at once in a vector preallocated using `VecSetPreallocationCOO()`

257:   Collective

259:   Input Parameters:
260: + x     - vector being set
261: . coo_v - the value array
262: - imode - the insert mode

264:   Level: beginner

266:   Note:
267:   This and `VecSetPreallocationCOO() or ``VecSetPreallocationCOOLocal()` provide an alternative API to using `VecSetValues()` to provide vector values.

269:   This API is particularly efficient for use on GPUs.

271:   The values must follow the order of the indices prescribed with `VecSetPreallocationCOO()` or `VecSetPreallocationCOOLocal()`.
272:   When repeated entries are specified in the COO indices the `coo_v` values are first properly summed, regardless of the value of `imode`.
273:   The imode flag indicates if `coo_v` must be added to the current values of the vector (`ADD_VALUES`) or overwritten (`INSERT_VALUES`).
274:   `VecAssemblyBegin()` and `VecAssemblyEnd()` do not need to be called after this routine. It automatically handles the assembly process.

276: .seealso: [](ch_vectors), `Vec`, `VecSetPreallocationCOO()`, `VecSetPreallocationCOOLocal()`, `VecSetValues()`
277: @*/
278: PetscErrorCode VecSetValuesCOO(Vec x, const PetscScalar coo_v[], InsertMode imode)
279: {
280:   PetscFunctionBegin;
284:   PetscCall(PetscLogEventBegin(VEC_SetValuesCOO, x, 0, 0, 0));
285:   if (x->ops->setvaluescoo) {
286:     PetscUseTypeMethod(x, setvaluescoo, coo_v, imode);
287:     PetscCall(PetscObjectStateIncrease((PetscObject)x));
288:   } else {
289:     IS              is_coo_i;
290:     const PetscInt *coo_i;
291:     PetscInt        ncoo;
292:     PetscMemType    mtype;

294:     PetscCall(PetscGetMemType(coo_v, &mtype));
295:     PetscCheck(mtype == PETSC_MEMTYPE_HOST, PetscObjectComm((PetscObject)x), PETSC_ERR_ARG_WRONG, "The basic VecSetValuesCOO() only supports v[] on host");
296:     PetscCall(PetscObjectQuery((PetscObject)x, "__PETSc_coo_i", (PetscObject *)&is_coo_i));
297:     PetscCheck(is_coo_i, PetscObjectComm((PetscObject)x), PETSC_ERR_COR, "Missing coo_i IS");
298:     PetscCall(ISGetLocalSize(is_coo_i, &ncoo));
299:     PetscCall(ISGetIndices(is_coo_i, &coo_i));
300:     if (imode != ADD_VALUES) PetscCall(VecZeroEntries(x));
301:     PetscCall(VecSetValues(x, ncoo, coo_i, coo_v, ADD_VALUES));
302:     PetscCall(ISRestoreIndices(is_coo_i, &coo_i));
303:     PetscCall(VecAssemblyBegin(x));
304:     PetscCall(VecAssemblyEnd(x));
305:   }
306:   PetscCall(PetscLogEventEnd(VEC_SetValuesCOO, x, 0, 0, 0));
307:   PetscFunctionReturn(PETSC_SUCCESS);
308: }

310: static PetscErrorCode VecPointwiseApply_Private(Vec w, Vec x, Vec y, PetscDeviceContext dctx, PetscLogEvent event, const char async_name[], PetscErrorCode (*const pointwise_op)(Vec, Vec, Vec))
311: {
312:   PetscErrorCode (*async_fn)(Vec, Vec, Vec, PetscDeviceContext) = NULL;

314:   PetscFunctionBegin;
321:   PetscCheckSameTypeAndComm(x, 2, y, 3);
322:   PetscCheckSameTypeAndComm(y, 3, w, 1);
323:   VecCheckSameSize(w, 1, x, 2);
324:   VecCheckSameSize(w, 1, y, 3);
325:   VecCheckAssembled(x);
326:   VecCheckAssembled(y);
327:   PetscCall(VecSetErrorIfLocked(w, 1));

330:   if (dctx) PetscCall(PetscObjectQueryFunction((PetscObject)w, async_name, &async_fn));
331:   if (event) PetscCall(PetscLogEventBegin(event, x, y, w, 0));
332:   if (async_fn) {
333:     PetscCall((*async_fn)(w, x, y, dctx));
334:   } else {
335:     PetscCall((*pointwise_op)(w, x, y));
336:   }
337:   if (event) PetscCall(PetscLogEventEnd(event, x, y, w, 0));
338:   PetscCall(PetscObjectStateIncrease((PetscObject)w));
339:   PetscFunctionReturn(PETSC_SUCCESS);
340: }

342: PetscErrorCode VecPointwiseMaxAsync_Private(Vec w, Vec x, Vec y, PetscDeviceContext dctx)
343: {
344:   PetscFunctionBegin;
345:   // REVIEW ME: no log event?
346:   PetscCall(VecPointwiseApply_Private(w, x, y, dctx, 0, VecAsyncFnName(PointwiseMax), w->ops->pointwisemax));
347:   PetscFunctionReturn(PETSC_SUCCESS);
348: }

350: /*@
351:   VecPointwiseMax - Computes the component-wise maximum `w[i] = max(x[i], y[i])`.

353:   Logically Collective

355:   Input Parameters:
356: + x - the first input vector
357: - y - the second input vector

359:   Output Parameter:
360: . w - the result

362:   Level: advanced

364:   Notes:
365:   Any subset of the `x`, `y`, and `w` may be the same vector.

367:   For complex numbers compares only the real part

369: .seealso: [](ch_vectors), `Vec`, `VecPointwiseDivide()`, `VecPointwiseMult()`, `VecPointwiseMin()`, `VecPointwiseMaxAbs()`, `VecMaxPointwiseDivide()`
370: @*/
371: PetscErrorCode VecPointwiseMax(Vec w, Vec x, Vec y)
372: {
373:   PetscFunctionBegin;
374:   PetscCall(VecPointwiseMaxAsync_Private(w, x, y, NULL));
375:   PetscFunctionReturn(PETSC_SUCCESS);
376: }

378: PetscErrorCode VecPointwiseMinAsync_Private(Vec w, Vec x, Vec y, PetscDeviceContext dctx)
379: {
380:   PetscFunctionBegin;
381:   // REVIEW ME: no log event?
382:   PetscCall(VecPointwiseApply_Private(w, x, y, dctx, 0, VecAsyncFnName(PointwiseMin), w->ops->pointwisemin));
383:   PetscFunctionReturn(PETSC_SUCCESS);
384: }

386: /*@
387:   VecPointwiseMin - Computes the component-wise minimum `w[i] = min(x[i], y[i])`.

389:   Logically Collective

391:   Input Parameters:
392: + x - the first input vector
393: - y - the second input vector

395:   Output Parameter:
396: . w - the result

398:   Level: advanced

400:   Notes:
401:   Any subset of the `x`, `y`, and `w` may be the same vector.

403:   For complex numbers compares only the real part

405: .seealso: [](ch_vectors), `Vec`, `VecPointwiseDivide()`, `VecPointwiseMult()`, `VecPointwiseMaxAbs()`, `VecMaxPointwiseDivide()`
406: @*/
407: PetscErrorCode VecPointwiseMin(Vec w, Vec x, Vec y)
408: {
409:   PetscFunctionBegin;
410:   PetscCall(VecPointwiseMinAsync_Private(w, x, y, NULL));
411:   PetscFunctionReturn(PETSC_SUCCESS);
412: }

414: PetscErrorCode VecPointwiseMaxAbsAsync_Private(Vec w, Vec x, Vec y, PetscDeviceContext dctx)
415: {
416:   PetscFunctionBegin;
417:   // REVIEW ME: no log event?
418:   PetscCall(VecPointwiseApply_Private(w, x, y, dctx, 0, VecAsyncFnName(PointwiseMaxAbs), w->ops->pointwisemaxabs));
419:   PetscFunctionReturn(PETSC_SUCCESS);
420: }

422: /*@
423:   VecPointwiseMaxAbs - Computes the component-wise maximum of the absolute values `w[i] = max(abs(x[i]), abs(y[i]))`.

425:   Logically Collective

427:   Input Parameters:
428: + x - the first input vector
429: - y - the second input vector

431:   Output Parameter:
432: . w - the result

434:   Level: advanced

436:   Notes:
437:   Any subset of the `x`, `y`, and `w` may be the same vector.

439: .seealso: [](ch_vectors), `Vec`, `VecPointwiseDivide()`, `VecPointwiseMult()`, `VecPointwiseMin()`, `VecPointwiseMax()`, `VecMaxPointwiseDivide()`
440: @*/
441: PetscErrorCode VecPointwiseMaxAbs(Vec w, Vec x, Vec y)
442: {
443:   PetscFunctionBegin;
444:   PetscCall(VecPointwiseMaxAbsAsync_Private(w, x, y, NULL));
445:   PetscFunctionReturn(PETSC_SUCCESS);
446: }

448: PetscErrorCode VecPointwiseDivideAsync_Private(Vec w, Vec x, Vec y, PetscDeviceContext dctx)
449: {
450:   PetscFunctionBegin;
451:   // REVIEW ME: no log event?
452:   PetscCall(VecPointwiseApply_Private(w, x, y, dctx, 0, VecAsyncFnName(PointwiseDivide), w->ops->pointwisedivide));
453:   PetscFunctionReturn(PETSC_SUCCESS);
454: }

456: /*@
457:   VecPointwiseDivide - Computes the component-wise division `w[i] = x[i] / y[i]`.

459:   Logically Collective

461:   Input Parameters:
462: + x - the numerator vector
463: - y - the denominator vector

465:   Output Parameter:
466: . w - the result

468:   Level: advanced

470:   Note:
471:   Any subset of the `x`, `y`, and `w` may be the same vector.

473: .seealso: [](ch_vectors), `Vec`, `VecPointwiseMult()`, `VecPointwiseMax()`, `VecPointwiseMin()`, `VecPointwiseMaxAbs()`, `VecMaxPointwiseDivide()`
474: @*/
475: PetscErrorCode VecPointwiseDivide(Vec w, Vec x, Vec y)
476: {
477:   PetscFunctionBegin;
478:   PetscCall(VecPointwiseDivideAsync_Private(w, x, y, NULL));
479:   PetscFunctionReturn(PETSC_SUCCESS);
480: }

482: PetscErrorCode VecPointwiseMultAsync_Private(Vec w, Vec x, Vec y, PetscDeviceContext dctx)
483: {
484:   PetscFunctionBegin;
486:   PetscCall(VecPointwiseApply_Private(w, x, y, dctx, VEC_PointwiseMult, VecAsyncFnName(PointwiseMult), w->ops->pointwisemult));
487:   PetscFunctionReturn(PETSC_SUCCESS);
488: }

490: /*@
491:   VecPointwiseMult - Computes the component-wise multiplication `w[i] = x[i] * y[i]`.

493:   Logically Collective

495:   Input Parameters:
496: + x - the first vector
497: - y - the second vector

499:   Output Parameter:
500: . w - the result

502:   Level: advanced

504:   Note:
505:   Any subset of the `x`, `y`, and `w` may be the same vector.

507: .seealso: [](ch_vectors), `Vec`, `VecPointwiseDivide()`, `VecPointwiseMax()`, `VecPointwiseMin()`, `VecPointwiseMaxAbs()`, `VecMaxPointwiseDivide()`
508: @*/
509: PetscErrorCode VecPointwiseMult(Vec w, Vec x, Vec y)
510: {
511:   PetscFunctionBegin;
512:   PetscCall(VecPointwiseMultAsync_Private(w, x, y, NULL));
513:   PetscFunctionReturn(PETSC_SUCCESS);
514: }

516: /*@
517:   VecDuplicate - Creates a new vector of the same type as an existing vector.

519:   Collective

521:   Input Parameter:
522: . v - a vector to mimic

524:   Output Parameter:
525: . newv - location to put new vector

527:   Level: beginner

529:   Notes:
530:   `VecDuplicate()` DOES NOT COPY the vector entries, but rather allocates storage
531:   for the new vector.  Use `VecCopy()` to copy a vector.

533:   Use `VecDestroy()` to free the space. Use `VecDuplicateVecs()` to get several
534:   vectors.

536: .seealso: [](ch_vectors), `Vec`, `VecDestroy()`, `VecDuplicateVecs()`, `VecCreate()`, `VecCopy()`
537: @*/
538: PetscErrorCode VecDuplicate(Vec v, Vec *newv)
539: {
540:   PetscFunctionBegin;
542:   PetscAssertPointer(newv, 2);
544:   PetscUseTypeMethod(v, duplicate, newv);
545: #if PetscDefined(HAVE_DEVICE)
546:   if (v->boundtocpu && v->bindingpropagates) {
547:     PetscCall(VecSetBindingPropagates(*newv, PETSC_TRUE));
548:     PetscCall(VecBindToCPU(*newv, PETSC_TRUE));
549:   }
550: #endif
551:   PetscCall(PetscObjectStateIncrease((PetscObject)*newv));
552:   PetscFunctionReturn(PETSC_SUCCESS);
553: }

555: /*@C
556:   VecDestroy - Destroys a vector.

558:   Collective

560:   Input Parameter:
561: . v - the vector

563:   Level: beginner

565: .seealso: [](ch_vectors), `Vec`, `VecCreate()`, `VecDuplicate()`, `VecDestroyVecs()`
566: @*/
567: PetscErrorCode VecDestroy(Vec *v)
568: {
569:   PetscFunctionBegin;
570:   PetscAssertPointer(v, 1);
571:   if (!*v) PetscFunctionReturn(PETSC_SUCCESS);
573:   if (--((PetscObject)*v)->refct > 0) {
574:     *v = NULL;
575:     PetscFunctionReturn(PETSC_SUCCESS);
576:   }

578:   PetscCall(PetscObjectSAWsViewOff((PetscObject)*v));
579:   /* destroy the internal part */
580:   PetscTryTypeMethod(*v, destroy);
581:   PetscCall(PetscFree((*v)->defaultrandtype));
582:   /* destroy the external/common part */
583:   PetscCall(PetscLayoutDestroy(&(*v)->map));
584:   PetscCall(PetscHeaderDestroy(v));
585:   PetscFunctionReturn(PETSC_SUCCESS);
586: }

588: /*@C
589:   VecDuplicateVecs - Creates several vectors of the same type as an existing vector.

591:   Collective

593:   Input Parameters:
594: + m - the number of vectors to obtain
595: - v - a vector to mimic

597:   Output Parameter:
598: . V - location to put pointer to array of vectors

600:   Level: intermediate

602:   Note:
603:   Use `VecDestroyVecs()` to free the space. Use `VecDuplicate()` to form a single
604:   vector.

606:   Fortran Notes:
607:   The Fortran interface is slightly different from that given below, it
608:   requires one to pass in `V` a `Vec` array of size at least `m`.
609:   See the [](ch_fortran) for details.

611: .seealso: [](ch_vectors), `Vec`, [](ch_fortran), `VecDestroyVecs()`, `VecDuplicate()`, `VecCreate()`, `VecDuplicateVecsF90()`
612: @*/
613: PetscErrorCode VecDuplicateVecs(Vec v, PetscInt m, Vec *V[])
614: {
615:   PetscFunctionBegin;
617:   PetscAssertPointer(V, 3);
619:   PetscUseTypeMethod(v, duplicatevecs, m, V);
620: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
621:   if (v->boundtocpu && v->bindingpropagates) {
622:     PetscInt i;

624:     for (i = 0; i < m; i++) {
625:       /* Since ops->duplicatevecs might itself propagate the value of boundtocpu,
626:        * avoid unnecessary overhead by only calling VecBindToCPU() if the vector isn't already bound. */
627:       if (!(*V)[i]->boundtocpu) {
628:         PetscCall(VecSetBindingPropagates((*V)[i], PETSC_TRUE));
629:         PetscCall(VecBindToCPU((*V)[i], PETSC_TRUE));
630:       }
631:     }
632:   }
633: #endif
634:   PetscFunctionReturn(PETSC_SUCCESS);
635: }

637: /*@C
638:   VecDestroyVecs - Frees a block of vectors obtained with `VecDuplicateVecs()`.

640:   Collective

642:   Input Parameters:
643: + m  - the number of vectors previously obtained, if zero no vectors are destroyed
644: - vv - pointer to pointer to array of vector pointers, if `NULL` no vectors are destroyed

646:   Level: intermediate

648:   Fortran Notes:
649:   The Fortran interface is slightly different from that given below.
650:   See the [](ch_fortran) for details.

652: .seealso: [](ch_vectors), `Vec`, [](ch_fortran), `VecDuplicateVecs()`, `VecDestroyVecsf90()`
653: @*/
654: PetscErrorCode VecDestroyVecs(PetscInt m, Vec *vv[])
655: {
656:   PetscFunctionBegin;
657:   PetscAssertPointer(vv, 2);
658:   PetscCheck(m >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Trying to destroy negative number of vectors %" PetscInt_FMT, m);
659:   if (!m || !*vv) {
660:     *vv = NULL;
661:     PetscFunctionReturn(PETSC_SUCCESS);
662:   }
665:   PetscCall((*(**vv)->ops->destroyvecs)(m, *vv));
666:   *vv = NULL;
667:   PetscFunctionReturn(PETSC_SUCCESS);
668: }

670: /*@C
671:   VecViewFromOptions - View a vector based on values in the options database

673:   Collective

675:   Input Parameters:
676: + A    - the vector
677: . obj  - Optional object that provides the options prefix for this viewing
678: - name - command line option

680:   Level: intermediate

682:   Note:
683:   See `PetscObjectViewFromOptions()` to see the `PetscViewer` and PetscViewerFormat` available

685: .seealso: [](ch_vectors), `Vec`, `VecView`, `PetscObjectViewFromOptions()`, `VecCreate()`
686: @*/
687: PetscErrorCode VecViewFromOptions(Vec A, PetscObject obj, const char name[])
688: {
689:   PetscFunctionBegin;
691:   PetscCall(PetscObjectViewFromOptions((PetscObject)A, obj, name));
692:   PetscFunctionReturn(PETSC_SUCCESS);
693: }

695: /*@C
696:   VecView - Views a vector object.

698:   Collective

700:   Input Parameters:
701: + vec    - the vector
702: - viewer - an optional `PetscViewer` visualization context

704:   Level: beginner

706:   Notes:
707:   The available visualization contexts include
708: +     `PETSC_VIEWER_STDOUT_SELF` - for sequential vectors
709: .     `PETSC_VIEWER_STDOUT_WORLD` - for parallel vectors created on `PETSC_COMM_WORLD`
710: -     `PETSC_VIEWER_STDOUT`_(comm) - for parallel vectors created on MPI communicator comm

712:   You can change the format the vector is printed using the
713:   option `PetscViewerPushFormat()`.

715:   The user can open alternative viewers with
716: +    `PetscViewerASCIIOpen()` - Outputs vector to a specified file
717: .    `PetscViewerBinaryOpen()` - Outputs vector in binary to a
718:   specified file; corresponding input uses `VecLoad()`
719: .    `PetscViewerDrawOpen()` - Outputs vector to an X window display
720: .    `PetscViewerSocketOpen()` - Outputs vector to Socket viewer
721: -    `PetscViewerHDF5Open()` - Outputs vector to HDF5 file viewer

723:   The user can call `PetscViewerPushFormat()` to specify the output
724:   format of ASCII printed objects (when using `PETSC_VIEWER_STDOUT_SELF`,
725:   `PETSC_VIEWER_STDOUT_WORLD` and `PetscViewerASCIIOpen()`).  Available formats include
726: +    `PETSC_VIEWER_DEFAULT` - default, prints vector contents
727: .    `PETSC_VIEWER_ASCII_MATLAB` - prints vector contents in MATLAB format
728: .    `PETSC_VIEWER_ASCII_INDEX` - prints vector contents, including indices of vector elements
729: -    `PETSC_VIEWER_ASCII_COMMON` - prints vector contents, using a
730:   format common among all vector types

732:   You can pass any number of vector objects, or other PETSc objects to the same viewer.

734:   In the debugger you can do call `VecView`(v,0) to display the vector. (The same holds for any PETSc object viewer).

736:   Notes for binary viewer:
737:   If you pass multiple vectors to a binary viewer you can read them back in the same order
738:   with `VecLoad()`.

740:   If the blocksize of the vector is greater than one then you must provide a unique prefix to
741:   the vector with `PetscObjectSetOptionsPrefix`((`PetscObject`)vec,"uniqueprefix"); BEFORE calling `VecView()` on the
742:   vector to be stored and then set that same unique prefix on the vector that you pass to `VecLoad()`. The blocksize
743:   information is stored in an ASCII file with the same name as the binary file plus a ".info" appended to the
744:   filename. If you copy the binary file, make sure you copy the associated .info file with it.

746:   See the manual page for `VecLoad()` on the exact format the binary viewer stores
747:   the values in the file.

749:   Notes for HDF5 Viewer:
750:   The name of the `Vec` (given with `PetscObjectSetName()` is the name that is used
751:   for the object in the HDF5 file. If you wish to store the same Vec into multiple
752:   datasets in the same file (typically with different values), you must change its
753:   name each time before calling the `VecView()`. To load the same vector,
754:   the name of the Vec object passed to `VecLoad()` must be the same.

756:   If the block size of the vector is greater than 1 then it is used as the first dimension in the HDF5 array.
757:   If the function `PetscViewerHDF5SetBaseDimension2()`is called then even if the block size is one it will
758:   be used as the first dimension in the HDF5 array (that is the HDF5 array will always be two dimensional)
759:   See also `PetscViewerHDF5SetTimestep()` which adds an additional complication to reading and writing `Vec`
760:   with the HDF5 viewer.

762: .seealso: [](ch_vectors), `Vec`, `VecViewFromOptions()`, `PetscViewerASCIIOpen()`, `PetscViewerDrawOpen()`, `PetscDrawLGCreate()`,
763:           `PetscViewerSocketOpen()`, `PetscViewerBinaryOpen()`, `VecLoad()`, `PetscViewerCreate()`,
764:           `PetscRealView()`, `PetscScalarView()`, `PetscIntView()`, `PetscViewerHDF5SetTimestep()`
765: @*/
766: PetscErrorCode VecView(Vec vec, PetscViewer viewer)
767: {
768:   PetscBool         iascii;
769:   PetscViewerFormat format;
770:   PetscMPIInt       size;

772:   PetscFunctionBegin;
775:   VecCheckAssembled(vec);
776:   if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)vec), &viewer));
778:   PetscCall(PetscViewerGetFormat(viewer, &format));
779:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)vec), &size));
780:   if (size == 1 && format == PETSC_VIEWER_LOAD_BALANCE) PetscFunctionReturn(PETSC_SUCCESS);

782:   PetscCheck(!vec->stash.n && !vec->bstash.n, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call VecAssemblyBegin/End() before viewing this vector");

784:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
785:   if (iascii) {
786:     PetscInt rows, bs;

788:     PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)vec, viewer));
789:     if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
790:       PetscCall(PetscViewerASCIIPushTab(viewer));
791:       PetscCall(VecGetSize(vec, &rows));
792:       PetscCall(VecGetBlockSize(vec, &bs));
793:       if (bs != 1) {
794:         PetscCall(PetscViewerASCIIPrintf(viewer, "length=%" PetscInt_FMT ", bs=%" PetscInt_FMT "\n", rows, bs));
795:       } else {
796:         PetscCall(PetscViewerASCIIPrintf(viewer, "length=%" PetscInt_FMT "\n", rows));
797:       }
798:       PetscCall(PetscViewerASCIIPopTab(viewer));
799:     }
800:   }
801:   PetscCall(VecLockReadPush(vec));
802:   PetscCall(PetscLogEventBegin(VEC_View, vec, viewer, 0, 0));
803:   if ((format == PETSC_VIEWER_NATIVE || format == PETSC_VIEWER_LOAD_BALANCE) && vec->ops->viewnative) {
804:     PetscUseTypeMethod(vec, viewnative, viewer);
805:   } else {
806:     PetscUseTypeMethod(vec, view, viewer);
807:   }
808:   PetscCall(VecLockReadPop(vec));
809:   PetscCall(PetscLogEventEnd(VEC_View, vec, viewer, 0, 0));
810:   PetscFunctionReturn(PETSC_SUCCESS);
811: }

813: #if defined(PETSC_USE_DEBUG)
814: #include <../src/sys/totalview/tv_data_display.h>
815: PETSC_UNUSED static int TV_display_type(const struct _p_Vec *v)
816: {
817:   const PetscScalar *values;
818:   char               type[32];

820:   TV_add_row("Local rows", "int", &v->map->n);
821:   TV_add_row("Global rows", "int", &v->map->N);
822:   TV_add_row("Typename", TV_ascii_string_type, ((PetscObject)v)->type_name);
823:   PetscCall(VecGetArrayRead((Vec)v, &values));
824:   PetscCall(PetscSNPrintf(type, 32, "double[%" PetscInt_FMT "]", v->map->n));
825:   TV_add_row("values", type, values);
826:   PetscCall(VecRestoreArrayRead((Vec)v, &values));
827:   return TV_format_OK;
828: }
829: #endif

831: /*@C
832:   VecViewNative - Views a vector object with the original type specific viewer

834:   Collective

836:   Input Parameters:
837: + vec    - the vector
838: - viewer - an optional `PetscViewer` visualization context

840:   Level: developer

842:   Note:
843:   This can be used with, for example, vectors obtained with `DMCreateGlobalVector()` for a `DMDA` to display the vector
844:   in the PETSc storage format (each MPI process values follow the previous MPI processes) instead of the "natural" grid
845:   ordering.

847: .seealso: [](ch_vectors), `Vec`, `PetscViewerASCIIOpen()`, `PetscViewerDrawOpen()`, `PetscDrawLGCreate()`, `VecView()`
848:           `PetscViewerSocketOpen()`, `PetscViewerBinaryOpen()`, `VecLoad()`, `PetscViewerCreate()`,
849:           `PetscRealView()`, `PetscScalarView()`, `PetscIntView()`, `PetscViewerHDF5SetTimestep()`
850: @*/
851: PetscErrorCode VecViewNative(Vec vec, PetscViewer viewer)
852: {
853:   PetscFunctionBegin;
856:   if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)vec), &viewer));
858:   PetscUseTypeMethod(vec, viewnative, viewer);
859:   PetscFunctionReturn(PETSC_SUCCESS);
860: }

862: /*@
863:   VecGetSize - Returns the global number of elements of the vector.

865:   Not Collective

867:   Input Parameter:
868: . x - the vector

870:   Output Parameter:
871: . size - the global length of the vector

873:   Level: beginner

875: .seealso: [](ch_vectors), `Vec`, `VecGetLocalSize()`
876: @*/
877: PetscErrorCode VecGetSize(Vec x, PetscInt *size)
878: {
879:   PetscFunctionBegin;
881:   PetscAssertPointer(size, 2);
883:   PetscUseTypeMethod(x, getsize, size);
884:   PetscFunctionReturn(PETSC_SUCCESS);
885: }

887: /*@
888:   VecGetLocalSize - Returns the number of elements of the vector stored
889:   in local memory (that is on this MPI process)

891:   Not Collective

893:   Input Parameter:
894: . x - the vector

896:   Output Parameter:
897: . size - the length of the local piece of the vector

899:   Level: beginner

901: .seealso: [](ch_vectors), `Vec`, `VecGetSize()`
902: @*/
903: PetscErrorCode VecGetLocalSize(Vec x, PetscInt *size)
904: {
905:   PetscFunctionBegin;
907:   PetscAssertPointer(size, 2);
909:   PetscUseTypeMethod(x, getlocalsize, size);
910:   PetscFunctionReturn(PETSC_SUCCESS);
911: }

913: /*@C
914:   VecGetOwnershipRange - Returns the range of indices owned by
915:   this process. The vector is laid out with the
916:   first n1 elements on the first processor, next n2 elements on the
917:   second, etc.  For certain parallel layouts this range may not be
918:   well defined.

920:   Not Collective

922:   Input Parameter:
923: . x - the vector

925:   Output Parameters:
926: + low  - the first local element, pass in `NULL` if not interested
927: - high - one more than the last local element, pass in `NULL` if not interested

929:   Level: beginner

931:   Note:
932:   The high argument is one more than the last element stored locally.

934:   Fortran Notes:
935:   `PETSC_NULL_INTEGER` should be used instead of NULL

937: .seealso: [](ch_vectors), `Vec`, `MatGetOwnershipRange()`, `MatGetOwnershipRanges()`, `VecGetOwnershipRanges()`
938: @*/
939: PetscErrorCode VecGetOwnershipRange(Vec x, PetscInt *low, PetscInt *high)
940: {
941:   PetscFunctionBegin;
944:   if (low) PetscAssertPointer(low, 2);
945:   if (high) PetscAssertPointer(high, 3);
946:   if (low) *low = x->map->rstart;
947:   if (high) *high = x->map->rend;
948:   PetscFunctionReturn(PETSC_SUCCESS);
949: }

951: /*@C
952:   VecGetOwnershipRanges - Returns the range of indices owned by EACH processor,
953:   The vector is laid out with the
954:   first n1 elements on the first processor, next n2 elements on the
955:   second, etc.  For certain parallel layouts this range may not be
956:   well defined.

958:   Not Collective

960:   Input Parameter:
961: . x - the vector

963:   Output Parameter:
964: . ranges - array of length size+1 with the start and end+1 for each process

966:   Level: beginner

968:   Notes:
969:   The high argument is one more than the last element stored locally.

971:   If the ranges are used after all vectors that share the ranges has been destroyed then the program will crash accessing ranges[].

973:   Fortran Notes:
974:   You must PASS in an array of length size+1

976: .seealso: [](ch_vectors), `Vec`, `MatGetOwnershipRange()`, `MatGetOwnershipRanges()`, `VecGetOwnershipRange()`
977: @*/
978: PetscErrorCode VecGetOwnershipRanges(Vec x, const PetscInt *ranges[])
979: {
980:   PetscFunctionBegin;
983:   PetscCall(PetscLayoutGetRanges(x->map, ranges));
984:   PetscFunctionReturn(PETSC_SUCCESS);
985: }

987: // PetscClangLinter pragma disable: -fdoc-section-header-unknown
988: /*@
989:   VecSetOption - Sets an option for controlling a vector's behavior.

991:   Collective

993:   Input Parameters:
994: + x    - the vector
995: . op   - the option
996: - flag - turn the option on or off

998:   Supported Options:
999: + `VEC_IGNORE_OFF_PROC_ENTRIES` - which causes `VecSetValues()` to ignore
1000:           entries destined to be stored on a separate processor. This can be used
1001:           to eliminate the global reduction in the `VecAssemblyBegin()` if you know
1002:           that you have only used `VecSetValues()` to set local elements
1003: . `VEC_IGNORE_NEGATIVE_INDICES` - which means you can pass negative indices
1004:           in ix in calls to `VecSetValues()` or `VecGetValues()`. These rows are simply
1005:           ignored.
1006: - `VEC_SUBSET_OFF_PROC_ENTRIES` - which causes `VecAssemblyBegin()` to assume that the off-process
1007:           entries will always be a subset (possibly equal) of the off-process entries set on the
1008:           first assembly which had a true `VEC_SUBSET_OFF_PROC_ENTRIES` and the vector has not
1009:           changed this flag afterwards. If this assembly is not such first assembly, then this
1010:           assembly can reuse the communication pattern setup in that first assembly, thus avoiding
1011:           a global reduction. Subsequent assemblies setting off-process values should use the same
1012:           InsertMode as the first assembly.

1014:   Level: intermediate

1016:   Developer Notes:
1017:   The `InsertMode` restriction could be removed by packing the stash messages out of place.

1019: .seealso: [](ch_vectors), `Vec`, `VecSetValues()`
1020: @*/
1021: PetscErrorCode VecSetOption(Vec x, VecOption op, PetscBool flag)
1022: {
1023:   PetscFunctionBegin;
1026:   PetscTryTypeMethod(x, setoption, op, flag);
1027:   PetscFunctionReturn(PETSC_SUCCESS);
1028: }

1030: /* Default routines for obtaining and releasing; */
1031: /* may be used by any implementation */
1032: PetscErrorCode VecDuplicateVecs_Default(Vec w, PetscInt m, Vec *V[])
1033: {
1034:   PetscFunctionBegin;
1035:   PetscCheck(m > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "m must be > 0: m = %" PetscInt_FMT, m);
1036:   PetscCall(PetscMalloc1(m, V));
1037:   for (PetscInt i = 0; i < m; i++) PetscCall(VecDuplicate(w, *V + i));
1038:   PetscFunctionReturn(PETSC_SUCCESS);
1039: }

1041: PetscErrorCode VecDestroyVecs_Default(PetscInt m, Vec v[])
1042: {
1043:   PetscInt i;

1045:   PetscFunctionBegin;
1046:   PetscAssertPointer(v, 2);
1047:   for (i = 0; i < m; i++) PetscCall(VecDestroy(&v[i]));
1048:   PetscCall(PetscFree(v));
1049:   PetscFunctionReturn(PETSC_SUCCESS);
1050: }

1052: /*@
1053:   VecResetArray - Resets a vector to use its default memory. Call this
1054:   after the use of `VecPlaceArray()`.

1056:   Not Collective

1058:   Input Parameter:
1059: . vec - the vector

1061:   Level: developer

1063: .seealso: [](ch_vectors), `Vec`, `VecGetArray()`, `VecRestoreArray()`, `VecReplaceArray()`, `VecPlaceArray()`
1064: @*/
1065: PetscErrorCode VecResetArray(Vec vec)
1066: {
1067:   PetscFunctionBegin;
1070:   PetscUseTypeMethod(vec, resetarray);
1071:   PetscCall(PetscObjectStateIncrease((PetscObject)vec));
1072:   PetscFunctionReturn(PETSC_SUCCESS);
1073: }

1075: /*@C
1076:   VecLoad - Loads a vector that has been stored in binary or HDF5 format
1077:   with `VecView()`.

1079:   Collective

1081:   Input Parameters:
1082: + vec    - the newly loaded vector, this needs to have been created with `VecCreate()` or
1083:            some related function before the call to `VecLoad()`.
1084: - viewer - binary file viewer, obtained from `PetscViewerBinaryOpen()` or
1085:            HDF5 file viewer, obtained from `PetscViewerHDF5Open()`

1087:   Level: intermediate

1089:   Notes:
1090:   Defaults to the standard `VECSEQ` or `VECMPI`, if you want some other type of `Vec` call `VecSetFromOptions()`
1091:   before calling this.

1093:   The input file must contain the full global vector, as
1094:   written by the routine `VecView()`.

1096:   If the type or size of `vec` is not set before a call to `VecLoad()`, PETSc
1097:   sets the type and the local and global sizes based on the vector it is reading in. If type and/or
1098:   sizes are already set, then the same are used.

1100:   If using the binary viewer and the blocksize of the vector is greater than one then you must provide a unique prefix to
1101:   the vector with `PetscObjectSetOptionsPrefix`((`PetscObject`)vec,"uniqueprefix"); BEFORE calling `VecView()` on the
1102:   vector to be stored and then set that same unique prefix on the vector that you pass to VecLoad(). The blocksize
1103:   information is stored in an ASCII file with the same name as the binary file plus a ".info" appended to the
1104:   filename. If you copy the binary file, make sure you copy the associated .info file with it.

1106:   If using HDF5, you must assign the `Vec` the same name as was used in the Vec
1107:   that was stored in the file using `PetscObjectSetName(). Otherwise you will
1108:   get the error message: "Cannot H5DOpen2() with `Vec` name NAMEOFOBJECT".

1110:   If the HDF5 file contains a two dimensional array the first dimension is treated as the block size
1111:   in loading the vector. Hence, for example, using MATLAB notation h5create('vector.dat','/Test_Vec',[27 1]);
1112:   will load a vector of size 27 and block size 27 thus resulting in all 27 entries being on the first process of
1113:   vectors communicator and the rest of the processes having zero entries

1115:   Notes for advanced users when using the binary viewer:
1116:   Most users should not need to know the details of the binary storage
1117:   format, since `VecLoad()` and `VecView()` completely hide these details.
1118:   But for anyone who's interested, the standard binary vector storage
1119:   format is
1120: .vb
1121:      PetscInt    VEC_FILE_CLASSID
1122:      PetscInt    number of rows
1123:      PetscScalar *values of all entries
1124: .ve

1126:   In addition, PETSc automatically uses byte swapping to work on all machines; the files
1127:   are written ALWAYS using big-endian ordering. On small-endian machines the numbers
1128:   are converted to the small-endian format when they are read in from the file.
1129:   See PetscBinaryRead() and PetscBinaryWrite() to see how this may be done.

1131: .seealso: [](ch_vectors), `Vec`, `PetscViewerBinaryOpen()`, `VecView()`, `MatLoad()`
1132: @*/
1133: PetscErrorCode VecLoad(Vec vec, PetscViewer viewer)
1134: {
1135:   PetscBool         isbinary, ishdf5, isadios, isexodusii;
1136:   PetscViewerFormat format;

1138:   PetscFunctionBegin;
1141:   PetscCheckSameComm(vec, 1, viewer, 2);
1142:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
1143:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
1144:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERADIOS, &isadios));
1145:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWEREXODUSII, &isexodusii));
1146:   PetscCheck(isbinary || ishdf5 || isadios || isexodusii, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid viewer; open viewer with PetscViewerBinaryOpen()");

1148:   PetscCall(VecSetErrorIfLocked(vec, 1));
1149:   if (!((PetscObject)vec)->type_name && !vec->ops->create) PetscCall(VecSetType(vec, VECSTANDARD));
1150:   PetscCall(PetscLogEventBegin(VEC_Load, viewer, 0, 0, 0));
1151:   PetscCall(PetscViewerGetFormat(viewer, &format));
1152:   if (format == PETSC_VIEWER_NATIVE && vec->ops->loadnative) {
1153:     PetscUseTypeMethod(vec, loadnative, viewer);
1154:   } else {
1155:     PetscUseTypeMethod(vec, load, viewer);
1156:   }
1157:   PetscCall(PetscLogEventEnd(VEC_Load, viewer, 0, 0, 0));
1158:   PetscFunctionReturn(PETSC_SUCCESS);
1159: }

1161: /*@
1162:   VecReciprocal - Replaces each component of a vector by its reciprocal.

1164:   Logically Collective

1166:   Input Parameter:
1167: . vec - the vector

1169:   Output Parameter:
1170: . vec - the vector reciprocal

1172:   Level: intermediate

1174:   Note:
1175:   Vector entries with value 0.0 are not changed

1177: .seealso: [](ch_vectors), `Vec`, `VecLog()`, `VecExp()`, `VecSqrtAbs()`
1178: @*/
1179: PetscErrorCode VecReciprocal(Vec vec)
1180: {
1181:   PetscFunctionBegin;
1182:   PetscCall(VecReciprocalAsync_Private(vec, NULL));
1183:   PetscFunctionReturn(PETSC_SUCCESS);
1184: }

1186: /*@C
1187:   VecSetOperation - Allows the user to override a particular vector operation.

1189:   Logically Collective; No Fortran Support

1191:   Input Parameters:
1192: + vec - The vector to modify
1193: . op  - The name of the operation
1194: - f   - The function that provides the operation.

1196:   Notes:
1197:   `f` may be `NULL` to remove the operation from `vec`. Depending on the operation this may be
1198:   allowed, however some always expect a valid function. In these cases an error will be raised
1199:   when calling the interface routine in question.

1201:   See `VecOperation` for an up-to-date list of override-able operations. The operations listed
1202:   there have the form `VECOP_<OPERATION>`, where `<OPERATION>` is the suffix (in all capital
1203:   letters) of the public interface routine (e.g., `VecView()` -> `VECOP_VIEW`).

1205:   Overriding a particular `Vec`'s operation has no affect on any other `Vec`s past, present,
1206:   or future. The user should also note that overriding a method is "destructive"; the previous
1207:   method is not retained in any way.

1209:   Level: advanced

1211:   Example Usage:
1212: .vb
1213:   // some new VecView() implementation, must have the same signature as the function it seeks
1214:   // to replace
1215:   PetscErrorCode UserVecView(Vec x, PetscViewer viewer)
1216:   {
1217:     PetscFunctionBeginUser;
1218:     // ...
1219:     PetscFunctionReturn(PETSC_SUCCESS);
1220:   }

1222:   // Create a VECMPI which has a pre-defined VecView() implementation
1223:   VecCreateMPI(comm, n, N, &x);
1224:   // Calls the VECMPI implementation for VecView()
1225:   VecView(x, viewer);

1227:   VecSetOperation(x, VECOP_VIEW, (void (*)(void))UserVecView);
1228:   // Now calls UserVecView()
1229:   VecView(x, viewer);
1230: .ve

1232: .seealso: [](ch_vectors), `Vec`, `VecCreate()`, `MatShellSetOperation()`
1233: @*/
1234: PetscErrorCode VecSetOperation(Vec vec, VecOperation op, void (*f)(void))
1235: {
1236:   PetscFunctionBegin;
1238:   if (op == VECOP_VIEW && !vec->ops->viewnative) {
1239:     vec->ops->viewnative = vec->ops->view;
1240:   } else if (op == VECOP_LOAD && !vec->ops->loadnative) {
1241:     vec->ops->loadnative = vec->ops->load;
1242:   }
1243:   ((void (**)(void))vec->ops)[(int)op] = f;
1244:   PetscFunctionReturn(PETSC_SUCCESS);
1245: }

1247: /*@
1248:   VecStashSetInitialSize - sets the sizes of the vec-stash, that is
1249:   used during the assembly process to store values that belong to
1250:   other processors.

1252:   Not Collective, different processes can have different size stashes

1254:   Input Parameters:
1255: + vec   - the vector
1256: . size  - the initial size of the stash.
1257: - bsize - the initial size of the block-stash(if used).

1259:   Options Database Keys:
1260: + -vecstash_initial_size <size> or <size0,size1,...sizep-1>           - set initial size
1261: - -vecstash_block_initial_size <bsize> or <bsize0,bsize1,...bsizep-1> - set initial block size

1263:   Level: intermediate

1265:   Notes:
1266:   The block-stash is used for values set with `VecSetValuesBlocked()` while
1267:   the stash is used for values set with `VecSetValues()`

1269:   Run with the option -info and look for output of the form
1270:   VecAssemblyBegin_MPIXXX:Stash has MM entries, uses nn mallocs.
1271:   to determine the appropriate value, MM, to use for size and
1272:   VecAssemblyBegin_MPIXXX:Block-Stash has BMM entries, uses nn mallocs.
1273:   to determine the value, BMM to use for bsize

1275:   PETSc attempts to smartly manage the stash size so there is little likelihood setting a
1276:   a specific value here will affect performance

1278: .seealso: [](ch_vectors), `Vec`, `VecSetBlockSize()`, `VecSetValues()`, `VecSetValuesBlocked()`, `VecStashView()`
1279: @*/
1280: PetscErrorCode VecStashSetInitialSize(Vec vec, PetscInt size, PetscInt bsize)
1281: {
1282:   PetscFunctionBegin;
1284:   PetscCall(VecStashSetInitialSize_Private(&vec->stash, size));
1285:   PetscCall(VecStashSetInitialSize_Private(&vec->bstash, bsize));
1286:   PetscFunctionReturn(PETSC_SUCCESS);
1287: }

1289: /*@
1290:   VecSetRandom - Sets all components of a vector to random numbers.

1292:   Logically Collective

1294:   Input Parameters:
1295: + x    - the vector
1296: - rctx - the random number context, formed by `PetscRandomCreate()`, or use `NULL` and it will create one internally.

1298:   Output Parameter:
1299: . x - the vector

1301:   Example of Usage:
1302: .vb
1303:      PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
1304:      VecSetRandom(x,rctx);
1305:      PetscRandomDestroy(&rctx);
1306: .ve

1308:   Level: intermediate

1310: .seealso: [](ch_vectors), `Vec`, `VecSet()`, `VecSetValues()`, `PetscRandomCreate()`, `PetscRandomDestroy()`
1311: @*/
1312: PetscErrorCode VecSetRandom(Vec x, PetscRandom rctx)
1313: {
1314:   PetscRandom randObj = NULL;

1316:   PetscFunctionBegin;
1320:   VecCheckAssembled(x);
1321:   PetscCall(VecSetErrorIfLocked(x, 1));

1323:   if (!rctx) {
1324:     PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)x), &randObj));
1325:     PetscCall(PetscRandomSetType(randObj, x->defaultrandtype));
1326:     PetscCall(PetscRandomSetFromOptions(randObj));
1327:     rctx = randObj;
1328:   }

1330:   PetscCall(PetscLogEventBegin(VEC_SetRandom, x, rctx, 0, 0));
1331:   PetscUseTypeMethod(x, setrandom, rctx);
1332:   PetscCall(PetscLogEventEnd(VEC_SetRandom, x, rctx, 0, 0));

1334:   PetscCall(PetscRandomDestroy(&randObj));
1335:   PetscCall(PetscObjectStateIncrease((PetscObject)x));
1336:   PetscFunctionReturn(PETSC_SUCCESS);
1337: }

1339: /*@
1340:   VecZeroEntries - puts a `0.0` in each element of a vector

1342:   Logically Collective

1344:   Input Parameter:
1345: . vec - The vector

1347:   Level: beginner

1349:   Note:
1350:   If the norm of the vector is known to be zero then this skips the unneeded zeroing process

1352: .seealso: [](ch_vectors), `Vec`, `VecCreate()`, `VecSetOptionsPrefix()`, `VecSet()`, `VecSetValues()`
1353: @*/
1354: PetscErrorCode VecZeroEntries(Vec vec)
1355: {
1356:   PetscFunctionBegin;
1357:   PetscCall(VecSet(vec, 0));
1358:   PetscFunctionReturn(PETSC_SUCCESS);
1359: }

1361: /*
1362:   VecSetTypeFromOptions_Private - Sets the type of vector from user options. Defaults to a PETSc sequential vector on one
1363:   processor and a PETSc MPI vector on more than one processor.

1365:   Collective

1367:   Input Parameter:
1368: . vec - The vector

1370:   Level: intermediate

1372: .seealso: [](ch_vectors), `Vec`, `VecSetFromOptions()`, `VecSetType()`
1373: */
1374: static PetscErrorCode VecSetTypeFromOptions_Private(Vec vec, PetscOptionItems *PetscOptionsObject)
1375: {
1376:   PetscBool   opt;
1377:   VecType     defaultType;
1378:   char        typeName[256];
1379:   PetscMPIInt size;

1381:   PetscFunctionBegin;
1382:   if (((PetscObject)vec)->type_name) defaultType = ((PetscObject)vec)->type_name;
1383:   else {
1384:     PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)vec), &size));
1385:     if (size > 1) defaultType = VECMPI;
1386:     else defaultType = VECSEQ;
1387:   }

1389:   PetscCall(VecRegisterAll());
1390:   PetscCall(PetscOptionsFList("-vec_type", "Vector type", "VecSetType", VecList, defaultType, typeName, 256, &opt));
1391:   if (opt) {
1392:     PetscCall(VecSetType(vec, typeName));
1393:   } else {
1394:     PetscCall(VecSetType(vec, defaultType));
1395:   }
1396:   PetscFunctionReturn(PETSC_SUCCESS);
1397: }

1399: /*@
1400:   VecSetFromOptions - Configures the vector from the options database.

1402:   Collective

1404:   Input Parameter:
1405: . vec - The vector

1407:   Level: beginner

1409:   Notes:
1410:   To see all options, run your program with the -help option.

1412:   Must be called after `VecCreate()` but before the vector is used.

1414: .seealso: [](ch_vectors), `Vec`, `VecCreate()`, `VecSetOptionsPrefix()`
1415: @*/
1416: PetscErrorCode VecSetFromOptions(Vec vec)
1417: {
1418:   PetscBool flg;
1419:   PetscInt  bind_below = 0;

1421:   PetscFunctionBegin;

1424:   PetscObjectOptionsBegin((PetscObject)vec);
1425:   /* Handle vector type options */
1426:   PetscCall(VecSetTypeFromOptions_Private(vec, PetscOptionsObject));

1428:   /* Handle specific vector options */
1429:   PetscTryTypeMethod(vec, setfromoptions, PetscOptionsObject);

1431:   /* Bind to CPU if below a user-specified size threshold.
1432:    * This perhaps belongs in the options for the GPU Vec types, but VecBindToCPU() does nothing when called on non-GPU types,
1433:    * and putting it here makes is more maintainable than duplicating this for all. */
1434:   PetscCall(PetscOptionsInt("-vec_bind_below", "Set the size threshold (in local entries) below which the Vec is bound to the CPU", "VecBindToCPU", bind_below, &bind_below, &flg));
1435:   if (flg && vec->map->n < bind_below) PetscCall(VecBindToCPU(vec, PETSC_TRUE));

1437:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
1438:   PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)vec, PetscOptionsObject));
1439:   PetscOptionsEnd();
1440:   PetscFunctionReturn(PETSC_SUCCESS);
1441: }

1443: /*@
1444:   VecSetSizes - Sets the local and global sizes, and checks to determine compatibility of the sizes

1446:   Collective

1448:   Input Parameters:
1449: + v - the vector
1450: . n - the local size (or `PETSC_DECIDE` to have it set)
1451: - N - the global size (or `PETSC_DETERMINE` to have it set)

1453:   Level: intermediate

1455:   Notes:
1456:   `N` cannot be `PETSC_DETERMINE` if `n` is `PETSC_DECIDE`

1458:   If one processor calls this with `N` of `PETSC_DETERMINE` then all processors must, otherwise the program will hang.

1460: .seealso: [](ch_vectors), `Vec`, `VecGetSize()`, `PetscSplitOwnership()`
1461: @*/
1462: PetscErrorCode VecSetSizes(Vec v, PetscInt n, PetscInt N)
1463: {
1464:   PetscFunctionBegin;
1466:   if (N >= 0) {
1468:     PetscCheck(n <= N, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Local size %" PetscInt_FMT " cannot be larger than global size %" PetscInt_FMT, n, N);
1469:   }
1470:   PetscCheck(!(v->map->n >= 0 || v->map->N >= 0) || !(v->map->n != n || v->map->N != N), PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot change/reset vector sizes to %" PetscInt_FMT " local %" PetscInt_FMT " global after previously setting them to %" PetscInt_FMT " local %" PetscInt_FMT " global", n, N,
1471:              v->map->n, v->map->N);
1472:   v->map->n = n;
1473:   v->map->N = N;
1474:   PetscTryTypeMethod(v, create);
1475:   v->ops->create = NULL;
1476:   PetscFunctionReturn(PETSC_SUCCESS);
1477: }

1479: /*@
1480:   VecSetBlockSize - Sets the block size for future calls to `VecSetValuesBlocked()`
1481:   and `VecSetValuesBlockedLocal()`.

1483:   Logically Collective

1485:   Input Parameters:
1486: + v  - the vector
1487: - bs - the blocksize

1489:   Level: advanced

1491:   Note:
1492:   All vectors obtained by `VecDuplicate()` inherit the same blocksize.

1494:   Vectors obtained with `DMCreateGlobalVector()` and `DMCreateLocalVector()` generally already have a blocksize set based on the state of the `DM`

1496: .seealso: [](ch_vectors), `Vec`, `VecSetValuesBlocked()`, `VecSetLocalToGlobalMapping()`, `VecGetBlockSize()`
1497: @*/
1498: PetscErrorCode VecSetBlockSize(Vec v, PetscInt bs)
1499: {
1500:   PetscFunctionBegin;
1503:   PetscCall(PetscLayoutSetBlockSize(v->map, bs));
1504:   v->bstash.bs = bs; /* use the same blocksize for the vec's block-stash */
1505:   PetscFunctionReturn(PETSC_SUCCESS);
1506: }

1508: /*@
1509:   VecGetBlockSize - Gets the blocksize for the vector, i.e. what is used for `VecSetValuesBlocked()`
1510:   and `VecSetValuesBlockedLocal()`.

1512:   Not Collective

1514:   Input Parameter:
1515: . v - the vector

1517:   Output Parameter:
1518: . bs - the blocksize

1520:   Level: advanced

1522:   Note:
1523:   All vectors obtained by `VecDuplicate()` inherit the same blocksize.

1525: .seealso: [](ch_vectors), `Vec`, `VecSetValuesBlocked()`, `VecSetLocalToGlobalMapping()`, `VecSetBlockSize()`
1526: @*/
1527: PetscErrorCode VecGetBlockSize(Vec v, PetscInt *bs)
1528: {
1529:   PetscFunctionBegin;
1531:   PetscAssertPointer(bs, 2);
1532:   PetscCall(PetscLayoutGetBlockSize(v->map, bs));
1533:   PetscFunctionReturn(PETSC_SUCCESS);
1534: }

1536: /*@C
1537:   VecSetOptionsPrefix - Sets the prefix used for searching for all
1538:   `Vec` options in the database.

1540:   Logically Collective

1542:   Input Parameters:
1543: + v      - the `Vec` context
1544: - prefix - the prefix to prepend to all option names

1546:   Level: advanced

1548:   Note:
1549:   A hyphen (-) must NOT be given at the beginning of the prefix name.
1550:   The first character of all runtime options is AUTOMATICALLY the hyphen.

1552: .seealso: [](ch_vectors), `Vec`, `VecSetFromOptions()`
1553: @*/
1554: PetscErrorCode VecSetOptionsPrefix(Vec v, const char prefix[])
1555: {
1556:   PetscFunctionBegin;
1558:   PetscCall(PetscObjectSetOptionsPrefix((PetscObject)v, prefix));
1559:   PetscFunctionReturn(PETSC_SUCCESS);
1560: }

1562: /*@C
1563:   VecAppendOptionsPrefix - Appends to the prefix used for searching for all
1564:   `Vec` options in the database.

1566:   Logically Collective

1568:   Input Parameters:
1569: + v      - the `Vec` context
1570: - prefix - the prefix to prepend to all option names

1572:   Level: advanced

1574:   Note:
1575:   A hyphen (-) must NOT be given at the beginning of the prefix name.
1576:   The first character of all runtime options is AUTOMATICALLY the hyphen.

1578: .seealso: [](ch_vectors), `Vec`, `VecGetOptionsPrefix()`
1579: @*/
1580: PetscErrorCode VecAppendOptionsPrefix(Vec v, const char prefix[])
1581: {
1582:   PetscFunctionBegin;
1584:   PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)v, prefix));
1585:   PetscFunctionReturn(PETSC_SUCCESS);
1586: }

1588: /*@C
1589:   VecGetOptionsPrefix - Sets the prefix used for searching for all
1590:   Vec options in the database.

1592:   Not Collective

1594:   Input Parameter:
1595: . v - the `Vec` context

1597:   Output Parameter:
1598: . prefix - pointer to the prefix string used

1600:   Level: advanced

1602:   Fortran Notes:
1603:   The user must pass in a string `prefix` of
1604:   sufficient length to hold the prefix.

1606: .seealso: [](ch_vectors), `Vec`, `VecAppendOptionsPrefix()`
1607: @*/
1608: PetscErrorCode VecGetOptionsPrefix(Vec v, const char *prefix[])
1609: {
1610:   PetscFunctionBegin;
1612:   PetscCall(PetscObjectGetOptionsPrefix((PetscObject)v, prefix));
1613:   PetscFunctionReturn(PETSC_SUCCESS);
1614: }

1616: /*@
1617:   VecSetUp - Sets up the internal vector data structures for the later use.

1619:   Collective

1621:   Input Parameter:
1622: . v - the `Vec` context

1624:   Level: advanced

1626:   Notes:
1627:   For basic use of the `Vec` classes the user need not explicitly call
1628:   `VecSetUp()`, since these actions will happen automatically.

1630: .seealso: [](ch_vectors), `Vec`, `VecCreate()`, `VecDestroy()`
1631: @*/
1632: PetscErrorCode VecSetUp(Vec v)
1633: {
1634:   PetscMPIInt size;

1636:   PetscFunctionBegin;
1638:   PetscCheck(v->map->n >= 0 || v->map->N >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Sizes not set");
1639:   if (!((PetscObject)v)->type_name) {
1640:     PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)v), &size));
1641:     if (size == 1) {
1642:       PetscCall(VecSetType(v, VECSEQ));
1643:     } else {
1644:       PetscCall(VecSetType(v, VECMPI));
1645:     }
1646:   }
1647:   PetscFunctionReturn(PETSC_SUCCESS);
1648: }

1650: /*
1651:     These currently expose the PetscScalar/PetscReal in updating the
1652:     cached norm. If we push those down into the implementation these
1653:     will become independent of PetscScalar/PetscReal
1654: */

1656: PetscErrorCode VecCopyAsync_Private(Vec x, Vec y, PetscDeviceContext dctx)
1657: {
1658:   PetscBool flgs[4];
1659:   PetscReal norms[4] = {0.0, 0.0, 0.0, 0.0};

1661:   PetscFunctionBegin;
1666:   if (x == y) PetscFunctionReturn(PETSC_SUCCESS);
1667:   VecCheckSameLocalSize(x, 1, y, 2);
1668:   VecCheckAssembled(x);
1669:   PetscCall(VecSetErrorIfLocked(y, 2));

1671: #if !defined(PETSC_USE_MIXED_PRECISION)
1672:   for (PetscInt i = 0; i < 4; i++) PetscCall(PetscObjectComposedDataGetReal((PetscObject)x, NormIds[i], norms[i], flgs[i]));
1673: #endif

1675:   PetscCall(PetscLogEventBegin(VEC_Copy, x, y, 0, 0));
1676: #if defined(PETSC_USE_MIXED_PRECISION)
1677:   extern PetscErrorCode VecGetArray(Vec, double **);
1678:   extern PetscErrorCode VecRestoreArray(Vec, double **);
1679:   extern PetscErrorCode VecGetArray(Vec, float **);
1680:   extern PetscErrorCode VecRestoreArray(Vec, float **);
1681:   extern PetscErrorCode VecGetArrayRead(Vec, const double **);
1682:   extern PetscErrorCode VecRestoreArrayRead(Vec, const double **);
1683:   extern PetscErrorCode VecGetArrayRead(Vec, const float **);
1684:   extern PetscErrorCode VecRestoreArrayRead(Vec, const float **);
1685:   if ((((PetscObject)x)->precision == PETSC_PRECISION_SINGLE) && (((PetscObject)y)->precision == PETSC_PRECISION_DOUBLE)) {
1686:     PetscInt     i, n;
1687:     const float *xx;
1688:     double      *yy;
1689:     PetscCall(VecGetArrayRead(x, &xx));
1690:     PetscCall(VecGetArray(y, &yy));
1691:     PetscCall(VecGetLocalSize(x, &n));
1692:     for (i = 0; i < n; i++) yy[i] = xx[i];
1693:     PetscCall(VecRestoreArrayRead(x, &xx));
1694:     PetscCall(VecRestoreArray(y, &yy));
1695:   } else if ((((PetscObject)x)->precision == PETSC_PRECISION_DOUBLE) && (((PetscObject)y)->precision == PETSC_PRECISION_SINGLE)) {
1696:     PetscInt      i, n;
1697:     float        *yy;
1698:     const double *xx;
1699:     PetscCall(VecGetArrayRead(x, &xx));
1700:     PetscCall(VecGetArray(y, &yy));
1701:     PetscCall(VecGetLocalSize(x, &n));
1702:     for (i = 0; i < n; i++) yy[i] = (float)xx[i];
1703:     PetscCall(VecRestoreArrayRead(x, &xx));
1704:     PetscCall(VecRestoreArray(y, &yy));
1705:   } else PetscUseTypeMethod(x, copy, y);
1706: #else
1707:   VecMethodDispatch(x, dctx, VecAsyncFnName(Copy), copy, (Vec, Vec, PetscDeviceContext), y);
1708: #endif

1710:   PetscCall(PetscObjectStateIncrease((PetscObject)y));
1711: #if !defined(PETSC_USE_MIXED_PRECISION)
1712:   for (PetscInt i = 0; i < 4; i++) {
1713:     if (flgs[i]) PetscCall(PetscObjectComposedDataSetReal((PetscObject)y, NormIds[i], norms[i]));
1714:   }
1715: #endif

1717:   PetscCall(PetscLogEventEnd(VEC_Copy, x, y, 0, 0));
1718:   PetscFunctionReturn(PETSC_SUCCESS);
1719: }

1721: /*@
1722:   VecCopy - Copies a vector `y = x`

1724:   Logically Collective

1726:   Input Parameter:
1727: . x - the vector

1729:   Output Parameter:
1730: . y - the copy

1732:   Level: beginner

1734:   Note:
1735:   For default parallel PETSc vectors, both `x` and `y` must be distributed in
1736:   the same manner; local copies are done.

1738:   Developer Notes:
1739:   `PetscCheckSameTypeAndComm`(x,1,y,2) is not used on these vectors because we allow one
1740:   of the vectors to be sequential and one to be parallel so long as both have the same
1741:   local sizes. This is used in some internal functions in PETSc.

1743: .seealso: [](ch_vectors), `Vec`, `VecDuplicate()`
1744: @*/
1745: PetscErrorCode VecCopy(Vec x, Vec y)
1746: {
1747:   PetscFunctionBegin;
1748:   PetscCall(VecCopyAsync_Private(x, y, NULL));
1749:   PetscFunctionReturn(PETSC_SUCCESS);
1750: }

1752: PetscErrorCode VecSwapAsync_Private(Vec x, Vec y, PetscDeviceContext dctx)
1753: {
1754:   PetscReal normxs[4], normys[4];
1755:   PetscBool flgxs[4], flgys[4];

1757:   PetscFunctionBegin;
1762:   PetscCheckSameTypeAndComm(x, 1, y, 2);
1763:   VecCheckSameSize(x, 1, y, 2);
1764:   VecCheckAssembled(x);
1765:   VecCheckAssembled(y);
1766:   PetscCall(VecSetErrorIfLocked(x, 1));
1767:   PetscCall(VecSetErrorIfLocked(y, 2));

1769:   for (PetscInt i = 0; i < 4; i++) {
1770:     PetscCall(PetscObjectComposedDataGetReal((PetscObject)x, NormIds[i], normxs[i], flgxs[i]));
1771:     PetscCall(PetscObjectComposedDataGetReal((PetscObject)y, NormIds[i], normys[i], flgys[i]));
1772:   }

1774:   PetscCall(PetscLogEventBegin(VEC_Swap, x, y, 0, 0));
1775:   VecMethodDispatch(x, dctx, VecAsyncFnName(Swap), swap, (Vec, Vec, PetscDeviceContext), y);
1776:   PetscCall(PetscLogEventEnd(VEC_Swap, x, y, 0, 0));

1778:   PetscCall(PetscObjectStateIncrease((PetscObject)x));
1779:   PetscCall(PetscObjectStateIncrease((PetscObject)y));
1780:   for (PetscInt i = 0; i < 4; i++) {
1781:     if (flgxs[i]) PetscCall(PetscObjectComposedDataSetReal((PetscObject)y, NormIds[i], normxs[i]));
1782:     if (flgys[i]) PetscCall(PetscObjectComposedDataSetReal((PetscObject)x, NormIds[i], normys[i]));
1783:   }
1784:   PetscFunctionReturn(PETSC_SUCCESS);
1785: }
1786: /*@
1787:   VecSwap - Swaps the values between two vectors, `x` and `y`.

1789:   Logically Collective

1791:   Input Parameters:
1792: + x - the first vector
1793: - y - the second vector

1795:   Level: advanced

1797: .seealso: [](ch_vectors), `Vec`, `VecSet()`
1798: @*/
1799: PetscErrorCode VecSwap(Vec x, Vec y)
1800: {
1801:   PetscFunctionBegin;
1802:   PetscCall(VecSwapAsync_Private(x, y, NULL));
1803:   PetscFunctionReturn(PETSC_SUCCESS);
1804: }

1806: /*@C
1807:   VecStashViewFromOptions - Processes command line options to determine if/how a `VecStash` object is to be viewed.

1809:   Collective

1811:   Input Parameters:
1812: + obj        - the `Vec` containing a stash
1813: . bobj       - optional other object that provides the prefix
1814: - optionname - option to activate viewing

1816:   Level: intermediate

1818:   Developer Notes:
1819:   This cannot use `PetscObjectViewFromOptions()` because it takes a `Vec` as an argument but does not use `VecView()`

1821: .seealso: [](ch_vectors), `Vec`, `VecStashSetInitialSize()`
1822: @*/
1823: PetscErrorCode VecStashViewFromOptions(Vec obj, PetscObject bobj, const char optionname[])
1824: {
1825:   PetscViewer       viewer;
1826:   PetscBool         flg;
1827:   PetscViewerFormat format;
1828:   char             *prefix;

1830:   PetscFunctionBegin;
1831:   prefix = bobj ? bobj->prefix : ((PetscObject)obj)->prefix;
1832:   PetscCall(PetscOptionsGetViewer(PetscObjectComm((PetscObject)obj), ((PetscObject)obj)->options, prefix, optionname, &viewer, &format, &flg));
1833:   if (flg) {
1834:     PetscCall(PetscViewerPushFormat(viewer, format));
1835:     PetscCall(VecStashView(obj, viewer));
1836:     PetscCall(PetscViewerPopFormat(viewer));
1837:     PetscCall(PetscOptionsRestoreViewer(&viewer));
1838:   }
1839:   PetscFunctionReturn(PETSC_SUCCESS);
1840: }

1842: /*@
1843:   VecStashView - Prints the entries in the vector stash and block stash.

1845:   Collective

1847:   Input Parameters:
1848: + v      - the vector
1849: - viewer - the viewer

1851:   Level: advanced

1853: .seealso: [](ch_vectors), `Vec`, `VecSetBlockSize()`, `VecSetValues()`, `VecSetValuesBlocked()`
1854: @*/
1855: PetscErrorCode VecStashView(Vec v, PetscViewer viewer)
1856: {
1857:   PetscMPIInt rank;
1858:   PetscInt    i, j;
1859:   PetscBool   match;
1860:   VecStash   *s;
1861:   PetscScalar val;

1863:   PetscFunctionBegin;
1866:   PetscCheckSameComm(v, 1, viewer, 2);

1868:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &match));
1869:   PetscCheck(match, PETSC_COMM_SELF, PETSC_ERR_SUP, "Stash viewer only works with ASCII viewer not %s", ((PetscObject)v)->type_name);
1870:   PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
1871:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)v), &rank));
1872:   s = &v->bstash;

1874:   /* print block stash */
1875:   PetscCall(PetscViewerASCIIPushSynchronized(viewer));
1876:   PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d]Vector Block stash size %" PetscInt_FMT " block size %" PetscInt_FMT "\n", rank, s->n, s->bs));
1877:   for (i = 0; i < s->n; i++) {
1878:     PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] Element %" PetscInt_FMT " ", rank, s->idx[i]));
1879:     for (j = 0; j < s->bs; j++) {
1880:       val = s->array[i * s->bs + j];
1881: #if defined(PETSC_USE_COMPLEX)
1882:       PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "(%18.16e %18.16e) ", (double)PetscRealPart(val), (double)PetscImaginaryPart(val)));
1883: #else
1884:       PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "%18.16e ", (double)val));
1885: #endif
1886:     }
1887:     PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "\n"));
1888:   }
1889:   PetscCall(PetscViewerFlush(viewer));

1891:   s = &v->stash;

1893:   /* print basic stash */
1894:   PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d]Vector stash size %" PetscInt_FMT "\n", rank, s->n));
1895:   for (i = 0; i < s->n; i++) {
1896:     val = s->array[i];
1897: #if defined(PETSC_USE_COMPLEX)
1898:     PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] Element %" PetscInt_FMT " (%18.16e %18.16e) ", rank, s->idx[i], (double)PetscRealPart(val), (double)PetscImaginaryPart(val)));
1899: #else
1900:     PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "[%d] Element %" PetscInt_FMT " %18.16e\n", rank, s->idx[i], (double)val));
1901: #endif
1902:   }
1903:   PetscCall(PetscViewerFlush(viewer));
1904:   PetscCall(PetscViewerASCIIPopSynchronized(viewer));
1905:   PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
1906:   PetscFunctionReturn(PETSC_SUCCESS);
1907: }

1909: PetscErrorCode PetscOptionsGetVec(PetscOptions options, const char prefix[], const char key[], Vec v, PetscBool *set)
1910: {
1911:   PetscInt     i, N, rstart, rend;
1912:   PetscScalar *xx;
1913:   PetscReal   *xreal;
1914:   PetscBool    iset;

1916:   PetscFunctionBegin;
1917:   PetscCall(VecGetOwnershipRange(v, &rstart, &rend));
1918:   PetscCall(VecGetSize(v, &N));
1919:   PetscCall(PetscCalloc1(N, &xreal));
1920:   PetscCall(PetscOptionsGetRealArray(options, prefix, key, xreal, &N, &iset));
1921:   if (iset) {
1922:     PetscCall(VecGetArray(v, &xx));
1923:     for (i = rstart; i < rend; i++) xx[i - rstart] = xreal[i];
1924:     PetscCall(VecRestoreArray(v, &xx));
1925:   }
1926:   PetscCall(PetscFree(xreal));
1927:   if (set) *set = iset;
1928:   PetscFunctionReturn(PETSC_SUCCESS);
1929: }

1931: /*@
1932:   VecGetLayout - get `PetscLayout` describing a vector layout

1934:   Not Collective

1936:   Input Parameter:
1937: . x - the vector

1939:   Output Parameter:
1940: . map - the layout

1942:   Level: developer

1944:   Note:
1945:   The layout determines what vector elements are contained on each MPI process

1947: .seealso: [](ch_vectors), `PetscLayout`, `Vec`, `VecGetSizes()`, `VecGetOwnershipRange()`, `VecGetOwnershipRanges()`
1948: @*/
1949: PetscErrorCode VecGetLayout(Vec x, PetscLayout *map)
1950: {
1951:   PetscFunctionBegin;
1953:   PetscAssertPointer(map, 2);
1954:   *map = x->map;
1955:   PetscFunctionReturn(PETSC_SUCCESS);
1956: }

1958: /*@
1959:   VecSetLayout - set `PetscLayout` describing vector layout

1961:   Not Collective

1963:   Input Parameters:
1964: + x   - the vector
1965: - map - the layout

1967:   Level: developer

1969:   Note:
1970:   It is normally only valid to replace the layout with a layout known to be equivalent.

1972: .seealso: [](ch_vectors), `Vec`, `PetscLayout`, `VecGetLayout()`, `VecGetSizes()`, `VecGetOwnershipRange()`, `VecGetOwnershipRanges()`
1973: @*/
1974: PetscErrorCode VecSetLayout(Vec x, PetscLayout map)
1975: {
1976:   PetscFunctionBegin;
1978:   PetscCall(PetscLayoutReference(map, &x->map));
1979:   PetscFunctionReturn(PETSC_SUCCESS);
1980: }

1982: PetscErrorCode VecSetInf(Vec xin)
1983: {
1984:   // use of variables one and zero over just doing 1.0/0.0 is deliberate. MSVC complains that
1985:   // we are dividing by zero in the latter case (ostensibly because dividing by 0 is UB, but
1986:   // only for *integers* not floats).
1987:   const PetscScalar one = 1.0, zero = 0.0;
1988:   PetscScalar       inf;

1990:   PetscFunctionBegin;
1991:   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
1992:   inf = one / zero;
1993:   PetscCall(PetscFPTrapPop());
1994:   if (xin->ops->set) {
1995:     PetscUseTypeMethod(xin, set, inf);
1996:   } else {
1997:     PetscInt     n;
1998:     PetscScalar *xx;

2000:     PetscCall(VecGetLocalSize(xin, &n));
2001:     PetscCall(VecGetArrayWrite(xin, &xx));
2002:     for (PetscInt i = 0; i < n; ++i) xx[i] = inf;
2003:     PetscCall(VecRestoreArrayWrite(xin, &xx));
2004:   }
2005:   PetscFunctionReturn(PETSC_SUCCESS);
2006: }

2008: /*@
2009:   VecBindToCPU - marks a vector to temporarily stay on the CPU and perform computations on the CPU

2011:   Logically collective

2013:   Input Parameters:
2014: + v   - the vector
2015: - flg - bind to the CPU if value of `PETSC_TRUE`

2017:   Level: intermediate

2019: .seealso: [](ch_vectors), `Vec`, `VecBoundToCPU()`
2020: @*/
2021: PetscErrorCode VecBindToCPU(Vec v, PetscBool flg)
2022: {
2023:   PetscFunctionBegin;
2026: #if defined(PETSC_HAVE_DEVICE)
2027:   if (v->boundtocpu == flg) PetscFunctionReturn(PETSC_SUCCESS);
2028:   v->boundtocpu = flg;
2029:   PetscTryTypeMethod(v, bindtocpu, flg);
2030: #endif
2031:   PetscFunctionReturn(PETSC_SUCCESS);
2032: }

2034: /*@
2035:   VecBoundToCPU - query if a vector is bound to the CPU

2037:   Not collective

2039:   Input Parameter:
2040: . v - the vector

2042:   Output Parameter:
2043: . flg - the logical flag

2045:   Level: intermediate

2047: .seealso: [](ch_vectors), `Vec`, `VecBindToCPU()`
2048: @*/
2049: PetscErrorCode VecBoundToCPU(Vec v, PetscBool *flg)
2050: {
2051:   PetscFunctionBegin;
2053:   PetscAssertPointer(flg, 2);
2054: #if defined(PETSC_HAVE_DEVICE)
2055:   *flg = v->boundtocpu;
2056: #else
2057:   *flg = PETSC_TRUE;
2058: #endif
2059:   PetscFunctionReturn(PETSC_SUCCESS);
2060: }

2062: /*@
2063:   VecSetBindingPropagates - Sets whether the state of being bound to the CPU for a GPU vector type propagates to child and some other associated objects

2065:   Input Parameters:
2066: + v   - the vector
2067: - flg - flag indicating whether the boundtocpu flag should be propagated

2069:   Level: developer

2071:   Notes:
2072:   If the value of flg is set to true, then `VecDuplicate()` and `VecDuplicateVecs()` will bind created vectors to GPU if the input vector is bound to the CPU.
2073:   The created vectors will also have their bindingpropagates flag set to true.

2075:   Developer Notes:
2076:   If a `DMDA` has the `-dm_bind_below option` set to true, then vectors created by `DMCreateGlobalVector()` will have `VecSetBindingPropagates()` called on them to
2077:   set their bindingpropagates flag to true.

2079: .seealso: [](ch_vectors), `Vec`, `MatSetBindingPropagates()`, `VecGetBindingPropagates()`
2080: @*/
2081: PetscErrorCode VecSetBindingPropagates(Vec v, PetscBool flg)
2082: {
2083:   PetscFunctionBegin;
2085: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
2086:   v->bindingpropagates = flg;
2087: #endif
2088:   PetscFunctionReturn(PETSC_SUCCESS);
2089: }

2091: /*@
2092:   VecGetBindingPropagates - Gets whether the state of being bound to the CPU for a GPU vector type propagates to child and some other associated objects

2094:   Input Parameter:
2095: . v - the vector

2097:   Output Parameter:
2098: . flg - flag indicating whether the boundtocpu flag will be propagated

2100:   Level: developer

2102: .seealso: [](ch_vectors), `Vec`, `VecSetBindingPropagates()`
2103: @*/
2104: PetscErrorCode VecGetBindingPropagates(Vec v, PetscBool *flg)
2105: {
2106:   PetscFunctionBegin;
2108:   PetscAssertPointer(flg, 2);
2109: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
2110:   *flg = v->bindingpropagates;
2111: #else
2112:   *flg = PETSC_FALSE;
2113: #endif
2114:   PetscFunctionReturn(PETSC_SUCCESS);
2115: }

2117: /*@C
2118:   VecSetPinnedMemoryMin - Set the minimum data size for which pinned memory will be used for host (CPU) allocations.

2120:   Logically Collective

2122:   Input Parameters:
2123: + v      - the vector
2124: - mbytes - minimum data size in bytes

2126:   Options Database Key:
2127: . -vec_pinned_memory_min <size> - minimum size (in bytes) for an allocation to use pinned memory on host.

2129:   Level: developer

2131:   Note:
2132:   Specifying -1 ensures that pinned memory will never be used.

2134: .seealso: [](ch_vectors), `Vec`, `VecGetPinnedMemoryMin()`
2135: @*/
2136: PetscErrorCode VecSetPinnedMemoryMin(Vec v, size_t mbytes)
2137: {
2138:   PetscFunctionBegin;
2140: #if PetscDefined(HAVE_DEVICE)
2141:   v->minimum_bytes_pinned_memory = mbytes;
2142: #endif
2143:   PetscFunctionReturn(PETSC_SUCCESS);
2144: }

2146: /*@C
2147:   VecGetPinnedMemoryMin - Get the minimum data size for which pinned memory will be used for host (CPU) allocations.

2149:   Logically Collective

2151:   Input Parameter:
2152: . v - the vector

2154:   Output Parameter:
2155: . mbytes - minimum data size in bytes

2157:   Level: developer

2159: .seealso: [](ch_vectors), `Vec`, `VecSetPinnedMemoryMin()`
2160: @*/
2161: PetscErrorCode VecGetPinnedMemoryMin(Vec v, size_t *mbytes)
2162: {
2163:   PetscFunctionBegin;
2165:   PetscAssertPointer(mbytes, 2);
2166: #if PetscDefined(HAVE_DEVICE)
2167:   *mbytes = v->minimum_bytes_pinned_memory;
2168: #endif
2169:   PetscFunctionReturn(PETSC_SUCCESS);
2170: }

2172: /*@
2173:   VecGetOffloadMask - Get the offload mask of a `Vec`

2175:   Not Collective

2177:   Input Parameter:
2178: . v - the vector

2180:   Output Parameter:
2181: . mask - corresponding `PetscOffloadMask` enum value.

2183:   Level: intermediate

2185: .seealso: [](ch_vectors), `Vec`, `VecCreateSeqCUDA()`, `VecCreateSeqViennaCL()`, `VecGetArray()`, `VecGetType()`
2186: @*/
2187: PetscErrorCode VecGetOffloadMask(Vec v, PetscOffloadMask *mask)
2188: {
2189:   PetscFunctionBegin;
2191:   PetscAssertPointer(mask, 2);
2192:   *mask = v->offloadmask;
2193:   PetscFunctionReturn(PETSC_SUCCESS);
2194: }

2196: #if !defined(PETSC_HAVE_VIENNACL)
2197: PETSC_EXTERN PetscErrorCode VecViennaCLGetCLContext(Vec v, PETSC_UINTPTR_T *ctx)
2198: {
2199:   SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "PETSc must be configured with --with-opencl to get a Vec's cl_context");
2200: }

2202: PETSC_EXTERN PetscErrorCode VecViennaCLGetCLQueue(Vec v, PETSC_UINTPTR_T *queue)
2203: {
2204:   SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "PETSc must be configured with --with-opencl to get a Vec's cl_command_queue");
2205: }

2207: PETSC_EXTERN PetscErrorCode VecViennaCLGetCLMem(Vec v, PETSC_UINTPTR_T *queue)
2208: {
2209:   SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "PETSc must be configured with --with-opencl to get a Vec's cl_mem");
2210: }

2212: PETSC_EXTERN PetscErrorCode VecViennaCLGetCLMemRead(Vec v, PETSC_UINTPTR_T *queue)
2213: {
2214:   SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "PETSc must be configured with --with-opencl to get a Vec's cl_mem");
2215: }

2217: PETSC_EXTERN PetscErrorCode VecViennaCLGetCLMemWrite(Vec v, PETSC_UINTPTR_T *queue)
2218: {
2219:   SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "PETSc must be configured with --with-opencl to get a Vec's cl_mem");
2220: }

2222: PETSC_EXTERN PetscErrorCode VecViennaCLRestoreCLMemWrite(Vec v)
2223: {
2224:   SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "PETSc must be configured with --with-opencl to restore a Vec's cl_mem");
2225: }
2226: #endif

2228: static PetscErrorCode VecErrorWeightedNorms_Basic(Vec U, Vec Y, Vec E, NormType wnormtype, PetscReal atol, Vec vatol, PetscReal rtol, Vec vrtol, PetscReal ignore_max, PetscReal *norm, PetscInt *norm_loc, PetscReal *norma, PetscInt *norma_loc, PetscReal *normr, PetscInt *normr_loc)
2229: {
2230:   const PetscScalar *u, *y;
2231:   const PetscScalar *atola = NULL, *rtola = NULL, *erra = NULL;
2232:   PetscInt           n, n_loc = 0, na_loc = 0, nr_loc = 0;
2233:   PetscReal          nrm = 0, nrma = 0, nrmr = 0, err_loc[6];

2235:   PetscFunctionBegin;
2236: #define SkipSmallValue(a, b, tol) \
2237:   if (PetscAbsScalar(a) < tol || PetscAbsScalar(b) < tol) continue

2239:   PetscCall(VecGetLocalSize(U, &n));
2240:   PetscCall(VecGetArrayRead(U, &u));
2241:   PetscCall(VecGetArrayRead(Y, &y));
2242:   if (E) PetscCall(VecGetArrayRead(E, &erra));
2243:   if (vatol) PetscCall(VecGetArrayRead(vatol, &atola));
2244:   if (vrtol) PetscCall(VecGetArrayRead(vrtol, &rtola));
2245:   for (PetscInt i = 0; i < n; i++) {
2246:     PetscReal err, tol, tola, tolr;

2248:     SkipSmallValue(y[i], u[i], ignore_max);
2249:     atol = atola ? PetscRealPart(atola[i]) : atol;
2250:     rtol = rtola ? PetscRealPart(rtola[i]) : rtol;
2251:     err  = erra ? PetscAbsScalar(erra[i]) : PetscAbsScalar(y[i] - u[i]);
2252:     tola = atol;
2253:     tolr = rtol * PetscMax(PetscAbsScalar(u[i]), PetscAbsScalar(y[i]));
2254:     tol  = tola + tolr;
2255:     if (tola > 0.) {
2256:       if (wnormtype == NORM_INFINITY) nrma = PetscMax(nrma, err / tola);
2257:       else nrma += PetscSqr(err / tola);
2258:       na_loc++;
2259:     }
2260:     if (tolr > 0.) {
2261:       if (wnormtype == NORM_INFINITY) nrmr = PetscMax(nrmr, err / tolr);
2262:       else nrmr += PetscSqr(err / tolr);
2263:       nr_loc++;
2264:     }
2265:     if (tol > 0.) {
2266:       if (wnormtype == NORM_INFINITY) nrm = PetscMax(nrm, err / tol);
2267:       else nrm += PetscSqr(err / tol);
2268:       n_loc++;
2269:     }
2270:   }
2271:   if (E) PetscCall(VecRestoreArrayRead(E, &erra));
2272:   if (vatol) PetscCall(VecRestoreArrayRead(vatol, &atola));
2273:   if (vrtol) PetscCall(VecRestoreArrayRead(vrtol, &rtola));
2274:   PetscCall(VecRestoreArrayRead(U, &u));
2275:   PetscCall(VecRestoreArrayRead(Y, &y));
2276: #undef SkipSmallValue

2278:   err_loc[0] = nrm;
2279:   err_loc[1] = nrma;
2280:   err_loc[2] = nrmr;
2281:   err_loc[3] = (PetscReal)n_loc;
2282:   err_loc[4] = (PetscReal)na_loc;
2283:   err_loc[5] = (PetscReal)nr_loc;
2284:   if (wnormtype == NORM_2) {
2285:     PetscCall(MPIU_Allreduce(MPI_IN_PLACE, err_loc, 6, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)U)));
2286:   } else {
2287:     PetscCall(MPIU_Allreduce(MPI_IN_PLACE, err_loc, 3, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)U)));
2288:     PetscCall(MPIU_Allreduce(MPI_IN_PLACE, err_loc + 3, 3, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)U)));
2289:   }
2290:   if (wnormtype == NORM_2) {
2291:     *norm  = PetscSqrtReal(err_loc[0]);
2292:     *norma = PetscSqrtReal(err_loc[1]);
2293:     *normr = PetscSqrtReal(err_loc[2]);
2294:   } else {
2295:     *norm  = err_loc[0];
2296:     *norma = err_loc[1];
2297:     *normr = err_loc[2];
2298:   }
2299:   *norm_loc  = (PetscInt)err_loc[3];
2300:   *norma_loc = (PetscInt)err_loc[4];
2301:   *normr_loc = (PetscInt)err_loc[5];
2302:   PetscFunctionReturn(PETSC_SUCCESS);
2303: }

2305: /*@
2306:   VecErrorWeightedNorms - compute a weighted norm of the difference between two vectors

2308:   Collective

2310:   Input Parameters:
2311: + U          - first vector to be compared
2312: . Y          - second vector to be compared
2313: . E          - optional third vector representing the error (if not provided, the error is ||U-Y||)
2314: . wnormtype  - norm type
2315: . atol       - scalar for absolute tolerance
2316: . vatol      - vector representing per-entry absolute tolerances (can be ``NULL``)
2317: . rtol       - scalar for relative tolerance
2318: . vrtol      - vector representing per-entry relative tolerances (can be ``NULL``)
2319: - ignore_max - ignore values smaller then this value in absolute terms.

2321:   Output Parameters:
2322: + norm      - weighted norm
2323: . norm_loc  - number of vector locations used for the weighted norm
2324: . norma     - weighted norm based on the absolute tolerance
2325: . norma_loc - number of vector locations used for the absolute weighted norm
2326: . normr     - weighted norm based on the relative tolerance
2327: - normr_loc - number of vector locations used for the relative weighted norm

2329:   Level: developer

2331:   Notes:
2332:   This is primarily used for computing weighted local truncation errors in ``TS``.

2334: .seealso: [](ch_vectors), `Vec`, `NormType`, ``TSErrorWeightedNorm()``, ``TSErrorWeightedENorm()``
2335: @*/
2336: PetscErrorCode VecErrorWeightedNorms(Vec U, Vec Y, Vec E, NormType wnormtype, PetscReal atol, Vec vatol, PetscReal rtol, Vec vrtol, PetscReal ignore_max, PetscReal *norm, PetscInt *norm_loc, PetscReal *norma, PetscInt *norma_loc, PetscReal *normr, PetscInt *normr_loc)
2337: {
2338:   PetscFunctionBegin;
2343:   if (E) {
2346:   }
2349:   if (vatol) {
2352:   }
2354:   if (vrtol) {
2357:   }
2359:   PetscAssertPointer(norm, 10);
2360:   PetscAssertPointer(norm_loc, 11);
2361:   PetscAssertPointer(norma, 12);
2362:   PetscAssertPointer(norma_loc, 13);
2363:   PetscAssertPointer(normr, 14);
2364:   PetscAssertPointer(normr_loc, 15);
2365:   PetscCheck(wnormtype == NORM_2 || wnormtype == NORM_INFINITY, PetscObjectComm((PetscObject)U), PETSC_ERR_SUP, "No support for norm type %s", NormTypes[wnormtype]);

2367:   /* There are potentially 5 vectors involved, some of them may happen to be of different type or bound to cpu.
2368:      Here we check that they all implement the same operation and call it if so.
2369:      Otherwise, we call the _Basic implementation that always works (provided VecGetArrayRead is implemented). */
2370:   PetscBool sameop = (PetscBool)(U->ops->errorwnorm && U->ops->errorwnorm == Y->ops->errorwnorm);
2371:   if (sameop && E) sameop = (PetscBool)(U->ops->errorwnorm == E->ops->errorwnorm);
2372:   if (sameop && vatol) sameop = (PetscBool)(U->ops->errorwnorm == vatol->ops->errorwnorm);
2373:   if (sameop && vrtol) sameop = (PetscBool)(U->ops->errorwnorm == vrtol->ops->errorwnorm);
2374:   if (sameop) PetscUseTypeMethod(U, errorwnorm, Y, E, wnormtype, atol, vatol, rtol, vrtol, ignore_max, norm, norm_loc, norma, norma_loc, normr, normr_loc);
2375:   else PetscCall(VecErrorWeightedNorms_Basic(U, Y, E, wnormtype, atol, vatol, rtol, vrtol, ignore_max, norm, norm_loc, norma, norma_loc, normr, normr_loc));
2376:   PetscFunctionReturn(PETSC_SUCCESS);
2377: }