Actual source code: vinv.c

  1: /*
  2:      Some useful vector utility functions.
  3: */
  4: #include <../src/vec/vec/impls/mpi/pvecimpl.h>

  6: /*@
  7:   VecStrideSet - Sets a subvector of a vector defined
  8:   by a starting point and a stride with a given value

 10:   Logically Collective

 12:   Input Parameters:
 13: + v     - the vector
 14: . start - starting point of the subvector (defined by a stride)
 15: - s     - value to set for each entry in that subvector

 17:   Level: advanced

 19:   Notes:
 20:   One must call `VecSetBlockSize()` before this routine to set the stride
 21:   information, or use a vector created from a multicomponent `DMDA`.

 23:   This will only work if the desire subvector is a stride subvector

 25: .seealso: `Vec`, `VecNorm()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`, `VecStrideScale()`
 26: @*/
 27: PetscErrorCode VecStrideSet(Vec v, PetscInt start, PetscScalar s)
 28: {
 29:   PetscInt     i, n, bs;
 30:   PetscScalar *x;

 32:   PetscFunctionBegin;
 35:   PetscCall(VecGetLocalSize(v, &n));
 36:   PetscCall(VecGetBlockSize(v, &bs));
 37:   PetscCheck(start >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative start %" PetscInt_FMT, start);
 38:   PetscCheck(start < bs, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Start of stride subvector (%" PetscInt_FMT ") is too large for stride. Have you set the vector blocksize (%" PetscInt_FMT ") correctly with VecSetBlockSize()?", start, bs);
 39:   PetscCall(VecGetArray(v, &x));
 40:   for (i = start; i < n; i += bs) x[i] = s;
 41:   PetscCall(VecRestoreArray(v, &x));
 42:   PetscFunctionReturn(PETSC_SUCCESS);
 43: }

 45: /*@
 46:   VecStrideScale - Scales a subvector of a vector defined
 47:   by a starting point and a stride.

 49:   Logically Collective

 51:   Input Parameters:
 52: + v     - the vector
 53: . start - starting point of the subvector (defined by a stride)
 54: - scale - value to multiply each subvector entry by

 56:   Level: advanced

 58:   Notes:
 59:   One must call `VecSetBlockSize()` before this routine to set the stride
 60:   information, or use a vector created from a multicomponent `DMDA`.

 62:   This will only work if the desire subvector is a stride subvector

 64: .seealso: `Vec`, `VecNorm()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`
 65: @*/
 66: PetscErrorCode VecStrideScale(Vec v, PetscInt start, PetscScalar scale)
 67: {
 68:   PetscInt     i, n, bs;
 69:   PetscScalar *x;

 71:   PetscFunctionBegin;
 75:   PetscCall(VecGetLocalSize(v, &n));
 76:   PetscCall(VecGetBlockSize(v, &bs));
 77:   PetscCheck(start >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative start %" PetscInt_FMT, start);
 78:   PetscCheck(start < bs, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Start of stride subvector (%" PetscInt_FMT ") is too large for stride. Have you set the vector blocksize (%" PetscInt_FMT ") correctly with VecSetBlockSize()?", start, bs);
 79:   PetscCall(VecGetArray(v, &x));
 80:   for (i = start; i < n; i += bs) x[i] *= scale;
 81:   PetscCall(VecRestoreArray(v, &x));
 82:   PetscFunctionReturn(PETSC_SUCCESS);
 83: }

 85: /*@
 86:   VecStrideNorm - Computes the norm of subvector of a vector defined
 87:   by a starting point and a stride.

 89:   Collective

 91:   Input Parameters:
 92: + v     - the vector
 93: . start - starting point of the subvector (defined by a stride)
 94: - ntype - type of norm, one of `NORM_1`, `NORM_2`, `NORM_INFINITY`

 96:   Output Parameter:
 97: . nrm - the norm

 99:   Level: advanced

101:   Notes:
102:   One must call `VecSetBlockSize()` before this routine to set the stride
103:   information, or use a vector created from a multicomponent `DMDA`.

105:   If x is the array representing the vector x then this computes the norm
106:   of the array (x[start],x[start+stride],x[start+2*stride], ....)

108:   This is useful for computing, say the norm of the pressure variable when
109:   the pressure is stored (interlaced) with other variables, say density etc.

111:   This will only work if the desire subvector is a stride subvector

113: .seealso: `Vec`, `VecNorm()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`
114: @*/
115: PetscErrorCode VecStrideNorm(Vec v, PetscInt start, NormType ntype, PetscReal *nrm)
116: {
117:   PetscInt           i, n, bs;
118:   const PetscScalar *x;
119:   PetscReal          tnorm;

121:   PetscFunctionBegin;
125:   PetscAssertPointer(nrm, 4);
126:   PetscCall(VecGetLocalSize(v, &n));
127:   PetscCall(VecGetBlockSize(v, &bs));
128:   PetscCheck(start >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative start %" PetscInt_FMT, start);
129:   PetscCheck(start < bs, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Start of stride subvector (%" PetscInt_FMT ") is too large for stride. Have you set the vector blocksize (%" PetscInt_FMT ") correctly with VecSetBlockSize()?", start, bs);
130:   PetscCall(VecGetArrayRead(v, &x));
131:   if (ntype == NORM_2) {
132:     PetscScalar sum = 0.0;
133:     for (i = start; i < n; i += bs) sum += x[i] * (PetscConj(x[i]));
134:     tnorm = PetscRealPart(sum);
135:     PetscCall(MPIU_Allreduce(&tnorm, nrm, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)v)));
136:     *nrm = PetscSqrtReal(*nrm);
137:   } else if (ntype == NORM_1) {
138:     tnorm = 0.0;
139:     for (i = start; i < n; i += bs) tnorm += PetscAbsScalar(x[i]);
140:     PetscCall(MPIU_Allreduce(&tnorm, nrm, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)v)));
141:   } else if (ntype == NORM_INFINITY) {
142:     tnorm = 0.0;
143:     for (i = start; i < n; i += bs) {
144:       if (PetscAbsScalar(x[i]) > tnorm) tnorm = PetscAbsScalar(x[i]);
145:     }
146:     PetscCall(MPIU_Allreduce(&tnorm, nrm, 1, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)v)));
147:   } else SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown norm type");
148:   PetscCall(VecRestoreArrayRead(v, &x));
149:   PetscFunctionReturn(PETSC_SUCCESS);
150: }

152: /*@
153:   VecStrideMax - Computes the maximum of subvector of a vector defined
154:   by a starting point and a stride and optionally its location.

156:   Collective

158:   Input Parameters:
159: + v     - the vector
160: - start - starting point of the subvector (defined by a stride)

162:   Output Parameters:
163: + idex - the location where the maximum occurred  (pass `NULL` if not required)
164: - nrm  - the maximum value in the subvector

166:   Level: advanced

168:   Notes:
169:   One must call `VecSetBlockSize()` before this routine to set the stride
170:   information, or use a vector created from a multicomponent `DMDA`.

172:   If xa is the array representing the vector x, then this computes the max
173:   of the array (xa[start],xa[start+stride],xa[start+2*stride], ....)

175:   This is useful for computing, say the maximum of the pressure variable when
176:   the pressure is stored (interlaced) with other variables, e.g., density, etc.
177:   This will only work if the desire subvector is a stride subvector.

179: .seealso: `Vec`, `VecMax()`, `VecStrideNorm()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMin()`
180: @*/
181: PetscErrorCode VecStrideMax(Vec v, PetscInt start, PetscInt *idex, PetscReal *nrm)
182: {
183:   PetscInt           i, n, bs, id = -1;
184:   const PetscScalar *x;
185:   PetscReal          max = PETSC_MIN_REAL;

187:   PetscFunctionBegin;
190:   PetscAssertPointer(nrm, 4);
191:   PetscCall(VecGetLocalSize(v, &n));
192:   PetscCall(VecGetBlockSize(v, &bs));
193:   PetscCheck(start >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative start %" PetscInt_FMT, start);
194:   PetscCheck(start < bs, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Start of stride subvector (%" PetscInt_FMT ") is too large for stride. Have you set the vector blocksize (%" PetscInt_FMT ") correctly with VecSetBlockSize()?", start, bs);
195:   PetscCall(VecGetArrayRead(v, &x));
196:   for (i = start; i < n; i += bs) {
197:     if (PetscRealPart(x[i]) > max) {
198:       max = PetscRealPart(x[i]);
199:       id  = i;
200:     }
201:   }
202:   PetscCall(VecRestoreArrayRead(v, &x));
203: #if defined(PETSC_HAVE_MPIUNI)
204:   *nrm = max;
205:   if (idex) *idex = id;
206: #else
207:   if (!idex) {
208:     PetscCall(MPIU_Allreduce(&max, nrm, 1, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)v)));
209:   } else {
210:     struct {
211:       PetscReal v;
212:       PetscInt  i;
213:     } in, out;
214:     PetscInt rstart;

216:     PetscCall(VecGetOwnershipRange(v, &rstart, NULL));
217:     in.v = max;
218:     in.i = rstart + id;
219:     PetscCall(MPIU_Allreduce(&in, &out, 1, MPIU_REAL_INT, MPIU_MAXLOC, PetscObjectComm((PetscObject)v)));
220:     *nrm  = out.v;
221:     *idex = out.i;
222:   }
223: #endif
224:   PetscFunctionReturn(PETSC_SUCCESS);
225: }

227: /*@
228:   VecStrideMin - Computes the minimum of subvector of a vector defined
229:   by a starting point and a stride and optionally its location.

231:   Collective

233:   Input Parameters:
234: + v     - the vector
235: - start - starting point of the subvector (defined by a stride)

237:   Output Parameters:
238: + idex - the location where the minimum occurred. (pass `NULL` if not required)
239: - nrm  - the minimum value in the subvector

241:   Level: advanced

243:   Notes:
244:   One must call `VecSetBlockSize()` before this routine to set the stride
245:   information, or use a vector created from a multicomponent `DMDA`.

247:   If xa is the array representing the vector x, then this computes the min
248:   of the array (xa[start],xa[start+stride],xa[start+2*stride], ....)

250:   This is useful for computing, say the minimum of the pressure variable when
251:   the pressure is stored (interlaced) with other variables, e.g., density, etc.
252:   This will only work if the desire subvector is a stride subvector.

254: .seealso: `Vec`, `VecMin()`, `VecStrideNorm()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMax()`
255: @*/
256: PetscErrorCode VecStrideMin(Vec v, PetscInt start, PetscInt *idex, PetscReal *nrm)
257: {
258:   PetscInt           i, n, bs, id = -1;
259:   const PetscScalar *x;
260:   PetscReal          min = PETSC_MAX_REAL;

262:   PetscFunctionBegin;
265:   PetscAssertPointer(nrm, 4);
266:   PetscCall(VecGetLocalSize(v, &n));
267:   PetscCall(VecGetBlockSize(v, &bs));
268:   PetscCheck(start >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative start %" PetscInt_FMT, start);
269:   PetscCheck(start < bs, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Start of stride subvector (%" PetscInt_FMT ") is too large for stride. Have you set the vector blocksize (%" PetscInt_FMT ") correctly with VecSetBlockSize()?", start, bs);
270:   PetscCall(VecGetArrayRead(v, &x));
271:   for (i = start; i < n; i += bs) {
272:     if (PetscRealPart(x[i]) < min) {
273:       min = PetscRealPart(x[i]);
274:       id  = i;
275:     }
276:   }
277:   PetscCall(VecRestoreArrayRead(v, &x));
278: #if defined(PETSC_HAVE_MPIUNI)
279:   *nrm = min;
280:   if (idex) *idex = id;
281: #else
282:   if (!idex) {
283:     PetscCall(MPIU_Allreduce(&min, nrm, 1, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)v)));
284:   } else {
285:     struct {
286:       PetscReal v;
287:       PetscInt  i;
288:     } in, out;
289:     PetscInt rstart;

291:     PetscCall(VecGetOwnershipRange(v, &rstart, NULL));
292:     in.v = min;
293:     in.i = rstart + id;
294:     PetscCall(MPIU_Allreduce(&in, &out, 1, MPIU_REAL_INT, MPIU_MINLOC, PetscObjectComm((PetscObject)v)));
295:     *nrm  = out.v;
296:     *idex = out.i;
297:   }
298: #endif
299:   PetscFunctionReturn(PETSC_SUCCESS);
300: }

302: /*@
303:   VecStrideSum - Computes the sum of subvector of a vector defined
304:   by a starting point and a stride.

306:   Collective

308:   Input Parameters:
309: + v     - the vector
310: - start - starting point of the subvector (defined by a stride)

312:   Output Parameter:
313: . sum - the sum

315:   Level: advanced

317:   Notes:
318:   One must call `VecSetBlockSize()` before this routine to set the stride
319:   information, or use a vector created from a multicomponent `DMDA`.

321:   If x is the array representing the vector x then this computes the sum
322:   of the array (x[start],x[start+stride],x[start+2*stride], ....)

324: .seealso: `Vec`, `VecSum()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`
325: @*/
326: PetscErrorCode VecStrideSum(Vec v, PetscInt start, PetscScalar *sum)
327: {
328:   PetscInt           i, n, bs;
329:   const PetscScalar *x;
330:   PetscScalar        local_sum = 0.0;

332:   PetscFunctionBegin;
335:   PetscAssertPointer(sum, 3);
336:   PetscCall(VecGetLocalSize(v, &n));
337:   PetscCall(VecGetBlockSize(v, &bs));
338:   PetscCheck(start >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative start %" PetscInt_FMT, start);
339:   PetscCheck(start < bs, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Start of stride subvector (%" PetscInt_FMT ") is too large for stride. Have you set the vector blocksize (%" PetscInt_FMT ") correctly with VecSetBlockSize()?", start, bs);
340:   PetscCall(VecGetArrayRead(v, &x));
341:   for (i = start; i < n; i += bs) local_sum += x[i];
342:   PetscCall(MPIU_Allreduce(&local_sum, sum, 1, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)v)));
343:   PetscCall(VecRestoreArrayRead(v, &x));
344:   PetscFunctionReturn(PETSC_SUCCESS);
345: }

347: /*@
348:   VecStrideScaleAll - Scales the subvectors of a vector defined
349:   by a starting point and a stride.

351:   Logically Collective

353:   Input Parameters:
354: + v      - the vector
355: - scales - values to multiply each subvector entry by

357:   Level: advanced

359:   Notes:
360:   One must call `VecSetBlockSize()` before this routine to set the stride
361:   information, or use a vector created from a multicomponent `DMDA`.

363:   The dimension of scales must be the same as the vector block size

365: .seealso: `Vec`, `VecNorm()`, `VecStrideScale()`, `VecScale()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`
366: @*/
367: PetscErrorCode VecStrideScaleAll(Vec v, const PetscScalar *scales)
368: {
369:   PetscInt     i, j, n, bs;
370:   PetscScalar *x;

372:   PetscFunctionBegin;
374:   PetscAssertPointer(scales, 2);
375:   PetscCall(VecGetLocalSize(v, &n));
376:   PetscCall(VecGetBlockSize(v, &bs));
377:   PetscCall(VecGetArray(v, &x));
378:   /* need to provide optimized code for each bs */
379:   for (i = 0; i < n; i += bs) {
380:     for (j = 0; j < bs; j++) x[i + j] *= scales[j];
381:   }
382:   PetscCall(VecRestoreArray(v, &x));
383:   PetscFunctionReturn(PETSC_SUCCESS);
384: }

386: /*@
387:   VecStrideNormAll - Computes the norms of subvectors of a vector defined
388:   by a starting point and a stride.

390:   Collective

392:   Input Parameters:
393: + v     - the vector
394: - ntype - type of norm, one of `NORM_1`, `NORM_2`, `NORM_INFINITY`

396:   Output Parameter:
397: . nrm - the norms

399:   Level: advanced

401:   Notes:
402:   One must call `VecSetBlockSize()` before this routine to set the stride
403:   information, or use a vector created from a multicomponent `DMDA`.

405:   If x is the array representing the vector x then this computes the norm
406:   of the array (x[start],x[start+stride],x[start+2*stride], ....) for each start < stride

408:   The dimension of nrm must be the same as the vector block size

410:   This will only work if the desire subvector is a stride subvector

412: .seealso: `Vec`, `VecNorm()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`
413: @*/
414: PetscErrorCode VecStrideNormAll(Vec v, NormType ntype, PetscReal nrm[])
415: {
416:   PetscInt           i, j, n, bs;
417:   const PetscScalar *x;
418:   PetscReal          tnorm[128];
419:   MPI_Comm           comm;

421:   PetscFunctionBegin;
424:   PetscAssertPointer(nrm, 3);
425:   PetscCall(VecGetLocalSize(v, &n));
426:   PetscCall(VecGetArrayRead(v, &x));
427:   PetscCall(PetscObjectGetComm((PetscObject)v, &comm));

429:   PetscCall(VecGetBlockSize(v, &bs));
430:   PetscCheck(bs <= 128, comm, PETSC_ERR_SUP, "Currently supports only blocksize up to 128");

432:   if (ntype == NORM_2) {
433:     PetscScalar sum[128];
434:     for (j = 0; j < bs; j++) sum[j] = 0.0;
435:     for (i = 0; i < n; i += bs) {
436:       for (j = 0; j < bs; j++) sum[j] += x[i + j] * (PetscConj(x[i + j]));
437:     }
438:     for (j = 0; j < bs; j++) tnorm[j] = PetscRealPart(sum[j]);

440:     PetscCall(MPIU_Allreduce(tnorm, nrm, bs, MPIU_REAL, MPIU_SUM, comm));
441:     for (j = 0; j < bs; j++) nrm[j] = PetscSqrtReal(nrm[j]);
442:   } else if (ntype == NORM_1) {
443:     for (j = 0; j < bs; j++) tnorm[j] = 0.0;

445:     for (i = 0; i < n; i += bs) {
446:       for (j = 0; j < bs; j++) tnorm[j] += PetscAbsScalar(x[i + j]);
447:     }

449:     PetscCall(MPIU_Allreduce(tnorm, nrm, bs, MPIU_REAL, MPIU_SUM, comm));
450:   } else if (ntype == NORM_INFINITY) {
451:     PetscReal tmp;
452:     for (j = 0; j < bs; j++) tnorm[j] = 0.0;

454:     for (i = 0; i < n; i += bs) {
455:       for (j = 0; j < bs; j++) {
456:         if ((tmp = PetscAbsScalar(x[i + j])) > tnorm[j]) tnorm[j] = tmp;
457:         /* check special case of tmp == NaN */
458:         if (tmp != tmp) {
459:           tnorm[j] = tmp;
460:           break;
461:         }
462:       }
463:     }
464:     PetscCall(MPIU_Allreduce(tnorm, nrm, bs, MPIU_REAL, MPIU_MAX, comm));
465:   } else SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown norm type");
466:   PetscCall(VecRestoreArrayRead(v, &x));
467:   PetscFunctionReturn(PETSC_SUCCESS);
468: }

470: /*@
471:   VecStrideMaxAll - Computes the maximums of subvectors of a vector defined
472:   by a starting point and a stride and optionally its location.

474:   Collective

476:   Input Parameter:
477: . v - the vector

479:   Output Parameters:
480: + idex - the location where the maximum occurred (not supported, pass `NULL`,
481:            if you need this, send mail to petsc-maint@mcs.anl.gov to request it)
482: - nrm  - the maximum values of each subvector

484:   Level: advanced

486:   Notes:
487:   One must call `VecSetBlockSize()` before this routine to set the stride
488:   information, or use a vector created from a multicomponent `DMDA`.

490:   The dimension of nrm must be the same as the vector block size

492: .seealso: `Vec`, `VecMax()`, `VecStrideNorm()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMin()`
493: @*/
494: PetscErrorCode VecStrideMaxAll(Vec v, PetscInt idex[], PetscReal nrm[])
495: {
496:   PetscInt           i, j, n, bs;
497:   const PetscScalar *x;
498:   PetscReal          max[128], tmp;
499:   MPI_Comm           comm;

501:   PetscFunctionBegin;
503:   PetscAssertPointer(nrm, 3);
504:   PetscCheck(!idex, PETSC_COMM_SELF, PETSC_ERR_SUP, "No support yet for returning index; send mail to petsc-maint@mcs.anl.gov asking for it");
505:   PetscCall(VecGetLocalSize(v, &n));
506:   PetscCall(VecGetArrayRead(v, &x));
507:   PetscCall(PetscObjectGetComm((PetscObject)v, &comm));

509:   PetscCall(VecGetBlockSize(v, &bs));
510:   PetscCheck(bs <= 128, comm, PETSC_ERR_SUP, "Currently supports only blocksize up to 128");

512:   if (!n) {
513:     for (j = 0; j < bs; j++) max[j] = PETSC_MIN_REAL;
514:   } else {
515:     for (j = 0; j < bs; j++) max[j] = PetscRealPart(x[j]);

517:     for (i = bs; i < n; i += bs) {
518:       for (j = 0; j < bs; j++) {
519:         if ((tmp = PetscRealPart(x[i + j])) > max[j]) max[j] = tmp;
520:       }
521:     }
522:   }
523:   PetscCall(MPIU_Allreduce(max, nrm, bs, MPIU_REAL, MPIU_MAX, comm));

525:   PetscCall(VecRestoreArrayRead(v, &x));
526:   PetscFunctionReturn(PETSC_SUCCESS);
527: }

529: /*@
530:   VecStrideMinAll - Computes the minimum of subvector of a vector defined
531:   by a starting point and a stride and optionally its location.

533:   Collective

535:   Input Parameter:
536: . v - the vector

538:   Output Parameters:
539: + idex - the location where the minimum occurred (not supported, pass `NULL`,
540:            if you need this, send mail to petsc-maint@mcs.anl.gov to request it)
541: - nrm  - the minimums of each subvector

543:   Level: advanced

545:   Notes:
546:   One must call `VecSetBlockSize()` before this routine to set the stride
547:   information, or use a vector created from a multicomponent `DMDA`.

549:   The dimension of `nrm` must be the same as the vector block size

551: .seealso: `Vec`, `VecMin()`, `VecStrideNorm()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMax()`
552: @*/
553: PetscErrorCode VecStrideMinAll(Vec v, PetscInt idex[], PetscReal nrm[])
554: {
555:   PetscInt           i, n, bs, j;
556:   const PetscScalar *x;
557:   PetscReal          min[128], tmp;
558:   MPI_Comm           comm;

560:   PetscFunctionBegin;
562:   PetscAssertPointer(nrm, 3);
563:   PetscCheck(!idex, PETSC_COMM_SELF, PETSC_ERR_SUP, "No support yet for returning index; send mail to petsc-maint@mcs.anl.gov asking for it");
564:   PetscCall(VecGetLocalSize(v, &n));
565:   PetscCall(VecGetArrayRead(v, &x));
566:   PetscCall(PetscObjectGetComm((PetscObject)v, &comm));

568:   PetscCall(VecGetBlockSize(v, &bs));
569:   PetscCheck(bs <= 128, comm, PETSC_ERR_SUP, "Currently supports only blocksize up to 128");

571:   if (!n) {
572:     for (j = 0; j < bs; j++) min[j] = PETSC_MAX_REAL;
573:   } else {
574:     for (j = 0; j < bs; j++) min[j] = PetscRealPart(x[j]);

576:     for (i = bs; i < n; i += bs) {
577:       for (j = 0; j < bs; j++) {
578:         if ((tmp = PetscRealPart(x[i + j])) < min[j]) min[j] = tmp;
579:       }
580:     }
581:   }
582:   PetscCall(MPIU_Allreduce(min, nrm, bs, MPIU_REAL, MPIU_MIN, comm));

584:   PetscCall(VecRestoreArrayRead(v, &x));
585:   PetscFunctionReturn(PETSC_SUCCESS);
586: }

588: /*@
589:   VecStrideSumAll - Computes the sums of subvectors of a vector defined by a stride.

591:   Collective

593:   Input Parameter:
594: . v - the vector

596:   Output Parameter:
597: . sums - the sums

599:   Level: advanced

601:   Notes:
602:   One must call `VecSetBlockSize()` before this routine to set the stride
603:   information, or use a vector created from a multicomponent `DMDA`.

605:   If x is the array representing the vector x then this computes the sum
606:   of the array (x[start],x[start+stride],x[start+2*stride], ....) for each start < stride

608: .seealso: `Vec`, `VecSum()`, `VecStrideGather()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`
609: @*/
610: PetscErrorCode VecStrideSumAll(Vec v, PetscScalar sums[])
611: {
612:   PetscInt           i, j, n, bs;
613:   const PetscScalar *x;
614:   PetscScalar        local_sums[128];
615:   MPI_Comm           comm;

617:   PetscFunctionBegin;
619:   PetscAssertPointer(sums, 2);
620:   PetscCall(VecGetLocalSize(v, &n));
621:   PetscCall(VecGetArrayRead(v, &x));
622:   PetscCall(PetscObjectGetComm((PetscObject)v, &comm));

624:   PetscCall(VecGetBlockSize(v, &bs));
625:   PetscCheck(bs <= 128, comm, PETSC_ERR_SUP, "Currently supports only blocksize up to 128");

627:   for (j = 0; j < bs; j++) local_sums[j] = 0.0;
628:   for (i = 0; i < n; i += bs) {
629:     for (j = 0; j < bs; j++) local_sums[j] += x[i + j];
630:   }
631:   PetscCall(MPIU_Allreduce(local_sums, sums, bs, MPIU_SCALAR, MPIU_SUM, comm));

633:   PetscCall(VecRestoreArrayRead(v, &x));
634:   PetscFunctionReturn(PETSC_SUCCESS);
635: }

637: /*----------------------------------------------------------------------------------------------*/
638: /*@
639:   VecStrideGatherAll - Gathers all the single components from a multi-component vector into
640:   separate vectors.

642:   Collective

644:   Input Parameters:
645: + v    - the vector
646: - addv - one of `ADD_VALUES`, `INSERT_VALUES`, `MAX_VALUES`

648:   Output Parameter:
649: . s - the location where the subvectors are stored

651:   Level: advanced

653:   Notes:
654:   One must call `VecSetBlockSize()` before this routine to set the stride
655:   information, or use a vector created from a multicomponent `DMDA`.

657:   If x is the array representing the vector x then this gathers
658:   the arrays (x[start],x[start+stride],x[start+2*stride], ....)
659:   for start=0,1,2,...bs-1

661:   The parallel layout of the vector and the subvector must be the same;
662:   i.e., nlocal of v = stride*(nlocal of s)

664:   Not optimized; could be easily

666: .seealso: `Vec`, `VecStrideNorm()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`, `VecStrideGather()`,
667:           `VecStrideScatterAll()`
668: @*/
669: PetscErrorCode VecStrideGatherAll(Vec v, Vec s[], InsertMode addv)
670: {
671:   PetscInt           i, n, n2, bs, j, k, *bss = NULL, nv, jj, nvc;
672:   PetscScalar      **y;
673:   const PetscScalar *x;

675:   PetscFunctionBegin;
677:   PetscAssertPointer(s, 2);
679:   PetscCall(VecGetLocalSize(v, &n));
680:   PetscCall(VecGetLocalSize(s[0], &n2));
681:   PetscCall(VecGetArrayRead(v, &x));
682:   PetscCall(VecGetBlockSize(v, &bs));
683:   PetscCheck(bs > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Input vector does not have a valid blocksize set");

685:   PetscCall(PetscMalloc2(bs, &y, bs, &bss));
686:   nv  = 0;
687:   nvc = 0;
688:   for (i = 0; i < bs; i++) {
689:     PetscCall(VecGetBlockSize(s[i], &bss[i]));
690:     if (bss[i] < 1) bss[i] = 1; /* if user never set it then assume 1  Re: [PETSC #8241] VecStrideGatherAll */
691:     PetscCall(VecGetArray(s[i], &y[i]));
692:     nvc += bss[i];
693:     nv++;
694:     PetscCheck(nvc <= bs, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Number of subvectors in subvectors > number of vectors in main vector");
695:     if (nvc == bs) break;
696:   }

698:   n = n / bs;

700:   jj = 0;
701:   if (addv == INSERT_VALUES) {
702:     for (j = 0; j < nv; j++) {
703:       for (k = 0; k < bss[j]; k++) {
704:         for (i = 0; i < n; i++) y[j][i * bss[j] + k] = x[bs * i + jj + k];
705:       }
706:       jj += bss[j];
707:     }
708:   } else if (addv == ADD_VALUES) {
709:     for (j = 0; j < nv; j++) {
710:       for (k = 0; k < bss[j]; k++) {
711:         for (i = 0; i < n; i++) y[j][i * bss[j] + k] += x[bs * i + jj + k];
712:       }
713:       jj += bss[j];
714:     }
715: #if !defined(PETSC_USE_COMPLEX)
716:   } else if (addv == MAX_VALUES) {
717:     for (j = 0; j < nv; j++) {
718:       for (k = 0; k < bss[j]; k++) {
719:         for (i = 0; i < n; i++) y[j][i * bss[j] + k] = PetscMax(y[j][i * bss[j] + k], x[bs * i + jj + k]);
720:       }
721:       jj += bss[j];
722:     }
723: #endif
724:   } else SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown norm type");

726:   PetscCall(VecRestoreArrayRead(v, &x));
727:   for (i = 0; i < nv; i++) PetscCall(VecRestoreArray(s[i], &y[i]));

729:   PetscCall(PetscFree2(y, bss));
730:   PetscFunctionReturn(PETSC_SUCCESS);
731: }

733: /*@
734:   VecStrideScatterAll - Scatters all the single components from separate vectors into
735:   a multi-component vector.

737:   Collective

739:   Input Parameters:
740: + s    - the location where the subvectors are stored
741: - addv - one of `ADD_VALUES`, `INSERT_VALUES`, `MAX_VALUES`

743:   Output Parameter:
744: . v - the multicomponent vector

746:   Level: advanced

748:   Notes:
749:   One must call `VecSetBlockSize()` before this routine to set the stride
750:   information, or use a vector created from a multicomponent `DMDA`.

752:   The parallel layout of the vector and the subvector must be the same;
753:   i.e., nlocal of v = stride*(nlocal of s)

755:   Not optimized; could be easily

757: .seealso: `Vec`, `VecStrideNorm()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`, `VecStrideGather()`,

759: @*/
760: PetscErrorCode VecStrideScatterAll(Vec s[], Vec v, InsertMode addv)
761: {
762:   PetscInt            i, n, n2, bs, j, jj, k, *bss = NULL, nv, nvc;
763:   PetscScalar        *x;
764:   PetscScalar const **y;

766:   PetscFunctionBegin;
768:   PetscAssertPointer(s, 1);
770:   PetscCall(VecGetLocalSize(v, &n));
771:   PetscCall(VecGetLocalSize(s[0], &n2));
772:   PetscCall(VecGetArray(v, &x));
773:   PetscCall(VecGetBlockSize(v, &bs));
774:   PetscCheck(bs > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Input vector does not have a valid blocksize set");

776:   PetscCall(PetscMalloc2(bs, (PetscScalar ***)&y, bs, &bss));
777:   nv  = 0;
778:   nvc = 0;
779:   for (i = 0; i < bs; i++) {
780:     PetscCall(VecGetBlockSize(s[i], &bss[i]));
781:     if (bss[i] < 1) bss[i] = 1; /* if user never set it then assume 1  Re: [PETSC #8241] VecStrideGatherAll */
782:     PetscCall(VecGetArrayRead(s[i], &y[i]));
783:     nvc += bss[i];
784:     nv++;
785:     PetscCheck(nvc <= bs, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Number of subvectors in subvectors > number of vectors in main vector");
786:     if (nvc == bs) break;
787:   }

789:   n = n / bs;

791:   jj = 0;
792:   if (addv == INSERT_VALUES) {
793:     for (j = 0; j < nv; j++) {
794:       for (k = 0; k < bss[j]; k++) {
795:         for (i = 0; i < n; i++) x[bs * i + jj + k] = y[j][i * bss[j] + k];
796:       }
797:       jj += bss[j];
798:     }
799:   } else if (addv == ADD_VALUES) {
800:     for (j = 0; j < nv; j++) {
801:       for (k = 0; k < bss[j]; k++) {
802:         for (i = 0; i < n; i++) x[bs * i + jj + k] += y[j][i * bss[j] + k];
803:       }
804:       jj += bss[j];
805:     }
806: #if !defined(PETSC_USE_COMPLEX)
807:   } else if (addv == MAX_VALUES) {
808:     for (j = 0; j < nv; j++) {
809:       for (k = 0; k < bss[j]; k++) {
810:         for (i = 0; i < n; i++) x[bs * i + jj + k] = PetscMax(x[bs * i + jj + k], y[j][i * bss[j] + k]);
811:       }
812:       jj += bss[j];
813:     }
814: #endif
815:   } else SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown norm type");

817:   PetscCall(VecRestoreArray(v, &x));
818:   for (i = 0; i < nv; i++) PetscCall(VecRestoreArrayRead(s[i], &y[i]));
819:   PetscCall(PetscFree2(*(PetscScalar ***)&y, bss));
820:   PetscFunctionReturn(PETSC_SUCCESS);
821: }

823: /*@
824:   VecStrideGather - Gathers a single component from a multi-component vector into
825:   another vector.

827:   Collective

829:   Input Parameters:
830: + v     - the vector
831: . start - starting point of the subvector (defined by a stride)
832: - addv  - one of `ADD_VALUES`, `INSERT_VALUES`, `MAX_VALUES`

834:   Output Parameter:
835: . s - the location where the subvector is stored

837:   Level: advanced

839:   Notes:
840:   One must call `VecSetBlockSize()` before this routine to set the stride
841:   information, or use a vector created from a multicomponent `DMDA`.

843:   If x is the array representing the vector x then this gathers
844:   the array (x[start],x[start+stride],x[start+2*stride], ....)

846:   The parallel layout of the vector and the subvector must be the same;
847:   i.e., nlocal of v = stride*(nlocal of s)

849:   Not optimized; could be easily

851: .seealso: `Vec`, `VecStrideNorm()`, `VecStrideScatter()`, `VecStrideMin()`, `VecStrideMax()`, `VecStrideGatherAll()`,
852:           `VecStrideScatterAll()`
853: @*/
854: PetscErrorCode VecStrideGather(Vec v, PetscInt start, Vec s, InsertMode addv)
855: {
856:   PetscFunctionBegin;
860:   PetscCheck(start >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative start %" PetscInt_FMT, start);
861:   PetscCheck(start < v->map->bs, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Start of stride subvector (%" PetscInt_FMT ") is too large for stride. Have you set the vector blocksize (%" PetscInt_FMT ") correctly with VecSetBlockSize()?", start,
862:              v->map->bs);
863:   PetscUseTypeMethod(v, stridegather, start, s, addv);
864:   PetscFunctionReturn(PETSC_SUCCESS);
865: }

867: /*@
868:   VecStrideScatter - Scatters a single component from a vector into a multi-component vector.

870:   Collective

872:   Input Parameters:
873: + s     - the single-component vector
874: . start - starting point of the subvector (defined by a stride)
875: - addv  - one of `ADD_VALUES`, `INSERT_VALUES`, `MAX_VALUES`

877:   Output Parameter:
878: . v - the location where the subvector is scattered (the multi-component vector)

880:   Level: advanced

882:   Notes:
883:   One must call `VecSetBlockSize()` on the multi-component vector before this
884:   routine to set the stride  information, or use a vector created from a multicomponent `DMDA`.

886:   The parallel layout of the vector and the subvector must be the same;
887:   i.e., nlocal of v = stride*(nlocal of s)

889:   Not optimized; could be easily

891: .seealso: `Vec`, `VecStrideNorm()`, `VecStrideGather()`, `VecStrideMin()`, `VecStrideMax()`, `VecStrideGatherAll()`,
892:           `VecStrideScatterAll()`, `VecStrideSubSetScatter()`, `VecStrideSubSetGather()`
893: @*/
894: PetscErrorCode VecStrideScatter(Vec s, PetscInt start, Vec v, InsertMode addv)
895: {
896:   PetscFunctionBegin;
900:   PetscCheck(start >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative start %" PetscInt_FMT, start);
901:   PetscCheck(start < v->map->bs, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Start of stride subvector (%" PetscInt_FMT ") is too large for stride. Have you set the vector blocksize (%" PetscInt_FMT ") correctly with VecSetBlockSize()?", start,
902:              v->map->bs);
903:   PetscCall((*v->ops->stridescatter)(s, start, v, addv));
904:   PetscFunctionReturn(PETSC_SUCCESS);
905: }

907: /*@
908:   VecStrideSubSetGather - Gathers a subset of components from a multi-component vector into
909:   another vector.

911:   Collective

913:   Input Parameters:
914: + v    - the vector
915: . nidx - the number of indices
916: . idxv - the indices of the components 0 <= idxv[0] ...idxv[nidx-1] < bs(v), they need not be sorted
917: . idxs - the indices of the components 0 <= idxs[0] ...idxs[nidx-1] < bs(s), they need not be sorted, may be null if nidx == bs(s) or is `PETSC_DETERMINE`
918: - addv - one of `ADD_VALUES`, `INSERT_VALUES`, `MAX_VALUES`

920:   Output Parameter:
921: . s - the location where the subvector is stored

923:   Level: advanced

925:   Notes:
926:   One must call `VecSetBlockSize()` on both vectors before this routine to set the stride
927:   information, or use a vector created from a multicomponent `DMDA`.

929:   The parallel layout of the vector and the subvector must be the same;

931:   Not optimized; could be easily

933: .seealso: `Vec`, `VecStrideNorm()`, `VecStrideScatter()`, `VecStrideGather()`, `VecStrideSubSetScatter()`, `VecStrideMin()`, `VecStrideMax()`, `VecStrideGatherAll()`,
934:           `VecStrideScatterAll()`
935: @*/
936: PetscErrorCode VecStrideSubSetGather(Vec v, PetscInt nidx, const PetscInt idxv[], const PetscInt idxs[], Vec s, InsertMode addv)
937: {
938:   PetscFunctionBegin;
941:   if (nidx == PETSC_DETERMINE) nidx = s->map->bs;
942:   PetscUseTypeMethod(v, stridesubsetgather, nidx, idxv, idxs, s, addv);
943:   PetscFunctionReturn(PETSC_SUCCESS);
944: }

946: /*@
947:   VecStrideSubSetScatter - Scatters components from a vector into a subset of components of a multi-component vector.

949:   Collective

951:   Input Parameters:
952: + s    - the smaller-component vector
953: . nidx - the number of indices in idx
954: . idxs - the indices of the components in the smaller-component vector, 0 <= idxs[0] ...idxs[nidx-1] < bs(s) they need not be sorted, may be null if nidx == bs(s) or is `PETSC_DETERMINE`
955: . idxv - the indices of the components in the larger-component vector, 0 <= idx[0] ...idx[nidx-1] < bs(v) they need not be sorted
956: - addv - one of `ADD_VALUES`, `INSERT_VALUES`, `MAX_VALUES`

958:   Output Parameter:
959: . v - the location where the subvector is into scattered (the multi-component vector)

961:   Level: advanced

963:   Notes:
964:   One must call `VecSetBlockSize()` on the vectors before this
965:   routine to set the stride  information, or use a vector created from a multicomponent `DMDA`.

967:   The parallel layout of the vector and the subvector must be the same;

969:   Not optimized; could be easily

971: .seealso: `Vec`, `VecStrideNorm()`, `VecStrideGather()`, `VecStrideSubSetGather()`, `VecStrideMin()`, `VecStrideMax()`, `VecStrideGatherAll()`,
972:           `VecStrideScatterAll()`
973: @*/
974: PetscErrorCode VecStrideSubSetScatter(Vec s, PetscInt nidx, const PetscInt idxs[], const PetscInt idxv[], Vec v, InsertMode addv)
975: {
976:   PetscFunctionBegin;
979:   if (nidx == PETSC_DETERMINE) nidx = s->map->bs;
980:   PetscCall((*v->ops->stridesubsetscatter)(s, nidx, idxs, idxv, v, addv));
981:   PetscFunctionReturn(PETSC_SUCCESS);
982: }

984: PetscErrorCode VecStrideGather_Default(Vec v, PetscInt start, Vec s, InsertMode addv)
985: {
986:   PetscInt           i, n, bs, ns;
987:   const PetscScalar *x;
988:   PetscScalar       *y;

990:   PetscFunctionBegin;
991:   PetscCall(VecGetLocalSize(v, &n));
992:   PetscCall(VecGetLocalSize(s, &ns));
993:   PetscCall(VecGetArrayRead(v, &x));
994:   PetscCall(VecGetArray(s, &y));

996:   bs = v->map->bs;
997:   PetscCheck(n == ns * bs, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Subvector length * blocksize %" PetscInt_FMT " not correct for gather from original vector %" PetscInt_FMT, ns * bs, n);
998:   x += start;
999:   n = n / bs;

1001:   if (addv == INSERT_VALUES) {
1002:     for (i = 0; i < n; i++) y[i] = x[bs * i];
1003:   } else if (addv == ADD_VALUES) {
1004:     for (i = 0; i < n; i++) y[i] += x[bs * i];
1005: #if !defined(PETSC_USE_COMPLEX)
1006:   } else if (addv == MAX_VALUES) {
1007:     for (i = 0; i < n; i++) y[i] = PetscMax(y[i], x[bs * i]);
1008: #endif
1009:   } else SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown norm type");

1011:   PetscCall(VecRestoreArrayRead(v, &x));
1012:   PetscCall(VecRestoreArray(s, &y));
1013:   PetscFunctionReturn(PETSC_SUCCESS);
1014: }

1016: PetscErrorCode VecStrideScatter_Default(Vec s, PetscInt start, Vec v, InsertMode addv)
1017: {
1018:   PetscInt           i, n, bs, ns;
1019:   PetscScalar       *x;
1020:   const PetscScalar *y;

1022:   PetscFunctionBegin;
1023:   PetscCall(VecGetLocalSize(v, &n));
1024:   PetscCall(VecGetLocalSize(s, &ns));
1025:   PetscCall(VecGetArray(v, &x));
1026:   PetscCall(VecGetArrayRead(s, &y));

1028:   bs = v->map->bs;
1029:   PetscCheck(n == ns * bs, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Subvector length * blocksize %" PetscInt_FMT " not correct for scatter to multicomponent vector %" PetscInt_FMT, ns * bs, n);
1030:   x += start;
1031:   n = n / bs;

1033:   if (addv == INSERT_VALUES) {
1034:     for (i = 0; i < n; i++) x[bs * i] = y[i];
1035:   } else if (addv == ADD_VALUES) {
1036:     for (i = 0; i < n; i++) x[bs * i] += y[i];
1037: #if !defined(PETSC_USE_COMPLEX)
1038:   } else if (addv == MAX_VALUES) {
1039:     for (i = 0; i < n; i++) x[bs * i] = PetscMax(y[i], x[bs * i]);
1040: #endif
1041:   } else SETERRQ(PetscObjectComm((PetscObject)s), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown norm type");

1043:   PetscCall(VecRestoreArray(v, &x));
1044:   PetscCall(VecRestoreArrayRead(s, &y));
1045:   PetscFunctionReturn(PETSC_SUCCESS);
1046: }

1048: PetscErrorCode VecStrideSubSetGather_Default(Vec v, PetscInt nidx, const PetscInt idxv[], const PetscInt idxs[], Vec s, InsertMode addv)
1049: {
1050:   PetscInt           i, j, n, bs, bss, ns;
1051:   const PetscScalar *x;
1052:   PetscScalar       *y;

1054:   PetscFunctionBegin;
1055:   PetscCall(VecGetLocalSize(v, &n));
1056:   PetscCall(VecGetLocalSize(s, &ns));
1057:   PetscCall(VecGetArrayRead(v, &x));
1058:   PetscCall(VecGetArray(s, &y));

1060:   bs  = v->map->bs;
1061:   bss = s->map->bs;
1062:   n   = n / bs;

1064:   if (PetscDefined(USE_DEBUG)) {
1065:     PetscCheck(n == ns / bss, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Incompatible layout of vectors");
1066:     for (j = 0; j < nidx; j++) {
1067:       PetscCheck(idxv[j] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "idx[%" PetscInt_FMT "] %" PetscInt_FMT " is negative", j, idxv[j]);
1068:       PetscCheck(idxv[j] < bs, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "idx[%" PetscInt_FMT "] %" PetscInt_FMT " is greater than or equal to vector blocksize %" PetscInt_FMT, j, idxv[j], bs);
1069:     }
1070:     PetscCheck(idxs || bss == nidx, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must provide idxs when not gathering into all locations");
1071:   }

1073:   if (addv == INSERT_VALUES) {
1074:     if (!idxs) {
1075:       for (i = 0; i < n; i++) {
1076:         for (j = 0; j < bss; j++) y[bss * i + j] = x[bs * i + idxv[j]];
1077:       }
1078:     } else {
1079:       for (i = 0; i < n; i++) {
1080:         for (j = 0; j < bss; j++) y[bss * i + idxs[j]] = x[bs * i + idxv[j]];
1081:       }
1082:     }
1083:   } else if (addv == ADD_VALUES) {
1084:     if (!idxs) {
1085:       for (i = 0; i < n; i++) {
1086:         for (j = 0; j < bss; j++) y[bss * i + j] += x[bs * i + idxv[j]];
1087:       }
1088:     } else {
1089:       for (i = 0; i < n; i++) {
1090:         for (j = 0; j < bss; j++) y[bss * i + idxs[j]] += x[bs * i + idxv[j]];
1091:       }
1092:     }
1093: #if !defined(PETSC_USE_COMPLEX)
1094:   } else if (addv == MAX_VALUES) {
1095:     if (!idxs) {
1096:       for (i = 0; i < n; i++) {
1097:         for (j = 0; j < bss; j++) y[bss * i + j] = PetscMax(y[bss * i + j], x[bs * i + idxv[j]]);
1098:       }
1099:     } else {
1100:       for (i = 0; i < n; i++) {
1101:         for (j = 0; j < bss; j++) y[bss * i + idxs[j]] = PetscMax(y[bss * i + idxs[j]], x[bs * i + idxv[j]]);
1102:       }
1103:     }
1104: #endif
1105:   } else SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown norm type");

1107:   PetscCall(VecRestoreArrayRead(v, &x));
1108:   PetscCall(VecRestoreArray(s, &y));
1109:   PetscFunctionReturn(PETSC_SUCCESS);
1110: }

1112: PetscErrorCode VecStrideSubSetScatter_Default(Vec s, PetscInt nidx, const PetscInt idxs[], const PetscInt idxv[], Vec v, InsertMode addv)
1113: {
1114:   PetscInt           j, i, n, bs, ns, bss;
1115:   PetscScalar       *x;
1116:   const PetscScalar *y;

1118:   PetscFunctionBegin;
1119:   PetscCall(VecGetLocalSize(v, &n));
1120:   PetscCall(VecGetLocalSize(s, &ns));
1121:   PetscCall(VecGetArray(v, &x));
1122:   PetscCall(VecGetArrayRead(s, &y));

1124:   bs  = v->map->bs;
1125:   bss = s->map->bs;
1126:   n   = n / bs;

1128:   if (PetscDefined(USE_DEBUG)) {
1129:     PetscCheck(n == ns / bss, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Incompatible layout of vectors");
1130:     for (j = 0; j < bss; j++) {
1131:       if (idxs) {
1132:         PetscCheck(idxs[j] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "idx[%" PetscInt_FMT "] %" PetscInt_FMT " is negative", j, idxs[j]);
1133:         PetscCheck(idxs[j] < bs, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "idx[%" PetscInt_FMT "] %" PetscInt_FMT " is greater than or equal to vector blocksize %" PetscInt_FMT, j, idxs[j], bs);
1134:       }
1135:     }
1136:     PetscCheck(idxs || bss == nidx, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must provide idxs when not scattering from all locations");
1137:   }

1139:   if (addv == INSERT_VALUES) {
1140:     if (!idxs) {
1141:       for (i = 0; i < n; i++) {
1142:         for (j = 0; j < bss; j++) x[bs * i + idxv[j]] = y[bss * i + j];
1143:       }
1144:     } else {
1145:       for (i = 0; i < n; i++) {
1146:         for (j = 0; j < bss; j++) x[bs * i + idxv[j]] = y[bss * i + idxs[j]];
1147:       }
1148:     }
1149:   } else if (addv == ADD_VALUES) {
1150:     if (!idxs) {
1151:       for (i = 0; i < n; i++) {
1152:         for (j = 0; j < bss; j++) x[bs * i + idxv[j]] += y[bss * i + j];
1153:       }
1154:     } else {
1155:       for (i = 0; i < n; i++) {
1156:         for (j = 0; j < bss; j++) x[bs * i + idxv[j]] += y[bss * i + idxs[j]];
1157:       }
1158:     }
1159: #if !defined(PETSC_USE_COMPLEX)
1160:   } else if (addv == MAX_VALUES) {
1161:     if (!idxs) {
1162:       for (i = 0; i < n; i++) {
1163:         for (j = 0; j < bss; j++) x[bs * i + idxv[j]] = PetscMax(y[bss * i + j], x[bs * i + idxv[j]]);
1164:       }
1165:     } else {
1166:       for (i = 0; i < n; i++) {
1167:         for (j = 0; j < bss; j++) x[bs * i + idxv[j]] = PetscMax(y[bss * i + idxs[j]], x[bs * i + idxv[j]]);
1168:       }
1169:     }
1170: #endif
1171:   } else SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown norm type");

1173:   PetscCall(VecRestoreArray(v, &x));
1174:   PetscCall(VecRestoreArrayRead(s, &y));
1175:   PetscFunctionReturn(PETSC_SUCCESS);
1176: }

1178: static PetscErrorCode VecApplyUnary_Private(Vec v, PetscDeviceContext dctx, const char async_name[], PetscErrorCode (*unary_op)(Vec), PetscScalar (*UnaryFunc)(PetscScalar))
1179: {
1180:   PetscFunctionBegin;
1182:   PetscCall(VecSetErrorIfLocked(v, 1));
1183:   if (dctx) {
1184:     PetscErrorCode (*unary_op_async)(Vec, PetscDeviceContext);

1186:     PetscCall(PetscObjectQueryFunction((PetscObject)v, async_name, &unary_op_async));
1187:     if (unary_op_async) {
1188:       PetscCall((*unary_op_async)(v, dctx));
1189:       PetscFunctionReturn(PETSC_SUCCESS);
1190:     }
1191:   }
1192:   if (unary_op) {
1194:     PetscCall((*unary_op)(v));
1195:   } else {
1196:     PetscInt     n;
1197:     PetscScalar *x;

1200:     PetscCall(VecGetLocalSize(v, &n));
1201:     PetscCall(VecGetArray(v, &x));
1202:     for (PetscInt i = 0; i < n; ++i) x[i] = UnaryFunc(x[i]);
1203:     PetscCall(VecRestoreArray(v, &x));
1204:   }
1205:   PetscFunctionReturn(PETSC_SUCCESS);
1206: }

1208: static PetscScalar ScalarReciprocal_Function(PetscScalar x)
1209: {
1210:   const PetscScalar zero = 0.0;

1212:   return x == zero ? zero : ((PetscScalar)1.0) / x;
1213: }

1215: PetscErrorCode VecReciprocalAsync_Private(Vec v, PetscDeviceContext dctx)
1216: {
1217:   PetscFunctionBegin;
1218:   PetscCall(VecApplyUnary_Private(v, dctx, VecAsyncFnName(Reciprocal), v->ops->reciprocal, ScalarReciprocal_Function));
1219:   PetscFunctionReturn(PETSC_SUCCESS);
1220: }

1222: PetscErrorCode VecReciprocal_Default(Vec v)
1223: {
1224:   PetscFunctionBegin;
1225:   PetscCall(VecApplyUnary_Private(v, NULL, NULL, NULL, ScalarReciprocal_Function));
1226:   PetscFunctionReturn(PETSC_SUCCESS);
1227: }

1229: static PetscScalar ScalarExp_Function(PetscScalar x)
1230: {
1231:   return PetscExpScalar(x);
1232: }

1234: PetscErrorCode VecExpAsync_Private(Vec v, PetscDeviceContext dctx)
1235: {
1236:   PetscFunctionBegin;
1238:   PetscCall(VecApplyUnary_Private(v, dctx, VecAsyncFnName(Exp), v->ops->exp, ScalarExp_Function));
1239:   PetscFunctionReturn(PETSC_SUCCESS);
1240: }

1242: /*@
1243:   VecExp - Replaces each component of a vector by e^x_i

1245:   Not Collective

1247:   Input Parameter:
1248: . v - The vector

1250:   Output Parameter:
1251: . v - The vector of exponents

1253:   Level: beginner

1255: .seealso: `Vec`, `VecLog()`, `VecAbs()`, `VecSqrtAbs()`, `VecReciprocal()`

1257: @*/
1258: PetscErrorCode VecExp(Vec v)
1259: {
1260:   PetscFunctionBegin;
1261:   PetscCall(VecExpAsync_Private(v, NULL));
1262:   PetscFunctionReturn(PETSC_SUCCESS);
1263: }

1265: static PetscScalar ScalarLog_Function(PetscScalar x)
1266: {
1267:   return PetscLogScalar(x);
1268: }

1270: PetscErrorCode VecLogAsync_Private(Vec v, PetscDeviceContext dctx)
1271: {
1272:   PetscFunctionBegin;
1274:   PetscCall(VecApplyUnary_Private(v, dctx, VecAsyncFnName(Log), v->ops->log, ScalarLog_Function));
1275:   PetscFunctionReturn(PETSC_SUCCESS);
1276: }

1278: /*@
1279:   VecLog - Replaces each component of a vector by log(x_i), the natural logarithm

1281:   Not Collective

1283:   Input Parameter:
1284: . v - The vector

1286:   Output Parameter:
1287: . v - The vector of logs

1289:   Level: beginner

1291: .seealso: `Vec`, `VecExp()`, `VecAbs()`, `VecSqrtAbs()`, `VecReciprocal()`

1293: @*/
1294: PetscErrorCode VecLog(Vec v)
1295: {
1296:   PetscFunctionBegin;
1297:   PetscCall(VecLogAsync_Private(v, NULL));
1298:   PetscFunctionReturn(PETSC_SUCCESS);
1299: }

1301: static PetscScalar ScalarAbs_Function(PetscScalar x)
1302: {
1303:   return PetscAbsScalar(x);
1304: }

1306: PetscErrorCode VecAbsAsync_Private(Vec v, PetscDeviceContext dctx)
1307: {
1308:   PetscFunctionBegin;
1310:   PetscCall(VecApplyUnary_Private(v, dctx, VecAsyncFnName(Abs), v->ops->abs, ScalarAbs_Function));
1311:   PetscFunctionReturn(PETSC_SUCCESS);
1312: }

1314: /*@
1315:   VecAbs - Replaces every element in a vector with its absolute value.

1317:   Logically Collective

1319:   Input Parameter:
1320: . v - the vector

1322:   Level: intermediate

1324: .seealso: `Vec`, `VecExp()`, `VecSqrtAbs()`, `VecReciprocal()`, `VecLog()`
1325: @*/
1326: PetscErrorCode VecAbs(Vec v)
1327: {
1328:   PetscFunctionBegin;
1329:   PetscCall(VecAbsAsync_Private(v, NULL));
1330:   PetscFunctionReturn(PETSC_SUCCESS);
1331: }

1333: static PetscScalar ScalarConjugate_Function(PetscScalar x)
1334: {
1335:   return PetscConj(x);
1336: }

1338: PetscErrorCode VecConjugateAsync_Private(Vec v, PetscDeviceContext dctx)
1339: {
1340:   PetscFunctionBegin;
1342:   if (PetscDefined(USE_COMPLEX)) PetscCall(VecApplyUnary_Private(v, dctx, VecAsyncFnName(Conjugate), v->ops->conjugate, ScalarConjugate_Function));
1343:   PetscFunctionReturn(PETSC_SUCCESS);
1344: }

1346: /*@
1347:   VecConjugate - Conjugates a vector. That is, replace every entry in a vector with its complex conjugate

1349:   Logically Collective

1351:   Input Parameter:
1352: . x - the vector

1354:   Level: intermediate

1356: .seealso: [](ch_vectors), `Vec`, `VecSet()`
1357: @*/
1358: PetscErrorCode VecConjugate(Vec x)
1359: {
1360:   PetscFunctionBegin;
1361:   PetscCall(VecConjugateAsync_Private(x, NULL));
1362:   PetscFunctionReturn(PETSC_SUCCESS);
1363: }

1365: static PetscScalar ScalarSqrtAbs_Function(PetscScalar x)
1366: {
1367:   return PetscSqrtScalar(ScalarAbs_Function(x));
1368: }

1370: PetscErrorCode VecSqrtAbsAsync_Private(Vec v, PetscDeviceContext dctx)
1371: {
1372:   PetscFunctionBegin;
1374:   PetscCall(VecApplyUnary_Private(v, dctx, VecAsyncFnName(SqrtAbs), v->ops->sqrt, ScalarSqrtAbs_Function));
1375:   PetscFunctionReturn(PETSC_SUCCESS);
1376: }

1378: /*@
1379:   VecSqrtAbs - Replaces each component of a vector by the square root of its magnitude.

1381:   Not Collective

1383:   Input Parameter:
1384: . v - The vector

1386:   Level: beginner

1388:   Note:
1389:   The actual function is sqrt(|x_i|)

1391: .seealso: `Vec`, `VecLog()`, `VecExp()`, `VecReciprocal()`, `VecAbs()`

1393: @*/
1394: PetscErrorCode VecSqrtAbs(Vec v)
1395: {
1396:   PetscFunctionBegin;
1397:   PetscCall(VecSqrtAbsAsync_Private(v, NULL));
1398:   PetscFunctionReturn(PETSC_SUCCESS);
1399: }

1401: static PetscScalar ScalarImaginaryPart_Function(PetscScalar x)
1402: {
1403:   const PetscReal imag = PetscImaginaryPart(x);

1405: #if PetscDefined(USE_COMPLEX)
1406:   return PetscCMPLX(imag, 0.0);
1407: #else
1408:   return imag;
1409: #endif
1410: }

1412: /*@
1413:   VecImaginaryPart - Replaces a complex vector with its imginary part

1415:   Collective

1417:   Input Parameter:
1418: . v - the vector

1420:   Level: beginner

1422: .seealso: `Vec`, `VecNorm()`, `VecRealPart()`
1423: @*/
1424: PetscErrorCode VecImaginaryPart(Vec v)
1425: {
1426:   PetscFunctionBegin;
1428:   PetscCall(VecApplyUnary_Private(v, NULL, NULL, NULL, ScalarImaginaryPart_Function));
1429:   PetscFunctionReturn(PETSC_SUCCESS);
1430: }

1432: static PetscScalar ScalarRealPart_Function(PetscScalar x)
1433: {
1434:   const PetscReal real = PetscRealPart(x);

1436: #if PetscDefined(USE_COMPLEX)
1437:   return PetscCMPLX(real, 0.0);
1438: #else
1439:   return real;
1440: #endif
1441: }

1443: /*@
1444:   VecRealPart - Replaces a complex vector with its real part

1446:   Collective

1448:   Input Parameter:
1449: . v - the vector

1451:   Level: beginner

1453: .seealso: `Vec`, `VecNorm()`, `VecImaginaryPart()`
1454: @*/
1455: PetscErrorCode VecRealPart(Vec v)
1456: {
1457:   PetscFunctionBegin;
1459:   PetscCall(VecApplyUnary_Private(v, NULL, NULL, NULL, ScalarRealPart_Function));
1460:   PetscFunctionReturn(PETSC_SUCCESS);
1461: }

1463: /*@
1464:   VecDotNorm2 - computes the inner product of two vectors and the 2-norm squared of the second vector

1466:   Collective

1468:   Input Parameters:
1469: + s - first vector
1470: - t - second vector

1472:   Output Parameters:
1473: + dp - s'conj(t)
1474: - nm - t'conj(t)

1476:   Level: advanced

1478:   Note:
1479:   conj(x) is the complex conjugate of x when x is complex

1481: .seealso: `Vec`, `VecDot()`, `VecNorm()`, `VecDotBegin()`, `VecNormBegin()`, `VecDotEnd()`, `VecNormEnd()`

1483: @*/
1484: PetscErrorCode VecDotNorm2(Vec s, Vec t, PetscScalar *dp, PetscReal *nm)
1485: {
1486:   PetscScalar work[] = {0.0, 0.0};

1488:   PetscFunctionBegin;
1491:   PetscAssertPointer(dp, 3);
1492:   PetscAssertPointer(nm, 4);
1495:   PetscCheckSameTypeAndComm(s, 1, t, 2);
1496:   PetscCheck(s->map->N == t->map->N, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Incompatible vector global lengths");
1497:   PetscCheck(s->map->n == t->map->n, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Incompatible vector local lengths");

1499:   PetscCall(PetscLogEventBegin(VEC_DotNorm2, s, t, 0, 0));
1500:   if (s->ops->dotnorm2) {
1501:     PetscUseTypeMethod(s, dotnorm2, t, work, work + 1);
1502:   } else {
1503:     const PetscScalar *sx, *tx;
1504:     PetscInt           n;

1506:     PetscCall(VecGetLocalSize(s, &n));
1507:     PetscCall(VecGetArrayRead(s, &sx));
1508:     PetscCall(VecGetArrayRead(t, &tx));
1509:     for (PetscInt i = 0; i < n; ++i) {
1510:       const PetscScalar txconj = PetscConj(tx[i]);

1512:       work[0] += sx[i] * txconj;
1513:       work[1] += tx[i] * txconj;
1514:     }
1515:     PetscCall(VecRestoreArrayRead(t, &tx));
1516:     PetscCall(VecRestoreArrayRead(s, &sx));
1517:     PetscCall(MPIU_Allreduce(MPI_IN_PLACE, work, 2, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)s)));
1518:     PetscCall(PetscLogFlops(4.0 * n));
1519:   }
1520:   PetscCall(PetscLogEventEnd(VEC_DotNorm2, s, t, 0, 0));
1521:   *dp = work[0];
1522:   *nm = PetscRealPart(work[1]);
1523:   PetscFunctionReturn(PETSC_SUCCESS);
1524: }

1526: /*@
1527:   VecSum - Computes the sum of all the components of a vector.

1529:   Collective

1531:   Input Parameter:
1532: . v - the vector

1534:   Output Parameter:
1535: . sum - the result

1537:   Level: beginner

1539: .seealso: `Vec`, `VecMean()`, `VecNorm()`
1540: @*/
1541: PetscErrorCode VecSum(Vec v, PetscScalar *sum)
1542: {
1543:   PetscScalar tmp = 0.0;

1545:   PetscFunctionBegin;
1547:   PetscAssertPointer(sum, 2);
1548:   if (v->ops->sum) {
1549:     PetscUseTypeMethod(v, sum, &tmp);
1550:   } else {
1551:     const PetscScalar *x;
1552:     PetscInt           n;

1554:     PetscCall(VecGetLocalSize(v, &n));
1555:     PetscCall(VecGetArrayRead(v, &x));
1556:     for (PetscInt i = 0; i < n; ++i) tmp += x[i];
1557:     PetscCall(VecRestoreArrayRead(v, &x));
1558:   }
1559:   PetscCall(MPIU_Allreduce(MPI_IN_PLACE, &tmp, 1, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)v)));
1560:   *sum = tmp;
1561:   PetscFunctionReturn(PETSC_SUCCESS);
1562: }

1564: /*@
1565:   VecMean - Computes the arithmetic mean of all the components of a vector.

1567:   Collective

1569:   Input Parameter:
1570: . v - the vector

1572:   Output Parameter:
1573: . mean - the result

1575:   Level: beginner

1577: .seealso: `Vec`, `VecSum()`, `VecNorm()`
1578: @*/
1579: PetscErrorCode VecMean(Vec v, PetscScalar *mean)
1580: {
1581:   PetscInt n;

1583:   PetscFunctionBegin;
1585:   PetscAssertPointer(mean, 2);
1586:   PetscCall(VecGetSize(v, &n));
1587:   PetscCall(VecSum(v, mean));
1588:   *mean /= n;
1589:   PetscFunctionReturn(PETSC_SUCCESS);
1590: }

1592: PetscErrorCode VecShiftAsync_Private(Vec v, PetscScalar shift, PetscDeviceContext dctx)
1593: {
1594:   PetscErrorCode (*shift_async)(Vec, PetscScalar, PetscDeviceContext) = NULL;

1596:   PetscFunctionBegin;
1597:   if (dctx) {
1598:     PetscErrorCode (*shift_async)(Vec, PetscScalar, PetscDeviceContext);

1600:     PetscCall(PetscObjectQueryFunction((PetscObject)v, VecAsyncFnName(Shift), &shift_async));
1601:   }
1602:   if (shift_async) {
1603:     PetscCall((*shift_async)(v, shift, dctx));
1604:   } else if (v->ops->shift) {
1605:     PetscUseTypeMethod(v, shift, shift);
1606:   } else {
1607:     PetscInt     n;
1608:     PetscScalar *x;

1610:     PetscCall(VecGetLocalSize(v, &n));
1611:     PetscCall(VecGetArray(v, &x));
1612:     for (PetscInt i = 0; i < n; ++i) x[i] += shift;
1613:     PetscCall(VecRestoreArray(v, &x));
1614:     PetscCall(PetscLogFlops(n));
1615:   }
1616:   PetscFunctionReturn(PETSC_SUCCESS);
1617: }

1619: /*@
1620:   VecShift - Shifts all of the components of a vector by computing
1621:   `x[i] = x[i] + shift`.

1623:   Logically Collective

1625:   Input Parameters:
1626: + v     - the vector
1627: - shift - the shift

1629:   Level: intermediate

1631: .seealso: `Vec`, `VecISShift()`
1632: @*/
1633: PetscErrorCode VecShift(Vec v, PetscScalar shift)
1634: {
1635:   PetscFunctionBegin;
1638:   PetscCall(VecSetErrorIfLocked(v, 1));
1639:   if (shift == (PetscScalar)0.0) PetscFunctionReturn(PETSC_SUCCESS);
1640:   PetscCall(PetscLogEventBegin(VEC_Shift, v, 0, 0, 0));
1641:   PetscCall(VecShiftAsync_Private(v, shift, NULL));
1642:   PetscCall(PetscLogEventEnd(VEC_Shift, v, 0, 0, 0));
1643:   PetscFunctionReturn(PETSC_SUCCESS);
1644: }

1646: /*@
1647:   VecPermute - Permutes a vector in place using the given ordering.

1649:   Input Parameters:
1650: + x   - The vector
1651: . row - The ordering
1652: - inv - The flag for inverting the permutation

1654:   Level: beginner

1656:   Note:
1657:   This function does not yet support parallel Index Sets with non-local permutations

1659: .seealso: `Vec`, `MatPermute()`
1660: @*/
1661: PetscErrorCode VecPermute(Vec x, IS row, PetscBool inv)
1662: {
1663:   PetscScalar    *array, *newArray;
1664:   const PetscInt *idx;
1665:   PetscInt        i, rstart, rend;

1667:   PetscFunctionBegin;
1670:   PetscCall(VecSetErrorIfLocked(x, 1));
1671:   PetscCall(VecGetOwnershipRange(x, &rstart, &rend));
1672:   PetscCall(ISGetIndices(row, &idx));
1673:   PetscCall(VecGetArray(x, &array));
1674:   PetscCall(PetscMalloc1(x->map->n, &newArray));
1675:   PetscCall(PetscArraycpy(newArray, array, x->map->n));
1676:   if (PetscDefined(USE_DEBUG)) {
1677:     for (i = 0; i < x->map->n; i++) PetscCheck(!(idx[i] < rstart) && !(idx[i] >= rend), PETSC_COMM_SELF, PETSC_ERR_ARG_CORRUPT, "Permutation index %" PetscInt_FMT " is out of bounds: %" PetscInt_FMT, i, idx[i]);
1678:   }
1679:   if (!inv) {
1680:     for (i = 0; i < x->map->n; i++) array[i] = newArray[idx[i] - rstart];
1681:   } else {
1682:     for (i = 0; i < x->map->n; i++) array[idx[i] - rstart] = newArray[i];
1683:   }
1684:   PetscCall(VecRestoreArray(x, &array));
1685:   PetscCall(ISRestoreIndices(row, &idx));
1686:   PetscCall(PetscFree(newArray));
1687:   PetscFunctionReturn(PETSC_SUCCESS);
1688: }

1690: /*@
1691:   VecEqual - Compares two vectors. Returns true if the two vectors are either pointing to the same memory buffer,
1692:   or if the two vectors have the same local and global layout as well as bitwise equality of all entries.
1693:   Does NOT take round-off errors into account.

1695:   Collective

1697:   Input Parameters:
1698: + vec1 - the first vector
1699: - vec2 - the second vector

1701:   Output Parameter:
1702: . flg - `PETSC_TRUE` if the vectors are equal; `PETSC_FALSE` otherwise.

1704:   Level: intermediate

1706: .seealso: `Vec`
1707: @*/
1708: PetscErrorCode VecEqual(Vec vec1, Vec vec2, PetscBool *flg)
1709: {
1710:   const PetscScalar *v1, *v2;
1711:   PetscInt           n1, n2, N1, N2;
1712:   PetscBool          flg1;

1714:   PetscFunctionBegin;
1717:   PetscAssertPointer(flg, 3);
1718:   if (vec1 == vec2) *flg = PETSC_TRUE;
1719:   else {
1720:     PetscCall(VecGetSize(vec1, &N1));
1721:     PetscCall(VecGetSize(vec2, &N2));
1722:     if (N1 != N2) flg1 = PETSC_FALSE;
1723:     else {
1724:       PetscCall(VecGetLocalSize(vec1, &n1));
1725:       PetscCall(VecGetLocalSize(vec2, &n2));
1726:       if (n1 != n2) flg1 = PETSC_FALSE;
1727:       else {
1728:         PetscCall(VecGetArrayRead(vec1, &v1));
1729:         PetscCall(VecGetArrayRead(vec2, &v2));
1730:         PetscCall(PetscArraycmp(v1, v2, n1, &flg1));
1731:         PetscCall(VecRestoreArrayRead(vec1, &v1));
1732:         PetscCall(VecRestoreArrayRead(vec2, &v2));
1733:       }
1734:     }
1735:     /* combine results from all processors */
1736:     PetscCall(MPIU_Allreduce(&flg1, flg, 1, MPIU_BOOL, MPI_MIN, PetscObjectComm((PetscObject)vec1)));
1737:   }
1738:   PetscFunctionReturn(PETSC_SUCCESS);
1739: }

1741: /*@
1742:   VecUniqueEntries - Compute the number of unique entries, and those entries

1744:   Collective

1746:   Input Parameter:
1747: . vec - the vector

1749:   Output Parameters:
1750: + n - The number of unique entries
1751: - e - The entries

1753:   Level: intermediate

1755: .seealso: `Vec`
1756: @*/
1757: PetscErrorCode VecUniqueEntries(Vec vec, PetscInt *n, PetscScalar **e)
1758: {
1759:   const PetscScalar *v;
1760:   PetscScalar       *tmp, *vals;
1761:   PetscMPIInt       *N, *displs, l;
1762:   PetscInt           ng, m, i, j, p;
1763:   PetscMPIInt        size;

1765:   PetscFunctionBegin;
1767:   PetscAssertPointer(n, 2);
1768:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)vec), &size));
1769:   PetscCall(VecGetLocalSize(vec, &m));
1770:   PetscCall(VecGetArrayRead(vec, &v));
1771:   PetscCall(PetscMalloc2(m, &tmp, size, &N));
1772:   for (i = 0, j = 0, l = 0; i < m; ++i) {
1773:     /* Can speed this up with sorting */
1774:     for (j = 0; j < l; ++j) {
1775:       if (v[i] == tmp[j]) break;
1776:     }
1777:     if (j == l) {
1778:       tmp[j] = v[i];
1779:       ++l;
1780:     }
1781:   }
1782:   PetscCall(VecRestoreArrayRead(vec, &v));
1783:   /* Gather serial results */
1784:   PetscCallMPI(MPI_Allgather(&l, 1, MPI_INT, N, 1, MPI_INT, PetscObjectComm((PetscObject)vec)));
1785:   for (p = 0, ng = 0; p < size; ++p) ng += N[p];
1786:   PetscCall(PetscMalloc2(ng, &vals, size + 1, &displs));
1787:   for (p = 1, displs[0] = 0; p <= size; ++p) displs[p] = displs[p - 1] + N[p - 1];
1788:   PetscCallMPI(MPI_Allgatherv(tmp, l, MPIU_SCALAR, vals, N, displs, MPIU_SCALAR, PetscObjectComm((PetscObject)vec)));
1789:   /* Find unique entries */
1790: #ifdef PETSC_USE_COMPLEX
1791:   SETERRQ(PetscObjectComm((PetscObject)vec), PETSC_ERR_SUP, "Does not work with complex numbers");
1792: #else
1793:   *n = displs[size];
1794:   PetscCall(PetscSortRemoveDupsReal(n, (PetscReal *)vals));
1795:   if (e) {
1796:     PetscAssertPointer(e, 3);
1797:     PetscCall(PetscMalloc1(*n, e));
1798:     for (i = 0; i < *n; ++i) (*e)[i] = vals[i];
1799:   }
1800:   PetscCall(PetscFree2(vals, displs));
1801:   PetscCall(PetscFree2(tmp, N));
1802:   PetscFunctionReturn(PETSC_SUCCESS);
1803: #endif
1804: }